注：本文有對應的演示視頻，必看！見https://www.bilibili.com/video/av33733749/

凡是按照本文的詳細說明還繪制不出來的人，一定要精確嚴格重復視頻里的操作！應先把視頻里的示例體系重復一遍，連這個例子都重復不出來者注意看本文文末的文字

使用Multiwfn+VMD快速地繪制靜電勢著色的分子范德華表面圖和分子間穿透圖

文/Sobereva @北京科音

First release: 2018-Oct-13  Last update: 2023-Jul-2

1 前言

靜電勢對于考察分子間靜電相互作用、預測反應位點、預測分子性質等方面有重要意義，被廣為使用，相關學習資料和筆者之前寫過的全部跟靜電勢有關的博文見《靜電勢與平均局部離子化能綜述合集》（http://bbs.keinsci.com/thread-219-1-1.html），免費的波函數分析程序Multiwfn（主頁為http://sobereva.com/multiwfn）是分析分子靜電勢最靈活、最強大的工具。很久以前，筆者寫過一篇文章《使用Multiwfn結合VMD分析和繪制分子表面靜電勢分布》（http://sobereva.com/196，以下簡稱“196”）介紹了怎么將Multiwfn和免費的可視化程序VMD相結合繪制靜電勢著色的分子范德華表面圖以及分子間范德華表面穿透圖，但是筆者在大量網上答疑中感到很多接受能力較差的Multiwfn用戶還是不能很好按照博文里詳細講述的步驟做出圖來，另一方面那篇博文的過程涉及很多手動操作過程，對于手慢的用戶繪圖效率較低。考慮到這些，筆者這次專門寫一篇文章，介紹通過利用批處理文件和繪圖腳本極其便利地繪制上述圖形，將步驟最大程度簡化。相信此文會令讀者深感Multiwfn+VMD繪圖又快又便利，初學者也能很容易地獲得非常漂亮的圖像。此文并不能完全取代之前那篇博文，那篇文章里對細節講述很詳細，如果能把那篇文章完整閱讀、深刻領會，讀者就可以根據實際要求恰當調節選項獲得更好效果。

本文涉及兩種展現分子表面靜電勢分布的做法，第一種是196文中已經詳細交代的，用Multiwfn的主功能12對靜電勢做定量分子表面分析導出表面頂點文件vtx.pdb，載入到VMD中顯示出來，并根據vtx.pdb的B因子字段數據（對應靜電勢）來著色。第二種是類似《基于Multiwfn產生的cube文件在VMD和GaussView中繪制填色等值面圖的方法》（http://sobereva.com/402）文中的做法，用Multiwfn產生電子密度和靜電勢的cube文件，載入到VMD里，將靜電勢的數據以不同顏色映射到電子密度的等值面上。本文說的范德華表面是指Bader的定義，即把電子密度0.001 a.u.等值面作為范德華表面。

讀者請務必使用Multiwfn官網上的最新版本。不了解Multiwfn的話參看入門貼《Multiwfn入門tips》（http://sobereva.com/167）和《Multiwfn FAQ》（http://sobereva.com/452）。筆者使用的VMD是1.9.3版，可在http://www.ks.uiuc.edu/Research/vmd/免費下載。這里筆者假定用戶使用的是Windows系統。

本文的做法絕對不僅限于Gaussian用戶，GAMESS-US、ORCA、Molpro等絕大部分主流程序的用戶都可以按照本文的做法繪制，只要提供帶有波函數信息的文件即可，參考《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://sobereva.com/379）。

2 本文涉及的文件和準備工作

此文用到的所有輸入文件都可以在此下載：http://sobereva.com/attach/443/file.rar。本文用到的.bat、.txt、.vmd文件在目前最新版Multiwfn的examples\drawESP目錄下都可以找到，這些文件在這里介紹一下。

以下三個Windows下的批處理文件，在繪圖前讀者必須將它們放到Multiwfn所在目錄下。
ESPiso.bat：用于產生電子密度和靜電勢的cube文件，即density[序號].cub和ESP[序號].cub
ESPpt.bat：用于產生分子結構和表面頂點的pdb文件，即mol[序號].pdb和vtx[序號].pdb
ESPext.bat：在ESPpt.bat基礎上，還會額外產生分子表面靜電勢極值點的pdb文件（surfanalysis.pdb）

