令Gaussian 16中SCF未收斂到默認收斂限也能繼續做后續計算的方法
令Gaussian 16中SCF未收斂到默認收斂限也能繼續做后續計算的方法
文/Sobereva@北京科音 2021-Nov-17
注：本文內容對目前我用過的Gaussian 16 A.03、B.01、C.01都適用，對未來的版本不一定。
Gaussian程序從09版開始默認的收斂限成為了tight，相當于SCF=conver=8關鍵詞設的收斂限，這在所有量子化學程序里幾乎是最嚴的。Gaussian程序從16版開始加入了一個規則，就是如果SCF過程中如果SCF收斂精度沒達到默認的收斂限，而且你的任務涉及到計算能量的導數，如幾何優化、振動分析，或者用的是后HF、TDDFT等涉及到多行列式/多組態的方法，在SCF模塊運行完之后就會有下面的提示，導致程序以報錯方式終止
SCF Error SCF Error SCF Error SCF Error SCF Error SCF Error SCF Error SCF Error
ERROR!!!!
SCF has not converged. Gradients and post-SCF results would be GARBAGE!!
Gaussian 16做這么一個檢查本身沒什么問題，因為如果SCF收斂精度如果不夠高，則上述那些任務、方法的計算精度可能會比較糟糕，而且這也避免了一些菜鳥被其它菜鳥誤導而亂用IOp(5/13=1)還不知道自己檢查收斂情況。Gaussian 16的這個收斂精度檢查設計的最變態的地方是沒法關閉，程序手冊里和IOp文檔里對關閉方法都只字未提，這給許多用戶的實際研究帶來極大的不便。比如很多SCF特別難收斂的體系，如過渡金屬團簇、用小核贗勢的鑭系錒系配合物（《使用Gaussian做鑭系金屬配合物的量子化學計算》http://www.shanxitv.org/581里提到了）、加較大外電場時，本來SCF收斂到默認的很嚴的收斂限就普遍較難，而對于幾何優化任務，由于初始搭的結構往往不夠理想（偏離極小點較遠）導致電子結構更復雜，SCF收斂到默認收斂限的難度就更大，很多情況下即便使用《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61）里的諸多做法嘗試很久也沒法解決。這種情況，一個常見解決策略是先用相對較松的SCF收斂限比如SCF=conver=6做粗略的幾何優化。當體系結構離極小點結構比較近的時候，SCF就相對容易收斂到默認收斂限了，因此之后可以再用默認的收斂限進一步優化（如果你是老司機而且永不會忘記檢查SCF收斂情況的話，IOp(5/13=1)也不是一定不能用）。然而由于Gaussian 16上述變態設定，導致這重要的技巧都無法實現了，令我感到Gaussian開發者嚴重脫離廣大群眾的應用場景。有的時候我為了用上述技巧我都不得不刻意改用Gaussian 09。還有時候，要對體系做TDDFT、雙雜化、CCSD(T)等計算，但個別情況SCF收斂到SCF=conver=8就是特別困難，而用#P監控SCF過程發現收斂到SCF=conver=7的程度是可以達到的，此時做這些計算的精度并沒什么明顯問題，可Gaussian 16來了個一刀切，這種情況下就是死活不讓你做SCF之后的計算，真是特別荒誕！
今天在CCL上有人分享了一個解決上述問題的奇技淫巧，令SCF沒收斂到默認收斂限也能做后續計算，這對Gaussian用戶極為重要，我覺得很有必要在這里具體說一下。不了解下述的Gaussian的Link、IOp的話看《Gaussian的Link、IOp與非標準計算路徑》（http://www.shanxitv.org/57）。
從Gaussian 16開始，在調用Link 701、Link 801模塊開始時都會調用SCFChk做一下SCF收斂性的檢查，沒收斂到默認收斂限就會出現前述報錯提示并終止任務。Link 701是計算單電子積分的一階或二階導數的模塊，所有算能量的導數的任務都會經歷它。Link 801是初始化雙電子積分變換，TDDFT、雙雜化、后HF等涉及到多行列式/組態函數的任務都會經歷它。這個Gaussian 16引入的新設計導致了前述問題。
實際上，可以通過官方IOp手冊（http://gaussian.com/iops/）里都只字未提的選項來關閉SCF收斂性檢查（估計CCL上那個人是看到了一般用戶都拿不到的Gaussian 16源代碼分析出來的）。對于需要經歷Link 701的任務，可以用IOp(7/127=-99)來關閉，對于需要經歷Link 801的任務，可以通過IOp(8/117=-99)來關閉。實際上不是必須=-99，只要是比-100更正的值都可以。
因此，在Gaussian 16中用比如# B3LYP/6-31G* SCF=conver=6 opt，由于負責SCF的Link 502模塊收斂到較松的收斂限就結束了，因此在計算幾何優化要用的受力時由于發現SCF沒收斂到默認的SCF=conver=8的程度，會報錯并出現前述提示。而如果你寫# B3LYP/6-31G* SCF=conver=6 opt IOp(7/127=-99)，就可以讓幾何優化進行下去，和Gaussian 09的情況一樣。
如果你嫌上述IOp不好記，每次寫的時候也麻煩，有一個技巧可用：如《Gaussian的安裝方法及運行時的相關問題》（http://www.shanxitv.org/439）所述，可以在Default.Rou（Windows版）或Default.Route（Linux版）里面寫上默認的計算資源設置，實際上默認的關鍵詞也可以寫進去。如果你在里面加上一行-#- IOp(8/117=-99) IOp(7/127=-99)，則這兩個IOp設置會對之后所有計算都會默認啟用，以后就再也不會遇到前述的煩人的問題了！
還有一點要注意，根據體系的不同，大多數時候Gaussian默認在SCF計算前期用精度較低的方式算電子積分（SCF一開始會看到Integral accuracy reduced to 1.0D-05 until final iterations.的提示。若強行要求總是這么做可以寫varacc關鍵詞），此時SCF=conver=x的設置其實是不起實際作用的，雖然在SCF開始之前顯示的收斂限和自己設的一致，但實際還是會迭代到滿足tight收斂標準才結束。如果想讓SCF=conver=x總是如實生效，應當在SCF里同時寫上novaracc。即前例建議改為# B3LYP/6-31G* SCF(conver=6,novaracc) opt IOp(7/127=-99)，可確保每一輪如實迭代到conver=6的程度就視為收斂。如果你想讓novaracc默認啟用的話，可以在Default.Rou/Route里的-#-后面加上scf=novaracc。
最后對初學者們強調一點，看本文絕對別斷章取義！絕對不要把SCF=conver=6時幾何優化的結果當做最終用的結構！SCF=conver=6情況下的優化只能算得上預優化，實際發表用的結構要在默認的SCF收斂限下優化得到！