通過獨立梯度模型(IGM)考察分子內和分子間的弱相互作用

后記1：筆者后來又寫了一篇文章介紹怎么對周期性體系做IGM分析，見《使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用》（http://www.shanxitv.org/588）。

后記2：2022年我提出了叫做IGMH的改進版IGM方法，圖像效果顯著好于IGM，還可以完全避免IGM圖中誤導性的著色出現，而且物理意義明顯比IGM更強。請讀者們務必閱讀IGMH的全面介紹文章《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621），其中還有十分詳細的用Multiwfn做IGMH分析的教程。由于IGMH的計算基于真實密度而非像IGM那樣基于粗糙的準分子密度計算，因此IGMH的耗時高于IGM。但對于幾百個原子體系，在一般計算條件下在Multiwfn程序里做IGMH分析并無壓力。因此但凡不是特別巨大的體系，一律強烈建議用IGMH代替IGM！！！

通過獨立梯度模型(IGM)考察分子間弱相互作用

文/Sobereva @北京科音

First release: 2018-Apr-5  Last update: 2022-Feb-14

0 前言

早在2010年，筆者寫了一篇《使用Multiwfn圖形化研究弱相互作用》（http://www.shanxitv.org/68），文中詳細介紹了楊偉濤等人通過約化密度梯度(RDG)圖形化考察弱相互作用的方法以及在Multiwfn中的實現。這種分析方法在文獻中也常被叫做NCI（Noncovalent interaction）方法。此文貼出來后，不斷有大量使用Multiwfn做RDG分析的文章涌現，而RDG分析方法也被很多研究者進一步發展，例如：
(1)IRI（相互作用區域指示函數）方法，由我于2021年提出。它可以通過一個函數同時完美展現體系中各種類型的相互作用，完全替代了RDG的用處，強烈建議一讀：《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）
(2)aRDG方法，用分析動力學過程中的弱相互作用，見《使用Multiwfn研究分子動力學中的弱相互作用》（http://www.shanxitv.org/186）
(3)將RDG與ELF聯用，同時考察弱相互作用和共價相互作用，例子見《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）。但有了IRI就沒必要用這個了
(4)對RDG等值面內部空間的實空間函數進行積分，從而能夠定量討論相互作用強度。見Multiwfn手冊4.200.14節的例子

2017年，在Phys. Chem. Chem. Phys., 19, 17928中新提出一種叫做獨立梯度模型（Independent Gradient Model, IGM）的方法。這個方法受到RDG方法的很大啟發，但是底層思想不同。和RDG一樣，IGM也是重在展現弱相互作用區域及其特征，但是有很多優勢，簡要列舉如下，后文會詳細說：
(a)給出的信息會明確劃分成片段間和片段內兩套數據，因此想考察分子間相互作用時，不會被分子內相互作用所干擾。（實際上，用RDG方法時，如果結合Multiwfn手冊4.13.4節示例的格點數據屏蔽功能，其實也能達到相同目的）
(b)計算時只依賴于原子坐標，而且計算很快。因為此方法計算時只需要自由原子密度，并不需要體系的波函數信息，而自由原子密度是內置于IGM分析程序當中的。（實際上，RDG也可以基于準分子(promolecular)密度來近似地快速計算）
(c)IGM方法給出的等值面圖比RDG更為飽滿一些，不像RDG圖在格點稀疏的時候容易出現難看的鋸齒和窟窿
(d)可以考察每一對原子間的相互作用程度
(e)可以定量給出每個原子、每對原子對片段間相互作用的影響程度，而且可以利用VMD通過著色直觀地展現
(f)在分子間相互作用區域，IGM方法里定義的δg函數的大小與相互作用強度有正相關性

本文的目的一方面是介紹IGM方法，一方面通過實際體系示例怎么通過Multiwfn做IGM分析。閱讀本文前，強烈建議仔細閱讀前述博文，至少弄清楚RDG分析的基本原理以及在Multiwfn中的操作，否則在理解本文內容時可能遇到困難。

Multiwfn可以免費在http://www.shanxitv.org/multiwfn下載，讀者務必使用目前官網上的最新版本。如果對Multiwfn不了解，看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。除本文提到的方法外，Multiwfn還支持一大票弱相互作用分析方法，可謂是弱相互作用分析的武器庫，看《Multiwfn支持的弱相互作用的分析方法概覽》（http://www.shanxitv.org/252）。

使用Multiwfn做IGM分析發表文章時，除了必須要引用Multiwfn程序在啟動時提示的Multiwfn的原文外，還應當引用以下兩篇文章：
(1) Phys. Chem. Chem. Phys., 19, 17928 (2017)：IGM方法的原文
(2) J. Comput. Chem., 43, 539 (2022) DOI: 10.1002/jcc.26812：里面對IGM方法、我對IGM的擴展以及在Multiwfn中的具體實現都有充分的介紹

此外，非常推薦閱讀《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）介紹的筆者的綜述文章，其中對IGM在內的各種弱相互作用可視化分析方法有非常全面系統的介紹并給了大量應用例子，也推薦在引用上述文章時一起引用。

1 IGM方法的原理

IGM的原文Phys. Chem. Chem. Phys., 19, 17928 (2017)里有此方法的介紹，后來我在IGMH的原文J. Comput. Chem. (2022) DOI: 10.1002/jcc.26812里對此方法有更全面、系統的介紹，這里主要解釋一下IGM最核心的思想。我們先看一個很簡單的情況，氫分子。在氫分子的軸線上，H1和H2的電子密度分布如下所示，這里用的原子的電子密度是原子在孤立狀態下計算的