ESPiso.txt、ESPpt.txt、ESPext.txt分別被上述三個.bat文件所利用，里面每一行記錄了Multiwfn中每一步要輸入的命令，在繪圖前讀者必須將它們放到Multiwfn所在目錄下。

運行以上批處理文件調用Multiwfn計算前，讀者要將被計算的分子的fch文件改名為1.fch放到當前目錄下。雙擊批處理文件時，就會自動調用Multiwfn對當前目錄下的1.fch文件進行計算并得到繪圖所需的文件，并且會自動把結果文件拷到VMD目錄下。在使用批處理文件前，用戶必須自己編輯.bat文件，把里面默認的VMD目錄D:\study\VMD193改成自己機子的實際VMD目錄，如果VMD的路徑里含有空格，則路徑必須通過雙引號擴住，如"D:\study\VMD 193"。

如果當前目錄下還有2.fch、3.fch、4.fch，那么它們也會被這些批處理文件計算并產生相應的結果文件并自動拷到VMD目錄下，這用于繪制范德華表面穿透圖。

ESPiso.bat批處理腳本的運作原理在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://sobereva.com/612）有非常詳細易懂的介紹，強烈建議感興趣的讀者看此文。

以下三個.vmd文件都是VMD的繪圖腳本，在使用VMD繪圖前讀者必須將之拷到VMD目錄下。
ESPiso.vmd和ESPiso2.vmd：其中ESPiso2.vmd用于讀取VMD目錄下density1.cub、ESP1.cub、density2.cub、ESP2.cub...繪制成靜電勢填色等值面圖。讀者可以自行編輯這個文件修改里面默認的色彩刻度下限和上限，以及設置文件序號最多讀取到多少，文件中的#開頭的注釋行都寫得非常明白。ESPiso.vmd和ESPiso2.vmd的差異在于前者只是用于繪制單個分子的，而且一些作圖設定不一樣。
ESPpt.vmd和ESPpt2.vmd：其中ESPpt2.vmd用于讀取VMD目錄下mol1.pdb, vtx1.pdb, mol2.pdb, vtx2.pdb...繪制靜電勢著色的分子表面頂點圖。讀者可以自行編輯這個文件修改里面默認的色彩刻度下限和上限、點的尺寸、文件序號最多考慮到多少。ESPpt.vmd和ESPpt2.vmd的差異在于前者只是用于繪制單個分子的，而且一些作圖設定不一樣。
ESPext.vmd：用于讀取VMD目錄下的surfanalysis.pdb，將其中記錄的分子表面靜電勢極值點繪制成圓球，里面默認用的顏色都可以自己改

最后，為了繪圖時敲的命令更簡短，讀者應當用文本編輯器打開VMD目錄下的vmd.rc文件（對Windows版VMD而言），在最后添加以下內容：
proc iso {} {source ESPiso.vmd}
proc iso2 {} {source ESPiso2.vmd}
proc pt {} {source ESPpt.vmd}
proc pt2 {} {source ESPpt2.vmd}
proc ext {} {source ESPext.vmd}
之后在VMD的文本窗口里只要輸入比如pt，就等價于輸入source ESPpt.vmd來執行這個作圖腳本了。注意運行這些命令前用戶在VMD里最好沒做過任何其它操作，否則應該重啟VMD，否則可能看不到預期的效果。

總結一下，繪圖之前用戶要做的事包括：
(1)把examples\drawESP目錄下的.txt和.bat文件拷到Multiwfn可執行文件所在目錄
(2)把.bat文件里的VMD路徑改成實際路徑
(3)把examples\drawESP目錄下的.vmd拷到VMD路徑下
(4)在vmd.rc里加入上述proc語句

按照以上說明配置后，就可以開始做下面的例子了。如果沒成功重復出文中的圖，說明要么沒一個字一個字看上面的配置方式說明，要么沒一個字一個字看下面的操作步驟。如果繪制某個體系的時候VMD里顯示出了之前繪制過的其它體系，應當自行把VMD目錄下相應的其它體系對應的pdb或cube文件刪除。如果下文的透明的等值面顯示得很詭異很難看，可能是顯卡驅動和VMD不兼容而沒能打開GLSL，嘗試更新驅動或者用其它機子。

