RESP2原子電荷的思想以及在Multiwfn中的計算

RESP2原子電荷的思想以及在Multiwfn中的計算

First release: 2020-Feb-1 Last update: 2022-Aug-6

1 前言

RESP是幾乎最適合用于分子動力學模擬的原子電荷計算方法，在《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）中筆者對RESP電荷有非常詳細的介紹。RESP電荷的具體數值明顯依賴于計算電子波函數用的量子化學方法和基組，以及對溶劑效應的考慮。有不少人關心到底什么計算級別最適合計算RESP電荷，筆者在上文中已經做了一些說明，在本文筆者將對這個問題做更具體的討論，并且介紹RESP電荷的變體RESP2電荷的思想，以及如何用Multiwfn計算RESP2電荷。目前已經有人使用Multiwfn計算了RESP2電荷并發表了文章，比如Nature Energy (2023) DOI: 10.1038/s41560-023-01275-y、Chem. Phys. Lett., 754, 137707 (2020)。

使用本文介紹的方法用Multiwfn計算RESP2電荷時請別忘了恰當引用Multiwfn，引用方式在《Multiwfn FAQ》（http://www.shanxitv.org/452）一開始就說明了。

2 關于溶劑對溶質電荷分布極化效應的描述

在談RESP2電荷之前，先說一些相關知識。

在溶劑環境下，溶質的電荷分布會被溶劑環境所極化。量子化學計算中溶劑環境通常通過PCM等隱式溶劑模型表現，溶劑對溶質極化的作用關鍵在于溶劑的介電常數。缺乏溶劑模型相關知識的話看《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）。在《談談如何衡量分子的極性》（http://www.shanxitv.org/518）中筆者提到，偶極矩是衡量小分子極性的常用指標之一。下圖展示了甲醛在不同介電常數的溶劑下的偶極矩，計算在B3LYP/aug-cc-pVTZ下進行，溶劑模型用Gaussian默認的IEFPCM，結構是在B3LYP/def-TZVP級別下在真空中優化的。圖中對應真空的情況（介電常數=1）是去掉溶劑模型后計算的結果。

可見，溶劑的介電常數越大，溶質在溶劑中就被極化得越厲害、極性比真空下就大得越多（但極化程度和介電常數不是線性關系，從上圖也能看出來）。水是最常見的溶劑，常溫下介電常數是78.3。

對于目前最常用的固定電荷力場（相對于可極化力場而言），在凝聚相下模擬時，原子電荷分布應當能等效地體現出實際環境產生的極化效果。例如常用的SPC/E水模型對應的偶極矩2.35 Debye就明顯高于水的氣相實驗偶極矩1.85 Debye。

以往的AMBER力場中大多數版本，包括2019年提出來的AMBER19SB，都是在HF/6-31G*級別下算的RESP電荷。稍有量子化學常識的人都知道HF/6-31G*非常垃圾、無法忍受。之所以HF/6-31G*下的RESP電荷還被沿用至今，是因為HF完全不能考慮電子庫侖相關，這使得波函數質量很爛、體系的極性整體被一定程度高估，而這誤差恰巧能一定程度表現凝聚相中溶劑對溶質的極化作用。另一方面，主要描述蛋白質和核酸的AMBER力場大多數版本都是結合著HF/6-31G*的RESP電荷擬合的力場參數，參數與電荷間有一定耦合（尤其是氨基酸的參數，稍微變變就可能導致長時間跑出來的蛋白質構象改變不少，開發者不敢輕易動），因此搞力場的人即便明知道HF/6-31G*組合如今看起來很渣，也還繼續容忍著這種陋習。

顯然，如果用像樣的DFT泛函結合隱式溶劑模型來體現溶劑的極化作用，比起靠HF/6-31G*氣相計算這種“以錯誤的方式得到湊合能用的結果”要優雅、理想得多，也顯然是未來發展新力場時勢必要用的做法。

在J. Chem. Inf. Model., 60, 249 (2020)一文中，作者專門考慮了不同方式計算的一批很小的有機分子的偶極矩，如下所示。其中橫坐標是氣相實驗偶極矩，縱坐標是不同方式計算出的偶極矩。（注：DFT結合圖中這么大的帶彌散函數的3-zeta基組計算的氣相偶極矩可以和實驗值對得很好，所以與氣相實驗值的偏差純粹體現的是隱式溶劑模型的影響）