從上圖可看出，在兩個原子間，兩個原子的梯度方向正好是相反的，這個區域H1的電子密度的梯度為負（即往左邊密度增大），而H2的電子密度的梯度為正（即往右邊密度增大）。因此，在原子間，兩個原子的密度梯度是相互抵消的，而且在兩原子正中央正好精確抵消，導致準分子密度（原子自由狀態密度疊加構建的密度）的梯度為0，此處也正是與準分子密度相對應的鍵臨界點(BCP)的位置。

一般方式計算的準分子密度的梯度等于各個原子密度梯度的加和。而IGM型密度梯度，則是把每個原子的梯度先取絕對值再加和。二者求差，就得到了所謂的δg函數。對上面的氫分子體系，圖示如下

圖中黑線是一般方式計算的準分子密度梯度，在原子核和原子連線中點都為0。綠線是IGM型密度梯度，由于計算時對兩個原子密度梯度都取了絕對值，所以不會相互抵消，在原子之間其數值也很大，它也是普通方式計算的密度梯度的上限。兩種密度梯度的差值，即δg，對應圖中深藍色曲線，可見利用δg這個函數可以明確地把原子間相互作用區域展現出來。而且，原子間相互作用越強，相互作用區域的δg就會越大。這是因為如果原子間相互作用較弱，那么優化后的結構中兩個原子距離就會較遠，它們相交疊的區域將是各自密度梯度已經比較小的部分，此時IGM型密度梯度和普通密度梯度相差也就沒那么大了。

之所以當前這個方法叫做IGM（獨立梯度模型），是因為常規方式計算準分子密度梯度時，相當于原子密度梯度之間發生了“干涉”，影響了“相位”。而在計算IGM型密度梯度時，純粹是原子的密度梯度大小的簡單疊加，不存在干涉現象，因此這些原子的密度梯度就彼此“獨立”了。

δg是個三維實空間函數，有不同的考察方式，如果圖形化展現的話，一般是像RDG那樣繪制等值面圖。比如對上面的氫分子體系，如果繪制δg=0.2的等值面圖，就肯定會看到兩個氫之間出現了等值面，因為上圖中對應Y=0.2的虛線和δg曲線在黑色箭頭標注的位置出現了相交。

δg對應于Multiwfn中第22號實空間函數。用主功能4對其繪制填色平面圖的話，就可以看到下面的效果。下圖是GC堿基二聚體平面上的情況

 

由圖可見，所有成鍵的原子間δg都較大，最大值處都超過了圖中色彩刻度上限0.3因此顯示為白色。而弱相互作用區域，比如N-H...N的氫鍵區域，只有一塊深藍色，說明這樣的地方δg相對較小，作用較弱，但終究比沒有相互作用的地方δg值要大。IGM原文里通過測試還發現，二聚體的結合能，和兩個單體間對應弱相互作用的BCP位置上的δg值有很好的正相關性。因此δg是個很有用的函數。

δg可以劃分為用于展現片段內相互作用的δg_intra和展現片段間相互作用的δg_inter。在說這個之前，這里先把δg的比較嚴謹的數學定義給出來。下式中i是原子序號，▽ρ是梯度矢量，abs(▽ρ)代表對里面▽ρ矢量的每個分量都取絕對值（還是保持矢量形式），| |代表對矢量取模。

 

IGM原文中計算δg_intra和δg_inter只考慮到兩個片段的情況，但筆者將之推廣到了多個片段的情況，可表達為：

 

式中A是片段編號。計算時不要求所有片段的并集必須對應整個體系，但不對應時，在計算δg的時候必須只考慮這些片段涉及的原子。只要想明白了δg是怎么定義的，就很容易理解δg_inter為什么是這樣定義的。

由上式可見，δg_intra+δg_inter=δg。δg衡量的是體系中所有原子間相互作用，稍加變通就成了只體現片段間相互作用的δg_inter，再把這部分從δg中扣除就成了描述片段內相互作用的δg_intra。IGM的這個思想相當不錯。

對δg_intra和δg_inter分別繪制等值面圖，就可以清楚看到分子內和分子間存在相互作用的區域。還可以像RDG分析時那樣，把sign(lambda2)rho函數以不同顏色投影到這兩個函數的等值面上，從而能夠清楚地判斷相互作用區域是吸引還是互斥作用，以及強度如何。另外，觀看δg_intra等值面的時候還可以像RDG分析那樣，把電子密度較大的區域的δg_intra值設為0，這樣等值面就只體現出分子內弱相互作用區域了，而化學鍵區域就都被屏蔽掉了。

還可以定義原子對δg指數，用于要考察相互作用的兩個片段間的每對原子的相互作用大小。這個指標在IGM原文里沒有提出，是我把IGM實現進Multiwfn的時候順帶想出來的。計算某兩個原子間的這個指數，就是把這兩個原子各自作為一個片段來計算它們之間的δg_inter函數，然后再對這個函數在全空間進行積分，數值越大顯然說明兩個原子間作用越強。然后，可以再計算原子δg指數，就是把這個原子與另一個片段中所有原子的原子對δg指數進行加和。之后，可以在VMD里將每個原子按照其δg指數進行著色，哪些原子對片段間相互作用起主要貢獻就一目了然了。