3 繪制單個分子的分子表面靜電勢圖例子

下面以乙酰胺為例，分別繪制兩種不同風格的分子表面靜電勢圖，二者是并列關系，你覺得顯示效果哪種好就用哪種。將Multiwfn的examples目錄下的CH3CONH2.fch拷到Multiwfn可執行文件所在目錄，改名為1.fch。

我們先通過表面頂點著色方式展現分子范德華表面靜電勢分布。雙擊ESPpt.bat，很快就運算完畢。然后啟動VMD，在文本窗口輸入pt，立刻看到下圖

ESPpt.vmd腳本里默認色彩刻度是-50kcal/mol ~ 50 kcal/mol，對應顏色按照“藍-白-紅”過渡。如果想在VMD的界面里直接修改顯示設定，諸如點的尺寸、色彩刻度范圍等，見196一文。

接下來我們通過電子密度等值面著色展現分子范德華表面靜電勢分布。雙擊ESPiso.bat，也是很快就算完了。啟動VMD（如果之前VMD已經打開了，關了它再重啟之），然后在VMD文本窗口直接輸入iso即可（在此之前絕對不要先運行pt命令！），看到的圖像將很接近于下圖。為了顯示效果更好，筆者此處使用了VMD自帶的Tachyon渲染器渲染，步驟在《用Multiwfn+VMD做RDG分析時的一些要點和常見問題》（http://sobereva.com/291）文中說明了，得到的圖像比在圖形窗口直接看到的更有光澤感而且有抗鋸齒效果。

最后，再看怎么把范德華表面靜電勢極值點位置顯示出來。雙擊ESPext.bat開始計算，它所干的事和ESPpt.bat一樣，只不過還額外產生一個記錄極值點的surfanalysis.pdb文件。按照如上任意一種方法顯示出靜電勢圖后，在文本窗口輸入ext即會看到極值點也顯示出來了。黃球對應靜電勢極大點，青球對應靜電勢極小點。以下左圖是先輸入pt再輸入ext的效果，右圖是先輸入iso再輸入ext然后用Tachyon渲染后的效果（之后在Graphics - Representation里把Selected Molecule切換為density1.cub，把Material選項材質改為Transparent，此時背面的極值點才能看到。如果你想默認就用這種材質，把ESPiso.vmd里的$id EdgyGlass替換為$id Transparent即可）

表面靜電勢極值點上的靜電勢的具體數值可以通過ps來標出來，在196文中詳細說了。基本過程是先點擊VMD的OpenGL圖形窗口將之激活，點擊鍵盤上的0進入VMD的查詢模式，然后點擊一個表面極值點對應的圓球（必須點正中心。如果目前同時顯示了表面頂點則先關閉顯示表面頂點免得妨礙選取），文本窗口里如果提示比如index 1，那么把這個數加上1就是它在surfanalysis.pdb中的實際序號，要加1是因為VMD的index編號是從0開始計的。要查看index 1的表面極值點對應的數值，就把surfanalysis.pdb里開頭的注釋行內容刪掉，找到第2行，其倒數第2列的數值就是靜電勢的數值了。對當前體系，第2行內容如下
HETATM    2  C   MOL A   1      -0.426  -2.953   0.281  1.00 44.12          C
故44.12 kcal/mol就是index 1這個點的靜電勢值。注意對于帶電體系，由于kcal/mol單位下的靜電勢數值可能超過了pdb格式的B因子這一列的可記錄范圍，此時Multiwfn自動會改用eV為單位記錄，看一眼pdb文件開頭提示用的是什么單位就知道了。順帶一提，在surfanalysis.pdb中靜電勢極大點和極小點分別用C和O原子表示。
還有一種更方便的查詢極值點靜電勢數值的方法：比如某個極值點是index 3，就在VMD窗口中輸入[atomselect top "index 3"] get beta即可顯示出數值。

4 繪制靜電勢著色的分子間范德華表面穿透圖

