能量最低反应路径

前提条件：绘制简单分子的可视化建模工具（Sketching simple molecules Visualizer）、绘制卟啉分子的可视化建模工具（Sketching a porphyrin Visualizer）教程

背景

在化学反应建模中，需要一种有效的方法识别和验证反应过渡态，这是理论化学中普遍存在的问题。一条化学反应路径通常有多个步骤，需要正确识别反应模拟路径中的中间态。

介绍

本教程介绍FlexTS模块的主要工作模式，使用系谱聚类法研究分子开关萘酞菁的互变异构反应势垒(Liljeroth, Repp, Meyer, 2007)。首先，使用Materials Studio中的DFTB+模块搜索整个反应路径，该模块表明互变异构反应中存在两个独立等效步骤。然后使用DMol3模块精确确定其中一个势垒的过渡态结构，并得到整个开关的能量分布。

本教程包括如下部分：

开始

建立计算所需的分子构型

利用DFTB+计算反应路径并对分析计算结果

利用DMol3模块进行过渡态优化

注意：为了和本教程中的参数保持一致，可以使用Settings Organizer对话框将项目中所有参数都设置为BIOVIA的默认值。

1、开始

首先启动Materials Studio并创建一个新项目。

打开New Project对话框，输入项目名MEP，单击OK按钮。

新项目将以MEP为项目名显示于Project Explorer中。

2、建立计算所需的分子构型

在本教程的这一部分中，将在两个不同的3D原子结构文档中构建反应物和生成物的结构。第一步是打开一个新的3D原子结构文档并绘制反应物萘酞菁的结构。

打开一个新的3D原子结构文档，将其命名为reactant_start.xsd。从卟啉核心开始，绘制出上图所示的反应物分子，中心H原子分列左右两侧。单击Clean

按钮并保存文档。确认该分子有86个原子，其化学式为C52H30N4。

提示：在绘制碳六元环时按住ALT键可以绘制苯环。

打开一个新的3D原子结构文档，将其命名为reactant_start.xsd。从卟啉核心开始，绘制出上图所示的反应物分子，中心H原子分列左右两侧。单击Clean

按钮并保存文档。确认该分子有86个原子，其化学式为C52H30N4。

提示：在绘制碳六元环时按住ALT键可以绘制苯环。

需通过运行DFTB+模块的结构能量最小化计算，确保反应物的起始几何结构合理。

打开DFTB+ Calculation对话框，将任务Task更改为几何优化GeometryOptimization，然后单击Run按钮。优化完成后，旋转生成的分子以确保所有原子处于同一平面上。

将生成的reactant_start.xsd文件复制到项目根目录下，并将其重命名为reactant.xsd。

接下来，通过修改反应物建立生成物模型，这样做的目的是便于在建立过渡态搜索路径时，反应物和生成物的原子可以一一匹配，从而创建合理的反应物到生成物的轨迹文件。

复制reactant.xsd文件的副本，并将新文件命名为product.xsd。在Atoms & Bonds工具栏中，单击Bonds按钮后的下拉菜单箭头，选择Monitor Bonding

将中心处的H原子移动到另一组相对位置（顶部和底部），如上图所示。选中两个位于中心的H原子，单击Clean按钮

并保存文档。

提示：为了使反应路径计算过程的可视化效果最佳，请保持实时检测成键功能打开。

下一步是创建一个初始轨迹，使反应物和生成物中所有原子正确匹配。

打开reactant.xsd和product.xsd两个文档，并关闭所有其他窗口。在菜单栏中选择Window | Tile Vertically，使得两个结构并列显示。

打开Tools | Reaction Preview对话框，将反应物Reactant和生成物Product设置为打开的两种结构文件。单击Match...按钮打开Find Equivalent Atoms对话框。

下一步是确保反应路径的化学意义的关键。然而，由于生成物的结构是对反应物的结构修改后创建的，两种结构中所有原子的顺序都是相同的。因此可以直接实现反应物和生成物原子之间的正确匹配，并创建从反应物指向生成物的初始猜测轨迹。

在Find Equivalent Atoms对话框中，选择reactant.xsd文件中第一个未匹配原子和product.xsd文件中相应的第一个未匹配原子。单击Set Match按钮。等效原子分析会自动选择下两个未匹配的原子。一直单击Set Match按钮，直到反应物和生成物中所有原子都完全匹配为止。关闭Find Equivalent Atoms对话框。

在Reaction Preview对话框中，选中叠加结构Superimpose structures复选框，然后单击Preview按钮。

将新创建的轨迹文件重命名为naphthalocyanine.xtd，并选择Monitor Bonding选项

保存轨迹文件并关闭Reaction Preview对话框。

此时需要将轨迹文件重复播放几次，以确保所有原子已经正确匹配。如果不进行该操作，有时可能会观测原子在空间中不规则运动，与反应物或生成物没有任何联系。

提示：如果在播放轨迹时观测到不正常的原子运动，请返回Reaction Preview | Match...，确认已经在反应物和生成物中正确选择了等效原子进行匹配。

