Gaussian中幾何優化收斂后Freq時出現NO或虛頻的原因和解決方法

Gaussian中幾何優化收斂后Freq時出現NO或虛頻的原因和解決方法

First release: 2015-Jan-16 Last update: 2023-Feb-18

1 前言

眾所周知，做振動分析之前必須先在相同級別下優化。總有Gaussian使用者問，優化收斂了，也顯示四個YES了，怎么做freq之后四個收斂判斷標準卻有的出現NO，甚至還出現虛頻了？此文件就專門說一下。Gaussian中opt=TS方式的搜索過渡態的算法和幾何優化很相似，因此下面的討論對于過渡態的優化也完全適用，只不過要消的虛頻顯然不是感興趣的那個（看振動動畫可判斷），而是額外存在的多余的虛頻。

特別要強調的是，絕對不存在消不掉的虛頻，筆者從來沒見過例外。不要在網上問我怎么按照此文的做法還是沒消掉虛頻，消不掉一定是沒認真領會此文并充分嘗試。

2 關于幾何優化收斂的判斷和出現虛頻的原因

牛頓法是高效、最常用的尋找多元函數駐點的方法，在量化研究中，這個方法被廣泛用來搜索過渡態和優化幾何結構。然而，牛頓法的每一步需要計算一次Hessian矩陣（更具體來說是其逆矩陣），也就是能量對幾何坐標的二階導數矩陣，如果進一步考慮質量權重后就是所謂的力常數矩陣。對于較高計算級別或大體系，計算一次Hessian矩陣是相當耗時的，另外很多方法在Gaussian里沒有二階解析導數（如CCSD、MP4等），對于這樣的方法計算大體系的Hessian矩陣更是耗時甚巨甚至無法承受。因此，為了節省優化時間，一般量子化學程序默認用的都是準牛頓法，也就是每一步不精確計算Hessian，而是只計算受力，通過受力和上一步的Hessian矩陣來近似得到這一步的Hessian矩陣。雖然準牛頓法由于用的是近似的Hessian矩陣，比起基于精確Hessian的牛頓法需要更多步數才能找到過渡態/極小點，但由于每一步的耗時低得多（比如準牛頓法的10步耗時可能只相當于牛頓法的1、2步），所以整體還是比牛頓法便宜得多得多，所以是大多數量化程序默認的優化方法。當Gaussian發現當前結構的受力以及位移（下一步結構相當于當前結構的變化）都收斂了，優化就停止。受力和位移都按照最大值和方均根（RMS）來檢測，因此共有四個判斷標準，即四個都YES時就宣告收斂。

做頻率計算，則必須先精確計算一次Hessian矩陣（或者讀取精確的Hessian矩陣），才能經過一番處理得到頻率，參見《基于fch中的Hessian矩陣計算振動頻率的簡單程序Hess2freq》（http://www.shanxitv.org/328）里的介紹。freq任務之后，Gaussian會自動根據當前結構的受力和現成的精確Hessian矩陣也用上述四個判斷標準考察一下當前結構是否真的準確收斂到了極小點。

正因為幾何優化默認用的是近似的Hessian矩陣，而頻率計算用的是精確Hessian矩陣，因此判斷收斂時，在后兩項，即最大位移和均方根位移上的判斷會存在差異（受力那兩項是一致的）。常見的情況是優化的最后一步四個標準都顯示YES了，但是到了freq任務時，基于精確Hessian矩陣判斷則發現后兩個并非都YES，甚至出現了虛頻，這說明優化的精度不是很高。特別是虛頻如果較大的話，更是體現當前結構偏離極小點可能挺明顯。

下面是個典型的例子，opt最后的輸出為
Maximum Force            0.000401     0.000450     YES
RMS     Force            0.000197     0.000300     YES
Maximum Displacement     0.001791     0.001800     YES
RMS     Displacement     0.000951     0.001200     YES
Freq最后的輸出為
Maximum Force            0.000401     0.000450     YES
RMS     Force            0.000197     0.000300     YES
Maximum Displacement     0.002231     0.001800     NO
RMS     Displacement     0.001325     0.001200     NO

那么，opt最后都YES了但freq最后沒有都YES，結果還能用么？筆者認為，如果數值比默認情況的收斂限高不太多，不超過2倍，且沒有虛頻出現，一般還是可以接受的。比如上面的情況就基本可以接受。但是，如果你對頻率/結構的準確度要求高，或者都出現了虛頻，則應當做更精確的優化。

3 提升優化和數值精度解決虛頻的做法

下面介紹讓優化更精確來避免freq最后出現NO乃至虛頻的做法。