繪制這種圖實際上就是把每個單體的靜電勢著色的分子表面圖顯示到一起，從圖中可以一目了然看到哪些地方范德華表面被互相穿透，由此可以討論相互作用強度，還可以通過顏色討論靜電相互作用特征（一般是呈現靜電勢正負互補的，在《靜電效應主導了氫氣、氮氣二聚體的構型》http://sobereva.com/209一文有詳細討論）。在產生輸入文件過程中一定要記住，只需要對復合物進行優化，而單體的坐標必須直接從優化后的復合物的坐標里直接摳出來（比如可以用gview打開優化后的輸出文件，把其中某個單體復制到一個新窗口里，保存成gjf文件），之后對各個單體做單點計算任務得到含有波函數信息的文件即可。對于Gaussian，一定要記得計算單體的時候要加nosymm關鍵詞，否則由于Gaussian自動會把分子調整到標準朝向，導致坐標和復合物不對應。詳見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》（http://sobereva.com/297）。

3.1 Guanine-Cytosine堿基對

這個例子的Guanine和Cytosine單體的Gaussian輸入文件，以及計算產生的.fch文件都提供在本文文件包里的GC目錄下了。將1.fch（對應Cytosine）和2.fch（對應Guanine）都拷到Multiwfn目錄下，雙擊ESPpt.bat，就會開始計算，輸出文件被自動拷到了VMD目錄下。然后啟動VMD，在文本窗口輸入pt2命令，得到基于表面頂點繪制的范德華表面穿透圖

由圖可見在形成氫鍵區域，范德華表面穿透很明顯，而且是紅色和藍色相互穿透，體現靜電勢的互補特征，也反映出氫鍵的靜電吸引相互作用的本質。在VMD里按鍵盤上的2，然后點擊恰當的兩個表面頂點，就可以通過測量頂點的距離來考察穿透距離，這在196文中已經說了，請仔細閱讀（由于顯示出來的距離值的文字是白色的，當前背景也是白色的，可能難以觀察。只要把ESPpt2.vmd里的color Display Background white刪掉再作圖，或者VMD文本窗口輸入color Display Background black命令，就可以切換到黑背景。為了便于點擊恰當的表面頂點，建議一次只顯示一個單體的表面頂點，即在VMD main窗口里雙擊vtx1.pdb或vtx2.pdb的單體的"D"標記將之取消顯示，并且恰當旋轉和縮放視角以便于點擊頂點）

下面基于靜電勢著色的電子密度等值面再來繪制穿透圖。雙擊ESPiso.bat，需要運算一會兒，算完后啟動VMD，輸入iso2命令，即得到下圖（對于這種重疊圖不建議用Tachyon渲染，否則等值面交界處會有明顯痕跡，不好看）

3.2 水四聚體

本文文件包里的water_tetramer里面的complex.gjf是含有優化過的四聚體的坐標的Gausian輸入文件，1/2/3/4.gjf是基于其中每個單體坐標構建的Gaussian輸入文件，計算產生的相應的.fch文件也都提供了。將這些fch文件拷到Multiwfn目錄下，然后按照上一節的做法操作，得到如下兩種圖

左邊的圖比較素雅清爽，右邊的圖比較濃郁。當然繪制左圖時也可以把表面頂點尺寸加大，縮窄色彩刻度上下限來讓顏色更艷一些，這些在196文章里都有涉及。

5 關于調節材質和色彩變化

對于有的復雜體系，在默認的EdgyGlass材質下可能繪圖效果不是很理想，比如下面這個團簇體系，畫面看起來很亂，不容易看清楚體系表面靜電勢特征：

這個問題主要是透明度太高了。只要改一下材質即可。做法是進入Graphics - Materials，選擇當前用的材質EdgyGlass，然后調節各個滑條，特別是Opacity（不透明度）。我們把材質改成下面的情況

此時看到的圖像如下所示，看著清楚多了。

如果你不喜歡作圖腳本用的藍-白-紅方式的色彩變化，可以進入Graphics - Colors - Color Scale，用Method下拉框將色彩變化方式設成別的。比如有很多文章用的是用紅色代表負值部分，藍色代表正值部分，為了與這種習俗對應，你就可以把Method設為RWB（Red-White-Blue）。

6 關于顯示色彩刻度軸

