使用Multiwfn做Hirshfeld surface分析直觀展現分子晶體和復合物中的相互作用

使用Multiwfn做Hirshfeld surface分析直觀展現分子晶體和復合物中的相互作用

First release: 2024-Feb-16 Last update: 2024-Feb-17

0 前言

Hirshfeld surface分析（以下簡稱HS分析）是展現分子晶體中一個或多個分子與周圍分子間相互作用的方法，它也同樣可以用于展現孤立體系中特定片段與周圍原子間的相互作用。HS分析的原理容易理解，圖像比較直觀，在分子晶體的研究領域已經用得非常普遍。Multiwfn的主功能12中很早以前就已經實現了HS分析，如今已經被不少文章使用。本文將結合許多例子，專門詳細具體講解一下Multiwfn做HS分析的各方面細節、操作和技巧。

本文的讀者請務必使用2024-Feb-16及以后更新的Multiwfn版本，否則與本文所述情況會有很多不同。Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者請參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。使用Multiwfn做HS分析在寫文章時請按照Multiwfn啟動時的提示對程序進行恰當的引用。

很值得一提的是，筆者提出的IGMH方法也非常適合展現分子晶體中的分子間相互作用，HS分析與IGMH分析的展現形式有明顯區別且有一定程度的互補性，二者亦可以同時使用以提供更多視角。請閱讀《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）和《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）提到的綜述了解IGMH的相關知識，里面也專門有IGMH用于分子晶體的實際例子。

1 Hirshfeld surface分析的基本思想

在筆者講授的量子化學波函數分析與Multiwfn程序培訓班（http://www.keinsci.com/workshop/WFN_content.html）和Multiwfn手冊3.15.5節對HS分析的原理有很具體、詳細的講解，下文只是把原理的最關鍵的部分簡要介紹一下，對于正確做HS分析基本夠了。

HS分析最早由Spackman等人于Chem. Phys. Lett., 267, 215 (1997)提出，后來又得到了發展，Acta Cryst., B60, 627 (2004)和CrystEngComm, 11, 19 (2009)是其兩篇綜述文章。HS分析關鍵思想是對分子晶體中的特定分子構造出Hirshfeld surface，它相當于這個分子在分子環境中的表面，然后再將一些有特殊意義的實空間函數映射到這個表面上，由此可以對分子晶體中的分子間相互作用特征進行考察。

簡要說一下Hirshfeld surface是怎么定義的。Hirshfeld在Theoret. Chim. Acta (Berl.), 44, 129 (1977)最早提出了一種定義化學體系中原子空間的方式，它給每個原子定義了Hirshfeld權重函數來描述這個原子在三維空間中各個位置所占權重，數值從0到1平滑變化，0和1分別對應于此位置完全不屬于和完全屬于這個原子。每個位置所有原子的權重函數加和為1。這是一種典型的模糊式原子空間定義方式。具體定義細節見《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）的2.5節和Multiwfn手冊3.9.1節，看完了就會知道產生Hirshfeld權重需要有各個原子的坐標，以及體系中的各種元素的原子在孤立狀態下的球對稱化的電子密度。將一個分子中所有原子的Hirshfeld權重函數加和就定義了這個分子的權重，因此分子晶體中各個分子都有各自的權重函數。某個分子的Hirshfeld surface就對應于它的權重函數數值為0.5的等值面。例如下圖的曲線展現了某個平面上各個分子的Hirshfeld surface對應的輪廓。可見，Hirshfeld surface算是分子環境中各個分子與其它分子接觸面的一種定義方式。下文所說的“表面”都是指Hirshfeld surface。

現實當中構造Hirshfeld surface是利用我在J. Mol. Graph. Model., 38, 314 (2012)介紹的Marching Tetrahedron或類似的算法實現的。這個表面被描述為大量小三角形的集合，每個三角形由三個表面頂點構成。由于構造Hirshfeld surface對應的等值面用的算法和《使用Multiwfn的定量分子表面分析功能預測反應位點、分析分子間相互作用》（http://www.shanxitv.org/159）涉及的構造電子密度等值面對應的分子范德華表面的算法在本質上相同，所以Multiwfn中HS分析也是在主功能12（定量分子表面分析）中實現的。

實際中做HS分析的時候是選取一個（也允許是多個）感興趣的分子并對它構造Hirshfeld surface進行分析，這個分子在下文管它叫“中心分子”，它周圍的分子也可以稱為“環境分子”。

HS分析定義了一些映射到Hirshfeld surface上的三維實空間函數，常見的有：
? d_i：表面內部（interior）的原子（也即中心分子的原子）到當前點的最近距離
? d_e：表面外部（exterior）的原子（也即周圍分子的原子）到當前點的最近距離
? d_norm：歸一化的（normalized）分子間接觸距離。由d_i、d_e和原子范德華半徑計算出來。數值越小體現表面上此處附近的內、外原子有越近的接觸，暗示此處的相互作用越強
? shape index（形狀指數）：數值在[-1,1]范圍，越負說明此處表面越凹，反之越凸
? curvedness（曲度）：在[-4,0.4]范圍，-4對應表面此處完全平坦，越正越凸，0對應單位球面的曲度

