使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷
使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷
Using Multiwfn to calculate Hirshfeld(-I), CM5 and MBIS atomic charges for periodic systems
文/Sobereva@北京科音 2024-Jun-8
1 前言
原子電荷（atomic charge）即原子所帶的凈電荷，對于討論化學體系中原子的帶電狀態有重要意義，作為原子的一種描述符它也和原子的很多屬性有密切聯系。強烈建議閱讀《一篇深入淺出、完整全面介紹原子電荷的綜述文章已發表！》（http://www.shanxitv.org/714）里提到的筆者寫的原子電荷的綜述以全面了解原子電荷。
強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）支持非常豐富的原子電荷計算方法，其中多數已經支持周期性體系。Multiwfn支持的方法中有一類是基于模糊式原子空間劃分，也就是每個原子對應一個平滑的權重函數，權重函數與體系的電子密度的乘積的全空間積分就是這個原子帶的電子數（原子的布居數），原子的核電荷數減去它就是原子電荷。Multiwfn支持的這類方法包括Hirshfeld、Hirshfeld-I、MBIS、Becke，以及對Hirshfeld電荷做后校正得到的ADCH和CM5電荷。從2024-May-25更新的Multiwfn版本開始，Hirshfeld、CM5、Hirshfeld-I、MBIS都已經支持了周期性體系，本文將具體介紹怎么用Multiwfn與非常流行、高效且免費的第一性原理程序CP2K相結合非常方便快速地計算這些原子電荷。Multiwfn計算這些電荷的功能是普適的，并不限于結合CP2K，也可以基于其它第一性原理程序產生的體系的電子密度的cube文件來算，還支持VASP產生的記錄電子密度的CHGCAR文件。
注：Multiwfn更老的一些版本也支持基于周期性波函數計算Hirshfeld和CM5電荷，但那時候用的算法對于周期性體系效率極低，本文介紹的做法與之有天壤之別。
下文第2節簡要介紹一下Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷的基本特點，第3節介紹用Multiwfn+CP2K計算它們的具體方法，第4節給出具體例子，同時還會介紹怎么計算片段電荷。使用Multiwfn計算原子電荷在發表文章時記得需要按Multiwfn啟動時的提示恰當引用Multiwfn。
Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，如果對Multiwfn不了解的話建議看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。如果你用CP2K創建給Multiwfn算原子電荷用的輸入文件的話，筆者假定你對CP2K的使用已經有基本且正確的了解，不具備這些知識的話非常建議通過北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）完整、系統地學習一遍。本文的例子利用到Multiwfn創建CP2K的輸入文件，相關介紹見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。
2 Hirshfeld、Hirshfeld-I、CM5和MBIS等原子電荷簡介
Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷在Multiwfn手冊3.9節的相應小節里都有很詳細的介紹，因此這里不介紹細節，只是簡單說一下它們的特點。
Theor. Chim. Acta (Berl.), 44, 129 (1977)中提出的Hirshfeld電荷是被廣為使用的原子電荷，計算方式簡單，物理意義清楚，在量子化學研究孤立體系和第一性原理研究周期性體系的領域中都用得很多。Hirshfeld電荷為人詬病的地方是原子電荷數值整體嚴重偏小，對偶極矩、靜電勢重現性很差，這在我寫的《原子電荷計算方法的對比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）里的對比測試中有充分體現。
JCP, 126, 144111 (2007)提出的Hirshfeld-I是在Hirshfeld基礎上的改進，通過迭代方式不斷更新原子權重函數直到收斂，這使得原子權重函數所描述的原子空間可以響應實際化學環境。其原子電荷數值顯著大于Hirshfeld，對靜電勢的重現性也有明顯改進。代價是Hirshfeld-I需要做迭代，往往幾十輪，計算耗時明顯高于Hirshfeld，而且需要事先提供體系中各種元素的各個氧化態的原子徑向密度信息，實現起來較麻煩（為了實現Hirshfeld-I，我計算了360多個原子的徑向密度文件并在Multiwfn中自帶）。
一定要注意，雖然CP2K程序自己也支持Hirshfeld-I電荷的計算，但我測試發現其結果完全不靠譜，我檢查了其源代碼也確認了其實現方式明顯不符合Hirshfeld-I的標準定義，因此切勿用CP2K自身的Hirshfeld-I電荷計算功能。CP2K算的Hirshfeld電荷在結合SHAPE_FUNCTION DENSITY選項時的結果不那么離譜，但結果和Multiwfn算的也存在不可忽視的差異，應當以Multiwfn的結果為準。
JCTC, 8, 527 (2012)中提出的CM5電荷和我在J. Theor. Comput. Chem., 11, 163 (2012)中提出的ADCH電荷一樣都是對Hirshfeld電荷的后校正，對大部分情況都能使得它的原子電荷變得更大、更符合化學直覺，對靜電勢的重現性也明顯變得更好。ADCH沒有CM5那么強的經驗性、不牽扯一堆經驗參數，因而原理更為理想、更推薦使用。CM5最初是用于分子體系的，如今在應用于晶體方面也已經有了一定探索，例如JCTC, 16, 5884 (2020)將CM5與其它一些原子電荷對比考察了它們對一些固體的電荷分布的描述。CM5還被用于一些計算化學領域的方法，比如《比SMD算溶解自由能更好的溶劑模型uESE的使用》（http://www.shanxitv.org/593）里介紹的uESE溶劑模型依賴于CM5電荷，CM5乘以1.2后適合結合OPLS-AA力場做動力學模擬，見《計算適用于OPLS-AA力場做模擬的1.2*CM5原子電荷的懶人腳本》（http://www.shanxitv.org/585）。
ADCH電荷在分子體系的研究方面已被廣為使用，之前有人問我是否會把ADCH電荷擴展到周期性體系。我目前沒這個打算，因為ADCH電荷的思想對周期性體系往往不適用。ADCH方法將原子偶極矩展開為周圍原子的校正電荷，但對于諸如NaCl這樣的體系，每個原子所處的環境是中心對稱的，因此原子偶極矩為0、校正電荷都為0，故ADCH電荷會與Hirshfeld電荷完全相同。
MBIS (Minimal Basis Iterative Stockholder)電荷于JCTC, 12, 3894 (2016)中提出。類似于Hirshfeld-I，此方法也是通過迭代方式（因此耗時較高）優化原子空間直到收斂。MBIS比Hirshfeld-I的一大好處是不需要利用各個元素各個氧化態的原子徑向密度信息，因此實現起來簡單得多，也更優雅。原文里的測試體現出MBIS比Hirshfeld-I整體更具有優勢。但MBIS經常收斂較慢，對于某些體系甚至需要200多輪才能充分收斂。目前Multiwfn里MBIS最高支持到Rn元素。
根據我的測試，對于大多數固體體系，原子電荷大小的關系是MBIS和Hirshfeld-I最大、誰大不一定，CM5大小介中，Hirshfeld電荷最小。如果你對數值的絕對大小不那么在乎，只是想體現相對大小關系的話，最“樸素”且很便宜的Hirshfeld也可以用。但如果對定量數值關心的話，建議用其它的。CM5、Hirshfeld-I和MBIS用于晶體哪個最理想沒有定論，都可以算一算看看，如果傾向于得到較大數值的結果可以優先考慮Hirshfeld-I和MBIS。如果其中有的原子電荷明顯不合理那就棄掉換別的方法，比如O的電荷如果算出來-2.34那就不應該用，因為顯然O最多只能帶-2電荷。
特別要強調的是，絕對不要把原子電荷與氧化態搞混，也切勿用原子電荷判斷氧化態！氧化態怎么正確地判斷看《使用Multiwfn結合CP2K計算晶體中原子的氧化態》（http://www.shanxitv.org/711）和《使用Multiwfn通過LOBA方法計算氧化態》（http://www.shanxitv.org/362）。
還有很多其它原子電荷計算方法，諸如Mulliken電荷、Lowdin電荷、SCPA電荷、AIM電荷（被個別文章稱為Bader電荷）等，Multiwfn也都支持將它們用于周期性體系，這不屬于本文的范疇。Mulliken電荷結合CP2K里常用的MOLOPT系列和pob系列基組用于周期性體系也說得過去。雖然它是基于希爾伯特空間劃分的原子電荷中最原始、思想最樸素的一種，但實際結果一般還算定性正確，而且計算耗時極低。被用得很多的AIM電荷實際上很糟糕，往往出現違背化學常識的結果，這在前述的《原子電荷計算方法的對比》里已經充分體現了，不建議用。
3 Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷在Multiwfn中的計算方法
3.1 CP2K用戶的情況
首先說對于CP2K用戶怎么計算。Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷都是完全基于電子密度計算的，不直接牽扯到波函數，在Multiwfn中計算它們有兩種可以用的輸入文件：
(1)用CP2K產生的體系的電子密度cube文件作為輸入文件。這種情況在CP2K計算時可以考慮k點，是我最推薦的，耗時非常低。如果不了解cube文件的格式的話，參看《Gaussian型cube文件簡介及讀、寫方法和簡單應用》（http://www.shanxitv.org/125）。注意CP2K有bug（起碼直到我撰文時最新的CP2K 2024.1為止），產生的cube文件里的原子有效電荷都為0，這會導致Multiwfn無法正確計算原子電荷。這里說的原子的有效核電荷是指當前基組描述的原子的電子數，對于贗勢基組它對應價電子數。為此，要么用文本編輯器手動改cube文件中每個原子的有效核電荷，要么把cube文件第一行寫成各個元素的有效核電荷，例如N 5 B 3 Ti 12，這代表當前體系中的N、B、Ti的有效核電荷數分別為5、3、12
(2)用CP2K產生的記錄了體系波函數的molden文件作為輸入文件。詳見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）。這個方法的缺點在于沒法考慮k點，對于小晶胞體系必須先擴胞到足夠大，從而在計算時只考慮gamma點，這顯然使得在CP2K中和Multiwfn中的計算耗時都會比用(1)的方式高很多
用Multiwfn載入上述輸入文件其中一種后，進入主功能7，然后選擇相應選項即可計算Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷。對于molden文件當輸入文件時，Multiwfn還會讓你輸入計算用的格點間距，格點間距越小耗時越高而計算精度越高，通常用0.2 Bohr就可以，是精度和耗時的較好權衡。而使用cube文件做為輸入文件時，計算原子電荷用的格點間距與此文件的格點間距直接對應，CP2K的CUTOFF設得越大則產生的cube文件里的格點間距會越小。
由于Multiwfn計算上述電荷用的是均勻分布的格點，不可能準確積分較重原子的內核區域電子密度，因此上述輸入文件應當是CP2K做GPW計算得到的，也需要用贗勢基組，此時cube文件記錄的電子密度或者molden文件里記錄的波函數只對應于價電子。
順帶一提，如果你要計算Mulliken電荷，用Multiwfn載入上述的molden文件，然后進主功能7，依次選擇5、1即可得迅速到結果。Multiwfn中計算AIM電荷的方法筆者以后會另文說明。
3.2 其它程序用戶的情況
如果你是VASP用戶，可以用記錄晶胞中價電子密度的CHGCAR文件給Multiwfn用于計算前述原子電荷，等同于3.1節說的第2類輸入文件。只要文件名里包含CHGCAR字樣就會被Multiwfn當做是CHGCAR文件來載入。注意為了讓Multiwfn得知各類原子的有效核電荷數，CHGCAR文件的第一行需要以Nval為開頭，后面依次寫上各類元素的名字和有效核電荷，并以空格分隔。比如寫為Nval Na 1 Cl 7，這代表Na和Cl的有效核電荷數分別為1和7。
如果你是其它第一性原理程序的用戶，若程序可以產生記錄晶胞中價電子密度的cube文件（即等同于3.1節說的第2類輸入文件），也都可以用于給Multiwfn計算前述的原子電荷。
之后在Multiwfn中的計算過程和CP2K用戶沒任何區別。
4 實例
下面給出Multiwfn結合CP2K計算Hirshfeld、Hirshfeld-I、CM5和MBIS原子電荷的具體例子。相關的文件都可以在http://www.shanxitv.org/attach/712/file.zip中得到（4.2、4.3和4.5節的cube/molden文件除外，因為文件過大）。
4.1 TiO2晶體
此例對TiO2的金紅石晶型計算前述的原子電荷。它的實驗晶體結構的cif文件是本文文件包里的TiO2-Rutile.cif，結構圖如下所示
先用CP2K計算出這個晶胞的電子密度的cube文件，用于之后Multiwfn做原子電荷計算。首先用Multiwfn創建CP2K輸入文件，啟動Multiwfn并載入TiO2-Rutile.cif，然后輸入
cp2k //創建CP2K的輸入文件
[回車] //產生的文件名用默認的TiO2-Rutile.inp
-3 //產生cube文件
1 //電子密度
8 //設置k點
4,4,6 //與各個晶胞邊長相乘皆為18埃左右，夠用了
0 //產生輸入文件，計算級別為默認的PBE/DZVP-MOLOPT-SR-GTH
用CP2K運行剛得到的TiO2-Rutile.inp，馬上就算完了，得到了TiO2-Rutile-ELECTRON_DENSITY-1_0.cube文件，在本文的文件包里已經提供了。用文本編輯器打開此文件，把第一行改為Ti 12 O 6并保存。這代表Ti和O的有效核電荷分別為12和6，直接對應于輸入文件里它們分別用的DZVP-MOLOPT-SR-GTH-q12和DZVP-MOLOPT-SR-GTH-q6基組名里面的q后面的值。
啟動Multiwfn，輸入
TiO2-Rutile-ELECTRON_DENSITY-1_0.cube //寫此文件的實際路徑
y //從第一行載入各個元素有效核電荷
7 //布居分析與原子電荷計算
1 //Hirshfeld電荷
立刻看到如下Hirshfeld電荷計算結果。由于格點間距不是無窮小，因此等價的原子算出來的原子電荷也會有輕微的數值差異，可以自己手動取平均。如果用更大的CUTOFF，則cube文件的格點間距會更小，不等價性會更低。
Final atomic charges:
Atom 1(Ti): 0.59349459
Atom 2(O ): -0.29461903
Atom 3(O ): -0.29472601
Atom 4(Ti): 0.58363165
Atom 5(O ): -0.29385236
Atom 6(O ): -0.29385236
現在Multiwfn還問你是否把原子電荷導出為chg文件，這里選n不導出。
接下來計算CM5電荷。在Multiwfn主功能7里選擇16，馬上就得到了結果
---------- CM5 charges ----------
Atom: 1Ti CM5 charge: 1.348393 Hirshfeld charge: 0.593495
Atom: 2O CM5 charge: -0.672093 Hirshfeld charge: -0.294619
Atom: 3O CM5 charge: -0.672191 Hirshfeld charge: -0.294726
Atom: 4Ti CM5 charge: 1.338580 Hirshfeld charge: 0.583632
Atom: 5O CM5 charge: -0.671306 Hirshfeld charge: -0.293852
Atom: 6O CM5 charge: -0.671306 Hirshfeld charge: -0.293852
Summing up all CM5 charges: 0.00007649
可見CM5電荷的大小明顯大于Hirshfeld電荷，看起來更合理，而Hirshfeld電荷值明顯偏小。
下面再演示計算Hirshfeld-I電荷。把Multiwfn自帶的examples目錄下的atomrad目錄挪到當前目錄下，這樣Multiwfn計算Hirshfeld-I電荷時就會自動用此目錄下的各種元素各個價態的徑向電子密度。之后在主功能7里選15，再選1用默認的設置進行Hirshfeld-I電荷計算，經過16輪迭代得到了結果：
Final atomic charges:
Atom 1(Ti): 3.01789464
Atom 2(O ): -1.50888824
Atom 3(O ): -1.50886131
Atom 4(Ti): 3.01787620
Atom 5(O ): -1.50897240
Atom 6(O ): -1.50897240
可見Hirshfeld-I電荷的大小不僅顯著大于Hirshfeld電荷，比CM5電荷還要大很多。雖然數值很大，但完全在合理范圍內，畢竟都沒有大過原子的氧化態。
接下來再計算MBIS電荷。在主功能7里選擇20，再選1用默認的設置進行MBIS電荷計算，經過19輪迭代得到了結果：
Final atomic charges:
Atom 1(Ti): 1.73473636
Atom 2(O ): -0.86733169
Atom 3(O ): -0.86731253
Atom 4(Ti): 1.73467041
Atom 5(O ): -0.86734303
Atom 6(O ): -0.86734303
對于當前例子，MBIS的數值介于CM5和Hirshfeld-I之間。但也有很多反例，MBIS也經常比Hirshfeld-I電荷還大不少，后文有例子。
以上原子電荷的計算順序是隨意的，不要誤以為必須按上述順序計算，你想算哪個就算哪個即可。
4.2 AlN晶體
這一節計算AlN晶體的原子電荷。完全可以按照和上一節相同的做法基于電子密度的cube文件來計算，但本節演示一下基于超胞的molden文件的計算過程。AlN原胞的cif文件是本文文件包里的AlN.cif，將之載入Multiwfn，然后輸入以下命令創建CP2K輸入文件
cp2k //創建CP2K輸入文件
AlN_443.inp //產生的輸入文件名
-11 //進入幾何操作界面
19 //構造超胞
4 //a方向擴胞成原本的4倍
4 //b方向擴胞成原本的4倍
3 //c方向擴胞成原本的3倍
-10 //返回
-2 //要求產生molden文件
0 //產生輸入文件
用CP2K運行剛產生的AlN_443.inp，得到AlN_443-MOS-1_0.molden。在此文件開頭插入以下內容定義晶胞和各個元素的價電子數
[Cell]
12.44400000 0.00000000 0.00000000
-6.22200000 10.77682012 0.00000000
0.00000000 0.00000000 14.93400000
[Nval]
Al 3
N 5
啟動Multiwfn，載入AlN_443-MOS-1_0.molden，然后輸入
7 //布居分析與原子電荷計算
16 //CM5電荷
[回車] //使用均勻分布的格點計算
[回車] //使用默認的0.2 Bohr格點間距。格點間距越小結果越精確但計算越耗時
在i9-13980HX CPU上24核并行，50秒就算完了，結果為
---------- CM5 charges ----------
Atom: 1Al CM5 charge: 0.467945 Hirshfeld charge: 0.302813
Atom: 2N CM5 charge: -0.467945 Hirshfeld charge: -0.302812
Atom: 3Al CM5 charge: 0.467945 Hirshfeld charge: 0.302812
Atom: 4N CM5 charge: -0.467946 Hirshfeld charge: -0.302814
...略
Al的電負性很小，N的電負性很大，因此二者的原子電荷理應相差懸殊。當前的Hirshfeld電荷明顯偏小。CM5電荷的數值雖然更大，但還是顯得有些偏低。
下面再計算Hirshfeld-I電荷。在主功能7里輸入
15 //Hirshfeld-I
1 //開始計算
[回車] //使用默認的0.2 Bohr格點間距
經過37輪收斂，結果如下。可見Hirshfeld-I電荷的數量級比較合理，體現出AlN中成鍵的極強離子性
Atom 1(Al): 1.80317384
Atom 2(N ): -1.80317397
Atom 3(Al): 1.80317257
Atom 4(N ): -1.80317539
...略
最后再計算下MBIS電荷。在主功能7里輸入
20 //MBIS
1 //開始計算
[回車] //使用默認的0.2 Bohr格點間距
對此體系MBIS收斂很慢，經過190輪才達到收斂，在i9-13980HX上24核并行算了5分鐘，在雙路7R32 96核并行下算了不到一分鐘。結果如下。可見對此體系，MBIS電荷的數值很大，和氧化態幾乎正好是一樣的。所以，如前所述，若你就是想算出來比較大的原子電荷，可以先考慮MBIS，但也有一些例外，后面會提到。
Atom 1(Al): 2.99599978
Atom 2(N ): -2.99601488
Atom 3(Al): 2.99600129
Atom 4(N ): -2.99599935
4.3 MOF-5晶體
MOF-5晶體含有Zn、O、C、H。其cif文件是本文文件包里的MOF-5.cif。由于此體系的晶胞邊長約26埃，因此計算電子密度格點數據時不用考慮k點。按照3.1節的做法，基于PBE/DZVP-MOLOPT-SR-GTH級別的電子密度cube文件（第一行改為Zn 12 O 6 C 4 H 1），用Multiwfn計算各種原子電荷，結果如下。此體系中O有兩類，分別給出，C則處于較多不同的化學環境，其電荷就不給出了。
Hirshfeld：
Zn=0.372 O=-0.332/-0.200 H=0.056
CM5：
Zn=0.644 O=-0.584/-0.305 H=0.116
Hirshfeld-I（36輪收斂，2*7R32機子上96核并行4分多鐘算完）
Zn=1.568 O=-1.917/-0.746 H=0.139
MBIS（39輪收斂，2*7R32機子上96核并行3分多鐘算完）
Zn=1.546 O=-1.662/-0.782 H=0.190
對這個體系，MBIS和Hirshfeld-I相符得很好，而且Zn的電荷大小很合理。Hirshfeld電荷再次數值最小、CM5比之大一些。
4.4 其它：MoS2、BaTiO3、MgF2、NaCl
注意MBIS電荷并非總是較大。例如《使用Multiwfn結合CP2K計算晶體中原子的氧化態》（http://www.shanxitv.org/711）里的單層MoS2的例子，在PBE/DZVP-MOLOPT-SR-GTH級別的波函數下，Mo的電荷的計算結果是Hirshfeld=0.330，CM5=0.827，Hirshfeld-I=1.987，MBIS=-0.059，Mulliken=-0.771。對這個體系來說CM5和Hirshfeld-I顯得靠譜，MBIS電荷接近0有點說不通，Mulliken電荷明顯有誤導性。
BaTiO3在PBE/DZVP-MOLOPT-SR-GTH級別下的原子電荷計算結果如下。Hirshfeld還是嚴重偏小。Hirshfeld-I大得離譜，Ba的原子電荷都超過氧化態+2了。相較之下CM5和MBIS的結果比較靠譜，正好二者算的O的原子電荷差不多，差異主要在Ba和Ti的電荷的相對大小上。
Atom: 1Ba CM5 charge: 1.103253 Hirshfeld charge: 0.404357
Atom: 2Ti CM5 charge: 1.200467 Hirshfeld charge: 0.537212
Atom: 3O CM5 charge: -0.767963 Hirshfeld charge: -0.310211
Atom: 4O CM5 charge: -0.767892 Hirshfeld charge: -0.315693
Atom: 5O CM5 charge: -0.767892 Hirshfeld charge: -0.315693
Hirshfeld-I
Atom 1(Ba): 2.37169239
Atom 2(Ti): 3.02878350
Atom 3(O ): -1.78799416
Atom 4(O ): -1.80625420
Atom 5(O ): -1.80625434
MBIS：
Atom 1(Ba): 0.57391985
Atom 2(Ti): 1.69416825
Atom 3(O ): -0.74678425
Atom 4(O ): -0.76066530
Atom 5(O ): -0.76066537
再看兩個典型離子化合物的情況，都在PBE/DZVP-MOLOPT-SR-GTH級別下計算。對NaCl計算的Na的原子電荷為：Hirshfeld=0.217，CM5=0.443，Hirshfeld-I：1.053，MBIS：1.043，Mulliken：0.596
對MgF2計算的Mg的原子電荷為：
Hirshfeld=0.414，CM5=0.828，Hirshfeld-I：2.048，MBIS：1.888，Mulliken：1.439
可見對于強離子性化合物，Hirshfeld-I和MBIS電荷都在氧化態附近，并且由于方法定義的缺陷，往往電荷的大小還輕微大于氧化態的大小（但如果用大核贗勢，對Na和Mg只描述最外層的3s價電子，則不會有這個問題）。其它原子電荷的大小關系為Mulliken>CM5>Hirshfeld。
以MgF2為例，這里順帶一提CP2K算的Hirshfeld-I電荷明顯不對，千萬不要用。在CP2K輸入文件的&DFT/&PRINT里加上以下內容要求輸出Hirshfeld-I電荷后，CP2K直接給出的Mg的Hirshfeld-I電荷為0.522，和Multiwfn算的相距甚遠。
&HIRSHFELD
SHAPE_FUNCTION DENSITY
SELF_CONSISTENT T
&END HIRSHFELD
4.5 NaCl板吸附CO的片段電荷
NaCl板吸附CO是《使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線》（http://www.shanxitv.org/638）里的一個例子體系，這里基于這個模型計算一下原子電荷和CO片段的電荷。CP2K優化這個體系后產生的NaCl_CO-1.restart文件在本文的文件包里提供了，將之載入Multiwfn，然后輸入
cp2k //創建CP2K輸入文件
NaCl_CO.inp //產生的輸入文件名
-7 //設置周期性
XY
-3 //要求產生cube文件
1 //電子密度
-2 //產生molden文件（用于之后算Mulliken電荷的目的）
0 //產生輸入文件
計算完畢后，把NaCl_CO-ELECTRON_DENSITY-1_0.cube文件第一行改為Na 9 Cl 7 C 4 O 6。先用Hirshfeld方法各個原子的電荷和CO的片段電荷。啟動Multiwfn，載入此cube文件，然后輸入
7 //布居分析與原子電荷計算
-1 //定義片段
109,110 //CO的原子序號
1 //Hirshfeld電荷
結果為
...略
Atom 109(C ): 0.13622804
Atom 110(O ): -0.03514734
Fragment charge: 0.10108070
Fragment population: 9.89891930
即CO的Hirshfeld片段電荷為0.101。
之后再選擇CM5方法，結果為
...略
Atom: 109C CM5 charge: 0.147306 Hirshfeld charge: 0.136228
Atom: 110O CM5 charge: -0.081972 Hirshfeld charge: -0.035147
Summing up all CM5 charges: 0.00510030
Fragment charge: 0.06533361
之后再選擇Hirshfeld-I方法，結果為
Atom 109(C ): 0.09790054
Atom 110(O ): -0.08744136
Fragment charge: 0.01045918
之后再選擇MBIS方法，結果為
...略
Atom 109(C ): 0.19642171
Atom 110(O ): -0.15937966
Fragment charge: 0.03704205
再來算Mulliken電荷。用Multiwfn載入加入了恰當的[Cell]和[Nval]信息的CP2K計算后產生的NaCl_CO-MOS-1_0.molden文件，輸入
7 //布居分析與原子電荷計算
-1 //定義片段
109,110 //CO的原子序號
5 //Mulliken電荷
1 //在屏幕上輸出結果
結果為
...略
Atom 109(C ) Population: 3.60671601 Net charge: 0.39328399
Atom 110(O ) Population: 6.23484007 Net charge: -0.23484007
Total net charge: -0.00000086
Fragment charge: 0.15844392
可見不同方法計算的結果在定量上有一定差異，但都有共性，即CO整體帶少量正電，即電子從CO往NaCl板上轉移了一些。并且C和O分別帶正電和負電，且原子電荷的絕對值都是C大于O。
5 總結
本文簡要介紹了四種知名的基于模糊式原子空間劃分計算原子電荷的方法的特點，包括Hirshfeld、Hirshfeld-I、CM5和MBIS，并介紹了它們在Multiwfn中的計算方法，并且給了很多具體例子，以使得讀者能了解如何使用Multiwfn結合CP2K或其它第一性原理程序很容易地計算它們。同時還說明了怎么直接計算片段電荷，這是定量討論體系中兩部分間電子轉移的最嚴格的方式。推薦將這些原子電荷應用于實際問題的研究當中以考察固體與表面體系的電荷分布情況。記得屆時應恰當引用Multiwfn啟動時提示的Multiwfn原文以及相應的原子電荷的原文。
通過本文的一些體系的對比可以明顯看到，Hirshfeld電荷原理上最簡單，但明顯缺點是整體嚴重偏小，CM5通常比它大不少，而且電荷的穩定性不錯。雖然結合化學直覺來看，CM5電荷的絕對大小仍然往往能偏小，但至少也不會明顯離譜，總比Hirshfeld電荷更建議使用。Hirshfeld-I和MBIS電荷普遍比較大，誰更大不一定，至少幾乎總比CM5還大不少，從絕對大小來看往往比CM5更合理，但穩健性比CM5弱一些，對極個別體系可能會表現得不符合常識。對于離子性特征很強的化合物，Hirshfeld-I和MBIS電荷往往接近氧化態，由于方法并不完美，有時候原子電荷甚至比氧化態還要大一點點。當然，沒有任何原子電荷計算方法是絕對完美的，本文介紹的這些原子電荷在實際當中都可以用，可以根據它們的特點和實際結果決定選用，在橫向對比時方法必須統一。
順帶一提，對于計算固體表面或者多孔物質的原子電荷以用于基于力場的動力學模擬的目的，我目前最推薦CP2K算REPEAT電荷，北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里講“電子結構的分析”的部分做了專門的介紹。但REPEAT電荷需要在真空區定義擬合點，因而無法用于諸如本文所示的AlN、TiO2等致密的固體，而且對于遠離真空區的原子也沒法擬合出質量靠譜的原子電荷。所以REPEAT電荷有特定用處，但不是像本文介紹的那些原子電荷一樣屬于普適性的原子電荷計算方法。