如果想把色彩刻度軸顯示出來，一個做法是進入Graphics - Colors，選擇Color Scale標簽頁，然后把色彩刻度條截圖截出來，用photoshop旋轉、拉伸一下，然后把上下限數值標注上去即可。

另一個做法是先把靜電勢圖顯示出來，在VMD里選擇Extensions - Visualization - Color Scale Bar，Autoscale選On，Use Molecule旁邊的項目選density1.cub（假設你是通過iso命令來顯示靜電勢著色的分子表面的），Use Representation選Isosurface。Color of labels選black，Label format選Decimal，然后點Draw Color Scale Bar。然后在VMD Main窗口里，把Color Scale Bar以外的項目的F字母都雙擊點成黑色將它們固定住，然后雙擊Color Scale Bar旁邊的F將之變成紅色使之可以移動。然后點擊VMD的OpenGL圖形窗口將之激活，點t鍵進入視角平移模式，按住鼠標在畫面里拖動，把色彩刻度條挪到合適的位置。之后再按r恢復默認的旋轉模式，在VMD Main窗口里把Color Scale Bar的F雙擊變成黑色使之固定住，而其它體系的F雙擊變成紅色使之解除固定。最后效果如下所示

可以在圖上標注一下或者在圖例里說明一下刻度的單位。對于通過iso命令顯示靜電勢著色等值面的情況，由于靜電勢cube文件里的單位是a.u.，所以刻度條上的單位也是a.u.，也可以自行ps成對應的其它單位的數值。如果你是用pt命令顯示的，如前所述，單位可能是kcal/mol也可能是eV。

如果你想讓iso、iso2命令產生的表面靜電勢圖以eV為單位，按下面的做法（必須用2021-Aug-17及以后更新的Multiwfn）：
把examples\drawESP目錄下的ESPiso_eV.bat和ESPiso_eV.txt拷到當前目錄下，用它們代替前文的ESPiso.bat和ESPiso.txt進行計算。把前述的ESPiso.vmd和ESPiso2.vmd拷到VMD目錄下之后，把里面的set colorlow -0.8和set colorhigh 0.8前面的#去掉使之生效。這樣得到的靜電勢cub文件就是以eV為單位記錄的了，默認色彩刻度下限和上限將是-0.8 eV和0.8 eV，用上述的Color Scale Bar繪制的色彩刻度軸也是以eV為單位的了。PS：如果想以kcal/mol為單位，把ESPiso_eV.txt里的a.u.到eV的轉換系數27.2114改成627.51，并恰當設置ESPiso.vmd和ESPiso2.vmd里的默認的色彩刻度上下限即可。

7 關于色彩刻度范圍

我建議讀者自行根據實際情況在VMD的Graphics - Representation的Trajectory標簽頁的文本框里調節色彩刻度的下限和上限，使得色彩可以盡可能充分地體現分子表面上不同區域靜電勢的差異。比如你感覺紅、藍顏色太淡，則應當把色彩刻度范圍區間設小一點。如果有一大片地方全是相同的紅色或相同的藍色，因而區分不開差異，則色彩刻度范圍需要設寬一些。需要考慮實際效果反復調節到最理想。

對于離子體系，VMD作圖腳本里默認的色彩刻度范圍是非改不可的。比如下圖的陽離子體系，分子表面靜電勢處處為正，而且數值很大，如果用默認色彩刻度的話分子表面顯然就只有一種顏色。對于這種情況，大家應啟動Multiwfn，載入輸入文件，進入主功能12，選0，即做分子表面上的靜電勢的定量分析，然后按照下圖的示意，把色彩刻度的下限和上限分別設成Multiwfn輸出的分子表面靜電勢最小值和最大值（最大值和最小值前頭都有星號提示），此時這個陽離子分子表面不同區域靜電勢的差異就得以體現了。注：對于用ESPext.bat繪圖的情況，從Multiwfn窗口里取a.u.為單位的值。而對于用ESPpt.bat繪圖的情況，單位可能是kcal/mol（中性體系一般如此）也可能是eV（離子體系一般如此），用文本編輯器打開自動被挪到VMD目錄下的vtx1.pdb文件，看一下開頭的注釋行的提示就知道應當讀哪個單位的。

