过渡态、反应路径的计算方法及相关问题

1、过渡态、反应路径的计算方法及相关问题SoberevaDepartment of Chemistry, University of Science and Technology Beijing, Beijing100083, China前言：本文主要介绍过渡态、反应路径的计算方法，并讨论相关问题。由于这类算法极多， 可以互相组合， 限于精力不可能面面俱到展开， 所以只介绍常用， 或者实用价值有限但有启 发性的方法。 文中图片来自相关文献， 做了一定修改。 由于本文作为帖子发布， 文中无法插 入复杂公式，故文中尽量将公式转化为文字描述并加以解释，这样必然不如公式形式严谨， 而且过于复杂的公式只能略

2、过， 但我想这样做的好处是更易把握方法的梗概， 有兴趣可以进 一步阅读原文了解细节。 对于 Gaussian 中可以实现的方法， 文中对其在 Gaussian 中的使用 进行了一些讨论，希望能纠正一些网上流传的误区。虽然绝大多数人不专门研究计算方法， 其中很多方法也不会用到，但多了解一下对开阔思路是很有好处的。文中指的“反应”包括构象变化、 异构化、单分子反应等任何涉及到过渡态的变化过程。 “反 应物”与“产物”泛指这些过程的初态和末态。 “优化”若未注明， 包括优化至极小点和优 化至过渡态。 势能面是高维的， 但为了直观以及表述方便， 文中一般用二维势能面模型来讨 论，应推广至高维情况。限于

3、纯文本格式，向量、矩阵无法加粗表示，但容易自行判断。目录：1. 过渡态2. 过渡态搜索算法2.1 基于初猜结构的算法2.1.1 牛顿 - 拉弗森法 (Newton-Raphson,NR) 与准牛顿法 (quasi-Newton,QN)2.1.2 AH 方法 (augmented Hessian)2.1.2.1 RFO 法 (Rational Function Optimization，有理函数优化 )2.1.2.2 P-RFO法 (Partitioned-RFO)2.1.2.3 QA 法 (Quadratic Approximation ，二次逼近 )2.1.2.4 TRIM 法 (trust

4、-region image minimization，置信区域镜像最小化 )2.1.2.5 在高斯中的常见问题2.1.3 GDIIS 法 (Geometry Direct Inversion in the Iterative Subspace)2.1.4 梯度模优化 (gradient norm minimization)2.1.5 Dimer 方法2.2 基于反应物与产物结构的算法2.2.1 同步转变方法 (synchronous transit,ST)2.2.2 STQN 方法 (Combined Synchronous Transit and Quasi-Newton Methods)2

5、.2.3 赝坐标法 (pseudo reaction coordinate)224 DHS方法(Dewar-Healy-Stewart ，亦称 Saddle 方法)与 LTP 方法(Line-Then-Plane)2.2.5 Ridge 方法2.2.6 Step-and-Slide 方法2.2.7 M ii ller -Brown 方法2.2.8 CI-NEB 、ANEBA方法2.3 基于反应物结构的算法2.3.1 最缓上升法 (least steep ascent,shallowest ascent)2.3.2 本征向量 / 本征 值跟踪法 (eigenvector/eigenvalue f

6、ollowing,EF。 也称 modewalking/mode following/Walking up valleys)2.3.3 ARTn(activation-relaxation technique nouveau)2.3.4 梯度极值法 (Gradient extremal,GE)2.3.5 约化梯度跟踪 (reduced gradient following,RGF)2.3.6 等势面搜索法 (Isopotenial Searching)2.3.7 球形优化 (Sphere optimization)2.4 全势能面扫描3. 过渡态相关问题3.1 无过渡态的反应途径 (barri

7、erless reaction pathways)3.2 Hammond-Leffler 假设3.2 对称性问题3.3 溶剂效应3.4 计算过渡态的建议流程4. 内禀反应坐标 (intrinsic reaction coordinate,IRC)5.IRC 算法5.1 最陡下降法 (Steepest descent)5.2 IMK 方法 (Ishida-Morokuma-Kormornicki)5.3 Mi ller -Brown 方法5.4 GS(Gonzalez-Schlegel) 方法6.chain-of-states 方法6.1 Drag method 方法6.2 PEB 方法 (pl

8、ain elastic band)6.3 Elber-Karplus 方法6.4 SPW 方法(Self-Penalty Walk)6.5 LUP 方法(Locally Updated planes)6.6 NEB 方法(Nudged Elastic Band)6.7 DNEB 方法 (Double Nudged Elastic Band)6.8 String 方法6.9 Simplified String 方法6.10 寻找过渡态的 chain-of-state 方法6.10.1 CI-NEB 方法6.10.2 ANEBA 方法 (adaptive nudged elastic band a

9、pproach)6.11 chain-of-states 方法的一些特点6.12 高斯中 opt 关键字的 path=M 方法6.13 CPK 方法 (Conjugate Peak Refinement)1. 过渡态过渡态结构指的是势能面上反应路径 上的能量最高点，它通过 最小能量路径 （minimum energy path,MEP） 连接着反应物和产物的结构（如果是多步反应的机理，则这里所指反应物 或产物包括中间体） 。对于多分子之间的反应， 更确切来讲过渡态结构连接的是它们由无穷 远接近后因为范德华力和静电力形成的复合物结构， 以及反应完毕但尚未无限远离时的复合 物结构。 确定过渡态有助

