解決SCF不收斂問題的方法

筆者注：讀者不要說諸如“此文我看了，還是不收斂”、“我用了文中的方法還是不收斂”這類話，然后問我怎么解決，筆者極度反感這種說法。絕對沒有哪個方法能100%保證解決SCF不收斂。本文絕對沒有一條餿主意，所有Gaussian支持的真正有可能能解決SCF不收斂的做法在本文都已經全面列舉了，沒有任何遺漏，沒有可補充的。倘若本文的做法都做了嘗試還沒解決，那也別指望有任何其它方法能解決。才隨便試了文中一個做法，或者以不正當的方式（沒有真正了解原理的基礎上）瞎試幾個就說沒解決問題，這根本什么也說明不了。要想解決問題，必須認真閱讀本文，所有有可能對當前問題有用的方法都依次嘗試，有必要的時候幾個一起結合使用。本文列舉的做法已經把最壞的情況都考慮了，所以出現所謂的“看了本文還是沒解決問題”的情況一定是讀者還沒真正仔細看、認真領會、充分嘗試。

1 前言

量子化學計算離不開SCF（自洽場）迭代，如半經驗方法、HF、DFT等。在SCF迭代中，由Fock矩陣F對角化獲得新的系數矩陣C和軌道能{ε}，然后構造密度矩陣D=C'C'^(T)，其中C'為不含虛軌道的C矩陣，再由D構造新的Fock矩陣，反復進行直到收斂，可以寫為F_(1)->C_(1)->D_(1)->F_(2)->C_(2)->D_(2)...。收斂判據不是唯一的，比如Gaussian中用的判據是當前步與上一步的密度矩陣元變化量的最大值和方均根(RMS)以及能量的變化，當數值都小于一定范圍就認為已經收斂了。默認判據和設定下多數情況在25步以內能達到收斂。但慢收斂甚至完全不能收斂的情況是經常會遇到的。常見的是迭代后期能量隨迭代呈現震蕩，直到達到默認的最大步數仍未收斂。也可能能量雖然震蕩但總趨勢是慢慢降低的，一直迭代下去能收斂，但震蕩行為明顯拖慢了收斂速度。也可能震蕩的規律性不顯著，迭代過程能量曲線看上去有隨機性，但就是很難達到收斂限。不收斂、難收斂情況的出現有很多數值巧合因素，但有些情況不好收斂是眾所周知的，例如：基組含彌散函數、體系處于明顯的非平衡構型、HOMO-LUMO能隙較小（過渡金屬化合物、成鍵方式古怪等靜態相關較強體系中容易出現此情況）、限制性開殼層(RO)計算、明尼蘇達系列泛函。

人們提出過一些方案加速收斂或試圖解決不收斂，如DIIS、阻尼方法、溫度展寬、能級移動、二次收斂方法等。主流的量化程序，如Gaussian、ORCA中默認就會使用一部分這樣的幫助收斂。而這些主流的量化程序為了加速SCF計算耗時，會引入一些數值近似，比如一開始先默認用較低的DFT積分格點、忽略數值較小的積分、Incremental fock方式加速Fock矩陣的構建等等，這些近似有時候會阻礙收斂。

Gaussian中對應SCF不收斂的錯誤提示是L502出錯，并伴隨Convergence criterion not met和Convergence failure -- run terminated.的具體提示（但從G16 B.01開始，僅提示Convergence criterion not met而不報錯中斷）。SCF不收斂怎么辦是初學者最最最常問的問題，本文第2節將解決不收斂問題的可行辦法進行了匯總，并給出在Gaussian中對應的關鍵詞。缺乏常識的讀者千萬不要輕易聽信計算化學公社論壇（http://bbs.keinsci.com）以外的網上其它一些地方（尤其是中文的）的關于解決SCF不收斂的說法，那里面的做法很多都是嚴重誤人子弟，白耽誤工夫還不有效解決問題。如果讀者對其中的DIIS、能級移動、阻尼、溫度展寬和二次收斂方法的原理感興趣，可以閱讀本文第3節的詳細介紹。

順帶一提，對于Gaussian用戶，如果憑經驗預感當前體系SCF不容易收斂，或者算的體系比較大，或者遇到SCF不收斂后加了本文提到的幫助收斂的關鍵詞重算，一定要用#P，這樣SCF迭代過程的每一輪信息才會輸出出來，才能夠了解當前已經迭代到了第幾輪，能量變化（Delta-E）、密度矩陣元最大變化（MaxDP）和密度矩陣平均變化（RMSDP）是多少，當前離收斂限還有多遠（收斂限會在SCF開始之前明確提示）以及迭代過程是否有收斂的趨勢。如果沒用#P，顯然用戶會被完全蒙在鼓里。

關于解決CP2K程序中SCF不收斂問題，筆者另有一篇文章：《解決CP2K中SCF不收斂的方法》（http://www.shanxitv.org/665），其中有些解決思想和本文是共通的。

2 Gaussian中解決SCF不收斂問題的常用做法

絕對沒有哪個或哪些關鍵詞只要加上就一定能解不收斂問題，解決不收斂問題需要在理解關鍵詞含義的基礎上反復嘗試。以下方法越靠前的往往越值得優先考慮，大多數可以結合在一起使用以得到更好效果。

