使用Multiwfn+VMD以原子著色方式表現原子電荷、自旋布居、電荷轉移、簡縮福井函數
使用Multiwfn+VMD以原子著色方式表現原子電荷、自旋布居、電荷轉移、簡縮福井函數
文/Sobereva @北京科音 2018-Jun-3
0 前言
在圖像中表現原子電荷分布、原子自旋布居的時候通常都是將數值標在原子上面，但當原子很多的時候，就會導致圖像上的數字密密麻麻，難以辨認，而且也不好直觀去考察。有一種做法可以解決這個問題，就是通過原子的顏色來表現其原子電荷、自旋布居等屬性的正負和大小，這樣非常直觀。在VMD程序中支持這種著色方式，在筆者之前一個文章中也利用過這種做法，見《將GROMACS的原子電荷信息讀入VMD的方法》（http://www.shanxitv.org/365）。在筆者的《通過獨立梯度模型(IGM)考察分子間弱相互作用》（http://www.shanxitv.org/407）一文中也利用了原子著色方法，表現了各個原子對指定的弱相互作用的貢獻。
Multiwfn程序主功能7在計算完原子電荷后可以導出記錄原子坐標和原子電荷的.chg格式的文件，它和常見的.xyz格式的唯一區別是每個原子后面增加了一列，用于記錄原子電荷。對于2018年6月3日及以后更新的Multiwfn程序，如果以這種.chg文件作為輸入文件，那么在主功能100的子功能2當中可以導出.pqr格式。pqr格式和pdb格式關鍵差異在于多增加了兩列，一列記錄原子電荷，另一列記錄原子半徑。.pqr是可以直接載入到VMD程序中的，而且載入后可以直接根據原子的電荷值來著色。本文就舉幾個例子完整說一下操作，并且把這種做法用于展示原子自旋布居、原子電荷的變化、簡縮福井函數上。
Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn上免費下載，本文用的是2018年6月3日更新的3.6(dev)版。如果對Multiwfn不了解，看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn波函數分析程序的意義、功能與用途》（http://www.shanxitv.org/184）。如果不了解原子電荷，看《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）。本文用的VMD是1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免費下載。
本文涉及到的各種文件可以在這里下載：http://www.shanxitv.org/attach/425/file.rar
1 用著色方式展現Au20團簇的原子電荷
此例我們用原子著色方式展示一下正四面體的20個Au組成的Au20團簇的原子電荷分布，這里用筆者提出的普適性和合理性良好的ADCH原子電荷，原文見http://dx.doi.org/10.1142/S0219633612500113。啟動Multiwfn，依次輸入
Au20.fch
7 //布居分析
11 //ADCH電荷
1 //使用內置的球對稱化密度計算ADCH電荷
算完之后，程序問你是否在當前目錄下導出Au20.chg文件，選擇y。之后，重啟Multiwfn（也可以退回到主菜單然后選隱藏選項-11來重新進入選擇輸入文件的界面），載入Au20.chg。之后進入主功能100，選擇子功能2，此時的選項1對應于輸出.pqr文件，選擇這個選項，然后輸入Au20.pqr。此時Au20.pqr就在當前目錄下產生了。此文件可以用文本編輯器打開，其倒數第一列是原子的Bondi范德華半徑，倒數第二列就是.chg文件里記錄的ADCH原子電荷了。
啟動VMD，把Au20.pqr拖入VMD Main窗口載入，進入Graphics - Representation，把Coloring Method設為Charge（按照原子電荷那一列著色），把Drawing Method設為VDW，此時看到各個原子以不同顏色展示了，如下所示。
接下來我們調節選項來改進圖像效果：
(a)關閉坐標軸：選擇Display-Axes-Off
(b)把背景設為白的：在VMD的命令行窗口輸入color Display Background white
(c)修改色彩刻度變化方式：默認是RWB方式色彩變化，即紅-白-藍，我們選擇Graphics-Colors-Color scale，Method設為BWR，這樣原子電荷從最小變化到最大就對應于藍-白-紅方式的色彩變化了
(d)修改圓球尺寸和材質：在Representation界面里把Sphere Scale設為0.5，然后把Material設為Glossy，這比默認的材質更有光澤感
(e)修改色彩刻度：在Representation界面的Trajectory標簽頁里把Color Scale Data Range下限和上限分別設為-0.15和0.15然后按回車（默認的色彩刻度是卡著最小值和最大值設定的，即-0.06~0.11，此時顏色變化過大，不那么柔和。我們把數值范圍設大一點，而且讓上限和下限絕對值相同，這樣色彩更舒服，而且還便于直接通過偏藍還是偏紅來判斷電荷是正值還是負值）
(f)增加鍵連顯示：在Representation界面里點擊Create Rep，Coloring Method設Color ID，旁邊的框選32 orange，Material設AOEdgy，Drawing Method設DynamicBonds，把Distance Cutoff設3.6，Bond Radius設0.2
此時效果應該已經很好了，應確保已經開了GLSL（Display-Rendermode-GLSL），這樣效果才能充分展現。如果要效果更好，可以用內置的Tachyon渲染器渲染，見此文的說明《用Multiwfn+VMD做RDG分析時的一些要點和常見問題》（http://www.shanxitv.org/291）。Tachyon渲染后的圖像如下
從圖中可見處于每個面中間和頂角的金帶微量正電荷，其它原子帶微量負電荷。由于這個體系中每個原子帶電量很小，因此不同原子電荷計算方法算出的結果正負號可能不同，沒必要太糾結于此。
值得一提的是，平時常用的gview也能根據原子電荷進行著色。載入Gaussian輸出文件或者fch文件，選Results - Charge distribution，按照如下設置即可顯示（當前顯示的是NPA電荷）
然而，跟VMD的效果一比，gview的效果簡直丑爆了。不僅色彩變化方式沒法改，而且模型很難看、缺乏靈活的可調節設定。
2 用著色方式展現丁烷雙自由基的自旋布居
如果不了解自旋布居，看《談談自旋密度、自旋布居以及在Multiwfn中的繪制和計算》（http://www.shanxitv.org/353）。這一節我們通過原子著色方式直觀展現丁烷雙自由基的自旋密度在各個原子上的凈分布量。這個體系在以前的博文中出現過多次，見《CASSCF計算雙自由基以及雙自由基特征的計算》（http://www.shanxitv.org/264）、《談談片段組合波函數與自旋極化單重態》（http://www.shanxitv.org/82）、《在Multiwfn中基于fch產生自然軌道的方法與激發態波函數、自旋自然軌道分析實例》（http://www.shanxitv.org/403）。本文用的丁烷雙自由基的fch文件是用對稱破缺DFT計算得到的。我們首先計算這個體系的原子自旋布居，這里就用常用的Becke方法計算。
啟動Multiwfn，依次輸入
C4H8.fch
15 //模糊空間分析
1 //對原子空間積分。默認是Becke方式劃分原子空間
5 //積分自旋密度
然后按照Multiwfn手冊5.4節的方法，把屏幕上輸出的自旋布居拷出來：
0.81112088
0.00391514
0.00391514
-0.04754640
0.01069745
0.01069745
0.04754640
-0.01069745
-0.01069745
-0.81112088
-0.00391514
-0.00391514
之后，我們回到主菜單，進入主功能100的子功能2，選擇導出.xyz文件。手動編輯此文件，把開頭兩行去掉（原子坐標和注釋行），然后利用Ultraedit等專業文本編輯器的列模式，把上面的自旋布居粘到.xyz文件最后一列，之后把文件后綴改為.chg。處理好的文件就是本文的文件包里的C4H8.chg。
把C4H8.chg載入Multiwfn，按照上一節的做法轉換成C4H8.pqr。此時此文件里原子電荷那里一列就對應于原子自旋布居了。（值得一說的是，xyz和chg都是自由格式，手動編輯方便，而pqr文件是固定格式，手動編輯不便，所以此例才先手動編輯得到chg文件再轉換成pqr文件）
把C4H8.pqr拖到VMD里，Drawing Method用CPK，Material用Diffuse，還是根據Charge進行著色，把色彩刻度范圍從默認的-0.81~0.81改小到-0.2~0.2，這樣可以讓自旋布居絕對值比較小的原子也能稍微有一點顏色（否則中間兩個原子由于自旋布居太小，完全是白的）。此時看到的圖像如下
此圖展現出未成對電子主要分布在兩端的碳原子上，而且自旋方向相反。而中間的碳原子也帶微量單電子，所以不是純白色。
從上圖中看不出各個原子是什么元素，如果通過顏色又想體現出元素，又想展現出原子屬性值，那么可以用兩個Rep。我們把目前已經有個一個Rep刪掉，即點擊Delete Rep按鈕，然后點Create Rep建立一個新Rep，用Charge著色，用VDW風格顯示，材質設Transparent，Sphere Scale設0.3，色彩范圍還用-0.2~0.2。然后再點擊Create Rep，用Name著色，用Licorice風格顯示，材質用Opaque。然后進入Graphics-Material，把Transparent的Opacity設0.7降低其透明度。此時用自帶的Tachyon渲染器渲染出的圖如下
可見元素和原子的自旋布居大小同時展現了出來，互不干擾。
3 根據原子電荷體現示意鳥嘌呤-胞嘧啶堿基對中的靜電相互作用以及電荷轉移
本文文件包中GC目錄中的GC.xyz是鳥嘌呤(G)-胞嘧啶(C)堿基對的結構，兩個堿基之間通過氫鍵維持二聚體的穩定結合。在《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）中提到，氫鍵的主要本質就是靜電相互作用，多數情況靜電相互作用又可以通過原子電荷之間的庫侖作用近似描述，本例的體系就是這種情況。本例將對GC中的兩個堿基都用原子著色方式展現其原子電荷，使得靜電相互作用可以較直觀地考察。
本例繪制過程不能像本文第1節那樣直接對二聚體計算原子電荷然后繪圖。因為兩個單體結合在一起的過程中，原子電荷也會發生變化。對二聚體計算的原子電荷，相當于是已經發生了電荷轉移、極化之后的情況了，而不是體現兩個單體原本的電荷分布。而要想描述單體間的靜電相互作用，我們應當用單體原本的電荷分布來表現。這類似于筆者之前在《使用Multiwfn結合VMD分析和繪制分子表面靜電勢分布》（http://www.shanxitv.org/196）中通過靜電勢填色的范德華表面穿透圖展現水二聚體中的靜電相互作用，當時就是對兩個單體分別計算產生的范德華表面，而非是對二聚體進行的計算。
我們首先把G、C的坐標從xyz文件中提出來，改寫成.gjf文件，關鍵詞用# M062X/6-311G** nosymm。注意這里nosymm很重要，不寫這個的話，單體計算時Gaussian會自動平移和旋轉體系坐標到標準朝向，這樣得到的單體的fch文件里的坐標就和復合物不同了，之后繪圖會比較麻煩。關于nosymm此文有更多討論：《談談Gaussian中的對稱性與nosymm關鍵詞的使用》（http://www.shanxitv.org/297）。
對G.gjf和C.gjf用Gaussian計算，得到G.fchk和C.fchk，然后再用Multiwfn計算，產生記錄了ADCH電荷的G.chg和C.chg，之后再轉換為G.pqr和C.pqr。之后，把這倆文件拖進VMD，在Representation界面里把顯示設定都設為相同，都用CPK風格，用charge著色，色彩刻度都設成-0.5~0.5，用BWR色彩過渡方式，材質用EdgyShiny，看到的圖像如下
從圖中的原子顏色可見，兩個單體是以原子電荷正負互補的方式結合的，即帶正電的氫（粉色）與帶負電的氮或氧原子（藍色）緊密接觸，這種正負電荷間的靜電吸引正是氫鍵作用的主要本質。
接下來，我們通過原子著色直觀考察一下形成二聚體時各個原子電荷的變化情況，這要計算二聚體狀態的原子電荷與單體的原子電荷的差值。筆者發現用ADCH電荷（原子偶極矩校正的Hirshfeld電荷）計算這個體系的二聚過程的原子電荷變化不太理想，和密度差圖差異較大，原因是因為Hirshfeld權重函數本身不太理想（過于平滑、隨徑向距離收斂太慢），如果改用權重函數收斂更快的Becke劃分就沒這個問題了。因此下面使用原子偶極矩校正的Becke電荷考察當前問題，計算這種電荷就是在Multiwfn主功能7里面選子功能10。
對GC.fch進行計算，得到GC.chg，然后把GC的原子電荷，以及對G和C分別計算的原子電荷都拷到Excel里（文件包里的charges.xls），令GC的電荷減去G和C的原子電荷，從而得到形成二聚體時原子電荷的變化量。然后把GC.chg復制成GCdiff.chg，把其中的原子電荷替換為剛才算出來的原子電荷的變化量。之后把GCdiff.chg轉化為GCdiff.pqr，用VMD作圖，色彩刻度用-0.1~0.1，色彩刻度用BWR，得到下圖
紅色和藍色的原子分別對應于形成二聚體后原子電荷增加和減小的原子。可見形成二聚體后氫鍵中的氫原子失去了電子，原子電荷變得更正，因此呈紅色；氫鍵受體重原子呈藍色，說明原子電荷有所減小，明顯得了電子。氫鍵給體重原子也同樣得了一定電子，但得的量很少，因此呈淡藍色。
可以用Multiwfn繪制此二聚體的片段密度差圖，方法見《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113），所得圖像如下，紅色實線和藍色虛線分別對應形成二聚體時電子密度增加和減小的部分。
可見密度差圖和上面根據原子偶極矩校正的Becke電荷的結論一致，在氫鍵的氫附近確實電子密度減小了。
4 對原子根據簡縮福井函數著色
之前筆者在《基于Multiwfn產生的cube文件在VMD和GaussView中繪制填色等值面圖的方法》（http://www.shanxitv.org/402）中舉了一個繪制暈苯的福井函數f-的例子，沒看過的一定要看一下。這回我們計算簡縮福井函數f-，并用原子著色方式展現。不了解福井函數者參看“概念密度泛函綜述和重要文獻合集”（http://bbs.keinsci.com/thread-384-1-1.html）。
暈苯的N電子和N-1電子狀態的fch文件都在本文的文件包里的coronene目錄下，對這兩個fch文件都計算Hirshfeld原子電荷。用這個電荷是因為在筆者寫的《親電取代反應中活性位點預測方法的比較》（物理化學學報, 30, 628 (2014) http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）等文章中都體現Hirshfeld電荷在各種原子電荷中相對來說較適合討論反應位點問題。然后把N和N-1態的原子電荷放到Excel里，令N-1態的原子電荷減去N態的原子電荷，即得到簡縮福井函數f-，將其數值寫入到對應于此體系的chg文件的最后一列（即coronene.chg），再轉換為coronene.pqr。用VMD作圖，以CPK顯示，色彩刻度用-0.07~0.07，色彩變化用BWR，得到下圖
圖中越紅的原子，擁有越大的簡縮福井函數f-，因此越容易發生親電反應。可以將此圖和《基于Multiwfn產生的cube文件在VMD和GaussView中繪制填色等值面圖的方法》里的f-等值面圖以及把f-投影到分子表面上的圖對比，會看到傳達的信息是相同的，只不過表現形式不同而已。
注：明顯更簡單、省事的計算簡縮福井函數的做法是用Multiwfn的專用功能，敲幾下鍵盤就能直接算出來所有類型的簡縮福井函數和簡縮雙描述符，見《使用Multiwfn超級方便地計算出概念密度泛函理論中定義的各種量》（http://www.shanxitv.org/484）。
最后值得一提的是，本文這種對原子著色表現原子屬性的做法是極度普適的，絕不僅僅限于文中討論的情況。比如還可以用這種方法表現Multiwfn計算出來的原子躍遷電荷、原子對分子軌道（或NTO等其它軌道）的貢獻、原子對空穴/電子分布的貢獻、原子偶極矩的大小、原子的電子動能（原子空間內對動能密度的積分）、原子的源函數值、原子的簡縮雙描述符、電子激發前后電荷的變化、外加電場前后電荷的變化等等。而其它程序計算的原子極化率、原子受力大小等也都可以用這種方式展現。