3、利用DFTB+计算反应路径并对分析计算结果

本教程的下一步骤是利用DFTB+模块计算分子开关的能量分布。

可以通过最小能量路径Minimum Energy Path任务，在DFTB+和DMol3计算对话框中使用FlexTS模块。FlexTS模块主要依赖于应力高度精确的计算结果，因此必须首先减小DFTB+的自洽收敛阈值，并增加迭代次数。

打开DFTB+ Calculation对话框，选择Electronic选项卡，然后单击More...按钮打开DFTB+ Electronic Options对话框。将SCC收敛容限SCC tolerance更改为1e-10，最大SCC循环次数Max. SCC cycles更改为100。关闭对话框。

在DFTB+ Calculation对话框中，选择Setup选项卡并将任务Task更改为最小能量路径Minimum Energy Path。点击More...按钮打开DFTB+ Minimum Energy Path对话框。

Setup选项卡包含FlexTS运行的基本信息，包括三种不同的模式：

全路径（Full Path）包括五个阶段：

1. 反应物和生成物结构的几何优化（可选）

2. 轻推弹性带计算，以估计过渡态（TS）的位置

3. 使用完全基于能量和应力的混合特征向量追踪方法优化过渡态(Kumeda et al. 2001)

4. 在过渡向量的正方向和负方向上移动过渡态，然后进行优化，以确定对应于该反应步骤的局部极小值

5. 尝试使用Dijkstra算法通过独立步骤连接反应物和生成物状态

该过程描述了一个单路径循环(Carr et al, 2005)。在本教程的DFTB+部分中，将使用两个连续的路径循环来研究萘酞菁分子开关。

轻推弹性带（Nudged Elastic Band）在反应的初始和最终状态之间进行插值，并使用双轻推弹性带方法(Trygubenko and Wales, 2004)优化整个反应路径，无需进一步优化过渡态。可用该方法研究势垒较小的反应，获得能量变化的粗略初步概况，或研究无势垒反应。

TS路径（TS Path）从预估的单一过渡态结构开始，仅运行全路径过程的阶段（3）和（4），以优化过渡态以及相应的反应物、生成物结构。在本教程的DMol3部分中使用此模式可以更准确地估计反应势垒。

Setup选项卡中的其他选项：

路径循环次数

是否优化反应物和生成物

确定局部极小值或过渡态时是否进行几何构型对比

FlexTS计算中轻推弹性带方法的基本参数

通常，对于全路径计算，选择自动确定NEB参数Determine the NEB parameters automatically。将基于预定义的构型密度和最大构型帧数生成NEB参数。对于某些情况，或者对于专用的轻推弹性带计算，也可以手动设置参数。

弹簧常数Spring constant控制弹性带的强度。在某些情况下，可以应用此选项通过允许体系更快地弛豫到实际最小能量路径，来加快计算速度。但是，确切的弹簧常数值在某种程度上取决于体系结构。

设置运行模式Run mode为全路径Full Path。设置路径循环次数Number of pathcycles为2。

选择DFTB+ Minimum Energy Path对话框的Advanced选项卡。

此选项卡上的大多数选项都是FlexTS模块运行时具体性质的阈值。

关闭DFTB+ Minimum Energy Path对话框。在DFTB+ Calculation对话框中，选择Job Control选项卡。选择合适的网关位置Gateway location和运行计算使用的CPU核数number of cores，单击Run按钮运行计算。关闭DFTB+ Calculation对话框。

此计算任务需要一些时间才能运行完成。计算完成后，可以分析计算结果。

在Project Explorer中，展开naphthalocyanine DFTB+ Minimum Energy Path文件夹并查看其中包含的文件。

该教程案例中最重要的图表是naphthalocyanine Connected Path.xcd。该曲线图是FlexTS模块成功确定过渡路径中连接反应物和生成物所需所有步骤后得到的完整路径的能量分布。

从曲线图中可以看到，该路径有两个单独的反应步骤，以及存在能量最小值的中间态。可以打开相应的轨迹文件naphthalocyanine Connected Path.xtd以查看路径。分子内核中的氢原子不是以协调一致的方式跳跃运动，而是一步一步地移动。这种运动是对称的，因为两个单独的势垒具有相同的形状（但方向相反），并且具有匹配的反应物、生成物和能量。可以使用这种对称性来简化后续步骤。

注意：路径长度可能取决于初始结构和最终结构的构型，并且结果可能与上图有所不同。

从曲线图中可以看到，该路径有两个单独的反应步骤，以及存在能量最小值的中间态。可以打开相应的轨迹文件naphthalocyanine Connected Path.xtd以查看路径。分子内核中的氢原子不是以协调一致的方式跳跃运动，而是一步一步地移动。这种运动是对称的，因为两个单独的势垒具有相同的形状（但方向相反），并且具有匹配的反应物、生成物和能量。可以使用这种对称性来简化后续步骤。

注意：路径长度可能取决于初始结构和最终结构的构型，并且结果可能与上图有所不同。

