本文是「第一性原理的微观计算模拟」系列的第 9 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
拖曳法与 NEB 方法
体系由一个状态跃迁至另一个状态的过程中，需要克服某种形式的能量势垒。设一个体系的自由度为 $N$，则体系的位置由一个 $N$ 维矢量描述，也即体系处在一个 $N$ 维空间中。相邻的两个稳态体系的坐标分别为 $\boldsymbol{R}_1^{N}$ 与 $\boldsymbol{R}_2^{N}$，连接这两个稳态的路径有无限多条。因此，跃迁路径特指最小能量路径（minimum energy path，MEP）。沿着这条路径前进，体系只需要越过最低的势垒就可以完成跃迁。而 MEP 的最高点是体系的一个一阶鞍点。在该点处，能量沿着最小能量路径方向达到极大值，而沿其他任何方向均是极小值。这意味着此处的声子谱只包含一个虚频的振动模式。该模式中原子沿着最小能量路径振动。相应地，过渡态（transition state，TS）特指最小能量路径上的这个一阶鞍点。图 3.16 为典型的一阶鞍点与最小能量路径示意图，图中圆点代表微动弹性带（nudged elastic band，NEB）方法中的映像。
图 3.16 一阶鞍点与最小能量路径示意图
过渡态是材料老化、变化以及化学反应过程中非常重要的一个概念，它直接反映了原子尺度上微观过程发生的路径与难易程度。第 8 章将要介绍的动态蒙特卡罗等方法也需要以过渡态（TS）相对于初态（IS）或者末态（FS）的能量变化作为参数进行大尺度的模拟。在前面的讨论中我们已经给出了原子受力的表达式，这使得直接对跃迁／反应路径进行原子模拟成为可能。但是，因为跃迁／反应的初态及末态均为稳态，也即在体系能量极小值处，即使将体系人为进行偏移，经过弛豫后体系也必然会自动回复到稳态。因此，必须加入额外的限制条件才能完成对跃迁／反应路径的模拟。
最为简单和直观的做法是拖曳法（drag method）。顾名思义，就是将体系由初态拖曳 $m$ 步至末态，产生 $m$ 个复制体系。之后，对每一个复制体系，固定沿着拖曳路径的自由度，而弛豫其他 $N-1$ 个自由度，寻找最小值。最后，取其中能量最高的一个复制体系所处的位置作为过渡态，将其能量作为势垒。拖曳法中每个复制体系均是独立弛豫的，因此对于某些原子数目较多的体系所需要的计算资源并不是太多。对于特定的反应路径，例如间隙原子在表面或块体内的迁移，拖曳法也往往能给出比较合理的结果。但是，这种方法最大的缺点就是可信度较低，而且失败率比较高，对很多路径都无法给出正确的反应路径。设体系初态为 $\boldsymbol{R}_{\mathrm{Ini}}^{N}$，末态为 $\boldsymbol{R}_{\mathrm{Fin}}^{N}$，则一般情况下拖曳法给出的初始路径是连接两者的直线 $\boldsymbol{V}^{N}=\boldsymbol{R}_{\mathrm{Ini}}^{N}-\boldsymbol{R}_{\mathrm{Fin}}^{N}$。因为 $m$ 个复制体系彼此独立，所以拖曳法实质上找到的是 $m$ 个彼此平行且垂直于 $\boldsymbol{V}^{N}$ 的超平面中的能量最小值。这无法保证找到的极大值在真正的一阶鞍点附近。实际上，如果鞍点处虚频对应的振动方向与拖曳路径夹角较大，拖曳法将无法给出正确的最小能量路径\cite{schwartz2000theoretical}。下面将要介绍的 NEB 方法可以很好地弥补拖曳法的这些缺陷。
Mills、Jónsson 等人提出的 NEB 方法可以给出含有多个鞍点的最小能量路径\cite{mills1994quantum,mills1995reversible,berne1998classical}。设体系在初态和末态之间移动，每个位置称为一个映像，现在假设这些映像同时出现在反应路径上，彼此之间由刚度系数为 $k$ 的弹簧连接。处于两个端点的初态和末态固定不动，其他位于中间的映像可以放开所有的自由度进行弛豫。与拖曳法不同，NEB 方法中的映像彼此之间通过弹簧耦合，而且参与弛豫的映像由于受到弹簧的阻力不会滑落回端点。在这种情况下，映像 $i$ 的受力相当于是
$$ \boldsymbol{F}_i=-\boldsymbol{\nabla}E(\boldsymbol{R}_i^{N})+k(\boldsymbol{R}_{i+1}^{N}-\boldsymbol{R}_i^{N})-k(\boldsymbol{R}_i^{N}-\boldsymbol{R}_{i-1}^{N}) \tag{3.579} $$
但是实践表明，应用上述方程时经常会出现两个问题。首先，充分弛豫后，能量较高的鞍点附近映像分布非常稀疏，而靠近端点的能量较低处映像分布比较集中，从而导致更有物理意义的鞍点附近的最小能量路径分辨率比较低。这种现象称为映像滑落（down-sliding），其原因是各个映像上的真实受力（例如 Hellmann-Feynman 力）沿弹簧方向的分量倾向于将各映像推向能量极小的端点处。其次，当 MEP 曲率较大时，因为弹簧将在垂直于相邻映像连线的方向上提供额外的力，所以该区域的映像会偏离实际的最小能量路径，而按照较为平直的路径分布。这种现象称为截弯（corner-cutting）。
为了解决上述两个问题，需要对映像上的受力进行投影。每个映像的受力按照径向和法向分为两部分：在每个映像上都可以定义一个单位超正切矢量 $\hat{\boldsymbol{\tau}}_i$ 作为径向，这个方向上的受力由连接两者的弹簧决定；垂直于该连线的方向为法向，沿法向的受力由该映像所处的势能面在该方向上的梯度决定。因此，在 NEB 方法中，映像 $i$ 的受力为
$$ \boldsymbol{F}_i=\boldsymbol{F}_i^{\perp}+\boldsymbol{F}_i^{s,\parallel} \tag{3.580} $$
其中
$$ \boldsymbol{F}_i^{\perp}=-\boldsymbol{\nabla}E(\boldsymbol{R}_i^{N})+\boldsymbol{\nabla}E(\boldsymbol{R}_i^{N})\cdot\hat{\boldsymbol{\tau}}_i\hat{\boldsymbol{\tau}}_i \tag{3.581} $$
$$ \boldsymbol{F}_i^{s,\parallel}=k(|\,\boldsymbol{R}_{i+1}^{N}-\boldsymbol{R}_i^{N}\,|-|\,\boldsymbol{R}_i^{N}-\boldsymbol{R}_{i-1}^{N}\,|)\hat{\boldsymbol{\tau}}_i \tag{3.582} $$
显然，定义映像上的径向矢量 $\hat{\boldsymbol{\tau}}_i$ 对最后的结果会有很大的影响。通常情况下可以通过与第 $i$ 个映像相连的 $i-1$ 和 $i+1$ 两个映像确定 $\hat{\boldsymbol{\tau}}_i$：
$$ \hat{\boldsymbol{\tau}}_i=\frac{\boldsymbol{R}_i^{N}-\boldsymbol{R}_{i-1}^{N}}{|\,\boldsymbol{R}_i^{N}-\boldsymbol{R}_{i-1}^{N}\,|}+\frac{\boldsymbol{R}_{i+1}^{N}-\boldsymbol{R}_i^{N}}{|\,\boldsymbol{R}_{i+1}^{N}-\boldsymbol{R}_i^{N}\,|},\quad\hat{\boldsymbol{\tau}}_i=\frac{\boldsymbol{\tau}_i}{|\,\hat{\boldsymbol{\tau}}_i\,|} $$
但是对于部分原子成键方向性较强的体系（如 Si 晶体或者 Ir(111)-CH4 体系等），这样给出的 $\hat{\boldsymbol{\tau}}_i$ 往往会导致 NEB 模拟不收敛\cite{henkelman2000improved}。Henkelman 和 Jónsson 为解决这个困难提出了改进的径向矢量\cite{henkelman2000improved}：
$$ \boldsymbol\tau_i=\begin{cases} \boldsymbol R_{i+1}^N-\boldsymbol R_i^N,&E_{i+1}\gt E_i\gt E_{i-1},\\ \boldsymbol R_i^N-\boldsymbol R_{i-1}^N,&E_{i+1}\lt E_i\lt E_{i-1}. \end{cases} \tag{3.583} $$
如果映像 $i$ 处于能量极值处，则
$$ \boldsymbol\tau_i=\begin{cases} (\boldsymbol R_{i+1}^N-\boldsymbol R_i^N)\Delta E_i^{\max} +(\boldsymbol R_i^N-\boldsymbol R_{i-1}^N)\Delta E_i^{\min},&E_{i+1}\gt E_{i-1},\\ (\boldsymbol R_{i+1}^N-\boldsymbol R_i^N)\Delta E_i^{\min} +(\boldsymbol R_i^N-\boldsymbol R_{i-1}^N)\Delta E_i^{\max},&E_{i+1}\le E_{i-1}. \end{cases} \tag{3.584} $$
式中
$$ \Delta E_i^{\max}=\max(|\,E_{i+1}-E_i\,|,|\,E_{i-1}-E_i\,|) \tag{3.585} $$
$$ \Delta E_i^{\min}=\min(|\,E_{i+1}-E_i\,|,|\,E_{i-1}-E_i\,|) \tag{3.586} $$
最后再将 $\tau_i$ 归一化：$\hat{\tau}_i=\tau_i/|\,\tau_i\,|$。利用这种改进的径向矢量定义，在包含足够映像数目的条件下，NEB 方法可以在大多数情况下收敛。
式（3.580）表明，在 NEB 方法中，体系受力并不等于能量对位置的导数的负值。这一事实使得 NEB 方法中所采用的优化算法不同于 1.3 节中所介绍的方法。在附录 A.7 节中，我们将简要介绍两种 NEB 常用的优化方法。除此之外，也可以采用最速下降法、共轭梯度法或者拟牛顿法。与 1.3 节中给出的算法不同，在进行一维搜索时，NEB 中采用的优化方法不寻找能量最低值，而是利用牛顿方向找到受力为零的点：沿优化方向取两点 1 和 2，各自计算映像 $i$ 的受力 $\boldsymbol{F}_{i,1}^{N}$ 和 $\boldsymbol{F}_{i,2}^{N}$，然后利用有限差分计算 $\boldsymbol{F}_i^{N}$ 的导数。关于这部分内容更详细的讨论可参阅文献\cite{sheppard2008optimization}。
近年来，基于 NEB 的其他寻找过渡态的方法，如 CI-NEB 方法、DNEB 方法等最近有了新的发展。我们在这里不详细介绍，有兴趣的读者可参阅文献\cite{henkelman2000climbing,trygubenko2004doubly}。
Dimer 方法
拖曳法和 NEB 方法要求预先知道初态和末态。如果采用通常的线性插值方法来寻找过渡态的初始路径，在一定意义上相当于预设一种跃迁机制，然后利用 NEB 方法确定这种跃迁机制的势垒。但是，这种寻找方法并不能保证找到的势垒是最低的。例如，W 原子在 W(001) 表面上跃迁，势垒最低的路径对应于吸附的 W 原子陷入表面，形成表面间隙原子，然后沿[100]或者[010]方向移动，若干步之后再转变为吸附 W 原子。这种跃迁机制称为表面挤列（crowdion）机制，如图 3.17(b) 所示。而此前认为的跳跃（hopping）机制是吸附 W 原子“跃过”表面原子到达下一个吸附位（所对应的势垒要高 0.7 eV\cite{chen2013biaxial}），如图 3.17 所示。对于在 Al(001) 表面上的自扩散也有类似的发现\cite{feibelman1990diffusion}。更为严谨的方法是只从已知的初态出发，通过一定的算法自动地寻找所有可能的跃迁路径。这种方法就是 Henkelman 等人提出的 Dimer 方法\cite{henkelman1999dimer}。
图 3.17 W 原子在 W(001) 表面上跃迁的两条典型路径
（a）跳跃机制；（b）表面挤列机制
在 Dimer 方法中仍然需要两个映像，但是只有其中一个映像要求是稳态（作为初态），另一个映像可以通过给初态的原子坐标加上一个方向随机的微扰而生成，或者由初态出发，在有限温度下按照分子动力学原理（将在第 6 章中详细讨论）生成一条轨迹，由其中某一时刻的即时构型给出。这两个映像组成一个偶矩（dimer），这个偶极矩在高维势能面中通过转动以及平动等运动方式经过鞍点，到达另一个势阱处（末态）。因为所加的微扰不同，偶极矩也会通过不同的途径到达多个终态。在尝试次数足够多的情况下，Dimer 方法可以找到连接给定初态的所有跃迁途径（更准确地说是势垒最低的跃迁路径方向上几个 $k_{\mathrm{B}}T$ 范围内的鞍点）。这无疑是 Dimer 方法非常有吸引力的一个优点。
图 3.18(a) 为偶极矩示意图。两个映像的位置、能量和受力分别为 $\boldsymbol{R}_1$、$E_1$、$\boldsymbol{F}_1$ 和 $\boldsymbol{R}_2$、$E_2$、$\boldsymbol{F}_2$。单位矢量 $\hat{\boldsymbol{N}}$ 由 $\boldsymbol{R}_2$ 指向 $\boldsymbol{R}_1$。该偶极矩的中点为 $\boldsymbol{R}$。因此有
图 3.18 Dimer 方法的应用
（a）偶极矩示意图以及受力的投影；（b）偶极矩的转动步骤
$$ \boldsymbol{R}_1=\boldsymbol{R}+\Delta R\hat{\boldsymbol{N}},\quad\boldsymbol{R}_2=\boldsymbol{R}-\Delta R\hat{\boldsymbol{N}},\quad\Delta R=\frac{1}{2}(\boldsymbol{R}_1-\boldsymbol{R}_2)\cdot\hat{\boldsymbol{N}} $$
注意：实际上偶极矩是在一个 $3N$ 维的空间里定义的，而非图 3.18 所示的仅由三维空间中的两个点确定。体系总能量 $E=E_1+E_2$。偶极矩中心的能量为 $E_0$，受力为 $\boldsymbol{F}_{\mathrm{R}}$，定义 $\boldsymbol{F}_{\mathrm{R}}=(\boldsymbol{F}_1+\boldsymbol{F}_2)/2$。由此可以通过有限差分以及定义计算此处势能面的曲率：
$$ C=\frac{(\boldsymbol{F}_2-\boldsymbol{F}_1)\cdot\hat{\boldsymbol{N}}}{2\Delta R}=\frac{E-2E_0}{(\Delta R)^{2}} \tag{3.587} $$
在 Dimer 方法中，每一步分为两个部分，即先转动偶极矩，使得其平行于势垒最低的跃迁路径，之后平移偶极矩，使其沿跃迁路径到达鞍点。这两部分均需采用优化算法。下面分别进行讨论。
1. 转动
由式（3.587）可知，偶极矩的总能量 $E$ 和势能面曲率 $C$ 呈线性关系。因此在限定偶极矩仅做转动的条件下，$E$ 有最小值意味着偶极矩平行于势能增加最缓慢的方向，即势垒最低的跃迁路径。首先定义垂直于偶极矩方向的受力（转动力）$\boldsymbol{F}^{\perp}$：
$$ \boldsymbol{F}^{\perp}=\boldsymbol{F}_1^{\perp}-\boldsymbol{F}_2^{\perp} \tag{3.588} $$
式中
$$ \boldsymbol{F}_i^{\perp}=\boldsymbol{F}_i-(\boldsymbol{F}_i\cdot\hat{\boldsymbol{N}})\hat{\boldsymbol{N}},\quad i=1,2 $$
设 $\boldsymbol{F}^{\perp}$ 作用在 $\boldsymbol{R}_1$ 上，且定义一个与 $\boldsymbol{F}^{\perp}$ 平行的单位矢量 $\hat{\boldsymbol{\varTheta}}$。$\hat{\boldsymbol{\varTheta}}$ 和 $\hat{\boldsymbol{N}}$ 组成了展开偶极矩转动平面 $S$ 上的一组正交基矢。如图 3.18(b) 所示，将偶极矩转动一个小角度 $\mathrm{d}\theta$，则
$$ \begin{cases} \boldsymbol{R}_1^{*}=\boldsymbol{R}+(\hat{\boldsymbol{N}}\cos(\mathrm{d}\theta)+\hat{\boldsymbol{\varTheta}}\sin(\mathrm{d}\theta))\Delta R\\ \boldsymbol{R}_2^{*}=\boldsymbol{R}-(\hat{\boldsymbol{N}}\cos(\mathrm{d}\theta)+\hat{\boldsymbol{\varTheta}}\sin(\mathrm{d}\theta))\Delta R \end{cases} \tag{3.589} $$
重新计算偶极矩中两个映像的受力 $\boldsymbol{F}_1^{*}$ 和 $\boldsymbol{F}_2^{*}$。偶极矩中心受力为 $\boldsymbol{F}^{*}=\boldsymbol{F}_1^{*}-\boldsymbol{F}_2^{*}$。此时，转动步骤的任务是利用上面已得到的信息，在 $\hat{\boldsymbol{\varTheta}}$ 和 $\hat{\boldsymbol{N}}$ 展开的平面 $S$ 内寻找转角 $\Delta\theta$，使得 $|\,\boldsymbol{F}^{\perp}\,|=0$。将势能面在平面 $S$ 内展开为二次函数 $U$，有
$$ U=E_0-(F_xx+F_yy)+\frac{1}{2}(c_xx^{2}+c_yy^{2}) \tag{3.590} $$
式中：$F_x$ 和 $F_y$ 分别为 $-\partial U/\partial x$ 和 $-\partial U/\partial y$；$c_x$ 和 $c_y$ 为沿两个不同方向的曲率。由此可将偶极矩的能量 $E$ 表示为转角 $\theta$ 的函数：
$$ E(\theta)=2E_0+\Delta R^{2}[c_x\cos^{2}(\theta-\theta_0)+c_y\sin^{2}(\theta-\theta_0)] $$
$$ =2E_0+\frac{\Delta R^{2}}{2}\{(c_x-c_y)\cos[2(\theta-\theta_0)]+(c_x+c_y)\} \tag{3.591} $$
式中一次项因为映像 1、2 关于中点 $\boldsymbol{R}$ 对称而相互抵消；$\theta_0$ 是一个常数。引入标量转动力 $F$：
$$ F=\frac{\boldsymbol{F}^{\perp}\cdot\hat{\boldsymbol{\varTheta}}}{\Delta R} \tag{3.592} $$
可以证明下述关系式\cite{henkelman1999dimer,heyden2005efficient}成立：
$$ F=-\frac{1}{\Delta R^{2}}\frac{\partial E}{\partial\theta}=A\sin[2(\theta-\theta_0)] \tag{3.593} $$
式中：$A$ 为未知常数，但是实际应用中并不需要知道 $A$ 的具体值。式（3.593）表明，$\theta_0$ 正是在平面 $S$ 内使得 $F=0$ 所需转动的角度。为了得出 $\theta_0$ 的具体计算公式，进一步计算
$$ F'=\frac{\mathrm{d}F}{\mathrm{d}\theta}=2A\cos[2(\theta-\theta_0)] \tag{3.594} $$
由此可得，如果在 $\theta=0$ 处的 $F$ 和 $F'$ 均已知（分别标注为 $F_0$ 和 $F_0'$），则使得偶极矩在平面 $S$ 内平行于曲率最小的方向所需的角度 $\Delta\theta$ 为
$$ \Delta\theta=\theta_0=-\frac{1}{2}\arctan\left(\frac{2F_0}{F_0'}\right) \tag{3.595} $$
如前所述，我们已经有了 $\theta=0$ 处的 $\boldsymbol{F}_i$、$\hat{\boldsymbol{\varTheta}}$ 和 $\theta=\mathrm{d}\theta$ 处的 $\boldsymbol{F}_i^{*}$、$\hat{\boldsymbol{\varTheta}}^{*}$，因此可以求得 $\theta=\mathrm{d}\theta/2$ 处的 $F$ 和 $F'$\cite{heyden2005efficient}：
$$ F_{\mathrm{d}\theta/2}=\frac{\boldsymbol{F}^{*}\cdot\hat{\boldsymbol{\varTheta}}^{*}+\boldsymbol{F}\cdot\hat{\boldsymbol{\varTheta}}}{2} \tag{3.596} $$
$$ F_{\mathrm{d}\theta/2}'=\frac{\boldsymbol{F}^{*}\cdot\hat{\boldsymbol{\varTheta}}^{*}-\boldsymbol{F}\cdot\hat{\boldsymbol{\varTheta}}}{\mathrm{d}\theta} \tag{3.597} $$
因此从偶极矩 $\boldsymbol{R}^{*}$ 出发，所需转动的角度为
$$ \Delta\theta=-\frac{1}{2}\arctan\left(\frac{2F_{\mathrm{d}\theta/2}}{F_{\mathrm{d}\theta/2}'}\right)-\frac{\mathrm{d}\theta}{2} \tag{3.598} $$
其与式（3.595）略有不同。最后的结果如图 3.18(b) 所示。
对上述转动步骤还可以做进一步的改进。当 Dimer 方法完成一次迭代后，下一次的转动平面由新得到的 $\hat{\boldsymbol{N}}^{**}$ 及 $\hat{\boldsymbol{\varTheta}}^{**}$ 展开，如图 3.18(b) 所示。因为每一步的 $\hat{\boldsymbol{\varTheta}}$ 均与当前偶极矩受力在法向上的投影 $\boldsymbol{F}^{\perp}$ 平行，所以相当于按照最速下降法来更新转动平面，并以此来寻找势能面上曲率最小的方向。按照 1.3.2 节中的讨论，利用共轭梯度法构建新的搜索方向会得到更高的效率。在 Dimer 方法中，也可以按照共轭梯度法来构建新的转动平面。但是因为转动步骤中采用的是偶极矩受力在法向上的投影，所以不能简单地采用式（1.86）构造共轭方向。Henkelman 等人指出，可以按下式构造新的共轭方向\cite{henkelman1999dimer}：
$$ \boldsymbol{d}_i^{\perp}=\boldsymbol{F}_i^{\perp}+\beta\,|\,\boldsymbol{d}_{i-1}^{\perp}\,|\,\hat{\boldsymbol{\varTheta}}_{i-1}^{**} \tag{3.599} $$
式中
$$ \beta=\frac{\boldsymbol{F}_i^{\perp}\cdot\boldsymbol{F}_i^{\perp}}{\boldsymbol{F}_{i-1}^{\perp}\cdot\boldsymbol{F}_{i-1}^{\perp}} $$
2. 平移
转动步骤完毕后，需要进行一次平移，即将偶极矩按照当前的取向平移，直至其到达鞍点处。与转动时步骤相同，在平移时也需要对偶极矩受力进行投影等操作，以阻止其沿势能面滑落至极小值处。因此，在平移中，使用如下力来平移偶极矩：
$$ \boldsymbol{F}^{\dagger}=\begin{cases}-(\boldsymbol{F}_R\cdot\hat{\boldsymbol{N}})\hat{\boldsymbol{N}},&C\gt 0\\\boldsymbol{F}_R-2(\boldsymbol{F}_R\cdot\hat{\boldsymbol{N}})\hat{\boldsymbol{N}},&C\leqslant0\end{cases} \tag{3.600} $$
式中：$C$ 是势能面在当前处的曲率。如果 $C\gt 0$，则偶极矩仍处于势阱（稳态）附近，将偶极矩沿连线的受力反向，使得偶极矩加速从势阱中逸出；如果 $C\leqslant0$，则偶极矩在鞍点附近，$\boldsymbol{F}^{\dagger}$ 指向鞍点，将驱使偶极矩平移至所求的鞍点处。偶极矩所移动的距离 $\Delta x$ 由式（1.88）给出。因此，平移时要求进行两次力的计算，分别在 $\Delta x=0$ 以及 $\Delta x=\delta x$ 处。利用有限差分，可以计算 $\delta x/2$ 处的力以及曲率，因此，平移过程中偶极矩沿 $\boldsymbol{F}^{\dagger}$ 方向移动的距离 $\Delta x$ 为
$$ \Delta x=-\frac{F^{\dagger}}{C^{\dagger}}=-\frac{(\boldsymbol{F}^{\dagger}\,|_{\Delta x=\delta x}+\boldsymbol{F}^{\dagger}\,|_{\Delta x=0})/2}{(\boldsymbol{F}^{\dagger}\,|_{\Delta x=\delta x}-\boldsymbol{F}^{\dagger}\,|_{\Delta x=0})/\delta x}+\frac{\delta x}{2} \tag{3.601} $$
为了控制算法的稳定性，平移时一般会设定一个上限 $\Delta x_{\max}$，若由方程（3.601）计算出的 $\Delta x\gt \Delta x_{\max}$，则强迫偶极矩沿 $\boldsymbol{F}^{\dagger}$ 方向仅移动 $\Delta x_{\max}$。
为了改进 Dimer 方法的收敛性，近年来对转动和平移的步骤和计算式有一些改进，具体请参考文献\cite{heyden2005efficient,kastner2008superlinearly}。