DRUGAI
今天为大家介绍的是来自美国的Teresa Head-Gordon团队与Samuel M. Blau团队发表的一篇论文。识别过渡态——反应物和产物极小值之间的势能面鞍点——对于预测动力学障碍和理解化学反应机制至关重要。在本研究中，作者训练了一个完全可微分的等变神经网络势能，NewtonNet，针对数千个有机反应，并推导出分析性海森矩阵。通过将计算成本相对于密度泛函理论（DFT）ab initio源降低几个数量级，作者能够在每一步的鞍点优化中使用学习到的海森矩阵。作者展示了完全机器学习（ML）海森矩阵在240个未见过的有机反应中稳健地找到过渡态，即使初始猜测结构的质量降低，同时将收敛所需的优化步骤数相比于准牛顿DFT和ML方法减少了2-3倍。所有数据生成、NewtonNet模型和ML过渡态寻找方法均可在自动化工作流程中获取。
计算识别过渡态（TS）在量子力学势能面（PES）上对于预测反应障碍和理解化学反应性至关重要。过渡态是一级鞍点，定位它们在高维PES上尤其具有挑战性。使用梯度信息的二阶方法（如海森矩阵）可以更稳健地识别局部极小值，但计算海森矩阵的代价高昂，几乎所有的过渡态优化方法都依赖于近似海森矩阵，通常称为准牛顿（QN）方法。
近年来，深度学习模型为获取和应用海森矩阵提供了新思路。作者对等变消息传递神经网络进行微调得到了NewtonNet，训练数据包含约1000万个通过NEB方法生成的气相有机反应配置。尽管训练数据中没有海森矩阵，但可以通过反向传播推断海森矩阵。作者将从NewtonNet模型获得的海森矩阵应用于240个独立的有机反应的过渡态优化中。结果表明，与近似海森矩阵方法相比，搜索步骤减少了2-3倍，并且提升了效率和稳健性。作者的方法计算海森矩阵的速度比对应的ab initio计算快1000倍，并在寻找过渡态方面比使用ML或DFT PES的QN方法更具稳健性。这种结合更高效的计算、更少依赖于良好初始猜测的需求以及对未见反应的强大收敛性，拓展了使用完整海森矩阵进行过渡态优化的机会。
学习预测DFT海森矩阵
图 1
图1a展示了NewtonNet模型，通过转化和聚合原子特征ai，预测DFT计算的分子能量E，这些特征从空间邻居aj和原子间距离Rij中获取局部化学环境信息。分子能量E相对于原子位置Ri进行微分，以预测原子力Fi或梯度gi，并且可以进行二次自动微分以获得海森矩阵Hij。作者展示了在小有机分子和甲烷、氢燃烧等多种化学体系中，其能量和力的预测精度优异。
作者的机器学习模型在ANI-1数据集上进行了预训练，包含超过2000万个小有机分子的非平衡构象。图1b表明，原始ANI-1数据集主要由接近平衡的几何构型组成，反应路径在反应的亚稳态附近显著不足。因此，预训练的机器学习模型在反应物和产物状态下能准确预测能量和力，但在过渡态附近表现较差（图1c, d）。幸运的是，T1x数据集是一组包含9644740个通过NEB方法生成的分子构象的基准数据，它涵盖10073个有机反应。这组数据更好地代表了整个反应路径，这使得能够微调预训练模型。微调后的模型在过渡态附近能量和力的预测精度提高了一个数量级。
由于海森矩阵计算和存储成本高，训练数据集中不包含ab initio海森矩阵的参考样本。尽管缺乏这样的训练实例，机器学习模型能平滑地预测原子力，强烈表明可以在不进行显式训练的情况下实现海森矩阵预测。基于这一假设，可以使用模型预测的梯度，通过在每个笛卡尔坐标轴上进行有限差分来实现海森矩阵的估计。然而，分析梯度更具成本效益，只要神经网络至少是二次可微的。NewtonNet模型采用了Sigmoid线性单元（SiLU）和多项式截断函数，因此是C2连续的。利用神经网络中的自动微分，可以通过反向传播解析地获得力和海森矩阵。
图 2
借助于微调后深度学习势能面的光滑性，图2显示出与预训练模型相比，分子过渡态的海森矩阵预测与参考DFT模型的准确性较高。海森矩阵预测在负特征值和正特征值区域均量化准确，未见测试反应的特征值均方根误差（RMSE）为318 kcal mol−1 Å−2，特征向量余弦相似度（MCS）为0.828。在微调后，这些值显著改善。主要误差来自正特征值空间，相对于DFT，海森矩阵特征值的低估常为20%。作者通过从ANI-1x数据集中选择1232469个化学式相同但化学键更为压缩的分子来增强T1x数据集，从而进一步提高海森矩阵预测的准确性，特征值RMSE降至267 kcal mol−1 Å−2，特征向量MCS升至0.839（图2c）。这一改进不仅归因于数据量的增加，而是通过建立更平衡的数据集有效减轻了正特征值的低估。
利用NewtonNet的过渡态优化
经过微调的NewtonNet模型用于预测过渡态（TS）特性，随后被应用于涉及新反应的实际TS优化场景，这些反应独立于扩展的ANI-1x/T1x训练和测试数据。这些反应包括氢迁移反应、内环化和外环化、广义Korcek步骤2反应、逆烯反应以及反向1,2和1,3插入。作者将微调的NewtonNet模型与Sella连接，其中Sella是一种先进的开源TS几何优化器。在Sella中，笛卡尔坐标与冗余内部坐标之间的转换是自动处理的，并且海森矩阵被迭代对角化，以便在几何鞍点优化中使用最左侧的特征向量。为了开始对240个Sella基准反应的TS优化，作者使用KinBot生成初始猜测，利用反应模板，其中每个模板还定义了给定反应的预期反应物和产物末态。
图 3
图3显示了在每次优化步骤中提供完整显式海森矩阵显著提高了优化效率，这在相对于DFT时是可以承受的。在图3a中，作者发现收敛到一个TS所需的步骤数量可以减少至QN方法的2倍。这一趋势显著非线性，对于需要更多优化步骤的困难任务，完整海森矩阵优化更具优势。如果考虑QN近似失效时初始海森矩阵构建和海森矩阵重构的迭代对角化步骤，当考虑这些梯度调用时，所需步骤减少接近3倍。更高的优化效率可以源于两个因素：每个优化步骤的置信度增加以增大步长，以及更优的整体优化路径。如图3b所示，这两个因素都提高了使用完整解析海森矩阵的TS优化效率。
图 4
接下来，作者通过比较KinBot的预期反应与TS优化后IRC预测的反应，考虑使用完整解析海森矩阵与近似海森矩阵（使用ML或DFT）的TS优化鲁棒性。如图4a所示，NewtonNet的完整或QN海森矩阵在实现2端匹配（意图的反应物和产物末点）时的表现优于DFT的QN海森矩阵。在某些情况下，仅找到1端匹配，因为预测的产物在中性和闭壳层势能面上比KinBot预期的产物更稳定且更具化学合理性。作者还对图4a中发现的TS类型进行了表征，其中反应物和产物的图形在化学反应中是不等价的，而同构末点则表明构象TS；无论如何，两者都收敛到一阶鞍点。
生成的初始猜测结构的不可靠性可能导致较差的收敛或不准确的预测。因此，作者考虑一种鲁棒性度量，即TS优化必须能够从不良起始结构中恢复，作者通过系统地向KinBot生成的猜测结构引入随机结构扰动进行分析。图4b演示了NewtonNet的完整海森矩阵优化的鲁棒性，在噪声水平增加至2-3 pm时保持一致的性能，甚至显示出轻微的改善。相比之下，使用近似QN海森矩阵的性能迅速下降，即使使用DFT，突显出准确海森矩阵在鲁棒TS优化中的重要性。
鲁棒性的另一个重要指标是NewtonNet预测的TS结构是否在DFT势能面上具有负特征值和/或改善了DFT势能面上的TS优化结果。在图4c中，作者比较了使用完整海森矩阵和QN海森矩阵优化的NewtonNet TS鞍点结构的振动频率，并使用这些结构作为输入计算DFT海森矩阵的频率。然后，作者识别与NewtonNet的负频率模式相对应的DFT频率模式。在约96%的反应中，无论是完整海森矩阵还是QN海森矩阵，NewtonNet都预测出高精度的虚频。作者进一步使用完整海森矩阵ML TS结构作为DFT表面的初始猜测重新优化了所有240个反应（图4d）。与使用原始KinBot猜测结构开始的优化结果相比，作者观察到在DFT重新优化下，2端匹配和化学转化都有所增加。因此，从ML优化的TS结构开始，作者在DFT表面上发现整体改进的解决方案。
讨论
作者提出了一种高度通用的机器学习方法，可以基于能量和梯度数据预测ab initio海森矩阵，且仅需满足二阶可微性。尽管先前研究发现缺乏海森矩阵训练数据使得预测效果差，但研究表明，经过良好训练的机器学习模型能够有效、准确地预测反应系统的海森矩阵。作者的方法依赖于高质量的势能面和深度等变消息传递神经网络的数学关系。通过使用NewtonNet，海森矩阵的计算速度比DFT快三个数量级，且在过渡态优化中所需步骤比准牛顿方法少2-3倍。此方法可推广至多个化学和材料科学领域，未来可结合更高理论水平的数据集，实现更精确的预测。
编译 | 于洲
审稿 | 王梓旭
参考资料
Yuan E C Y, Kumar A, Guan X, et al. Analytical ab initio hessian from a deep learning potential for transition state optimization[J]. Nature Communications, 2024, 15(1): 8865.