下面我們用Multiwfn通過IGM方法研究一系列體系，使讀者了解IGM分析在Multiwfn中的用法以及IGM方法的價值。IGM方法靈活程度很高，讀者務必舉一反三。本文用到的所有文件都可以在此下載：http://www.shanxitv.org/attach/407/file.rar。本文使用的VMD是1.9.3版，可在http://www.ks.uiuc.edu/Research/vmd/免費下載。由于本文討論的都是基于promolecular密度的IGM分析，所以輸入文件只需要包含結構信息就行了，常見的記錄分子結構的.xyz、.pdb、.mol等格式，以及記錄波函數信息的.fch、.wfn、.molden等格式都可以用。Multiwfn支持的輸入文件中哪些能提供結構信息見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。

注：下文的例子的數據和圖像是基于較早的Multiwfn版本得到的，后來Multiwfn在IGM分析方面又做了一些改進，故目前版本的結果與下文例子會有些許不同，以最新版本為準。

要注意，和RDG方法一樣，所有用IGM分析的體系的幾何結構一定要經過在恰當的級別下優化，否則可能會給出錯誤的分析結果。上面提供的文件包里的結構都是經過優化的。用高精度晶體衍射得到的結構也可以，但是如果考察的相互作用牽扯到氫原子，則應當固定重原子而把氫原子位置進行優化（因為氫原子位置是沒法準確測定的），詳見《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）。

2 實例：繪制甲酸二聚體的δg圖像

此例演示一下怎么將IGM方法中定義的δg函數繪制成平面圖和等值面圖。本例使用的是文件包里的aceticacid_dimer.xyz，這是甲酸二聚體，所有原子都在Z=0的XY平面上。

首先繪制體系平面上的填色圖。啟動Multiwfn，依次輸入
aceticacid_dimer.xyz
4   //繪制平面圖
22   //δg函數
1   //填色圖
[回車]  //用默認的格點設定
1  //XY平面
0  //Z=0
從屏幕上直接顯示的圖像看不出什么特征來，這是因為默認的色彩刻度對當前體系并不合適。屏幕輸出的信息中可以找到這樣的提示，表明當前平面上函數最小值為0，最大值為0.54956：
The minimum of data:  0.000000000000000E+000
The maximum of data:  0.549560853101836
因此，為了讓不同區域的函數值能夠被不同顏色區分開，我們需要調整色彩刻度，讓下限和上限大體與函數值最小值和最大值分別相仿佛。于是關閉圖像，輸入
1  //修改色彩刻度
0,0.5  //設定色彩刻度范圍為0~0.5
4  //允許顯示原子標簽
1  //使用紅色標簽
此時再輸入-1，重新顯示出的圖像如下

此圖體把體系中的共價鍵和單體之間的兩條氫鍵都展現出來了。由于共價鍵明顯比氫鍵強，所以相應區域的δg數值明顯比氫鍵區域的要大得多。δg函數在某種程度上與DORI有相同的價值，即可以體現出體系中任何類型的相互作用，但是δg還可以把強度差異體現出來，而DORI不反映這一點。

接下來我們繪制一下δg函數的等值面圖。選-5退到主菜單，然后依次輸入
5  //計算格點數據
22   //δg函數
2  //中等質量格點
-1  //觀看等值面
在不同的等值面數值(isovalue)下看到的等值面如下所示

 

可見在isovalue=0.15時，只有化學鍵區域被展現了出來。而當isovalue減小到0.05時，氫鍵作用區域才出現。讀者可以嘗試把isovalue逐漸由大到小調節，看等值面怎么變化，看到的現象是越弱的相互作用對應的等值面出現得越晚。因此，當我們只對化學鍵感興趣而不想討論弱相互作用的時候，只需要繪制isovalue較大時的δg等值面即可。

3 實例：用IGM考察苯酚二聚體中的相互作用

上面的例子僅僅是前戲，當前這個例子才展示Multiwfn中專門的IGM分析模塊的用法。IGM模塊在使用時需要定義片段，Multiwfn允許定義無數個片段。當定義n個片段時，n個片段內的相互作用和n個片段間的相互作用可以分別展示互不干擾，而且程序還允許對每一對片段之間的相互作用進行細節分析。所有定義的片段的并集并不要求是整個體系，比如對10聚體團簇我們可以只把其中離得近的三個分子定義成三個片段來分析它們的相互作用，而其它單體在計算時會被忽略。一個片段并不需要對應一個分子，完全可以將一個大分子劃分成兩個甚至多個片段。用IGM模塊時即便你并不打算分離片段間和片段內相互作用，也至少要定義一個片段來將你感興趣的區域的原子納入進去。

本例我們將把苯酚二聚體的兩個單體各定義為一個片段來分析片段間的相互作用。啟動Multiwfn，載入苯酚二聚體的結構文件phenoldimer.xyz，然后依次輸入
20   //圖形化分析弱相互作用
10   //IGM分析
2    //定義兩個片段
1-13   //片段1包含的原子。即第一個苯酚的原子序號范圍
14-26    //片段2包含的原子。即第二個苯酚的原子序號范圍
2    //中等質量格點（約512000個點）
然后程序開始計算格點數據，很快就算完了。

