在本实例中，你将学习使用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}]$

作为一个模型体系，我们将研究 Pt 吸附原子在 Pt（100）表面的扩散。对于吸附原子在表面的扩散有两种众所周知的机制。第一种是“跳跃”机制，吸附原子直接移动到一个邻近的空位。第二种是“交换”机制，吸附原子与一个表面原子交换位置然后将其放置在一个邻近的空位。

对于这个系统我们将要学到几个很有趣的事情。首先是尽管 DFT 和经典势通常对平衡态性质算的不错，经典势在这里很合适，在鞍点附近的吻合度却通常是很糟糕的。通常鞍点描述了断键和成键步骤，这通常在经典势中很难被描述，尤其是像经常被用来描述 Si 和 C 键的一些使用表观角度项的势。其次是初看上去是一个简单的单步反应有时实际上包含了很多步。在所选的初始和最终几何构型间遇到一个相对不稳定的反应中间体是很常见的。一个必须理解的要点是，HTST 在形式上只描述简单的单步反应。

在 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计算中的扩散吸附原子。

通过点击 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。

双击模块栏中的 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 计算则会有一个鞍点的最高能量构象。

双击模块栏中的 Analysis ‣ HTSTEvent 并在脚本栏中双击它。有两个选项来确定 HTST 速率方程的前因子（见介绍部分的第一个方程）。

计算指前因子的第一种方法是计算系统中原子的一个子集的动力学矩阵。包含在计算中的这些原子同最近邻原子（使用共价半径来决定）一起沿着反应坐标移动超过 Minimum Displacement 参数。第二种方法是简单地设置前因子为一个已知量。

改变默认输出文件为 pt hop.nc。然后保存脚本使用 atkpython 运行或者将其送到 Job Manager 运行。

计算可能会花费几个小时，这取决于你的硬件和所使用的MPI进程个数。

计算完成后，文件 pt hop.nc 应包含初始和收敛的 NEB 构型，以及 HTSTEvent。点击第二个 NEB 构型并点击右侧面板的 Movie Tool。

势能曲线表明能垒为 1.3 eV。反应的动画展示了预期的跳跃过程。在初始态和末态之间有一个单能垒。关闭 Movie Tool 窗口。选择 LabFloor 中的 HTSTEvent 对象并点击右侧面板的 HTST Rates 插件。

显示了关于反应的信息如能垒和指前因子。可以通过改变温度值来计算任意温度下的速率，默认温度为 300K。

正反应和逆反应的能垒和前因子在理论上应该相同，然而在实际过程中是不同的。这是由于 NEB 端点几何构型的微小不同（它们每一个都是从一个略微不同的开始结构被最小化的）。而且，由于这是一个 ATK-DFT 计算，由于 eggbox 效应，平移对称性没有被完美的保持。这个效应描述了用以表示电荷密度和势的实空间网格与原子位置之间的关系。这些网格点会随着不同原子而向不同方向对齐，从而打破了平移对称性。在本实例中对于总能影响非常小（<1 meV），但是对于前因子影响很大，变化达到 1.6 倍。可以通过增加在 ATK-DFT 计算器设置中的 Density mesh cut-off 来减小这些误差，但是这样会增加模拟的计算成本，而且对于得到有意义的结果并不是必要的。

我们也可以画一个 Arrhenius 图来显示 HTST 速率随着温度的变化情况。一个 Arrhenius 图是 $\log({\rm rate})$ 与 $1/T$ 的关系。当与温度的关系以这种方式绘出，图应该总是直线。直线的斜率是 $-\Delta E/k_{\rm B}$ 而 y 轴截距是指前因子。从图可以很容易地说明指前因子是 HTST 预测的最大可能速率。

Pt 的跳跃扩散反应速率在室温下非常低。反应之间的平均时间（速率的倒数）是在世纪量级的。然而，在 600 K 我们可以期望反应每秒发生几十次。

Pt 吸附原子在 Pt（100）面扩散的交换机制是一个很有趣的学习案例，因为它不像跳跃机制那样通过一个单元反应来进行。计算的设置细节在另一个实例中给出，但是计算结果以及如何理解 HTST 中的反应将在这里阐明。 收敛的 CI-NEB 计算显示了如下的反应坐标图表：

这里你可以看到有两个鞍点因而有两个反应步。第一个能垒大约为 1.4 eV，第二个能垒约为 0.35 eV。对于这个反应的 HTST 速率计算将会涉及将这个计算分成两个分开的 NEB 构型以确定每一步反应的速率常数。

根据你试图想要构建的模型，忽略中间体态通常会是一个合理的近似。举例来说，在600 K 初始反应步的速率大约为 1.7 Hz，而第二步的反应速率大约为 1GHz（$10^{12}\exp[-0.35 eV/ 600 K k_{\rm B}]$）。所以，忽略反应中间体可能是合理的。

更细致的研究应当运行一个动力学 Monte Carlo 模拟来确定系统的动力学。只要涉及的反应步超过一个，反应动力学就不会严格的服从一阶速率方程，尽管它通常是一个有效的近似。

• [1] George H. Vineyard. Frequency factors and isotope effects in solid state rate processes. Journal of Physics and Chemistry of Solids, 3(1):121 – 127, 1957. doi:10.1016/0022-3697(57)90059-8.
• [2] Henry Eyring. The activated complex in chemical reactions. The Journal of Chemical Physics, 3(2):107–115, 1935. doi:10.1063/1.1749604.
• [3] Arthur F. Voter and Jimmie D. Doll. Transition state theory description of surface self‐diffusion: Comparison with classical trajectory results. The Journal of Chemical Physics, 80(11):5832–5838, 1984. doi:10.1063/1.446610.
• 原文：http://docs.quantumwise.com/tutorials/neb_htst.html
• 本文翻译：王吉章