3.1 幾何優化時用更嚴的收斂限

在opt里面寫上tight可以將默認的優化收斂限設嚴一個數量級以上。由于極小點定位得更精準了，出現虛頻的可能性自然大大降低了。但代價是達到收斂限的難度大大增加，因此優化達到收斂所需的步數會大大增加、耗時大大提升，還增加了因為震蕩而一直達不到較嚴格的收斂限的概率。

注意如果寫了opt=tight freq，則freq最后也會用tight標準來判斷是否收斂，而實際上freq最后只要能滿足默認收斂限判斷標準就夠了，不需要非得滿足tight那么嚴苛的判斷標準。

值得一提的是優化過程中即便位移沒收斂，但只要受力小于收斂限的100倍也自動被Gaussian算作收斂。我認為這主要是考慮到勢能面非常平緩的柔性大分子或分子團簇，相對于其體系的尺度，把結構收斂到那么精確意義也不大。這種情況下，連opt最后都有NO，顯然freq最后也肯定有NO。若對結構精度要求高的話可以把優化的收斂限設嚴。

3.2 提升DFT積分格點精度

對于DFT計算要考慮交換-相關泛函的積分格點質量問題。如果不了解什么是DFT的積分格點，可參看《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69）。如果積分格點質量不夠理想，DFT計算的受力和Hessian的精度就不夠高，這往往是造成虛頻出現的原因。此時即便用了tight收斂限也解決不了問題，而且改用tight后可能還極難收斂、容易震蕩。如果你是Gaussian 09的用戶，用opt=tight時應同時搭配int=ultrafine使用比默認的int=fine更高精度的積分格點。從Gaussian 16開始int=ultrafine已成為默認，因此就不需要手動寫了。

對積分格點的要求與使用的泛函、體系都有密切聯系。一些明尼蘇達系列泛函（比如M06-2X等，但不是所有）對于積分格點精度要求比B3LYP、PBE0、TPSS等其它常見泛函高得多。因此比如M06-2X在int=ultrafine下依然有時會因為積分格點精度不夠而導致很微小的虛頻的出現，此時可嘗試更高檔次的積分格點int=superfine，但這非常昂貴。

對于DFT計算極個別色散作用主導的弱相互作用體系，有時候更需要用int=superfine。比如筆者曾經在ωB97XD/def2-TZVP(-f)下計算C30嵌套結構，老是有二點幾波數的微小虛頻死活消不掉。最終有一次我雖然沒設更嚴收斂限，但是用了int=superfine，可算把虛頻給搞掉了。但這種情況極度罕見，切勿視為普遍情況。只是說對于某些弱相互作用體系，哪怕不是明尼蘇達系列泛函，當你碰到極微小虛頻且實在干不掉時，還是值得得試一把int=superfine。

3.3 使用精確的Hessian矩陣

在opt里寫上calcall，那么幾何優化過程中每一步都會精確計算Hessian，自然最終能優化得比基于近似Hessian的準牛頓法更準確。此時幾何優化和freq所用的Hessian都是一致的，因此對收斂的判斷結果也是完全相同的。也就是說，只要opt最后是YES則freq最后也必然是YES，并且多數情況可以確保無虛頻（值得一提的是，用calcall的時候，優化任務結束后會自動做振動分析，因此完全沒有必要再單獨寫freq關鍵詞要求做振動分析了）。

不過，每一步優化都精確計算Hessian太耗時。如果你之前已經在默認情況下優化過，可以取最后的結構重新優化，并在opt里寫上calcfc，這樣只在優化第一步的時候精確計算Hessian，而之后還是近似估計Hessian，這樣收斂后也往往可以避免freq之后出現NO和虛頻。但是如果初始結構離實際極小點太遠則這么做無效，因為等到結構收斂時Hessian矩陣可能又偏離精確Hessian比較遠了。

從Gaussian 16開始，opt里加入了一個很有用的選項recalc。比如opt=recalc=n就代表第一步精確計算一次Hessian矩陣，之后每n步重新計算一次Hessian矩陣。這種做法的效果和耗時都介于calcfc和calcall之間。對計算Hessian成本不很高的情況，n取3~5左右比較合適，對計算Hessian成本較高的情況，n取5~10比較合適，能達到效果和耗時的權衡。讀者請隨機應變。

4 其它導致虛頻的因素和解決方法

4.1 結構對稱性問題

