計算適用于OPLS-AA力場做模擬的1.2*CM5原子電荷的懶人腳本
計算適用于OPLS-AA力場做模擬的1.2*CM5原子電荷的懶人腳本
文/Sobereva@北京科音 2021-Jan-28
根據J. Phys. Chem. B., 121, 3864 (2017)中的測試,基于OPLS-AA力場的模擬很適合結合1.2*CM5原子電荷。雖然用1.2*CM5的時候水合自由能的計算精度稍遜于文中測試的另一種原子電荷1.14*CM1A-LBCC(因為這種原子電荷本來就是針對水合自由能訓練的校正參數),但在計算其它屬性上,如蒸發焓、密度,用1.2*CM5時的誤差都更小,而且經驗性更低、更有普適性。所以當大家用OPLS-AA力場研究新的小分子的時候,我比較推薦用1.2*CM5電荷。1.2*CM5電荷就是在Truhlar等人提出的原始的CM5電荷基礎上乘上1.2,這相當于增大了原子電荷的數量級,等效地體現出了溶劑環境對溶質的極化作用。
筆者之前在《計算RESP原子電荷的超級懶人腳本(一行命令就算出結果)》(http://www.shanxitv.org/476)和《RESP2原子電荷的思想以及在Multiwfn中的計算》(http://www.shanxitv.org/531)中分別給出了超級懶人計算RESP和RESP2電荷的Linux腳本,完全不會用Gaussian的人都能輕松計算,已經有不少人在用。最近碰到完全不會用Gaussian的人試圖計算1.2*CM5電荷向我求助,我遂又寫了個類似的計算1.2*CM5電荷的懶人Linux腳本,在這里說一下。
這個腳本是Multiwfn文件包里的examples\scripts\1.2CM5.sh,在2021-Jan-28及之后更新的Multiwfn中才有,Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn免費下載。
此腳本使用很簡單:先確保Gaussian在當前機子里已經裝了,見《Gaussian的安裝方法及運行時的相關問題》(http://www.shanxitv.org/439),也確保將Multiwfn按照手冊2.1.2節的說明在機子里裝了。之后假設1.2CM5.sh和一個結構文件phenol.xyz都在當前目錄下,只需要運行./1.2CM5.sh phenol.xyz,之后在屏幕上就會看到腳本的運行過程:
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
Running optimization task via Gaussian...
Done!
Running formchk...
Running Multiwfn...
Finished! The optimized atomic coordinates with 1.2*CM5 charges (the last column) have been exported to phenol.chg in current folder
最終得到的phenol.chg用文本編輯器打開后可見,其中2、3、4列是優化后的XYZ坐標(埃),最后一列就是1.2*CM5電荷了。
腳本原理:這個腳本會將輸入文件里的結構作為初猜結構用Gaussian在B3LYP-D3(BJ)/def2-SVP下做幾何優化,然后調用formchk將得到的chk轉化為fchk文件,然后調用Multiwfn計算CM5電荷,最后給出1.2*CM5電荷。
幾點相關事項:
? 用的輸入文件可以是任意含有結構信息而且Multiwfn支持的文件格式,比如pdb/mol/mol2/xyz/fch/gjf/wfn等等等等,見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/379)。
? 默認情況下電荷當成0,自旋多重度當成1。如果是比如不帶電的二重態體系,應當寫比如./1.2CM5.sh phenol.xyz 0 2。
? 如果運行腳本時提示沒可執行權限,先運行chmod +x ./1.2CM5.sh
? 此腳本調用的是Gaussian 16,如果當前機子里裝的是諸如Gaussian 09,需要將文件的第22行改為Gaussian=g09。
? 從腳本里的關鍵詞可見,計算是在真空下進行的,這是刻意而為之,不要自行改成在溶劑模型下進行。
? 如果腳本運行失敗,有這么幾種可能:(1)Gaussian根本沒恰當安裝,應當確保能手動調用Gaussian運行一個最簡單的計算任務 (2)Multiwfn沒按照手冊2.1.2節恰當安裝。其中有些圖形庫不方便裝的話,可以用Multiwfn的noGUI版,此時不需要裝圖形庫 (3)Gaussian計算失敗,可能是版本太老而不支持em=GD3BJ關鍵詞,可將腳本中此關鍵詞去掉。也可能幾何優化或SCF不收斂,注意看計算中途產生的gau.out文件內容判斷原因。
? 如果你手頭已經有Gaussian或其它程序產生的Multiwfn支持的波函數文件了(fch/wfn/wfx/molden/mwfn等),就沒必要再用這個腳本了。直接載入Multiwfn,依次輸入7、16、1就能得到CM5電荷,再手動乘上1.2即可。
? 前述J. Phys. Chem. B.文章中只是考慮了中性體系,沒有考慮帶電體系。帶電體系用什么原子電荷結合OPLS-AA沒有統一說法,一般情況下不能用1.2*CM5電荷,因為此時總電荷都不是整數了。筆者建議此時用Multiwfn算RESP2電荷,雖然與OPLS-AA兼容性沒有充分測試,但原理上問題不大。
? 如果使用此腳本計算1.2*CM5電荷,請按照Multiwfn程序包里How to cite Multiwfn.pdf文檔的說明恰當引用Multiwfn。
如果你沒買Gaussian,也可以用免費的ORCA量子化學程序結合Multiwfn算1.2*CM5原子電荷,筆者也提供了相應的傻瓜式腳本,見《ORCA結合Multiwfn計算RESP、RESP2和1.2*CM5原子電荷的懶人腳本》(http://www.shanxitv.org/637)。