是否對各種體系總是適合使用表面靜電勢最負和最正值分別作為色彩刻度下限和上限呢？答案是否定的。雖然這么做能確保分子表面不同數值靜電勢的位置能對應不同顏色，但這往往不是什么好的選擇，比如：
(1)上、下限是個零碎的小數，令色彩刻度軸顯得不夠工整
(2)上、下限不對稱，導致白色的區域不是正好落在靜電勢為0的地方從而有礙判斷。當然并不需要非得上下限對稱，只不過不對稱的時候需要自己確保色彩變化的中間顏色正好對應靜電勢零點，可以通過Graphics - Colors - Color Scale里的Midpoint進行調節。
(3)視覺效果未必好。比如分子某處靜電勢特別負，若為了照顧這個位置而設置下限，可能導致靜電勢比較負的區域在著色上不容易區分差異。
(4)橫向對比一批分子表面靜電勢時，如果色彩刻度范圍用得不統一，會造成誤判。

8 其它細節

如果你目前還沒根據自己機子實際情況修改Multiwfn的settings.ini文件里控制并行核數的nthreads參數，記得一定要將之改成CPU的物理核心數，要不然沒法發揮CPU的實際計算能力。

本文的腳本的運行機制在《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://sobereva.com/612）中做了十分詳細的介紹。希望大家充分理解腳本中的語句的內容，以便根據自己實際情況修改，從而用起來更方便，比如可以在ESPpt.bat和ESPiso.bat倒數第二行分別加上del D:\study\VMD193\*.pdb和del D:\study\VMD193\*.cub來使得運行時自動把VMD目錄下之前的pdb和cub文件刪除。

本文的.bat批處理文件是Windows下的，如果想在Linux下也類似地通過腳本快速繪制，看Multiwfn手冊4.A.13節的第9部分（此外，也可以手動在性能較好的Linux服務器上把VMD腳本繪圖所需的文件算出來，然后自行拷到Windows的VMD目錄下再用VMD的腳本繪圖）。

ESPpt.txt里的0.15是Multiwfn中定量表面分析時的格點間距，如果改大這個值，表面頂點就會變得更稀疏，同時由于表面頂點數變少了，計算耗時也會降低。ESPiso.txt里面的命令對應于計算電子密度和靜電勢格點數據時分別用high quality grid和low quality grid（后者低于前者是為了節約靜電勢計算時間），對于一般情況這個設置已經效果較理想了，但如果體系涵蓋的空間范圍特別大（比如400個原子以上），必須提升格點質量增加計算的格點數這樣作圖效果才好（當然代價是耗時顯著增加），否則等值面可能很不平滑而且對應靜電勢的顏色很模糊。具體做法是把ESPiso.txt內容改為
5
1
4
0.4
0
5
12
2
2

本文最多考慮到了四聚體，更多聚體的穿透圖也可以類似地繪制，比如繪制8聚體，就把ESPiso.bat和ESPpt.bat里計算命令的序號上限從4擴展到8，并且把ESPiso2.vmd和ESPpt2.vmd里的nsystem參數從4改成8，則運行這些繪圖腳本時VMD目錄下的后綴從1到8號的文件都會被載入和繪制。

常有人問為什么他通過ext命令繪制出的分子表面靜電勢極值點的數目和他自己在Multiwfn里用主功能12做靜電勢的定量分子表面分析給出的不符，這是因為ESPext.bat腳本把格點間距設成了0.15 Bohr，而Multiwfn默認用的是0.25 Bohr。格點間距越小表面頂點確定得越準確。顯然不符的時候，通過ESPext.bat給出的更準。如果想要二者結果一樣，在Multiwfn主功能12里進行計算前先選選項3修改格點間距為0.15 Bohr即可。

對于巨型體系的分子表面靜電勢的繪制，DFT方法優化并產生波函數文件的過程可能比較費勁，但可以使用xtb程序在很便宜的GFN-xTB級別（類似于半經驗的DFT方法）下做優化并產生molden格式的波函數文件，可以直接作為上文提及的腳本的輸入文件。這在《巨大體系的范德華表面靜電勢圖的快速繪制方法》（http://sobereva.com/481）一文中介紹了。

