关键词：反应力场训练集，键扫描，参数拟合

软件版本：AMS2024.101以上

概述：本文通过一个简单例子展示拟合ReaxFF力场的基本流程。以水分子中H–OH键解离过程的DFT扫描作为训练集，对现有力场中的一些参数进行优化、微调。简单起见，本教程中的训练集只包含了一个作业，后面的其他高阶例子，将从各种角度讲解力场训练涉及的问题。

注意：真正的训练，需要包含大量训练集，以充分描述所训练参数，而不是仅仅一两个Job。

打开AMSjobs面板，新建任务AMSinput，创建水分子（建模教程参考链接）。注意：输入的初始结构尽量合理，然后对结构进行优化。优化时可采用低水平DFT方法或者分子力场方法预先优化，以避免显然不合理的结果，然后再用较为精确的方法优化，结构优化教程参考优化分子的几何结构、以聚乙烯为例，演示低维周期性体系的结构优化、Quantum ESPRESSO结构优化、晶格常数优化。

如下图所示，任务类型选择PES Scan，电荷、自旋等信息根据分子结构选择，在这里我们做的是中性水分子，计算方法和基组根据体系选择，训练集一般不需要太高精度。

设置结构扫描：按住shift+鼠标左键分别选中O、H原子，Model → Geometry Constrains and PES scan，点击distance前面的+，即扫描O-H键长，扫描点数设置为11，扫描范围为0.85埃 ~ 1.35埃。保存作业（例如名为H2O_pesscan）并运行。

注意：在这里希望大家一定要养成检查作业设置情况和运行情况的习惯，了解每一个参数设置的原因，否则会浪费大量的计算资源和时间

更详细的灵活的势能面扫描教程参考：

• 周期体系：BAND对周期性体系进行势能面扫描(得到鞍点的初始猜测)
• 分子体系：势能面扫描PES
• 分子体系：PES多条件同步扫描

SCM → Movie：

在这里需要检查能量曲线是否合理，错误的训练集可能会导致训练集的污染。在键解离过程中如果存在自旋翻转的话，需要特别注意，只保留正确的、需要的那一段能量曲线。双击纵坐标可以修改纵坐标的单位。

在AMSjobs中选择完成的训练集任务H2O_pesscan，Jobs→Add to ParAMS（或快捷键Ctrl-T），将打开ParAMS界面和导入对话框如下所示： 在对话框中，ResultsImporter 选择 Add Single Job，Task (for new job) 选择 PESScan，将Properties 设置为 pes（这表示将这个作业的每个PES收敛点都会被读入），能量单位可选kcal/mol（根据个人习惯），然后击“Import”即可导入。这种方式导入训练集，会得到一条完整的能量曲线，是将这条能量曲线当作一个整体去作为训练集或者参考集的。如果把PES中的各个点作为single point导入，也是可以的，如此就会与其他单点训练集一视同仁对待，但这样会失去了训练能量曲线、能垒的计算特性。

导入后，可以看到相关训练集结果如下图。其中All列出各种条目，有训练集也有其他信息条目，不过我们只需要关心训练集，因此为了显示的更简略，实际上也可以切换到Trainning Set，这样就只列出下图红色框中的一条了。它包含的就是前面PES作业这个训练集（具体信息包括扫描的键、采样的点的数量等）。Detail那一列显示H2O_pesscan，relative_to=2，其中relative_to_2表示势能面扫描的能量以第3个点为参考（0表示第1个点），因为该点的能量最低（具体可以看上面的能量曲线图）。双击该行即显示下图中第二个窗口，可以看到具体能量数值，其中第三个值为0，其他点的值大于第三个点的值。

如果要将参考点改为第1个点，可以将relative_to_2改为relative_to_0，点击OK，左侧菜单栏Training set → Get Data From Ref Jobs。更新后，重新双击该训练集的各个点的能量就变成了以第一个点为参考（第1个点的能量成为0），不过能量曲线并没有发生变化，因此对拟合结果不会有影响。

还可以看到权重和sigma参数：