用過RDG方法分析弱相互作用的人肯定對RDG vs sign(lambda2)rho的散點圖不陌生，類似地，我們可以繪制δg或δg_intra或δg_inter vs sign(lambda2)rho的散點圖來考察體系中的相互作用。我們先繪制δg vs sign(lambda2)rho的散點圖。在后處理界面選擇-1，然后選1，就得到了相應的圖像
 

按照討論RDG散點圖的方式討論此圖，我們看到圖中sign(lambda2)rho顯著大于0的區域有一大堆點，因此可以說體系里存在位阻作用。在sign(lambda2)rho約-0.04的位置有一個峰，由于這個位置電子密度不很大，但又不很接近于0，因此可以初步判斷應該是對應苯酚二聚體之間的氫鍵。在sign(lambda2)rho很負的區域有大面積的峰，由于這些位置電子密度較大，因此可以認為是對應化學鍵區域的點。若我們在后處理界面選擇-1后選擇1，得到的將是體現片段間作用的δg_inter vs. sign(lambda2)rho散點圖；如果選擇2，得到的將是體現片段內作用的δg_intra vs. sign(lambda2)rho散點圖。我們上面展示的圖相當于這兩種圖的疊加。有興趣者也可以選擇-1后選擇4，會把δg_intra和δg_inter分別用黑色和紅色的點同時顯示在同一張圖上便于分別討論。
 

然后我們來看看等值面圖。后處理菜單選4 Show isosurface of grid data，然后我們可以選擇觀看δg_inter、δg_intra以及δg的等值面圖，前兩者如下所示

 

由圖可見，利用IGM方法，可以非常好地把片段內和片段間的相互作用分離開。而用RDG方法的時候，為了達到這個目的，需要通過Multiwfn的主功能13來試圖對RDG格點數據進行屏蔽（見Multiwfn手冊4.13.4節的實例），這麻煩得多，有時候還得來回試參數。

接下來我們利用VMD來繪制sign(lambda2)rho填色的δg_inter等值面圖。在后處理菜單中選擇3來把δg_inter、δg_intr、δg、sign(lambda2)rho的格點數據分別導出為當前目錄下的dg_inter.cub、dg_intra.cub、dg.cub和sl2r.cub。將dg_inter.cub和sl2r.cub都拷到VMD目錄下，把Multiwfn自帶的examples目錄下的VMD作圖腳本IGM_inter.vmd也拷到VMD目錄下，然后啟動VMD，在命令行窗口輸入source IGM_inter.vmd，就立馬看到以下圖像

 

這個IGM_inter.vmd腳本和Multiwfn自帶的繪制RDG填色圖的腳本RDGfill.vmd非常類似。不過由于δg_inter和RDG的函數特征還是存在差異的，所以IGM_inter.vmd腳本默認用的色彩刻度是-0.05~0.05，與RDGfill.vmd用的不同，不過對圖像的解讀還是和RDG填色圖完全一樣。上圖當中形成氫鍵的區域等值面顏色是明顯的藍色，因此應當認為確實形成了氫鍵。而在旁邊的區域還有大面積綠色扁片，體現相應弱相互作用區域電子密度很小，數值接近0，因此作用強度較弱，應當解釋為色散作用。對此體系的IGM分析結論和RDG方法是完全一樣的。不過，做過RDG分析的人會明顯感覺IGM的等值面和RDG等值面的特點還是有一定區別的，主要差異在于RDG等值面整體比較薄，而IGM等值面整體比較豐滿。豐滿帶來的好處是等值面邊緣比較平滑，也不像RDG等值面那樣往往出現斑駁難看的窟窿。用IGM的時候，使用比RDG分析時明顯更少的格點數就可以達到與之同檔次的圖像效果。

繪制δg_intra的填色圖的做法是把dg_intra.cub、sl2r.cub、examples\IGM_intra.vmd都拷到VMD目錄下，然后在VMD里執行source IGM_intra.vmd。圖片就不在這里展示了。

下面我們把兩個苯酚之間的相互作用分解成原子的貢獻來考察。在Multiwfn的IGM模塊的后處理菜單中選6 Evaluate contribution of atomic pairs and atoms to interfragment interaction，然后選High quality積分格點，算完后，當前目錄下就出現了atmdg.txt。文件中先給出這樣的信息
 Atom dg index in fragment  1
  Atom    9 :    0.468547
 Atom   10 :    0.279359
 Atom    2 :    0.237590
 Atom    7 :    0.216155
 Atom    1 :    0.182576
...略
 Atom dg index in fragment  2
  Atom   20 :    0.427777
 Atom   15 :    0.277831
 Atom   14 :    0.273219
 Atom   22 :    0.243161
...略

這就是之前提到的原子δg指數，體現了兩個苯酚中各個原子對苯酚間相互作用貢獻的程度，由大到小排序。通過這樣信息，可以很方便而且定量地考察原子對感興趣的弱相互作用的貢獻程度。

再往下是原子對δg指數，體現兩個片段間的原子對對總相互作用的貢獻
  Atom pair dg index (zero terms are not shown)
    9   20 :    0.215895
    7   20 :    0.079773
    9   14 :    0.078497
...略

上面給出的定量數據靠不靠譜？我們對照前面給出的帶著原子序號的δg_inter圖，可以看到δg指數較大的原子確實都離δg_inter等值面較近。原子對δg指數中9-20是最大的，這倆原子確實緊緊夾著對應氫鍵的δg_inter等值面。7-20、9-14等原子對的距離也都不很遠，而且連線穿過δg_inter等值面，這是它們的原子對δg指數不太小的原因。可見，程序給出的原子和原子對δg指數都是對考察弱相互作用很有用的指標。