d_norm經常用來對Hirshfeld surface進行著色來直觀展現分子環境中的分子間的相互作用，d_norm數值較小的區域對應于較強烈的分子間相互作用。筆者發現用電子密度來著色也很有意義，Hirshfeld surface上電子密度越大的地方體現相互作用越強，效果比用d_norm著色時明顯更好，色彩過渡更為平滑，物理意義也更強，也沒有范德華半徑選取的任意性。由于對大體系做量子化學計算產生比較精確的電子密度比較耗時，因此只需要用準分子近似的電子密度（promolecular density）對Hirshfeld surface著色就夠了，它直接由各個原子孤立狀態的電子密度疊加得到，耗時極低，我實測和使用量子化學計算的電子密度著色的效果差不多。

由于各種元素的孤立狀態的電子密度以及范德華半徑在Multiwfn中是內置的，故產生Hirshfeld surface以及計算以上提及的各種函數只需要原子坐標和元素信息就夠了。由于HS分析依賴的信息非常簡單，不牽扯基于波函數的計算，因此耗時極低，用起來也很方便。

HS分析還經常繪制指紋圖（fingerprint map），是把d_i和d_e作為散點圖的橫軸和縱軸，然后把Hirshfeld surface上的各個頂點根據其d_i和d_e的數值繪制在指紋圖上作為一個個小點。根據指紋圖上散點分布的位置可以對中心和周圍分子間的相互作用特征進行討論，后文會結合具體例子來講。

HS分析中還經常做局部接觸（local contact）分析。完整的HS分析描述了中心分子的所有原子和周圍分子的所有原子間的相互作用，而局部接觸分析可以指定在HS分析中只考慮中心分子的哪些原子和周圍分子的哪些原子的相互作用。比如可以了解在整個Hirshfeld surface中體現中心分子的氧和周圍分子的H之間的相互作用的區域的位置和面積。

Multiwfn還支持Becke surface分析，是我自己提出的概念。它和HS分析的唯一差別是用Becke權重函數而非Hirshfeld權重函數的0.5等值面來定義表面。二者實際效果差別不大，一般沒必要用Becke surface，但它的一個特殊好處是允許表面出現在電子密度為0的區域，此處沒法定義Hirshfeld權重并構造Hirshfeld surface，詳見本文的第5節。Becke權重函數的具體定義方式參見《密度泛函計算中的格點積分方法》（http://www.shanxitv.org/69），它基于原子坐標和原子共價半徑得到。出現在相互作用的原子間的Becke surface會離半徑較小的原子較近、離半徑較大的原子較遠，Hirshfeld surface也有這樣的特點，這是由于定義它的準分子密度分布特征所自然而然帶來的。

2 Multiwfn的HS分析的功能

Multiwfn中做HS分析需要提供含有原子信息的文件作為輸入文件，如.pdb、.xyz、.mwfn、.cif、.fch、.mol2、.gjf等等等等，詳見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。

HS分析大多研究的是分子晶體，cif是最常用的記錄晶體結構的格式。通常不能載入cif文件后上來就做HS分析，因為HS分析一般需要提供一個中心分子+環境分子的簇模型，這樣才能靠HS分析考察中心分子與環境分子的相互作用，而cif文件記錄的是晶胞里各個原子的坐標，分子往往是被截斷的，中心分子或環境分子一般都不完整。因此首先需要用《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）中介紹的自動挖團簇的功能構造簇模型。如果你研究的不是分子晶體的情況，就是比如分子二聚體中兩個分子間的相互作用，那么就不牽扯挖簇的過程了，直接提供含有二聚體坐標信息的文件當輸入文件就行了。

HS分析在Multiwfn主功能12實現。進入這個功能后，首先選1把定義表面的方式切換為Hirshfeld surface。此時被映射到表面的函數會自動改為電子密度（對于輸入文件沒有波函數信息的情況具體是指準分子電子密度），你也可以選2把被映射的函數改成其它的。之后選0就開始計算了，Hirshfeld surface會被構造出來，組成它的所有表面頂點上的被映射的函數值會被計算出來，并顯示出Hirshfeld surface的面積和包圍的體積。然后會看到后處理菜單，里面有豐富的選項，提示得都很清楚，Multiwfn手冊3.15.5節都有解釋。利用后處理菜單的選項，可以導出用于在VMD程序中繪制著色的Hirshfeld surface圖要用的.cub文件，還可以繪制指紋圖、做局部接觸分析，這些在后文的例子里都有體現。

后文的例子涉及到的各種文件可以在http://www.shanxitv.org/attach/701/file.zip中獲得。

3 用Multiwfn對NAOB晶體做Hirshfeld surface分析實例

這一節以(Z)-4-((2-nitrophenyl)amino)-4-oxobut-2-enoic acid (NAOB)晶體為例做HS分析。NAOB的分子結構如下，其晶體在DOI: 10.1007/s13738-023-02904-9里進行了研究，文章的補充材料里給了它的cif文件的信息，是本文文件包里的NAOB.cif。

3.0 準備工作：構造簇模型

NAOB.cif對應的晶體結構如下，可見連一個完整的分子都沒有，因此在進行HS分析之前，我們顯然得先構造出中心分子+周圍分子的團簇結構才行。