• 权重I：每个训练集都可以有自己的权重。可查看Error function公式（这是拟合的非常重要的公式，也叫损失函数Loss function，我们优化力场的核心就是得到这个函数的全局最小值）和ReaxFF手册，其中权重对应不同单位的训练集有所不同，如能量一般选1.0，键长一般选0.1，键角一般选5.0。代表了能量键长键角的归一化精度为1.0kcal/mol, 0.1埃和5.0度，在其核心意义将有物理意义的单位变为无量纲误差值，在这里权重设置越小，代表它的误差占比越高。
• 权重II：每个训练集的权重假设为1，如果有10个类似的训练集，那么就相当于在这方面的权重是10了。所以训练集的条目数量分布情况，也有权重的效果。例如，你的Training Set列表里面，有1000条数据，Angle有500条，其他的是Energy、Bond之类（简单假定每条的Weight都相同），那么训练的时候，就有有一半的需求都在满足角度的准确性上，对其他方面的要求就被降低了，这样是非常容易导致“过拟合”（过拟合的概念，参考ParAMS高阶、重要教程：如何选择训练集，第9节）的。
• Sigma参数：含义见下文

菜单栏Parameter → ReaxFF，将会导入一个力场模板，所有值都需要优化。

这里我们从Water2017.ff作为初始力场去优化，因此菜单栏Parameter → Load ReaxFF ForceField→Water2017.ff，如下图所示：

窗口底部显示该力场的各个参数的值，以及训练过程中需要填写的最大值和最小值，训练时力场参数将在这个区间内调整，以使得Loss function达到全局最低。

在这种情况下，我们将优化ReaxFF力场的所有O-H键参数，起始点（初始参数）为Water2017.ff（DOI:10.1021/acs.jpcb/7b02548）。

• 选择初始参数对拟合结果十分重要，虽然反应力场具有较强的物理意义，可迁移性也非常强，但是选择合适的力场参数依然十分关键，对拟合结果和分子动力模拟结果影响非常大，首先反应力场参数主要由三个大部分组成，即Combustion（主要针对燃烧爆炸分解）,Water（主要针对液相，蛋白质等），以及一些其它体系。想了解更多的可以去看DOI: 10.1038/npjcompumats.2015.11。
• 此外还要根据你的模拟体系选择相应的元素，如果有很多套力场的话，最好能够查看原始文献的训练集有哪些，测试哪个力场对于你的体系模拟结果最可靠，对于常见元素CHON体系，Combustion和Water的参数有明显差异，尤其是对于弱相互作用，因此，如果做有机液相反应，Water2017.ff是比较好的，因为在训练集中增加了密度以及径向分布函数等数据集，对于非键相互作用的描述是做了修正。

参数中我们需要优化的是OH键参数，因此这里下半窗口的Block栏选择BND（键参数，Bond缩写，类似地ATM是原子参数，ANG是键角参数，HBD是氢键参数），Category栏选择Standard，如此列表中的参数将只显示BND相关的标准参数。在Active栏，将H-O相关的参数全部选中（Shift+第一个和最后一个，然后勾选，即可全部选中）。

通常情况下，你不会优化所有的键参数，而是只优化其中的几个（你可以尝试不同的组合）。在这里，我们优化所有键参数仅用于演示。例如，优化HO键的π键参数（第二、第三个参数）并不是很有用！因为OH键只有sigma键，只有双键和三键才有π键，ReaxFF的函数设计非常巧妙，具有一定的物理意义。AMS提供了参数选择的一个工具，为了保持本文逻辑流畅性，这里先略去不讲，整个过程演示完毕后，最后一节中介绍了一个辅助工具，充分理解后，可以灵活使用该工具。

• 在这里有需要深入学习的，建议了解所有的参数意义和参数类型分类，查看ReaxFF手册，我们需要针对性的对于哪些参数进行优化，了解参数意义才能更好地根据自己体系的内容拟合特定组合的参数，例如Na+和Na原子参数如何区分，共轭体系和超价分子SF6需要拟合哪些参数，因为目前ReaxFF(除了eReaxFF)没有电子的概念，因此需要解决上述问题需要对ReaxFF函数形式和参数有更深入的了解，建议阅读原始手册、讲义！