• (1) 對于M06、M06-2X等明尼蘇達系列泛函可嘗試增加泛函積分格點的精度，對于其它泛函這個做法見效幾率則比較有限。G09默認是int=fine，相當于(75,302)，可以提升到int=ultrafine，相當于(99,590)，耗時也會增加，收斂的幾率也會增加，結果精度也會提升。
  注：從G16開始，默認就是int=ultrafine了。
• (2) Gaussian會自動在計算初期時降低積分計算精度以加快計算，但有可能因此阻礙收斂。對于使用了彌散函數的情況（其它情況不算）出現不收斂時建議先嘗試SCF=NoVarAcc來避免因自動降低積分精度而導致的這個問題。
• (3) 用int=acc2e=12可以增加積分精度，對因為使用大量彌散函數（其它情況不算）導致的大體系難收斂問題奏效幾率較大。G09默認的是int=acc2e=10。
  注：從G16開始，默認就是int=acc2e=12了。
• (4) Gaussian默認使用Incremental Fock方式以近似方式構建Fock矩陣來顯著節約迭代過程每一步的時間，但可能因此阻礙收斂。可以用SCF=noincfock避免這個做法。
  Hint：含有很大量彌散函數的基組（如aug-cc-pVTZ及更大）計算中、大體系時，往往很難收斂，建議SCF(novaracc,noincfock) int=acc2e=12三管齊下，解決問題的幾率較大，而耗時增加不算很大。但這樣并非總是奏效，把彌散去掉后將收斂的波函數作為初猜也是很有效的，見下文的(8)。如果你的基組里彌散函數不多，或者根本都沒加彌散函數，靠SCF(novaracc,noincfock) int=acc2e=12解決不收斂的概率不太高，不奏效時別忘了嘗試本文其它方法。
• (5) 用能級移動法提升虛軌道能量以增大HOMO-LUMO gap，避免虛軌道和占據軌道間過度混合：SCF=vshift=x。x一般為300~500。這只影響收斂過程，而絕對不影響任何最終計算結果，包括軌道能級。對于HOMO-LUMO gap較小情況，常見的比如含過渡金屬的體系，值得嘗試此做法。
• (6) 用SCF=conver=N關鍵詞改變收斂限，這代表令密度矩陣RMS收斂限為1E-N，密度矩陣最大變化和能量收斂限都為1E(-N+2)。G09/16對單點計算默認是SCF=tight，相當于SCF=conver=8，然而此時對密度矩陣的收斂要求有點太嚴了，往往不容易達到，而且達到的時候能量通常都早已經收斂到極高精度了。對于計算單點能的目的可放心降到conver=6，相當于把收斂標準放寬100倍，通常收斂時能量變化已經非常小了。但是做幾何優化、振動分析的時候不建議降低默認的SCF收斂限，否則可能結果不準確，還可能阻礙幾何優化的收斂。
  不過，對于優化來說，有些體系的初猜結構可能構建得和實際結構偏差較大，導致優化初期SCF收斂難，此時可降低SCF收斂限到=6（但不要設得更松）使結構粗略收斂，再在最后的結構下用默認的SCF收斂限進一步嚴格優化，是很有幫助的做法。
  注1：大多時候Gaussian會在計算初期默認用相對較低的精度算積分，到后期才還原為全精度，相當于啟用varacc，這時候SCF=conver=N關鍵詞并不生效（我認為這是bug）。用novaracc關閉此做法時才能確保SCF=conver=N肯定能實際生效。
  注2：如果你用G16做幾何優化、振動分析等任務，或使用TDDFT、雙雜化、后HF等方法，放寬SCF收斂限的時候必須加上特殊的IOp，否則會因為通不過收斂精度檢測而遇到報錯，參見《令Gaussian 16中SCF未收斂到默認收斂限也能繼續做后續計算的方法》（http://www.shanxitv.org/625）。
  
• (7) 若某個泛函下計算不收斂，可以嘗試其它泛函，如果能收斂，再用guess=read讀取其收斂的波函數作為初猜。比如M06-2X等明尼蘇達系列泛函SCF難度往往高于其它泛函，碰到難收斂的時候，可以試試B3LYP等泛函，如果發現能收斂，將之作為M06-2X的初猜，解決問題的幾率不小。有一個值得一提的經驗是通常HF成份越高的泛函往往越容易收斂，特別是對于含過渡金屬的體系來說。
• (8) 小基組通常比大基組更容易收斂。直接用大基組難收斂但小基組能收斂時可將小基組收斂的波函數作為大基組計算時的初猜。比如原本打算用的def2-TZVP收斂不了，但發現def2-SVP能收斂，于是可嘗試用后者收斂的波函數作為def2-TZVP計算時的初猜。再比如，眾所周知加彌散函數后收斂會變得更困難，因此比如用aug-cc-pVTZ計算時發現不收斂，應立刻想到試試cc-pVTZ能否收斂，若能則將其收斂的波函數當aug-cc-pVTZ的初猜。但兩個基組尺寸差異太大的話這么做也沒什么意義，比如拿STO-3G收斂的波函數當def2-TZVP的初猜對收斂并不會帶來明顯的益處。對于特別難收斂的情況，不嫌麻煩的話可以嘗試逐級提升基組檔次的做法，步步為營，如STO-3G->3-21G->6-31G*->6-311G**->def2-TZVP。
• (9) 改變默認的初猜方法，可使用比如guess=huckel（指的是擴展Huckel方法）、INDO（實指INDO、CNDO、Huckel的混合）等關鍵詞。默認的是使用Harris泛函的DFT方法的結果作為初猜，是通過自由原子密度疊加作為分子電子密度來構建KS算符，然后將變分得到的軌道用做后續計算的初猜，這步只做一次而不像普通DFT方法繼續做迭代。core的初猜是過于粗糙的，說明初始的Fock矩陣元只含單電子項而不含雙電子項，即密度矩陣用空矩陣，這默認被用作一些半經驗方法計算的初猜。擴展Huckel計算也不需要迭代，也默認被用作一些半經驗方法的初猜。INDO由于是迭代的方法，給高等級計算提供初猜波函數之前自身也需要初猜（擴展Huckel方法），本身亦有不收斂的可能。guess=AM1一般不能直接使用。
  