從上圖可以得到以下結論
(1)HF/6-31G*傾向于高估體系極性（擬合的直線斜率大于1），但數據點分布比較散，因此高估程度不一，甚至對于個別體系反倒還低估了偶極矩（散點在對角線下方）。
(2)如果HF計算時改用較大基組，則線性相關性比起用6-31G*好得多，即數據質量更好了，而對極性的平均高估程度倒是和6-31G*相仿佛。
(3)常用的B3LYP泛函結合較大基組并考慮溶劑效應時，哪怕只用介電常數比較低的苯作為溶劑（介電常數為2.2），體系的極性也比HF在真空下計算的更高。若改用極性很大的水作為溶劑，則極化程度還能再提高一倍。

從上述對比可以看出，使用HF/6-31G*算RESP電荷的做法絕對應當被摒棄，因為其數據質量太低，溶劑的等效極化效應表現的可靠程度太差（數據點太分散，時高時低）。再者，實際的溶劑只會對直接暴露在溶劑環境的位于溶質表面的原子的電荷分布產生明顯的極化，并不會顯著影響溶質內部區域的電荷分布，而在HF下溶質所有區域的極性都傾向于被高估，這和實際不符。另外，HF/6-31G*下計算RESP電荷完全體現不出不同溶劑對溶質極化程度的差異，這又是基本原理上的一個硬傷。

由于從上圖可見，DFT結合隱式溶劑模型計算出的偶極矩分布沒有那么分散，因此以這種組合在RESP電荷計算時考慮溶劑效應在理論上是非常理想的。然而要注意的是，對于RESP電荷計算目的，并非體系實際處于什么溶劑環境，就應當選擇什么溶劑用于溶劑模型。有以下幾點值得關注：
(1)基于HF/6-31G*的RESP電荷下用AMBER、GAFF力場得到的水環境下的模擬結果通常還不錯。而從上圖來看，如果DFT計算時用隱式溶劑模型表現水，那么溶質的極性就顯得有點高過頭了，整體遠高于HF/6-31G*的情況，因此從直覺來看此時的結果有可能更差。
(2)AMBER03力場也是用的RESP電荷，且計算RESP電荷時是在看起來比較理想的B3LYP/cc-pVTZ結合IEFPCM溶劑模型下進行的，但是設的介電常數是4。之所以設得這么小是避免過度極化，否則最終會高估靜電相互作用，導致出現氫鍵過強等問題。而且從前面的圖來看，哪怕用苯這么低極性的溶劑，產生的極化效應也并不小，比HF/6-31G*高估極性的程度甚至整體還高。這體現出刻意用較小的介電常數結合隱式溶劑模型來計算水環境中模擬用的RESP電荷有一定合理性。
(3)在J. Chem. Inf. Model., 60, 249 (2020)文中提到，對于高極性溶劑環境下的模擬，用于計算原子電荷的體系波函數中，體系的被極化程度理應低于實際溶劑中的情況，這用于體現體系電荷分布從孤立狀態被極化到溶劑環境中的狀態過程中的能量消耗。

3 IPolQ與RESP2電荷

J. Phys. Chem. B, 117, 2328 (2013)中提出了IPolQ原子電荷，它是將氣相和顯式水下算的RESP電荷取平均。后來J Comput Aided Mol Des, 28, 277 (2014)中提出的IPolQ-mod將IPolQ的計算簡化，變成了將氣相和隱式水下算的RESP電荷取平均，用的是MP2/aug-cc-pVDZ級別（此級別對電荷分布的描述和上圖中B3LYP/aug-cc-pV(T+d)Z其實差不多）。測試發現用IPolQ-mod電荷算水合自由能的結果和基于HF/6-31G*下RESP電荷計算的結果差不多。因此，IPolQ-mod這種在靠譜計算級別下，將氣相和隱式溶劑下算的RESP電荷取平均的做法是值得提倡的方法，因為雖然結果看起來并不比用基于HF/6-31G*的RESP電荷算的有顯著優勢，但做法明顯嚴格、合理得多得多，還能有效避免因為用了垃圾的HF/6-31G*被他人批評。