值得一提的是，也不要過度解釋上面這兩個指數、夸大這兩個指數的意義，這倆指標和衡量相互作用程度最關鍵的相互作用能并沒有直接的內在聯系。如果你真的需要準確討論原子和原子對對弱相互作用能的貢獻，比較簡單的做法是利用分子力場來計算，見《使用Multiwfn做基于分子力場的能量分解分析》（http://sobereva.com/442）。有的原子電荷計算方法如RESP、ADCH需要波函數信息，如果你只能提供結構信息，則可以用Multiwfn里的電負性均衡方法(EEM)計算原子電荷，此方法計算飛快，幾百個原子體系瞬間就能算完，默認參數下得到的EEM電荷與B3LYP/6-31G*下的CHELPG電荷相仿佛，很適合描述原子間靜電相互作用。

輸出完atmdg.txt后，Multiwfn還問你是否把atmdg.pdb輸出到當前目錄下。這個pdb文件里beta值和占有度（occupancy）字段分別記錄的是原子δg指數的10倍和百分比原子δg指數。如果你將這個文件載入VMD，并且按照beta或occupancy屬性進行著色的話，那么不同原子對弱相互作用的貢獻將可以通過顏色一目了然地考察。這里我們輸入y，讓程序輸出此文件，然后將文件載入VMD，進入Graphics - Representation，Drawing method設CPK，Coloring Method設Beta，此時看到圖形窗口中不同原子已經根據δg指數著色了。默認的色彩刻度是程序根據當前情況自動判斷的，想改的話可以點擊Trajectory標簽頁然后調整Color Scale Data Range。我們把色彩刻度上限從自動判斷的值改為更大的值1，使得色彩變化更柔和一些，然后看到的就是下圖的效果。VMD默認的色彩刻度是紅-白-藍（RWB）方式變化，因此下圖中顏色越紅的原子δg指數越小，越白的原子δg指數越大。如果想調整色彩刻度定義就進入Graphics - Colors - Color Scale來修改。

 

根據原子δg指數著色的分子結構圖和δg_inter等值面圖也可以繪制到一起。即首先按照前述方法把填色的δg_inter等值面圖用VMD腳本顯示出來，再載入atmdg.pdb，然后以上述方式修改顯示設定，得到根據原子δg指數著色的CPK風格的分子結構圖，并且把顯示δg_inter等值面圖的那個ID里面的顯示分子結構的那個表示(Rep)給刪除掉。最后效果將會是下面這樣：

 

之所以上面的圖對原子δg指數的著色和之前不一樣了，是因為VMD里色彩刻度只能用一種。由于IGM_inter.vmd腳本里面有語句把色彩刻度變化設成了藍-綠-紅(BGR)，因此此時的圖像中，原子δg指數越大的原子顯得越綠，越小的顯得越藍。繪制上圖時我還在Graphics-Representation里把顯示等值面的Material設成了AOShiny，這樣填色效果比默認的Opaque材質更不容易受到光線的影響。

Multiwfn在輸出atmdg.txt的同時還輸出我提出的IBSIW（intrinsic bond strength index for weak interaction）指數，這借鑒了《Multiwfn支持的分析化學鍵的方法一覽》（http://sobereva.com/471）中介紹的用于衡量化學鍵強度的IBSI指數的思想。IBSIW相當于將原子對δg指數除以兩個原子間以埃為單位的距離再乘以100。經過筆者粗略測試，IBSIW比起單純的原子對δg指數在衡量原子間弱相互作用強度上往往更合理一些。Multiwfn把IBSIW輸出到當前目錄下的IBSIW.txt中，其中不僅有每對原子的數值，也有把各個原子涉及到的IBSIW都加和的數值。

4 實例：考察主-客體復合物中的弱相互作用

這個例子我們考察一個主-客體復合物之間的弱相互作用。結構文件是本文提供的文件包里的host-guest.xyz。

我們要把主體和客體分子各定義為一個片段，所以我們得知道兩個分子各自的原子序號范圍。如何快速得到一個大體系中某個分子的序號范圍？這可以利用gview6的片段選擇功能，這里示例一下。gview6不支持xyz文件，但支持pdb文件，所以我們先得轉換一下格式，用Multiwfn就可以做這個轉換。啟動Multiwfn，載入host-guest.xyz，然后依次輸入
100  //主功能100
2   //導出文件
1   //導出為pdb格式
C:\host-guest.pdb  //導出的路徑

用gview6載入此文件，然后在一個客體原子上點右鍵，選擇Select Fragment of Atom xxx，此時與這個原子有鍵連關系的所有原子都會被選中成為黃色，然后選主菜單中的Tools-Atom Selection，在新窗口中會看到原子的編號范圍，是127-153，把這行字符復制下來，之后就可以直接粘貼到Multiwfn窗口中了（對有的體系，選中的區域的原子序號是不連續的，gview會給出比如1-5,8,20-30這樣的字符串，這種字符串也同樣可以直接粘貼到Multiwfn里定義片段）。顯然，當前體系里其它原子，即1-126，就是主體分子的序號范圍了。