為了構造簇模型，啟動Multiwfn，然后載入NAOB.cif，之后輸入
300 //主功能300
7 //幾何相關操作
25 //構造“中心分子+臨近一圈分子”的團簇
1 //當前晶胞里1號原子所在的分子將被作為中心分子（由于NAOB晶體里所有分子都是等價的，所以當前隨便輸入一個原子序號即可）
[回車] //代表若一個周圍分子與中心分子間最近原子對距離小于這倆原子的范德華半徑和的1.2倍，則這個周圍分子就被整個納入團簇

Multiwfn瞬間就構造出了團簇，從屏幕上的提示可看到這個簇有375個原子，并且屏幕上還巨貼心地把中心分子中的原子序號給了出來，此例為1-10,17-19,22-25,205-210,214,215。把這個序號記下來，之后HS分析時要用到。

現在可以選當前菜單中的選項0看一眼新構造出的團簇是什么樣，如下所示，可見非常理想，確實是中心分子被周圍一層分子所圍繞（為了中心分子看得清楚，在Multiwfn圖形界面的菜單欄里選了Other settings - Set atom highlighting，然后輸入了前述的中心分子里的原子序號，使中心分子用青色高亮了）。

點圖形界面右上角的RETURN按鈕關閉圖形窗口，然后選擇-2 Output system to .pdb file并輸入NAOB_cluster.pdb，以將當前簇結構導出為當前目錄下的這個文件。這個pdb文件在本文的文件包里也提供了。

注：如《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）所強調的，由于X光衍射實驗一般難以確定氫原子的準確位置，因此原理上做HS分析之前最好先固定重原子而優化一下所有氫原子的位置。用免費高效的CP2K程序對晶體結構優化氫原子位置然后再用Multiwfn摳團簇，或是先摳團簇再用Gaussian等量子化學程序優化氫，都是可以的。

3.1 繪制電子密度著色的Hirshfeld surface圖

這一節演示繪制基于準分子近似的電子密度著色的Hirshfeld surface圖，這是HS分析最重要的一種圖。除了Multiwfn外還會用到非常流行的VMD可視化程序，可以在http://www.ks.uiuc.edu/Research/vmd/下載，使用筆者現在用的VMD 1.9.3版肯定沒問題，用其它版本不保證能按照本文的例子正常作圖。

啟動Multiwfn，載入上一節產生的NAOB_cluster.pdb，然后輸入
12 //定量分子表面分析
1 //選擇定義表面的方式
5 //Hirshfeld surface
1-10,17-19,22-25,205-210,214,215 //中心分子的原子序號
0 //開始分析

從屏幕上可以看到許多信息，以下兩條是值得注意的，第一個是Hirshfeld surface包圍的體積，相當于分子晶體中屬于這個中心分子的體積，第二個是Hirshfeld surface的面積
Volume: 1644.06256 Bohr^3 ( 243.62494 Angstrom^3)
Overall surface area: 859.01840 Bohr^2 ( 240.54965 Angstrom^2)

現在看到了后處理菜單。選擇-2在當前目錄下導出surf.cub，它是中心分子的Hirshfeld權重的格點數據。再選13，Multiwfn會計算被映射的函數的格點數據并導出為當前目錄下的mapfunc.cub。在這個界面里如果選擇-3也可以直接在Multiwfn里預覽Hirshfeld surface，但沒有著色效果。

將剛剛產生的surf.cub和mapfunc.cub以及Multiwfn自帶的examples\scripts\目錄下的hirsh_rho.vmd作圖腳本一起拷到VMD目錄下（即啟動VMD后在其文本窗口里運行pwd命令看到的目錄）。啟動VMD，在文本窗口里輸入source hirsh_rho.vmd來執行作圖腳本，然后就看到了下圖。

上圖顯示出了中心分子的Hirshfeld surface，并根據電子密度進行了著色。通過hirsh_rho.vmd腳本里的mol scaleminmax top 1 0.0 0.015這條命令可知，默認用的色彩刻度范圍的下限和上限分別為0和0.15。根據腳本里color scale method BWR這行命令可知，當前用的色彩變化是藍-白-紅，色彩變化示意圖在Graphics - Colors - Color Scale里可以看到。因此上圖中越紅的地方就是Hirshfeld surface上電子密度越大的地方，無疑對應于越強的相互作用。上圖中有一塊非常大的紅色，這對應于中心分子羧基的O-H作為氫鍵給體、周圍一個分子的O作為氫鍵受體形成的顯著的氫鍵作用區域。它的旁邊還有一小塊淡紅色區域，對應于中心分子的羧基氧作為氫鍵受體與周圍分子的C-H形成的弱氫鍵。上圖中還有一些發白的區域，在VMD里旋轉圖像仔細觀察的話可以看到對應的是比較遠距離的C-H...O-N很弱氫鍵，以及pi-pi堆積和普通色散吸引作用顯著的區域。這些特征區域都可以自行用powerpoint之類畫個箭頭標注在圖上便于讀者看清楚。上圖還有很多偏藍色的區域，這些地方電子密度接近0，因此不牽扯任何值得一提的分子間相互作用。特別藍的地方也往往對應于分子晶體中的孔洞區域，這種地方電子密度自然特別低，非常建議感興趣的讀者按照《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》（http://www.shanxitv.org/617）介紹的方法作圖考察。