提示：在查看轨迹和其他结果时，需为每个文档启用实时检测价键功能Monitor Bonding

naphthalocyanine Results.std数据表包含确定过渡态所需的所有数据和结构。本数据表的All Segments表单的每行为FlexTS模块计算出的一个势垒。Connected Path表单为连接反应物和生成物的实际路径中涉及的各个步骤。在该例中的计算中，两个表单包含相同的结果，但顺序可能不同。

结果数据表还包含一个Computational Settings列，其中以脚本格式包含对话框上的所有非默认设置，因此可以重新创建计算。当在单个数据表中存储多个最小能量路径Minimum Energy Path计算结果时，此列将尤为重要，因为它分别存储每个计算得到的势垒所用的不同计算模块、电子结构设置、电荷和自旋设置等。

打开Connected Path表单，可以看到以kcal/mol为单位列出的正向和反向势垒。每行以一个收集文档开始，该文档包含根据相对能量排序的反应物、过渡态和生成物结构。可以使用这些文档进行进一步分析。例如作为精确计算势垒的输入文件，或作为DMol3或CASTEP中动力学任务的输入文件。

确定DFTB+计算的正向和反向势垒。

反应物或生成物与过渡态TS之间的能量差约为22 kcal/mol。两种过渡态TS与中间态之间的能量差约为15.5 kcal/mol。

naphthalocyanine MEP Report.txt文档为FlexTS模块运行的文本输出。查看本文档时，可以识别单个路径循环的各个部分，并发现在本例中，计算需要两个路径循环才能完成。

如果需要，可以查看曲线图连贯性以帮助解决计算问题。首先，NEB计算必须成功收敛到指定的阈值。可以通过在文档中搜索NEB部分标题中的术语每帧构型的能量Energy per image对其加以验证。

成功完成NEB计算后，FlexTS将自动优化候选的过渡态构型。可以通过搜索优化循环起始Beginning of optimization cycle来找到此部分。每个循环包括两个步骤：首先收敛Hessian矩阵的特征值最小的特征向量的方向（例如，过渡态的上坡向量），然后跟随该向量直到它不再相关。另外要验证每个新方向的收敛特征值，必须为负值。

最后，每次FlexTS计算都会返回许多轨迹和链接图表。结果数据表的AllSegments选项卡的每一行对应一条轨迹和一个图表，若确定了连接的路径，则会返回另一个图表和轨迹。轨迹对应于过渡态下的两个下坡几何优化，取自全路径任务的第4阶段。

4、利用DMol3模块进行过渡态优化

打开数据表naphthalocyanine DFTB+ Minimum Energy Path\naphthalocyanine Results.std，双击Connected Path表单第一行中的结构文档。

该收集文件包含所有可作为DMol3最小能量路径计算起点的相关信息，可以根据指定的反应物和生成物的物理体系运行完整的路径循环，在该特定体系中重现整个反应步骤。在这种情况下，应使用DFTB+计算中找到的TS作为TS路径计算的初始猜测构型。

打开DMol3 Calculation对话框，将任务Task更改为最小能量路径MinimumEnergy Path，计算精度Quality更改为Fine。选择密度泛函，将Theory Level更改为GGA，Non Local Functional更改为PBE。

在Electronic选项卡上，单击More...按钮打开DMol3 Electronic Options对话框。将SCF收敛容限SCF tolerance降低至1e-8，并将最大SCF循环次数Max. SCFcycles增加至100。关闭DMol3 Electronic Options对话框，然后选择DMol3 Calculation对话框的Setup选项卡。

点击More...打开DMol3 Minimum Energy Path对话框。

DMol3最小能量路径计算任务的原理与DFTB+中的名称一致的任务相同。主要区别在于，输入数据的能量单位为Ha。单位将转换为kcal/mol，即FlexTS模块内部和Materials Studio通常使用的单位。

在本例中，使用DFTB+中找到的TS对DMol3过渡态搜索进行初始化。要实现这一目标，只有一种选择。

在DMol3 Minimum Energy Path对话框中，将运行模式Run mode更改为TS路径TS Path，然后关闭该对话框。

在DMol3 Calculation对话框的Job Control选项卡上，选择合适的网关和多个处理器，然后单击Run按钮。

这个计算需要一些时间才能完成。

DMol3最小能量路径Minimum Energy Path任务的工作方式与DMol3中的其他任务不同。该任务需要运行FlexTS模块，FlexTS又需要与DMol3计算引擎协同工作。这意味着没有可以直接查看或修改的输入文件。相反，DMol3输入文件是在FlexTS需要时动态生成的。

计算完成后，下载并查看结果。可以在初始DFTB+计算中找到与单个势垒之一类似的结果。现在，在反应物或生成物和过渡态TS之间有一个更精确的能量势垒，约为13 kcal/mol，在中间最小值和过渡态之间约为6 kcal/mol。数据表文档的每一行还包含用于此计算的参数设置可供以后参考。

最终结果数据表中的收集文件可用于进一步处理，例如在DMol3动力学计算中，为Cantera中的反应机理生成输入数据。