10、于了解反应机理， 以及通过势垒高度计算反应速率。 一般来讲， 势 垒小于 21kcal/mol 就可以在室温下发生。在势能面上，过渡态结构的能量对坐标的一阶导数为0，只有在反应坐标方向上曲率（对坐标二阶导数） 为负， 而其它方向上皆为正， 是能量面上的一阶鞍点。 过渡态结构的能量二阶 导数矩阵（ Hessian 矩阵）的本征值仅有一个负值，这个负值也就是过渡态拥有唯一虚频的 来源。 若将分子振动简化成谐振子模型， 这个负值便是频率公式中的力常数， 开根号后即得 虚数。分子构象转变、 化学反应过程中往往都有过渡态的存在， 即这个过程在势能面上的运动往往 都会经历满足上述条件的一点。 化学反应的过

11、渡态更确切应当成为“反应过渡态”。 需要注 意的是化学反应未必都经历过渡态结构。由于过渡态结构存在时间极短， 所以很难通过实验方法获得， 直到飞秒脉冲激光光谱的出现 才使检验反应机理为可能。 计算化学方法在目前是预测过渡态的最有力武器， 尽管计算上仍 有一些困难，比如其附近势能面相对于平衡结构更为平坦得多、低水平方法难以准确描述、 难以预测过渡态结构、缺乏绝对可靠的方法（如优化到能量极小点可用的最陡下降法）等。搜索过渡态的算法一般结合从头算、DFT方法，在半经验、或者小基组条件下，难以像描述平衡结构一样正确描述过渡态结构， 使得计算尺度受到了限制。 结合分子力场可以描述构象 变化的过渡态， 但

12、不适用描述反应过渡态， 因为大部分分子力场的势函数不允许分子拓扑结 构的改变，虽然也有一些力场如ReaxFF可以支持，有的力场还有对应的过渡态原子类型，但目前来看适用面仍然较窄，而且不够精确，尽管更为快速。注：严格来说， “过渡结构”是指势能面上反应路径上的能量最高点， 而“过渡态”是指自 由能面上反应路径上的能量最高点， 由于自由能变主要贡献自势能部分， 所以多数情况二者 结构近似一致。 但随着温度升高， 往往熵变的贡献导致自由能面与势能面形状发生明显偏离， 从而导致过渡结构与过渡态明显偏离， 两个词就不能混用了。 但本文不涉及相关问题， 故文 中过渡态、过渡结构一律指势能面上反应路径上的能

13、量最高点。2. 过渡态搜索算法2.1 基于初猜结构的算法2.1.1 牛顿 - 拉弗森法 (Newton-Raphson,NR) 与准牛顿法 (quasi-Newton,QN)NR法是寻找函数一阶导数为0 (驻点)位置的方法。通过对能量函数的泰勒级数的二阶近似展开，然后使用稳态条件 dE/dr=0 ，可导出步进公式： 下一步的坐标向量 = 当前坐标向量 - 能量一阶导数向量 * Hessian矩阵的逆矩阵。在势能面上以NR法最终找到的结果是与初猜位置 Hessian 矩阵本征值正负号一致、 离初猜结构最近的驻点， 由于能量极小点、 过渡态和 高阶鞍点的能量一阶导数皆为0,故都可以用NR法寻找。对

14、于纯二次形函数 NR法仅需一步即可找到正确位置，而势能面远比之复杂，所以需要反复 走步直至收敛。也因为势能面这个特点，为了改进优化，实际应用中NR法一般还结合线搜索步(line search)，对于优化至极小点，就是找当前点与NR法算出来的下一点的连线上的能量极小点作为实际下一步结构； 若优化至过渡态, 且连线方向主要指向过渡态, 则找的是 连线上能量极大点, 若主要指向其它方向则找连线的能量极小点, 若指向二者程度均等则一 般不做线搜索。 由于精确的线搜索很花时间, 所以一般只是在连线的当前位置附近计算几个 点的能量,以高阶多项式拟和后取其最小 / 最大点。NR法每一步需要计算 Hessia

15、n矩阵并且求其逆，所以十分昂贵。QN法与NR法的走步原理一 样,但 Hessian 矩阵最初是用低级或经验方法猜出来的, 每一步优化中通过当前及前一步的 梯度和坐标对 Hessian 矩阵逆矩阵逐渐修正。 由于只需计算一阶导数, 即便 Hessian 矩阵不 准确造成所需收敛步数增加，但一般仍比NR法速度快得多。QN法泛指基于此原理的一类方法,常用的是 BFGS(Broyden Fletcher Goldfarb Shanno) ,此法对 Hessian 的修正保持其 对 称 性 和 正 定 性 , 最 适 合 几 何 优 化 , 但 显 然 不 能 用 于 找 过 渡 态 。 还 有 DFP