要注意中心分子與各個方向的周圍分子都有相互作用，光是靠一張圖的話很難展現完整，因此文章中可以多給幾張圖展現不同視角的Hirshfeld surface圖。

下面再說一下怎么改進作圖效果。Hirshfeld surface圖的效果受到多方面影響：
(1)光源。可以通過VMD的Display菜單里的Lighting選項設置打開哪些光源。如果選Mouse - Move light，然后在圖形窗口中拖動，還可以移動特點光源的位置。
(2)材質。hirsh_rho.vmd默認對等值面用Translucent材質，可以自行在Graphics - Materials界面里對這個材質的具體定義進行調節。
(3)Graphics - Representation界面里的作圖設置。里面可以創建更多Rep，每個Rep都可以獨立設置顏色和材質，并且通過選擇語句可以定義各個rep顯示的原子，不懂選擇語句怎么寫的話參考《VMD里原子選擇語句的語法和例子》（http://www.shanxitv.org/504）。特別值得一提的是，當前的圖中每個分子都有一個獨立的residue編號并被分子中所有原子所共享，因此可以利用residue選擇特定分子。若想查詢某個分子的residue號，就選Mouse - Query，然后點擊這個分子上任意一個原子的正中央，在VMD的文本窗口中就能看到它的residue號了。

為了讓上面例子的圖像效果更好，我在VMD的Graphics - Representation里做了些修改：把第1個rep設為了residue 0專門用于顯示中心分子，用CPK方式顯示并把Bond Radius設為了0.5。點擊Create Rep按鈕新增一個Rep，選擇語句用residue 7 11 5使得三個與中心分子作用顯著的分子顯示出來，Material設Edgy，Drawing Method用Licorice，把Bond Radius設成0.1。點擊用來顯示Isosurface的那個Rep，再選擇Trajectory標簽頁，在里面把色彩刻度上限設為0.012（注意文本框只顯示兩位小數，輸入0.012后按回車會如實設為0.012），使得著色的色彩顯得更鮮明。在Display菜單里把所有四個光源都打開。最后選File - Render，選擇Tachyon(Internal)進行渲染。之后在圖像上展現特征作用的地方通過powerpoint進行標注，并且把Graphics - Color - Color Scale里顯示的色彩刻度條挪到圖上并適當拉伸、標記上色彩刻度上下限。最終得到下圖，可見對相互作用展現得非常直觀清楚。

如果要繪制d_norm著色的Hirshfeld surface圖，與上面的例子只有兩個差別
(1)進入主功能12并選擇用Hirshfeld surface方式定義分子表面后，選擇2 Select mapped function，再選d_norm。然后再開始分析
(2)作圖使用Multiwfn目錄下examples\scripts\里的hirsh_dnorm.vmd代替hirsh_rho.vmd
由于d_norm著色的圖明顯不如電子密度著色的圖色彩變化平滑，在很多地方有顏色的突越不好看，故這里就不多說了。

Multiwfn還可以給出Hirshfeld surface上被映射的函數（當前為電子密度）的極大點的位置和數值，便于定量對比討論。在Multiwfn后處理菜單中選8 Export all surface vertices and surface extrema as vtx.pqr and extrema.pqr，此時當前目錄下就出現了vtx.pqr和extrema.pqr，分別記錄了所有表面頂點和表面極值點的坐標和被映射的函數數值。如屏幕上的提示所示，extrema.pqr里碳和氧原子分別用于記錄表面極大點和極小點，此文件在本文的文件包里也提供了。按前文用hirsh_rho.vmd作圖后，將extrema.pqr載入VMD，將其顯示方式設為VDW并把Sphere Scale設為0.2，顏色用黃色，此時看到下圖，每個黃色小球都對應Hirshfeld surface上電子密度極大點位置。

若想查詢表面極值點處的電子密度數值，就在VMD里選Mouse - Query，然后點擊要考察的黃球的正中心，文本窗口就出現了它的index號，從0開始記。上圖箭頭所指的那個對應O-H...O作用的極值點的index為9，相當于從1開始記的編號為10。打開extrema.pqr，找到對應10號碳的下面這一行，倒數第3列的數值0.04228468就是此處的電子密度了，單位為a.u.。
HETATM 10 C MOL A 1 -1.214 3.606 5.014 0.04228468 1.0000 C
以類似的方式查詢旁邊那個C-H...O氫鍵對應的Hirshfeld surface上的電子密度極大點，數值為0.01093053，可見作用顯著弱于O-H...O。

3.2 繪制指紋圖

這一節繪制HS分析中很常見的指紋圖。先按上一節的過程做HS分析并進入到主功能12的后處理菜單，然后輸入
20 //指紋圖與局部接觸分析
0 //開始分析
計算很快就完成了，然后看到一個新的后處理菜單，里面的選項不言自明。直接選0就會在屏幕上顯示完整的指紋圖，選1就會把指紋圖保存為當前目錄下的pdf文件（程序用pdf格式是因為它可以無損縮放、線條平滑），此pdf文件已經提供在了本文的文件包里（NAOB_HS.pdf）。當前的圖像如下

