在计算化学的世界里，科学家们一直在追求一种完美的“替身”——它能够像高精度的量子化学计算方法（如DFT）那样精确地描述分子中原子的相互作用，又能像传统分子力学力场那样高效，从而让我们能够在计算机中进行长时间的分子动力学(Molecular Dynamics, MD)模拟，探索药物如何与靶点结合、材料如何在原子层面变化。这种理想的“替身”被称为机器学习势能函数(Machine-Learned Potentials, MLP)。近二十年来，MLP发展迅速，已成为连接高精度计算与高效模拟的桥梁。
然而，这个“替身”有一个致命的弱点：不稳定。尤其是在高温模拟中，当原子运动加剧、分子构象发生剧烈变化时，许多先进的MLP往往会“崩溃”，导致模拟失败。这就像一个导航软件，虽然在熟悉的道路上表现优异，一旦驶入陌生或崎岖的地形就可能给出错误指令，导致“车祸”。这种不稳定性源于模型在“外推”时的不可靠性，即当模拟进入训练数据未充分覆盖的构象空间时，模型的预测会偏离物理实际，产生荒谬的原子力和能量，最终导致分子结构爆炸（原子间距离过大）或内爆（原子间距离过小）。尽管模型在静态测试集上可以达到极低的平均绝对误差(MAE)，但这并不能保证其在动态模拟中的稳健性。因此，开发能够在高温等极端条件下保持长期稳定的MLP，是计算化学和材料科学领域一个亟待解决的核心挑战。
在此背景下，一项发表在《Communications Chemistry》上的研究取得了突破性进展。该研究团队并未采用主流MLP中将总能量分割为无物理意义的“位点能量”的做法，而是转向了基于量子力学基本原理的“物理信息”原子能量。他们利用高斯过程回归(Gaussian Process Regression, GPR)框架，构建了基于“相互作用量子原子”(Interacting Quantum Atoms, IQA)能量的原子能量模型。IQA能量源于量子理论中的原子(Quantum Theory of Atoms in Molecules, QTAIM)，它将分子的总能量严格分解为每个拓扑原子的贡献，这些能量具有清晰的物理意义（如原子内能、原子间库仑和交换相关能）。研究的关键创新在于发现了高斯过程中一个常被忽视的参数——先验均值函数(prior mean function, m)——对模型稳健性有着决定性影响。传统上，这个函数常被简单地设为零。而本研究通过系统性地将先验均值函数向高原子能量状态偏移，成功地“训练”出了具有非凡“恢复力”的模型。这些模型能够在外推至非物理构型（如被故意拉伸或压缩的化学键）时，自动预测出将分子“拉回”合理状态的恢复力(restoring forces)，从而防止模拟崩溃。
为了验证这一理念，研究人员对四种结构灵活的有机分子——肽封端的甘氨酸(GLY)、丝氨酸(SER)、丙二醛(MAL)和阿司匹林(ASP)——进行了深入研究。他们利用良好调温的元动力学(well-tempered metadynamics)采样构象空间，生成包含1000个构象的训练集，并为每个构象计算了基于密度泛函理论(DFT)的IQA原子能量。随后，他们训练了五组不同的GPR模型，其区别仅在于先验均值函数m的定义（从MF1到MF5，MF5指向最高的能量状态）。所有模型在静态测试集上都表现出相近且合理的预测精度，MAE接近或优于1 kcal/mol的化学精度门槛。然而，当这些模型被部署到NVT（恒粒子数、恒体积、恒温）分子动力学模拟中时，其表现天差地别。
在高达1000K的温度下进行的稳健性测试结果令人震惊。基于MF1和MF2（均值接近或略高于训练集平均能量）的模型极其脆弱，模拟平均在1.5皮秒(ps)内就会崩溃。基于MF3（均值设为训练集最大原子能量）的模型稳定性有所提升，但最高也仅能维持约45 ps。而基于MF4和MF5（均值在最大原子能量基础上进一步增加）的模型则表现出了“近乎无限”的稳定性。特别是MF5模型，在所有四个分子、所有测试温度（300K, 500K, 800K, 1000K）下，均成功完成了预设的1000 ps（1纳秒）模拟，稳健性得分达到满分1000 ps。更进一步的挑战是，研究人员使用MF5模型，从50个不同的起始构型出发，在500K下进行了长达10纳秒的扩展模拟。所有200个模拟（4个分子 × 50次重复）全部成功，累计模拟时间达到了惊人的0.5微秒(μs)，即20亿个模拟步长。
为了探究MF5模型如此稳健的根源，研究人员设计了一系列巧妙的“压力测试”。他们构建了能量极高、键长被严重扭曲的非物理起始构型。例如，将一个丝氨酸分子的O-H键、N-H键和C=O键拉伸至1.5 Å，并将末端C-C和C-N键拉伸至2.0 Å。然后，他们分别用基于MIN（均值设为训练集最小原子能量）、MF1和MF5的模型在500K下进行模拟。结果，只有MF5模型能够成功地从这种畸变状态中“恢复”过来，在1皮秒内将键长调整至合理范围，并稳定运行1纳秒。而MIN和MF1模型则在模拟初期就因预测的原子力异常而迅速崩溃。对丙二醛的深入分析显示，MF5模型在模拟初期能够预测出方向正确的、强大的恢复力矢量，这些力指向压缩过短的化学键、拉伸过长的化学键，从而快速地将系统从高能畸变态“松弛”到一个合理的构型，这个过程类似于一个阻尼谐振子。这种预测恢复力的内在能力，是模型实现超常稳健性的关键。
此外，研究还证明，这些稳健的原子能量模型不仅能用于稳定的动力学模拟，还能准确地重现势能面上的已知稳定点（驻点）。以丙氨酸二肽(AD)为例，利用基于7000个构象训练的MF5模型进行几何结构优化，能够成功地将其不同起始构型优化到正确的C5、C7ax和C7eq构象，得到的相对能量和关键结构参数（如二面角ϕ和ψ）与高精度量子化学计算(GAUSSIAN16)的结果高度一致，同时计算成本降低了约200倍。
本研究主要应用了以下几项关键技术方法：1. 基于量子理论中的原子(QTAIM)和相互作用量子原子(IQA)的能量分解方法，从DFT计算中获得具有物理意义的原子能量作为训练目标。2. 高斯过程回归(GPR)机器学习框架，用于构建原子能量预测模型，并创新性地系统调整了先验均值函数。3. 良好调温的元动力学(WTMetaD)模拟，用于有效采样柔性分子的构象空间，生成具有结构多样性的训练数据集。4. 基于多样性感知的子采样器(DAS)从大量构象中选取代表性训练子集。5. 使用研究团队自研的DL_FFLUX程序（基于DL_POLY4）进行分子动力学模拟，以及FEREBUS引擎进行高斯过程模型的训练和超参数优化。
研究结果
原子IQA能量模型的静态性能
所有基于不同先验均值函数（MF1-MF5）训练的高斯过程原子能量模型，在包含1000个构象的测试集上均表现出合理的预测精度。分子总能量的平均绝对误差(MAE)接近或优于1 kcal/mol的化学精度。元素级别的分析显示，非氢原子（重原子）的原子能量误差通常高于氢原子，这与重原子IQA能量的变化范围更大有关。有趣的是，尽管不同均值函数模型的静态MAE差异很小，但这并未预示它们在动态模拟中会有相似的稳健性。
模拟稳定性与稳健性测试
动态模拟结果是模型稳健性的终极试金石。在300K至1000K的NVT模拟中，基于不同先验均值函数的模型表现出了数量级级别的差异。MF1和MF2模型极其脆弱，平均模拟时间难以超过1.5 ps。MF3模型稳健性有所改善，但仍有限（最高~45 ps）。而MF4和MF5模型则展现出卓越的稳定性，其中MF5模型在所有测试系统和温度下均实现了1000 ps的完美稳健性得分。提高模拟温度有助于系统探索更广阔的构象空间，但即使在1000K下，基于FFLUX（MF5模型）的模拟也主要保持在训练数据所覆盖的构象空间内，并在偶尔超出时能通过预测恢复力来维持物理合理性。
均值函数与高能结构的松弛
为理解稳健性差异的根源，研究测试了模型从严重扭曲的高能起始构型中恢复的能力。对于键长被故意扭曲的丝氨酸和阿司匹林分子，只有基于MF5的模型能够成功地将系统松弛到合理状态并完成1纳秒的稳定模拟。而基于MIN（均值最低）和MF1的模型则迅速崩溃。这表明，将先验均值函数m向高能态偏移，赋予了模型强大的外推能力和系统松弛的倾向。
模拟温度与恢复力的预测
以丙二醛为例，分析了MF5模型在从畸变态开始模拟时预测的原子力。在模拟最初的几步，模型预测的原子力明确指向修正畸变键长的方向（压缩过短的键，拉伸过长的键）。这些恢复力迅速降低了系统的势能，并使键长回归合理值。恢复过程伴随着系统总能和最大原子力的阻尼振荡式下降。模拟温度越高，初始恢复力引发的振荡幅度越大，但模型均能成功完成松弛并进入稳定模拟阶段。
扩展的MD模拟与部署成本
使用MF5模型对四个分子进行了规模更大的测试：每个分子进行50次独立的500K NVT模拟，每次长达10纳秒。总计200次模拟全部成功，累计模拟时间达0.5微秒。在计算成本方面，使用与原子数相等的CPU核心数时，每个MD步长的耗时在2-4毫秒之间，其计算效率与甚至优于部分在GPU上运行的先进MLP。
超越原子IQA能量模型的稳健性：重现势能面驻点
研究还评估了稳健的原子能量模型再现势能面稳定点的能力。使用基于7000个丙氨酸二肽构象训练的模型进行几何优化，能够将从高能畸变态出发的起始结构成功优化到正确的C5、C7ax和C7eq构象。优化得到的结构能量、关键二面角（ϕ, ψ）和原子间距离与高精度量子化学优化(GAUSSIAN16)的结果高度吻合，相对能量误差很小，并能正确识别C7eq为全局能量最低点。同时，FFLUX模型优化成功率更高（能处理GAUSSIAN无法优化的畸变结构），且计算速度快约200倍。
结论与讨论
这项研究取得了里程碑式的成果，首次报道了能够在高达1000K的温度下实现长期稳定（微秒级）分子动力学模拟的物理信息高斯过程原子能量模型。其核心结论是：高斯过程回归模型中常被忽视的先验均值函数m，是决定其稳健性的关键控制参数。 研究发现，将m向高原子能量状态偏移（如本研究中的MF5），能够产生极其稳健的模型。这是因为在模型外推至训练数据未覆盖的区域（模拟中常见）时，高斯过程的后验均值预测会收敛于其先验均值m。一个被设为高能态的m，使得模型倾向于预测更高的原子能量，这在物理上对应于预测出阻止化学键断裂（爆炸）或过度压缩（内爆）的恢复力。这种与生俱来的、防止非物理构型出现的能力，是MF5模型实现“前所未有稳健性”的根本原因。
该研究的重大意义在于多方面的突破。首先，它解决了机器学习势能函数在高温模拟中稳定性差的长期难题，将稳定模拟的时间尺度推向了微秒级，温度范围扩展至1000K，为模拟高温过程（如材料烧结、化学反应）提供了新工具。其次，它揭示了先验知识融入机器学习模型的一个新维度——先验均值函数的设计，为开发更稳健的物理信息机器学习模型提供了关键理论指导。第三，该方法基于具有严格物理意义的IQA原子能量，继承了量子拓扑原子的可迁移性，与主流Behler-Parrinello类MLP中任意、无物理意义的位点能量划分有本质区别，预测结果更具物理可靠性和可解释性。最后，该模型仅需训练原子能量即可通过解析公式获得准确的原子力，既能进行稳定的动力学模拟，又能完成精确的几何优化，且计算效率高，展现了强大的实用潜力。
这项工作为未来挑战性的体相模拟、涉及反应的过程模拟以及更高精度的电子相关能计算奠定了基础。通过将量子力学基本原理作为归纳偏置深度融入机器学习框架，该研究成功地构建了兼具高精度、高稳健性、高效率和物理可解释性的下一代原子尺度模拟模型，是计算化学和机器学习交叉领域的一项重要进展。
打赏
