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

过渡态、反应路径的计算方法及相关问题SoberevaDepartment of Chemistry,University of Science and Technology Beijing,Beijing 100083,China 前言：本文主要介绍过渡态、反应路径的计算方法，并讨论相关问题。由于这类算法极多，可以互相组合，限于精力不可能面面俱到展开，所以只介绍常用，或者实用价值有限但有启发性的方法。文中图片来自相关文献，做了一定修改。由于本文作为帖子发布，文中无法插入复杂公式，故文中尽量将公式转化为文字描述并加以解释，这样必然不如公式形式严谨，而且过于复杂的公式只能略过，但我想这样做的好处是更易把握方法的梗概，有兴趣可以进一步阅读原文了解细节。对于 Gaussian 中可以实现的方法，文中对其在Gaussian 中的使用进行了一些讨论，希望能纠正一些网上流传的误区。虽然绝大多数人不专门研究计算方法，其中很多方法也不会用到，但多了解一下对开阔思路是很有好处的。文中指的“反应”包括构象变化、异构化、单分子反应等任何涉及到过渡态的变化过程。“反应物”与“产物”泛指这些过程的初态和末态。“优化”若未注明，包括优化至极小点和优化至过渡态。势能面是高维的，但为了直观以及表述方便，文中一般用二维势能面模型来讨论，应推广至高维情况。限于纯文本格式，向量、矩阵无法加粗表示，但容易自行判断。目录：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-region image minimization，置信区域镜像最小化)2.1.2.5 在高斯中的常见问题 2.1.3 GDIIS法(Geometry Direct Inversion in the Iterative Subspace)名师资料总结-精品资料欢迎下载-名师精心整理-第 1 页，共 30 页 -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.2.3 赝坐标法(pseudo reaction coordinate)2.2.4 DHS方法(Dewar-Healy-Stewart，亦称 Saddle 方法)与 LTP方法(Line-Then-Plane)2.2.5 Ridge方法 2.2.6 Step-and-Slide方法2.2.7 M ller-Brown 方法 2.2.8 CI-NEB、ANEBA 方法2.3 基于反应物结构的算法 2.3.1 最缓上升法(least steep ascent,shallowest ascent)2.3.2 本征向量/本征 值跟踪法(eigenvector/eigenvalue following,EF。也称mode walking/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 无过渡态的反应途径(barrierless 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 Mller-Brown 方法5.4 GS(Gonzalez-Schlegel)方法6.chain-of-states方法6.1 Drag method方法6.2 PEB 方法(plain 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方法名师资料总结-精品资料欢迎下载-名师精心整理-第 2 页，共 30 页 -6.10.1 CI-NEB方法 6.10.2 ANEBA方法(adaptive nudged elastic band approach)6.11 chain-of-states方法的一些特点6.12 高斯中 opt 关键字的path=M 方法6.13 CPK 方法(Conjugate Peak Refinement)1.过渡态过渡态结构指的是势能面上反应路径上的能量最高点，它通过最小能量路径(minimum energy path,MEP)连接着反应物和产物的结构（如果是多步反应的机理，则这里所指反应物或产物包括中间体）。对于多分子之间的反应，更确切来讲过渡态结构连接的是它们由无穷远接近后因为范德华力和静电力形成的复合物结构，以及反应完毕但尚未无限远离时的复合物结构。确定过渡态有助于了解反应机理，以及通过势垒高度计算反应速率。一般来讲，势垒小于 21kcal/mol就可以在室温下发生。在势能面上，过渡态结构的能量对坐标的一阶导数为0，只有在反应坐标方向上曲率（对坐标二阶导数）为负，而其它方向上皆为正，是能量面上的一阶鞍点。过渡态结构的能量二阶导数矩阵（Hessian 矩阵）的本征值仅有一个负值，这个负值也就是过渡态拥有唯一虚频的来源。若将分子振动简化成谐振子模型，这个负值便是频率公式中的力常数，开根号后即得虚数。分子构象转变、化学反应过程中往往都有过渡态的存在，即这个过程在势能面上的运动往往都会经历满足上述条件的一点。化学反应的过渡态更确切应当成为“反应过渡态”。需要注意的是化学反应未必都经历过渡态结构。由于过渡态结构存在时间极短，所以很难通过实验方法获得，直到飞秒脉冲激光光谱的出现才使检验反应机理为可能。计算化学方法在目前是预测过渡态的最有力武器，尽管计算上仍有一些困难，比如其附近势能面相对于平衡结构更为平坦得多、低水平方法难以准确描述、难以预测过渡态结构、缺乏绝对可靠的方法（如优化到能量极小点可用的最陡下降法）等。搜索过渡态的算法一般结合从头算、DFT方法，在半经验、或者小基组条件下，难以像描述平衡结构一样正确描述过渡态结构，使得计算尺度受到了限制。结合分子力场可以描述构象变化的过渡态，但不适用描述反应过渡态，因为大部分分子力场的势函数不允许分子拓扑结构的改变，虽然也有一些力场如ReaxFF可以支持，有的力场还有对应的过渡态原子类型，但目前来看适用面仍然较窄，而且不够精确，尽管更为快速。注：严格来说，“过渡结构”是指势能面上反应路径上的能量最高点，而“过渡态”是指自由能面上反应路径上的能量最高点，由于自由能变主要贡献自势能部分，所以多数情况二者结构近似一致。但随着温度升高，往往熵变的贡献导致自由能面与势能面形状发生明显偏离，从而导致过渡结构与过渡态明显偏离，两个词就不能混用了。但本文不涉及相关问题，故文中过渡态、过渡结构一律指势能面上反应路径上的能量最高点。2.过渡态搜索算法名师资料总结-精品资料欢迎下载-名师精心整理-第 3 页，共 30 页 -2.1 基于初猜结构的算法2.1.1 牛顿-拉弗森法(Newton-Raphson,NR)与准牛顿法(quasi-Newton,QN)NR法是寻找函数一阶导数为0（驻点）位置的方法。通过对能量函数的泰勒级数的二阶近似展开，然后使用稳态条件dE/dr=0，可导出步进公式：下一步的坐标向量=当前坐标向量-能量一阶导数向量*Hessian 矩阵的逆矩阵。在势能面上以NR法最终找到的结果是与初猜位置 Hessian 矩阵本征值正负号一致、离初猜结构最近的驻点，由于能量极小点、过渡态和高阶鞍点的能量一阶导数皆为0，故都可以用NR法寻找。对于纯二次形函数NR法仅需一步即可找到正确位置，而势能面远比之复杂，所以需要反复走步直至收敛。也因为势能面这个特点，为了改进优化，实际应用中NR法一般还结合线搜索步(line search)，对于优化至极小点，就是找当前点与NR法算出来的下一点的连线上的能量极小点作为实际下一步结构；若优化至过渡态，且连线方向主要指向过渡态，则找的是连线上能量极大点，若主要指向其它方向则找连线的能量极小点，若指向二者程度均等则一般不做线搜索。由于精确的线搜索很花时间，所以一般只是在连线的当前位置附近计算几个点的能量，以高阶多项式拟和后取其最小/最大点。NR法每一步需要计算Hessian 矩阵并且求其逆，所以十分昂贵。QN法与 NR法的走步原理一样，但 Hessian 矩阵最初是用低级或经验方法猜出来的，每一步优化中通过当前及前一步的梯度和坐标对Hessian 矩阵逆矩阵逐渐修正。由于只需计算一阶导数，即便 Hessian 矩阵不准确造成所需收敛步数增加，但一般仍比NR法速度快得多。QN法泛指基于此原理的一类方法，常用的是BFGS(Broyden Fletcher Goldfarb Shanno)，此法对Hessian 的修正保持其对 称 性 和 正 定 性，最 适 合 几 何 优 化，但 显 然 不 能 用 于 找 过 渡 态。还 有DFP(Davidon-Fletcher-Powell)，MS(Murtagh-Sargent，亦称symmetric rank 1，SR1)，PSB(Powell-symmetric-Broyden)。也有混合方法，如Bofill法是 PSB和 MS法对 Hessian修正量的权重线性组合，比二者独立使用更优，权重系数通过位移、梯度改变量和当前Hessian 计算得到，它对Hessian 的修正不强制正定，很适宜搜索过渡态。将 NR步进公式放到Hessian 本征向量空间下其意义更为明显（此时 Hessian 为对角矩阵），可看出在每个方向上的位移就是这个方向势能的负梯度除以对应的本征值，比如在 i 方向上的位移可写为X(i)=-g(i)/h(i)，在受力越大、越平坦的方向位移越大。每一步实际位移就是这些方向上位移的矢量和。对于寻找过渡态，因为虚频方向对应Hessian 本征值为负，使位移为受力相反方向，所以NR法在过渡态附近每一步都是使虚频方向能量升高，而在其它正交的方向朝着能量降低的方向位移，通过这个原理步进到过渡态。若有n 个虚频，则NR法就在 n 个方向升高能量而其它方向降低能量找到n 阶鞍点。由于NR 法的这个特点，为找到正确类型的驻点，初猜结构必须在目标结构的二次区域(quadratic region)内。所谓的二次区域，是指驻点附近保持Hessian 矩阵本征值符号不变的区域，它的形状可以用多变量的二次函数近似描述，例如二维势能面情况下这样的区域可以用 F(x,y)=A*x2+B*y2+C*x+D*y+E*x*y来近似描述。对于能量极小点，就是指初猜点在目标结构附近Hessian 矩阵为正定矩阵的范围；对于找过渡态，就需要初猜点在它附近含有且仅含有一个负本征值的范围内。并且这个范围内不能有其它同类驻点比目标结构距离初猜结构更近。NR 法方便之处是只需要提供一个初猜结构即可，但是由于过渡态二次区域很小名师资料总结-精品资料欢迎下载-名师精心整理-第 4 页，共 30 页 -（相对于能量极小点来讲），复杂反应过渡态又不容易估计，故对使用者的直觉和经验有一定要求，即便是老手，也往往需要反复尝试。NR 法对初猜结构比较敏感，离过渡态越近所需收敛步数越少，成功机率越高。模版法可以帮助给出合理的初猜，也就是如果已经知道其它机理相同的反应的过渡态结构，可以保持反应位点部分的结构不变而替换周围的原子，使之变成自己要研究的化合物反应的初猜结构。2.1.2 AH方法(Augmented Hessian)AH方法并不是独立的寻找过渡态的方法，而是通过修改原始Hessian 矩阵来调整NR法步进的长度和方向的一种方法。在NR 法的步进公式中加入了一个移位参数，式子变为X(i,)=-g(i)/(h(i)-)，NR法相当于 等于 0 时的特例。控制着每步步进距离，它与 h(i)的相对大小也控制着这个坐标上的步进方向。根据设定 方法的不同，常见的有RFO、P-RFO和 QA/TRIM。这些方法每一步也使用QN方法来快速地更新Hessian。下面提及的置信半径R(Trust radius)是指二阶泰勒级数展开这种近似的合理的区域，可以在优化过程中固定也可以动态改变，比如下一步位置的实际能量与使用二阶泰勒级数展开预测的能量符合较好则加大R，反之减小。优化的每一步移动距离不应超过R，否则可能进入二阶泰勒级数展开近似的失效区域，NR 法在势能面平坦的时候容易超过这个范围，应调整 避免。2.1.2.1 RFO法(Rational Function Optimization，有理函数优化)对能量函数根据有理近似展开，而不是 NR法的二阶泰勒级数近似展开，可推得与 AH方法形式相同的步进公式。确定其中 的公式是=(g(i)2/(h(i)-)，g(i)和 h(i)代表此方向的梯度和本征值，加和是对所有本征向量方向加和。通过迭代方法会解出N+1 个（N 代表势能面维数），将 按大小排列，则有(i)h(i)(i+1)。故选其中最小的 可使各个方向位移公式的(h(i)-)项皆为正，保证每步位移都向着极小点。选其中大于m个 Hessian 本征值的，将会在本征值最低的m个方向上沿其上的受力反方向位移提升能量，在其余N-m个方向上降低能量，由此确保优化到m阶鞍点，若m为 1 即用来找过渡态。所以用了这个方法寻找指定类型驻点不再像NR法对初猜位置Hessian 本征值符号有要求，而是直接通过选择 来设定向着何种鞍点位移。如果每步步长度超过了R，则乘以一个小于 1 的因子来减小步长。值得一提的是，与势能面维数N 的平方根近似成正比，随着体系尺度的增大，RFO的 对 NR法的二次近似就会趋现“校正过度”情况，产生大小不一致问题，可使用 SIRFO(Size independent RFO)方法解决，即 AH走步公式中的 改为/N0.5。2.1.2.2 P-RFO法(Partitioned-RFO)专用于优化过渡态，效率比RFO更高。RFO对所有方向的步进都使用同一个，而在 P-RFO中在指向过渡态的方向使用独立计算的(TS)，(TS)=g(TS)2/(h(TS)-(TS)，应选这个一元二次方程的最大的解，可保证在这个方向上升高能量。其余方向 的确定和RFO的公式一样，加和就不再包含指向过渡态的方向，并且选最小的 解以使这些方向能量降低。这里所谓指向过渡态的方向一般是指最低本征值的方向，在上述 RFO方法 m为 1 时也是如此假设（限于其形式RFO也只能用这最低模式），但有时会是其它的非最低的模式，P-RFO也可以将这样的模式作为指向过渡态的模式，见后文EF方法的讨论。2.1.2.3 QA法(Quadratic Approximation，二次逼近)名师资料总结-精品资料欢迎下载-名师精心整理-第 5 页，共 30 页 -确定 的公式是(X(i)2=(-g(i)/(h(i)-)2=R2，也就是说每一步移动的距离恰好是置信半径，这样步进速度较快。若优化到过渡态，计算 公式的加和中指向过渡态本征向量的那一项的 改为-，即 X(TS)=-g(TS)/(h(TS)+)。2.1.2.4 TRIM法(trust-region image minimization，置信区域镜像最小化)这个方法假设Hessian 本征值最小的方向的梯度和曲率符号与原本相反，而其它方向不变。经过这样的变化后原来的过渡态位置就成为了能量极小点(过渡态的 image)，这样就可以通过优化到极小点而得到过渡态。将TRIM的假设 g(TS)=-g(TS)，h(TS)=-h(TS)代入 AH方法的步进公式X(i,)=-g(i)/(h(i)-)，再使分子分母同乘以-1，可知在过渡态方向上的步进公式与其它方向区别仅在于反转了 的符号。又由于TRIM也是通过调整 使步进长度等于为置信半径，所以在公式的形式上与QA法找过渡态的公式完全一致，QA 与 TRIM可互为同义词。通过如上调整AH方法引入的 可使 NR法的步进更有效率、更稳定，还可以通过它改变步进公式在不同方向上的分母项符号，使优化过渡态的初猜点不限于过渡态的二次区域。可直接指定沿某个振动模式升高能量来找过渡态，即便当前点这个方向的Hessian 本征值可能是正值，例如从极小点开始跟踪至过渡态，见后文的EF方法。2.1.2.5 在高斯中的常见问题高斯中 opt=ts是使用 Berny 算法来找过渡态，需要提供一个初猜结构。Berny 默认的走步的方法是RFO/P-RFO(分别对于优化至极小值/鞍点），若加了Newton 选项，则走步基于NR法。每一步对 Hessian 矩阵的更新方法以UpdateMethod 选项指定，寻找极小点时默认用BFGS，找过渡态时默认用Bofill。Berny 算法还包括一些细节步骤在内，比如投影掉被冻结的变量、更新置信半径、设定了线搜索过程中几种方案等等，详见手册opt 关键字。使用了每步修正Hessian 的准牛顿法后，初猜的Hessian 矩阵质量明显影响结构收敛速度，它的不准确容易导致搜索过渡态失败（在高斯中默认使用价键力场得到Hessian）。这种情况需要昂贵的calcfc关键字以当前方法水平计算最初的Hessian 矩阵，若使用的方法在程序中支持解析二阶导数，速度会较好。或者用 readfc来读取包含了Hessian 矩阵信息的chk文件，可以先使用低水平方法进行简正振动分析得到chk 文件，再将之读入作为Hessian矩阵初猜，能够节约时间，但前提是此势能面对方法等级不敏感（一般如此）。使用了更准确的初猜后不仅可以增加找到过渡态的成功机率，还有助于在更短的优化步数内达到收敛标准。若使用calcall，则每一点都重新准确计算Hession，会更为可靠，但极为昂贵。高斯中 berny 方法寻找过渡态默认每步会检查Hessian 矩阵的本征值是否仅有一个为负，如果不符，就会提示“a wrong sign eigenvalue in hessian matrix”，经常一开始就报错，原因是初猜结构不符合这个条件，即便这个初猜通过berny 方法最终能够正确优化到过渡态，这时应加noeigentest选项避免本征值符号的检查，不符合要求也继续优化，但因此可能收敛到其它类型驻点。有时这种情况由初猜的Hessian 不准导致，可用calcfc解决。如果搜索的过渡态出现多个负本征值，可根据适当的虚频（高斯中以负数频率表示）振动方向调整结构以降低能量，直至剩下一个虚频，再重新优化。高斯中默认的置信半径为0.3 bohr，若优化中步长（RFO/P-RFO步）超过就会输出“Maximum step size(0.300)exceeded in Quadratic search”和“Step size scaled by xxx”，名师资料总结-精品资料欢迎下载-名师精心整理-第 6 页，共 30 页 -即乘以 xxx 调小步长至置信半径内。也可以使用iop(1/8=k)将置信半径改为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 与 DIIS 原理一致，但用于几何优化，这个方法趋于收敛到离初始位置最近的驻点，包括过渡态。下一步坐标X(new)=X-Hg，H 代表当前步的Hessian 逆矩阵，可见公式形式与 NR法是一致的，但是 X与 g 不再指当前步的坐标和梯度，而是由之前走过的点的坐标X(k)和 梯 度g(k)插 值 得 到 的，X=c(k)X(k)，g=c(k)g(k)，代 入 上 式 即X(new)=c(i)(X(i)-Hg(i)，其中是对之前全部走过的点加和。系数c(k)通过使误差向量 r 的模最小化得到，r=c(k)e(k)，并以c(k)=1为限制条件。e(k)常见有两种定义，一种是 e(k)=-Hg(k)，另一种更常用，是 e(k)=g(k)，可看出 GDIIS 利用的是已经搜索过的子空间中坐标与梯度的相关性，目的是估出梯度（即误差向量）的模尽可能小的坐标，这一点与后述的梯度模方法相似。此方法缺点是由于势能面复杂，步进中容易被拉到已经过的势能面的其它驻点而不能到达指定类型驻点，还容易走到类似肩膀形状的拐点，梯度虽小却不为0，由于不能达到收敛标准而反复在此处震荡。另外随优化步数增加，误差向量数目逐渐加大，会逐渐出现的误差向量之间的线性相关，导致伪收敛和数值不稳定问题。在改进的方法中将GDIIS 与更可靠的RFO方法结合，比较二者的步进方向和长度，并检验GDIIS 中的组合系数c，根据一定规则来决定每一步对之前走过的点的保留方式，必要时全部舍去而重新开始GDIIS。Gaussian 中用的这种改进的GDIIS 方法解决了上述问题同时提高了效率，速度等于或优于RFO方法，尤其是以低水平对势能面平坦的大体系优化时更为突出。GDIIS 计算量小，对 Hessian 矩阵很不敏感，可以在优化中不更新，也可以用 QN法更新来改善性能。此方法自 Gaussian98 起就是默认的半经验优化算法，其它方法下也可以用OPT的 gdiis关键词打开。2.1.4 梯度模优化(gradient norm minimization)势能面上的驻点，包括能量极小点、过渡态和高阶鞍点的势能梯度都为0，所以在相应于势能面的梯度模面上进行优化找到数值为0 的点，经过 Hessian 矩阵本征值符号的检验，就能得到过渡态。这相当于把搜索过渡态问题转化为了能量极小化问题，就有了更可靠的算法可用。（注：梯度模指的是势能梯度在各个维度分量平方和的平方根，即梯度大小的绝对值）。但是寻找数值更小点的优化方法比如最陡下降法只能找到离初始位置最近的极小点，若找到的梯度模面上的极小点数值大于0 则是势能面肩膀形拐点，没有什么用处，而这样的点收敛半径往往很大，例如图中在x=2 至 8 的区域内都会收敛到函数拐点，只有提供的初猜结构在x=1 和 9 附近很小的范围内才会收敛到过渡态，收敛半径太小，难以提供合理初猜。梯度模面上还多出一些极大点，如x=1.5 处，若使用收敛更快的NR法找极小点还容易收敛到这样没有意义的点上。基于这些原因，梯度模法很少使用。名师资料总结-精品资料欢迎下载-名师精心整理-第 7 页，共 30 页 -图 1 原函数与它的梯度模曲线。2.1.5 Dimer方法Dimer 方法是一种高效的定位过渡态的方法。这个方法定义了由两个点R1和 R2组成的一个Dimer，能量和所受势能力（由原始的势能面梯度造成受力，下同）分别为E1和 E2、F1 和F2。两个点间距为2R，R为定值。这两点的中间点为R，其受力 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、共轭梯度法。当 C0（极小点二次区域内），F 等于 F(R)平行于 Dimer 方向力的分量的负值，而没有垂直于Dimer 方向的力，促使 Dimer 尽快离开这个区域。由于 Dimer 的方向就是曲率最小的方向，在过渡态二次区域内就是指虚频方向，在Dimer 方法中F 的定义使这个方向以受力相反方向移动以升高能量，而其它方向顺着受力方向移动来最小化能量，可看出原理上与NR法相似。费时的计算Hessian 矩阵最小本征值以确定提升能量方向的过程被旋转 Dimer 这一步代替了，仅需要计算一阶导数。Dimer 法对初始位置要求很宽松，并不需要在过渡态二次区域内，若在极小点二次区域内就类似于后述的EF方法沿着最小振动模式爬坡。如果在高阶鞍点二次区域内，只在曲率最负的虚频方向沿着受力反向移动，在其它虚频方向上仍最小化能量，而不会像NR法收敛到高阶鞍点。名师资料总结-精品资料欢迎下载-名师精心整理-第 8 页，共 30 页 -图 2 右侧为 Dimer 法在 M ller-Brown 模型势上面搜索两个过渡态过程中Dimer 走的路径。势能面上往往有许多鞍点，Dimer 方法还可以做鞍点搜索。通过分子动力学方法给予Dimer一定动能，使之能够在势能面上广阔的区域内运动，根据一定标准提取轨迹中的一些点作为初猜，再执行标准Dimer 方法就可以得到许多不同的鞍点。Dimer 方法很适合双处理器并行，两个点的受力分别由两个处理器负责，速度可增加将近一倍。2.2 基于反应物与产物结构的算法2.2.1 同步转变方法(synchronous transit,ST)提供合理的初猜结构往往不易，ST方法可以只根据反应物和产物结构自动得到过渡态结构。“同步转变”这个名字强调的是反应路径上所有坐标一起变化，这是相对于后面提到的赝坐标法来说的（即只变化指定的坐标，尽管其它坐标优化后坐标也会变化）。ST 分为两种模型，最简单的就是LST 模型(Linear synchronous transit,线性同步转变)，这个方法假设反应过程中，反应物结构的每个坐标都是同步、线性地变化到产物结构。如果反应物、产物的坐标分别以向量A、B表示，则反应过程中的结构坐标可表示为(1-x)*A+x*B，x 由 0 逐渐变到 1 代表反应进度。注意LST并不是指反应中原子在真实空间上以直线运动，只有笛卡尔坐标下的LST才是如此，在内坐标下的LST，原子在真实空间中一般以弧线运动。以 LST的假设，反应路径在其所用坐标下的势能面图上可描述为一条直线，LST给出的过渡态就是这条直线上能量最高点（图3 的点 1）。LST的问题也很显著，其假设的坐标线性变化多数是错误的，绘制在势能面图上也多数不会是直线，故给出的过渡态也有较大偏差，容易带两个或多个虚频。比 LST更合理的是QST(quadratic synchronous transit,二次同步转变)，它假设反应路径在势能面上是一条二次曲线。QST在 LST得到的过渡态位置上，对 LST直线路径的垂直方向进行线搜索找到能量极小点A（图 3 的点 2）。QST给出的反应路径可以用经过反应物、A、产物的二次曲线来表示，如果这条路径上能量最大点的位置恰为A，则 A 就是 QST方法给出的过渡态；如果不是，则以最大点作为过渡态。若想结果更精确，可以再对这个最大点向垂直于路径的方向优化，再次得到A并检验，反复重复这个步骤，逐步找到能量更低、更准确名师资料总结-精品资料欢迎下载-名师精心整理-第 9 页，共 30 页 -的过渡态。QST方法在计算能力较低的年代曾是简单快速的获得过渡态和反应路径的方法，然而如今看来其结果是相当粗糙的，已极少单独使用，可以将其得到的过渡态作为AH法的初猜。图 3LST 与 QST方法示意图2.2.2 STQN 方法(Combined Synchronous Transit and Quasi-Newton Methods)STQN是 ST与 QN方法的结合（更准确地说是与EF法的结合）。但不要简单认为是按顺序独立执行这两步，即认为“先利用反应物和产物结构以ST方法得到粗糙过渡态，再以之作为初猜用 QN法精确寻找过渡态”是错误的。STQN 方法大意是：使结构从低能量的反应物出发，以 ST路径在当前位置切线为引导，沿着LST或 QST假设的反应路径行进（爬坡步），目的是使结构到达假设路径的能量最高处附近（真实过渡态二次区域附近）。当符合一定判据时就转换为 QN法寻找精确过渡态位置（EF步）。下面介绍具体步骤。先说明后面用到的切线的定义：STQN当中的 LST路径与前面ST部分介绍的LST 路径无异，都是直线，切线T 在优化中是不变的，就是反应物R 指向产物P 的单位向量。STQN 方法中的 QST路径定义与ST方法介绍的不同，走的不是二次曲线而是圆形的一段弧，如图 4 所示。这个圆弧经过R、P 以及优化中的当前步位置X，切线就是圆在X处的单位切线向量，圆弧和切线在每一步都是变化的。虽然QST路径比 LST更为合理，但对于QSTN方法，QST路径在收敛速度和成功机率上的优势并不显著。名师资料总结-精品资料欢迎下载-名师精心整理-第 10 页，共 30 页 -图 4STQN对 QST路径的定义STQN每一步执行内容如下：(1)首先重新计算或用QN法更新 Hessian。(2)按上述方法计算当前位置处的切线。(3)决定这一步是爬坡步还是EF步。如果是优化的第前两步，则一定认为是爬坡步，因为此时离过渡态区域还较远，应当先爬坡。如果是第3、4 步，则估算出在切线方向的位移，超过一定标准就是爬坡步，否则说明爬得差不多了就进入EF步找过渡态。如果是第 5 步之后，一般已离过渡态区域较近，故一定认为是EF步。(4)如果是爬坡步，则在切线方向上移动（将切线方向作为EF方法所跟踪的振动方向来计算位移大小）。如果是EF步，首先计算Hessian 各个本征向量的与切线重叠情况，如果有重叠大于0.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需额外输入猜测的过渡态，它直接作为X，一般比 QST2效果更好。对于经验不足的用户，用 STQN方法往往比只提供过渡态初猜的方法更为适合。注意产物和反应物应当使用同样方法同样基组进行优化，如果是多分子比如A+B=C+D 这样的反应，应当优化A和 B/C 和 D的复合物作为输入的产物/产物，而不是单独优化 A、B 然后拼到一起，因为形成范德华复合物后孤立的分子会有一定构象改变，能量也低于它们孤立状态的加和。2.2.3 赝坐标法(pseudo reaction coordinate)也称为坐标驱动法(Coordinate Driving)。这个方法在高斯中就是柔性扫描(Relaxed Scan)，即扫描一个变量，但每一步对其它变量自动进行优化，每一步得到的结构就是在这个变量为一定值情况下的最优结构。赝坐标法扫描的是反应物转变到产物过程中的关键坐标，比如扫描化学键断裂/生成反应中的键长。扫描的结果就是近似的IRC，可以再将能量最高点作为初猜找过渡态，或者用更小的步长再次扫描能量最高点附近找更精确的过渡态结构。这个找名师资料总结-精品资料欢迎下载-名师精心整理-第 11 页，共 30 页 -过渡态方法实际上用的是能量极小化优化过程，由于这样的算法比寻找过渡态的算法更为稳健，所以赝坐标法是颇可靠的，其它方法失败时可考虑这种方法。这个方法缺点是费时间，而且不适合通向过渡态路径中反应区域涉及多个坐标变化的反应过程，因为自定义扫描的内容很难全面、准确考虑到这些坐标变量的变化，结果难以说明问题，没有考虑进去的关键变量容易产生滞后效应(hysteresis effect)。比如乙烷由交叉构象变化到另一个交叉构象，需要经历重叠构象的过渡态，会涉及到三个HCCH二面角同时由60度变化到0 度，如果用赝坐标法只扫描其中一个HCCH 由 60 度变到 0 度，则每一步其它两个HCCH角一定会大于这个扫描的二面角，与实际不符。这是因为这两个角越小，分子的能量越高，每一步自动优化的时候它们更倾向于保持在大角度。最终到达过渡态时，所扫描的二面角到达了0 度，另外两个二面角却大于0 度，说明它们的运动比实际的过程滞后了。由于滞后效应，从反应物和产物两个方向扫描同相同的坐标，得到的路径也不同。上述简单的反应此方法滞后效应尚且严重，对于复杂变化，这种效应导致的问题更难以预测。故此方法确定的 IRC、过渡态不可靠，只建议对简单的反应使用这种方法，扫描变量的选择注意避免滞后效应。在高斯中此方法可以使用opt=modredundant 或 Opt=Z-matrix结合分子结构部分标记的扫描变量来实现。例如使用opt=modredundant并在分子结构末尾写上A 3 2 1 S 10 1.000000来指定 3 2 1原子组成的角度进行柔性扫描，共10 步，每步 1.0 度。如果不熟悉，也可以很方便地在GaussView 里的冗余坐标编辑器里面添加要柔性扫描的变量。如果只执行常规的某个变量的扫描，比如高斯中的scan 来找能量最高点作为初猜结构，对于简单体系可行，但对于复杂体系，这样忽略了此变量的变化导致分子其它部分结构的驰豫，如此得到的能量最高点作为过渡态初猜很不可靠，因为势能中掺入了不合理的结构造成的能量升高，使势能曲线形状改变。2.2.4 DHS 方法(Dewar-Healy-Stewart，亦称 Saddle 方法)与 LTP方法(Line-Then-Plane)DHS方法中第一步将反应和产物分别作为A点和 B点，确定哪个点能量低，比如A比 B低，就把 A点的结构向B点稍微做调整(5%)得到 A，然后限制变量空间中A 与 B的距离不变（即在超球面上）对A 进行优化得到A。将 A 与 B当作下一步的起始点A与 B，重复上述方法。这样反复进行迭代，若以序号n 代表第 n 次得到的 A 或 B，会依次得到例如A(1)、A(2)、B(1)、A(3