這個指紋圖里每一個小點都是一個Hirshfeld surface上的頂點，并且根據點在指紋圖上的分布密度進行了著色，越密的地方顏色越黃，越稀疏的地方顏色越紫。Multiwfn自動根據指紋圖上的最大分布密度設置色彩變化范圍的上限。Multiwfn畫指紋圖用的這種色彩過度方式稱為viridis，用Google搜圖功能一搜viridis就能找到相應的色條。

上圖在左下角有兩個顯著的尖兒（spike。兒化音明確體現這是個名詞），是圖像的特征區域，相當于“指紋”。靠左的那個尖兒的頂端大約是d_i≈0.7、d_e≈1.1位置，此處d_e顯著大于d_i，這是此體系具有氫鍵給體特征的體現（如上一節所示此體系的羧基氫確實是氫鍵給體）。為什么這個尖兒的存在能體現此分子存在氫鍵給體？因為這說明在Hirshfeld surface的這個區域，中心分子的原子到表面的最近距離（d_i）顯著小于周圍分子的原子到表面的最近距離（d_e），只有當氫鍵跨越這個區域，半徑很小的氫在表面內、半徑較大的重原子在表面外的時候才會出現這種情況，如下圖所示。

顯然不難理解前面的指紋圖中靠下的那個在d_e≈0.7、d_i≈1.1的尖兒體現的是中心分子作為氫鍵受體與周圍分子形成的氫鍵跨越了Hirshfeld surface。所以指紋圖說明當前研究的NAOB分子同時作為氫鍵給體和氫鍵受體。

指紋圖還蘊含了更多信息，但當前的指紋圖的中間區域密密麻麻一片，很難直接討論，進一步考察就需要利用下一節的局部接觸分析了。

順帶一提，如果你想改變指紋圖上的點的密度，可以在主功能12里做分析之前先選擇3 Spacing of grid points for generating molecular surface修改格點間距，格點間距越小則產生的Hirshfeld surface上的點就越多，指紋圖也會越密。

3.3 局部接觸分析與局部指紋圖

這一節演示一下怎么對中心分子特定的原子與周圍原子特定的原子進行局部接觸分析并繪制與之對應的指紋圖。作為例子，這里考察中心分子的氧與周圍分子的氫的局部接觸情況。

按照上一節說的進入到主功能12的后處理菜單中的20 Fingerprint plot and local contact analyses選項后，輸入以下內容
1 //設置內部原子范圍。默認是所有中心分子的原子都包括
[回車] //代表對原子序號不做限制
O //必須是氧元素
現在中心分子的所有氧原子都納入到了要考慮的內部原子范圍了。然后再輸入
2 //設置外部原子范圍。默認是所有周圍分子的原子都包括
[回車] //代表對原子序號不做限制
H //必須是氫元素
現在周圍分子的所有氫原子都納入到了外部原子要考慮的范圍了。

選擇0開始分析，算完后屏幕上顯示以下信息，告訴你當前考察的這種接觸面積是61.5 Angstrom^2，占Hirshfeld surface總面積的25.57%。
The area of the local contact surface is    61.504 Angstrom^2
The area of the total contact surface is   240.550 Angstrom^2
The local surface occupies   25.57% of the total surface

然后在后處理菜單選擇繪制指紋圖，看到下圖（對應本文文件包里的NAOB_HS_O-H.pdf）。此圖中只有對應于當前考察的局部接觸面上的頂點才被繪制為彩色，可見指紋圖中d_e≈0.7、d_i≈1.1的尖兒確實對應于內O與外H的接觸。

這種局部接觸還可以在立體結構圖上體現。在當前界面里選擇4 Export surface points to .pqr file in current folder，然后當前目錄下就產生了finger.pqr和finger_all.pqr，它們都提供在了本文的文件包里。前者記錄的是當前指定的局部接觸表面上的頂點，后者記錄的是完整的Hirshfeld surface上的頂點。它們可以用文本編輯器打開，可以看到每個表面頂點在pqr文件里用一個碳原子表示，倒數第三列記錄的是Charge屬性，當前被用來記錄表面頂點上被映射的函數數值，對當前來說就是準分子電子密度，單位為a.u.。

為了得到同時展現體系結構和局部接觸表面的圖，現在將NAOB_cluster.pdb載入VMD并恰當設置顯示方式，再把finger.pqr載入VMD，在Graphics - Representation里將它的Drawing Method設為Points并恰當設置Size，Coloring Method設為Charge（即根據Charge屬性著色），在Trajectory標簽頁里把色彩刻度設成與之前繪制等值面圖用的相同的0到0.012。確保Display - Rendermode已經設為了GLSL。在Graphics - Colors - Color Scale里把色彩刻度設為與之前相同的BWR，現在看到的圖如下，確實這些局部表面只對應中心分子的氧和周圍分子的氫的接觸。

也可以先用hirsh_rho.vmd照常繪制完整的Hirshfeld surface圖，然后再載入finger.pqr，把Drawing Method設為VDW并把Sphere Scale設為0.1，并且對Charge屬性用0到0.012色彩刻度范圍著色，此時看到的圖像如下，局部接觸部分以小球形式著重展現了。