16、(Davidon-Fletcher-Powell) , MS(Murtagh-Sargent ,亦称 symmetric rank 1, SR1),PSB(Powell-symmetric-Broyden)。也有混合方法，如 Bofill 法是 PSB和 MS法对 Hessian修正量的权重线性组合,比二者独立使用更优,权重系数通过位移、梯度改变量和当前 Hessian 计算得到,它对 Hessian 的修正不强制正定,很适宜搜索过渡态。将NR步进公式放到 Hessian本征向量空间下其意义更为明显(此时Hessian为对角矩阵)，可看出在每个方向上的位移就是这个方向势能的负梯度除以对应的本征

17、值,比如在 i 方向上的位移可写为 X(i)=-g(i)/h(i)，在受力越大、越平坦的方向位移越大。每一步实际位移就是这些方向上位移的矢量和。对于寻找过渡态,因为虚频方向对应Hessian 本征值为负,使位移为受力相反方向，所以NR法在过渡态附近每一步都是使虚频方向能量升高，而在其它正交的方向朝着能量降低的方向位移,通过这个原理步进到过渡态。若有n 个虚频,则NR法就在n个方向升高能量而其它方向降低能量找到n阶鞍点。由于 NR 法的这个特点,为找到正确类型的驻点,初猜结构必须在目标结构的二次区域 (quadratic region) 内。所谓的二次区域,是指驻点附近保持 Hessian 矩阵

18、本征值符号不变 的区域, 它的形状可以用多变量的二次函数近似描述, 例如二维势能面情况下这样的区域可 以用F(x,y)=A*xA2+B*yA2+C*x+D*y+E*x*y来近似描述。对于能量极小点，就是指初猜点在目标结构附近 Hessian 矩阵为正定矩阵的范围； 对于找过渡态, 就需要初猜点在它附近含有 且仅含有一个负本征值的范围内。 并且这个范围内不能有其它同类驻点比目标结构距离初猜 结构更近。NR法方便之处是只需要提供一个初猜结构即可，但是由于过渡态二次区域很小(相对于能量极小点来讲) ，复杂反应过渡态又不容易估计， 故对使用者的直觉和经验有一 定要求，即便是老手，也往往需要反复尝试。N

19、R 法对初猜结构比较敏感，离过渡态越近所需收敛步数越少，成功机率越高。模版法可以 帮助给出合理的初猜， 也就是如果已经知道其它机理相同的反应的过渡态结构， 可以保持反 应位点部分的结构不变而替换周围的原子，使之变成自己要研究的化合物反应的初猜结构。2.1.2 AH 方法 (Augmented Hessian)AH方法并不是独立的寻找过渡态的方法，而是通过修改原始 Hessian矩阵来调整NR法步进的长度和方向的一种方法。在NR法的步进公式中加入了一个移位参数入，式子变为 X(i,入)=-g(i)/(h(i)- 入)，NR法相当于 入等于0时的特例。入控制着每步步进距离， 它与h(i)的相对大小

