使用Multiwfn做周期性體系的atom-in-molecules (AIM)拓撲分析
使用Multiwfn做周期性體系的atom-in-molecules (AIM)拓撲分析
Using Multiwfn to perform atom-in-molecules (AIM) topology analysis for periodic systems
文/Sobereva@北京科音 2024-Jul-3
0 前言
Bader提出的非常知名的atom-in-molecules (AIM)理論在分析孤立體系和周期性體系的電子結構方面有非常廣泛的應用，相關資料匯總見《AIM學習資料和重要文獻合集》（http://bbs.keinsci.com/thread-362-1-1.html）。AIM分析框架中最常見的分析手段之一是電子密度的拓撲分析，也常被稱為AIM拓撲分析，例如根據這種拓撲分析得到的鍵臨界點（BCP）的位置可以用來討論成鍵和弱相互作用特征，參考《Multiwfn支持的分析化學鍵的方法一覽》（http://www.shanxitv.org/471）和《Multiwfn支持的弱相互作用的分析方法概覽》（http://www.shanxitv.org/252）里的相關介紹。非常強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）很早以前就已經支持了結合Gaussian、ORCA等量子化學程序產生的波函數文件對孤立體系做AIM拓撲分析，見《使用Multiwfn做拓撲分析以及計算孤對電子角度》（http://www.shanxitv.org/108）和Multiwfn手冊4.2節的大量例子。如今Multiwfn也已支持了對周期性體系做AIM拓撲分析，本文的目的是對此進行演示，以使得讀者可以輕松地舉一反三充分將這種手段應用于實際研究中。使用Multiwfn做AIM拓撲分析在發表文章時記得需要按照Multiwfn啟動時的提示恰當引用Multiwfn的原文。
讀者務必使用2024-Jun-16及以后更新的Multiwfn版本，否則情況與本文所述不符。Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載。不了解Multiwfn者建議閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。本文要用到非常流行、高效且免費的第一性原理程序CP2K，不熟悉者推薦通過“北京科音CP2K第一性原理計算培訓班”（http://www.keinsci.com/workshop/KFP_content.html）系統性學習。本文要使用Multiwfn創建CP2K輸入文件，這在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）里有簡要介紹。本文用的CP2K是2024.1版。
下面的例子涉及到的所有文件都可以在http://www.shanxitv.org/attach/717/file.rar里得到。
1 實例1：NaCl晶體
這一節對NaCl做AIM拓撲分析，要使用CP2K產生的molden文件作為Multiwfn的輸入文件。molden文件的產生方法在《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）里有詳細說明。由于考慮k點時CP2K沒法產生molden文件，因此必須使用足夠大的超胞模型并只考慮gamma點。本文文件包里NaCl.cif是NaCl單胞的晶體結構文件，基于它我們創建CP2K算單點任務的超胞的輸入文件，同時得到molden文件。
用Multiwfn載入NaCl.cif，然后輸入
cp2k //創建CP2K輸入文件的界面
[回車] //產生的輸入文件名用默認的NaCl.inp
-2 //要求產生molden文件
-11 //幾何操作界面
19 //構造超胞
2 //第1個方向平移復制為原先的2倍
2 //第2個方向平移復制為原先的2倍
2 //第3個方向平移復制為原先的2倍
-10 //返回
0 //產生輸入文件
現在當前目錄下就有了NaCl.inp，對應的是PBE/DZVP-MOLOPT-SR-GTH計算級別，對此體系是合適的。當前晶胞邊長為11.28埃，對于絕緣體來說只考慮gamma點基本夠用。用CP2K運行NaCl.inp，得到NaCl-MOS-1_0.molden。用文本編輯器打開NaCl.inp，將以下內容插入文件開頭來定義晶胞信息從而使得Multiwfn在分析時能考慮周期性，并且定義Na和Cl的價電子數分別為9和7（不了解為什么這么設置的人看前述的《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》）。
[Cell]
11.28112000 0.00000000 0.00000000
0.00000000 11.28112000 0.00000000
0.00000000 0.00000000 11.28112000
[Nval]
Na 9
Cl 7
下面開始做AIM拓撲分析。啟動Multiwfn，載入NaCl-MOS-1_0.molden，然后輸入
2 //拓撲分析（默認是對電子密度做拓撲分析）
2 //以每個原子核位置為初猜搜索臨界點
當前從屏幕上會看到找到了64個(3,-3)臨界點，分別對應64個原子核位置的電子密度極大點。之所以雖然當前用的是贗勢基組，在原子核位置還能找到電子密度極大點，是因為Multiwfn自動對這些原子指認了合適的EDF信息來描述內核電子，詳見《在贗勢下做波函數分析的一些說明》（http://www.shanxitv.org/156）。
然后再輸入
3 //以每兩個原子間連線中點作為搜索臨界點的初猜位置
4 //以每三個原子的中心位置作為搜索臨界點的初猜位置
0 //觀看結果
此時文本窗口顯示了各種臨界點的數目，如下所示，并且提示Poincare-Hopf關系已經滿足了，因此大概率所有臨界點都已經找全了。
The number of critical points of each type:
(3,-3): 64, (3,-1): 384, (3,+1): 384, (3,+3): 64
Poincare-Hopf verification: 64 - 384 + 384 - 64 = 0
Fine, Poincare-Hopf relationship is satisfied, all CPs may have been found
在Multiwfn的圖形窗口可以看到下圖，棕紅色大球是Na，綠色大球是Cl，橙色、黃色、綠色小球分別是電子密度的(3,-1)、(3,+1)和(3,+3)臨界點，而(3,-3)臨界點都被原子球擋住了。在Multiwfn手冊3.14節對Multiwfn的拓撲分析功能有完整、非常詳細的說明，建議閱讀。
下面再把AIM拓撲分析里常涉及的鍵徑產生出來。在圖形界面右上角點Return關閉之，然后選8產生(3,-3)和(3,-1)之間的拓撲路徑，也即鍵徑。在筆者的i9-13980HX CPU上24核并行很快就產生完了。之后再選0進入圖形界面觀看，看到下圖，可見所有鍵徑（棕色細線）都產生了出來，既連接相鄰的Na和Cl，也連接每一對臨近的兩個Cl。
默認情況下Multiwfn把處于晶胞邊緣的所有鏡像原子和臨界點也顯示了出來，如果不想顯示它們的話（也即只顯示晶胞中唯一的那些），在圖形界面頂端依次選擇Other settings - Toggle showing all boundary atoms和Toggle showing all boundary CPs and paths。在文本窗口也會看到當前設置的狀態。
下面來考察一下Na和Cl之間的BCP的屬性。用選項0的圖形界面右邊的選項把其它臨界點關閉顯示而只顯示(3,-1)臨界點，選擇CP labels復選框要求顯示臨界點序號，關閉鍵徑顯示并要求顯示原子，然后就能看到下圖，這是當前圖像的一個邊緣區域（點擊save picture按鈕可在當前目錄下得到高分辨率圖像文件，下圖是剪切出的一部分）
可見189、512、238等序號的臨界點都對應Na-Cl的BCP，考察哪個都可以，都是等價的。點擊Return按鈕關閉圖形窗口，然后選擇7進入臨界點屬性計算界面，然后再輸入臨界點序號189，此時就看到了這個臨界點處各種函數的值，如下所示。這些函數在Multiwfn手冊2.6節都有簡要介紹，且在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/workshop/WFN_content.html）里有十分全面深入詳細的講解。
CP type: (3,-1)
Density of all electrons: 0.1144132659E-01
Density of Alpha electrons: 0.5720663296E-02
Density of Beta electrons: 0.5720663296E-02
Spin density of electrons: 0.0000000000E+00
Lagrangian kinetic energy G(r): 0.1168592222E-01
G(r) in X,Y,Z: 0.7296255010E-03 0.1022667133E-01 0.7296253875E-03
Hamiltonian kinetic energy K(r): -0.3025433987E-02
Potential energy density V(r): -0.8660488230E-02
Energy density E(r) or H(r): 0.3025433987E-02
Laplacian of electron density: 0.5884598270E-01
Electron localization function (ELF): 0.1993303578E-01
Localized orbital locator (LOL): 0.1249064363E+00
...略
如我在全面提到的《Multiwfn支持的分析化學鍵的方法一覽》文中所述，BCP位置上電子密度拉普拉斯函數和能量密度都為正是離子鍵的典型特征。當前BCP這兩個值分別為0.0588 a.u.和0.00302 a.u.，都為正值，完全體現出Na-Cl是離子鍵的實事。此外，當前BCP的電子密度僅為0.01144 a.u.、ELF僅為0.0199，BCP處很小的電子密度和ELF值也是典型的離子鍵的普遍共性。
我在《使用Multiwfn+VMD快速地繪制高質量AIM拓撲分析圖（含視頻演示）》（http://www.shanxitv.org/445）介紹了一種極好的將Multiwfn和VMD聯用快速、理想地繪制AIM拓撲分析圖的做法，比在Multiwfn里直接顯示的效果更好，而且視角可以自由旋轉，沒看過的話務必先看一下。對于當前體系，也可以像此文一樣用Multiwfn文件包里的examples\scripts\AIM.bat和AIM.txt來作圖，得到的圖像如下。
上圖效果明顯不好，顯得殘缺不全，這是因為此時顯示的原子、臨界點、鍵徑都只有晶胞里唯一的那些，而處于棱和面上的邊界原子、臨界點、鍵徑的周期鏡像都沒有顯示出來。為了避免這個問題，對于當前體系有原子、臨界點、鍵徑處于棱和面上的情況，應當改用examples\scripts\目錄下的AIM+boundary.bat和AIM+boundary.txt來繪制，其用法與AIM.bat和AIM.txt的組合完全一致。只不過前者傳遞給Multiwfn的指令要求Multiwfn導出mol.pdb、CPs.pdb和paths.pdb文件時把邊界的鏡像也都納入進去。利用AIM+boundary.bat和AIM+boundary.txt來繪制的結果如下，可見此時的效果就很好了，和前面在Multiwfn里看到的完全相同。注意此時這三個pdb文件里只有序號靠前的那些才是晶胞中唯一的那些。
2 層狀COF
下面再舉個例子，用Multiwfn對一個層狀的共價有機框架化合物（COF）體系做AIM拓撲分析。此體系的cif文件是本文文件包里的COF-12000N2.cif，其結構圖如下所示。
由于此體系的晶胞尺寸不小，計算時可以直接對原胞做只考慮gamma點的計算來得到molden文件。如《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）所說，由于X光衍射測的晶體結構中氫的原子位置通常不準確，最好做分析前先凍結重原子而優化一下氫原子的坐標，但此例為了省事就跳過這一步了。如果你要用AIM拓撲分析考察氫鍵，由于氫的位置是重中之重，那顯然是絕對不能略過這一步的。
當前體系用PBE結合CP2K的GTH贗勢基組計算完全可以，而此例作為演示，這回使用PBE結合pob-TZVP-rev2全電子基組計算。啟動Multiwfn，載入COF-12000N2.cif，然后輸入
cp2k //創建CP2K輸入文件的界面
[回車] //產生的輸入文件名用默認的COF-12000N2.inp
-2 //要求產生molden文件
2 //選擇基組
15 //pob-TZVP-rev2
4 //開OT，對當前用的基組通常比對角化明顯更容易收斂
0 //產生輸入文件
用CP2K計算COF-12000N2.inp，然后將以下晶胞設置信息插入到計算得到的COF-12000N2-MOS-1_0.molden文件的開頭
[Cell]
22.55600000 0.00000000 0.00000000
-11.27800000 19.53406901 0.00000000
0.00000000 0.00000000 6.80000000
啟動Multiwfn，載入COF-12000N2-MOS-1_0.molden，然后輸入
2 //拓撲分析（默認是對電子密度做拓撲分析）
2 //以每個原子核位置為初猜位置搜索臨界點
3 //以每兩個原子間連線中點作為搜索臨界點的初猜位置
4 //以每三個原子的中心位置作為搜索臨界點的初猜位置
8 //產生(3,-3)到(3,-1)的拓撲路徑
0 //觀看結果
現在看到下圖。可見臨界點和鍵徑都成功地得到了，不僅有化學鍵對應的BCP和鍵徑，COF內氫鍵N-H...O作用以及層之間相互作用相對應的BCP和鍵徑都有。而且不僅有晶胞內兩層COF之間相互作用的BCP和鍵徑，當前晶胞里的COF和相鄰鏡像晶胞里的COF相互作用對應的BCP和鍵徑也都有，顯然周期性完全充分考慮了。
雖然Multiwfn的文本窗口提示當前并未符合Poincare-Hopf關系，但由于感興趣的臨界點和鍵徑都有了，所以就沒必要進一步搜索了。
再利用前述名為AIM+boundary.bat的Windows批處理腳本結合AIM+boundary.txt里記錄的指令，用Multiwfn+VMD繪制體系結構+臨界點+鍵徑圖。然后再在VMD的文本窗口里輸入pbc box命令顯示盒子，VMD main界面里選擇Display - Orthographic用正交視角。此時VMD顯示的圖像如下，可見非常理想！
順帶一提，對于當前體系要想完整展現層間相互作用的話，最理想的做法是用我提出的IGMH，見《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）和《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）。
3 對Ag晶胞基于電子密度格點數據做AIM拓撲分析
Multiwfn不僅可以基于CP2K產生的波函數文件做AIM拓撲分析，還可以基于從cube等文件里載入的電子密度的格點數據做拓撲分析，此時各個位置的電子密度、梯度和Hessian矩陣基于格點數據以B-spline方式插值得到。由于這個特性，Multiwfn做周期性體系的AIM拓撲分析不限于CP2K用戶，諸如從VASP產生的CHGCAR文件里讀取電子密度格點數據做AIM拓撲分析也可以。Multiwfn支持哪些格點數據文件格式見手冊2.5節的說明。而且CP2K用戶還可以在考慮k點的情況下對原胞計算得到電子密度cube文件結合Multiwfn做AIM拓撲分析，而不必在晶胞較小時非得弄成超胞模型再算。但基于電子密度格點數據做AIM拓撲分析的關鍵不足是對臨界點只能得到電子密度及其導數信息以及一些與之直接相關的量，而依賴于波函數才能計算的如動能密度、ELF等函數就沒法計算了。另外，如果電子密度格點數據是在贗勢下計算產生的，基于它插值出的電子密度函數中，原子核處將沒有電子密度極大點，也因此無法產生鍵徑。
此例演示用Multiwfn對Ag的單胞基于電子密度cube文件做AIM拓撲分析的過程。使用的級別是PBE/TZVPP-MOLOPT-PBE-GTH-q19，這個基組是CP2K 2024.1的data目錄下的BASIS_MOLOPT_UZH基組文件里的。之所以不用常用的BASIS_MOLOPT里的Ag的基組，在于那里面Ag的基組都是-q11的，即只描述11個價電子。使用贗勢的情況下做AIM拓撲分析時，被描述的價電子越多，即使用越小核贗勢，結果越準確、越可能接近全電子計算的結果。
啟動Multiwfn，載入本文文件包里的Ag單胞的晶體結構文件Ag.cif，然后輸入
cp2k //創建CP2K輸入文件的界面
[回車] //產生的輸入文件名用默認的Ag.inp
-3 //要求產生cube文件
1 //電子密度
6 //開smearing
8 //設置k點
10,10,10
0 //產生輸入文件
把Ag.inp里BASIS_SET_FILE_NAME設為BASIS_MOLOPT_UZH，POTENTIAL_FILE_NAME設為POTENTIAL_UZH。Ag的BASIS_SET設為TZVPP-MOLOPT-PBE-GTH-q19，POTENTIAL設為GTH-PBE-q19。然后用CP2K計算，之后會得到記錄晶胞內價電子密度的Ag-ELECTRON_DENSITY-1_0.cube。如果不了解cube格式的話，看《Gaussian型cube文件簡介及讀、寫方法和簡單應用》（http://www.shanxitv.org/125）。用文本編輯器打開這個cube文件，在第一行里寫上box2cell。這樣Multiwfn載入這個cube文件時，一旦發現第一行的內容包含box2cell，就會把格點數據的盒子信息轉化為晶胞信息，從而在計算的時候能夠考慮周期性（也可以不修改cube文件，而是載入cube文件后用隱藏的主功能1000里的子功能18來實現盒子信息向晶胞信息的轉化）。
把Multiwfn的settings.ini里的iuserfunc設為-3，使得user-defined function對應基于載入的格點數據通過B-spline插值的函數。然后啟動Multiwfn，載入Ag-ELECTRON_DENSITY-1_0.cube，之后輸入
2 //拓撲分析
-11 //修改被分析的函數
100 //User-defined function，當前對應于插值方式產生的價層電子密度
3 //將每兩個原子之間的中點作為搜索臨界點的初猜
然后進入選項0，并把原子標簽顯示出來、把鍵設細，就會看到下圖。可見臨近原子之間的BCP都已經找到了
下面再把更多的臨界點也找出來。當前情況不適合選擇選項4、5分別用3、4個原子中心作為初猜點來搜索臨界點。選項3用每兩個原子之間中點當做初猜點時是考慮循環相鄰晶胞中的鏡像原子的，而選項4、5則只會循環當前晶胞里那些唯一的原子，因此此例用選項4、5搜索后還是會遺漏許多重要的臨界點。當前最適合的方式是選擇選項6里的子選項-1，此時會在每個原子附近特定半徑內撒一批初猜點來搜索臨界點。這么選過之后，會發現找到了一大堆新的臨界點，進入選項0后看到的圖如下所示。下圖左側把所有臨界點都顯示了出來，可以看到離原子核較近的區域臨界點有不少，而且各種類型的都有。這是因為當前用了贗勢，由于電子密度不是從原子核開始向四周單調下降的，因此會在價層區域出現各種各樣的臨界點，這些不是我們感興趣的。把原子球設大掩蓋這些沒什么意義的臨界點，并且旋轉體系使得視角與晶軸平行后，就看到了下面的右圖。可見在原子間相互作用區域找到的臨界點的分布整齊、對稱，因此可以認為有意義的臨界點都找到了。
有一個概念叫平面性指數，在“量子化學波函數分析與Multiwfn程序培訓班”（http://www.keinsci.com/workshop/WFN_content.html）里有一頁幻燈片做了介紹，如下所示。我們這里來對Ag晶體計算一下這個平面性指數，看看是否和文章里說的情況吻合。注意這里說的平面性指數和《使用Multiwfn定量化和圖形化考察分子的平面性（planarity）》（http://www.shanxitv.org/618）里介紹的我提出的衡量體系幾何結構平面性的參數是兩碼事。
首先我們得確定要考察的籠臨界點也即(3,+3)的臨界點編號。下圖左邊只顯示了籠臨界點，可見當前體系有兩類籠臨界點，一種是諸如晶胞中心的那個（15號），它處于八面體中心。另一種是處于四面體中心的，比如63號臨界點。下圖右邊只顯示了所有BCP，此體系中所有的BCP都是等價的，我們考慮其中任意一個，比如1號。利用當前拓撲分析界面里的選項7依次對這三個臨界點的屬性進行考察，可以得到的這三個位置的基于格點數據插值產生的電子密度，也即User-defined real space function后面的值，結果如下（注意不要讀Multiwfn直接輸出的Density of electrons后面的值，那個是用孤立狀態原子的電子密度疊加得到的準分子密度）
CP(15): 0.1572533472E-01 a.u.
CP(63): 0.2212960297E-01 a.u.
CP(1): 0.2706526462E-01 a.u.
可見兩類籠臨界點里電子密度最小值是0.01572 a.u.，BCP處的電子密度為0.02706 a.u.，因此平面性指數為0.01572/0.02706=0.58，在0.5左右，符合幻燈片里說的“其它金屬及合金”的情況。
4 總結&其它
本文通過三個有代表性的例子詳細介紹了在Multiwfn中對周期性體系做電子密度的拓撲分析，也即一般說的AIM拓撲分析的完整過程，所有要注意的細節都做了解釋。通過例子可見，Multiwfn做周期性體系的拓撲分析的過程很簡單，并不比分析孤立體系更復雜。讀者如果實際操作一遍，還會感受到Multiwfn做這些分析的速度相當快。AIM拓撲分析對于考察周期性體系的電子結構有重要意義，推薦讀者在實際研究中運用，使文章增光添彩。
很值得一提的是，Multiwfn的拓撲分析功能極其普適，絕不僅限于電子密度的拓撲分析。對于Multiwfn自身可以直接計算的函數，以及載入的格點數據里記錄的函數，Multiwfn都可以做拓撲分析，對分子體系的例子參看比如《使用Multiwfn對靜電勢和范德華勢做拓撲分析精確得到極小點位置和數值》（http://www.shanxitv.org/645）對靜電勢、范德華勢的拓撲分析，以及《在Multiwfn中單獨考察pi電子結構特征》（http://www.shanxitv.org/432）里面對ELF-pi函數的拓撲分析。對周期性體系，Multiwfn也完全可以以相同的方式對任意函數做拓撲分析。
使用本文的做法做AIM拓撲分析在發表文章時必須按照Multiwfn啟動時的提示恰當引用Multiwfn，如果用于給別人代算也要明確告知對方這一點。