下面再做個演示，對中心分子的芳環和周圍分子的芳環之間考察局部接觸，展現它們之間的pi-pi堆積。在GaussView里，或者在VMD里使用我在《在VMD中顯示原子序號的方法》（http://www.shanxitv.org/197）中提供的腳本，都可以在它們載入NAOB_cluster.pdb后顯示出原子序號，可以看到中心分子的芳環的原子序號是7-9,205,207,209。與中心分子芳環較近因而有pi-pi堆積作用的周圍分子的芳環有兩個，原子序號為11,13,15,112,114,116,201-203,356-358。因此，進入前述的Multiwfn的20 Fingerprint plot and local contact analyses選項后，輸入
1 //設置內部原子范圍
7-9,205,207,209 //原子序號范圍
[回車] //對元素不做限制
2 //設置外部原子范圍
11,13,15,112,114,116,201-203,356-358 //原子序號范圍
[回車] //對元素不做限制
0 //開始分析
4 //導出finger.pqr和finger_all.pqr

然后按照上一節的方法，用hirsh_rho.vmd繪制出Hirshfeld surface，再將finger.pqr載入VMD并恰當設置顯示方式。顯示體系結構的rep的選擇語句寫residue 0 1 11，其中residue 0是中心分子，residue 1和11對應與它有pi-pi堆積作用的上、下兩個分子。然后在Display - Material里把Translucent材質的Opacity改為0.5以降低透明度，并打開Angle-Modulated Transparency選項使得等值面立體感更強一些。之后看到的圖如下，可見小圓球把中心分子芳環上下兩側的pi-pi堆積對應的接觸區域都清晰展示了出來。

在Multiwfn里顯示出相應的局部接觸表面的指紋圖，如下所示。由于碳的原子半徑不小，因此pi-pi堆積對應的局部表面與表面內和表面外的碳原子都有一定距離，而碳原子間距離若太遠也不會有pi-pi堆積效應，因此這些表面頂點在指紋圖中的位置是d_i和d_e都不大不小的區域。

3.4 統計不同元素間的接觸面積

如果想得到中心分子與周圍分子的每一對元素之間的局部接觸面積，雖然通過前面演示的局部接觸分析可以手動實現，但有多少種元素組合就得操作多少次，很麻煩。因此Multiwfn提供了一次性完成的功能。進入前述的Multiwfn的20 Fingerprint plot and local contact analyses選項后，不需要定義內部外部原子，直接選擇選項3 Calculate contact area between different elements，就可以馬上得到以下信息

Inside element, outside element, their contact area (Angstrom^2) and percentage (%)
H-H        47.330      19.676
H-C        10.423       4.333
H-N         0.521       0.216
H-O        56.618      23.537
C-H        12.821       5.330
C-C        22.304       9.272
C-N         2.614       1.087
C-O         4.552       1.892
N-H         0.876       0.364
N-C         2.546       1.059
N-N         0.464       0.193
N-O         1.910       0.794
O-H        61.504      25.568
O-C         4.280       1.779
O-N         1.816       0.755
O-O         9.970       4.145

The same as above, but do not distinguish inside and outside elements
H-H              47.330      19.676
H-C/C-H          23.244       9.663
H-N/N-H           1.396       0.581
H-O/O-H         118.121      49.105
C-C              22.304       9.272
C-N/N-C           5.160       2.145
C-O/O-C           8.831       3.671
N-N               0.464       0.193
N-O/O-N           3.727       1.549
O-O               9.970       4.145

Area of total contact surface is 240.550 Angstrom^2

可見以上信息包含了每一對元素的分析結果。例如O-H 61.504 25.568這一行代表中心分子的氧元素與周圍分子的氫元素之間的局部接觸面積是61.504 Angstrom^2，占總面積的25.568%，這和3.4節我們專門算的結果完全一致。而后面輸出的比如H-O/O-H 118.121 49.105代表不管O和H誰在表面內、誰在表面外，兩種元素之間的接觸面積總共是118.121 Angstrom^2，占總面積的49.105%。

為了看起來更直觀，可以把上面的三列形式的數據拷到txt文件里，然后拖到Origin中導入，再繪制餅形圖，如下所示，一目了然。圖中只有占比相對較大的幾種接觸直接標注了數值。

3.5 Hirshfeld surface等值面的截斷問題

Multiwfn在構造Hirshfeld surface之前，需要先對一個矩形區域（稱為“盒子”）中均勻分布的格點計算你設定的體系片段的Hirshfeld權重。盒子范圍是對你定義的片段往各個方向延展一定距離來確定的，延展距離是位于最邊界的原子的范德華半徑乘以一個系數得到的。默認的系數值是1.7，一般來說夠大了，但如果Hirshfeld surface延伸到距離當前片段較遠的地方，導致超過了盒子范圍，則Hirshfeld surface就會在相應地方被截斷，等值面在那個地方看起來就會有窟窿。此時如果你想要讓等值面完整、完全封閉，顯然就需要增大系數。這一節拿富勒烯晶體舉個例子。