对于ReaxFF，我们建议使用CMA-ES优化力场。这是一种在中心点（平均值）周围对Popsize参数进行采样的方法,第一个中心点是初始参数。Sigma参数影响采样分布的宽度（采样时距离中心点参数有多远）。在由Popsize参数评估组成的迭代之后，更新中心点和sigma参数。中心点将越来越接近最优化的力场参数，Sigma参数将随着优化的进行而降低。这样，采样分布在优化开始时非常宽，但随着时间的推移会变小。sigma参数达到阈值时，优化停止。初始设置很小的Sigma参数可以加快收敛速度，更快得到一个局部最优结果，但坏处是优化的全局性减弱，即，有可能有更好的力场在其他参数区域。

此外，优化方法除了最初的单参数扫描法（该方法严重依赖初始参数好坏，且忽略参数之间耦合作用，但是对于有经验的来讲不错）还有蒙特卡洛-模拟退火方法，有兴趣可以在英文教程中阅读相关文献内容。修改优化方法如下图，在Option → Optimizer中设置。典型参数设置如下：

其中Max loss function calls控制运行总的运行长度，可以酌情增大，优化结果如何，重点关注下文所示的几个曲线。

File → Save as保存文件为reaxff_water.params并运行。

运行过程中，在ParAMS右下角Graph中可以实时显示优化的情况如何。如下图所示：

分别为

• 图1：Error function误差值，误差值如果小于1kcal/mol就比较好了，当然也可以继续优化
• 图2：平均绝对误差（MAE）类似Error function
• 图3：势能面DFT训练集（即reference）的最优拟合结果情况（即所谓prediction）点与线越靠近，表示质量越好。

如果图3，PES改为H2O_pesscan_relative_to_2，即将ReaxFF力场的结果与之前DFT势能面扫描的结果进行对比，同时Data from确保三项都勾选，则如图4所示一样：显示训练前后、原始DFT训练集的对比，其中

• Traning:ΔE_prediction即训练后的结果
• Traning:ΔE_reference即DFT的结果
• Traning:(Initial)ΔE_prediction即初始力场的结果

注意此处第四图左上角是Best Training，指Traning:ΔE_prediction是训练后的最优结果，如果改为Latest Training，则Traning:ΔE_prediction是生成的最后一个力场扫描势能面的结果，一般关心Best Training即可。该势能曲线的数据位于*result/optimization/training_set_results/best/pes_predictions/h2o_pesscan_relative_to_2.txt。

优化好的ReaxFF力场位于：*result/optimization/training_set_results/best/ffield.ff，拷贝出去即可使用。或在ParAMS窗口 → File→ Open Optimized Engine in AMSinput → Best training，即可打开一个新的AMSinput窗口，已经设置了力场为该力场：

力场文件的第一行是注释，可以用一些有意义的内容来描述是谁训练的，以及训练数据的类型（体系、参考方法） 在ParAMS 界面中也可以直接更改标题Parameters→Edit ReaxFFHeader。修改完毕即可如上所说使用该力场。

建议对于新参数对应的分子进行优化和MD测试。如果模拟过程中出现的意外、不合理的反应过程，需要据此重新思考，优化训练集，进一步调整参数，不断迭代训练得到合理可靠的力场参数，这个过程十分的繁杂，希望大家保持良好心态和科学严谨的态度!

AMS2024还有一个敏感度测试的功能，实际上在帮助用户选择哪些参数需要优化方面，非常有用。该功能用来测试力场里面的各个参数对Error function的敏感情况，得到敏感度排序后，我们原则上只需要优化敏感参数，不优化不敏感参数。从而能够减少参数数目，达到减少计算量（计算量随参数数目呈指数增长）、更快收敛的效果。

在进行力场优化前，将Optimization改为Sensitivity，勾选Run sampling，Repeat calculation n times设置为5，Number of samples per repeat 选择500。保存job，提交运行任务。

如下图所示从结果我们可以看到不同的敏感度值，其中前三位的参数敏感度大于0.03,一般我们选择大于0.03参数进行优化会更加节约计算量。

Graph中有饼图反应损失函数等的分布情况。

如果训练集包含的内容与勾选的要优化的力场参数不匹配，例如参数里面只勾选了O-H键的参数，但是生成了键角、能量、电荷训练集，那么计算Sensitivity的时候就会报错。

感谢苏州大学刘越老师整理资料！