使用Multiwfn作電子密度差圖
使用Multiwfn作電子密度差圖
文/Sobereva @北京科音
First release: 2011-Dec-24 Last update: 2019-May-22
1 前言
在Multiwfn（http://www.shanxitv.org/multiwfn）里可以十分輕松地作密度差圖，是繪制密度差圖最方便的工具，沒有之一。在程序手冊里也已經給出幾個例子很清楚地介紹了繪制密度差圖的方法（4.4.7節的變形密度圖，4.5.4、4.13.6、4.18.3節的體系在不同狀態間的密度差圖、以及4.5.5、4.4.8節的片段密度差圖）。在此寫個專門帖子通過實例來說說繪制過程，會把每一步都做出詳細解釋。本文同時也會談到密度差圖的意義、用處。另外，對于作密度差圖常常被忽視的原子密度球對稱化問題也會著重說一下。
本文所說的密度都是指電子密度。利用和本文完全相同的方法，還可以作其它函數的差值圖，比如靜電勢差值圖、動能密度差值圖、電子定域化函數(ELF)差值圖等等，只要在選擇實空間函數時選擇對應的實空間函數就可以了。
如果讀者對Multiwfn不了解，建議參看這些文章：《Multiwfn入門tips》（http://www.shanxitv.org/167）、《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文只是簡單提及了特定情況下產生Multiwfn輸入文件的做法，更全面的介紹見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。利用類似作密度差的方法還可以得到（超）極化率密度，這對于討論體系不同區域對（超）極化率的貢獻、探究不同計算級別導致（超）極化率計算結果存在差異的原因特別有用，此問題另有專門的帖子，見《使用Multiwfn計算超極化率密度》（http://www.shanxitv.org/305）。
2 密度差圖的種類
“密度差”就是指一些體系在各自特定的狀態下彼此間的密度的差值，它是三維實空間函數。有三種情況常會涉及：
1 某個體系的某個狀態的密度減去這個體系在另外一個狀態的密度。例如，某分子的激發態的密度減去它在基態時的密度、某分子在有外加電場時的密度減去它在自由狀態下的密度。注意在計算這兩個狀態的密度時分子的笛卡爾坐標必須一致，不能在各自狀態下分別優化，而只能共用其中一個狀態的幾何構型。同時也不能允許量化程序為了能利用分子對稱性而自動旋轉、平移坐標。否則都導致兩個狀態下原子坐標沒法對應，密度差也就沒什么意義。
2 某個體系的密度減去組成它的各個片段的密度。例如，水分子二聚體的密度減去這兩個水分子都在孤立時的密度。同樣要注意，求各個片段的密度時，它們的笛卡爾坐標必須與在整個體系中的坐標一致。
3 某個體系的密度減去組成它的各個原子在自由狀態下的密度，這也被稱為變形密度(Deformation density)。這種密度差對于研究分子形成過程中電子密度的變化、探究化學鍵的本質很有用。
密度差必須作成圖才便于被直觀地考察，圖分為三類：
1 曲線圖：展現某條直線上密度差的數值，這種圖對于密度差的研究用得較少。
2 平面圖：能夠清楚地展現一個截面上各處的密度差數值，十分有用，其中最適合展現密度差的是等值線圖。平面圖的缺點是難以從一個截面上了解整個空間中密度差的數值分布。
3 等值面圖：在等值面上每個點的數值都等于人為指定的數值(這個值被稱為isovalue)，通過它可以很直觀地了解密度差在三維空間中的分布。但是等值面里面和外面的密度差數值沒法一下子看清楚，必須反復調節isovalue考察等值面的變化才能獲知。等值面圖和平面圖的用處是互補的。
3 Multiwfn中的自定義運算功能
Multiwfn可以直接計算密度并做曲線、平面、等值面圖，Multiwfn也提供了自定義運算(custom operation)功能，用于對多個波函數的密度、靜電勢等實空間函數進行相互運算，求密度差實際上只是Multiwfn的自定義運算功能的一個具體應用而已。
比如，啟動Multiwfn后載入的波函數是A.wfn，在自定義運算中設定有3個波函數將要依次對它進行運算，且算符和文件名分別寫
+,madoka.wfn
x,TMA.wfn
-,homura.wfn
然后，假設要計算的屬性選擇的是電子密度，那么Multiwfn會依次計算這四個波函數文件的電子密度格點數據，并且對這些格點數據做如下運算：(A.wfn的密度+madoka.wfn的密度) * TMA.wfn的密度 - homora.wfn的密度
在計算后三個波函數的格點數據時Multiwfn都會自動使用與計算A.wfn格點數據時相同的格點設定，無需用戶考慮格點位置的對應性問題。
Multiwfn可以對無限多個波函數進行自定義運算，支持的算符是+ - x /，對應加減乘除（乘號寫為字母x）。
下面將通過一些實例，說明怎么用自定義運算功能繪制各種類型密度差圖。
4 某體系某狀態減去此體系另一個狀態的密度差圖
4.1 施加電場后COCl2電子密度的變化的等值面圖
在外場作用下，電子密度分布會被極化。這里考察一下COCl2分子在C-O軸線方向上加上0.03a.u.均勻電場（1a.u.=51.421 V/Angstrom）后電子密度的變化。注意電場作用也會誘導體系幾何結構發生改變，但繪制密度差要求幾何結構保持不變，所以這里不考慮電場對分子結構的影響，也因此不在電場下重新優化結構。
首先生成COCl2在自由狀態下的波函數，將以下內容保存為Gaussian輸入文件并執行即可，wfn文件輸出路徑和文件名可以自定。這個結構是事先在b3lyp/6-31+G*下優化好的。
# b3lyp/6-31+G* out=wfn nosymm
[空行]
COCl2 opted
[空行]
0 1
C 0.00000000 0.00000000 0.49676300
Cl 0.00000000 1.46464400 -0.48297200
Cl 0.00000000 -1.46464400 -0.48297200
O 0.00000000 0.00000000 1.68005900
[空行]
C:\COCl2.wfn
算完后將這個輸入文件的route section部分末尾加上field=z+300關鍵詞，并將輸出的wfn文件名改為C:\COCl2_z+300.wfn，然后重新算一次以得到在Z方向上（C-O軸線方向）上加了0.03a.u.均勻電場后的波函數文件。
Gaussian會自動將體系坐標進行旋轉、平移使之處于標準朝向下，以便利用對稱性節約計算耗時。如果你對Gaussian缺乏經驗的話，在計算用于做密度差的波函數文件時建議寫上nosymm，并且使用笛卡爾坐標，這樣可避免Gaussian自動做這種調整，從而確保了分子在三維空間中的位置、朝向是固定的，和輸入文件嚴格一致。對于外加電場的計算，如果不寫nosymm，而且體系在輸入文件里的坐標不是處在標準朝向，那么計算時經過Gaussian自動調整體系朝向后，電場和體系的相對方向就亂套了。
我們先來作等值面圖。啟動Multiwfn，依次輸入
C:\COCl2_z+300.wfn
5 //生成格點文件
0 //設定自定義運算
1 //只有1個文件要對已讀入的波函數進行運算
-,C:\COCl2.wfn //運算符（減號）和文件名
1 //選電子密度函數。因此接下來得到的就是COCl2_z+300.wfn的密度減去COCl2.wfn的密度
2 //中等質量格點精度
現在Multiwfn開始依次計算這兩個波函數文件的密度的格點數據，之后自動相減。接下來有諸多選項可以后處理生成的格點數據，用選項2可以將格點數據輸出為cube文件，之后可以載入進VMD、GaussView、ChimeraX等軟件觀看等值面。Multiwfn也自帶了觀看等值面的功能，這里選擇選項-1來進入，如下圖所示。由于默認的isovalue為0.05，而密度差數據通常很小，所以要將isovalue值降低一些，比如降到0.005，才能比較清楚地看到密度差圖。isovalue的大小的設定沒絕對標準，只要設定一個值，能讓密度差特征充分清楚地展現出來就行了。設得太大時等值面會很小而看不清楚甚至消失不見，而設得太小會使等值面涵蓋的空間范圍太大，混在一起，看不出到底哪些區域的電子密度變化較大。
密度差數值與isovalue的符號相同和相反的等值面分別以綠色和藍色顯示。圖中isovalue被設為0.005，也就是說綠色等值面上密度差值都為+0.005，體現了加了電場后電子密度增加的區域；藍色的面上密度差都為-0.005，體現了加了電場后電子密度降低的區域。本例中外加電場矢量是由Z正方向指向負方向的，受其影響電子密度應趨向于向Z軸正方向轉移，圖中也驗證了這一點，可以看到主要發生極化的是C-O鍵的pi電子云。下個例子將繪制穿過C-O鍵的垂直于分子平面的截面上的密度差圖，這個截面對應于Y=0的XZ平面。
如果你想要更好的繪制效果，建議將導出的cube文件弄到VMD里顯示成等值面。在VMD里將cube文件里的格點數據顯示成等值面的操作可參考《使用Multiwfn結合VMD繪制自旋密度等值面圖》https://www.bilibili.com/video/av26312131）。
4.2 施加電場后COCl2電子密度的變化的等值線圖
如果你剛做完上例子的等值面圖，可以在格點數據后處理界面中選擇0從而直接退回到主菜單。如果直接做這個例子，那么還是先啟動Multiwfn，然后輸入C:\COCl2_z+300.wfn來載入之。之后依次輸入
4 //繪制平面圖
0 //設定自定義運算，后面幾步和4.1完全相同
1
-,C:\COCl2.wfn
1
2 //等值線圖
直接按回車使用推薦的格點設定，即平面圖中兩個方向都是200個點
2 //XZ平面
0 //表明這是Y=0的XZ平面
瞬間就彈出了下面的圖
實線代表加電場后電子密度增加的區域，虛線代表減小的區域。這幅圖比等值面圖展現的密度差更為清楚具體。在圖上點擊右鍵關閉它之后可以看到一個菜單，其中有豐富的選項，用于保存圖片、調節等值線設定、調節坐標軸刻度、選擇是否將等值線的數值標在圖中等等，請大家自行嘗試，在手冊對應的章節也都有細致的說明。從這個例子可見使用Multiwfn做密度差平面圖是極其容易的，而且圖像精細漂亮，可直接作為文獻的插圖。既不需要費時費事地去自行計算格點數據，也完全不需要在Sigmaplot、Origin里面折騰數據和調節設定。
在手冊的4.13.6節有個繪制聚炔烴在外加電場后電子密度變化的例子，建議看看，其中還額外繪制了電荷轉移曲線，對于定量分析密度差變化很有幫助。
4.3 電子激發前后的密度差圖
還是以COCl2為例子，這次來看看TDDFT計算得到的第一單重激發態(S1)相對于基態的電子密度的變化。如果你對激發態計算一無所知，建議參看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）。這里還是用4.1節的COCl2的輸入文件，但這次將route section改為
# B3LYP/6-31+G* out=wfn TD(root=1) density
并且將末尾定義的wfn文件名改為C:\COCl2_S1.wfn。用Gaussian計算之，COCl2_S1.wfn就將包含S1態的波函數信息。此例的root=1可以省，因為默認就是root=1。density關鍵詞一定要加，否則wfn文件包含的是基態波函數，而不是你指定的第root態的（但從后期的G09開始，density可以不用寫，因為out=wfn默認用了density）。
注：這里我們沒寫nosymm，是因為它對于當前情況沒什么用處，因為基態和激發態結構是相同的，Gaussian雖然會自動旋轉、平移體系，但這么自動調整之后基態和激發態的坐標依然是一致的。而且我們也沒加外電場，體系的朝向無所謂。如果寫了nosymm，雖然對結果也沒害處，但Gaussian沒法利用對稱性來節約計算耗時。
做等值面圖和做截面圖的方法和前兩節的過程完全一致，將C:\COCl2_z+300.wfn替換為C:\COCl2_S1.wfn就行了。得到的等值面圖如下所示，isovalue=0.03。可以看到在體系從S0垂直激發到S1后，碳和氧的px軌道上電子密度都有所增加，而且增加的區域并沒有相連，這是由于電子主要躍遷到了pi反鍵軌道所致。相應地，氧的py軌道（此例對應于氧的孤對電子）上電子密度有所減小。因此這種激發模式可指認為LP->pi*躍遷。
在手冊4.18.3節有個繪制D-A體系的電荷轉移激發態的密度差的例子，隨后還對密度差進一步轉化，得到了轉移距離、正負電荷區域重疊程度等定量指標，對于分析電子激發模式很有用，建議看看。
筆者還寫過一篇文章《使用Multiwfn計算激發態之間的密度差》（http://www.shanxitv.org/429）。如果你打算對好多個激發態之間求密度差，那么一個一個生成激發態的波函數文件又耗時又麻煩，而如此文所示的，使用Multiwfn可以一次性產生所有激發態的波函數文件，繪制密度差就省事多了。
4.4 計算Fukui函數預測反應位點
Fukui（福井）函數是常用的預測反應位點的實空間函數，尤其適合軟酸軟堿之間的反應（靜電勢的分析方法更適合用于分析硬酸硬堿反應）。Fukui函數有三類，分別用于預測親核、親電、自由基反應位點。預測親電反應的函數f- = ρ(N) - ρ(N-1)，其中ρ(N)是分子正常狀態下的密度函數，ρ(N-1)是分子少了一個電子，或者說分子的凈電荷為+1狀態下的密度。f-在哪個原子附近數值越大，說明親電反應越容易發生在哪個原子上。注意計算ρ(N-1)時幾何結構必須使用計算ρ(N)時的，這不僅是繪制密度差圖的要求，這也是Fukui函數的定義本身所要求的。
這里計算呋喃的f-函數的等值面圖。首先計算正常狀態下的波函數文件，Gaussian輸入文件如下，結構已在b3lyp/6-31G*下優化好。nosymm不是必需的，理由和上一節一樣。
# b3lyp/6-31g* out=wfn
[空行]
opted b3lyp/6-31G*
[空行]
0 1
C -0.66222300 -0.93877200 0.00000000
C 0.69784300 -0.97448500 0.00000000
C 1.13199200 0.39408700 0.00000000
C 0.00000000 1.14881300 0.00000000
O -1.10670000 0.35104000 0.00000000
H -1.43048300 -1.69690900 0.00000000
H 1.31938500 -1.85911500 0.00000000
H 2.14981500 0.75868100 0.00000000
H -0.19078800 2.21116400 0.00000000
[空行]
C:\furan.wfn
算完后把凈電荷和自旋多重度從0 1改成1 2，C:\furan.wfn改成C:\furan+1.wfn，重算一遍得到N-1電子時的波函數文件。
之后過程和4.1節完全相同，即在Multiwfn里依次輸入
C:\furan.wfn
5
0
1
-,C:\furan+1.wfn
1
2
算完后選-1預覽等值面，適當調節isovalue，以使得等值面區域只集中在個別原子上。當前體系isovalue設為0.015時比較清楚，如下圖所示，等值面幾乎只在氧的鄰位上出現，而且綠色區域顯著大于藍色區域，這顯示出氧鄰位碳附近的f-函數值非常正，而且明顯大于間位碳的，預測了呋喃的親電反應會發生在氧的鄰位碳上，這與實驗結論是一致的。
值得一提的是，Multiwfn有一個專門計算概念密度泛函框架里定義的包括福井函數在內的各種量的功能。用這個功能，你甚至都不用手動準備波函數文件了，只要提供一個含有優化后的結構的文件，在Multiwfn里敲擊下鍵盤就能直接給出各類福井函數，方便至極，請閱讀：《使用Multiwfn超級方便地計算出概念密度泛函理論中定義的各種量》（http://www.shanxitv.org/484）。
如果對于反應活性的預測方法感興趣，想了解更多理論知識和相應的計算流程，推薦閱讀幻燈片《反應位點的預測》（http://www.shanxitv.org/234）、Multiwfn手冊4.A.4節、《親電取代反應中活性位點預測方法的比較》（物理化學學報，30，628 (2014)，http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）、《概念密度泛函綜述和重要文獻合集》（http://bbs.keinsci.com/thread-384-1-1.html）。筆者還寫過一個相關文章《使用Multiwfn+VMD以原子著色方式表現原子電荷、自旋布居、電荷轉移、簡縮福井函數》（http://www.shanxitv.org/425）。
5 繪制某個體系與組成它的各個片段的密度差圖實例
這一節以水二聚體為例介紹怎么繪制某個體系與組成它的各個片段的密度差圖。二聚體的幾何結構直接從S22弱相互作用數據庫獲得，如下圖所示
首先獲得二聚體的波函數文件，輸入文件如下（關于計算電子密度的理論級別的討論見文末注1）。這個例子必須寫nosymm，否則由于Gaussian會自動調整坐標到標準朝向，將導致單體的坐標和它們在復合物中的坐標不對應，使密度差失去意義。
# B3LYP/6-311G** out=wfn nosymm
[空行]
From S22
[空行]
0 1
O -1.551007 -0.114520 0.000000
H -1.934259 0.762503 0.000000
H -0.599677 0.040712 0.000000
O 1.350625 0.111469 0.000000
H 1.680398 -0.373741 -0.758561
H 1.680398 -0.373741 0.758561
[空行]
C:\waterdimer.wfn
二聚體算完之后，將后三個原子從輸入文件中刪掉，并將輸出波函數文件名改為C:\water1.wfn，然后計算。算完后再將后三個原子恢復回來，而將前三個原子從輸入文件中刪掉，將輸出波函數文件名改為C:\water2.wfn，再次計算。此時water1.wfn和water2.wfn里面就分別是第一個水分子和第二個水分子的波函數信息。
啟動Multiwfn，依次輸入
C:\waterdimer.wfn
4 //繪制平面圖
0 //設定自定義運算
2 //有兩個波函數將從當前載入的波函數中減去
-,C:\water1.wfn
-,C:\water2.wfn
1 //電子密度函數。最終得到的密度就是waterdimer.wfn的密度減去water1.wfn的密度再減去water2.wfn的密度
2 //等值線圖
直接按回車使用推薦的格點設定
1 //用XY平面作為繪圖平面
0 //繪制Z=0的XY平面
馬上密度差圖就會蹦出來，如下圖所示。之后還可以關閉圖像，用后處理界面的選項改進圖像效果、修改等值線設定。
從這幅圖上能看出兩個水分子形成復合物時電子密度是如何轉移的，實線和虛線代表結合后密度增加和減少的區域。從圖中看到形成氫鍵的H和O的電子密度都顯著降低，這體現出形成二聚體時由于Pauli互斥作用令電子密度分布彼此排斥。同時，可以看到O-H之間密度變化有個截面，這也反映了O-H鍵的電子共享程度有一定弱化，使得O-H鍵被削弱（反映為鍵長增長，振動頻率降低）。
此例要繪制成等值面圖也非常容易。啟動Multiwfn，依次輸入
C:\waterdimer.wfn
5 //計算格點數據
0 //設定自定義運算
2 //有兩個波函數將從當前載入的波函數中減去
-,C:\water1.wfn
-,C:\water2.wfn
1
2 //中等質量格點
-1 //觀看等值面圖
將isovalue調至0.001，看到下圖。綠色和藍色分別對應密度差=0.001和-0.001的等值面，勾勒出電子密度增加和減少的主要區域。
手冊4.5.5節也有個繪制片段密度差的例子建議看看。
在Multiwfn中繪制體系與構成它的多個片段間的密度或者其它函數的差值圖也非常容易。在Multiwfn手冊4.4.8節的例子中展示了對水四聚體繪制密度差的等值線圖、等值面圖，以及填色圖形式的ELF差值圖的過程，所得圖像分別如下所示：
6 繪制體系與組成它的每個原子的密度差圖
做完上個例子后，一定已經知道怎么做體系與組成它的各原子的密度差圖了（下面稱為變形密度），只需要在填入每個片段波函數文件名的地方改為填入每個原子的波函數文件名就行了。然而，原子的數目通常很多，依次計算原子的波函數文件并加入自定義運算列表里是很費事的，而且還得手動去做原子密度的球對稱化（第8節將談到）。為了使變形密度圖能夠方便地繪制出來，Multiwfn提供了特殊的處理方法。有兩種情況，一種是讓Multiwfn直接調用Gaussian產生原子波函數文件，如下面6.1節所述；另一種方法是直接用Multiwfn自帶的原子波函數文件，如下面6.2節所述，這樣既不用自己算原子波函數文件也不用調用Gaussian了，最為方便。
6.1 自動調用Gaussian生成原子波函數文件
Multiwfn可以自動調用Gaussian生成原子波函數文件并作球對稱化，并自動將其中心挪到各個對應的原子位置上，求密度并作差值，完全不用人工干預。這里來繪制鎂卟吩平面的變形密度等值線圖。此體系的波函數文件MN.wfn是在B3LYP/6-31G*下生成的，結構如下
首先確保Gaussian環境變量已經設好，否則Multiwfn在調用Gaussian時會出錯。對于Win7，設定方法是進入“系統”-“高級系統設置”-“高級”標簽頁，點“環境變量”，在用戶變量中點擊“新建”，變量名填GAUSS_EXEDIR，變量值填Gaussian的安裝路徑，比如D:\study\g09w。設過一次之后以后就不用再重新設了。
啟動Multiwfn依次輸入
MN.wfn
4 //繪制平面圖
-2 //獲取變形屬性（注2）。這里的“屬性”是指接下來你所選的實空間函數，因為接下來要選電子密度函數，所以將得到變形密度
b3lyp/6-31G* //設定計算原子波函數所用的方法和基組
D:\study\g09w\g09.exe //Gaussian可執行文件的路徑。注意可執行文件名是g09.exe而不是g09w.exe。如果你懶得每次做變形密度圖時都重新輸入一遍路徑，那么就在settings.ini文件里將gaupath=后面空一格寫上這個路徑，并且這路徑兩邊要用雙引號擴住，比如gaupath= "D:\study\G09W\g09.exe"。當Multiwfn發現這個路徑正確，就不會再讓你手動輸入路徑了。
現在Multiwfn自動調用Gaussian計算每個原子的波函數文件，由于此體系只含C H N Mg這四種元素，所以會調用Gaussian四次。待算完后接著輸入
1 //電子密度函數
2 //等值線圖
敲回車，用默認格點設定
4 //由于MN.wfn中鎂卟吩不是恰好處在某個XY或YZ或XZ平面上，所以這次還是通過三個原子來定義平面。
19,22,23 //用19、22、23號原子定義繪圖平面。這三個原子位于鎂卟吩的三個頂角，因此做出來的圖會覆蓋整個鎂卟吩。（也可以像第5節那樣通過設定延展范圍來調節作圖區域大小）
很快就得到下圖
可以看到，在每個共價鍵之間都有明顯的正值區域，表明共價鍵的形成伴隨著電子密度在成鍵原子間聚集。N和Mg之間是離子鍵，如預期的，并沒有電子密度以兩原子之間為中心聚集。由N向著Mg方向凸起來的等值線區域表現的是N與C成鍵后，在那片區域凝聚起來的孤對電子。
6.2 使用自帶的原子波函數文件
Multiwfn自帶了前四個周期元素的原子波函數文件，在examples\atomwfn里面，如果你的體系只包含前四周期的元素，那么直接利用這些就可以了，最為省事方便。更詳細的說明見下一節。
為了利用這些自帶的原子波函數文件，要把atomwfn目錄拷貝到當前目錄下。然后，再去做6.1節的例子。當你選-2后，Multiwfn發現所有要用到的原子波函數文件在當前目錄下的atomwfn目錄中都找到了，于是就不再提示你用Gaussian來計算了。接下來的操作和前面一樣，依次選擇電子密度、設定作圖類型等等。
（注：所謂的“當前目錄”是調用Multiwfn時所在目錄。如果你是通過雙擊圖標來啟動Multiwfn，則當前目錄就是指Multiwfn可執行文件所在的目錄。如果是通過命令行方式調用，比如你是在/sob/love/Mio_Akiyama下調用的Multiwfn，那么/sob/love/Mio_Akiyama這個目錄就叫當前目錄）
如果你的體系包含第四周期以后的元素，比如Au，那么你就必須自己計算Au原子得到它的.wfn文件（Au.wfn），所用基組或贗勢和你計算分子時給Au用的應盡量相同。然后將Au.wfn放到當前目錄下的atomwfn目錄下，來讓Multiwfn自動使用之。（如果原子的密度分布不是球對稱的，建議自行對原子波函數文件做球對稱化，見第8節）
7 Multiwfn處理原子波函數的內部機制
本節介紹Multiwfn處理原子波函數的內部機制，以便用戶更靈活地、更方便地運用Multiwfn做變形密度圖。如果看了前面的介紹已經能夠完全滿足你的實際需求，得到了想要的結果，且懶得了解更多細節，此節和下一節都可以跳過去。
在用戶啟動Multiwfn后選擇了作圖的類型，并選了-2以獲取變形屬性時，Multiwfn會執行以下步驟：
1 刪除當前目錄下wfntmp文件夾然后重建此文件夾。此文件夾用于儲存計算原子波函數的臨時文件。
2 將當前目錄下atomwfn目錄下所有波函數文件拷貝到wfntmp文件夾內。
3 檢查wfntmp文件夾是否已經有了當前體系中所有元素的波函數文件。例如，如果發現此文件夾里有N .wfn文件（元素符號為單字符時文件名必須留一個空格），它就被認為是氮元素的波函數文件。元素的波函數文件里原子中心必須位于原點。
4 如果體系中所有元素的波函數文件都能在wfntmp文件夾內找到，就跳到第7步。如果少了任何一個元素的波函數文件，就讓用戶輸入計算它們的理論方法和基組。然后檢查settings.ini下gaupath的路徑是否正確，如果不正確就讓用戶輸入Gaussian的可執行文件路徑。
5 在wfntmp目錄下生成缺失的元素的Gaussian輸入文件，調用Gaussian執行之，以得到它們的波函數文件。Gaussian任務的輸出文件和波函數文件也都生成在wfntmp目錄下。
6 對上一步生成的元素波函數文件做球對稱化處理，這將會改寫元素的波函數文件。
7 將每個元素的波函數文件根據體系中相應元素原子的位置復制成原子的波函數文件。比如體系中3、6、7號原子都是氮原子，則程序會將N .wfn里的原子坐標分別移動到3、6、7號原子位置上并復制一份，分別命名為N 3.wfn、N 6.wfn和N 7.wfn。
8 將所有原子波函數加入自定義運算列表里，運算符都是負號。表明它們的屬性將依次從體系中減去。
如果覺得每次調用Gaussian重新計算各元素的波函數太麻煩，可以自己先手動算好，并做球對稱化處理，改成 [元素名].wfn 格式后放到atomwfn目錄下。由上述過程描述可知，Gaussian發現已經有了相應元素的波函數文件就不會再通過Gaussian計算了。如前所述，在Multiwfn的examples文件夾里就有一個atomwfn文件夾，里面存的是前四周期在ROHF（主族）或UB3LYP（過渡金屬）下結合6-31G*事先計算好的元素波函數文件（已球對稱化）。將這個atomwfn文件夾移動到當前目錄下，就可以省得再調用Gaussian了。之所以沒有自帶>=第五周期的元素，是因為這些重原子計算通常需要用贗勢，但是贗勢種類很多，有大核有小核，也可能有人用全電子基組，所以沒法統一提供原子波函數文件，而只有和實際體系所用的對應上，密度差的結果才有意義。
如果想直接得到特定方法和基組下已球對稱化好的元素波函數文件，目的是將之放到atomwfn目錄下以避免每次都調用Gaussian計算它們，那么可以隨便找個包含那些元素的波函數文件載入進Multiwfn，選完繪圖類型后選-2，使元素波函數文件在wfntmp下生成，然后將它們拷到atomwfn目錄下即可。另外，如果想一次就把前四個周期的元素波函數文件都生成出來，那么可以載入examples目錄下的genatmwfn.pdb文件，它包含了前4個周期所有元素，因此用-2選項時就會生成所有這些元素的波函數文件。
Multiwfn只能對前四個周期元素自動通過Gaussian計算波函數文件并球對稱化。在提示輸入理論方法和波函數時，如果只輸入基組，比如6-31G*，那么就會默認用ROHF計算主族元素，用UB3LYP計算過渡金屬元素，這也是推薦的方法。或許這樣的理論方法和你的體系波函數文件所用方法不符，但實際上沒什么問題，因為原子密度對理論方法依賴性并不大，對基組依賴性也不很大，所以examples\atomwfn目錄下預置的那些元素波函數文件就足以用于做大部分體系的密度差圖了，起碼定性是正確的。
如果同時輸入方法名和基組名，那么必須嚴格按照諸如B3LYP/6-31G**這樣的格式，不要在前面加上RO、R或者U，因為Multiwfn對于過渡金屬會自動加上U，對主族自動加上RO。方法名只能用HF或DFT（雙雜化泛函除外），Multiwfn的球對稱方法無法處理其它情況。如果在自動調用Gaussian計算波函數文件時出錯（現象是突然退出），那么應檢查wfntmp目錄下自動生成的Gaussian輸入文件的route section語法是否正確，以及Gaussian產生的out文件中是否出現錯誤，比如沒能收斂（建議換個泛函）。
如果你的體系是混合基組，尤其是常見的使用贗勢的情況，比如計算二茂鐵（鐵用lanl2TZ、其它原子用6-31G*），那么建議你自行在lanl2TZ下計算出鐵的波函數文件（最好也自行對之作球對稱化處理），改名為Fe.wfn后放到atomwfn文件下，然后按照正常步驟做變形密度圖，填入基組時填6-31G*。這樣Multiwfn就會調用Gaussian在6-31G*下生成Fe以外原子的波函數文件，而對Fe直接使用你提供的lanl2TZ下的波函數文件。
8 原子密度的球對稱化問題
原子在自由狀態下的密度是球對稱的。變形密度是實際體系與組成它的各個原子在自由狀態下密度的差值，這也就要求原子波函數文件對應的密度必須是球對稱的。然而，對一些元素的一些狀態，在一般量子化學計算中得到的密度卻不是球對稱的。比如基態為s2p2組態的碳原子，有兩個p軌道上占了單電子，還有一個p軌道是空的，所以空著的p軌道方向上電子密度低，碳原子的密度成了扁橢圓的（而在現實中，我們無法觀測到空p軌道朝向哪里，只能觀測到平均狀態--球對稱的密度）。原子密度不滿足球對稱會造成變形密度出現不合理性。比方說CH4的變形密度，本應該是有Td對稱性的，如果用了扁橢圓的碳的密度，就會導致變形密度的對稱性被破壞，缺失了物理意義。再比如，基態氧分子是s2p4組態，有一個雙占據p軌道和兩個單占據p軌道，因此氧的電子密度也不是球形的，而是長橢圓形。假設要算O2分子的變形密度，那么計算氧原子的密度時應當讓雙占據的p朝向什么方向？讓雙占據p軌道順著分子軸沒什么道理，畢竟兩個雙占據軌道在一起不會成鍵；然而如果將它沖著垂直于分子軸的方向，那么變形密度就不是本應有的D∞h對稱性了。
因此，生成變形密度時必須要求原子波函數對應的密度是球對稱的。對于第一、二、五主族元素，惰性氣體元素，半滿或滿占據的過渡金屬元素，一般量化算得的基態波函數的密度本身就是球對稱的，所以不必去做球對稱化處理。而其它元素就必須考慮如何獲得球對稱的密度。一般有這么幾種辦法：
1 用密度是球對稱的組態代替基態。比如C，可以用sp3代替sp2p2，由于C在分子中往往是sp3雜化，所以這么做是有道理的。然而對于基態組態為s2p1的B，就沒有合適的生成球對稱密度的組態了，如果把p軌道的那個電子激發到3s軌道上，顯然密度會太過于彌散，與真實很不符。所以這種球對稱方法并不普適。
2 用分數占據的DFT計算。在這種形式中，不要求每個軌道電子占據數都是整數，比如可以讓B的三個p軌道都占1/3個電子，這保證了密度是球對稱的。這種方法很好，但是分數占據的DFT計算并沒有被廣泛支持，Gaussian也不支持（但可以用NWChem，見《使用NWChem做分數占據數的DFT計算》http://www.shanxitv.org/363）
3 做MCSCF計算。最后得到的其中每個組態波函數的密度不是球對稱的，但MCSCF波函數的密度會是球對稱的，產生的自然軌道及占據數和分數占據的DFT的情況是相似的。然而做這種MCSCF會把問題復雜化。
4 用一些s型GTF或STO函數擬合非球對稱的原子密度，之后由這些s型函數重現出的密度就是球對稱的了。這種方法也相對麻煩，費時費事。
在Multiwfn里，對原子密度球對稱化的方式簡單高效，對不同類型元素用不同處理方法。對于第三主族，將限制性開殼層計算中那個單占據p軌道旋轉向其它兩個方向并且都復制一份，然后將這三個軌道的占據數都設成1/3；對于第四主族，提供了兩種球對稱化方法，默認的是使用sp3組態，另一種做法和處理第三主族方式類似，可以通過將settings.ini里的SpherIVgroup設為1來啟用；對于六、七主族，由于每個p軌道都是占據軌道（只是占據數不同），而占據軌道間形狀差異很小（遠小于占據軌道和空軌道的差異），因此直接將p軌道上全部電子在三個p軌道上平攤；對于Sc、Ti、V用的方法和第三主族類似，是將單占據的d軌道向其它方向旋轉并復制（Sc的4s軌道被極化得較厲害，所以也通過旋轉復制的策略削弱它造成的各向異性）；對于Fe、Co、Ni，將d軌道上beta電子都均攤到5個alpha的d軌道上。
注意在過渡金屬原子波函數計算時，所用的理論方法必須能給出正確的4s和3d能級順序，否則Multiwfn沒法正確判斷哪些是d軌而無法正確地球對稱化。一般的泛函如B3LYP沒問題，不能用HF，給出的順序是錯的。Multiwfn對第六、七主族和過渡金屬的處理實際上并沒有使之密度精確球對稱化，但是密度的各向已經已經被削弱到可以忽略的程度。如果不希望Multiwfn自動對元素的波函數文件進行球對稱化處理，可以將settings.ini里的ispheratm設為0來關閉之。
注1：關于應當用什么泛函生成弱相互作用體系的密度（假設都使用相同的可靠的幾何結構），我想在這里表達一下我的觀點。目前的泛函全都是近似的，而且不少都沒有清楚的物理意義。在完全真實的密度下用這些泛函得不到真實的能量，而能夠得到較好的能量/幾何結構的泛函產生的密度未必比那些能量/幾何結構預測得不好的泛函更好。究竟用什么泛函產生密度更合適，這沒有確切結論，如果想獲得在原理上很可靠的密度，那就去做高精度post-HF。但所幸，不同泛函之間產生的密度的差異遠不像能量/幾何結構差異那么大，尤其是對于弱相互作用體系來說。比如，就連公認對氫鍵等弱相互作用體系描述不好的B3LYP也能獲得和MP2很相近的密度（見約化密度梯度分析弱相互作用的原文JACS,132,6498）。所以，第五節就用B3LYP來獲得水二聚體的密度。M06-2X對弱相互作用的能量/結構預測上表現得比B3LYP好很多，會有不少人認為使用M06-2X產生密度更合理，表面上來看用它也的確不會有不妥之處，只不過并不能直接論證它的密度會比B3LYP更接近真實。而對于DFT-D類的泛函，對弱相互作用也普遍表現得較好，這是由于它們在原有泛函能量表達式中額外引入了經驗校正項的結果，而單電子有效勢算符并沒有改變。故產生的密度實際上和原有泛函是一致的。（有一些DFT-D泛函在引入經驗校正項的時候還重新擬合了原有泛函中的參數，這類泛函在原理上我不建議使用。因為經驗校正項不光描述遠程相關作用以改進弱相互作用預測能力，也會描述一定程度中程的相關作用，而原有泛函也描述了中程相關。因此重新擬合出的參數會弱化原有泛函展現的中程相關，這將可能導致產生的密度并沒有充分體現出電子相關效應的影響，還不如原有泛函的密度更可靠）
注2：如果選了-1，就會得到promolecule屬性，它是體系中每個原子在自由狀態下屬性的疊加。也就是說，變形屬性=體系實際屬性-promolecule屬性。如果實空間函數選的是電子密度，就將得到promolecule密度，它代表了分子結構剛剛形成，但是電子密度還沒來得及弛豫時的密度。生成promolecule屬性也需要生成原子波函數文件，和生成變形屬性經歷的流程是完全一致的。