过渡态搜索、自由能、活化能的计算

本教程软件版本AMS2020.101，以F₂与CH₄的反应，演示非周期性体系，也就是分子、团簇体系的过渡态搜索。

ADF搜索过渡态一般分为如下几个步骤：

优化反应物结构（可酌情跳过，或者用MOPAC简单优化得到大致结构即可，不需要很精确）；
猜测反应路径（使用PES Scan功能，让反应物通过某一个或两个参数的线性变化，到达产物端）；
计算PES Scan得到的最高点（即近似的过渡态或者称为“鞍点”）的频率；
使用Transition State（以后简称为TSS）功能搜索真正的精确的过渡态；
计算TSS得到的过渡态的频率，校验是不是有且仅有一个虚频，并且该虚频对应的振动模式分别向反应物、产物方向振动，同时得到Gibbs自由能；
从过渡态的虚频沿着正负方向振动一帧（或三帧）分别得到倾向反应物和产物的结构，分别优化之，分别得到反应物和产物的稳定结构，并与一开始设想的反应物、产物结构对比，从而确定是单步反应还是多步反应；
计算精确优化的反应物、产物的Gibbs自由能，得到活化能

第一步，优化反应物结构

反应物分子，各自进行大致的优化，参考优化分子的几何结构。这一步，不需要优化的很精确，因为接下来的计算，都不需要很高精度，只需要提高Transit State那一个步骤的精度，即可以完全弥补前面所有精度不足的问题。

第二步，猜测反应路径

反应物的进攻位点、相对位置对搜索过渡态有决定性影响。不同的进攻位点，有可能对应不同的反应。

这一步对化学直观经验和感觉的要求较高，主要的技术含量也在这里。猜测的好，计算能够少费很多力气，猜测的不好，就会导致无用的计算太多。
如果经验足够丰富的话，可以不用从完全的反应物一端开始计算。打个比方，$F_2$的F原子在反应物一端与C原子距离可能很远，例如4埃，逐渐靠近C，从而接近过渡态，如果用户大概估计出过渡态处F-C距离是1.5埃，那么只需要从F-C距离为1.8埃开始，扫描到1.3埃即可。而不需要从4埃扫描到1.1埃（估计的产物键长）。毕竟扫描的范围越大，需要扫描的点数就越多，计算也就越耗时间

在本例中，我们猜测F₂进攻CH₄的姿态如下图：

上图中，按住Shift键，选中F、C原子，右下角出现两个原子之间的距离，之所以这个距离是因为我们估计过渡态键长大约1.4~1.6埃，所以扫描范围略比这个范围大，确保过渡态键长在扫描范围内即可。因此我们将从该距离，扫描到1.3埃。如果缺乏相关经验，则应老老实实地确保扫描起点为反应物，终点为产物。

PES Scan的计算原理，参考势能面扫描PES。参数设置如下：

因为精度不需要高，所以这里我们选择计算量小的GGA，使用小基组DZP（较重的元素，可以使用TZP，特别重的元素，例如锕系可以使用TZ2P）、Frozen Core：Large。

按住Shift键，选中F、C两个原子，点击右边窗口“……(distance)”前面的加号：

从而可以设置F-C键长初始值为当前值（1.851埃），末值为1.3埃，从1.851到1.3之间扫描5个点应该就足够了：

点数越多，势能面越光滑，但计算量也越大，这里基本上大约每隔0.11埃左右就扫了一个点，已经非常光滑了。

扫描的参数可以是键长，也可以是键角、二面角，也可以多个参数配合。建议尽量只用一个“决定性”的参数即可，该参数的变化，决定了反应物到产物的变化。

保存任务，并提交作业。

计算过程中，为了节省计算时间，可以在计算过程中，在SCM - Movie中监控计算过程（如果在服务器计算，服务器没有GUI，可以将*.logifle下载到本地用SCM - Movie打开），往往并不需要算完，即可得到需要的结果，例如这里：