如果你要繪制1000原子以上的特大體系的分子表面靜電勢圖，像本文這樣基于波函數來計算肯定是算不動的，應當使用《基于原子電荷極快速繪制超大體系的分子表面靜電勢圖》（http://www.shanxitv.org/639）里的做法，用Multiwfn基于原子坐標和原子電荷來計算，雖然準確性會打一定折扣，但耗時超級低，對于幾千原子體系在一般個人電腦上也就是耗費幾分鐘的事。

ESPiso.bat里有個-ESPrhoiso 0.001的設置，如果你想了解這是什么含義的話，參看《巨幅降低Multiwfn結合VMD繪制分子表面靜電勢圖耗時的一個關鍵技巧》（http://sobereva.com/602）。

如果體系帶電荷特別多導致分子表面有的區域靜電勢數量級特別大，有可能導致ESPpt.bat涉及的vtx.pdb的B因子段落由于位數限制無法記錄實際的靜電勢數值。解決辦法是使用examples\drawESP目錄下帶_pqr后綴的那些文件代替前文提及的不帶_pqr后綴的那些文件來計算以及在VMD中作圖。帶_pqr后綴的腳本會調用Multiwfn產生記錄分子表面靜電勢的.pqr文件，其中charge那一列被用于記錄靜電勢數值，由于沒有位數限制，靜電勢無論多大都可以被記錄（注意，總是以a.u.為單位記錄）。

使用與本文類似的方法，還可以非常方便地繪制對于研究親電反應位點非常重要的平均局部離子化能(ALIE)著色的分子表面圖，并且ALIE表面極小點也可以一起顯示出來，繪制過程見《使用Multiwfn和VMD繪制平均局部離子化能(ALIE)著色的分子表面圖（含視頻演示）》（http://sobereva.com/514）。

9 總結

本文介紹的繪制分子表面靜電勢圖的方法算是目前最快、最好的繪制方法了，而且用到的Multiwfn和VMD都完全免費。相比之下，用GaussView繪制這類圖的操作步驟多于此文，耗時遠高于此文，效果也不如本文，而且程序是收費的而且挺貴，還沒法把靜電勢極值點顯示出來，也沒法繪制穿透圖，而且只能用fch/fchk作為輸入文件，顯然相對于本文介紹的方法來是不可取的。

如果你以本文方式繪制靜電勢圖發表文章，請務必記得引用Multiwfn程序一啟動時就顯示的程序原文，也請同時按照《Multiwfn使用的高效的靜電勢算法的介紹文章已于PCCP期刊發表！》（http://sobereva.com/614）文末說的引用Multiwfn計算靜電勢的算法在PCCP上發表的原文。如果你同時在圖中顯示了分子表面靜電勢極值點位置，請同時也引用筆者的文章J. Mol. Graph .Model., 38, 314 (2012) DOI: 10.1016/j.jmgm.2012.07.004，此文詳細介紹了Multiwfn中搜索分子表面靜電勢極值點的算法。

10 附：幾個文獻中的應用實例

下面給出幾個通過本文做法繪制的已發表的文章中的圖像作為實例。

筆者專門發表了一篇18碳環電子結構的分析文章，簡介見《一篇最全面、系統的研究新穎獨特的18碳環的理論文章》（http://sobereva.com/524），其中就利用了本文的方法對此體系表面靜電勢分布進行了繪制，并根據《使用Multiwfn結合VMD分析和繪制分子表面靜電勢分布》（http://sobereva.com/196）里的做法對靜電勢不同區間面積分布進行了考察，結果如下圖所示，可見效果極佳！

下面這張圖是筆者研究18碳環弱相互作用的文章Carbon, 171, 514 (2021)中的圖，左邊是18碳環二聚體的，右邊是18碳環在外部吸附水分子。可以見可以通過靜電勢互補原理能極好地解釋為什么二聚體是這種構型。強烈建議大家閱讀此文：《全面探究18碳環獨特的分子間相互作用與pi-pi堆積特征》（http://sobereva.com/572）。更多的筆者關于18碳環的研究工作看http://sobereva.com/carbon_ring.html。

下圖是《8字形雙環分子對18碳環的獨特吸附行為的量子化學、波函數分析與分子動力學研究》（http://www.shanxitv.org/674）介紹的筆者的Phys. Chem. Chem. Phys., 25, 16707 (2023)文中的圖，充分展現了18碳環吸附進OPP的大環部分后與其范德華表面的相互穿透，體現出這一對主、客體分子之間形狀極其般配。