• (10) 使用二次收斂方法（Quadratic convergence, QC）。這種方法收斂所需步數通常比默認情況更少，有一定可能性解決SCF不收斂。SCF=QC代表使用這種方法。G09中，SCF=XQC代表先用普通方式迭代到64輪（或者自己用scf=maxcyc設的值），如果不收斂再自動切換為QC。而G16中用SCF=XQC時，是用MaxConventionalCycles=N設定常規迭代多少輪不收斂才切換為XQC，N默認為32。
      筆者強烈不建議輕易嘗試QC來試圖解決SCF不收斂，更絕對不要把scf=qc或scf=xqc當做默認用的關鍵詞！！！（然而老有初學者總想著用QC，自取其辱！）用此做法前一定要充分嘗試本文說的一大堆其它方法，實在絕望時再試圖用QC。QC不能輕易用的原因有三：(1)QC方法耗時頗高，SCF每一輪的耗時是平常的好幾倍 (2)如果必須借助此方法SCF才能收斂的話，非常有可能最終并沒有收斂到最穩定的波函數，導致得到的能量高于此結構下當前級別的基態的真實值，因此結果無意義。用QC時建議同時加上stable關鍵詞檢測一下所得波函數穩定性 (3)用QC時有很高概率照樣最后還是不收斂，最終會提示L508報錯，這期間會做漫長的計算，導致花了大量時間做無用功。順帶一提，網上有人在低水平場所說QC方法是什么“強制收斂”，這完全不懂裝懂的人胡說八道以訛傳訛！QC哪可能保證最后一定能收斂？
• (11) 使用Fermi展寬：SCF=Fermi。此時默認也會擇情自動使用能級移動和阻尼。
• (12) DIIS是默認的加速SCF收斂的方法，使得收斂所需的步數普遍減少許多。但極個別情況也可能反倒導致不收斂，可嘗試關閉之：SCF=noDIIS。
• (13) 略微改變分子幾何結構，比如稍微縮短/伸長鍵長、改變鍵角，若能得到收斂的波函數，可作為原來的幾何構型下計算的初猜。Gaussian在幾何優化中一般會自動將當前步結構的收斂波函數作為計算下一步結構的波函數的初猜，與這個策略有相似之處。
• (14) 對于開殼層體系，可以先計算相應的閉殼層離子體系得到收斂的波函數，然后讀取其收斂的波函數作為初猜。注意，RO方式計算開殼層體系比非限制性計算(U)難收斂得多得多，如果RO收斂不了，先改成U再說！
• (15) 電子數少的時候往往更容易收斂，比如陽離子比陰離子一般容易收斂。所以可以先計算電離態，將得到收斂波函數用作初猜。
• (16) 對于溶劑下的計算若SCF不收斂，也可以試試其它溶劑下或者真空下或者其它溶劑模型下能否收斂，如果可以，將收斂的波函數當做原溶劑下計算的初猜。
  
• (17) 用SCF=maxcyc=N（等價于scfcyc=N）加大迭代次數上限到幾百、上千。這是解決SCF不收斂問題的最蠢的方法，但由于以訛傳訛，卻成了初學者們最愛用的做法。如果一個搞量化的碰見SCF不收斂時第一反應是使用這個辦法解決，或者看到某人的輸入文件里總是帶著比如SCF=maxcyc=500這樣的關鍵詞，那么他一定是個菜鳥！初學者總喜歡將N設得很大，這種做法>99.9%的情況下都絕對沒用。在Gaussian默認的128次迭代次數內若不能收斂，N加大到多大幾乎都是白搭，還白白浪費計算時間。僅當你觀察迭代過程能量變化，發現隨著迭代能量有緩慢降低，超過128步后繼續迭代下去有希望收斂，才有必要加大循環次數上限，但極少會有這種情況出現。其實稍微有點基本的邏輯就能明白，倘若迭代次數上限設到幾百上千就能解決不收斂，那量化程序干嘛不默認就設到幾百上千？
• (18) 上述做法若都不奏效，索性改用其它方法或基組，或嘗試其它程序來計算當前體系。如果發現其它程序用相同的級別計算后SCF能收斂，而你為了保持結果嚴格可比性又不想換程序，那么可以用Multiwfn載入那個程序產生的含有基函數信息的文件（如.molden，詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》http://www.shanxitv.org/379），然后進入主功能100的子功能2，選擇導出波函數為fch文件，之后再用Gaussian的unfchk工具將之轉換成chk文件，然后Gaussian計算時用guess=read從中讀取波函數作為初猜，此時有很大概率可以收斂（Gaussian收斂的波函數也可以反過來給其它程序用，見比如《將Gaussian等程序收斂的波函數作為ORCA的初猜波函數的方法》http://www.shanxitv.org/517）。

