1. 项目概述当机器学习势函数遇见分子晶体的“心跳”在材料科学和凝聚态物理领域理解材料的振动性质就好比在聆听材料的“心跳”。这阵心跳——晶格振动或声子——直接决定了材料的热导率、光学响应、电声耦合乃至化学反应活性。对于并苯Acenes这类经典的有机半导体分子晶体如萘、蒽、并四苯、并五苯其振动行为尤为复杂。分子内共价键的强相互作用与分子间微弱的范德华力交织在一起形成了从低频的分子平动/转动到高频的C-H伸缩振动的丰富频谱。更复杂的是当我们将一个“客人”分子如并五苯嵌入到“主人”晶体如萘中形成主客体系统时主客之间的振动耦合、能量转移路径会如何变化这些微观细节对于设计高性能有机光电器件、理解能量传递机制至关重要。传统上要精确计算这些振动性质尤其是超越简谐近似的非谐效应必须依赖第一性原理计算如密度泛函理论。然而DFT计算成本高昂进行一次完整的非谐振动谱分析往往需要数千甚至数万次能量和力的计算对于包含数十上百个原子的超胞体系这几乎是不可承受之重。这就引出了我们今天要深入探讨的核心工具机器学习势函数。MLIPs的核心思想很直观既然DFT计算本质上是输入原子坐标输出能量和原子受力那么我们可以用一个高效的机器学习模型来学习这个复杂的映射关系。一旦模型训练完成它就能以比DFT快数个数量级的速度提供近乎量子力学精度的力和能量从而使得长时间、大尺度的分子动力学模拟成为可能进而计算各种动态性质包括非谐振动光谱。但问题来了如何确保这个“黑箱”模型在广阔且复杂的构型空间尤其是我们关心的振动激发态中依然可靠如何用最少的高质量DFT数据训练出一个既能精准预测已知晶体又能可靠外推到未知主客体系统的通用势函数这正是我们这项工作的挑战与目标。我们选择了MACE架构作为基石。与一些早期MLIP模型不同MACE通过高阶张量消息传递和原子簇展开能够更有效地捕捉多体相互作用这对于描述分子晶体中方向性较强的相互作用至关重要。但仅有先进的架构还不够数据从哪里来、怎么选才是决定模型成败的关键。我们采用了委员会主动学习策略让模型自己“告诉”我们它在哪里还不确定从而引导DFT计算去探索那些最具有信息量的原子构型。最终我们仅用了总计数百次的DFT计算就构建出了一系列针对并苯家族的通用势函数G-MLIPs并成功将其应用于主客体系统的振动性质预测与误差量化。接下来我将拆解整个流程从模型架构选择、主动学习策略设计到振动性质计算与误差分析分享其中的核心思路、实操细节以及我们踩过的一些坑。2. 核心思路与方案选型为什么是MACE主动学习面对并苯晶体振动预测这个问题我们有几个核心需求高精度力误差需在几个meV/Å量级、高数据效率DFT计算成本必须严格控制、良好的可迁移性一个模型能用于多个相似体系并能外推到主客体复合物、以及可靠的不确定性估计我们需要知道预测结果的可信度。2.1 架构选型MACE为何脱颖而出在众多MLIP架构中如Behler-Parrinello神经网络、高斯近似势、DeepMD等我们选择了MACE。这背后有几个关键考量对等变性的严格保证MACE的架构设计严格遵循了物理系统的对称性包括平移、旋转和反演对称性。这意味着无论分子如何旋转模型预测的能量和力都保持不变。这对于振动分析至关重要因为我们需要在分子动力学模拟中采样大量随机取向的构型模型必须对这些对称变换具有不变性。高效的高阶相互作用建模MACE通过消息传递和原子簇展开可以系统地构建高阶原子环境描述符。分子晶体中的相互作用不仅仅是两两原子之间的还涉及三个、四个原子构成的角、二面角等。MACE能够有效地捕捉这些多体效应这对于准确描述分子内弯曲、扭转以及分子间复杂的范德华接触至关重要。与基础模型的衔接近年来出现了如MACE-OFF这样的“基础模型”它们在海量有机分子数据上预训练具有良好的迁移能力。我们的工作可以看作是在此基础上的“专业化”微调。我们验证了即使将MACE的最大角动量L设为0仅使用不变特征通过精心设计的训练也能达到与大型基础模型相媲美的精度这大大降低了模型复杂度和计算开销。注意模型选择不是越新越好。我们对比过VASP内置的on-the-fly MLFF一种基于高斯过程回归的方法。测试发现即使使用完全相同的数据集训练VASP MLFF在预测振动频率尤其是低频分子间模式时误差显著大于MACE见图6。这凸显了模型架构本身在捕捉特定物理相互作用如弱范德华力方面的内在能力差异。2.2 数据策略委员会主动学习的精妙之处有了好的模型还需要好的“饲料”——数据。主动学习的核心思想是让模型参与数据收集过程而不是随机或均匀地采样。我们采用的是委员会查询策略。具体流程如下初始化委员会我们不是训练一个模型而是同时训练8个结构相同但初始权重和训练/验证集划分不同的MACE模型组成一个“委员会”。启动循环首先用一个小的初始数据集例如100个构型训练委员会中的所有模型。探索与查询使用训练好的委员会模型对一个预先通过廉价方法如经典力场或基础模型生成的大规模构型池进行预测。对于池中的每个构型计算委员会成员间预测结果的标准差。这个标准差就是模型不确定性的量化指标。选择不确定性最高的N个例如25个构型认为这些是模型最“吃不准”、信息量最大的区域。标注与增强对这N个高不确定性构型进行昂贵的DFT计算获得精确的能量和力将它们加入训练集。迭代优化用增强后的训练集重新训练委员会模型然后回到步骤3。如此循环直到模型在目标性质对我们来说是Γ点声子频率上的误差收敛。这个策略的高明之处在于高效勘探它自动将DFT计算资源引导到势能面上模型尚未掌握的区域避免了在平滑、简单区域的无谓采样。不确定性量化委员会成员间的分歧直接提供了预测不确定性的估计这为后续的误差传播分析奠定了基础。停止准则我们以Γ点声子频率的预测精度作为主动学习的停止准则确保模型在振动性质这个核心任务上达到要求而不是单纯追求能量和力误差的最小化。2.3 训练目标与损失函数设计我们的损失函数是标准的多任务加权形式Loss w_energy * MSE(E_pred, E_DFT) w_force * MSE(F_pred, F_DFT)其中w_force力权重通常设置为远大于w_energy能量权重例如1000:1。这是因为在分子动力学中力的精度直接决定了轨迹的可靠性而振动频率对力的误差更为敏感。能量误差主要影响绝对结合能但对相对能量和力的影响较小。我们使用的关键超参数总结如下表超参数值说明委员会成员数量8用于不确定性估计和主动学习查询交互层数2消息传递的深度平衡表达能力和效率通道数256特征空间的维度最大角动量 L0仅使用标量不变特征简化模型关联阶数3描述原子环境的多体阶数截断半径 r_max6.0 Å原子相互作用的距离截断力权重 w_force1000损失函数中力的权重能量权重 w_energy10损失函数中能量的权重训练/验证集划分9:1用于监控过拟合实操心得r_max的选择需要谨慎。对于分子晶体6.0 Å足以覆盖第一壳层甚至部分第二壳层的分子间相互作用。设置过大会增加计算成本且可能引入噪声设置过小则无法正确描述长程作用。我们通过测试不同截断半径下对晶格常数和弹性常数的预测来辅助确定。3. 数据集构建与模型训练实战理论再好也需要落地。这部分将详细说明我们从数据准备到模型训练完成的每一步操作。3.1 参考计算与数据池生成一切机器学习的基础是高质量数据。我们的参考计算全部基于DFT使用FHI-aims和VASP软件采用PBE泛函加上MBD多体色散修正来处理关键的范德华相互作用。这一步必须保证一致性因为任何系统误差都会“教”给机器学习模型。数据池的构建策略因体系而异萘我们进行了最彻底的采样。对一个1×2×2的超胞在80K, 120K, 150K, 220K, 295K五个温度下分别进行了50 ps的从头算分子动力学模拟。这确保了训练数据覆盖了从低温到室温的整个热振动范围包含了丰富的非谐运动信息。蒽、并四苯、并五苯为了节省成本我们使用了预训练的通用力场模型MACE-OFF来生成初始构型。例如对蒽在100K和295K下进行AIMD模拟。MACE-OFF虽然精度不及目标DFT但其产生的构型分布是物理合理的为主动学习提供了优秀的初始探索空间。最终我们为每个晶体构建了庞大的未标记数据池见下表而主动学习仅从中挑选了数百个构型进行DFT计算。分子晶体数据池中的结构数量萘36,859蒽200,000并四苯183,123并五苯85,5783.2 主动学习迭代与模型泛化我们训练了四个逐步泛化的模型N-MLIP仅用萘的数据训练专用于萘晶体。G-MLIP 1用萘和蒽的数据训练。G-MLIP 2用萘、蒽和并四苯的数据训练。G-MLIP 3用全部四种并苯萘、蒽、并四苯、并五苯的数据训练。每个模型的训练结构数量如下表所示。可以看到即使是最复杂的G-MLIP 3总训练数据也仅为825个构型。MLIP萘蒽并四苯并五苯总计N-MLIP450---450G-MLIP 1450125--575G-MLIP 2450125125-700G-MLIP 3450125125125825训练过程监控在每次主动学习迭代后我们不仅看训练集和验证集上能量、力的均方根误差更关键的是监控Γ点声子频率的误差。下图图S.18展示了G-MLIP 3训练过程中声子频率最大误差和平均误差的下降情况。当误差曲线趋于平缓且新增数据不再显著改善性能时即可停止迭代。踩坑记录早期我们曾尝试仅以力和能量的验证集误差作为停止准则。结果发现模型在验证集上的力误差虽然很低2 meV/Å但预测某些特定振动模式尤其是低频的分子外模式时仍会出现较大偏差。这警示我们最终的评估标准必须与下游应用目标紧密对齐。对于振动性质预测声子频率误差是最直接的度量。3.3 不确定性校准与传播委员会模型给出的不确定性是原始的需要进行校准才能真实反映预测误差的分布。我们采用了一个缩放因子α来重新标定委员会的不确定性。具体操作是在一个独立的测试集上计算预测值与DFT参考值的差异并与委员会的标准差进行比较通过公式(D1)计算α。对于力我们得到的α约为2.8。这意味着原始的委员会不确定性被低估了需要放大2.8倍才能更好地匹配实际误差。校准后的不确定性可以通过公式(D3)传播到每个委员会成员预测的力上F_i(A) F(A) α * [F_i(A) - F(A)]其中F(A)是委员会平均力F_i(A)是第i个成员的预测力。用这些缩放后的力F_i(A)分别计算物理量如声子频率其成员间的分布就代表了该物理量的不确定性。4. 振动性质计算与误差分析全流程模型训练好后就进入了“用武之地”。我们主要计算了三类振动性质简谐声子谱、非谐振动态密度以及主客体系统的模式投影VDOS。4.1 简谐声子谱计算与误差评估简谐近似下声子频率可以通过计算动力学矩阵的本征值得到。动力学矩阵由原子间的力常数构成而力常数可以通过对超胞中的原子施加微小位移并计算Hellmann-Feynman力来获得有限位移法或直接使用密度泛函微扰理论。在我们的工作中使用phonopy软件包基于2×2×2的超胞用有限位移法计算力常数。关键步骤用训练好的MACE委员会模型为超胞中每个原子在正负方向上施加微小位移通常0.01 Å计算原子受力。由受力变化得到力常数矩阵。对角化动力学矩阵得到所有声子模式的频率和本征矢量极化向量。误差分析我们计算了每个声子模式的预测频率与DFT参考频率的绝对误差和百分比误差。图S.9展示了G-MLIP 3对四种晶体的预测误差。整体趋势是随着模型泛化能力的增强从N-MLIP到G-MLIP 3对于新加入的、更大的并苯分子如并四苯、并五苯预测精度显著提升。例如仅用萘训练的N-MLIP预测并五苯时平均误差高达7.09 cm⁻¹而G-MLIP 3将其降低到了2.85 cm⁻¹。一个有趣的发现是分子间振动模式频率通常低于150 cm⁻¹的误差普遍小于分子内模式。这是因为分子间模式主要由较弱的、相对简单的范德华相互作用主导而分子内模式涉及复杂的键弯曲、键伸缩和面外振动对势能面的细节更为敏感。尽管如此MACE模型在C-H伸缩模式~3100 cm⁻¹这种高频振动上也表现出了出色的精度平均误差在1-3 cm⁻¹以内这非常令人鼓舞。4.2 非谐振动态密度计算简谐近似忽略了声子-声子相互作用无法描述温度效应和谱线展宽。要获得真实的振动光谱必须进行非谐计算。我们通过分子动力学模拟来获取非谐振动态密度。操作流程平衡在目标温度如100K下使用Nose-Hoover控温器进行NVT系综模拟使系统达到热平衡。采样切换至NVE微正则系综进行生产模拟在此期间系统总能量守恒我记录下所有原子的速度随时间演化的轨迹。计算VDOS振动态密度通过速度自相关函数的傅里叶变换得到VDOS(ω) ∝ ∫ v(0)·v(t) exp(-iωt) dt其中v(0)·v(t)是速度自相关函数。我们对每个委员会成员用其缩放后的力分别进行多组独立的MD模拟例如100组每组20 ps计算其VDOS。误差估计总误差σ_total(ω)由两部分构成统计误差σ_stat(ω)源于有限的MD采样时间通过将多次独立MD轨迹的VDOS进行块平均来估计。委员会误差σ_com(ω)源于模型本身的不确定性通过计算不同委员会成员预测的VDOS之间的标准差来估计。 最终某个频率ω处的VDOS值可以表示为VDOS(ω) ± σ_total(ω)。图S.7展示了萘晶体VDOS的误差分布。可以看到在大部分频率区间统计误差是主要贡献但在一些特定的尖锐峰位置模型不确定性会显著增加提示这些区域的预测需要谨慎对待。4.3 主客体系统模式投影与耦合分析这是本项目最精彩的部分。我们将一个并五苯分子嵌入到萘晶体中构建了一个主客体系统。我们不仅计算了整体的VDOS更深入分析了每个简谐本征模式中主体和客体的贡献各是多少。技术核心模式投影VDOS。首先计算主体纯萘晶体的简谐声子本征矢e_j(s)其中s是模式索引j是原子索引。在主体-客体系统的MD模拟中记录所有原子的速度v_j(t)。对于第s个简谐模式计算其投影速度v_s(t) Σ_j √m_j e*_j(s) · v_j(t)这个投影速度可以理解为在时刻t整个系统的原子运动在简谐模式s这个“基矢”上的投影分量。将投影速度按原子归属主体或客体进行拆分v_s(t) v_s^host(t) v_s^guest(t)分别计算主体投影VDOS (C_s^host(ω))、客体投影VDOS (C_s^guest(ω))以及它们的交叉项VDOS (C_s^cross(ω))。总VDOS满足C_s(ω) C_s^host(ω) C_s^guest(ω) C_s^cross(ω)。交叉项C_s^cross(ω)的物理意义极其重要它直接量化了主体和客体振动模式之间的动态耦合强度。如果交叉项为零说明主体和客体的运动在该模式下完全独立交叉项越大说明能量在该模式下在主客体间交换越频繁。我们的分析图5清晰地揭示了一系列混合模式。例如一个频率为~128 cm⁻¹的模式其本征矢在孤立并五苯分子中对应一个~100 cm⁻¹的振动。在嵌入晶体后这个模式不仅客体分子在振动周围的主体萘分子也发生了显著的协同运动见图S.15中的Mode 7。交叉项分析表明这些模式存在强烈的耦合为振动能量从客体到主体的转移提供了直接的通道。4.4 振动寿命提取从非谐VDOS中我们还可以提取振动模式的寿命。对于客体主导的振动模式我们将其投影VDOS的谱峰用洛伦兹线型进行拟合。洛伦兹峰的半高宽FWHM单位cm⁻¹的倒数乘以一个常数因子即给出了该模式的振动寿命单位ps。我们对低频300 cm⁻¹的客体模式进行了分析得到的平均寿命约为3.5 ps图S.16。我们发现一个有趣的规律对于频率低于主体晶体声子截止频率两倍的振动模式其寿命较短而高于该截止频率的模式寿命则较长。这与Dlott等人早期的理论预测一致因为非谐弛豫需要满足能量和动量守恒频率过高的模式缺乏可耦合的低能声子通道从而弛豫更慢。5. 关键发现、局限性与未来展望通过这一系列工作我们得到了几个核心结论数据效率惊人基于MACE和委员会主动学习仅用825个DFT计算就构建了一个能同时精确描述四种并苯晶体及其主客体复合物振动性质的通用势函数G-MLIP 3。这比传统方法所需的数据量少了1-2个数量级。精度与泛化能力兼备G-MLIP 3在预测所有四种晶体简谐声子频率时平均绝对误差在3 cm⁻¹以内。更重要的是它能成功外推到训练集中从未出现过的“萘主体-并五苯客体”系统并精确预测其非谐振动谱和模式耦合。不确定性量化链条贯通我们建立了一套从力预测的不确定性到声子频率不确定性再到非谐VDOS不确定性的完整传播方法。这让我们不仅能给出预测值还能给出置信区间这对于基于模拟的材料设计至关重要。谐波不确定性可作为非谐学习的代理我们发现谐波声子频率的不确定性与从分子动力学得到的非谐性质的不定性存在较强的相关性。这意味着在主动学习过程中我们可以用计算成本低廉的谐波不确定性来指导数据采集以改进对非谐性质的预测这是一个非常实用的策略。当然这项工作也有其局限性长程相互作用当前模型的截断半径为6 Å没有显式包含长程的范德华相互作用。这对于预测晶格膨胀系数等对长程力敏感的性质可能不足。但在预测振动性质尤其是我们关注的频率范围内其精度已经足够。温度外推模型主要在室温及以下温度的数据上训练。将其应用于远高于训练温度的场景时需要谨慎评估其可靠性。应力预测如附录A.1所述我们使用的MACE版本0.3.10在预测应力张量上存在困难这影响了其在可变晶胞模拟中的应用。不过后续的MACE版本已在此方面做了改进。我个人在实际操作中的体会是机器学习势函数的研究正在从“追求通用”转向“追求可靠且高效的专业化”。像MACE-OFF这样的基础模型提供了一个强大的起点但针对特定性质如振动和特定材料体系如有机晶体进行有针对性的主动学习微调是达到应用级精度的关键。未来将这种“基础模型主动学习微调”的模式与更先进的不确定性量化、误差传播方法结合并集成到自动化计算工作流中将极大地加速功能材料的发现与设计进程。对于从事计算材料学的研究者而言掌握这套从数据生成、模型训练到性质分析与误差评估的全链条技能正变得越来越重要。