過渡態(tài)、反應路徑的計算方法及相關問題

過渡態(tài)、反應路徑的計算方法及相關問題SoberevaDepartmentofChemistry,UniversityofScienceandTechnologyBeijing,Beijing

100083,China前言：本文主要介紹過渡態(tài)、反應路徑的計算方法，并討論相關問題。由于這類算法極多，可以互相組合，限于精力不可能面面俱到展開，所以只介紹常用，或者實用價值有限但有啟發(fā)性的方法。文中圖片來自相關文獻，做了一定修改。由于本文作為帖子發(fā)布，文中無法插入復雜公式，故文中盡量將公式轉化為文字描述并加以解釋，這樣必然不如公式形式嚴謹，而且過于復雜的公式只能略過，但我想這樣做的好處是更易把握方法的梗概，有興趣可以進一步閱讀原文了解細節(jié)。對于Gaussian中可以實現(xiàn)的方法，文中對其在Gaussian中的使用進行了一些討論，希望能糾正一些網(wǎng)上流傳的誤區(qū)。雖然絕大多數(shù)人不專門研究計算方法，其中很多方法也不會用到，但多了解一下對開闊思路是很有好處的。文中指的“反應”包括構象變化、異構化、單分子反應等任何涉及到過渡態(tài)的變化過程?！胺磻铩迸c“產(chǎn)物”泛指這些過程的初態(tài)和末態(tài)?！皟?yōu)化”若未注明，包括優(yōu)化至極小點和優(yōu)化至過渡態(tài)。勢能面是高維的，但為了直觀以及表述方便，文中一般用二維勢能面模型來討論，應推廣至高維情況。限于純文本格式，向量、矩陣無法加粗表示，但容易自行判斷。目錄：過渡態(tài)過渡態(tài)搜索算法基于初猜結構的算法牛頓-拉弗森法(Newton-Raphson,NR)與準牛頓法(quasi-Newton,QN)AH方法(augmentedHessian)RFO法(RationalFunctionOptimization,有理函數(shù)優(yōu)化)P-RFO法(Partitioned-RFO)0人法(QuadraticApproximation,二次逼近)TRIM法(trust-regionimageminimization,置信區(qū)域鏡像最小化)在高斯中的常見問題GDIIS法(GeometryDirectInversionintheIterativeSubspace)梯度模優(yōu)化(gradientnormminimization)Dimer方法基于反應物與產(chǎn)物結構的算法同步轉變方法(synchronoustransit,ST)STQN方法(CombinedSynchronousTransitandQuasi-NewtonMethods)贋坐標法(pseudoreactioncoordinate)DHS方法(Dewar-Healy-Stewart,亦稱Saddle方法)與LTP方法(Line-Then-Plane)Ridge方法Step-and-Slide方法Muller-Brown方法CI-NEB、ANEBA方法基于反應物結構的算法最緩上升法(leaststeepascent,shallowestascent)本征向量/本征值跟蹤法(eigenvector/eigenvaluefollowing,EF。也稱modewalking/modefollowing/Walkingupvalleys)ARTn(activation-relaxationtechniquenouveau)梯度極值法(Gradientextremal,GE)約化梯度跟蹤(reducedgradientfollowing,RGF)等勢面搜索法(IsopotenialSearching)球形優(yōu)化(Sphereoptimization)全勢能面掃描過渡態(tài)相關問題無過渡態(tài)的反應途徑(barrierlessreactionpathways)Hammond-Leffler假設對稱性問題溶劑效應計算過渡態(tài)的建議流程內稟反應坐標(intrinsicreactioncoordinate,IRC)算法最陡下降法(Steepestdescent)IMK方法(Ishida-Morokuma-Kormornicki)Muller-Brown方法GS(Gonzalez-Schlegel)方法方法Dragmethod方法PEB方法(plainelasticband)Elber-Karplus方法SPW方法(Self-PenaltyWalk)LUP方法(LocallyUpdatedplanes)NEB方法(NudgedElasticBand)DNEB方法(DoubleNudgedElasticBand)String方法SimplifiedString方法尋找過渡態(tài)的chain-of-state方法CI-NEB方法ANEBA方法(adaptivenudgedelasticbandapproach)chain-of-states方法的一些特點高斯中opt關鍵字的path=M方法CPK方法(ConjugatePeakRefinement)1.過渡態(tài)過渡態(tài)結構指的是勢能面上反應路徑上的能量最高點，它通過最小能量路徑（minimumenergypath,MEP）連接著反應物和產(chǎn)物的結構（如果是多步反應的機理，則這里所指反應物或產(chǎn)物包括中間體）。對于多分子之間的反應，更確切來講過渡態(tài)結構連接的是它們由無窮遠接近后因為范德華力和靜電力形成的復合物結構，以及反應完畢但尚未無限遠離時的復合物結構。確定過渡態(tài)有助于了解反應機理，以及通過勢壘高度計算反應速率。一般來講，勢壘小于21kcal/mol就可以在室溫下發(fā)生。在勢能面上，過渡態(tài)結構的能量對坐標的一階導數(shù)為0，只有在反應坐標方向上曲率（對坐標二階導數(shù)）為負，而其它方向上皆為正，是能量面上的一階鞍點。過渡態(tài)結構的能量二階導數(shù)矩陣（Hessian矩陣）的本征值僅有一個負值，這個負值也就是過渡態(tài)擁有唯一虛頻的來源。若將分子振動簡化成諧振子模型，這個負值便是頻率公式中的力常數(shù)，開根號后即得虛數(shù)。分子構象轉變、化學反應過程中往往都有過渡態(tài)的存在，即這個過程在勢能面上的運動往往都會經(jīng)歷滿足上述條件的一點?；瘜W反應的過渡態(tài)更確切應當成為“反應過渡態(tài)”。需要注意的是化學反應未必都經(jīng)歷過渡態(tài)結構。由于過渡態(tài)結構存在時間極短，所以很難通過實驗方法獲得，直到飛秒脈沖激光光譜的出現(xiàn)才使檢驗反應機理為可能。計算化學方法在目前是預測過渡態(tài)的最有力武器，盡管計算上仍有一些困難，比如其附近勢能面相對于平衡結構更為平坦得多、低水平方法難以準確描述、難以預測過渡態(tài)結構、缺乏絕對可靠的方法（如優(yōu)化到能量極小點可用的最陡下降法）等。搜索過渡態(tài)的算法一般結合從頭算、DFT方法，在半經(jīng)驗、或者小基組條件下，難以像描述平衡結構一樣正確描述過渡態(tài)結構，使得計算尺度受到了限制。結合分子力場可以描述構象變化的過渡態(tài)，但不適用描述反應過渡態(tài)，因為大部分分子力場的勢函數(shù)不允許分子拓撲結構的改變，雖然也有一些力場如ReaxFF可以支持，有的力場還有對應的過渡態(tài)原子類型，但目前來看適用面仍然較窄，而且不夠精確，盡管更為快速。注：嚴格來說，“過渡結構”是指勢能面上反應路徑上的能量最高點，而“過渡態(tài)”是指自由能面上反應路徑上的能量最高點，由于自由能變主要貢獻自勢能部分，所以多數(shù)情況二者結構近似一致。但隨著溫度升高，往往熵變的貢獻導致自由能面與勢能面形狀發(fā)生明顯偏離，從而導致過渡結構與過渡態(tài)明顯偏離，兩個詞就不能混用了。但本文不涉及相關問題，故文中過渡態(tài)、過渡結構一律指勢能面上反應路徑上的能量最高點。2.過渡態(tài)搜索算法基于初猜結構的算法牛頓-拉弗森法(Newton-Raphson,NR)與準牛頓法(quasi-Newton,QN)NR法是尋找函數(shù)一階導數(shù)為0(駐點)位置的方法。通過對能量函數(shù)的泰勒級數(shù)的二階近似展開，然后使用穩(wěn)態(tài)條件dE/dr=0,可導出步進公式：下一步的坐標向量=當前坐標向量-能量一階導數(shù)向量*Hessian矩陣的逆矩陣。在勢能面上以NR法最終找到的結果是與初猜位置Hessian矩陣本征值正負號一致、離初猜結構最近的駐點，由于能量極小點、過渡態(tài)和高階鞍點的能量一階導數(shù)皆為0,故都可以用NR法尋找。對于純二次形函數(shù)NR法僅需一步即可找到正確位置，而勢能面遠比之復雜，所以需要反復走步直至收斂。也因為勢能面這個特點，為了改進優(yōu)化，實際應用中NR法一般還結合線搜索步(linesearch),對于優(yōu)化至極小點，就是找當前點與NR法算出來的下一點的連線上的能量極小點作為實際下一步結構；若優(yōu)化至過渡態(tài)，且連線方向主要指向過渡態(tài)，則找的是連線上能量極大點，若主要指向其它方向則找連線的能量極小點，若指向二者程度均等則一般不做線搜索。由于精確的線搜索很花時間，所以一般只是在連線的當前位置附近計算幾個點的能量，以高階多項式擬和后取其最小/最大點。NR法每一步需要計算Hessian矩陣并且求其逆，所以十分昂貴。QN法與NR法的走步原理一樣，但Hessian矩陣最初是用低級或經(jīng)驗方法猜出來的，每一步優(yōu)化中通過當前及前一步的梯度和坐標對Hessian矩陣逆矩陣逐漸修正。由于只需計算一階導數(shù)，即便Hessian矩陣不準確造成所需收斂步數(shù)增加，但一般仍比NR法速度快得多。QN法泛指基于此原理的一類方法，常用的是BFGS(BroydenFletcherGoldfarbShanno)，此法對Hessian的修正保持其對稱性和正定性，最適合幾何優(yōu)化，但顯然不能用于找過渡態(tài)。還有DFP(Davidon-Fletcher-Powell),MS(Murtagh-Sargent,亦稱symmetrierank1，SR1)，PSB(Powell-symmetric-Broyden)。也有混合方法，如Bofill法是PSB和MS法對Hessian修正量的權重線性組合，比二者獨立使用更優(yōu)，權重系數(shù)通過位移、梯度改變量和當前Hessian計算得到，它對Hessian的修正不強制正定，很適宜搜索過渡態(tài)。將NR步進公式放到Hessian本征向量空間下其意義更為明顯(此時Hessian為對角矩陣)，可看出在每個方向上的位移就是這個方向勢能的負梯度除以對應的本征值,比如在i方向上的位移可寫為AX(i)=-g(i)/h(i),在受力越大、越平坦的方向位移越大。每一步實際位移就是這些方向上位移的矢量和。對于尋找過渡態(tài)，因為虛頻方向對應Hessian本征值為負，使位移為受力相反方向，所以NR法在過渡態(tài)附近每一步都是使虛頻方向能量升高，而在其它正交的方向朝著能量降低的方向位移，通過這個原理步進到過渡態(tài)。若有n個虛頻，則NR法就在n個方向升高能量而其它方向降低能量找到n階鞍點。由于NR法的這個特點，為找到正確類型的駐點，初猜結構必須在目標結構的二次區(qū)域(quadraticregion)內。所謂的二次區(qū)域，是指駐點附近保持Hessian矩陣本征值符號不變的區(qū)域，它的形狀可以用多變量的二次函數(shù)近似描述，例如二維勢能面情況下這樣的區(qū)域可以用F(x,y)二A*x"2+B*y"2+C*x+D*y+E*x*y來近似描述。對于能量極小點，就是指初猜點在目標結構附近Hessian矩陣為正定矩陣的范圍；對于找過渡態(tài)，就需要初猜點在它附近含有且僅含有一個負本征值的范圍內。并且這個范圍內不能有其它同類駐點比目標結構距離初猜結構更近。NR法方便之處是只需要提供一個初猜結構即可，但是由于過渡態(tài)二次區(qū)域很小(相對于能量極小點來講)，復雜反應過渡態(tài)又不容易估計，故對使用者的直覺和經(jīng)驗有一定要求，即便是老手，也往往需要反復嘗試。NR法對初猜結構比較敏感，離過渡態(tài)越近所需收斂步數(shù)越少，成功機率越高。模版法可以幫助給出合理的初猜，也就是如果已經(jīng)知道其它機理相同的反應的過渡態(tài)結構，可以保持反應位點部分的結構不變而替換周圍的原子，使之變成自己要研究的化合物反應的初猜結構。AH方法(AugmentedHessian)AH方法并不是獨立的尋找過渡態(tài)的方法，而是通過修改原始Hessian矩陣來調整NR法步進的長度和方向的一種方法。在NR法的步進公式中加入了一個移位參數(shù)入，式子變?yōu)锳X(i,入)二-g(i)/(h(i)-入)，NR法相當于入等于0時的特例。入控制著每步步進距離，它與h(i)的相對大小也控制著這個坐標上的步進方向。根據(jù)設定入方法的不同，常見的有RFO、P-RFO和QA/TRIM。這些方法每一步也使用QN方法來快速地更新Hessian。下面提及的置信半徑R(Trustradius)是指二階泰勒級數(shù)展開這種近似的合理的區(qū)域，可以在優(yōu)化過程中固定也可以動態(tài)改變，比如下一步位置的實際能量與使用二階泰勒級數(shù)展開預測的能量符合較好則加大R,反之減小。優(yōu)化的每一步移動距離不應超過R,否則可能進入二階泰勒級數(shù)展開近似的失效區(qū)域，NR法在勢能面平坦的時候容易超過這個范圍，應調整入避免。RFO法(RationalFunctionOptimization,有理函數(shù)優(yōu)化)對能量函數(shù)根據(jù)有理近似展開，而不是NR法的二階泰勒級數(shù)近似展開，可推得與AH方法形式相同的步進公式。確定其中入的公式是入=E(g(i)"2/(h(i)-入))，g(i)和h(i)代表此方向的梯度和本征值，加和是對所有本征向量方向加和。通過迭代方法會解出N+1個入(N代表勢能面維數(shù))，將入按大小排列，則有入(i)Wh(i)W入(i+1)。故選其中最小的入可使各個方向位移公式的(h(i)-入)項皆為正，保證每步位移都向著極小點。選其中大于m個Hessian本征值的入，將會在本征值最低的m個方向上沿其上的受力反方向位移提升能量，在其余N-m個方向上降低能量，由此確保優(yōu)化到m階鞍點，若m為1即用來找過渡態(tài)。所以用了這個方法尋找指定類型駐點不再像NR法對初猜位置Hessian本征值符號有要求，而是直接通過選擇入來設定向著何種鞍點位移。如果每步步長度超過了R,則乘以一個小于1的因子來減小步長。值得一提的是，入與勢能面維數(shù)N的平方根近似成正比，隨著體系尺度的增大，RFO的入對NR法的二次近似就會趨現(xiàn)“校正過度”情況，產(chǎn)生大小不一致問題，可使用SIRFO(SizeindependentRFO)方法解決，即AH走步公式中的入改為入/N"。P-RFO法(Partitioned-RFO)專用于優(yōu)化過渡態(tài)，效率比RFO更高。RFO對所有方向的步進都使用同一個入，而在P-RFO中在指向過渡態(tài)的方向使用獨立計算的入(TS)，入(TS)=g(TS廠2/(h(TS)-入(TS)),應選這個一元二次方程的最大的解，可保證在這個方向上升高能量。其余方向入的確定和RFO的公式一樣，加和就不再包含指向過渡態(tài)的方向，并且選最小的入解以使這些方向能量降低。這里所謂指向過渡態(tài)的方向一般是指最低本征值的方向，在上述RFO方法m為1時也是如此假設邙限于其形式RFO也只能用這最低模式)，但有時會是其它的非最低的模式,P-RFO也可以將這樣的模式作為指向過渡態(tài)的模式，見后文EF方法的討論。0人法（QuadraticApproximation,二次逼近）確定入的公式是(AX(i))2=E(-g(i)/(h(i)-入))2=R2,也就是說每一步移動的距離恰好是置信半徑，這樣步進速度較快。若優(yōu)化到過渡態(tài)，計算入公式的加和中指向過渡態(tài)本征向量的那一項的入改為-入，即AX(TS)=-g(TS)/(h(TS)+入)。TRIM法(trust-regionimageminimization,置信區(qū)域鏡像最小化)這個方法假設Hessian本征值最小的方向的梯度和曲率符號與原本相反，而其它方向不變。經(jīng)過這樣的變化后原來的過渡態(tài)位置就成為了能量極小點(過渡態(tài)的image)，這樣就可以通過優(yōu)化到極小點而得到過渡態(tài)。將TRIM的假設g(TS)'=-g(TS)，h(TS)'=-h(TS)代入AH方法的步進公式AX(i,入)=-g(i)/(h(i)-入)，再使分子分母同乘以-1，可知在過渡態(tài)方向上的步進公式與其它方向區(qū)別僅在于反轉了入的符號。又由于TRIM也是通過調整入使步進長度等于為置信半徑，所以在公式的形式上與QA法找過渡態(tài)的公式完全一致，QA與TRIM可互為同義詞。通過如上調整AH方法引入的入可使NR法的步進更有效率、更穩(wěn)定，還可以通過它改變步進公式在不同方向上的分母項符號，使優(yōu)化過渡態(tài)的初猜點不限于過渡態(tài)的二次區(qū)域。可直接指定沿某個振動模式升高能量來找過渡態(tài)，即便當前點這個方向的Hessian本征值可能是正值，例如從極小點開始跟蹤至過渡態(tài)，見后文的EF方法。在高斯中的常見問題高斯中opt=ts是使用Berny算法來找過渡態(tài)，需要提供一個初猜結構。Berny默認的走步的方法是RFO/P-RF0(分別對于優(yōu)化至極小值/鞍點)，若加了Newton選項，則走步基于NR法。每一步對Hessian矩陣的更新方法以UpdateMethod選項指定，尋找極小點時默認用BFGS,找過渡態(tài)時默認用Bofill。Berny算法還包括一些細節(jié)步驟在內，比如投影掉被凍結的變量、更新置信半徑、設定了線搜索過程中幾種方案等等，詳見手冊。pt關鍵字。使用了每步修正Hessian的準牛頓法后，初猜的Hessian矩陣質量明顯影響結構收斂速度，它的不準確容易導致搜索過渡態(tài)失敗(在高斯中默認使用價鍵力場得到Hessian)。這種情況需要昂貴的calcfc關鍵字以當前方法水平計算最初的Hessian矩陣，若使用的方法在程序中支持解析二階導數(shù)，速度會較好?；蛘哂胷eadfc來讀取包含了Hessian矩陣信息的chk文件，可以先使用低水平方法進行簡正振動分析得到chk文件，再將之讀入作為Hessian矩陣初猜，能夠節(jié)約時間，但前提是此勢能面對方法等級不敏感(一般如此)。使用了更準確的初猜后不僅可以增加找到過渡態(tài)的成功機率，還有助于在更短的優(yōu)化步數(shù)內達到收斂標準。若使用calcall，則每一點都重新準確計算Hession，會更為可靠，但極為昂貴。高斯中berny方法尋找過渡態(tài)默認每步會檢查Hessian矩陣的本征值是否僅有一個為負，如果不符，就會提示“awrongsigneigenvalueinhessianmatrix”，經(jīng)常一開始就報錯,原因是初猜結構不符合這個條件，即便這個初猜通過berny方法最終能夠正確優(yōu)化到過渡態(tài)，這時應加noeigentest選項避免本征值符號的檢查，不符合要求也繼續(xù)優(yōu)化，但因此可能收斂到其它類型駐點。有時這種情況由初猜的Hessian不準導致，可用calcfc解決。如果搜索的過渡態(tài)出現(xiàn)多個負本征值，可根據(jù)適當?shù)奶擃l(高斯中以負數(shù)頻率表示)振動方向調整結構以降低能量，直至剩下一個虛頻，再重新優(yōu)化。高斯中默認的置信半徑為bohr，若優(yōu)化中步長(RF0/P-RF0步)超過就會輸出“Maximumstepsize( exceededinQuadraticsearch”和“Stepsizescaledbyxxx”，即乘以xxx調小步長至置信半徑內。也可以使用iop(l/8二k)將置信半徑改為k*bohr(lbohr=埃)，調大后往往可以顯著減少收斂步數(shù)，很適合勢能面平坦的大體系。注意并不是每一步的步長都固定為k*bohr，若沒超過置信半徑則步長并不因此改變。尋找極小點時默認為允許動態(tài)改變置信半徑，此時iop(1/8)設的就是最初的置信半徑，對于尋找過渡態(tài)默認為關閉此功能(相當于用了NoTrustUpdate),可以使用trustupdate關鍵字來打開這個功能。GDIIS法(GeometryDirectInversionintheIterativeSubspace)GDIIS與DIIS原理一致，但用于幾何優(yōu)化，這個方法趨于收斂到離初始位置最近的駐點，包括過渡態(tài)。下一步坐標X(new)二X"-H'g",H'代表當前步的Hessian逆矩陣，可見公式形式與NR法是一致的，但是X"與g"不再指當前步的坐標和梯度，而是由之前走過的點的坐標X(k)和梯度g(k)插值得到的，X"=Ec(k)X(k)，g"=Ec(k)g(k)，代入上式即X(new)=Ec(i)(X(i)-H'g(i))，其中工是對之前全部走過的點加和。系數(shù)c(k)通過使誤差向量r的模最小化得到，r=工c(k)e(k)，并以Ec(k)=1為限制條件。e(k)常見有兩種定義，一種是e(k)=-H'g(k)，另一種更常用，是e(k)=g(k)，可看出GDIIS利用的是已經(jīng)搜索過的子空間中坐標與梯度的相關性，目的是估出梯度(即誤差向量)的模盡可能小的坐標，這一點與后述的梯度模方法相似。此方法缺點是由于勢能面復雜，步進中容易被拉到已經(jīng)過的勢能面的其它駐點而不能到達指定類型駐點，還容易走到類似肩膀形狀的拐點，梯度雖小卻不為0，由于不能達到收斂標準而反復在此處震蕩。另外隨優(yōu)化步數(shù)增加，誤差向量數(shù)目逐漸加大，會逐漸出現(xiàn)的誤差向量之間的線性相關，導致偽收斂和數(shù)值不穩(wěn)定問題。在改進的方法中將GDIIS與更可靠的RFO方法結合，比較二者的步進方向和長度，并檢驗GDIIS中的組合系數(shù)c,根據(jù)一定規(guī)則來決定每一步對之前走過的點的保留方式，必要時全部舍去而重新開始GDIIS。Gaussian中用的這種改進的GDIIS方法解決了上述問題同時提高了效率，速度等于或優(yōu)于RFO方法，尤其是以低水平對勢能面平坦的大體系優(yōu)化時更為突出。GDIIS計算量小，對Hessian矩陣很不敏感，可以在優(yōu)化中不更新，也可以用QN法更新來改善性能。此方法自Gaussian98起就是默認的半經(jīng)驗優(yōu)化算法，其它方法下也可以用OPT的gdiis關鍵詞打開。梯度模優(yōu)化(gradientnormminimization)勢能面上的駐點，包括能量極小點、過渡態(tài)和高階鞍點的勢能梯度都為0，所以在相應于勢能面的梯度模面上進行優(yōu)化找到數(shù)值為0的點，經(jīng)過Hessian矩陣本征值符號的檢驗，就能得到過渡態(tài)。這相當于把搜索過渡態(tài)問題轉化為了能量極小化問題，就有了更可靠的算法可用。(注：梯度模指的是勢能梯度在各個維度分量平方和的平方根，即梯度大小的絕對值)。但是尋找數(shù)值更小點的優(yōu)化方法比如最陡下降法只能找到離初始位置最近的極小點，若找到的梯度模面上的極小點數(shù)值大于0則是勢能面肩膀形拐點，沒有什么用處，而這樣的點收斂半徑往往很大，例如圖中在x=2至8的區(qū)域內都會收斂到函數(shù)拐點，只有提供的初猜結構在x=1和9附近很小的范圍內才會收斂到過渡態(tài)，收斂半徑太小，難以提供合理初猜。梯度模面上還多出一些極大點，如x=處，若使用收斂更快的NR法找極小點還容易收斂到這樣沒有意義的點上?；谶@些原因，梯度模法很少使用。［圖1］原函數(shù)與它的梯度模曲線。Dimer方法Dimer方法是一種高效的定位過渡態(tài)的方法。這個方法定義了由兩個點R1和R2組成的一個Dimer,能量和所受勢能力（由原始的勢能面梯度造成受力，下同）分別為E1和E2、F1和F2。兩個點間距為2AR,AR為定值。這兩點的中間點為R,其受力F（R）=（Fl+F2）/2。Dimer的總能量為E=（El+E2）/2。這個方法的每一步包括平移Dimer和旋轉Dimer兩步。旋轉Dimer:保持Rl、R2中點位置R不變作為軸，旋轉Dimer直到總能量E最小。通過推導可知在旋轉過程中，E與R點在dimer方向（R1-R2方向）上的曲率關系C是線性的，即最小化E的過程就是最小化C的過程。所以每一步的Dimer方向都是曲率最小方向，當最終R收斂到過渡態(tài)位置時，Dimer就會平行于虛頻方向。平移Dimer：Dimer根據(jù)受力F'移動R的位置，結合不同方法有具體步進方式，如quick-win、共軛梯度法。當C<0（過渡態(tài)或高階鞍點的二次區(qū)域內），F(xiàn)'等于將F（R）平行于Dimer方向力的分量符號反轉；當C>0（極小點二次區(qū)域內），F(xiàn)'等于F（R）平行于Dimer方向力的分量的負值，而沒有垂直于Dimer方向的力，促使Dimer盡快離開這個區(qū)域。由于Dimer的方向就是曲率最小的方向，在過渡態(tài)二次區(qū)域內就是指虛頻方向，在Dimer方法中F'的定義使這個方向以受力相反方向移動以升高能量，而其它方向順著受力方向移動來最小化能量，可看出原理上與NR法相似。費時的計算Hessian矩陣最小本征值以確定提升能量方向的過程被旋轉Dimer這一步代替了，僅需要計算一階導數(shù)。Dimer法對初始位置要求很寬松，并不需要在過渡態(tài)二次區(qū)域內，若在極小點二次區(qū)域內就類似于后述的EF方法沿著最小振動模式爬坡。如果在高階鞍點二次區(qū)域內，只在曲率最負的虛頻方向沿著受力反向移動，在其它虛頻方向上仍最小化能量，而不會像NR法收斂到高階鞍點。［圖2］右側為Dimer法在MUller-Brown模型勢上面搜索兩個過渡態(tài)過程中Dimer走的路徑。勢能面上往往有許多鞍點，Dimer方法還可以做鞍點搜索。通過分子動力學方法給予Dimer一定動能，使之能夠在勢能面上廣闊的區(qū)域內運動，根據(jù)一定標準提取軌跡中的一些點作為初猜，再執(zhí)行標準Dimer方法就可以得到許多不同的鞍點。Dimer方法很適合雙處理器并行,兩個點的受力分別由兩個處理器負責，速度可增加將近一倍。基于反應物與產(chǎn)物結構的算法同步轉變方法（synchronoustransit,ST）提供合理的初猜結構往往不易，ST方法可以只根據(jù)反應物和產(chǎn)物結構自動得到過渡態(tài)結構。“同步轉變”這個名字強調的是反應路徑上所有坐標一起變化，這是相對于后面提到的贗坐標法來說的（即只變化指定的坐標，盡管其它坐標優(yōu)化后坐標也會變化）。ST分為兩種模型，最簡單的就是LST模型（Linearsynchronoustransit,線性同步轉變），這個方法假設反應過程中，反應物結構的每個坐標都是同步、線性地變化到產(chǎn)物結構。如果反應物、產(chǎn)物的坐標分別以向量A、B表示，則反應過程中的結構坐標可表示為（1-x）*A+x*B,x由0逐漸變到1代表反應進度。注意LST并不是指反應中原子在真實空間上以直線運動，只有笛卡爾坐標下的LST才是如此，在內坐標下的LST,原子在真實空間中一般以弧線運動。以LST的假設，反應路徑在其所用坐標下的勢能面圖上可描述為一條直線，LST給出的過渡態(tài)就是這條直線上能量最高點（圖3的點1）。LST的問題也很顯著，其假設的坐標線性變化多數(shù)是錯誤的，繪制在勢能面圖上也多數(shù)不會是直線，故給出的過渡態(tài)也有較大偏差，容易帶兩個或多個虛頻。比LST更合理的是QST（quadraticsynchronoustransit,二次同步轉變），它假設反應路徑在勢能面上是一條二次曲線。QST在LST得到的過渡態(tài)位置上，對LST直線路徑的垂直方向進行線搜索找到能量極小點A（圖3的點2）。QST給出的反應路徑可以用經(jīng)過反應物、A、產(chǎn)物的二次曲線來表示，如果這條路徑上能量最大點的位置恰為A，則A就是QST方法給出的過渡態(tài)；如果不是，則以最大點作為過渡態(tài)。若想結果更精確，可以再對這個最大點向垂直于路徑的方向優(yōu)化，再次得到A并檢驗，反復重復這個步驟，逐步找到能量更低、更準確的過渡態(tài)。QST方法在計算能力較低的年代曾是簡單快速的獲得過渡態(tài)和反應路徑的方法，然而如今看來其結果是相當粗糙的，已極少單獨使用，可以將其得到的過渡態(tài)作為AH法的初猜。[圖3]LST與QST方法示意圖STQN方法（CombinedSynchronousTransitandQuasi-NewtonMethods）STQN是ST與QN方法的結合（更準確地說是與EF法的結合）。但不要簡單認為是按順序獨立執(zhí)行這兩步，即認為“先利用反應物和產(chǎn)物結構以ST方法得到粗糙過渡態(tài)，再以之作為初猜用QN法精確尋找過渡態(tài)”是錯誤的oSTQN方法大意是：使結構從低能量的反應物出發(fā)，以ST路徑在當前位置切線為引導，沿著LST或QST假設的反應路徑行進（爬坡步），目的是使結構到達假設路徑的能量最高處附近（真實過渡態(tài)二次區(qū)域附近）。當符合一定判據(jù)時就轉換為QN法尋找精確過渡態(tài)位置（EF步）。下面介紹具體步驟。先說明后面用到的切線的定義：STQN當中的LST路徑與前面ST部分介紹的LST路徑無異，都是直線，切線T在優(yōu)化中是不變的，就是反應物R指向產(chǎn)物P的單位向量。STQN方法中的QST路徑定義與ST方法介紹的不同，走的不是二次曲線而是圓形的一段弧，如圖4所示。這個圓弧經(jīng)過R、P以及優(yōu)化中的當前步位置X，切線就是圓在X處的單位切線向量，圓弧和切線在每一步都是變化的。雖然QST路徑比LST更為合理，但對于QSTN方法，QST路徑在收斂速度和成功機率上的優(yōu)勢并不顯著。[圖4]STQN對QST路徑的定義STQN每一步執(zhí)行內容如下：（1）首先重新計算或用QN法更新Hessian。（2）按上述方法計算當前位置處的切線。（3）決定這一步是爬坡步還是EF步。如果是優(yōu)化的第前兩步，則一定認為是爬坡步，因為此時離過渡態(tài)區(qū)域還較遠，應當先爬坡。如果是第3、4步，則估算出在切線方向的位移，超過一定標準就是爬坡步，否則說明爬得差不多了就進入EF步找過渡態(tài)。如果是第5步之后，一般已離過渡態(tài)區(qū)域較近，故一定認為是EF步。（4）如果是爬坡步，則在切線方向上移動（將切線方向作為EF方法所跟蹤的振動方向來計算位移大?。Ｈ绻荅F步，首先計算Hessian各個本征向量的與切線重疊情況，如果有重疊大于的本征向量，則以EF法跟蹤本征值最大的本征向量來移動，相當于繼續(xù)向上爬。如果沒有大于的，就跟蹤最小本征值的本征向量移動來尋找過渡態(tài)。(5)步長長度若大于標準則調小，默認bohr。(6)根據(jù)預置受力、位移標準判斷是否已收斂，收斂則結束循環(huán)。注意，ST方法中具體包含LST和QST兩種方法，STQN也用到了LST和QST兩種反應路徑的假設。高斯中的LST方法指的是ST中的LST方法，而QST2/3指的是利用QST路徑假設的STQN方法，它們原理上截然不同，不要混淆。高斯中的QST2只需輸入反應物和產(chǎn)物結構，通過幾何方法估出STQN的初始步結構X。QST3需額外輸入猜測的過渡態(tài)，它直接作為X,一般比QST2效果更好。對于經(jīng)驗不足的用戶，用STQN方法往往比只提供過渡態(tài)初猜的方法更為適合。注意產(chǎn)物和反應物應當使用同樣方法同樣基組進行優(yōu)化，如果是多分子比如A+B=C+D這樣的反應，應當優(yōu)化A和B/C和D的復合物作為輸入的產(chǎn)物/產(chǎn)物，而不是單獨優(yōu)化A、B然后拼到一起，因為形成范德華復合物后孤立的分子會有一定構象改變，能量也低于它們孤立狀態(tài)的加和。贋坐標法(pseudoreactioncoordinate)也稱為坐標驅動法(CoordinateDriving)。這個方法在高斯中就是柔性掃描(RelaxedScan)，即掃描一個變量，但每一步對其它變量自動進行優(yōu)化，每一步得到的結構就是在這個變量為一定值情況下的最優(yōu)結構。贋坐標法掃描的是反應物轉變到產(chǎn)物過程中的關鍵坐標，比如掃描化學鍵斷裂/生成反應中的鍵長。掃描的結果就是近似的IRC,可以再將能量最高點作為初猜找過渡態(tài)，或者用更小的步長再次掃描能量最高點附近找更精確的過渡態(tài)結構。這個找過渡態(tài)方法實際上用的是能量極小化優(yōu)化過程，由于這樣的算法比尋找過渡態(tài)的算法更為穩(wěn)健，所以贋坐標法是頗可靠的，其它方法失敗時可考慮這種方法。這個方法缺點是費時間，而且不適合通向過渡態(tài)路徑中反應區(qū)域涉及多個坐標變化的反應過程，因為自定義掃描的內容很難全面、準確考慮到這些坐標變量的變化，結果難以說明問題，沒有考慮進去的關鍵變量容易產(chǎn)生滯后效應(hysteresiseffect)。比如乙烷由交叉構象變化到另一個交叉構象，需要經(jīng)歷重疊構象的過渡態(tài)，會涉及到三個HCCH二面角同時由60度變化到0度，如果用贋坐標法只掃描其中一個HCCH由60度變到0度，則每一步其它兩個HCCH角一定會大于這個掃描的二面角，與實際不符。這是因為這兩個角越小，分子的能量越高，每一步自動優(yōu)化的時候它們更傾向于保持在大角度。最終到達過渡態(tài)時，所掃描的二面角到達了0度，另外兩個二面角卻大于0度，說明它們的運動比實際的過程滯后了。由于滯后效應，從反應物和產(chǎn)物兩個方向掃描同相同的坐標，得到的路徑也不同。上述簡單的反應此方法滯后效應尚且嚴重，對于復雜變化，這種效應導致的問題更難以預測。故此方法確定的IRC、過渡態(tài)不可靠，只建議對簡單的反應使用這種方法，掃描變量的選擇注意避免滯后效應。在高斯中此方法可以使用opt=modredundant或Opt=Z-matrix結合分子結構部分標記的掃描變量來實現(xiàn)。例如使用opt=modredundant并在分子結構末尾寫上A321S10來指定321原子組成的角度進行柔性掃描，共10步,每步度。如果不熟悉，也可以很方便地在GaussView里的冗余坐標編輯器里面添加要柔性掃描的變量。如果只執(zhí)行常規(guī)的某個變量的掃描，比如高斯中的scan來找能量最高點作為初猜結構，對于簡單體系可行，但對于復雜體系，這樣忽略了此變量的變化導致分子其它部分結構的馳豫，如此得到的能量最高點作為過渡態(tài)初猜很不可靠，因為勢能中摻入了不合理的結構造成的能量升高，使勢能曲線形狀改變。DHS方法(Dewar-Healy-Stewart,亦稱Saddle方法)與LTP方法(Line-Then-Plane)DHS方法中第一步將反應和產(chǎn)物分別作為A點和B點，確定哪個點能量低，比如A比B低，就把A點的結構向B點稍微做調整(~5%)得到A'，然后限制變量空間中A'與B的距離不變(即在超球面上)對A'進行優(yōu)化得到A''。將A''與B當作下一步的起始點A與B,重復上述方法。這樣反復進行迭代，若以序號n代表第n次得到的A''或B'',會依次得到例如A''(1)、A''(2)、B''(l)、A''(3) 直到A與B十分接近時才停止迭代，此位置就是過渡態(tài)。將得到的全部A''(n)按序號n依次連接，B''(n)也按序號依次連接，再將序號最大的A''(n)與B''(n)連接，得到的就是近似的IRC。LTP與DHS方法基本一致，不同的是每步是在垂直于A'與B連線的超平面上優(yōu)化。DHS方法雖然可以很快地走到過渡態(tài)附近的位置，但是越往后每步的AB距離縮近也越少，故并不能有效率地貼近過渡態(tài)。然而每步的在連線上調整的距離不可過大，否則可能造成一側的點跨過過渡態(tài)勢壘跑到另一側得到錯誤結果。[圖5]DHS方法示意圖Ridge方法第一步時將反應物、產(chǎn)物作為A點和B點，在其LST的路徑上找到能量最大點C，然后在AC與BC直線上相距C為s的位置上分別設一點A'和B'，將A'與B'分別沿著此處勢能面負梯度優(yōu)化p距離，將得到的A''與B''作為下一步的A和B。反復進行這個步驟，收斂后C的位置就是過渡態(tài)位置。s和p是計算過程中動態(tài)調節(jié)的參數(shù)，對結果影響較大，它們應當隨C逐漸接近過渡態(tài)而減小，可設若當前步的C能量高于上一步的C,則減小p至原先一半；若s與p的比值大于某個數(shù)值，s也減半。Ridge方法的缺點是接近過渡態(tài)時效率較低，可以當C進入過渡態(tài)二次區(qū)域后改用QN法來加快收斂。也可以結合DIIS法，速度比原先有一半以上的提升，效率有時還高于基于二階導數(shù)的方法，而且在某些勢能面非常平坦的體系比二階導數(shù)方法更可靠。[圖6]Ridge方法示意圖Step-and-Slide方法使產(chǎn)物和反應物的結構同時順著LST描述的路徑相對移動(step步)，直到它們的能量都等于某個預先設定的能量，然后讓這兩個結構在它們當前所在的勢能等值面上滑動(slide步)，使二者結構在坐標空間中的距離最小。重復上述step和slide步驟，最終當兩個結構碰上，這個位置就是過渡態(tài)。[圖7]StepandSlide方法示意圖Muller-Brown方法見下文IRC算法相應部分CI-NEB、ANEBA方法見下文“尋找過渡態(tài)的chain-of-state方法”相應部分基于反應物結構的算法最緩上升法（leaststeepascent,shallowestascent）由反應物結構到達過渡態(tài)結構的過程是沿著勢能面最容易行進的路徑進行的（不考慮動力學問題），這個途徑一般比其它方向要緩和，所以由反應物結構開始，沿著勢能面最緩的方向逐漸往上爬，往往可以沿著MEP到達過渡態(tài)。但要注意這條路徑時常與從過渡態(tài)沿最陡下降路徑所走出的MEP并不一致，因此原理上此法不能保證一定能到達過渡態(tài)。圖8描述的是LEPS勢結合諧振勢的勢能面，最緩上升法所走的黑色粗曲線嚴重不符合實際MEP（黑點所示路徑），而且曲線是中斷的。此法也可能走到與此平衡結構相連的其它過渡態(tài)，而非預期的過渡態(tài)。還容易因為步長問題導致走到中途時跑到另外一條錯誤路徑上，雖然設小步長能得到解決，但是需要花費更長時間。因為種種問題，這個方法使用較少。［圖8］勢能面上最緩上升法所走的路徑（黑色粗曲線）本征向量/本征值跟蹤法（eigenvector/eigenvaluefollowing,EF。也稱modewalking/modefollowing/Walkingupvalleys）由于平衡結構越過勢壘發(fā)生反應的能量主要來自分子某振動模式提供的動能，考慮這一點，由平衡結構沿著此振動矢量方向步進，能夠找到過渡態(tài)，經(jīng)歷的路徑就是反應路徑。這種方法需要首先對平衡結構進行振動分析，由用戶最初指定一個可能指向過渡態(tài)的振動模式。因為平衡態(tài)通向過渡態(tài)路徑勢能面平緩，曲率（可視為振子力常數(shù)）一般小于其它方向，故一般跟蹤頻率最低的振動模式（高斯中默認）。每走一步后重新計算Hessian矩陣的本征值和本征向量，如果跟蹤的是本征值最低的模式，仍取本征值最小的本征向量繼續(xù)跟蹤；如果跟蹤的是其它振動模式，就取與上一步所跟蹤的向量重疊最大的向量繼續(xù)跟蹤。重復執(zhí)行，直到符合收斂標準為止。如果一個結構涉及到多個過渡態(tài)，則跟蹤不同的本征向量有可能得到不同的過渡態(tài)，即便所跟蹤的不是最低模式，當接近過渡態(tài)后也會成為最低的模式。此方法也可以直接由過渡態(tài)初猜結構開始跟蹤，或者說EF方法是一種不需要初猜在過渡態(tài)二次區(qū)域內的尋找過渡態(tài)的方法。由穩(wěn)定結構通過EF方法跟蹤至過渡態(tài)相對與直接給出初猜顯然更為費時，但對于不能預測過渡態(tài)結構的情況下往往是有用的。LMOD法搜索構象也是基于這一原理，不斷地根據(jù)低頻振動方向越過構象轉變的過渡態(tài)到達新的構象。最初的EF方法只是簡單地沿所跟蹤的振動模式移動來升高能量。高斯中opt=（EF,TS）方法還使結構同時在其余方向上沿能量更低的方向移動，其實它用的就已介紹的P-RF0法，所跟蹤的模式用獨立計算的入的最大解，其它的模式使用相同的另外計算的入的最小解。由于Berny方法尋找過渡態(tài)已經(jīng)包含了P-RF0步，所以EF方法實際上也已經(jīng)包含在內了，除非要用到跟蹤特定模式等功能時才有使用的必要。ARTn（activation-relaxationtechniquenouveau）此方法主要用于研究無序材料的在能量面上由極小點穿過過渡態(tài)到達其它極小點的過程，解決由于勢壘高而難以用MD和MC方法研究的問題。方法分兩步，（1）將初始結構由極小點位置激活并收斂到過渡態(tài)(activation步),(2)由過渡態(tài)通過常規(guī)的能量極小化算法尋找極小點(relaxation步)。(1)中的每一步中在任意方向上移動結構，然后在垂直于走過的路徑方向的超平面上做能量極小化，反復執(zhí)行，直到Hessian矩陣出現(xiàn)一個負本征值為止。之后進入收斂至鞍點的步驟，在最小本征值的方向上沿受力反方向移動，其余方向根據(jù)受力移動，最終將找到一階鞍點。由于大體系Hessian矩陣本征值求解困難，此方法中使用Lanczos算法快速求解最低本征值和本征向量。ART法可以獲得與初始極小點相連的許多過渡態(tài)。梯度極值法(Gradientextremal,GE)梯度極值路徑連接的是每一個等值線(高維情況為超曲面)上的梯度的模|g|為極大或極小值的點(相對于同一等值線上的其它點的梯度模來說)。因為勢能面的每一點的梯度垂直于此點等值線的切線，故梯度模極值點的位置相當于垂直于等值線方向上等值線間隔比處在相同等值線上相鄰的點更遠或更近。|g|的極值與g'2—致，設勢能函數(shù)為f,限制所在等值線能量為k，通過拉格朗日乘子法求g"2的極值▽[g"2-2入(f-k)]=O，可知梯度極值點的梯度方向等于此點Hessian矩陣某一本征向量。由于勢能面上每個駐點必有一條或多條梯度極值路徑通過而互相構成網(wǎng)絡(但任意駐點間不一定有梯度極值路徑直接相連)，所以系統(tǒng)地跟蹤梯度極值路徑是一種獲得勢能面上全部駐點的方法，目前已有幾種跟蹤算法，然而即便對于簡單體系，梯度極值路徑數(shù)目也極多，尤其是包含對稱性情況下。由極小點跟蹤梯度極值路徑也能夠用于尋找過渡態(tài)，但極小點未必與過渡態(tài)通過梯度極值路徑直連，且此方法并不能控制要尋找哪類駐點，故為了尋找過渡態(tài)可能需要從多個其它駐點跟蹤多個梯度極值路徑，計算量很大，所以單純?yōu)榱藢ふ疫^渡態(tài)而使用這種方法不切實際。[圖9]梯度極值路徑示意圖約化梯度跟蹤(reducedgradientfollowing,RGF)這個方法同梯度極值法一樣可以得到包括過渡態(tài)、極小點在內的各種駐點。設勢能面為N維，此方法將跟蹤N條路徑，其中第i條(i=1,2,3...N)路徑只有在第i維上梯度不為0,而其它N-1個維度上皆為0,故稱為約化梯度。這樣的路徑交匯的位置，就是所有維度上梯度皆為0的位置，即駐點。例如簡單的二維情況E(x,y)=x'3+y'3-6xy,跟蹤的RGF方程就是Ex(x,y)=3x"2-6y=0和Ey(x,y)=3y"2-6x=0,前者僅y方向梯度不為0,后者僅x方向梯度不為0，相交得到的駐點為一個一階鞍點和一個極小點。也可以使用原始坐標組合的正交坐標系，例如跟蹤僅x+y和僅x-y方向上梯度不為0的兩條路徑。[圖10]x"3+y"3-6xy面上約化梯度路徑示意圖跟蹤約化梯度的步進算法是第m點的坐標x(m+1)=x(m)+StL*x'(m)/|x'(m)|。StL是步長，x'(m)/|x'(m)|代表路徑切線方向單位向量。x'可以通過H'x'=0方程以QR分解法獲得，其中H'與Hessian矩陣唯一不同的是，若當前跟蹤的是僅第k維梯度不為0的約化梯度路徑，則H'沒有Hessian矩陣的第k行。一般起始步由某駐點開始，此步準確計算Hessian,步進過程中Hessian可用前述的DFP方法修正。每步檢驗所跟蹤方向上的朝向下一個駐點的牛頓步步長，若小于標準則停止，并且再精確計算一次Hessian以確認此駐點是什么類型。每次走步的結果如果在數(shù)值上與“僅某維度上梯度為0”條件符合較好，可以動態(tài)增加步長，類似AH法的置信半徑概念，如果相差較大，則調用校正步(后期方法將校正步合并入步進步，改善了效率和穩(wěn)定性)。這個方法計算量也很大，而且也無法指定要搜索的駐點的類型，所以不適合獨立用作尋找過渡態(tài)。等勢面搜索法(IsopotenialSearching)如果將反應物位置附近的勢能面比做一個湖，這個方法可以看作逐漸往湖里面灌水，由于過渡態(tài)能量比周圍地方更低，所以隨著水位(勢能)逐漸升高，水最先流出來的地方就是過渡態(tài)。繼續(xù)灌水，隨著水位繼續(xù)升高，還可以找到其它能量更高的過渡態(tài)。具體實現(xiàn)的方法是：首先最小化反應物的能量E0,在反應物位置附近設置一些測試點，可以隨機也可以根據(jù)經(jīng)驗設定，作為“水位”來檢測是否已到達過渡態(tài)能量。然后設定目標能量E(target),一般高于E0幾百KJ/mol。計算那些測試點的能量和勢能梯度，檢查其能量與E(target)的差的絕對值，若大于10KJ/mol,即沒達到目標水位，就讓它們沿著梯度方向行進以提升能量，之后再次檢查是否符合條件，直到小于10KJ/mol，即已到達目標水位，就對這些點進行人工的檢查，包括結構、成鍵分析等，考察在E(target)時是否已經(jīng)達到或超過了過渡態(tài)的能量。如果找到了過渡態(tài)，就調整這些點的位置繼續(xù)找別的過渡態(tài)；如果未找到，就提高E(target)并且調測試點整位置以增大找到過渡態(tài)的概率，然后再沿著梯度方向提升測試點的能量并進行接下來的檢測，反復如此。上述提到的“調整點的位置”有很多算法，但主要都是使那些測試點在垂直于梯度，即在等值面上移動。因為測試點無法密集覆蓋整個等勢面，受計算能力制約其數(shù)目有限的，很難有哪個點隨著E(target)的提升而移動后恰好落在過渡態(tài)的位置。直到E(target)提升到有測試點可判斷為過渡態(tài)時，其能量一般已高出實際過渡態(tài)很多。所以使用此方法得到的過渡態(tài)能量與初始點位置和調整點位置的算法都有很大關系，一般都顯著偏高，甚至不能找到過渡態(tài)，可嘗試以不同初始位置和調整算法重新執(zhí)行以改善結果。等勢面搜索法適合在只有反應物結構而難以預測過渡態(tài)和產(chǎn)物結構的情況下尋找過渡態(tài)，例如預測質譜中分子的可能裂解的方式，有時還可能找到全新未曾考慮到的反應機理。但是此方法的結果很粗糙，而且計算量極大，尤其是大分子的高維勢能面，有限的測試點很容易漏掉許多重要過渡態(tài)。球形優(yōu)化(Sphereoptimization)在幾何參數(shù)的變量空間上，以反應物或產(chǎn)物為中心，在不斷增加半徑的超球面上做能量最小化。將相鄰球面上得到的能量極小點相連接，就得到一條由反應物或產(chǎn)物為起點的低能量的路徑，可做為IRC(未必正確，考慮圖8的勢能面)，并由此找到過渡態(tài)。如果每個球面上可以找到多個極小點，則連接后有可能得到多條反應路徑。此法若以坐標驅動法為類比，此方法就是對幾何參數(shù)空間中反應物或產(chǎn)物結構代表的點的距離進行柔性掃描。［圖11］球形優(yōu)化示意圖全勢能面掃描當一切方法都不能找到過渡態(tài)，全勢能面掃描是最終途徑。由于掃描得到的勢能面格點是離散的，可通過插值提高格點密度以增加精度。得到勢能面后，就可以通過一些算法找到過渡態(tài)，例如用這些點擬合出解析表達式，然后用標準微分方法找過渡態(tài)。但全勢能面掃描極為昂貴，內坐標下需要計算X'(3N-6)次住代表每個變量掃描步數(shù))，只限于反應中僅涉及幾個自由度的勢能面掃描，往往不得不考慮更低級的方法如半經(jīng)驗或者分子力學，變量稍多的體系則完全不能實現(xiàn)。全勢能面掃描的結果還提供了過渡態(tài)位置以外結構的信息，例如可以用于研究反應路徑、用于構象搜索等。過渡態(tài)相關問題無過渡態(tài)的反應途徑(barrierlessreactionpathways)并非所有反應途徑都需要越過勢壘，這類反應在很低的溫度下就能發(fā)生，盲目找它們的過渡態(tài)是徒勞的。常見的包括自由基結合，比如甲基自由基結合為乙烷；自由基向烯烴加成，比如甲基自由基向乙烯加成成為丙基自由基；氣相離子向中性分子加成，比如叔碳陽離子向丙烯加成。等等。Hammond-Leffler假設過渡態(tài)在結構上一般會偏向反應物或者產(chǎn)物結構一邊。Hammond-Leffler假設對預測過渡態(tài)結構往哪個方向偏是很有用的，意思是反應過程中，如果兩個結構的能量差異不大，則它們的構型差異也不大。由此可知對于放熱反應，因為過渡態(tài)能量與反應物差異小，與產(chǎn)物差異大，故過渡態(tài)結構更偏向反應物，相反，吸熱反應的過渡態(tài)結構更偏向產(chǎn)物。所以初猜過渡態(tài)結構應考慮這一問題。對稱性問題如果已經(jīng)明確地知道過渡態(tài)是什么對稱性，而且對稱性高于平衡態(tài)對稱性，且可以確信在這個高對稱性下過渡態(tài)是能量最低點，則可以強行限制到這個對稱性之后進行幾何優(yōu)化，幾何優(yōu)化算法比尋找過渡態(tài)算法方法更可靠。比如F+CH3F—〉FCH3+F這個SN2反應，過渡態(tài)就是傘形翻轉的一刻，恰為高對稱性的D3h點群，而反應路徑上的其它結構對稱性都比它低，所以在D3h點群條件下優(yōu)化，得到的能量最低點就是過渡態(tài)。如果過渡態(tài)對稱性不確定，則找過渡態(tài)計算的時候不宜設任何對稱性，否則若默認保持了平衡態(tài)下的對稱性，得到的此對稱下的過渡態(tài)并不是真正的過渡態(tài)，容易得到二階或高階鞍點。溶劑效應計算凝聚態(tài)條件下過渡態(tài)的性質，必須考慮溶劑效應，它明顯改變了勢能面。一般對過渡態(tài)的結構影響較小，但對能量影響很大。有時溶劑效應也會改變反應途徑，或產(chǎn)生氣相條件下沒有的勢壘。溶劑條件下，上述尋找過渡態(tài)的方法依然適用。應注意涉及到與溶劑產(chǎn)生氫鍵等強相互作用的情況，隱式溶劑模型是不適合的，需要用顯式溶劑考察它對過渡態(tài)的影響，即在輸入文件中明確表達出溶劑分子。計算過渡態(tài)的建議流程直接用高水平方法計算過渡態(tài)往往比較花時間，可以使用逐漸提高方法等級的方法加速這一過程，一般建議是：執(zhí)行低水平的計算找過渡態(tài)，如半經(jīng)驗。將第1步得到的過渡態(tài)作為初猜，用高級別的方法找過渡態(tài)。在相同水平下對上一步找到的過渡態(tài)做振動分析，檢驗是否僅有一個虛頻，以及觀看其振動模式的動畫來考察振動方向是否連接反應物與產(chǎn)物結構。有必要時可以做IRC進一步檢驗。為獲得更精確的過渡態(tài)能量，可使用更高等級方法比如含電子相關的方法計算能量。內稟反應坐標(intrinsicreactioncoordinate,IRC)MEP指的是勢能面上，由一個點到達另一個點的能量最低的路徑，滿足最小作用原理。若質量權重坐標下的MEP連接的是反應物、過渡結構和產(chǎn)物，則稱為IRC。所謂質權坐標在笛卡兒坐標下即r(i,x)=sqrt(m(i))*R(i,x),m(i)為i原子質量，R(i,x)為i原子原始x方向坐標，同樣有r(i,y)、r(i,z)。IRC描述了原子核運動速度為無限小時，質權坐標下由過渡態(tài)沿著勢能負梯度方向行進的路徑(最陡下降路徑)，其中每一點的負梯度方向就是此處核的運動方向，在垂直于路徑方向上是能量極小點。注意質量權重和非權重坐標下的路徑是不一樣的。IRC可看作0K時的實際在化學反應中原子核所走的路徑，溫度較低時IRC也是一個很好的近似。但是當溫度較高，即核動能較大時，實際反應路徑將明顯偏離IRC，而趨于沿最短路徑變化，即便經(jīng)歷的是勢能面上能量較高的的路徑，這時就需要以動力學計算的平均軌跡來表征反應路徑。算法最陡下降法(Steepestdescent)最簡單的獲得IRC的方法就是固定步長的最陡下降法，由過渡態(tài)位置開始，每步沿著當前梯度方向行進一定距離直到反應物/產(chǎn)物位置，也稱Euler法。由于最陡下降法及下文的IMK、GS等方法第一步需要梯度，而過渡態(tài)位置梯度為0,所以第一步移動的方向沿著虛頻方向。最陡下降方法與IRC的本質相符，但是此法實際得到的路徑是一條在真實IRC附近反復震蕩的曲折路徑，而非應有的平滑路徑，對IRC描述不夠精確。雖然可以通過更小的步長得以一定程度的解決，但是太花時間，對于復雜的反應機理，需要更多的點。也可以通過RK4(四階Runga-Kutta)來走步，比上面的方法更穩(wěn)定、準確，但每步要需要算四個梯度，比較費時。IMK方法(Ishida-Morokuma-Kormornicki)它是最陡下降法的改進，解決其震蕩問題。首先計算起始點X(k)的梯度g(k)，獲得輔助點X'(k+l)=X(k)-g(k)*s,其中s為可調參數(shù)。然后計算此點梯度g'(k+l),在g(k)與-g'(k+l)方向的平分線上(紅線所示)進行線搜索，所得能量最小點即為X(k+1)，之后再將X(k+1)作為上述步驟的X(k)重復進行。整個過程類似先做最陡下降法，然后做校正。此方法仍然需要相對較小的步長，獲得較精確IRC所需計算的點數(shù)較多。[圖12]IMK方法示意圖Schmidt,Gordon，Dupuis改進了IMK的三個細節(jié)，使之更有效率、更穩(wěn)定。(1)將X'(k+1)的確定方式改為了X(k)-g(k)/|g(k)|*s，即每一步在負梯度方向上行進固定的s距離，與梯度大小不再有關。(2)線搜索步只需在平分線上額外計算一個點的能量即可，這個點和X'(k+1)點的能量以及g'(k+1)在此平分線上的投影三個條件作聯(lián)立方程即可解出曲線方程，減少了計算量。IMK原始方法則需要在平分線上額外計算兩個點的能量與X'(k+1)的能量一起擬和曲線方程。(3)第一步在過渡態(tài)位置的移動距離Aq如此確定:AE=k*(Aq"2)/2,k為虛頻對應的力常數(shù)，AE為降低能量的期望值(一般為hartree)，這樣可避免在虛頻很大的鞍點處第一步位移使能量降低過多。Muller-Brown方法這是通過球形限制性優(yōu)化找IRC的方法。首先將過渡態(tài)和能量極小點位置定義為P1和P2,由P1開始步進,當前步結構以Q(n)表示。每一步,在相距Q(n)為r距離的超球面上用simplex法優(yōu)化獲得能量極小點Q'(圖中綠點)，優(yōu)化的起始點是Q(n-1)Q(n)與Q(n)P2方向的平分線b上距Q(n)為r距離的位置S(紅點)。若Q(n)Q'與Q(n)P2的夾角較小，則Q'可當作是下一步位置Q(n+1)。如此反復，直到符合停止標準，比如下一步能量比當前更高(已走過頭了)、與P2距離已很近(如小于)、或者與P2方向偏離太大(P1與P2點通過此法無法找到IRC)。最終所得到全部結構點依次相連即為近似的IRC,減小步長r值可使結果更貼近實際IRC?；诖朔椒ㄒ部梢杂糜趯ふ疫^渡態(tài)，先將反應物和產(chǎn)物作為P1和P2,將二者距離的約2/3作為r,由其中一點在P1-P2連線上相距其r位置為初始位置進行球形優(yōu)化得到0點，在0與P1、0與P2上也如此獲得P1'與P2',根據(jù)P1、P1'、0、P2'、P2的能量及之間距離信息以一定規(guī)則確定其中哪兩個點作為下一步的P1和P2,確定新的P1和P2后重復上述步驟，直至P1與P2十分接近，即是過渡態(tài)。此方法計算IRC可以步長可設得稍大，第一步不需要費時的Hessian矩陣確定移動方向，缺點是獲得的路徑曲率容易有問題，對于曲率較大的反應路徑需要減小步長。[圖13]Muller-Brown方法示意圖GS(Gonzalez-Schlegel)方法這是目前很常用，也是Gaussian使用的方法，見圖14。首先計算起始點X(k)的梯度，沿其負方向行進s/2距離得到X'(k+1)點作為輔助點。在距X'(k+1)點距離為s/2的超球面上做限制性能量最小化，找到下一個點X(k+1)。因為這個點的負梯度(黑色箭頭)在弧方向上分量為0,故垂直于弧，即其梯度方向在X'(k+1)到X(k+1)的直線上。這必然可以得到一段用于描述IRC的圓弧(虛線)，它通過X(k)與X(K+1)點，且在此二點處圓弧的切線等于它們的梯度方向，這與IRC的特點一致，這段圓弧可以較好地(實線)。之后再將X(k+1)作為上述步驟的X(k)重復進行。GS方法對IRC描述得比較精確，在研究反應過程等問題中，由于對中間體結構精度有要求,GS是很好的選擇，而且用大步長可以得到與小步長相近的結果，優(yōu)于IMK、MUller-Brown等方法。若只想得到與過渡態(tài)相連的反應物和產(chǎn)物結構，或者粗略驗證預期的反應路徑，對IRC精度要求不咼，使用最陡下降法往往效率更咼，盡管GS可以用更大步長，但每步更花時間。[圖14]GS方法示意圖除上述外，IRC也可以通過已提及的EF、最緩上升法、球形優(yōu)化等方法得到，它們的好處是不需要事先知道過渡態(tài)的結構。贋坐標法除了簡單的反應以外，只能得到近似的IRC,由于結構的較小偏差會帶來能量的較大變化，容易引入滯后效應，所以這樣得到的勢能曲線難以說明問題。6.chain-of-states方法這類方法主要好處是只需要提供反應物和產(chǎn)物結構就能得到準確的反應路徑和過渡態(tài)。首先在二者結構之間以類似LST的方式線性、均勻地插入一批新的結構(使用內坐標更為適宜),一般為5~40個，每個結構就是勢能面上的一個點(稱為image)，并將相鄰的點以某種勢函數(shù)相連，這樣它們在勢能面上就如同組成了一條鏈子。對這些點在某些限制條件下優(yōu)化后，在勢能面上的分布描述的就是MEP,能量最高的結構就是近似的過渡態(tài)位置。Dragmethod方法這個方法最簡單，并不是嚴格的chain-of-states方法，因為每個結構點是獨立的。插入的結構所代表的點均勻分布在圖8所示的短虛線上，也可以在過渡態(tài)附近位置增加點的密度。每個點都在垂直于短虛線的超平面上優(yōu)化，在圖中就是指平行于長虛線方向優(yōu)化。這種方法一般是奏效的，但也很容易失效，圖8就是一例，優(yōu)化后點的分布近似于從產(chǎn)物和反應物用最緩上升法得到的路徑(黑色粗曲線)，不僅反應路徑錯誤，而且兩段不連接，與黑色小點所示的真實MEP相距甚遠(黑色點是用下文的NEB方法得到的)。目前基本不使用此方法。PEB方法(plainelasticband)這是下述Chain-of-state方法的基本形式。也是在反應物到產(chǎn)物之間插入一系列結構，共插入P-1個，反應物編號為0,產(chǎn)編號物為P。不同的是優(yōu)化不是對每個點孤立地優(yōu)化，而是優(yōu)化一個函數(shù)，每一步所有點一起運動。下文用工[i=l,P]X(i)符號代表由X(l)開始加和直到X(P)。PEB函數(shù)是這樣的：S(R(1),R⑵…R(P-1))=工[i=l,P-l]V(R(i))+工[i=l,P](k/2*(R(i)-R(i-l)廠2)。其中R(i)代表第i個點的勢能面上的坐標，V(R(i))是R(i)點的能量，k代表力常數(shù)。優(yōu)化過程中反應物R(0)和產(chǎn)物R(P)結構保持不變，優(yōu)化此函數(shù)相當于對一個N*(P-2)個原子的整體進行優(yōu)化，N為體系原子數(shù)。優(yōu)化過程中，式中的第一項目的是讓每個點盡量向著能量極小的位置移動。第二項相當于將相鄰點之間用自然長度為0、力常數(shù)為k的彈簧勢連了起來，目的是保持優(yōu)化中相鄰點之間距離均衡，避免過大。當只有第一項的時候，函數(shù)優(yōu)化后結構點都會跑到作為能量極小點的反應物和產(chǎn)物位置上去而無法描述MEP,這時必然會有一對兒相鄰結構點距離很大。當?shù)诙棾霈F(xiàn)后，由于此種情況下彈簧勢能很高,在優(yōu)化中不可能出現(xiàn)，從而避免了這個問題。dragmethod法在圖8中失敗的例子中，也有一對兒相鄰結構點距離太遠，所以也不會在PEB方法中出現(xiàn)。簡單來說，PEB方法就是保持相鄰結構點的間距盡量小的情況下，優(yōu)化每個結構點位置?？梢越票扔鞒稍趧菽苊娴哪Ｐ蜕?，將一串以彈簧相連的珠子，一邊掛在反應物位置，另一邊掛在產(chǎn)物位置，拉直之后松手，這串珠子受重力作用在模型上滾動，停下來后其形狀可當作MEP，最高的位置近似為過渡態(tài)。但是PEB方法的結果并不能很好描述MEP。圖15描述的是常見的A、B、C三原子反應的LEPS勢能面，B可與A或C成鍵，黑色弧線為NEB方法得到的較真實的MEP。左圖中，在過渡態(tài)附近PEB的結構點沒有貼近MEP，得到的過渡態(tài)能量過高，稱為corner-cutting問題。這是因為每點間的彈簧勢使這串珠子僵硬、不易彎曲，由圖15右圖可見，⑴朝R(i-1)與R(i+1)方向都會受到彈簧拉力，其合力牽引R(i)，使R(i-1)、R(i)、R(i+1)的弧度有減小趨勢。如果將彈簧力常數(shù)減小以減弱其效果，就會出現(xiàn)圖15中間的情況，雖然結構點貼近了MEP,但相鄰點間距沒有得到保持，過渡態(tài)附近解析度很低，錯過了真實過渡態(tài)，若以能量最高點作為過渡態(tài)則能量偏低，這稱為sliding-down問題?？梢姀椈闪Τ?shù)k的設定對PEB結果有很大影響，為權衡這兩個問題只能取折中的k，但結果仍不準確。[圖15]LEPS勢能面上不同k值的PEB結果Elber-Karplus方法與PEB函數(shù)定義相似。第一項定義為1/L*工[i=1,P-1](V(R(i))*d(i,i-1))，其中L為鏈子由0點到P-1點的總長，d(i,i+1)為R(i)與R(i+1)的距離，此項可視為所有插入點總能量除以點數(shù)，即插入點的平均能量。第二項為Y*E[i=1,P](d(i,i-1)-<d>)"2，其中〈d〉代表相鄰點的平均距離，是所有d(i,j)的RMS。此項相當于將彈簧自然長度設為了當前各個彈簧長度的平均值，由Y參數(shù)控制d(i,j)在平均值上下允許的波動的范圍。此方法最初被用于研究蛋白質體系的構象變化。SPW方法(Self-PenaltyWalk)在Elber-Karplus方法的基礎上增加了第三項互斥項，工[i=0,P-1]工[i二j+1,P-1]U(ij)，其中U(ij)=P*exp(-d(i,j)/(入*心〉))，<d〉定義同上。此項相當于全部點之間的“非鍵作用能U(ij)”之和，不再僅僅是相鄰點之間才有限制勢。任何點之間靠近都會造成能量升高，可以避免Elber-Karplus方法中出現(xiàn)的在能量極小點處結構點聚集、路徑自身交錯的問題，能夠使路徑充分地展開，確保過渡態(tài)區(qū)域有充足的采樣點。式中P和入都是可調參數(shù)來設定權重。此外相對與Elber-Karplus方法還考慮了笛卡兒坐標下投影掉整體運動的問題。LUP方法(LocallyUpdatedplanes)特點是優(yōu)化過程中，只允許每個結構點R(i)在垂直于R(i-1)R(i+1)向量的超平面上運動。由于每步優(yōu)化后R(i-1)與R(i+1)連線方向也會變化，故每隔一定步數(shù)重新計算這些向量，重新確定每個點允許移動的超平面。但是LUP缺點是結構點之間沒有以上述彈簧勢函數(shù)相連來保持間隔，容易造成結構點在路徑上分布不均勻，甚至不連續(xù)，還可能逐漸收斂至兩端的極小點。NEB方法(NudgedElasticBand)NEB方法集合了LUP與PEB方法的優(yōu)點，其函數(shù)形式基于PEB。從PEB方法的討論可以看出，彈簧勢是必須的，它平行于路徑切線(R(i)-R(i-1)與R(i+1)-R(i)矢量和的方向)的分量保證結構點均勻分布在MEP上來描述它；但其垂直于路徑的分量造成的弊端也很明顯，它改變了這個方向的實際的勢能面，優(yōu)化后得到的MEP'就與真實的MEP發(fā)生了偏差，造成corner-cutting問題。解決這個問題很簡單，在NEB中稱為nudge過程，即每個點在平行于路徑切線上的受力只等于彈簧力在這個方向分量，每個點在垂直于路徑切線方向的受力只等于勢能力在此方向上分量。這樣彈簧力垂直于路徑的分量就被投影掉了，而有用的平行于路徑的分量完全保留；勢能力在路徑方向上的分量也不會再對結構點分布的均勻性產(chǎn)生影響，被保留的它在垂直于路徑上的分量將會引導結構點地正確移動。這樣優(yōu)化收斂后結構點就能正確描述真實的MEP,矛盾得到解決。彈簧力常數(shù)的設定也比較隨意，不會再對結果產(chǎn)生明顯影響。但是當平行于路徑方向能量變化較快，垂直方向回復力較小的情況，NEB得到的路徑容易出現(xiàn)曲折，收斂也較慢，解決這一問題可以引入開關函數(shù)，即某點與兩個相鄰點之間形成的夾角越小，此點就引入更多的彈簧勢垂直于路徑的分量，使路徑不易彎曲而變得光滑，但也會帶來一定corner-cutting問題。也可以通過將路徑切線定義為每個點指向能量更高的相鄰點的方向來解決。DNEB方法（DoubleNudgedElasticBand）彈簧勢垂直于路徑的分量壞處是造成corner-cutting問題，好處是避免路徑卷曲。更具體來說，前者是由于它平行于勢能梯度方向的那個分量造成的，若只將這個分量投影掉，就可避免corner-cutting問題，而其余分量的力F（DNEB）仍可以避免路徑卷曲，這便是DNEB的主要思想。故DNEB與NEB的不同點就是DNEB保留了彈簧勢垂直于路徑的分量其中的垂直于勢能梯度的分量。DNEB的這個設定卻導致結構點不能精確收斂到MEP上。正確的MEP上的點在垂直于路徑方向上受勢能力一定為0,但是當用了DNEB方法后，若其中某一點處路徑是彎曲的，即彈簧力在垂直于路徑方向上有分量F而且此點勢能梯度方向不垂直于此點處路徑的切線，即F'不會被完全投影掉，F(xiàn)'力的分量F（DNEB）將繼續(xù)帶著這個點移動，也就是說結構點就不在正確的MEP上了。只有當結構點所處路徑恰為直線，即F'為0則不會有此問題。為了解決此問題有人將開關函數(shù)加入到DNEB，稱為swDNEB，當結果越接近收斂，即垂直于路徑的勢能力越小的時候，F(xiàn)（DNEB）也越小，以免它使結構點偏