用IOp(5/13=1)是最腦殘的“解決”不收斂的做法，文末的附錄專門對此進行了大批判。

值得注意的是，化學意義越強的結構、電子態通常越容易收斂。諸如成鍵方式普通的有機分子在平衡結構下計算基態時SCF極少會碰到不收斂的情況。遇到非常難收斂的情況應注意是否有以下問題：

需要注意的是，化學意義越強的結構、電子態通常越容易收斂。諸如成鍵方式普通的有機分子在平衡結構下計算基態時SCF極少會碰到不收斂的情況。遇到非常難收斂的情況一定要注意是否有以下問題：
1. 某些原子間距離太近
2. 用的結構相對于真實的結構有嚴重變形，此時電子結構很詭異
3. 搭建的結構異想天開（對體系特征缺乏基本的認識、缺乏結構化學基本常識），或者存在嚴重不合理性（如該有氫的地方少了氫，用簇模型計算時沒有對被截斷的共價鍵進行飽和）
4. 對于過渡金屬配合物，自旋多重度設定得不符合實際基態（非基態的情況往往更難收斂）
5. 關鍵詞寫錯了。諸如Gaussian里用贗勢卻寫gen而非genecp，導致設的贗勢根本沒讀入
6. 一些設定明顯不妥，比如凈電荷設的根本不對。或者有些設置就是會導致難收斂，比如考察外電場下的情況，電場越大通常SCF越難收斂，過大的電場甚至會導致電子脫離體系，根本就不可能收斂，參見《一篇文章深入揭示外電場對18碳環的超強調控作用》（http://www.shanxitv.org/570）里的討論

使用Stuttgart小核贗勢的時候，含有稀土元素的體系SCF收斂難度比起其它體系往往大得多，這主要是由于它的f電子所致，這時候可以考慮使用Stuttgart大核贗勢，這種贗勢對稀土元素通常有針對不同氧化態的版本，對應于不同f電子數，應當根據實際體系恰當選取，否則結果沒意義。如果算能量時對能量精度要求較高，或者要對稀土元素和周圍原子之間成鍵情況做波函數分析，則一般應使用小核Stuttgart贗勢，此時碰到SCF不收斂只能硬著頭皮試圖通過上述方式解決。相關信息見《在贗勢下做波函數分析的一些說明》（http://www.shanxitv.org/156）、《談談贗勢基組的選用》（http://www.shanxitv.org/373）、《詳解Gaussian中混合基組、自定義基組和贗勢基組的輸入》（http://www.shanxitv.org/60）。

另外，MCSCF（最常用的形式是CASSCF）也是需要做SCF迭代的方法，其SCF收斂的難度明顯高于HF/DFT，其不收斂的解決和上述方法既有一定共性也有明顯差異（比如vshift、QC等都沒法用）。在Gaussian中，碰到其不收斂的情況有以下辦法可以嘗試：
(1)檢查活性空間的設置是否合理，是否已經正確調換了軌道順序
(2)嘗試其它類型的初猜軌道
(3)嘗試不同基組。將小基組收斂的波函數做為中/大基組的初猜的這種做法對CASSCF尤為重要
(4)用SCF=conver放寬收斂限
(5)用SCF=maxcyc增大迭代次數上限
Gaussian在MCSCF這方面的收斂性做得明顯不及Molpro、Molcas等專長于MCSCF、多參考方法的程序。實在沒轍的話建議嘗試這倆程序，或者嘗試這方面相對遜色一些的ORCA也可以。

3 促進SCF收斂的一些方法的原理

3.1 阻尼方法(Damping)

設正常解得的第n步的密度矩陣為D_(n)，阻尼方法使實際用于構建第n+1步Fock矩陣F_(n+1)的密度矩陣變為D'_(n)，D'_(n)=w*D_(n-1)+(1-w)*D_(n)。w是權重系數，既可以設為常數也可以根據迭代過程動態調整。在迭代出現震蕩時，D(n),D(n+1),D(n+2)...的變化呈鋸齒狀，用w參數平均化后的密度矩陣D'_(n)代替D_(n)就削弱了當前步與上一步密度矩陣之間的差異，使密度矩陣隨迭代變化更為平滑，幫助收斂。Gaussian在迭代初期會使用動態阻尼，這個方法對收斂很有幫助，但也并非總能奏效，比如前線軌道能級密集的情況。

3.2 迭代子空間中直接求逆(DIIS,Direct Inversion in the Iterative Subspace)

