計算RESP原子電荷的超級懶人腳本(一行命令就算出結果)
計算RESP原子電荷的超級懶人腳本(一行命令就算出結果)
文/Sobereva@北京科音
First release: 2019-Apr-17 Last update: 2022-Aug-6
RESP電荷主要優點是特別適合對柔性分子做分子動力學計算。之前筆者寫過《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》(http://www.shanxitv.org/441)一文,極其詳細地介紹了RESP電荷,并結合實際例子充分展現了Multiwfn在計算RESP電荷上的強大性,可謂極度便利而且極度靈活,明顯比Antechamber等其它能計算RESP的程序好用、方便得多得多得多。但在網上的答疑過程中,筆者也看到很多要計算RESP電荷的人完全是量子化學外行,甚至GaussView都不會用,Gaussian的關鍵詞一點也不會寫,需要一個極度傻瓜化的“一鍵”工具來計算RESP電荷,為此筆者在這里介紹一個基于Gaussian和Multiwfn的僅僅需要寫一條命令就可以算出RESP電荷的Linux環境下的腳本,可謂將RESP電荷的計算極盡簡化,不可能做到更簡單了。
另外,筆者在《RESP2原子電荷的思想以及在Multiwfn中的計算》(http://www.shanxitv.org/531)中介紹了RESP2電荷,在文中也提供了和本文使用方式基本一樣的腳本RESP2.sh來專門用于一鍵計算RESP2電荷。如果你計算RESP電荷的目的是用于分子動力學模擬,筆者更建議用RESP2來算,有較大可能結果更好。
如果你沒買Gaussian,也可以用免費的ORCA量子化學程序結合Multiwfn算各種原子電荷,筆者也提供了相應的傻瓜式腳本,見《ORCA結合Multiwfn計算RESP、RESP2和1.2*CM5原子電荷的懶人腳本》(http://www.shanxitv.org/637)。
如果你想了解本文提供腳本原理是怎么回事,強烈建議仔細閱讀《詳談Multiwfn的命令行方式運行和批量運行的方法》(http://www.shanxitv.org/612),里面將一切細節都解釋得特別清楚,并能充分體會到寫腳本運行Multiwfn多么重要和便利。
1 準備工作
在機子里安裝Gaussian。不會裝的話仔細看《Gaussian的安裝方法及運行時的相關問題》(http://www.shanxitv.org/439)。去http://www.shanxitv.org/multiwfn下載Multiwfn最新版本。然后嚴格按照手冊2.1.2節說的方式安裝到Linux下面,從而能通過輸入Multiwfn命令啟動程序。
對Multiwfn目錄下的examples\RESP\目錄下的RESP.sh根據需要進行恰當的修改。此腳本默認調用的是g09,Gaussian16的用戶應當將其中的g09替換為g16。此腳本默認為在B3LYP-D3(BJ)/def2-SVP級別下對結構進行優化,然后在B3LYP-D3(BJ)/def2-TZVP級別下計算單點能同時產生分子表面靜電勢數據,然后調用formchk把chk文件轉化為fchk,最后調用Multiwfn產生RESP電荷。量化計算級別在這個腳本里都體現得很清楚,大家可以根據需要對關鍵詞進行修改。
此腳本默認的計算級別算出的RESP電荷的質量已經很好。幾何優化部分如果算不動或者想顯著降低耗時的話,可以把原有的關鍵詞改為# PM6D3 opt,即使用PM6-D3半經驗方法在真空中做優化,但這個級別只適合普通有機體系而且沒有體系局部顯離子特征的情況。單點任務部分算不動的話,可以把def2-TZVP降到6-311G**,如果是陰離子建議改為6-311+G**。有時候對柔性大體系的優化過程可能難收斂,可以嘗試把opt改為諸如opt(loose,gdiis)之類使得達到收斂更容易的關鍵詞,詳見《量子化學計算中幫助幾何優化收斂的常用方法》(http://www.shanxitv.org/164)。
注意Gaussian09 D.01版本之前不支持DFT-D3校正,因此用老版本者應當把腳本里em=GD3BJ關鍵詞取消,相關信息看《DFT-D色散校正的使用》(http://www.shanxitv.org/210)。
2 使用腳本
將Multiwfn目錄下的examples\RESP\目錄下的RESP.sh拷到當前目錄,使用chmod +x ./RESP.sh給它增加可執行權限。把含有分子結構信息的文件也拷到當前目錄,格式必須是Multiwfn認識的(可以是xyz/mol/mol2/pdb/fch/wfn/molden/gjf等等)。然后運行下面的命令:
./RESP.sh [文件名] [凈電荷] [自旋多重度] [溶劑名]
比如./RESP.sh H2O.xyz 0 1 ethanol命令將在IEFPCM隱式溶劑模型表現的乙醇環境中計算中性單重態水分子的RESP電荷。如果凈電荷和自旋多重度不寫,則分別默認為0和1。如果溶劑名不寫,默認為water,如果寫gas,則在真空下計算。支持的溶劑名見http://www.shanxitv.org/g09/k_scrf.htm的末尾。
計算完畢后,當前目錄下就有了和輸入文件同名的.chg文件。這是個文本文件,前四列是元素名和優化后的坐標,最后一列就是RESP電荷。以下是運算過程的輸出例子:
[root@192 other2]# ./RESP.sh H2O.xyz
Net charge was not defined. Default to 0
Spin multiplicity was not defined. Default to 1
Running optimization task via Gaussian...
Done!
Running single point task via Gaussian...
Done!
Running formchk...
Running Multiwfn...
Finished! The optimized atomic coordinates with RESP charges (the last column) have been exported to H2O.chg in current folder
如果調用Gaussian運行失敗,腳本就會自動停止,屆時大家請仔細檢查當前目錄下的gau.out文件判斷出錯原因。
計算的時候用的CPU核數取決于Default.Route文件里的-P-設置,這點在《Gaussian的安裝方法及運行時的相關問題》(http://www.shanxitv.org/439)里說了。
本文的做法只是以最一般的方式計算RESP電荷,如果需要更靈活的擬合方式,比如設置電荷等價性、設置某個片段總電荷為特定值等等,需要在Multiwfn的RESP模塊的界面里進行設置,在《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》(http://www.shanxitv.org/441)文中有充分、詳細的介紹。
如果你的輸入的結構文件里的結構就已經足夠好,不想讓腳本自動再做優化浪費時間,可以用examples\RESP\目錄下的RESP_noopt.sh代替前述的RESP.sh,二者用法完全一樣,只不過前者不做優化步驟。
如果你的體系里含有18號及以后的元素,由于Gaussian沒有內置擬合靜電勢用的原子半徑,因此無法通過以上腳本自動得到RESP電荷,而必須手動做Gaussian計算得到fch文件后自行利用Multiwfn按照http://www.shanxitv.org/441說的操作來計算RESP電荷(屆時也會提示缺半徑,但按照屏幕上的提示,用Multiwfn自動推薦的半徑即可)。
如果以本文的方式計算了RESP電荷,請別忘了引用Multiwfn原文。引用方式在《Multiwfn FAQ》(http://www.shanxitv.org/452)里寫明了。