我們開始做這個體系的IGM分析，在Multiwfn中載入host-guest.xyz，然后依次輸入
20
10
2
1-126
127-153
-10  //修改延展距離
0  //把默認的延展距離從默認的2 Bohr改為0
2  //中等質量格點

此例修改延展距離是因為根據體系結構可知，當前體系中弱相互作用區域肯定都在體系內部區域，因此計算IGM格點數據的時候用的盒子沒必要在體系四周往外擴一定距離。在計算的點數相同的情況下，盒子越小，則均勻分布在盒子里的格點之間的間隔就越小，圖像效果就會越好，等值面會越平滑。

本例我們就看一下填色的δg_inter等值面圖。按照上一節例子用IGM_inter.vmd腳本繪圖，我們發現圖上只出現了很小很零散的等值面，并沒有很好體現主-客體之間的相互作用。這是因為前面提到過，相互作用越弱的區域δg越小，由于當前體系主-客體之間是范德華作用，強度很弱，因此如果不把δg_inter等值面數值設小的話是基本看不到什么的。我們遂在Representation界面里把Isovalue改為較小的值，比如0.003，此時看到的圖像如下。作為對比，下圖也把同樣格點設定下的基于準分子密度的RDG填色等值面圖附上了（是用Multiwfn主功能20的子功能2來產生格點數據然后由VMD繪制的，具體步驟看手冊4.20.2節）。
 

圖像效果好壞高下立見！同樣格點間距情況下，IGM的圖圓潤豐滿富態得多，而RDG的圖則頗為毛糙，鋸齒明顯，顯得十分貧窮。要改善RDG圖的圖像效果，需要用明顯更小的格點間距，但代價就是要算的格點數將大為增加。雖然把RDG等值面數值稍微改大點可以讓等值面更平滑、連續一些，但代價就是周圍會出現很多多余的等值面，把圖像弄臟。

下面我們把δg_inter等值面圖和根據原子δg指數著色的結構圖繪制在一起。由于當前體系弱相互作用幾乎都是色散作用，根據sign(lambda2)rho填色后的顏色基本都是綠色，所以要不要填色效果無所謂，所以我們繪制時候只把δg_inter顯示成等值面而不填色，這樣在根據原子δg指數對結構著色的時候就可以隨意設定色彩刻度了。讓Multiwfn產生atmdg.pdb文件，由于體系比苯酚二聚體略大了一些，所以計算δg指數過程的耗時也高了不少。把dg_inter.cub和atmdg.pdb都拖入VMD，在VMD里經過一番調節，最后得到了下圖，效果相當不賴！圖中主體分子中越紅的地方對主-客體相互作用貢獻越大，看著非常直觀。

 

為了獲得上圖，在VMD中具體做的設定和調節如下
(a)在Graphics-Representation中，切換到對應atmdg.pdb的頁面，設立兩個Rep，Selected Atoms分別設為fragment 0和fragment 1，分別對應主體和客體分子（在VMD中，相互鍵連的每批原子都有一個獨立的fragment號）。選中主體分子對應的Rep，Drawing Method設Licorice，Material設EdgyShiny，Coloring Method設Beta，然后在Trajectory標簽頁里把色彩刻度下限和上限分別設為-0.2和0.2。之后選中客體分子對應的Rep，Drawing Method設CPK，Coloring Method用Name，Material設Glossy，然后把Sphere Scale設為1.2，Bond Radius設0.9。
(b)在Graphics-Representation中，切換到對應dg_inter.cub的頁面，把Drawing Method設為Isosurface，Isovalue設為0.003，Drawing設為Solid Surface，Show設為Isosurface。然后把Coloring Method設ColorID并且在旁邊的下拉框中選7 Green。最后把Material設為Edgy Glass。
(c)在Graphics-Color里選Color Scale標簽頁，把Method從默認的RWB改為BWR，這樣原子δg指數從0.0到0.2范圍的變化就通過由白變紅方式展現了。

5 實例：考察Actos四聚體中的相互作用

Actos是一種藥物分子，筆者過做分子動力學模擬了此體系在真空下的四聚體（共180個原子）。這里筆者隨便從中抽出一幀，用Grimme的xtb程序作了優化，得到的最終結構是本文的文件包里的4Actos.xyz。本節我們考察一下其中單體間的弱相互作用，由此進一步展示IGM分析的片段設定。

首先我們考察體系中所有單體間的弱相互作用對應的δg_inter填色等值面圖，因此我們要把每個單體各定義為一個片段。啟動Multiwfn，載入4Actos.xyz，然后依次輸入
20
10
4   //四個片段
1-45  //單體1的序號
46-90  //單體2的序號
91-135  //單體3的序號
136-180  //單體4的序號
3   //鑒于體系不小，為了保證等值面質量，這次我們用高質量格點
3   //導出cube文件
此時產生的dg_inter.cub和dg_intra.cub分別描述的是這個四聚體當中所有單體間的弱相互作用和所有單體內的弱相互作用。然后和之前的例子一樣用IGM_inter.vmd和IGM_intra.vmd繪制等值面圖，如下所示。左圖是腳本直接繪制出來的，右圖是為了便于辨別分子，把不同單體用不同顏色顯示了
 