這是由Pulay發展的基于外推的方法。常規的SCF迭代收斂緩慢，DIIS能明顯減少收斂所需迭代次數，加快SCF計算，是SCF計算中使用得十分普遍的方法，使多數分子都能在20步以內收斂。這個方法利用之前步的信息來估算出最好（最接近收斂，亦即“誤差”最小）的Fock矩陣。下一步的Fock矩陣F_(n+1)由之前步的Fock矩陣線性組合而成，F_(n+1)=∑[i]c(i)*F_(i)，c(i)是組合系數，i的加和從1到n，也可以設成比如從n-10到n。引入的每一步的F(i)都會引入誤差，以誤差矩陣Err_(i)表示，總誤差矩陣Err_tot=∑[i]c(i)Err_(i)，誤差函數ErrF(c)為Err_tot的模。若找到一套系數能讓誤差函數ErrF(c)最小，則組合出的F_(n+1)就是最佳、最接近收斂的Fock矩陣。只需通過令ErrF對每個c(i)求導得0，并且用拉格朗日乘子法將歸一化條件∑[i]c(i)=1作為限制，就能求解出各個組合系數。這等價于求解矩陣方程Ac=b，因此c=A^(-1)b，其中對矩陣A直接求逆的操作就是名字中Direct inversion的由來，iterative subspace就是指由之前迭代步的信息構成的子空間。若只利用當前步與此前的10步的信息構建F_(n+1)，則子空間的維度就是10，之前陳舊的信息就被丟掉了，這也是常見做法，因為那些信息對當前迭代的貢獻已經很小了，留著還會占內存。

有了上面的思路后還需要具體定義誤差矩陣。第n步迭代求解HF方程就是求解F_(n)C_(n)=SC_(n)E_(n)，E是對角矩陣，矩陣元為軌道能。顯然用上一步的系數矩陣不可能令此式成立，否則當C_(n-1)=C_(n)時就說明已收斂了，即F_(n)C_(n-1)-SC_(n-1)E_(n)≠0，此式左邊的結果正是誤差矩陣，若將式子右乘C_(n-1)^(T)S，則可整理為Err_(n-1)=F_(n)D_(n-1)S-SD_(n-1)F_(n)，這就是Pulay定義的誤差矩陣。

也有使用以能量描述誤差的相似方法，稱為EDIIS，E代表Energy，它是DIIS與后文中RCA的思想結合的方法。EDIIS中下一步密度矩陣D_(n+1)=∑[i]c(i)D_(i)，由于總能量是密度矩陣的函數，求出一套系數令總能量最小，就確定了D_(n+1)。注意EDIIS求解c(i)的時候限制它們必須都大于或等于0，而DIIS沒這個要求，所以EDIIS被稱為內插方法，而DIIS被稱為外推方法，但實際上DIIS迭代中既有內插也有外推的情況。EDIIS初猜的密度矩陣比DIIS更為隨意，都能很快收斂，但初猜的密度矩陣離收斂比較近時DIIS收斂更快。將EDIIS與DIIS組合使用，誤差矩陣最大值較大時使用EDIIS，較小時使用DIIS，可以成為兼具二者優點的既高效又穩健的方法，這也是Gaussian默認的。

3.3 能級提升法(Level Shifting)

也稱能級移動法。先介紹下SCF迭代過程的本質。設第F_(n)對角化得到的MO軌道波函數系為ψ_(n)（其中包括ψ_(n,1)、ψ_(n,2)...ψ_(n,N)，N為基函數數目，能量由低到高排序），其中占據的軌道部分叫ψocc_(n)，非占據的部分叫ψvir_(n)，下一步n+1的Fock矩陣元F_(n+1)依靠ψocc_(n)的信息構建。以ψ_(n)為基時，顯然F_(n)是對角矩陣，矩陣元就是各個軌道的能量。在收斂前以ψ_(n)為基時F_(n+1)不是對角矩陣，如果仍是對角矩陣，等于不用再做對角化，ψ_(n+1)就等于ψ_(n)，就說明已經收斂了。對F_(n+1)進行對角化所得到的變換矩陣X，就是描述在這一次迭代中新得到的ψ_(n+1)是如何由上一步的ψ_(n)轉換而來。由于ψ_(n)和ψ_(n+1)都是正交基，所以X為酉矩陣。X可看成由四個子方陣組成，分別是對角塊X[occ(n)_occ(n+1)]和X[vir(n)_vir(n+1)]，以及非對角塊X[occ(n)_vir(n+1)]和X[vir(n)_occ(n+1)]（二者互為轉置關系）。由于占據軌道波函數之間的酉變換不改變體系能量，而虛軌道對能量有沒有貢獻，所以X[occ(n)_occ(n+1)]和X[vir(n)_vir(n+1)]的數值對迭代能量變化并沒有貢獻，只有對應占據軌道與虛軌道混合的非對角塊X[vir(n)_occ(n+1)]才會使總能量變化。

當前步解出來的虛軌道與占據軌道能量差值越小，說明它們越相似，HF方程的近似解也越容易是由它們為基組合而成，所以下一步所得的占據軌道就會混進很多當前步的虛軌道成分，使總能量變化劇烈，當然總能量下降得快是好事，但也因此可能出現總能量不降反升的情況。若每步虛軌道與占據軌道能量差都不大，則每步它們都有大規模的混合，可能能量降低與升高會混雜出現而產生振蕩的情況，必然使收斂緩慢。Saunders等人提出的能級提升法就是人為地將每步得到的虛軌道能量提升（對應的虛軌道波函數也因此改變），這樣每步占據軌道與虛軌道的混合就會減小，如果提升得足夠多，可以證明一定能在迭代中令總能量不斷降低，而不會出現升高的情況，如此可以使能量最終收斂于極小點。雖然虛軌道能量提升越多，收斂越穩定，總能量變化越緩慢，越能保證總能量每步迭代都能降低，但如果原始的迭代過程本來就能讓每步的總能量都降低，提升虛軌道能量則會使降低的量也減小，導致收斂速度減慢，還有可能出現低于占據軌道能量的空軌道，形成所謂的“空穴”，這樣的非基態結果沒什么意義，所以能量也不能肆意提升。能量提升值的選取有一定任意性，對不同虛軌道可以統一提升同一個值也可以提升不同值。

