量子化學計算中幫助幾何優化收斂的常用方法
筆者注：讀者不要說諸如“此文我看了，還是不收斂”、“我用了文中的方法還是不收斂”這類話。絕對沒有哪個方法能100%解決幾何優化不收斂。本文絕對沒有一條餿主意，所有Gaussian支持的真正有可能能解決幾何優化不收斂的做法在本文都已經全面列舉了，沒有任何遺漏。倘若本文的做法都做了嘗試還沒解決，那也別指望有任何其它方法能解決。才隨便試了文中一個做法，或者以不正當的方式（沒有真正了解原理的基礎上）瞎試幾個就說沒解決問題，這根本什么也說明不了。要想解決問題，必須十分仔細閱讀本文，所有有可能對當前問題有用的方法都依次嘗試，有必要的時候幾個一起結合使用。本文列舉的做法已經把最壞的情況都考慮了，所以出現所謂的“看了本文還是沒解決問題”的情況一定是讀者還沒真正仔細看、仔細領會、充分嘗試。
量子化學計算中幫助幾何優化收斂的常用方法
文/Sobereva @北京科音
First release: 2012-Oct-13 Last update: 2022-May-5
0 前言
幾何優化就是尋找勢能面極小點結構的過程。而所謂幾何優化不收斂，就是優化過程始終，或者很難達到收斂要求，這種情況通常會伴隨著震蕩行為，即原子受力、幾何結構隨優化步數增加呈現一定周期性變化趨勢。顯然盲目增加優化步數的上限試圖解決是超級愚蠢的做法。解決這種問題必須在結合經驗和理論知識的前提下，通過考察實際收斂的趨勢，嘗試各種可能奏效處理辦法，本文就專門說說這個問題。更詳細、更具體的東西在筆者講授的北京科音（http://www.keinsci.com）的量子化學培訓班里會系統地講解。
經常有人問Gaussian程序做幾何優化為什么有時會出現下面這種報錯
Error termination request processed by link 9999.
Error termination via Lnk1e in d:\study\G09W\l9999.exe at Sat Jun 02 20:45:20 2018.
實際上往上看，就會看到更明確的錯誤提示
Optimization stopped.
-- Number of steps exceeded, NStep= [當前步數上限]
-- Flag reset to prevent archiving.
即達到了步數上限還沒收斂，于是優化就報錯結束了。這多數情況就是因為出現了上述震蕩問題而造成了難收斂所導致的。另外，還經常有初學者問，怎么優化任務算了非常非常久還沒算完，這很可能也是因為發生了震蕩，導致優化始終達不到收斂限所致（比如都跑了好幾百步了，還在震蕩中）。
如果想確認是否發生了震蕩，可以用GaussView載入輸出文件（打開文件的界面下方必須勾選上Read intermediate geometries），然后進入Results - Optimization，看能量和受力變化曲線。如果曲線到后期沒有整體降低、逐漸收斂到平坦的大趨勢，而是來回上下波動，并且觀看優化軌跡也發現結構在反復波動，這就說明震蕩了。幾何優化正在執行的過程中也可以隨時對收斂情況這樣進行監視。
說怎么解決幾何優化難收斂之前先說一下收斂標準。Gaussian中判斷幾何優化收斂有四個標準，在默認收斂設定下，這四個標準是：
最大受力<0.00045；方均根受力<0.00030；最大位移<0.00180；方均根位移<0.00120
當這四個標準都滿足了就宣告收斂。方均根受力/方均根位移體現的是體系中所有原子的平均受力/位移情況。另外，優化過程中只要受力小于預定的收斂限100倍，哪怕位移還沒低于收斂限，也算作已收斂。這主要考慮到勢能面非常非常緩的大的柔性分子，相對于這樣尺度的分子，幾何結構收斂到那么精確意義不大，放寬位移收斂限避免了收斂太慢。
也注意有時候優化任務中途出錯，不是因為幾何收斂方面的問題，而是因為優化當中某一步連能量計算都沒能正常完成，很常見的情況是由于SCF不收斂導致的，這種情況的解決見《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61）。幾何優化過程也可能朝著明顯錯誤的方向進行而導致難以收斂，這可能是理論方法、基組、電子態、計算模型的構建，以及其它某些設定不合理所致。這些方面和優化不收斂問題本身沒關系，所以不會在本文提到。
如果在原先設置下的優化任務最終沒能收斂，可以取優化過程中受力或能量最低的一個結構（可以用gview切換到相應的幀，保存成新的輸入文件。也可以用geom(check,step=n)從%chk指定的chk文件中讀取第n步的坐標），然后加上下文提到的關鍵詞接著進行優化。
下面，就詳細談談碰到幾何優化不收斂（通常由于上述震蕩問題導致）應當怎么解決。本文列舉各種有可能解決不收斂，也包括加速收斂的辦法。其中很多方法可以相互結合使用以達到更好的效果。本文雖然說的是Gaussian的情況，但很多方法在其它程序如ORCA中也可以類似地使用。對于優化過渡態過程出現難收斂的情況，在本文末尾有專門的提及。
1 嘗試不同的優化方法
優化幾何結構的方法有很多，以前我在《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）當中詳細介紹過的很多搜索過渡態的方法其實和搜索勢能面極小點（即幾何優化）的方法本質是一致的。對于從頭算，在Gaussian03中默認的是Berny方法，這方法實際上相當于在一般的RFO方法上添油加醋，而RFO法可理解為是Newton-Raphson方法（一種很常用的非線性優化方法）的改良。G03也支持GDIIS方法（改進版，與RFO方法相結合），它對半經驗方法優化是默認的，對于大體系、勢能面較平坦的情況收斂得一般比RFO要快（勢能面平坦表現出的特征是每一步受力小但位移大）。在G09中新支持了GEDIIS方法，并且作為了默認的優化方法，號稱是目前效率最高的優化方法，自稱比RFO或GDIIS在各種條件下都能更快收斂，但實際情況中往往并非如此。當使用某一種方法不收斂時，可以考慮用別的方法。這三種方法可以用opt=RFO、opt=GDIIS、opt=GEDIIS來選擇。個人強烈推薦對于弱相互作用體系、柔性大體系優化遇到不收斂時優先嘗試opt=GDIIS。
2 使用更好的Hessian矩陣
RFO、GDIIS、GEDIIS這三種優化方法在走步時都需要Hessian矩陣（質權坐標下也叫力常數矩陣），然而計算Hessian矩陣是很耗時的。因此，在默認情況下，Gaussian會通過價力場的方法近似估計出初始的Hessian矩陣，在每一步優化中只精確計算梯度，利用梯度對原先Hessian矩陣進行修正不斷得到新的Hessian矩陣。所以，在默認情況下，優化從頭到尾使用的Hessian矩陣都只是近似的（這也是為什么在優化出來的結構上做Freq往往會顯示還沒收斂，因為Freq用的是精確的Hessian矩陣）。當Hessian矩陣離精確值偏差較大，就會造成收斂緩慢，或者始終不收斂。為解決這個問題，可以用opt=calcfc，這會在優化的第一步使用精確計算的Hessian矩陣，但是仍可能后續優化過程中Hessian矩陣逐漸變得越來越不精確而依然收斂失敗。opt=calcall則不僅在第一步，在后續的每一步中也都精確計算Hessian矩陣，這使得很多優化失敗的情況都能得到解決，優化所需步數通常也會減少很多，而且能夠保證最終優化結果準確（因為最終判斷是否收斂時是基于精確的Hessian矩陣所得結果），但代價是每一步計算量會很大。
在G16中，增加了opt=recalc=n關鍵詞來設置每n步重新精確計算一次Hessian矩陣。顯然n越小，收斂的可能性越大，但耗時越高。n=1時對應calcall。實際中n建議設為3~5，以達到耗時和幫助收斂效果的權衡。
有時候，在當前計算級別下，哪怕只是calcfc都因為計算量太大（或程序限制）而難以做到，此時也可以考慮在更低級別下計算出Hessian矩陣作為高水平優化時的初始Hessian矩陣。具體來說，比如可以在低水平方法下做freq得到精確的Hessian矩陣，然后在高水平下優化時用opt=readfc來讀取它。
3 增加步數上限
這個方法是菜鳥最愛用的方法，但實際上幾乎不會起到用處！如果一個搞量化的碰見幾何優化不收斂時第一反應是使用這個辦法解決，那么他一定是個菜鳥！之所以這個做法這么流行，一方面是長期的以訛傳訛，另一方面是初學者們想當然地以為迭代次數越多自然越可能收斂。然而，絕大多數情況下，在默認步數下如果還不能收斂，那么多數情況是早已發生了震蕩，增大步數上限也無濟于事，只會白費功夫，糟蹋計算資源。
僅當初始結構離最終結構比較遠，特別是體系比較小（此時默認的步數上限較少）或者用了較小的步長上限時，且從gview顯示的優化過程受力和位移曲線上能看到結構變化和受力有整體降低趨勢，那么重算或續算時才有必要將步數上限設大，通過opt=maxcyc=N來指定上限為N步。注意實際上Gaussian程序在16版之前有內定步數上限（與坐標變量數正相關，也就是計算輸出中顯示的"maximum allowed number of steps"），即便N設得巨大，若超過了內定步數上限，N將自動設為內定步數上限值。
4 使用不同的坐標系和坐標定義方式
坐標系的選擇，以及坐標的定義方式，對于優化收斂速度乃至于能否收斂都有很大關系。各種坐標系各有所長，某種坐標系對于某個體系可能比較快就能收斂，但是換到另一個體系可能收斂會比較慢。當一種坐標系無法收斂時，可嘗試換另外的坐標系。一般來說，對于一個體系，哪種坐標設定下能夠使各個幾何變量之間耦合程度越弱，就越容易收斂。這容易理解，假設每個幾何變量之間都完全沒有耦合（Hessian矩陣非對角元為0），那么比如3N-6維下的幾何優化就等于3N-6個一維的優化，顯然一維優化是很容易完成的，所以整體優化也容易完成。如果各個坐標之間耦合比較強烈，每個坐標的微小改變又會導致其它坐標上受力的明顯改變，問題就比較復雜而難以拆分，收斂會變得相對困難，對優化算法的能力要求就越高。
在Gaussian中，支持以下三種坐標系下的優化。輸入文件里的坐標采用哪種坐標系記錄和實際優化中使用哪種坐標系完全無關。如果不注明優化過程所用的坐標系，都會自動生成冗余內坐標，并在冗余內坐標下優化。
(1) 笛卡爾坐標（Opt=Cartesian）：通常的優化不用笛卡爾坐標，因為對多數體系這種坐標下收斂都比在冗余內坐標下要慢。主要適合于含有較多分子的簇體系的優化。另外，當冗余內坐標下收斂過程因為坐標問題報錯時（比如構成二面角的四個原子中的三個在優化過程中排成了直線，此時無法定義二面角），改用此坐標系優化可以解決。注意，使用GDIIS時即便寫了cartesian，Gaussian依然會在默認的冗余內坐標下優化。
(2) 內坐標（Opt=Z-Matrix）：對于單分子體系，在Gview自動建立的內坐標系下優化一般比在笛卡爾坐標系下優化收斂會快。這不難理解。比如對應化學鍵長度的坐標之間，尤其是相隔較遠這樣的坐標之間，沒有明顯的耦合，一個鍵稍微伸縮一下不至于使得另一個鍵長坐標上的受力受到明顯影響。而對于笛卡爾坐標，一個原子朝笛卡爾坐標某個軸移動一下，這很可能使周圍很多原子受力發生不小的改變，因此耦合較大。但注意，如果內坐標建得很爛，胡亂定義，比如坐標跨著好幾個原子，則會造成耦合嚴重，反倒比笛卡爾坐標下收斂更慢。
內坐標不太適合有較多分子或原子的團簇體系，因為每個單體/原子之間內坐標怎么去鏈接比較合適往往說不清楚，很有可能會弄得很糟。內坐標尤為不適合環狀體系的優化（是指不利用虛原子進行特殊設定的情況下），這是因為內坐標對于環狀體系是以一個原子為起點繞一圈定義的，首尾雖然在幾何結構上是挨著的，但是在內坐標下卻沒有連著。這就導致了首尾原子間本該有的鍵長、鍵角、二面角項，要被其它各個內坐標所一起等效地“湊合”著描述，這勢必造成了那些內坐標之間強烈耦合。
內坐標對于高對稱性體系的優化特別有用。通過巧妙地自行設定要內坐標，并結合一些虛原子作為輔助，可以大大減少要優化的坐標數目，不僅易于收斂也顯著減少了計算梯度的耗時。對于一些昂貴的，而且只能通過有限差分方法計算梯度的理論方法來說，有時只有靠這種做法才能優化得動。
另外，Gaussian也支持內坐標和笛卡爾坐標的混合，可以結合二者所長。不過需要用到這種情況的場合不多，就不多提了。
(3) 冗余內坐標：N個原子的體系的結構至少由3N-6個坐標才能描述，也即內坐標。如果在內坐標基礎上，額外添加一些幾何變量，使總坐標數多于3N-6個，就叫冗余內坐標。這些多出來的變量雖然對于描述體系結構是多余（冗余）的，但是它們的存在對于幫助幾何優化收斂往往是很有意義的。在默認情況下，Gaussian程序會根據輸入的幾何坐標自動構建出冗余內坐標，構建的方法參見J.Comp.Chem,17,49。這樣自動構建的冗余內坐標中只要是離得近的原子之間都有鍵長項（并可能伴隨產生角度、二面角項）。對于單分子體系冗余內坐標下優化通常總比在內坐標下和笛卡爾坐標下優化收斂得要容易。特別是對于環狀體系，在這樣的冗余內坐標下優化明顯比在內坐標下好得多，這是由于環中首尾原子離得很近，因此在冗余內坐標下它們之間會有對應的鍵長、角度、二面角項，故而環上的每個原子在優化中都能夠同等地對待。
對于氫鍵，在自動設定的冗余內坐標中一般也會被認為是成鍵的。氫鍵長度在優化過程中有專們對應的鍵長項來描述，這對收斂是有益的。然而，弱一些的氫鍵（或者鹵鍵之類偏弱的相互作用），以及范德華作用的原子間由于相距較遠而超過了閾值，往往不會自動添加對應的鍵長項。為了幫助收斂或試圖解決不收斂，這時應當自行添加冗余內坐標讓它們之間存在鍵長項。在Gaussian中，可以用opt=modredundant，然后在分子坐標末尾空一行寫上要添加鍵長項的原子序號即可，例如
3 7
12 14
8 22
這樣，在優化過程中輸出的冗余內坐標的信息中，就會看到新增了R(3,7)、R(12,14)、R(8,22)鍵長項了。
如果是很多分子的團簇的優化，例如水簇的優化，若發生震蕩難以收斂，可以對所有H---O氫鍵都加上對應的鍵長項。據說給O---O也加上鍵長項的話還能令收斂變得更容易一些。
盡管冗余內坐標比起內坐標增加了一些變量，因此計算梯度、Hessian時原則上會略慢于用內坐標的情況，但是相對于減少了的優化步數，這點代價是不值得一提的。
5 調整步長上限
opt=maxstep=N（等價于IOp(1/8=N)）可以將步長上限（也稱置信半徑）調整到0.01N Bohr或弧度，默認是30。當按照優化方法判斷出的下一步幾何變量的變化超過這個數值，則實際的改變量會調節到恰好等于這個上限值（具體來說，是在置信半徑對應的球面上做最小化尋找最佳的位置）。如果優化發生震蕩，且每步位移不是很大而受力較大的話，往往表明優化過程陷在比較深、比較窄的勢阱里，此時可以降低步長上限，比如降低到N=5乃至N=3，以避免每次都走過頭而始終達不到勢能面極小點。減小步長的做法不要在優化過程中過早地使用，否則由于優化過程走得太慢，所需的優化步數會增加。另外，如果用了calcall，減小步長上限不一定有好處，因為在精確Hessian下進行優化，每一步走的是比較準確的，人為篡改了步長有可能還阻礙收斂。
注意，maxstep設的只是最初的步長上限。隨著優化的進行，步長上限其實是動態調整的。如果想一直用自己的設定，需要在opt中加上NoTrust關鍵詞。對于結構反復微小震蕩的情況，在使用maxstep減小步長上限的同時總建議同時寫上此關鍵詞。
6 調整對稱性
如果優化一開始的初猜結構具有較高對稱性，那么優化過程中這個對稱性多數情況下會一直保持。如果實際上體系本身勢能面極小點的結構并沒有那么高的對稱性，則很容易造成優化不能收斂，而且即便收斂了也會有虛頻。此時應當嘗試人為地隨意擾動一下初始結構，比如讓某個二面角稍微扭轉一點來破壞對稱性。
反之，如果一個體系看起來有對稱性，而且優化過程中也確實明顯看到體系逐漸變得對稱性越來越高，但是始終沒收斂，則建議將結構取出，在Gview里打開Edit-Point Group工具，將tolerance設低一些，使Gview能判斷出此分子應有的最高階點群，然后點Symmetrize，即調整結構使得體系的結構嚴格滿足這樣的點群。接下來再用Gaussian優化時，就會在較高階點群下優化，這大大節省了計算能量及其導數的時間，結構還會更精確（無對稱性下優化的結構雖然看上去有了對稱性，但由于數值上的微小偏差其坐標并沒有完全精確滿足對稱性），也由于要優化的自由度顯著變少了，使得優化過程容易完成。有時候在Gview里已經按照指定的點群對初始結構做了對稱化，但是Gaussian在算的時候卻認不出相應的點群，尤其是對于高階點群（比如Gaussian很難認出C60的Ih點群），這是由于程序自身毛病以及輸入文件中記錄精度的問題導致的。可嘗試在輸入文件中使用內坐標而非笛卡爾坐標記錄結構，也可以用symm=loose來放寬判斷點群的標準。
7 增加DFT泛函積分精度
DFT計算中交換相關泛函由于形式復雜，而且種類奇多，所以它涉及到的項不是利用解析積分方法精確計算而是利用數值格點積分方法來近似計算的，詳見《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）。積分格點數目越多，積分就越準確，DFT計算的能量、受力、Hessian矩陣就算得越準確。如果積分精度不夠高，不僅可能造成優化收斂困難，最終優化的結果也可能不很精確。Gaussian默認的泛函積分精度一般來說也夠了，但如果碰上不收斂（特別是對于明尼蘇達系列泛函，如M06-2X），可以試試提高積分格點數，比如使用int=ultrafine。注：從G16開始，默認就是int=ultrafine了。
8 改變收斂標準
改變收斂標準是沒辦法的辦法。如果當前研究對結構精度要求不需要那么高，比如只要有個大概結構就可以，或者是作為其它什么任務（比如動力學）的初始結構，若發現結構反復微小地震蕩，而且離收斂限也不遠了，索性就直接停掉優化，當做已經收斂就行了。在Gaussian中也可以用opt=loose放寬收斂限。反之，如果對幾何優化精度要求高，比如一些存在很弱相互作用的體系的研究（如H2二聚體），或是需要計算準確的振動頻率，那么應該將收斂限設嚴，用opt=tight乃至verytight。
9 嘗試不同的初始結構
比如，優化過程中某個二面角反復在兩個角度來回晃蕩就是不收斂，那么可以試試將初始結構中這個二面角設在這個兩個角度的平均值，以此作為初始結構（最好同時結合calcfc）。此做法解決問題的幾率不大，但是偶爾碰巧能解決，畢竟本身幾何優化這種迭代的數值過程能否收斂就是有很多巧合性。如果優化過程的結構變化已經向著不靠譜或者非預期的方向上走了，往往是因為初始結構明顯不合理所致，需要重新調整結構重新優化再試。
10 嘗試不同理論方法和基組
某個計算級別下如果優化很難收斂，可以嘗試換個相近的基組，或換個相近的理論方法（比如B3LYP換成B3PW91之類），如果優化能收斂，可將此結構作為原先級別下優化的初始結構（也可以將最后的波函數作為初猜，Hessian矩陣也一起讀入）。理論方法和基組不應和當前計算級別相差太多，否則勢能面極小點位置可能相差得不少，對解決收斂問題將沒多大用處。這種做法屬于撞大運，起效幾率不大，但實在沒轍了可以試試。或者，你的文章中所用的理論方法干脆就換成一個精度相仿佛，對當前體系又發現相對更容易收斂的方法。
如果主要目的是想節省整個優化過程的計算耗時，那么建議先在級別低一些的理論方法和基組下（或者用比較好的半經驗方法，或適宜的分子力場）在loose收斂限下做預優化。不過這對解決原先級別下不收斂的問題倒未必有多大幫助，因為低級別下與原先級別下勢能面的極小點結構往往有不可忽視的偏離。另外，較低的計算級別選擇得必須有意義，比如想在高精度下優化一個范德華復合物，事先用諸如B3LYP這種根本沒法合理描述色散作用的泛函預優化不僅白費時間，預優化出的結構還很可能比你最初搭的結構偏離極小點結構更多。
11 關于溶劑模型
溶劑環境下的幾何優化有沒有必要帶溶劑模型在《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）里我已經專門說過了，不再累述。這里特別強調的一點是，對于opt freq任務，不建議用常用的SMD溶劑模型，因為我和很多人都發現用SMD的時候幾何優化可能變得更困難，而且收斂后做振動分析時還容易出現虛頻，這大抵是SMD溶劑模型使得勢能面出現一定數值噪音所致。用Gaussian的scrf關鍵詞默認的IEFPCM溶劑模型則沒有顯著這種問題。但使用IEFPCM也不收斂的話，而且若當前體系的優化沒有絕對必要帶溶劑模型，也可以嘗試在真空下優化再試。
幾何優化和優化后計算單點兩個過程對溶劑效應的考慮、用的溶劑模型都沒有必要非得統一。在優化之后算單點能時可以放心用SMD。
具體怎么解決優化不收斂得看實際情況，特別是要用gview檢查優化過程的結構變化趨勢和能量、受力曲線，看是否最后開始震蕩了（如果能量/受力還有降低的大趨勢則應當繼續跑下去）。雖然上面討論了很多，但從實用性出發，對于大多數情況，對于初學者筆者推薦按照以下順序去嘗試解決震蕩問題：
1 對于DFT，特別是明尼蘇達系列泛函（如M062X、M06、M06L），首先加上int=ultrafine再試（從G16開始此為默認，不再需要嘗試此項）。如果用了SMD，改用IEFPCM
2 如果體系小，嘗試opt=calcall。如果體系大算得慢，此關鍵詞會令計算量進一步激增，實在沒轍再考慮這個，或者嘗試opt=recalc=3或5
3 嘗試opt=gdiis
4 嘗試opt(gdiis,maxstep=x,notrust)，x為3~5
5 如果上述做法都不靈，干脆改成loose收斂限，但如果之后做振動分析容易有虛頻
尋找過渡態，或者說優化過渡態，本身也是幾何優化過程，只不過是向勢能面一階鞍點優化，而不是向極小點優化。在Gaussian里一般用opt=TS方式基于一個初猜結構來尋找過渡態，通常用這樣的關鍵詞：opt(calcfc,noeigen,TS)。尋找過渡態也經常會遇到幾何優化不收斂問題。由于在主流量化程序里優化過渡態和優化極小點的算法基本沒有區別，因此上面所有解決優化極小點不收斂的方法，對于優化過渡態也完全適用。但優化過渡態對初始結構的敏感性遠遠高于優化極小點，因此優化過渡態過程幾何優化不收斂有很大可能是因為初猜的過渡態結構離實際過渡態結構不夠接近，導致最終被優化到沒有化學意義的結構所造成的。到底是由于初猜結構不好，還是由于最終發生了震蕩導致優化過渡態沒有收斂，通過gview查看優化過程的結構變化立刻就能知道。對于不是很簡單的化學反應，研究者一次就能給出能最終優化到過渡態的初猜結構并不容易，故失敗時不要輕易氣餒，要耐心反復嘗試。Gaussian還支持QST2方法優化過渡態，它是根據用戶提供的反應物和產物結構通過插值產生過渡態初猜結構。但這個初猜結構往往很不合理，明顯不如有化學直覺的用戶自己手動擺出來的，因此一般情況不建議用QST2。而當你用QST2時發現結構向著不靠譜的方向優化，那應當立馬放棄QST2而改用opt=TS。另外，優化過渡態對于Hessian矩陣的質量比起優化極小點敏感得多得多，因此解決優化過渡態不收斂的問題，在計算條件允許的情況下，建議優先考慮calcall或recalc=3（對于小體系，由于計算Hessian矩陣耗時不算特別高，對于不太好找過渡態的情況，為了減少反復嘗試的次數，建議總是用calcall或recalc）。