本文文件包里的C60.cif是C60富勒烯晶體結構文件。按照3.0和3.1節的做法摳團簇、繪制電子密度著色的Hirshfeld surface，會得到下圖。這里我在Graphics - Representation里選擇顯示等值面的那個Rep后在界面右下角Show旁邊的下拉框里選擇了Box+Isosurface，這樣除了等值面外，盒子范圍還會用細線同時顯示出來。由下圖可見，等值面上有個難看的窟窿，因為6個富勒烯之間有孔洞區域，Hirshfeld surface實際上會延伸到這里，但當前被尺寸有限的盒子截斷了。

當你發現存在這樣的窟窿，想通過增大延展距離來避免，就應當從主功能12的后處理菜單返回之前的菜單（即剛進入主功能12看到的菜單），然后輸入
4 //高級選項
1 //設置確定延展距離用的范德華半徑的倍數
2.3 //設一個比原本更大的值，數值可以反復嘗試。設得越大，盒子就越大，格點數就越多，計算耗時就越高、產生的cub文件尺寸就越大。由于之前的等值面只有一個小窟窿，所以盒子再大一點就夠了，因此當前嘗試的2.3比之前默認的1.7沒大太多
0 //返回
0 //開始計算

之后按照之前的步驟繪制Hirshfeld surface圖，如下所示，窟窿已經沒了。當前的Hirshfeld surface是個正十二面體，每個面都有塊白色且中間偏紅的區域，清楚地展現出中間的富勒烯和周圍富勒烯在此處有顯著的pi-pi堆積作用，由于這些區域離碳原子不遠因此電子密度不是很小。而表面上藍色部分都是晶體中的縫隙、孔洞區域，電子密度極低。

順帶一提，此體系的指紋圖如下所示。繪圖前用Multiwfn菜單里的選項3 Set range of axes適當加大了上限。可見散點的分布區域很窄，各處d_i和d_e數值都差不多大，涵蓋從中等到較大范圍。這體現出表面頂點基本都是位于中間分子和周圍分子之間的正中央區域，并且中間分子和周圍分子有的地方挨得近（d_i和d_e為中等大小），有的地方離得遠（d_i和d_e較大）

4 將Hirshfeld surface分析用于分子復合物

這一節示例HS分析用于展現孤立體系中的片段間相互作用。片段是指一批原子的集合，可以根據需要自由定義，比如既可以是分子復合物中的一個或多個分子，也可以是一個分子中的一個基團，也可以是一個過渡金屬配合物中的一個或多個配體，等等。對于孤立體系，由于感興趣的片段通常不是像晶體環境一樣被其它原子所完整包圍的，所以相應的Hirshfeld surface通常是開放的而不是完全封閉的。

筆者之前做了大量的和18碳環及其衍生物有關的研究工作，匯總見http://www.shanxitv.org/carbon_ring.html。其中研究了具有雙環特征的OPP分子與18碳環的結合問題，成果介紹見《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）和《理論設計新穎的基于18碳環構成的雙馬達超分子體系》（http://www.shanxitv.org/684）。研究中使用了IGMH方法直觀展現了OPP與18碳環之間的弱相互作用，此例通過繪制電子密度著色的Hirshfeld surface圖也來展現一下。

上述研究中使用DFT優化出來的OPP結合一個18碳環的結構文件是本文文件包里的C18-OPP.pdb。將之載入Multiwfn，然后進入主功能0，會發現此結構文件中18碳環是斜著的，而由于HS分析用的盒子的邊框總是平行于笛卡爾軸的，因此直接對它產生Hirshfeld surface的話表面開口的地方也將是斜著的，看著很惡心。因此此例做HS分析之前還需要先用《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）里介紹的功能令18碳環平行于XY平面，如下所示。

啟動Multiwfn，載入C18-OPP.pdb，然后輸入
300 //其它功能（Part 3）
7 //對當前體系做幾何操作
11 //令一批原子擬合的平面平行于某個笛卡爾平面
225-242 //18碳環的原子序號
1 //平行于XY平面
現在可以選0在圖形界面觀看當前結構，會發現確實18碳環已經在XY平面上了。點RETURN關閉圖形窗口，然后接著輸入
-10 //返回
0 //返回到主菜單
12 //定量分子表面分析
1 //設置表面的定義
5 //Hirshfeld surface
225-242 //18碳環的原子序號
4 //高級選項
1 //設置定義盒子延展距離用的范德華半徑的倍數
2.3 //比默認值更大，實測不這樣的話會導致不好看的截斷
0 //返回
0 //開始分析
-2 //導出surf.cub
13 //導出mapfunc.cub

現在基于surf.cub和mapfunc.cub用hirsh_rho.vmd照常繪制Hirshfeld surface圖。這里我把VMD的四個光源都打開了，18碳環用CPK方式顯示，雙環OPP分子用Edgy材質。為了著色效果更鮮明，色彩刻度上限設為了比默認更小的0.01。此時圖像如下所示，可見白色和偏紅的條狀區域把18碳環與OPP的pi-pi作用最強烈的區域展現得挺清楚，此處電子密度比周圍相對更大。由于OPP的大環不是像18碳環一樣理想的圓形，形狀沒有完美匹配，因此它們之間的相互作用明顯不是處處均勻的。通過Hirshfeld surface的顏色可見，在大環末端，特別是兩個大環連接區域，大環與18碳環的pi-pi作用都很弱、電子密度很低。當前的Hirshfeld surface在18碳環上、下方都是開口的，因為在這里被盒子邊界截斷了，顯然這樣的截斷是理應有的，因為上下沒有其它原子了，原理上就不可能是封閉的。