圖像把各處弱相互作用出現區域和特征都很好地展現了，但是看著終究有點亂。我們可以把弱相互作用一對一對地考察。比如我們只想把上圖中紅色和肉色分子之間的δg_inter等值面展現出來，那么我們可以選0退出IGM分析界面，然后依次輸入
10  //再次進入IGM分析模塊
2   //定義兩個片段
1-45  //紅色的分子對應的原子序號
135-180   //粉色的分子對應的原子序號
2   //中等質量格點（用高質量格點當然也可以。通過對比會看到，中等質量格點的等值面也不比之前用高質量格點的差多少，體現出IGM對格點質量依賴性較低的優點）

然后和之前介紹的一樣利用IGM_inter.vmd繪制填色等值面圖。在VMD里顯示出來后，我們把分子結構用兩個Rep來顯示，都用Licorice風格。第一個Rep的選擇范圍寫的是fragment 0 3，用來顯示我們感興趣的那兩個分子，Bond Radius設0.2；另一個Rep的選擇范圍寫的是not fragment 0 3，用來顯示另外兩個分子，并且為了讓它們存在感低一些，把Material設成Transparent使它們透明，且Bond Radius設為更小的0.1。此時看到的圖像如下：

 

圖像效果非常理想，僅有感興趣的兩個分子之間的弱相互作用被展現出來，其它區域的弱相互作用完全沒有出現在圖中擾亂視覺。IGM比RDG方法關鍵性的好處就是指哪打哪，而且打的地方還可以進一步分解成原子和原子對的貢獻。我們讓Multiwfn把當前情況對應的atmdg.pdb產生出來，拖到VMD里以Licorice風格顯示并用Beta字段著色，并把原先顯示fragment 0 3的那個Rep取消顯示，看到如下效果

 

當前色彩刻度變化是藍-綠-紅，因此對于圖中以不透明方式顯示的兩個分子，藍色代表原子的δg指數接近于0，故對兩個分子間的弱相互作用幾乎沒有貢獻；綠色代表原子的δg指數不小，明顯參與了弱相互作用；桔色甚至紅色代表相應原子的δg指數很大，仔細圖像的話，會發現原因是這些原子都參與了不止一對分子間弱相互作用，因此值得特別關注。

當前這個四聚體還可以以其它方式劃分片段。比如我們可以把分子1作為一個片段，而分子2、3、4作為第二個片段，那么IGM分析就能專門考察1與2-3-4三聚體之間的弱相互作用，而三聚體內部的弱相互作用則不會搗亂。我們也可以把其中兩個分子作為一個片段，另兩個分子作為另一個片段，這樣兩個二聚體之間的弱相互作用就會被展現，而二聚體內部的弱相互作用則不會呈現。IGM方法可謂極盡靈活。

6 實例：考察DNA中堿基之間的弱相互作用

前面的例子都是考察分子間弱相互作用，利用IGM也可以考察分子內弱相互作用，只要恰當定義片段即可。這里我們以考察DNA中堿基之間的弱相互作用作為例子。通過分子力場優化后的DNA片段是文件包里的DNA_single.pdb。為了看起來更清楚一點，這個DNA只保留了一條鏈。這個文件里有10個核苷酸，我們這里只考察靠中部位置的三個堿基之間的相互作用，即下圖中用球棍方式繪制的三個堿基

 

我們還是可以利用gview選中堿基部分，然后在Tools-Atom Selection里方便地看到編號。這三個堿基的原子編號范圍分別為107-120、139-152、171-184。啟動Multiwfn，載入DNA_single.pdb，然后依次輸入
20
10
3  //定義三個片段
107-120
139-152
171-184
10  //在圖形窗口中設定盒子
為了節約計算量而且省內存，我們在圖形窗口中通過調節盒子尺寸和位置，使得盒子差不多剛好套住可能出現弱相互作用的區域，如下圖這樣
 

關閉圖形窗口，程序就開始計算格點數據了。之后把格點數據導出，用IGM_inter.vmd腳本在VMD里作圖，然后在Graphics-Representation里把等值面數值改為0.005，看到的圖像效果如下

 

可見圖像很干凈，等值面很平滑，很好地展現出了堿基之間的pi-pi堆積作用。

7 考察BCP位置的δg衡量相互作用強度

AIM理論中的鍵臨界點(BCP)一般被認為是原子間相互作用路徑上最有代表性的點，通過考察這個點的實空間函數可以用來表征相互作用強度、特征。IGM原文表1中指出BCP位置的δg和相互作用強度有很好正相關性。我們這里對HF線型四聚體當中三個氫鍵對應的BCP的δg進行一下考察。如果對Multiwfn做拓撲分析搜索臨界點不熟悉，可以參看《使用Multiwfn做拓撲分析以及計算孤對電子角度》（http://sobereva.com/108）。

本文的文件包里4HF.wfn是基于M062X/TZVP下進行結構優化，然后在M062X/def2-TZVP下產生的波函數文件。我們先啟動Multiwfn，載入4HF.wfn，然后依次輸入
2  //拓撲分析
3  //以每兩個原子連線中點為初猜搜索臨界點，一般能夠找到所有的BCP
0  //顯示分子結構和臨界點

 

由圖可見對應三個氫鍵的三個BCP的編號分別是2、4、6。關閉圖形窗口，進入選項7，然后輸入要查看屬性的BCP序號。我們先輸入2，從屏幕上的數據我們可以直接讀到包括δg在內的各種實空間函數值，部分貼在這里：
Density of all electrons:  0.3332066650E-01
Potential energy density V(r): -0.3698198590E-01
Delta_g:  0.7562586676E-01