其中红色曲线是能量变化曲线，整个扫描过程，我们只关注收敛的点，也就是上图中圆圈所指示的前4个点。可以看到前4个点，随着F-C键长缩短，能量逐渐升高，第5个点还没有收敛，但是看到当前能量已经低于第4个点了，收敛后，能量必定更低。

因此，可以确定第4个点，就是我们想要的——过渡态结构的初始猜测。

因此基于这个结构，进行第3步计算。可以在logfile中，直接复制这个结构对应的坐标，新建Input窗口，然后ctrl v粘贴到该Input窗口，或者在Movie窗口，鼠标拖动到该帧，File - Save Geometry，保存为*.xyz文件。新建Input窗口，File - Import Coordinates导入该*.xyz文件。然后设置第三步的参数，如下。

第三步，计算频率

参数设置：

这一步计算频率的作用有两个：

通过虚频（频率为负数）的情况，看第2步得到的这个结构，是否确实接近过渡态：如果频率最大、强度最大的虚频，其振动方式是分别朝着反应物和产物变化，那么这就是“过渡态（或者说鞍点）”的很好初猜；
得到该结构的能量梯度（保存在*.results/adf.rkf中），第4步将读取该文件，促进过渡态搜索更快收敛。

得到的虚频，可以用SCM LOGO > Spectra打开（在adfjobs窗口，选中任务列表中该相，变成灰色，然后点击SCM LOGO > Spectra），如下图：

下面的表格列出了振动的频率、强度，对于虚频（频率为负），点击该行，即可看到对应的振动模式动画。进而检查是否符合过渡态特征。如果最大、最强的这个虚频确实满足过渡态特征，我们就可以继续基于该分子结构真正开始搜索过渡态，也就是下一步。

如果有多个虚频，只需要看波数最大、强度明显最大的虚频是否满足过渡态特征即可，有其他小虚频没有关系
如果遇到两个虚频强度接近，则表示该结构非常劣质，用于过渡态初始猜测将很难收敛到精确过渡态结构。如何消除不需要的那个峰呢？参考如何消除虚频

4，过渡态搜索

这一步计算，仍然采用第三步的分子结构，参数设置：

前一步计算生成了一个*.t21文件，这里如果能够在Details > Files (Restart)读取该文件，则能大大提供过渡态搜索的效率、成功率：

计算收敛之后，*.logfile显示最大梯度（constrained gradient max）接近为0：

<Nov13-2020> <18:53:35>  energy change                      -0.00000680     0.00007000    T
<Nov13-2020> <18:53:35>  constrained gradient max            0.00052115     0.00100000    T
<Nov13-2020> <18:53:35>  constrained gradient rms            0.00019580     0.00066667    T
<Nov13-2020> <18:53:35>  gradient max                        0.00052115
<Nov13-2020> <18:53:35>  gradient rms                        0.00019580
<Nov13-2020> <18:53:35>  cart. step max                      0.00222910     0.01000000    T
<Nov13-2020> <18:53:35>  cart. step rms                      0.00067977     0.00666667    T

同时检查结构，确实是一个非常不稳定的状态，因此可以基本确定这是过渡态。但完全的确认过渡态，还需要再次计算频率，来进行验证。

5，虚频（频率为负数）验证，计算Gibbs自由能

基于第4步最后收敛的结构，计算频率（注意整个过程最好保持相同的基组、泛函、积分精度、冻芯情况）

设置热力学性质计算参数——温度、压强，从而得到该分子结构的Gibbs自由能，用于计算活化能：

计算结束后，也如同如第3步一样观察虚频的情况，发现只有一个虚频，并且振动模式比较符合我们的期望（分别朝反应物、产物振动）。即确认这确实是一个过渡态结构。