將第n步虛軌道的能量值提升x也就是使得E_(n)對應虛軌道能量的矩陣元都增加了x成為E'_(n)，因此需要修改F_(n)成為F'_(n)，這樣通過求解F'_(n)C'_(n)=SC'_(n)E'_(n)就能得到所期望的E'_(n)以及相應的軌道波函數C'_(n)。用于限制性閉殼層SCF，可以推出F'_(n)=F_(n)+x*(S-0.5*SD_(n)S)。

在Gaussian中使用SCF=Vshift=n關鍵詞可以設定將虛軌道能量提升n*0.001hartree，一般n設為幾百，若仍不收斂可嘗試提高更多，這個方法對解決不收斂問題很有效，可以嘗試多提升一些。如果體系本身就容易收斂則不要用這個方法，會使收斂更慢。注意能級移動只是幫助收斂的方法，在最終輸出的能量中會從所得虛軌道能量中減去這個值，修正回實際的虛軌道能量，所以對結果沒影響。

PS：這個問題可以通過二階微擾理論做更細致的分析。將f_(n+1)與f_(n)的差作為微擾項f~，對第i個占據軌道的一級校正波函數δψ(i)=∑[j∈vir]H(i,j)/(e(i)-e(j))*ψ_(n,j)，其中H(i,j)=<ψ_(n,i)|f~|ψ_(n,j)>，j為虛軌道，e是相應軌道的能量。所以ψ_(n+1,i)≈ψ_(n,i)+δψ(i)，因此軌道i能量變化為δe(i)=∑[j∈vir]H(i,j)^2/(e(i)-e(j))，體系總能量的變化δE=∑[i∈occ]δe(i)。可見軌道能級差e(i)-e(j)越小，從虛軌道j混進占據軌道i的量越多，δE的值也會越負。

當F~較大時，由于高階項不可忽略，加上后可能實際的δE是正值，對收斂不利，如果將H(i,j)乘上一個數值較小的阻尼因子，令高階項可忽略，就能讓δE仍是負值，這是上述阻尼方法另外形式的表述。

能級移動法則相應于將e(i)-e(j)項替換為e(i)-e(j)-x。在一般情況下由于占據軌道能量低于虛軌道，因此e(i)-e(j)總是負值，此時能量總是降低的，這沒問題。但迭代過程中aufbau定律（即電子按軌道能量從低往上排）并不總是遵守的，這種情況經常出現在過渡金屬絡合物的體系。原因是波函數離收斂比較遠時，分子的平均勢場與真實勢場偏差較大（“真實”是指相對當前步更真實的下一次迭代的平均勢場），此時SCF給出的軌道能量因此與真實偏差較大，甚至順序是倒轉的，這樣程序根據此時SCF軌道能量由低到高判斷出的軌道占據順序就與真實順序不符。當順序相反時，e(i)-e(j)就成為正值，可能δE也會成為正值，導致迭代后總能量不降反升，這光靠乘上阻尼因子已無法解決，此時若用較大的能級提升量x使e(i)-e(j)-x項為負，就能使能量仍然降低。若e(i)-e(j)本來為負值，從中減去了x就會使其絕對值變大，造成δE絕對值減小，使收斂更慢。

3.4 直接最小化(Direct minimization)

這是一類方法而不特指某種具體方法。總能量E是系數矩陣C（或密度矩陣D）的函數，可寫作E(C)，變化軌道組合系數會引起總能量發生變化，所以通過SCF尋找總能量最低的波函數實際上可視為常規的非線性優化問題。這十分類似于分子幾何結構優化，都是在能量面上尋找極小點。幾何優化是不斷調整結構找到總能量最低的構型；而SCF過程則是在特定構型下，不斷優化各軌道組合系數找到總能量最低的體系波函數。所以幾何優化常用的最陡下降法、共軛梯度法、牛頓法、準牛頓法等都可以應用到求解最優電子波函數問題上，MCSCF也利用了牛頓法。而上述的DIIS、EDIIS則反過來也被用在幾何優化上，分別稱GDIIS和GEDIIS，G代表Geometry。使用直接最小化方法每步不再通過對角化Fock矩陣得到新的系數，而是根據不同方法的定義來不斷調整系數。直接最小化能保證一定收斂到一個能量極小點，對于DIIS收斂困難的體系很管用。