然后再次進入選項7，輸入4，得到對應中間氫鍵的BCP的性質
Density of all electrons:  0.3714890342E-01
Potential energy density V(r): -0.4270021580E-01
Delta_g:  0.8258706914E-01

之后再查看6號BCP的性質
Density of all electrons:  0.3110958174E-01
Potential energy density V(r): -0.3333690144E-01
Delta_g:  0.6951460311E-01

由數據可見，δg的大小順序是BCP4 > BCP2 > BCP6，因此體現出氫鍵強度應該是F7-H8...F1 > F1-H2...F3 > F5-H6...F7。對于同類鍵來說，BCP電子密度越大、勢能密度越負則相互作用強度越大。由上面給出的數據可見，通過δg體現的強度順序和電子密度、勢能密度體現的強度順序是完全一致的。而且如果考察氫鍵鍵長的話，根據鍵長越短鍵強越大的這個多數情況奏效的規律，結論也是和δg判斷的是一致的。

8 總結&其它

本文介紹的IGM方法是對流行的RDG（NCI）方法的一個重要的補充。IGM的優點主要有以下幾條
(1)可以很自由地定義片段，沒納入片段的部分不會被考慮，而且片段內和片段間的相互作用可以明確地分離展現，這使得圖像可控性很強，想展現什么就展現什么，不像RDG方法那樣各種相互作用同時出現在圖上，還得想辦法對格點數據做屏蔽或者ps來消除不感興趣又擾亂視覺的部分。
(2)可以通過計算原子和原子對δg指數定量反映原子和原子對對片段間相互作用起到的影響程度，而且原子的貢獻還可以通過著色一目了然地在VMD中展現出來，便于識別關鍵原子。
(3)IGM等值面比RDG等值面平滑豐滿得多，因此用較少數目的格點就可以達到滿意的圖像效果，計算耗時和內存消耗量因此較低。
(4)不同區域的δg的大小可以體現相互作用強度，因此可以通過調節δg等值面數值來決定是只顯示較強的相互作用，還是所有相互作用都顯示出來。另外，還可以通過考察BCP處的δg值來定量衡量特定相互作用的強度。

但是，IGM也并非能完全代替RDG，有以下兩點原因
(1)IGM雖然有基于波函數的版本，但是定義得比較抽象，效果也沒比基于promolecular密度的IGM方法更好，因此Multiwfn也沒打算支持之。對于對精度要求較高，特別是相互作用區域的實際密度可能與promolecular密度相差較大時，RDG還是應當優先考慮。
(2)IGM不適合用來考察體系中所有弱相互作用。因為考察這個需要把整個體系（或某個局部區域的原子）定義為一個片段，然后考察這個片段的δg_intra。但是，這樣呈現出的圖像要么只顯示化學鍵區域，要么顯示的是化學鍵+弱相互作用區域。如果想只把弱相互作用區域顯示出來，就需要根據sign(lambda2)rho的范圍將sign(lambda2)rho很負的區域進行屏蔽（用IGM后處理菜單中的5 Screen delta_g_intra at high density region實現），但實際發現這么的做效果并不如RDG好，往往是isovalue調到能令弱相互作用等值面能較好顯示時，其它區域的一些雜碎等值面又出來了。另一種做法是直接計算整體的δg函數格點數據，再用Multiwfn的主功能13對距離原子較近區域的δg函數的格點數據進行屏蔽，從而只保留弱相互作用區域，但這樣終究稍微麻煩點，不如直接用RDG。

在Multiwfn手冊的4.20.10節里還有關于IGM方法的更多例子，建議大家閱讀本文后閱讀。另外，在《18碳環（cyclo[18]carbon）與石墨烯的相互作用：基于簇模型的研究一例》（http://sobereva.com/615）介紹的筆者的一篇文章中還用IGM考察了18碳環與石墨烯之間的相互作用，效果很好，建議一看。

IGM的作者在2020年基于δg還定義了一個叫IBSI（內稟鍵強指數）的量，用于衡量化學鍵強度，Multiwfn也可以計算。關于這個方法請在《Multiwfn支持的分析化學鍵的方法一覽》（http://sobereva.com/471）搜索IBSI查看介紹。

還值得一說的是，如果你用的輸入文件包含波函數信息（如.wfn、.fch、.molden文件），則進行IGM分析的時候程序會讓你選擇sign(lambda2)rho的類型，第一種是基于實際電子密度計算的，第二種是基于promolecular密度近似計算的。原理上來說，用前者時散點圖、等值面圖的填色效果結果比后者會更有意義（但文中涉及的腳本里的色彩刻度在某些時候可能需要進行一些修改），而耗時也會高得多。當輸入文件只含結構信息時，sign(lambda2)rho總是基于promolecular密度近似計算的。

之前在http://sobereva.com/399里介紹過怎么繪制有填色效果的RDG散點圖，這很便于將散點圖上的位置與填色等值面圖相互對照，對于IGM散點圖也完全可以繪制類似的效果，只不過需要用的gnuplot程序的作圖腳本有所不同。具體做法見手冊4.20.10.1節末尾。效果如下所示

在http://bbs.keinsci.com/forum.php?mod=redirect&goto=findpost&ptid=9472&pid=149548還有人分享了基于Python的繪制這種圖的腳本，這樣就不需要裝gnuplot了。