使用Multiwfn可視化分子孔洞并計算孔洞體積
注1：本文的方法適合展現分子內部的孔洞。如果你是希望展現周期性動力學模擬或者晶胞當中的自由區域（未被原子占據的區域）并計算其體積，應當使用此文的方法：《使用Multiwfn圖形化展示分子動力學模擬體系中的孔洞、自由區域》（http://www.shanxitv.org/539）和《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》（http://www.shanxitv.org/617）。
注2：Multiwfn還有計算分子孔洞直徑的功能，見《使用Multiwfn計算分子和晶體中孔洞的直徑》（http://www.shanxitv.org/643）。
使用Multiwfn可視化分子孔洞并計算孔洞體積
文/Sobereva @北京科音
First release: 2018-Feb-28 Last update: 2021-Jul-26
1 前言
前些天有個人在思想家公社QQ群里問這種圖怎么繪制，怎么計算孔洞體積
其實有很多程序和在線服務器都是專門用來顯示孔洞、口袋并計算其體積的，但是這些程序絕大部分都是給生物分子用的，而對于小分子體系，要么完全不適用要么不好用。而利用十分靈活的Multiwfn，可以很容易地顯示出分子孔洞并計算其體積。為了實現這個目的時更方便，筆者在原先Multiwfn基礎上又稍微做了擴展。本文使用的是2021-Jul-26更新的Multiwfn。本文會結合幾個例子講解原理以及說明具體怎么操作。
Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，如果對Multiwfn不了解，看《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。使用本文功能計算孔洞體積發表文章時務必按照程序啟動時顯示的要求恰當引用Multiwfn。
2 原理
何以定義分子孔洞？一個簡單的做法是，把每個原子在自由狀態的電子密度進行疊加成為promolecular density（準分子密度），相當于原子出現在分子中對應的位置，但電子密度還沒有因為成鍵而分布發生弛豫的狀態。準分子密度的某個等值面內的區域我們認為是分子內的區域，而這個等值面外部的區域就可以作為分子孔洞了。有些孔洞是封閉的，比如富勒烯內部區域；但有些孔洞是與外界連通的，比如納米管，對于這種情況，計算準分子密度的時候必須恰當設定盒子（即格點數據計算范圍），讓這個盒子恰好括住孔洞出現的區域，此時既在盒子內部，又在準分子密度等值面外的區域就是孔洞區域了。
周期表里幾乎所有元素在自由狀態的球對稱化的密度都是內置于Multiwfn的，因此Multiwfn可以很迅速、方便地構建準分子密度。計算準分子密度不需要波函數信息，所以輸入文件就用含有結構信息的文件即可，比如.pdb、.xyz、.mol這些常用的記錄分子結構的格式都是Multiwfn支持的。
上面是基本原理，更具體到程序實現時是這樣的：進入Multiwfn的域分析功能（主功能200里的子功能14），將被計算的實空間函數設為準分子密度。然后恰當設定域的定義，令域對應于準分子密度小于某個值（如<0.001或0.0001 a.u.）的空間范圍。然后選擇計算格點數據并做域分析，期間要恰當設定盒子范圍，選擇恰當的格點間距（盒子尺寸相同時，格點間距越小則點數越多，耗時越高、耗內存越多，而體積計算結果也越準確）。程序會計算均勻分布在盒子內每個格點的準分子密度，然后通過特殊的方法將每一批空間上相互連通的格點都賦予唯一的編號，也因此一個編號對應于一個域。之后，在后處理菜單中，就可以選擇查看各個域的分布，找到對應于孔洞的那個域的編號，之后就可以計算相應的域的體積。在Multiwfn里可以通過小圓點疊加的方式顯示指定的域，如果想要更好的繪制效果，還可以把域導出成為.cub文件，其中處于域內部的格點數值為1，外部為0（盒子邊界的點也被設為0以保證等值面看起來總是封閉的），然后可以載入VMD、gview等程序來觀看等值面，把等值面數值設為0~1之間的值比如0.5即可令等值面勾勒出孔洞范圍。
注意孔洞的圖形和體積是沒法唯一定義的，對于Multiwfn這種考察孔洞的方式，結果會受到定義域所用的準分子密度數值的影響。下面是alpha環糊精體系準分子密度=0.001和0.0001 a.u.時候的等值面圖，盒子被設為恰好緊貼著分子邊緣
一般情況為了方便、快速，基于上述準分子密度計算孔洞體積。但如果你能得到體系的波函數文件，想計算孔洞更為準確，不怕花費更多耗時，那么可以把含有波函數信息的文件（如.wfn、.fch等）載入Multiwfn，按照與本文一致的操作即可得到實際密度下的孔洞體積。
3 例：alpha環糊精
本例將計算alpha環糊精中間的孔洞區域體積并繪制孔洞圖形，結構文件在此下載alpha-cyclodextrin.pdb，這是從CCDC結構數據庫中提取的環糊精結構。啟動Multiwfn，載入pdb文件，然后輸入
200 //主功能200
14 //域分析
2 //設定用于定義域的實空間函數
1 //對于只含有分子結構信息而沒有波函數信息的輸入文件，函數1對應于準分子密度
3 //設定域的定義規則
<0.0001 //將電子密度小于0.0001 a.u.的區域作為域
1 //計算格點數據并產生域
-10 //修改延展距離
0 //延展距離被設為0，此時盒子會緊貼著分子邊界設定
1 //低質量格點。如屏幕提示所示，這個模式使用間距為0.2 Bohr的格點，格點數共978588。這樣的精度對當前體系已經足夠用了
接下來我們可以選3來在圖形界面里觀看各個域。在右下角列表里選擇哪個域，哪個域對應的區域就會以小圓球疊加方式顯示出來。我們挨個看后發現對應孔洞的那個域的編號是4號，如下所示，孔洞圖形可以很好地反映當前分子實際的孔洞特征。
然后介紹下怎么把域導出來放到VMD程序里觀看。VMD是最強大、最靈活的化學體系可視化程序，可在http://www.ks.uiuc.edu/Research/vmd/免費下載。在剛才的域分析后處理界面里選擇10，然后輸入4，則4號域就被導出為了當前目錄下的domain.cub。將分子.pdb文件和domain.cub都拖到VMD里，進入"Graphics"-"Representation"，Selected Molecule那一欄切換到對應分子結構的項，把Drawing method設為CPK。然后把Selected Molecule那一欄切換到對應.cub文件的項，把Drawing method設為Isosurface，Draw設為Solid Surface，Show設為Isosurface，Isovalue設為0.5。Coloring Method設為Color ID，假設我們要用黃色就在旁邊選擇4 Yellow。然后在文本窗口輸入color Display Background white把背景變成白的，此時圖像效果如下：
稍微修改作圖設定，弄成本文開頭的圖那種風格也很容易。
另外，我們還可以用后處理界面的選項11把指定的域的邊緣格點導出為當前目錄下的domain.pdb文件，其中每個粒子對應一個邊緣格點。把這個pdb拖進VMD程序，Drawing Method設為sphere，即可在圖形窗口看到各個邊緣格點。之后還可以按鍵盤上的2鍵，然后點擊兩個粒子，這兩個粒子的連線，以及它們的間距，就會直接顯示出來，因此可以十分方便地測量孔洞的各個地方的尺寸，例如下圖。為了看著清楚，這里選擇了Display - Orthographic以使用正交視角。另外，注意在Graphics - Labels里可以設定原子標簽、連線是否顯示。
當前的pdb文件里分子朝向剛剛好，正好是孔洞方向沖著分子Z軸，因此按照上文計算方式可以直接得到對應于孔洞的域。但是，如果pdb文件里這個分子是斜著的，本例的做法就不行了。對這種情況，應當先在gview或VMD等程序里把分子朝向調正，從而能夠在Multiwfn里定義矩形盒子來恰好框住孔洞區域。
4 例：富勒烯
下面我們來計算富勒烯C60的內部空腔體積。如上文提到的，計算結果和空腔圖形也是受到計算時選擇的準分子密度閾值影響的。為了幫助選擇合適的閾值，我們可以先利用主功能5看看等值面圖形。C60的pdb文件可在這里下載C60.pdb。啟動Multiwfn，載入C60.pdb，然后依次輸入
5 //計算格點數據
1 //準分子密度
-10
0 //不留延展距離，而讓盒子緊貼體系邊緣
1 //低質量格點，對于粗略考察足夠了（注意此處主功能5里的格點質量的“低/中/高”和域分析模塊的這個概念不同。前者的“低/中/高”對應的是不同的總格點數，而后者對應的是不同格點間距）
-1 //觀看等值面
在界面里選擇Show data range可以把盒子范圍顯示出來。我們隨便考察幾個等值面數值，圖像分別如下
可見，準分子密度等值面為0.1的時候里外連通，此設定下沒法劃分出對應空腔的域。而0.01和0.001時里外都斷開了，因此都可以計算空腔體積，而不同設定下結果必定不同，等值面數值越小時空腔對應的域越小。取什么數值up to you，但一定要在文章中說清楚空腔是怎么定義的。
下面同上一節一樣計算域，先回到主界面，然后依次輸入
200
14 //域分析
2 //設定用于定義域的實空間函數
1 //準分子密度
3 //設定域的定義規則
<0.001 //將小于0.001 a.u.電子密度的區域作為域
1 //計算格點數據并產生域
-10 //修改延展距離
0 //延展距離設為0
1 //低質量格點
輸出信息如下：
根據前面0.001 a.u.下的準分子密度等值面圖，從屏幕輸出的信息上一眼就能看出5號域肯定是對應富勒烯的空腔的，因為體積7.08 ?^3比其它域都大得多。5號域的圖形如下。由于空腔接近圓形，我們可以推算出空腔半徑為1.19埃。
5 例：修飾的富勒烯
這是我隨便畫的一個體系
對這個體系，如果我們要計算富勒烯里的空腔體積，讓盒子緊貼著整個體系的邊緣定義就不合適了，因為由于修飾上去的基團很大，很多格點將會分布在無關區域，白浪費計算時間。對這種只需要考察某個局部區域空腔特征的情況，我們應當令盒子只包圍感興趣的區域。
在Multiwfn設定格點的界面里，提供了很多種方式定義格點。對于當前情況，使用圖形界面來調整盒子中心位置、盒子尺寸是最方便的，這對應10 Set box of grid data visually using a GUI window選項。進入之后會看到下面的界面。在其中可以拖動相應滑框來移動盒子中心X/Y/Z位置、調整盒子X/Y/Z尺寸，以及調節格點間距。當前狀況對應的格點數實時顯示在窗口右下角。把盒子設成下圖這樣就是比較合適的了，正好框住富勒烯部分。
各項都設好后，就可以選Return，程序馬上就按照當前的設定計算格點數據、產生域了，其它的就沒什么好說的了。