有時候虛頻是由于優化所用的初始結構對稱性太高導致的。比如聯苯在基態下實際上兩個苯環間是有一定二面角的，但優化聯苯所用的初始結構如果是純平面的，則優化過程就會一直保持平面狀態，最終也得到平面構型，顯然會發現有對應于兩個平面彼此間扭轉的虛頻。碰到這種情況，最佳的解決方法是在GaussView觀看振動模式的界面中選中虛頻模式，選上Manual Displacement并拉動滑條來沿著虛頻模式略微調整結構，然后點save structure得到新的結構，再拿這個新結構保存新的輸入文件，重新優化和做振動分析后通常就已經沒虛頻了。有時候對于大體系，也可能因為某些局部區域的初始對稱性太高，一直維持到最后而出現虛頻，這時候也應按照虛頻調節結構。對于其它情況，這么按照虛頻調節結構的做法解決虛頻的幾率比較有限，雖然也可以嘗試一下。而對于有多個虛頻的情況，想通過這種做法一次性就把所有虛頻都解決沒太大可能。我經常看到網上有很多人但凡只要看到虛頻就總是妄圖通過這種按照虛頻調結構的辦法試圖消掉，這是明顯不對的！！！真正應該優先嘗試的是本文前面說的做法！

4.2 SMD隱式溶劑模型問題

SMD隱式模型流行度很高，在《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）我專門介紹過。我強烈不建議對opt和freq任務用SMD模型，因為此模型下很容易由于數值噪音而產生虛假的虛頻，而且還可能由此造成幾何優化難收斂。因此當優化和振動分析必須在溶劑模型下做時，強烈建議用Gaussian的scrf關鍵詞默認的IEFPCM溶劑模型。從幾何優化結果來講，SMD相較IEFPCM也沒有絲毫額外的好處，如327號博文所述，SMD比IEFPCM長處在于額外多考慮了非極性部分，對計算溶劑環境下的能量很重要，而非極性部分對優化的結構影響甚微。所以，opt freq在IEFPCM下做，而之后算能量時再在SMD下做，是比較得當的做法。

4.3 IOp定義泛函的問題

在Gaussian中可以自定義泛函，見《Gaussian中非內置的理論方法和泛函的用法》（http://www.shanxitv.org/344）和《在Gaussian中自定義范圍分離泛函的方法》（http://www.shanxitv.org/550），這需要用到IOp關鍵詞。很多人用自定義泛函的時候還是按習慣寫上opt freq關鍵詞，這是不行的。要知道，IOp是沒法傳遞到下一個任務的，結合opt freq用的時候只對opt有效，而對freq無效，自然會導致振動分析和優化時用的泛函定義不同、對應的勢能面不同，故很可能因此出虛頻。所以用自定義泛函的時候opt和freq必須分開做。

4.4 程序bug問題

Gaussian從G09 D.01版開始加入DFT-D3，相關介紹見《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）、《DFT-D色散校正的使用》（http://www.shanxitv.org/210）、《亂談DFT-D）《http://www.shanxitv.org/83》。DFT-D3有零阻尼和BJ阻尼兩個版本，分別用em=gd3和em=gd3bj關鍵詞使用。后者最常用，因為算相互作用能的精度整體更好一點。BJ阻尼版本在G09 D.01里的解析二階導數存在bug，會導致你結合DFT-D3(BJ)做優化完全精確的情況下，結合DFT-D3(BJ)做的振動分析也可能出現虛頻。因此解決辦法是用G09 E.01及以后的版本，或者用零阻尼版本DFT-D3(0)做opt和freq。

PS1：準牛頓法、牛頓法在不同量子化學程序里具體實現不同，有很多數值方面的技巧和額外考慮，本文中沒有細談。想了解更多的話可以看看《過渡態、反應路徑的計算方法及相關問題》（http://www.shanxitv.org/44）以及《幾何優化、過渡態搜索、IRC綜述與原文合集》（http://bbs.keinsci.com/thread-105-1-1.html）當中的詳細討論。Gaussian做opt默認用的方法并不是簡簡單單的準牛頓法，而是牽扯到很多細節，在手冊的opt部分里有專門的描述。

PS2：有人問，為什么有的時候優化時明明用了calcall，但是基于其結構再做freq任務的時候顯示的判斷標準卻依然和優化最后一步不同，甚至又出現了NO。可能原因有三 (1)優化時用的級別、數值設定（如DFT積分格點）等因素和幾何優化時不完全一致 (2)用GaussView讀取優化的輸出文件，然后又保存成freq任務的輸入文件的時候，由于小數位數有限，因此造成了一點數值誤差 (3)幾何優化的時候用了GDIIS（默認的GEDIIS也有一定GDIIS的成份），由于這種方法預測下一步位移的時候還會參考之前步的信息，因此會和freq時基于牛頓法判斷出的位移有所不同。