这是本文档旧的修订版！
用简谐近似下的过渡态理论研究反应速率
介绍
在本实例中，你将学习使用VNL利用NEB方法和简谐近似下的过渡态理论（HTST）来计算反应速率。这两项技术的有力结合可以让研究者对固态反应建立模型（或者更一般的说任何小概率事件动力学）。典型地，固态中的反应很慢以至于无法通过分子动力学（MD）模拟进行高效的研究。HTST，从另一方面，运用统计力学方法计算反应速率。所以，与MD模拟相比，HTST计算反应速率需要的计算量与反应时间尺度无关。
理论
HTST[1]源自于一个众所周知的理论——过渡态理论（TST）[2][3]。TST原理是通过定义过渡态（一个超平面）将反应物与产物分开。反应速率是这样计算的，相对于整个反应物状态找到系统处于过渡态的概率，乘上在过渡态处从反应区到产物区的穿过速率。总的来说，TST很难使用，因为它需要对过渡态超平面进行定义。而且，计算在过渡态的概率需要对过渡态区和反应区完整的配分函数进行估计，这在使用密度泛函理论时难以计算，即使对于速度很快的经典势方法也是个挑战。
TST可以通过两个简化假定成为一个更有实用价值的理论。第一个假定是势能面上的鞍点代表了分隔反应物和生成物的瓶颈区。这解决了在TST中确定超平面这个问题（将其定义为鞍点处垂直于不稳定方向的超平面）。第二个假定是反应物构型附近和鞍点附近的势能面是局域简谐的（二次型的）。这使得这些点的势能面可以用二阶泰勒级数展开，并且得出的在TST中所需的构型积分是显式表达式。这两个假定的结合被称为HTST。 HTST对于很多固态反应来说是一个很好的近似，因为能量位垒通常比系统中的平均动能$k_{\rm B} T$ 要大。这很重要，因为它确保了系统在反应事件之间达到局部平衡。如果能垒相对热能较小，那么对于轨迹就会有动态关联。这意味着轨迹的历史是决定反应速率的一个很重要的因素，因此需要使用更为复杂的统计模型。另一个原因是固态反应之所以能被HTST很好的描述是因为简谐近似通常凑效。接近鞍点势能迅速升高并产生一个小的瓶颈区，然而，在某些“软”系统（比如蛋白质结构变化）势能在鞍点周围区域相对平缓，从而不能用简谐展开对鞍点进行表征。 使用HTST的第一步是去寻找所研究的反应对应的鞍点。寻找鞍点可以使用很多种方法，但本实例我们将用NEB。NEB优化可以实现某一个构象收敛于鞍点（而不是将所有构象沿着反应坐标均匀隔开）。这被称作climbing image NEB (CI-NEB)，并且一个优化的CI-NEB应该有一个鞍点作为其最高能量构象。 得到鞍点后可以通过以下公式来计算反应速率： $$k_{\rm HTST} = \frac{\prod_i^{3N} \nu_i^{\rm R}}{\prod_i^{3N-1} \nu_i^{\rm S}}\, \exp \left[ -\left(E^{\rm S} - E^{\rm R}\right)/k_{\rm B} T \right]$$ 其中N是原子数；$\nu^{\rm R}_i$ 和$\nu^{\rm S}_i$ 是在反应物和鞍点稳定的（实际值）简正模频率（鞍点的$3N-1$是由于在鞍点存在一个虚频）；$E^{\rm R}$和$E^{\rm S}$ 是最小点和鞍点的能量；$k_{\rm B}$是玻耳兹曼常数；$T$是温度。 HTST的函数形式与经验得出的阿累尼乌斯反应速率指数方程相似： $k_{\rm Arrhenius} = A \exp \left[ -E_{\rm a} / k_{\rm B} T \right].$ 指前因子，A，在阿累尼乌斯速率方程对应最小点振动模式的乘积比上在HTST鞍点的稳定振动模式的乘积。这可以被解读为一个“尝试频率”，也就是，系统沿反应方向每秒振动的次数。 活化能，$E_{\rm a}$，在阿累尼乌斯速率方程中被定义为在HTST中鞍点和最小点之间的能量差。指数项代表了沿着反应方向的振动克服位垒并继续到产物态的玻尔兹曼概率。 所以，HTST速率方程可被解读为： $k_{\rm HTST} = [\text{Attempt Frequency}] \cdot[\text{Probability of Success Per Attempt}]$
注意！
在多数体系中，指前因子在不同反应不会显著变化（大约10倍的变化）。然而，指数项（处于过渡态的系统玻尔兹曼概率）可以轻松跨越10个数量级。这是因为多数固态反应是热激发过程，它们的速率主要取决于能量位垒。
尽管前因子对不同机制之间速率的差别贡献很小，计算它比计算能量位垒要花费更多的计算资源。
所以推荐对系统进行测试以判断是否有必要每次都计算前因子。在有轻原子（元素周期表前两行）参与的反应中，前因子通常是在1014Hz的量级，而对重金属原子它可能会小至1011Hz。这些频率与这些系统的特有振动周期（主要取决于原子质量）有关。
Pt吸附原子在Pt（100）表面扩散的建模
作为一个模型体系，我们将研究Pt吸附原子在Pt（100）表面的扩散。对于吸附原子在表面的扩散有两种众所周知的机制。第一种是“跳跃”机制，吸附原子直接移动到一个邻近的空位。第二种是“交换”机制，吸附原子与一个表面原子交换位置然后将其放置在一个邻近的空位。
对于这个系统我们将要学到几个很有趣的事情。首先是尽管DFT和经典势通常对平衡态性质算的不错，经典势在这里很合适，在鞍点附近的吻合度却通常是很糟糕的。通常鞍点描述了断键和成键步骤，这通常在经典势中很难被描述，尤其是像经常被用来描述Si和C键的一些使用表观角度项的势。其次是初看上去是一个简单的单步反应有时实际上包含了很多步。在所选的初始和最终几何构型间遇到一个相对不稳定的反应中间体是很常见的。一个必须理解的要点是， HTST在形式上只描述简单的单步反应。
构造Pt（100）表面
在Builder 中点击Add ‣ From Database然后搜索“platinum”。含有铂的结构列表将会出现，滚动并寻找“Platinum”然后双击它。 点击Builders ‣ Surface (Cleave)然后设置密勒指数为（100）并点击Next。创建一个2*2的表面晶格并点击Next。然后选择Out-of-plane cell vector 为 Non-periodic and normal (slab) ，设置厚度为5层，点击Finish。
为了之后的可视化，最好将平板模型在x-y平面居中。点击Coordinate Tools ‣ Center，点掉C方向并点击Apply。 在第五层删除多余原子，只保留一个原子。剩下的原子将会是NEB计算中的扩散吸附原子。
构造NEB构型
通过点击Stash中的Copy按钮将构型拷贝。双击新的构型。选择吸附原子并点击工具栏中的move tool将其拖到一个邻近的空位。如果你将视角方向对着xy平面，这将很容易。
点击Builders ‣ Nudged Elastic Band并将stash中的两个构型拖到initial和 final放置区。点击Constraints 并选择最底层。点击Selection中的Add tag，然后将约束类型从None改为Fixed并点击OK。将NEB方法从Linear改为 Image Dependent Pair Potential。在本实例中，5个构象就可以很有效了，因为此反应坐标相对简单。设置Maximum distance为1 Å并点击Create。 右击Stash中的NEB构型并选择Send To ‣ Script Generator。
优化NEB与HTST速率计算
双击模块栏中的New Calculator来为脚本添加一个计算工具。双击这个新加的New Calculator 。改变k-point sampling为3*3*1并将交换关联函数设为GGA。点击OK。
双击模块栏里的Optimization ‣ OptimizeGeometry并在脚本栏中双击它。改变Force Tolerance为0.01 eV/Å。然后点击Add Constraints，选择最底部层的原子，从Selection中点击Add tag ，将约束类型从None 改为Fixed。选中Climbing Image Method和Optimize end points勾选框。
运行CI-NEB计算是很重要的。如果使用一个常规的NEB优化，构象将会在反应坐标中均匀隔开。这意味着很有可能在鞍点上没有构象。而一个收敛的CI-NEB计算应该有一个鞍点作为最高能量构象。
注意！
在本实例中，我们在一步之内优化NEB端点和中间体构象。这对于本例（端点接近于它们的优化结构）来说很有效，然而，将这个过程分为三步通常更好：
- 优化初始端点。
- 优化最终端点。
- 从优化的端点创建NEB构型并对其优化。
这个方法允许你去分析并核实优化的端点来确保几何优化不会显著地改变它们。
双击模块栏中的Analysis ‣ HTSTEvent并在脚本栏中双击它。有两个选项来确定HTST速率方程的前因子（见介绍部分的第一个方程）。
计算指前因子的第一种方法是计算系统中原子的一个子集的动力学矩阵。包含在计算中的这些原子同最近邻原子（使用共价半径来决定）一起沿着反应坐标移动超过Minimum Displacement参数。第二种方法是简单地设置前因子为一个已知量。
改变默认输出文件为 pt hop.nc。然后保存脚本使用atkpython运行或者将其送到Job Manager运行。
计算可能会花费几个小时，这取决于你的硬件和所使用的MPI进程个数。
注意！
收敛的CI-NEB计算的最高能量构象可能不是一个一阶鞍点。这可以通过计算鞍点的频率来检测。日志里将会写入一个报错信息，计算的频率也会显示出来。出现这种情况有以下几个原因：
- NEB优化没有收敛。考虑增加在优化里的最大步数，并检查输入几何构型来确保初始猜测是合理的。
- NEB优化也许收敛了，但是force tolerance太大了。尝试将force tolerance变小。对于很多体系 $10^{-2}eV/Å$ 是一个很好的值，但是有时候需要使用 $10^{-3}eV/Å$。
- 如果鞍点的虚频仍然存在，那么可能是NEB收敛于一个二阶鞍点（也就是一个拥有两个虚频的鞍点）。如果是这样，你需要为NEB构型建一个不同的初始猜测来看它是否会收敛于正确的一阶鞍点。
结果分析
计算完成后，文件pt hop.nc 应包含初始和收敛的NEB构型，以及HTSTEvent。点击第二个NEB构型并点击右侧面板的Movie Tool。
势能曲线表明能垒为1.3eV。反应的动画展示了预期的跳跃过程。在初始态和末态之间有一个单能垒。关闭Movie Tool窗口。选择LabFloor中的HTSTEvent对象并点击右侧面板的HTST Rates插件。
显示了关于反应的信息如能垒和指前因子。可以通过改变温度值来计算任意温度下的速率，默认温度为300K。
正反应和逆反应的能垒和前因子在理论上应该相同，然而在实际过程中是不同的。这是由于NEB端点几何构型的微小不同（它们每一个都是从一个略微不同的开始结构被最小化的）。而且，由于这是一个ATK-DFT计算，由于eggbox效应，平移对称性没有被完美的保持。这个效应描述了用以表示电荷密度和势的实空间网格与原子位置之间的关系。这些网格点会随着不同原子而向不同方向对齐，从而打破了平移对称性。在本实例中对于总能影响非常小（<1meV），但是对于前因子影响很大，变化达到1.6倍。可以通过增加在ATK-DFT计算器设置中的Density mesh cut-off 来减小这些误差，但是这样会增加模拟的计算成本，而且对于得到有意义的结果并不是必要的。
我们也可以画一个Arrhenius图来显示HTST速率随着温度的变化情况。一个Arrhenius图是$\log({\rm rate})$与$1/T$的关系。当与温度的关系以这种方式绘出，图应该总是直线。直线的斜率是$-\Delta E/k_{\rm B}$而y轴截距是指前因子。从图可以很容易地说明指前因子是HTST预测的最大可能速率。
Pt的跳跃扩散反应速率在室温下非常低。反应之间的平均时间（速率的倒数）是在世纪量级的。然而，在600K我们可以期望反应每秒发生几十次。