在2019年的https://doi.org/10.26434/chemrxiv.10072799.v1一文中，作者將IPolQ-mod的思想進一步廣義化，提出了RESP2電荷（PS：這名字起得真夠吸引眼球的）。它的定義是：q_RESP2 = (1-δ)*q_gas + δ*q_aqu。其中q_gas和q_aqu分別是氣相和水模型下計算的RESP電荷，δ是可調參數。諸如δ=0.5的時候可以叫RESP2(0.5)。實際上RESP2(0.5)就相當于IPolQ-mod電荷。

2020-Apr-29補充：RESP2的文章已發表于Communications Chemistry, 3, 44 (2020)，地址：https://www.nature.com/articles/s42004-020-0291-4。

注：RESP2的文章里用的計算級別是PW6B95/aug-cc-pV(D+d)Z，在我來看用這個并沒有什么特別的道理。PW6B95是個非主流泛函，對于RESP電荷計算來說不會比常用的B3LYP有任何優勢，絕對沒必要刻意用。而用aug-cc-pV(D+d)Z這種基組其實還不如用def2-TZVP，因為對中性體系，與其加彌散函數還不如先從2-zeta升到3-zeta，而且def2-TZVP的f極化函數也對于更好地描述電荷分布有不可忽視的幫助（尤其是雜原子）。如果是陰離子，個人建議用def2-TZVP加彌散的版本ma-TZVP，詳見http://www.shanxitv.org/509。

RESP2這篇文章基于非主流的SMIRNOFF力場，對一批有機體系模擬了密度、蒸發焓、介電常數、水合自由能，根據文中的數據可以發現以下情況（以下RESP代表HF/6-31G*氣相算的RESP電荷，RESP2代表PW6B95/aug-cc-pV(D+d)Z算的RESP2電荷）：
(1)不重新優化LJ參數時，基于RESP2(0.6)電荷算的介電常數比RESP好很多，其它性質則半斤八兩。如果加大δ，可以令介電常數和水合自由能的模擬精度比RESP更有優勢（δ=1時最理想）。如果減小δ，可以令密度的精度更有優勢（δ=0.2時比較合適）
(2)如果重新優化LJ參數，δ在0.5~0.7時整體誤差（蒸發焓、密度、介電常數、水和自由能的平均誤差）最低，盡管不同的屬性對于δ的傾向性有所不同。重新優化LJ參數時，RESP2(0.6)比RESP除了算蒸發焓以外都更好，說明用RESP2電荷的話比基于RESP電荷明顯有更大的余地去進一步優化LJ參數來改進結果，這點從下圖可以非常清楚地看出。圖中中紅線是基于RESP電荷算的各種性質的平均誤差，藍線曲線是對RESP2電荷掃描δ的情況，每次計算時都重新優化了LJ參數。

RESP2這篇文章進一步體現，在像樣的計算級別下，將氣相和溶劑模型下的RESP電荷相混合，比起用垃圾的HF/6-31G*在氣相下算RESP電荷更妥，后者儼然已經過時，是時候拋棄它了。

此文中發現RESP2計算性質整體誤差最小時δ是0.6，和IPolQ電荷對應的δ=0.5非常接近，而且RESP2(0.5)和RESP2(0.6)的誤差差不多，因此δ=0.5應當認為是一個比較普適的參數。RESP2文章里沒有測試Amber03那種算RESP電荷的做法，即基于介電常數為4的溶劑用像樣級別算RESP電荷，這有點可惜。

RESP2原文里的RESP2電荷的定義在我來看有個不妥之處是不管體系處在什么外環境下，在溶劑模型下計算時都是用的水，這明顯不科學。一個溶質在低極性的溶劑比如氯仿中以及在高極性的溶劑比如水當中的電荷分布差異是很大的，而隱式溶劑模型明明能考慮不同溶劑何故不如實地考慮？如果此文里諸如計算呋喃的密度時就用呋喃作為溶劑來計算RESP2公式中溶劑下的RESP電荷，估計會使文中的除了水合自由能以外的性質都得到不小改進。

4 到底該怎么計算適合AMBER/GAFF力場模擬的電荷？