SCM - Output - Other properties - Statistical Thermal Analysis，即可看到Gibbs自由能为-566.06kcal/mol：

     Temp                                                        Transl      Rotat     Vibrat     Total
     ----                                                        ------      -----     ------     -----
 
     298.15   Entropy (cal/mol-K):                               37.883     22.667      7.632     68.182
              Nuclear Internal Energy (kcal/mol):                 0.889      0.889     26.892     28.669
              Constant Volume Heat Capacity (cal/mol-K):          2.981      2.981      9.094     15.055
          (c) Constant Volume Heat Capacity (cal/mol-K):          2.981      2.981      8.872     14.834
 
       Summary of energy terms
                                                   hartree              eV         kcal/mol           kJ/mol
                                      --------------------     -----------       ----------      -----------
      Energy from Engine:               -0.916314999066729        -24.9342          -575.00         -2405.78
      Nuclear Internal Energy:           0.045687613803229          1.2432            28.67           119.95
  (c) Nuclear Internal Energy:           0.045553626970109          1.2396            28.59           119.60
      Internal Energy U:                -0.870627385263500        -23.6910          -546.33         -2285.83
      pV/n = RT:                         0.000944186013486          0.0257             0.59             2.48
      Enthalpy H:                       -0.869683199250014        -23.6653          -545.73         -2283.35
      -T*S:                             -0.032395578008292         -0.8815           -20.33           -85.05
  (c) -T*S:                             -0.032380708844105         -0.8811           -20.32           -85.02
      Gibbs free energy:                -0.902078777258306        -24.5468          -566.06         -2368.41

6，提取虚频的正方向和反方向振动的第3帧，保存二者的几何结构

在Spectra窗口点击该唯一虚频，Play - Open Mode in AMSMovie可以将该振动模式的动画单独窗口显示出来。将结构从第0帧开始，手动地向右移动3帧，之后暂停在该帧，然后File-Save Geometry保存结构；以及从第0帧开始像向左移动3帧，大致到达录像的倒数第3帧）之后暂停在该帧，然后File-Save Geometry保存结构。

两个结构与“过渡态”相比，分别靠近反应物和产物方向。对这两个结构，分别进行几何结构优化（对应下载的文件中的Reactant和Product），分别得到反应物和产物的结构。

如果确实是我们所预想的反应物和产物，那么就表示这个过渡态是该反应物、产物之间的过渡态，否则意味着该反应可能不是一步完成的，中间可能还有其他单步反应。

优化结束后，使用各自优化好的分子结构，重新计算在指定温度、压强下的热力学性质，从而得到Gibbs自由能。将上面过渡态得到的Gibbs自由能与反应物的Gibbs自由能相减，即得到活化能。

Gibbs自由能的计算，参考：气相分子的频率、红外IR、零点能、转动能量、转动惯量、熵、焓、热容与Gibbs自由能，液相自由能与焓的计算

活化能

简而言之就是过渡态的“能量”-反应物的“能量”。

如果活化能是应用于厄伦尼乌斯公式，用来计算反应速率，那么“能量”一般采用计算结束之后得到的Total Bonding energy（*.logfile尾部的“Bond Energy”、*.out文件中“Total Bonding Energy”），加上零点能（SCM → Output → Other Properties → Zero-Point Energy ，计算频率默认就可以得到该值）。

如果考虑外部条件，则可以用Gibbs自由能的计算时，得到的焓（Enthalpy H）来相减。

如果用Gibbs自由能相减，这个能量差值，是不能用于厄伦尼乌斯公式的。

总结

通过ADF的过渡态搜索，能够搜索到复杂的反应过程。也能够顺利地在搜索过程中，发现可能存在的中间反应过程。ADF搜索过渡态的成功率几乎可以达到100%。
注意在计算能量的时候，需要用相同的参数（基组、泛函、数值精度等）。因此，一般在优化好反应物、产物、过渡态之后，可以统一的比较高精度的参数计算一下三个结构的单点能或Gibbs自由能作为三者的能量。
用户也可以基于收敛的过渡态、过渡态结构的频率计算生成的*.results/adf.rkf文件，提高精度重新再搜索一次过渡态，这样会非常快地收敛。在统一精度下计算反应物、产物、过渡态的能量。
即使初始的反应路径猜测的不是很准确，甚至是错误的，最后也能通过这个过程，得到正确的反应路径。

AMS软件提供免费试用，试用申请方式参见：AMS免费试用

目录