直接最小化方法用在求解波函數上雖然要做一些特殊考慮，但大體思路都是一樣的，常見方法包括下列：
3.4.1 最陡下降法：每一步按照當前位置的能量負梯度方向走步，也就是能量下降最快的方向，直到走到這個方向能量最低點。最陡下降法往往在初期能量降低得較快，但由于每步之間方向是正交的，進入能量面呈山谷形狀的區域后容易反復震蕩，收斂緩慢。如果將步幅乘上一個調節因子來減小，避免走到此方向能量最低點，對避免震蕩有好處，稱調整步幅的最陡下降法(scaled steepest descent)。
3.4.2 共軛梯度法：令每次最陡下降法移動的方向通過前一步的方向進行校正，可緩解最陡下降法的震蕩問題，加快收斂。
3.4.3 牛頓法：令E(C)對C進行Taylor展開，舍去三階及以上項作為近似，并根據dE/dC=0的條件，推得C_(n+1)=C_(n)-E(C)'|C_(n) * (E(C)''|C_(n))^(-1)，這里|C_(n)代表將C_(n)代入左邊的變量。顯然對于純二次型區域（是指可以用多變量的二次函數準確描述的區域，三階及以上導數為0）這個表達式是準確的，可以一步達到極小點。在接近收斂時能量面可近似視為二次區域，所以牛頓法收斂很快，但是在優化初期效率很低，一方面是二次區域的近似不再那么合理，相對于DIIS迭代次數可能需要更多；另一方面是牛頓法需要計算E(C)的二階導數，每一步運算耗時很長。
3.4.4 準牛頓法：使用估算方法快速得到E(C)''而不是直接計算，比牛頓法大大節約時間，效率甚至接近DIIS。
3.4.5 增強型牛頓法(Augmented Newton)：單純的牛頓法并不穩定，需要特殊的處理，其中通過修改原始E(C)''來調整步進長度和方向的方法稱為增強型牛頓法。

往往不將C或D作為E的變量，而是將前后兩次MO的變換矩陣X作為變量，好處是MO是正交歸一集，相互獨立利于優化。前面已說到，在對能量有影響的只是X的非對角塊，所以更具體來說E是X[vir(n)_occ(n+1)]的函數。對E(X)的求導可寫為以MO為基的單電子和雙電子積分，由于需要從AO的雙電子積分變換過來，所以很耗時。

將D作為變量時，根據密度矩陣的性質可將能量寫為E(D)=Tr(DF)=Tr(Dh)+Tr(DG(D))，其中h是和G是Fock矩陣中單、雙電子部分。D并不能自由變化，由于MO需要正交歸一，D需要受制于等冪條件DSD=D，否則結果無意義。然而在優化中滿足這個要求并不容易，一種方法是通過3*D^2-2*D^3運算將不等冪的D純化。等冪條件限制對應于軌道占據數為0或者1，而RCA(relaxed constraints algorithms)過程允許在優化過程中放松這個要求，占據數可以在0至1之間，相當于將等冪條件放松為DSD<=D，而到收斂時自動恢復等冪條件，RCA當中最簡單的實現為ODA(Optimal Damping Algorithm)，但速度不如DIIS快。

下面介紹ODA的大致過程。在ODA迭代過程中的任意一段可以這么表示：...D_(n)~-->F_(n+1)~-->D_(n+1)-->D_(n+1)~-->F_(n+2)~-->D_(n+2)...。其中D_(n+1)~=D_(n)~+λ(D_(n+1)-D_(n)~)，λ是通過令能量E(D_(n+1)~)在λ∈[0,1]區間內取得極小值來得到的。由D_(n)~構建F_(n+1)~，再由其本征向量矩陣C_(n+1)構建D_(n+1)，如上方法計算D_(n+1)~，就形成了一次迭代。每一步可以視作由當前位置D_(n)~向D_(n+1)-D_(n)~方向走到能量最低點，這個方向正是能量下降最快的方向。

在Gaussian中可以用SCF選項中的SD和SSD來使用最陡下降法和調整步幅的最陡下降法，但容易出錯，不建議使用。SCF=QC是線搜索法和加入了阻尼的準牛頓法的結合，稱二次收斂法(quadratically convergent)，雖然往往能解決收斂問題，但速度非常慢，另外可能在收斂末期出現振蕩而無法收斂。建議嘗試SCF=XQC，在迭代初期仍用效率高的DIIS方法，難以收斂時才轉為速度慢但收斂性好的SCF=QC，使收斂性和速度都較好，但往往在DIIS的階段就已經報錯停止。

3.5 溫度展寬

也叫分數占據法(FON, Fractional occupation number)，在Gaussian中通過SCF=Fermi使用，默認是關閉的。這個方法的特點是使軌道占據數n出現非整數，其數值根據費米-狄拉克分布計算，即n(i)=1/(1+exp((e(i)-e(F))/kT))，e(F)是費米能級，T是虛構的溫度，密度矩陣D(i,j)=2*∑[k]n(k)C(i,k)C(j,k)，由于不再分清軌道是否占據，k要對所有軌道加和。此方法用的費米-狄拉克分布函數并沒有實際的物理意義，只是由于此方法意圖將一部分占據軌道的電子挪到虛軌道上而借用了這個函數的形式。這樣的結果就是每一步占據軌道和虛軌道之間混合加劇，前面也已提到，這種混合是使能量降低的驅動力，所以此方法可以加快收斂速度。實際應用中，此方法對于容易收斂的體系收斂速度沒什么影響，對一些收斂慢的體系可以大大加快收斂速度，有時也能解決個別體系不能收斂的問題。溫度展寬可以結合DIIS、阻尼一起使用以得到更好效果。