上面說了一大堆，到底怎樣計算最佳的用于AMBER/GAFF力場模擬的原子電荷？我個人目前推薦這樣計算RESP2(0.5)電荷：
q_RESP2(0.5) = 0.5*q_gas + 0.5*q_solv
其中q_gas是氣相下的RESP電荷，q_solv是用PCM模型表現實際溶劑時的RESP電荷（注：G09、G16默認的IEFPCM是PCM最佳的實現形式）。幾何優化使用B3LYP-D3(BJ)搭配6-311G**（或更好一點點的def-TZVP），優化時在真空下即可，但如果是帶電體系或局部帶顯著電荷的體系則應當在溶劑模型表現實際溶劑環境下進行優化。之后產生用于計算q_gas和q_solv的波函數文件時的單點任務使用B3LYP-D3(BJ)/def2-TZVP，如果是陰離子則基組改為ma-TZVP。只要有個基本像樣的服務器，這樣的計算級別用Gaussian+Multiwfn算上百個原子的RESP2電荷完全沒有壓力。但如果體系更大，或者要算一大批體系，或者機子實在爛，那么優化可以用6-31G**，單點可以用def-TZVP，再低就沒法忍受了。

以上述做法計算電荷，從任何角度都說得通。但可能有人會嘰嘰歪歪說，對之前AMBER力場的大多數版本，以及專門用于有機小分子的GAFF力場，官方都是HF/6-31G*下算的RESP電荷，你這樣拿RESP2(0.5)電荷搭配會不會有不兼容性？這完全不必顧慮。前文已經說了，原理較嚴格的RESP2(0.5)的電荷實際表現不會比HF/6-31G*的RESP電荷差，很多情況還更好，而且前文里IPolQ-mod文章就是基于GAFF力場做的對比，SMIRNOFF力場的LJ參數也都是來自AMBER/GAFF的。

5 使用Multiwfn計算RESP2電荷實例

下面就通過實例完整演示一下使用Gaussian結合Multiwfn計算RESP2(0.5)電荷的過程，以計算乙醇環境中的甲醛為例。Gaussian使用16 A.03版，Multiwfn使用2020-Jun-30更新的3.7(dev)版。如果還沒讀過《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》（http://www.shanxitv.org/441）一文的話請先閱讀，里面有很多關鍵性知識應當了解。對Multiwfn不了解者強烈建議閱讀《Multiwfn FAQ》（http://www.shanxitv.org/452）。下面用到的文件都可以在http://www.shanxitv.org/attach/531/file.rar里找到。

將以下內容存為Gaussian輸入文件然后運行。初始結構是隨便畫的。此任務分三步，會先在B3LYP-D3(BJ)/def-TZVP下進行優化，然后基于優化過的結構利用B3LYP-D3(BJ)/def2-TZVP下在氣相和IEFPCM模型表現的乙醇環境下分別算單點，最終在C:\下得到SP_gas.chk和SP_solv.chk。讀者應擇情對內容進行恰當修改（諸如G09 D.01之前不支持DFT-D3校正、G09 C.01之前不支持%oldchk。當前體系用DFT-D3校正其實完全沒必要，加這個只是為了增加這個模板文件的普適性而已）

%chk=C:\opt.chk
# B3LYP/TZVP em=GD3BJ opt

niconiconi

0 1
C                  0.00000000    0.00000000    0.52887991
H                  0.00000000    0.93775230    1.12379107
O                  0.00000000    0.00000000   -0.67757652
H                  0.00000000   -0.93775230    1.12379107

--link1--
%oldchk=C:\opt.chk
%chk=C:\SP_gas.chk
# B3LYP/def2TZVP em=GD3BJ geom=allcheck

--link1--
%oldchk=C:\opt.chk
%chk=C:\SP_solv.chk
# B3LYP/def2TZVP em=GD3BJ scrf=solvent=ethanol geom=allcheck
<--空行
<--空行

算完之后，把SP_gas.chk和SP_solv.chk都轉換為fch文件，不知道怎么做的話看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。

先計算氣相下的RESP電荷。啟動Multiwfn后，輸入
SP_gas.fch //填實際的路徑
7 //原子電荷計算
18 //RESP電荷
1 //開始標準兩步式RESP電荷擬合
結果為
Center Charge
1(C ) 0.4195430529
2(H ) -0.0050085243
3(O ) -0.4095260044
4(H ) -0.0050085243