下面這張圖是Comput. Biol. Chem., 101, 107786 (2022)中用Multiwfn+VMD畫的Hirshfeld surface圖，展現了配體和附近的蛋白質的氨基殘基間的相互作用情況，效果不錯，作用位點能看得比較清楚。

值得一提的是HS分析不僅限于用于考察弱相互作用，讓Hirshfeld surface跨越化學鍵也完全可以。在《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）介紹的我的綜述文章中給出了一個九并苯共價二聚體的例子，這里我也繪制了其電子密度著色的Hirshfeld surface，其中一個九并苯被定義為了片段，色彩刻度用0-0.01，圖像如下所示。可見Hirshfeld surface正好嚴格平行于兩個九并苯、精確處于二者之間。在兩個九并苯之間形成共價鍵的區域及附近電子密度很大，顏色都為紅色，而其它白色區域體現了九并苯之間的明顯的pi-pi堆積作用。

5 Becke surface分析一例：DNA

距離原子非常遠的地方電子密度精確為0，這樣的地方是沒法計算Hirshfeld權重的，Hirshfeld surface也沒法在這個區域出現，這是HS分析的一個缺點。如果極個別情況下被研究的片段與其它部分的接觸面就是會牽扯到距離所有原子都很遠的區域，那么這個時候就必須用我提出的Becke surface代替Hirshfeld surface了，因為在電子密度為0的地方也是可以計算片段的Becke權重的。除了這種極個別情況外沒有用Becke surface的必要，圖像效果不會更好，而且對大體系耗時更高。

這一節就以DNA雙螺旋結構為例繪制做Becke surface分析。Multiwfn自帶的examples目錄下的DNA.pdb是一個DNA片段的結構文件，其中一條鏈的原子序號是1-319，它將被定義為片段。

啟動Multiwfn，載入DNA.pdb，然后輸入
12 //定量分子表面分析
1 //選擇表面的定義
6 //Becke surface
1-319 //第一個片段的原子序號
4 //高級選項
1 //設置定義盒子延展距離用的范德華半徑的倍數
0 //倍數設為0，使得盒子緊貼著邊界原子。因為感興趣的相互作用區域是在兩條鏈間，這樣減小盒子尺寸可以避免Becke surface過大、延伸到不感興趣的區域去
0 //返回
0 //開始計算
在筆者的i9-13980HX筆記本上，用16核并行，花了15分鐘。嫌慢的話可以用核很多的服務器。對于預覽目的，也可以把格點間距設大到比如0.5，耗時只有原本的約八分之一，但此時圖像會明顯更為粗糙。

導出mapfunc.cub和surf.cub并照常用hirsh_rho.vmd繪圖，得到的圖像如下所示，可見Becke surface把兩條DNA鏈間的接觸面非常直觀地展現了出來。在每一層的兩個堿基之間都有兩塊紅色區域，體現出兩個典型氫鍵的存在導致相應地方電子密度相對較大。仔細看的話，會發現有的兩塊紅色區域旁邊還有一小塊白色區域，這對應于C-H...O很弱的氫鍵。Becke surface的其它區域都是藍色，說明兩條鏈在其它地方并沒有值得一提的相互作用（注意不能描述為“沒有相互作用”，用詞必須嚴謹。重原子間的色散作用能即便到了10埃也沒完全衰減到0，兩條鏈的骨架之間的色散作用對相互作用能的貢獻不可完全忽略，只不過強度屬于“不值得專門一提”的程度）

按照HS分析的做法也可以顯示Becke surface的指紋圖，如下所示，可見左下方存在兩個尖兒，體現出每條DNA鏈既作為氫鍵給體也作為氫鍵受體和另一條鏈形成了氫鍵，和HS分析能展現的信息基本一致。但尖兒的具體位置不可能十分一致，畢竟Hirshfeld和Becke權重函數的定義方式差異顯著。

6 總結

本文介紹了HS分析的基本思想，并通過大量例子非常詳細介紹了Multiwfn做HS分析的方法和很多要點。可見HS分析用起來靈活方便，可以較好地直觀展現化學體系的特定片段與周圍其它原子之間的相互作用。同樣的目的用Multiwfn做基于電子波函數的IGMH分析也同樣可以達到，而且IGMH原理更為嚴格，通過sign(lambda2)rho函數對等值面著色還能區分相互作用類型，還能給出各個原子的貢獻量。HS分析可以作為IGMH分析的展現形式的補充，并且由于HS分析只依賴于原子坐標而且計算量很低，因此可以快速地用于很大體系，也不需要事先做量子化學計算產生電子波函數。注：只有原子坐標信息時也可以用Multiwfn做IGM分析，耗時也極低，但圖像效果明顯不及IGMH，參見《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）。

最后提醒一下，使用Multiwfn做HS分析或Becke surface分析時請按Multiwfn啟動時的提示恰當引用Multiwfn原文，對于給別人代算的目的也需要明確告知對方這一點。