技術標籤：gaussian 和gaussview

筆者經常碰到小夥伴在用Gaussian軟體計算涉及自由基的反應時，不清楚何時該加關鍵詞guess=mix，何時不該加；也可能會有師兄/老師這樣告訴新手：碰到自由基一律用guess(mix,always)。前者可能量化基礎不紮實，碰到這類問題不懂；後者則可能缺乏實際計算經驗。趁假期有空，正好寫上一篇，詳細解釋一下。當然，筆者寫的絕對不是標準答案，只能力求合理性和正確性，僅供對這個問題不清楚的小夥伴們參考。

為方便起見，本文僅討論RHF和UHF，所有內容對DFT同樣適用。對UHF計算不熟悉的新手可以閱讀本公眾號發過的軟體教程《用Gaussian做UHF計算》。本文中所涉及計算皆使用G16 A.03。

首先給出簡單的結論：(1)guess=mix只在自旋極化單重態(即使用UHF方法在單重態下做計算，發現有嚴重自旋汙染)時需要考慮加，其他情況(如二重態、三重態)等無需考慮這個問題；(2)always表示在結構優化的每一步中都執行guess=mix。顧名思義，這隻在結構優化中可能有用，而在單點計算中無需加、加了也沒用。

接下來我們一一解釋。關於自旋極化單重態和guess=mix的含義，Sob老師的博文《談談片段組合波函式與自旋極化單重態》(http://sobereva.com/82)講得十分詳細，強烈推薦經常做此類計算、但又還沒看過博文的小夥伴仔細閱讀。這裡再舉一個簡單的單點計算例項：在UHF/STO-3G水平下計算鍵長為2.0 Å的H₂

(1) #p RHF/STO-3G nosymm(2) #p UHF/STO-3G nosymm(3)#pUHF/STO-3Gnosymmguess=mix

這樣的UHF計算和結果顯然沒有太大價值，多花了計算時間卻得到與RHF一樣的解。而guess=mix的作用就是混合HOMO與LUMO軌道，使UHF的初始軌道與RHF不同，這樣迭代就有可能收斂到更低的解上。

寫法(3)得到的結果為−0.937213 a.u.， = 0.946，能量明顯低許多，代價是自旋不嚴格等於零，即不再是純態。開啟其log檔案搜尋Initial guess，會發現波函式初猜就已經不是純態 = 1.0。

mix具體的做法並沒有統一規定，可以是隻混合alpha列的HOMO與LUMO軌道，也可以是alpha、beta兩列的HOMO與LUMO軌道各自混合、但混合方式不一樣。從Sob博文中展示的例子可以看出高斯是採用後者。當然，我們也可以自己手動操作一番，加深理解。例如，將下列計算產生的fchk檔案開啟

h2_uhf_only.gjf檔案內容：

%chk=h2_uhf_only.chk#p UHF/STO-3G nosymm guess=(only,save)Title Card0 1H   0.0   0.0   0.0H   0.0   0.0   2.0

h2_uhf_only.fchk檔案部分內容(修改前)：

Alpha MO coefficients                      R   N=           4  6.68386924E-01  6.68386924E-01  7.53443037E-01 -7.53443037E-01Beta MO coefficients                       R   N=           46.68386924E-016.68386924E-017.53443037E-01-7.53443037E-01

h2_uhf_only.fchk檔案部分內容(修改後)

Alpha MO coefficients                      R   N=           4  1.00538561E+00 -6.01437543E-02 -6.01437543E-02  1.00538561E+00Beta MO coefficients                       R   N=           46.68386924E-016.68386924E-017.53443037E-01-7.53443037E-01

容易看出，上述改動僅是將alpha一列軌道的HOMO與LUMO軌道混合了一下。更具體地說，就是按照如下公式進行混合(此即等比例混合)：

再使用命令

unfchk h2_uhf_only.fchk h2_uhf_only.chk

將修改後的fchk轉回chk檔案。將gjf檔案內關鍵詞guess=(only,save)改為guess=read，重新提交計算，即可得出較低的UHF能量−0.937213 a.u.。當然，等比例混合只是一種做法，也可以用 (根號3)/2、1/2這種搭配，保持軌道的正交歸一性即可。高斯內部原始碼實現的就是這種功能，只不過我們舉了個最簡單的例子，用計算器就能算。

但若混合比例很離譜，比如在這個例子裡直接把HOMO、LUMO軌道互換(這相當於把波函式直接置於鞍點上)

h2_uhf_only.fchk檔案部分內容(再一次修改)：

Alpha MO coefficients                      R   N=           4  7.53443037E-01 -7.53443037E-01  6.68386924E-01  6.68386924E-01Beta MO coefficients                       R   N=           46.68386924E-016.68386924E-017.53443037E-01-7.53443037E-01

再轉化為chk檔案，使用guess=read提交計算，也能立即收斂出一個能量−0.665399 a.u.，但明顯高很多，說明這種混合不合理。

我們還可以舉出很多guess=mix有效的例子，如單重態O₂分子、單重態卡賓、並苯體系(並的環越多，自旋汙染越大)等等。

高斯的guess=mix雖然很傻瓜、黑箱化，很多時候很好用，但它不能保證沒有其他更低的UHF解，也不能保證收斂到一個穩定的UHF波函式，少數時候甚至連對稱性破缺波函式都得不到。這些情況筆者在平時計算中也見過幾例：

(1)O₂在鍵長1.25 Å時，僅使用

#p UHF/cc-pVDZ nosymm guess=mix

得到的電子能量為−149.589479 a.u.，對應=0.627，確實是對稱破缺的，有自旋汙染。然而如果加上關鍵詞stable=opt算，會發現輸出檔案內提示內部不穩定性

The wave function has an internal instability

並自動使用二階優化方法繼續優化至穩定，能量為−149.596725 a.u.，對應=1.015，能量更低，自旋偏離純態更多。實際上在這個例子裡，O-O鍵長大於1.16 Å的幾乎所有點都會有這個問題。

(2)自由基引發劑tBuOOH(叔丁基過氧化氫)分子，在M062X/6-31G(d)水平下優化得平衡幾何結構，然後將O-O鍵長拉至1.90 Å，僅使用關鍵詞

#p UM062X/6-31G(d) nosymm guess=mix

算單點得到的電子能量為−308.590951 a.u.，對應 = 0，顯然是收斂到了RM062X的解上。但是一旦加上stable=opt，會發現能量降為−308.594232 a.u.， = 0.439。

所以筆者的建議是，用guess=mix時總是加上stable=opt，即檢驗波函式穩定性，若發現不穩定則優化至穩定。

如果僅寫stable，表示只檢驗波函式穩定性，若穩定則結束；若不穩定只會輸出不穩定資訊然後結束，不會繼續優化至穩定。這時候還要手動補交guess=read stable=opt任務，相當於波函式穩定性檢驗了兩次，多費了些時間，因此不如直接用stable=opt簡單。

有三點要注意：(1)高斯中stable=opt不能與結構優化，限制性優化或IRC等一起使用(除非自己寫指令碼實現，但可以預見計算量會很大)；(2)stable=opt只能保證得到穩定的波函式，對一些極特殊情況仍不能保證沒有其他更低的UHF解，這一點以後再展開；(3)沒有opt=stable這種東西。

諸如二重態、三重態這些情況，alpha、beta兩列初始軌道可以一模一樣，但由於兩列軌道中的電子數不一樣，對應的密度矩陣也不一樣，因此不存在上述問題，無需guess=mix。

Reference

1. 《談談片段組合波函式與自旋極化單重態》http://sobereva.com/82

2. https://chemistry.stackexchange.com/questions/42005/initial-guess-for-unrestricted-hartree-fock-calculation