下面這張圖右邊是J. Mol. Struct.期刊上的一篇研究共晶的文章（DOI: 10.1016/j.molstruc.2020.128480）中對晶體局部的四個分子按照本文的做法繪制的靜電勢穿透圖，可見將相互作用展現得非常清楚直觀

下面是Ind. Eng. Chem. Res. (2020)（DOI: 10.1021/acs.iecr.0c01451）文章中使用本文方法繪制的圖，體現的是1-乙基-3-甲基咪唑鎓氯化物和咪唑選擇性吸附SO2，給出了六種構型。可見圖中非常清晰地展現出氯離子與SO2靜電勢為正的硫區域接觸，其它分子的靜電勢為正的部分與SO2靜電勢為負的氧區域接觸，很好地體現出通過靜電作用使得SO2被牢固吸附。

附：圖像無法繪制出來的排查方法

老有人問怎么按照此文的方法繪制不出來、VMD文本窗口提示找不到mol1.pdb、運行iso后提示Unable to load file 'density1.cub'...之類，幾乎是周經問題，在這里明確統一回復一遍

1 首先必須確保用的是Multiwfn官網上的最新版本。VMD最好用1.9.3版（至少這個版本肯定沒問題。1.9.4的非正式版絕對不要用！）
2 確保ESPiso.bat或ESPpt.bat里的VMD路徑設對了，否則腳本調用Multiwfn運行完產生的pdb、cub文件沒法被自動挪到VMD目錄下，顯然運行iso等命令時VMD沒法加載文件。如果VMD的安裝路徑里有空格，兩邊必須用雙引號括起來。注：不要把VMD安裝到默認的路徑下，路徑又長又有空格，初學者老是容易寫不對！而且可能操作系統或某些防護軟件的安全設置會妨礙在其目錄下正常創建文件。筆者強烈建議把VMD安裝到簡短的比如D:\study\VMD193這樣的目錄中
3 確保上述bat文件里給Multiwfn用的輸入文件的路徑無誤。如果fch文件是在當前目錄下的，直接寫文件名就可以。如果你實際的文件是fchk后綴的，顯然要么把實際文件的后綴手動改成fch，要么把bat文件里的輸入文件的后綴改成fchk
4 確保用的輸入文件的文件格式合理，諸如用.xyz、.pdb之類根本沒有波函數信息的格式顯然不行。仔細看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）了解什么格式包含波函數信息，用諸如.molden、.gms、mwfn、.wfn、.wfx等文件格式也都是可以的
5 確保文中的.vmd后綴的腳本文件都放到了VMD默認的目錄下（具體來說就是啟動VMD后，在文本窗口里輸入pwd并按回車，顯示出的目錄就是.vmd腳本文件該放的目錄。正常情況就是VMD可執行文件所在的目錄）
6 確保改過vmd.rc后保存了此文件，并且確保VMD是在保存此文件后才啟動的，否則不生效
7 如果你要讓Multiwfn調用cubegen做靜電勢計算且當前體系特別大，應確保已經設了GAUSS_MEMDEF環境變量，否則超過默認分配的內存上限時cubegen會崩潰。關于這個環境變量的含義和設置，參考《巨大體系的范德華表面靜電勢圖的快速繪制方法》（http://www.shanxitv.org/481）
8 死活搞不明白原因的話，仔細從頭閱讀《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）直到介紹完ESPiso.bat的地方為止，由此搞清楚bat批處理文件的運作原理（極其簡單！）。如果還沒醒悟過來問題所在，在當前目錄下下打開cmd窗口（不知道做怎么看http://bbs.keinsci.com/thread-22940-1-1.html），輸入bat腳本名字運行之，從提示上判斷哪里有問題（如果提示找不到文件，明擺著就是路徑寫錯了，沒別的原因！如果還不懂路徑怎么寫錯了，就找個身邊有最基本計算機常識的人給你看）。如果自行不會判斷、死活無法解決，把cmd窗口輸出的信息連同當前目錄的窗口截圖貼到Multiwfn的官方論壇上，筆者會回復