e(F)和T都是可變參數，影響FON的效果。e(F)一般設為HOMO與LUMO能量的平均值，實際上設在HOMO之下LUMO之上都可以。T越高，在費米能級附近的分數占據的軌道越多，軌道混合程度越大；T越低，則越趨于還原為整數占據。對于HOMO-LUMO能隙大的分子，T應設大來以起到明顯的加速效果；而能隙小的分子由于軌道混合已比較重，不需要太高的T。由于單slater行列式要求軌道占據數必須是整數，T在迭代完成前必須回歸為0以滿足整數占據條件，因此T是動態變化的。在迭代開始，視分子的能隙或DIIS誤差矩陣或將溫度設在1000~2000K范圍內，之后逐步降低T。可以線性降低，比如每步100K，也可以將T作為DIIS誤差矩陣的函數，這樣隨迭代進行DIIS誤差矩陣逐漸減小最終趨近0，T也隨之逐漸降至0，對迭代的影響也逐漸減弱，這樣在迭代完成時占據數就自動恢復為整數了。


另外，還有一種通過Car-Parrinello動力學來獲得基態收斂波函數的方法，當其它各種解決波函數不收斂的方法都不靈的時候，這個方法就是最后選擇。但這種方法不能在Gaussian中實現。也就是一開始將電子溫度設得較高，在動力學過程中慢慢降溫，降溫到0時就得到了收斂的基態波函數，也就是電子退火過程。

---

附：對Gaussian中IOp(5/13=1)的大批判

通過IOp(5/13=1)來“解決”SCF不收斂是計算化學界最最最最腦殘的行為，沒有之一！我不吝嗇任何惡毒的詞語來批判用5/13=1的做法。灑家一看到一些人輸入文件里的5/13=1就特別來氣，簡直就像看到大街上有一坨惡翔！

5/13=1意思是，SCF到達指定圈數不收斂也不報錯，繼續做后面的事。為什么用這個IOp腦殘，稍微有一點點量子化學理論常識的人都自然明白。這哪是解決不收斂啊，分明是對不收斂視而不見！完全就是飲鴆止渴！如果到最大圈數時，恰好已經離收斂限不遠了（并且收斂限本身設得并不太松），那么結果可以用，但這必須自己查閱輸出文件來確認這一點；然而，初學者才不知道查閱輸出文件來檢查收斂到什么程度了呢！如果到達最大圈數時能量還在明顯震蕩，由于用了5/13=1又沒報錯，初學者就會傻乎乎地直接用最后給出的能量，然而這能量根本就不靠譜啊！！！用這樣的結果害人害己，輕則結果不準確，重則結果定性錯誤！可以說，對于不知道自行檢查SCF收斂程度的初學者用5/13=1，完全就是對數據不負責任，甚至個別情況可以說是數據造假！

很多對理論知識一無所知的初學者發現平時不收斂會報錯的任務一加5/13=1就不報錯了，于是就傻乎乎地、無條件地所有任務都加上5/13=1，甚至還讓別人加5/13=1，別人甚至連5/13=1都不知道是干什么的，被以訛傳訛也傻乎乎地用，這種膚淺至極的做學問的態度真是不忍直視！！！不懂某個關鍵詞是干什么的話就不要加這個關鍵詞，不要盲目地效仿別人（尤其是網上大量計算化學的討論都是錯的），不要被以訛傳訛！

在http://bbs.keinsci.com/forum.php?mod=viewthread&tid=3344帖子中Gaussian官方客服也已經明確表態絕不應當在任何實際計算中用5/13=1。

順帶一提，有人根據自己的一些計算經驗，認為5/13=1對于解決幾何優化過程的SCF不收斂是有益的。這種觀點的原理是：幾何優化初期由于結構不太好，SCF不容易收斂，于是加了5/13=1讓程序不報錯；等到優化末期，結構收斂到差不多極小點結構附近了，SCF通常也就容易收斂了，此時雖然帶著5/13=1，但是照樣也能收斂到了標準收斂限，大不了人工檢查一下是否SCF收斂了。實際上，即便是出于這種考慮，用5/13=1也是不對的！這里給個例子：用5_13=1有害的例子.rar
從例子可見，用了5/13=1時，由于波函數未收斂從而進一步導致受力計算爛，而又不試圖通過恰當解決SCF不收斂的做法去解決，導致結構也反復震蕩難以收斂。而按照本文做法，加了輔助SCF收斂的關鍵詞后，每一步SCF都收斂了，結構也順利收斂了。

在我來看，能用5/13=1的僅限一種情況：你已經是量化的老司機，很清楚自己在干什么，在對某個初始結構偏離極小點可能較多的體系做幾何優化時，發現在優化的第一步SCF不收斂，而且你已經毫無遺漏、費勁千辛萬苦去嘗試了本文所有辦法，但就是死活達不到SCF收斂限，在絕望之際才可以用一次5/13=1。別忘了，如果優化最終順利結束了，一定要檢查最后一步的SCF收斂情況，看是否達到了收斂限，并且確保沒虛頻，還要用stable關鍵詞做波函數穩定性檢測確保收斂到的波函數是穩定的。一般研究極少遇到這種情況，逼不得已這樣被迫用5/13=1的往往也就是過渡金屬團簇、用小核贗勢計算鑭系錒系體系（見《使用Gaussian做鑭系金屬配合物的量子化學計算》http://www.shanxitv.org/581）等SCF特別難收斂的情況。注意從G16開始，起碼截止到目前最新的G16 C.01來說，如果SCF沒收斂到默認的收斂精度，幾何優化一步都不給你做，所以用這個做法只能用之前的Gaussian版本。