20、也控制着这个坐标上的步进方向。根据设定入方法的不同，常见的有RFO P-RF0和QA/TRIM。这些方法每一步也使用QN方法来快速地更新 Hessian。下面提及的置信半径 R(Trust radius) 是指二阶泰勒级数展开这种近似的合理的区域，可以 在优化过程中固定也可以动态改变，比如下一步位置的实际能量与使用二阶泰勒级数展开预测的能量符合较好则加大R,反之减小。优化的每一步移动距离不应超过R,否则可能进入二阶泰勒级数展开近似的失效区域，NR法在势能面平坦的时候容易超过这个范围，应调整入避免。2.1.2.1 RFO 法(Rational Function Optimization，有理函数

21、优化 )对能量函数根据有理近似展开， 而不是NR法的二阶泰勒级数近似展开， 可推得与AH方法形 式相同的步进公式。确定其中 入的公式是 入=E ( g(i)A2/(h(i)-入)，g(i)和h(i)代表此方向的梯度和本征值，加和是对所有本征向量方向加和。通过迭代方法会解出N+1个入(N代表势能面维数)，将 入按大小排列，则有入(i) < h(i) w入(i+1)。故选其中最小的入可使各个方向位移公式的(h(i)-入)项皆为正，保证每步位移都向着极小点。选其中大于 m个Hessian本征值的入，将会在本征值最低的 m个方向上沿其上的受力反方向位移提升能 量，在其余N-m个方向上降低能量，由

22、此确保优化到m阶鞍点，若m为1即用来找过渡态。所以用了这个方法寻找指定类型驻点不再像NR法对初猜位置 Hessian本征值符号有要求，而是直接通过选择入来设定向着何种鞍点位移。如果每步步长度超过了R,则乘以一个小于1的因子来减小步长。值得一提的是，入与势能面维数 N的平方根近似成正比，随着体系尺度的增大，RFO的入对NR法的二次近似就会趋现“校正过度”情况，产生大小不一致 问题，可使用SIRFO(Size independent RFO)方法解决，即AH走步公式中的 入改为入/20.5。2.1.2.2 P-RFO 法(Partitio ned-RFO)专用于优化过渡态，效率比RFO更高。RFO

23、对所有方向的步进都使用同一个入，而在P-RFO中在指向过渡态的方向使用独立计算的入(TS),入(TS)=g(TS)A2/(h(TS)- 入(TS)，应选这个一元二次方程的最大的解，可保证在这个方向上升高能量。其余方向入的确定和RFO的公式一样，加和就不再包含指向过渡态的方向，并且选最小的 入解以使这些方向能量降低。这里所谓指向过渡态的方向一般是指最低本征值的方向，在上述RFO方法m为1时也是如此假设(限于其形式 RFO也只能用这最低模式)，但有时会是其它的非最低的模式，P-RFO也可以将这样的模式作为指向过渡态的模式，见后文EF方法的讨论。2.1.2.3 QA 法（Quadratic Appr

24、oximation，二次逼近）确定入的公式是( X(i)A2=刀(-g(i)/(h(i)- 入)A2=RA2，也就是说每一步移动的距离 恰好是置信半径，这样步进速度较快。若优化到过渡态，计算入公式的加和中指向过渡态本征向量的那一项的 入改为-入，即 X(TS)=-g(TS)/(h(TS)+ 入)。2.1.2.4 TRIM 法 (trust-region image minimization，置信区域镜像最小化 )这个方法假设 Hessian 本征值最小的方向的梯度和曲率符号与原本相反，而其它方向不变。 经过这样的变化后原来的过渡态位置就成为了能量极小点 ( 过渡态的 image) ，这样就可以

25、通 过优化到极小点而得到过渡态。将TRIM的假设g(TS)'=-g(TS) , h(TS)'=-h(TS)代入AH方法的步进公式 X(i,入)=-g(i)/(h(i)- 入)，再使分子分母同乘以-1，可知在过渡态方向上 的步进公式与其它方向区别仅在于反转了入的符号。又由于TRIM也是通过调整入使步进长度等于为置信半径，所以在公式的形式上与QA法找过渡态的公式完全一致，QA与TRIM可互为同义词。通过如上调整 AH方法引入的 入可使NR法的步进更有效率、更稳定，还可以通过它改变步 进公式在不同方向上的分母项符号， 使优化过渡态的初猜点不限于过渡态的二次区域。 可直 接指定沿某个振

26、动模式升高能量来找过渡态， 即便当前点这个方向的Hessian 本征值可能是正值，例如从极小点开始跟踪至过渡态，见后文的 EF方法。2.1.2.5 在高斯中的常见问题高斯中 opt=ts 是使用 Berny 算法来找过渡态，需要提供一个初猜结构。 Berny 默认的走步 的方法是RFO/P-RFO分别对于优化至极小值/鞍点)，若加了 Newton选项，则走步基于 NR 法。每一步对 Hessian 矩阵的更新方法以 UpdateMethod 选项指定，寻找极小点时默认用 BFGS， 找过渡态时默认用 Bofill o Berny算法还包括一些细节步骤在内，比如投影掉被冻结的变量、更新置信半径、

27、设定了线搜索过程中几种方案等等，详见手册opt 关键字。使用了每步修正 Hessian 的准牛顿法后，初猜的 Hessian 矩阵质量明显影响结构收敛速度， 它的不准确容易导致搜索过渡态失败(在高斯中默认使用价键力场得到Hessian )。这种情况需要昂贵的 calcfc 关键字以当前方法水平计算最初的 Hessian 矩阵，若使用的方法在程 序中支持解析二阶导数， 速度会较好。 或者用 readfc 来读取包含了 Hessian 矩阵信息的 chk 文件，可以先使用低水平方法进行简正振动分析得到chk 文件，再将之读入作为 Hessian矩阵初猜，能够节约时间，但前提是此势能面对方法等级不敏

28、感(一般如此)。使用了更准 确的初猜后不仅可以增加找到过渡态的成功机率， 还有助于在更短的优化步数内达到收敛标 准。若使用 calcall ，则每一点都重新准确计算 Hession ，会更为可靠，但极为昂贵。高斯中 berny 方法寻找过渡态默认每步会检查 Hessian 矩阵的本征值是否仅有一个为负， 如 果不符，就会提示“ a wrong sign eigenvalue in hessian matrix ”，经常一开始就报错， 原因是初猜结构不符合这个条件， 即便这个初猜通过 berny 方法最终能够正确优化到过渡态， 这时应加 noeigentest 选项避免本征值符号的检查， 不符合

29、要求也继续优化，但因此可能收敛到其它类型驻点。有时这种情况由初猜的 Hessian 不准导致，可用 calcfc 解决。如果搜 索的过渡态出现多个负本征值， 可根据适当的虚频 (高斯中以负数频率表示)振动方向调整结构以降低能量，直至剩下一个虚频，再重新优化。高斯中默认的置信半径为 0.3 bohr,若优化中步长(RFO/P-RFO步)超过就会输出“ Maximumstep size (0.300) exceeded in Quadratic search ”和“ Step size scaled by xxx”， 即乘以 xxx 调小步长至置信半径内。 也可以使用 iop(1/8=k) 将置信

30、半径改为 k*0.01 bohr(1 bohr=0.5292 埃) ，调大后往往可以显著减少收敛步数，很适合势能面平坦的大体系。注意 并不是每一步的步长都固定为 k*0.01 bohr ，若没超过置信半径则步长并不因此改变。寻找 极小点时默认为允许动态改变置信半径，此时 iop(1/8) 设的就是最初的置信半径，对于寻 找过渡态默认为关闭此功能(相当于用了 NoTrustUpdate )，可以使用 trustupdate 关键字 来打开这个功能。2.1.3 GDIIS 法 (Geometry Direct Inversion in the Iterative Subspace)GDIIS 与

31、DIIS 原理一致，但用于几何优化，这个方法趋于收敛到离初始位置最近的驻点， 包括过渡态。下一步坐标 X(new)=X"-H'g" ， H' 代表当前步的 Hessian 逆矩阵，可见公式形 式与NR法是一致的，但是X"与g"不再指当前步的坐标和梯度，而是由之前走过的点的坐标 X(k)和梯度g(k)插值得到的，X"=Ec(k)X(k), g"=刀c(k)g(k),代入上式即X(new)=c(i)(X(i)-H'g(i)，其中刀是对之前全部走过的点加和。系数c(k)通过使误差向量r的模最小化得到，r=刀c(k)e

32、(k),并以刀c(k)=1为限制条件。e(k)常见有两种定义， 一种是 e(k)=-H'g(k) ，另一种更常用， 是 e(k)=g(k) ，可看出 GDIIS 利用的是已经搜索过的 子空间中坐标与梯度的相关性， 目的是估出梯度 (即误差向量) 的模尽可能小的坐标， 这一 点与后述的梯度模方法相似。此方法缺点是由于势能面复杂， 步进中容易被拉到已经过的势能面的其它驻点而不能到达指 定类型驻点，还容易走到类似肩膀形状的拐点，梯度虽小却不为0，由于不能达到收敛标准而反复在此处震荡。 另外随优化步数增加， 误差向量数目逐渐加大， 会逐渐出现的误差向量 之间的线性相关，导致伪收敛和数值不稳定问

33、题。在改进的方法中将GDIIS 与更可靠的 RFO方法结合，比较二者的步进方向和长度，并检验 GDIIS中的组合系数c,根据一定规则来决 定每一步对之前走过的点的保留方式， 必要时全部舍去而重新开始 GDIIS。 Gaussian 中用的 这种改进的GDIIS方法解决了上述问题同时提高了效率，速度等于或优于 RFO方法，尤其是以低水平对势能面平坦的大体系优化时更为突出。 GDIIS 计算量小, 对 Hessian 矩阵很不敏 感，可以在优化中不更新，也可以用QN法更新来改善性能。此方法自Gaussian98起就是默认的半经验优化算法，其它方法下也可以用OPT的gdiis关键词打开。2.1.4

34、梯度模优化 (gradient norm minimization) 势能面上的驻点,包括能量极小点、过渡态和高阶鞍点的势能梯度都为0,所以在相应于势能面的梯度模面上进行优化找到数值为 0 的点, 经过 Hessian 矩阵本征值符号的检验, 就能 得到过渡态。 这相当于把搜索过渡态问题转化为了能量极小化问题,就有了更可靠的算法可用。(注：梯度模指的是势能梯度在各个维度分量平方和的平方根,即梯度大小的绝对值) 。但是寻找数值更小点的优化方法比如最陡下降法只能找到离初始位置最近的极小点,若找到的梯度模面上的极小点数值大于 0 则是势能面肩膀形拐点, 没有什么用处, 而这样的点收敛 半径往往很大，

35、例如图中在x=2至8的区域内都会收敛到函数拐点，只有提供的初猜结构在x=1 和 9 附近很小的范围内才会收敛到过渡态,收敛半径太小,难以提供合理初猜。梯度模 面上还多出一些极大点，如x=1.5处，若使用收敛更快的 NR法找极小点还容易收敛到这样没有意义的点上。基于这些原因,梯度模法很少使用。2嚅>uo 一UUP 亠2.1.5 Dimer 方法Dimer方法是一种高效的定位过渡态的方法。这个方法定义了由两个点 R1和R2组成的一个Dimer，能量和所受势能力(由原始的势能面梯度造成受力，下同)分别为E1和E2、F1和F2。两个点间距为 2A R, A R为定值。这两点的中间点为 R,其受力

36、F(R)=(F1+F2)/2。Dimer 的总能量为E=(E1+E2)/2。这个方法的每一步包括平移 Dimer和旋转Dimer两步。旋转Dimer :保持R1、R2中点位置R不变作为轴，旋转 Dimer直到总能量E最小。通过推 导可知在旋转过程中， E与R点在dimer方向(R1-R2方向)上的曲率关系 C是线性的，即 最小化E的过程就是最小化 C的过程。所以每一步的Dimer方向都是曲率最小方向， 当最终 R收敛到过渡态位置时， Dimer就会平行于虚频方向。平移Dimer: Dimer根据受力F'移动R的位置，结合不同方法有具体步进方式，如quick-win、共轭梯度法。当C&l

37、t;0(过渡态或高阶鞍点的二次区域内)，F'等于将F(R)平行于Dimer方向力的分量符号反转；当C>0(极小点二次区域内)，F'等于F(R)平行于Dimer方向力的分量 的负值，而没有垂直于 Dimer方向的力，促使Dimer尽快离开这个区域。 由于Dimer的方向 就是曲率最小的方向，在过渡态二次区域内就是指虚频方向，在Dimer方法中F'的定义使这个方向以受力相反方向移动以升高能量，而其它方向顺着受力方向移动来最小化能量，可看出原理上与NR法相似。费时的计算 Hessian矩阵最小本征值以确定提升能量方向的过程 被旋转Dimer这一步代替了，仅需要计算一阶导

38、数。Dimer法对初始位置要求很宽松，并不需要在过渡态二次区域内，若在极小点二次区域内就类似于后述的EF方法沿着最小振动模式爬坡。如果在高阶鞍点二次区域内，只在曲率最负的虚频方向沿着受力反向移动，在其它虚频方向上仍最小化能量，而不会像NR法收敛到高阶鞍点。图2右侧为Dimer法在Mu Iler -Brown模型势上面搜索两个过渡态过程中Dimer走的路径。势能面上往往有许多鞍点，Dimer方法还可以做鞍点搜索。通过分子动力学方法给予Dimer一定动能，使之能够在势能面上广阔的区域内运动，根据一定标准提取轨迹中的一些点作为初猜，再执行标准Dimer方法就可以得到许多不同的鞍点。 Dimer方法很

39、适合双处理器并行， 两个点的受力分别由两个处理器负责，速度可增加将近一倍。2.2基于反应物与产物结构的算法2.2.1 同步转变方法（synchronous transit,ST）提供合理的初猜结构往往不易，ST方法可以只根据反应物和产物结构自动得到过渡态结构。“同步转变”这个名字强调的是反应路径上所有坐标一起变化，这是相对于后面提到的赝坐标法来说的（即只变化指定的坐标，尽管其它坐标优化后坐标也会变化）。ST分为两种模型，最简单的就是 LST模型（Linear synchronous transit,线性同步转变），这个方法假设反应过程中， 反应物结构的每个坐标都是同步、 线性地变化到产物结构。

40、 如果 反应物、产物的坐标分别以向量 A、B表示，则反应过程中的结构坐标可表示为（1-x）*A+x*B ，x由0逐渐变到1代表反应进度。注意 LST并不是指反应中原子在真实空间上以直线运动， 只有笛卡尔坐标下的 LST才是如此，在内坐标下的LST，原子在真实空间中一般以弧线运动。以LST的假设，反应路径在其所用坐标下的势能面图上可描述为一条直线，LST给出的过渡态就是这条直线上能量最高点（图 3的点1 ）。LST的问题也很显著，其假设的坐标线性变 化多数是错误的，绘制在势能面图上也多数不会是直线，故给出的过渡态也有较大偏差，容易带两个或多个虚频。比LST更合理的是QST（quadratic s

41、ynchronous transit,二次同步转变），它假设反应路径在势能面上是一条二次曲线。QST在 LST得到的过渡态位置上，对LST直线路径的垂直方向进行线搜索找到能量极小点A （图3的点2）。QST给出的反应路径可以用经过反应物、A产物的二次曲线来表示，如果这条路径上能量最大点的位置恰为A,则A就是QST方法给出的过渡态；如果不是，则以最大点作为过渡态。若想结果更精确，可以再对这个最大点向垂直于路径的方向优化，再次得到A并检验，反复重复这个步骤，逐步找到能量更低、 更准确的过渡态。QST方法在计算能力较低的年代曾是简单快速的获得过渡态和反应路径的方法，然而如今看来其结果是相当粗糙的，已

42、极少单独使用，可以将其得到的过渡态作为AH法的初猜。图3LST与QST方法示意图2.2.2 STQN 方法（Combi ned Syn chro nous Tran sit and Quasi-Newton Methods）STQN是 ST与QN方法的结合（更准确地说是与EF法的结合）。但不要简单认为是按顺序独立执行这两步，即认为“先利用反应物和产物结构以ST方法得到粗糙过渡态，再以之作为初猜用QN法精确寻找过渡态”是错误的。STQh方法大意是：使结构从低能量的反应物出发， 以ST路径在当前位置切线为引导，沿着LST或QST假设的反应路径行进（爬坡步），目的是使结构到达假设路径的能量最高处附近

43、（真实过渡态二次区域附近）。当符合一定判据时就转换为QN法寻找精确过渡态位置（EF步）。下面介绍具体步骤。先说明后面用到的切线的定义：STQN当中的LST路径与前面ST部分介绍的LST路径无异，都是直线，切线 T在优化中是不变的，就是反应物R指向产物P的单位向量。STQN方法中的QST路径定义与ST方法介绍的不同，走的不是二次曲线而是圆形的一段弧，如图4所示。这个圆弧经过 R、P以及优化中的当前步位置 X,切线就是圆在 X处的单位切线向量，圆弧 和切线在每一步都是变化的。虽然QST路径比LST更为合理，但对于 QSTN方法，QST路径在收敛速度和成功机率上的优势并不显著。图4STQN对QST路

44、径的定义STQN每一步执行内容如下：(1)首先重新计算或用 QN法更新Hessian。 (2)按上述方法计算 当前位置处的切线。(3)决定这一步是爬坡步还是EF步。如果是优化的第前两步，则一定认为是爬坡步，因为此时离过渡态区域还较远，应当先爬坡。如果是第3、4步，则估算出在切线方向的位移，超过一定标准就是爬坡步，否则说明爬得差不多了就进入EF步找过渡态。如果是第5步之后，一般已离过渡态区域较近，故一定认为是EF步。(4)如果是爬坡步，则在切线方向上移动(将切线方向作为EF方法所跟踪的振动方向来计算位移大小)。如果是EF步，首先计算 Hessian各个本征向量的与切线重叠情况，如果有重叠大于0.

45、8的本征向量，则以EF法跟踪本征值最大的本征向量来移动，相当于继续向上爬。如果没有大于0.8的，就跟踪最小本征值的本征向量移动来寻找过渡态。(5)步长长度若大于标准则调小，默认0.3 bohr 。 (6)根据预置受力、位移标准判断是否已收敛，收敛则结束循环。注意，ST方法中具体包含 LST和QST两种方法，STQN也用到了 LST和QST两种反应路径的 假设。高斯中的LST方法指的是ST中的LST方法，而QST2/3指的是利用 QST路径假设的 STQN方法，它们原理上截然不同，不要混淆。高斯中的QST2只需输入反应物和产物结构，通过几何方法估出 STQN的初始步结构X。QST3需额外输入猜测

46、的过渡态，它直接作为X，一般比QST2效果更好。对于经验不足的用户，用STQN方法往往比只提供过渡态初猜的方法 更为适合。注意产物和反应物应当使用同样方法同样基组进行优化，如果是多分子比如 A+B=C+D这样的反应，应当优化 A和B/C和D的复合物作为输入的产物 /产物，而不是单独 优化A、B然后拼到一起，因为形成范德华复合物后孤立的分子会有一定构象改变，能量也 低于它们孤立状态的加和。2.2.3 赝坐标法(pseudo reaction coord in ate) 也称为坐标驱动法(Coordinate Driving)。这个方法在高斯中就是柔性扫描(Relaxed Scan)，即扫描一个变

47、量，但每一步对其它变量自动进行优化，每一步得到的结构就是在这个变量为一定值情况下的最优结构。赝坐标法扫描的是反应物转变到产物过程中的关键坐标，比如扫描化学键断裂/生成反应中的键长。扫描的结果就是近似的IRC,可以再将能量最高点作为初猜找过渡态，或者用更小的步长再次扫描能量最高点附近找更精确的过渡态结构。这个找 过渡态方法实际上用的是能量极小化优化过程， 由于这样的算法比寻找过渡态的算法更为稳 健，所以赝坐标法是颇可靠的，其它方法失败时可考虑这种方法。这个方法缺点是费时间， 而且不适合通向过渡态路径中反应区域涉及多个坐标变化的反应过 程，因为自定义扫描的内容很难全面、 准确考虑到这些坐标变量的变

48、化， 结果难以说明问题， 没有考虑进去的关键变量容易产生滞后效应 (hysteresis effect) 。比如乙烷由交叉构象变 化到另一个交叉构象，需要经历重叠构象的过渡态，会涉及到三个HCCH二面角同时由 60度变化到0度，如果用赝坐标法只扫描其中一个HCCF由60度变到0度，则每一步其它两个HCCH角 一定会大于这个扫描的二面角，与实际不符。这是因为这两个角越小，分子的能量 越高， 每一步自动优化的时候它们更倾向于保持在大角度。 最终到达过渡态时， 所扫描的二 面角到达了 0 度， 另外两个二面角却大于 0 度，说明它们的运动比实际的过程滞后了。 由于 滞后效应， 从反应物和产物两个方向

49、扫描同相同的坐标， 得到的路径也不同。 上述简单的反 应此方法滞后效应尚且严重， 对于复杂变化， 这种效应导致的问题更难以预测。 故此方法确 定的IRC、过渡态不可靠，只建议对简单的反应使用这种方法，扫描变量的选择注意避免滞 后效应。在高斯中此方法可以使用 opt=modredundant 或 Opt=Z-matrix 结合分子结构部分标记的扫描 变量来实现。例如使用 opt=modredundant 并在分子结构末尾写上 A 3 2 1 S 10 1.000000 来指定 3 2 1 原子组成的角度进行柔性扫描，共 10 步，每步 1.0 度。如果不熟悉，也可以 很方便地在 GaussVie

50、w 里的冗余坐标编辑器里面添加要柔性扫描的变量。如果只执行常规的某个变量的扫描，比如高斯中的 scan 来找能量最高点作为初猜结构，对 于简单体系可行， 但对于复杂体系， 这样忽略了此变量的变化导致分子其它部分结构的驰豫， 如此得到的能量最高点作为过渡态初猜很不可靠， 因为势能中掺入了不合理的结构造成的能 量升高，使势能曲线形状改变。2.2.4 DHS 方法 (Dewar-Healy-Stewart ，亦称 Saddle 方法 )与 LTP 方法 (Line-Then-Plane)DHS方法中第一步将反应和产物分别作为A点和B点，确定哪个点能量低，比如A比B低就把A点的结构向B点稍微做调整(5

51、%)得到A，然后限制变量空间中 A'与B的距离不变(即 在超球面上)对 A'进行优化得到 A''。将A''与B当作下一步的起始点 A与B,重复上述方 法。这样反复进行迭代， 若以序号 n 代表第 n 次得到的 A'' 或 B'' ，会依次得到例如 A''(1)、A''(2) 、B''(1) 、A''(3)直到A与B十分接近时才停止迭代，此位置就是过渡态。将得到的全部 A''(n) 按序号 n 依次连接， B''(n) 也

52、按序号依次连接，再将序号最大的 A''(n) 与B”(n)连接，得到的就是近似的IRC。LTP与DHS方法基本一致，不同的是每步是在垂直于A'与B连线的超平面上优化。DHS方法虽然可以很快地走到过渡态附近的位置，但是越往后每步的 AB 距离缩近也越少，故并不能有效率地贴近过渡态。然而每步的在连线上调整的 距离不可过大，否则可能造成一侧的点跨过过渡态势垒跑到另一侧得到错误结果。225 Ridge 方法第一步时将反应物、 产物作为A点和B点，在其LST的路径上找到能量最大点C,然后在AC与BC直线上相距C为s的位置上分别设一点 A和B'，将A与B'分别沿着此

53、处势能面负梯 度优化p距离，将得到的 A"与B"作为下一步的 A和B。反复进行这个步骤，收敛后 C的 位置就是过渡态位置。s和p是计算过程中动态调节的参数，对结果影响较大，它们应当随 C逐渐接近过渡态而减小，可设若当前步的C能量高于上一步的 C,则减小p至原先一半；若s与p的比值大于某个数值，s也减半。Ridge方法的缺点是接近过渡态时效率较低，可 以当C进入过渡态二次区域后改用QN法来加快收敛。也可以结合DIIS法，速度比原先有一半以上的提升，效率有时还高于基于二阶导数的方法，而且在某些势能面非常平坦的体系比二阶导数方法更可靠。图6Ridge方法示意图226 Step-a

54、 nd-Slide方法使产物和反应物的结构同时顺着LST描述的路径相对移动（step步），直到它们的能量都等于某个预先设定的能量，然后让这两个结构在它们当前所在的势能等值面上滑动（slide 步），使二者结构在坐标空间中的距离最小。重复上述step和slide步骤，最终当两个结构碰上，这个位置就是过渡态。图7Step and Slide 方法示意图2.2.7 M ii ller -Brown 方法见下文IRC算法相应部分2.2.8 CI-NEB、ANEBA方法见下文"寻找过渡态的chai n-o f-state方法”相应部分2.3基于反应物结构的算法2.3.1 最缓上升法 （leas

55、t steep ascent,shallowest ascent）由反应物结构到达过渡态结构的过程是沿着势能面最容易行进的路径进行的（不考虑动力学问题），这个途径一般比其它方向要缓和，所以由反应物结构开始，沿着势能面最缓的方向逐渐往上爬，往往可以沿着MEP到达过渡态。但要注意这条路径时常与从过渡态沿最陡下降 路径所走出的 MEP并不一致，因此原理上此法不能保证一定能到达过渡态。图8描述的是LEPS势结合谐振势的势能面，最缓上升法所走的黑色粗曲线严重不符合实际MEFP黑点所示路径），而且曲线是中断的。 此法也可能走到与此平衡结构相连的其它过渡态，而非预期的过渡态。还容易因为步长问题导致走到中途时

56、跑到另外一条错误路径上，虽然设小步长能得到解决，但是需要花费更长时间。因为种种问题，这个方法使用较少。图8势能面上最缓上升法所走的路径（黑色粗曲线）2.3.2 本征向量 / 本征值跟踪法（eigenvector/eigenvaluefollowing,EF。也称 modewalki ng/mode followi ng/Walki ng up valleys）由于平衡结构越过势垒发生反应的能量主要来自分子某振动模式提供的动能，考虑这一点，由平衡结构沿着此振动矢量方向步进，能够找到过渡态，经历的路径就是反应路径。这种方法需要首先对平衡结构进行振动分析，由用户最初指定一个可能指向过渡态的振动模式。

57、因为平衡态通向过渡态路径势能面平缓，曲率（可视为振子力常数）一般小于其它方向，故一般跟踪频率最低的振动模式（高斯中默认）。每走一步后重新计算 Hessian矩阵的本征值和 本征向量，如果跟踪的是本征值最低的模式，仍取本征值最小的本征向量继续跟踪；如果跟踪的是其它振动模式，就取与上一步所跟踪的向量重叠最大的向量继续跟踪。重复执行，直到符合收敛标准为止。如果一个结构涉及到多个过渡态，则跟踪不同的本征向量有可能得到不同的过渡态，即便所跟踪的不是最低模式，当接近过渡态后也会成为最低的模式。此方法也可以直接由过渡态初猜结构开始跟踪，或者说 EF方法是一种不需要初猜在过渡态二次区域内的寻找过渡态的方 法。由稳定结构通过 EF方法跟踪至过渡态相对与直接给出初猜显然更为费时，但对于不能 预测过渡态结构的情况下往往是有用的。LMOD法搜索构象也是基于这一原理，不断地根据低频振动方向越过构象转变的过渡态到达新的构象。opt=(EF,TS)方法最初的EF方法只是简单地沿所跟踪的振动模式移动来升高能量。高斯中还使结构同时在其余方向上沿能量更低的方向移动，其实它用的就已介绍的 P-RFO法，所跟踪的模式用独立计算的入的最大解，其它的模式使用相同的另外计算的入的最小解。由于Berny方法寻找过渡态已经包含了P-RFO步，所以EF方法实际上也已经包含在内了，除非要用到跟踪特定