Targeted MD：强制反应

Targeted MD使用的Reaction Boost原理很简单：在分子动力学过程中，体系以少量Step逐渐推向某个已知的最终构型。这种“转向”，可以使用holonomic constraints（标准的Targeted MD，TMD），或者restraints，一些作者也称之为steered MD。AMS使用两种Restraints：“原子对”Pair约束、（质量加权）RMSD。RMSD约束适用于非常大的体系或Region，而Pair约束适用于小、中等尺寸的体系或Region。

这两种方法都允许用户选择一个Region，Restraints将限制在这个Region内，从而“转向”也就限定发生在这个Region内。Region是AMS用户惯用的原子集合，如果不熟悉的话，创建Region的方法参考：如何创建分区

效果展示：

H₂O + propene → 1-propanol → H₂O + propene → 2-propanol

此功能可用于：

• 迅速得到化学反应过程的采样，便于使用Active Learning训练机器学习势
• 通过分子动力学，快速地对过渡态得到合理的初始猜测（不保证是能量最低的路径）
• 适用于ADF、BAND、Quantum ESPRESSO、DFTB、ReaxFF、机器学习势的分子动力学

版本要求：AMS2024.101以上。本文使用“Pair”约束进行演示，该方法会先确定出从反应物到产物需要断裂和形成哪些键，然后对这些键施加移动约束，以逐渐断裂和形成这些键，因此这种分子动力学模拟可以翻越（可能很大的）反应能垒。

如果施加的约束足够温和，那么从反应物到产物的路径，更可能接近最小能量路径（不能保证）。要更接近最小能量路径，应使用专门的过渡态搜索方法例如NEB等。

• Bond Boost：如果在模拟过程中识别出特定的局部条件，则添加约束。使用这种方法时，用户不知道最终产物会是什么样子，目标是获得最终结构。
• 本文将介绍的Reaction Boost：约束是通过输入结构来定义。使用这种方法时，用户规定了最终产物的是什么样子，目标是看看中间会经历什么过程、什么中间结构。

该方法支持AMS的一切MD，本文以机器学习势AIMNet2-wB97MD3为例，以便观察轨迹过程中，计算能量的不确定性。该功能支持周期性结构，但本例以非周期性结构为例。

首先创建H₂O + propene模型，然后Ctrl a全选，Ctrl c复制，然后菜单栏Edit → New Molecule则创建了一个新的窗口。

注意：底部可以点击栏目切换到各个窗口，也可以点击下三角按钮，选择Multi View 2×2 Synced，这样多个窗口都合并到一个窗口中显示，鼠标可以选择不同窗口区域，然后编辑该区域的分子结构。

将底部新窗口的名字Mol-2可以改为State2，依次类推创建State3、State4窗口。将刚才复制的Mol-1的分子结构在各个新的窗口中粘贴，然后通过手动移动原子的位置，然后重新猜测键级（菜单栏Bonds → Guess bonds，当然不猜测键级也无所谓，计算过程并不读取键级信息）分别将其修改为如下结构：

• Mol-1: H₂O + propene
• state2: 1-propanol
• state3: H₂O + propene (identical to Mol-1)
• state4: 2-propanol

看起来是这样：

这样的话，整个体系在MD过程中，将首先转化为1-propanol，然后再转化回H₂O + propene，最后转化为2-propanol。在本教程中，我们在一次MD模拟中执行了3次转换。如果其中一个转换不成功，就会导致问题，因此通常对于Reaction Boost模拟，只进行一次转换（只有两个输入结构而不是四个）更加可靠。

可以在底部，单击下三角箭头，选择Single view，当然这并不影响模拟，只是个人习惯。

先按照常规分子动力学设置参数即可，具体可以参考：分子动力学：ReaxFF & 机器学习势。

注意：与常规MD不同，这里要求恒温器、恒压器的damping constant需要设置的非常小，例如5fs。

Model → MD → Reaction Boost，Reaction Boost有许多输入选项，以下是对最重要的选项的解释：

当Restraint type设置为时Pair，意味着应用约束来建立或断开单个键。这些限制作用于分子的内部坐标，输入结构如何相对于彼此旋转并不重要。输入结构各个原子间的距离（键长）很重要，因为将用于匹配目标结构中的键长。因此，建议使用与 MD 期间相同的计算引擎，进行Pre-Optimize是一个不错的方法（如果用ADF、Quantum ESPRESSO、BAND等，就没有Pre-Optimize按钮，单独优化后导入，就需要注意将原子的坐标顺序弄到和Mol-1一致，Model → Coordinates 中，选中原子可以用↑↓按钮调整顺序）。

对于约束类型 RMSD、Region的用法，请参阅其他Reaction Boost教程。

• 您可以通过单击+按钮和-按钮按钮，来添加或删除Target system。
• MD 期间访问Target system的顺序由上图面板中的顺序决定，而不是左侧分子模型区域中输入体系的顺序。
• MD 的初始结构未列在Target system中，其他所有体系均需要列入Target system列表。
• 这里的Target system是state2、state3、state4，初始结构是Mol-1。
• 这里只有2个分子，所以全局添加约束即可，没有使用Region功能。

详细说明，请参考手册：

1. Move type设置为LogForce，表示力将逐渐增加，力的大小取决于 MD 步长以及当前结构与Target system的距离。
2. Initial fraction设置为0.1，表示log force在每次的第1个Step为最终（最大）值的 10％。
3. Steps per target：指定到达每个Target system要执行多少个 MD Step，这里将其设置为1000。由于有 3 个Target system，Reaction Boost 总共花费 3000 MD Step。常规 MD 设置里面，步数大于3000，后续多余的Step，恢复正常MD，没有任何约束。
4. Reaction boost约束类型与设置：
   • Non-bonded restraints（可选）：通过设置Exponential，在初始结构中未成键的所有原子和Target system中未成键的所有原子之间添加排斥势，确保在 MD 期间不会形成不需要的新键。Epsilon和Sigma是关于排斥力的大小和衰减指数。
   • Bonded restraints（可选）：通过设置Harmonic，使用谐波势能使所有在初始结构中，确保在模拟过程中不应断裂的键不断裂。
   • Bond breaking restraints（如果至少 1 个键断裂，建议使用）：通过设置Erf 键断裂力，将键断裂力添加到应该断裂的键上，促成其断裂。这里我们施加的最大力是 0.1hartree/bohr$^2$。
   • Bond making restraints（如果至少形成 1 个键，建议使用）：通过设置 Erf 键形成力，将键形成力添加到应该形成的键上。在这种情况下，这里我们施加的最大力是 0.4hartree/bohr$^2$。

如果您预期的转变不成功，您可以考虑增加力常数、最大力或步数。保存并运行作业，如果您使用 ReaxFF 引擎，则只需几秒钟即可完成。

SCM → AMSmovie 中打开轨迹以可视化结果。它看起来应该像这样：

AMS2024 中的 AIMNet2 可以计算在 MD 模拟期间的能量不确定性，并可视化，前提是已经安装 AIMNet2。

请注意，Movie 显示的能量，包括约束能量和动能。通常，查看“Engin Energy”更有用，这是不同的计算引擎本身计算出的势能（不含约束能量和动能）。

点击能量曲线，del键删除，Graph → Engine Energy，在此图中，您可以看到计算出的能量不确定性是线周围的红色阴影区域。如果放大一点，会更容易看到：

您还可以在它自己的图表中绘制不确定性：

• Graph → Add graph
• Graph → Other… → Engine Energy U

这里，你可以清楚地看到，当分子分离时，能量不确定性更高。在转变为 1-propanol 的过程中获得的不确定性最高。