然后用同樣的方法，以SP_solv.fch作為輸入文件，得到乙醇下的RESP電荷：
Center       Charge
1(C )   0.4625344852
2(H )   0.0093031373
3(O ) -0.4811407598
4(H )   0.0093031373

將以上二者用諸如Excel取平均，就得到了RESP2(0.5)電荷，即
0.441038769
0.002147307
-0.445333382
0.002147307

6 通過腳本方便地調用Multiwfn計算RESP2電荷

將兩個文件分別計算再手動取平均是非常簡單的事，但有的人就是嫌麻煩，還有人誤以為Multiwfn計算電荷的過程沒法納入批處理什么的。為此，筆者提供了一個一鍵完成RESP2電荷計算的Linux腳本calcRESP.sh。此腳本在目前最新版本的Multiwfn文件包的examples\RESP2目錄下可以找到。使用前需要先按照手冊2.1.2節將Multiwfn正確安裝，使得程序可以通過輸入Multiwfn命令直接啟動。這個腳本三種用法示例：

對SP_gas.fch計算RESP電荷：./calcRESP.sh SP_gas.fch
基于SP_gas.fch和SP_solv.fch計算RESP2(0.5)電荷：./calcRESP.sh SP_gas.fch SP_solv.fch
基于SP_gas.fch和SP_solv.fch計算RESP2(0.7)電荷：./calcRESP.sh SP_gas.fch SP_solv.fch 0.7

運行過后，RESP電荷會輸出到當前目錄下的與輸入文件同名但是后綴為chg的文件中，而RESP2電荷會輸出為RESP2.chg。chg是Multiwfn的私有的記錄原子電荷的格式，其中2~4列是來自fch文件里的坐標，最后一列是算出來的電荷值。

此腳本的輸入文件不限于fch，用其它的Multiwfn可以讀取基函數信息的格式如.molden等都可以，見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。如果你已經將Multiwfn的settings.ini里的formchkpath設成了機子里實際的formchk路徑，也可以直接用chk作為輸入文件。

7 基于結構文件傻瓜式一鍵調用Gaussian和Multiwfn計算RESP2電荷

之前筆者針對RESP電荷寫過《計算RESP原子電荷的超級懶人腳本（一行命令就算出結果）》（http://www.shanxitv.org/476）一文，對其中介紹的腳本RESP.sh只需要提供含有結構信息的文件，腳本就會一條龍地調用Gaussian和Multiwfn完成RESP電荷計算，完全不懂Gaussian的人都能一鍵計算出RESP電荷，沒有讀過此文者請先閱讀。本文提供一個用法類似的腳本，但專門用于計算RESP2電荷，這個腳本就是Multiwfn文件包里的examples\RESP\RESP2.sh。用法示例如下：
對中性單重態分子計算用于水溶劑MD模擬的RESP2(0.5)電荷：./RESP2.sh H2O.pdb
對中性三重態分子計算用于水溶劑MD模擬的RESP2(0.5)電荷: ./RESP2.sh foo.xyz 0 3
對陰離子單重態分子計算用于乙醇溶劑MD模擬的RESP2(0.5)電荷 ./RESP2.sh nozomi.mol -1 1 ethanol
可見，電荷和自旋多重度分別默認為0和1，溶劑默認是水。δ用的是0.5，想改的話直接修改腳本即可。算完之后，當前目錄下gas.chg和solv.chg分別是氣相和溶劑下的chg文件，與輸入文件同名但后綴為chg的是RESP2(0.5)電荷的文件。

此腳本用的計算級別只要打開.sh文件一看便知，對幾何優化用的是B3LYP-D3(BJ)/def2-SVP，對單點用的是B3LYP-D3(BJ)/def2-TZVP。調用的是Gaussian 09，想改成調用Gaussian 16就把腳本中的g09改為g16。

如果你的輸入的結構文件里的結構就已經足夠好，不想讓腳本自動再做優化浪費時間，可以用examples\RESP\目錄下的RESP2_noopt.sh代替前述的RESP2.sh，二者用法完全一樣，只不過前者不做優化步驟。

如果你沒買Gaussian，也可以用免費的ORCA量子化學程序結合Multiwfn算RESP2原子電荷，筆者也提供了相應的傻瓜式腳本，見《ORCA結合Multiwfn計算RESP、RESP2和1.2*CM5原子電荷的懶人腳本》（http://www.shanxitv.org/637）。