詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件

詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件

First release: 2022-Aug-31 Last update: 2022-Dec-17

0 前言

強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）已經有很多分析功能支持免費又速度極快的CP2K程序產生的周期性波函數，諸如AIM拓撲分析、原子電荷計算、Mayer鍵級計算、軌道成分分析、繪制DOS圖、模擬隧道顯微鏡圖、空穴-電子分析、IRI/IGMH/NCI相互作用分析、計算ELF/LOL/密度差、軌道定域化，等等等等，目前已明確支持周期性體系的功能詳見最新版Multiwfn手冊2.9.2.2和2.9.2.3節的列表。我也寫過一些CP2K結合Multiwfn做波函數分析的博文，參看http://www.shanxitv.org/category/CP2K/。

用Multiwfn做周期性體系波函數分析需要CP2K產生記錄原子、基函數、軌道等信息的molden文件作為輸入文件，在本文將全面、詳細介紹如何使用CP2K產生這樣的文件。本文是對于當前Multiwfn最新版本而言的，對于CP2K是對>=8.1版的情況而言的，并且截止到本文最后更新時的最新CP2K版本來說都是適用的。如果你不了解Multiwfn，建議參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167），另建議參看《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）以對Multiwfn的輸入文件支持情況有更全面的了解。

本文涉及到的文件可以在http://www.shanxitv.org/attach/651/file.zip里得到，便于讀者對照。

1 產生molden文件

在CP2K輸入文件的&DFT字段中加入以下內容就可以在SCF收斂后在當前目錄下產生.molden為后綴的文件

&PRINT
&MO_MOLDEN
NDIGITS 9
&END MO_MOLDEN
&END PRINT

其中NDIGITS控制輸出的軌道展開系數。上例要求絕對值大于1E-9的展開系數才會被輸出。NDIGITS越大不僅輸出的系數越多，系數保留的有效位數也越多，也因此在Multiwfn里分析結果越準確，但文件也相應地越大。NDIGITS的設置不影響在Multiwfn里的分析耗時。NDIGITS 9對于波函數分析的目的足夠精確了。

在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）里專門介紹過使用Multiwfn創建CP2K輸入文件的功能，在此界面里選-2將狀態切換為Yes后再產生的輸入文件即包含以上內容。

對任何CP2K支持的半經驗級別的方法，比如GFN1-xTB、PM6、DFTB，CP2K都無法產生molden文件。順帶一提，Grimme的xtb程序倒是可以產生GFN-xTB級別下的波函數文件，可以用于Multiwfn做波函數分析，不過xtb主要面向的是孤立而非周期性體系。

考慮k點時也無法產生molden文件，因此晶胞必須足夠大，不夠大就需要擴成超胞來算。

對于MD等涉及結構變化的任務，加了以上字段后，默認每一步都輸出一個molden文件，可以通過&MO_MOLDEN中的&EACH子字段控制對各種任務輸出molden文件的頻率。利用這個特征，可以令分子動力學模擬過程中不斷產生molden文件，然后寫腳本批量調用Multiwfn做波函數分析，以考察電荷、成鍵等方面隨模擬過程的變化，還可以制作成動畫。參考《詳談Multiwfn的命令行方式運行和批量運行的方法》（http://www.shanxitv.org/612）、《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）、《制作動畫分析電子結構特征》（http://www.shanxitv.org/86）。

2 盒子信息的添加

molden文件沒有標準的字段用來記錄盒子（或稱晶胞）信息，因此如果直接把CP2K對周期性體系產生的molden文件當Multiwfn的輸入文件用，Multiwfn做的分析將無法考慮周期性，顯然對于周期性體系來說得到的結果是明顯錯誤的，至少是對于靠近盒子邊界的原子來說。

為了讓Multiwfn獲得盒子信息，需要用文本編輯器打開molden文件，在里開頭部分插入[Cell]字段，里面通過三行來分別定義盒子的三個平移矢量的X、Y、Z分量（以埃為單位），例如

[Molden Format]
[Cell]
16.99927817     0.00000000     0.00000000
-8.49963911    14.72180676     0.00000000
0.00000000     0.00000000    22.46050431
[Atoms] AU
Cl       1      17       3.212398       1.854679      23.537502
Fe       2      26       0.000000       0.000000      21.222101
略...

顯然可以直接把CP2K輸入文件或者restart文件里的&CELL中的A/B/C后面的信息直接粘貼到[Cell]字段里。

也可以通過盒子的三個邊長(a,b,c)和夾角(α,β,γ)來定義盒子信息，例如下面的寫法定義a = 15 ?, b = 13 ?, c = 18.5 ?, α= 90°, β= 90°, γ= 121.3°
[Cell]
15 13 18.5 90 90 121.3

此外，如果你不想修改molden文件就能給Multiwfn提供盒子信息，可以在當前目錄下創建一個叫[Cell].txt的文件，內容是本應寫在[Cell]下面的信息，即對應三個平移矢量的三行信息，或者對應六個晶胞參數的一行信息。當Multiwfn發現molden文件里沒有[Cell]字段而當前目錄下有[Cell].txt，就會問你是否從此文件里讀取（此特性從2022-Dec-17更新的Multiwfn才支持）。

3 有效核電荷數的添加

molden文件不記錄原子的有效核電荷數（即元素序數減去被贗化的內核電子數，也即在計算中被基函數所描述的原子在孤立狀態下的價電子數），因此當使用贗勢時，若用molden文件作為Multiwfn的輸入文件，原子電荷等涉及到核電荷數的分析功能的結果可能不正常。另外，Multiwfn內置了一套內核電子密度數據庫EDFlib，只有當Multiwfn知道了原子被贗化了的電子是多少時，才能自動取用恰當的內核電子密度信息，這帶來的好處是直接基于電子密度的分析結果將和使用全電子基組時基本一樣，比如AIM拓撲分析可以找到核臨界點、能產生出完整的鍵徑。為了令Multiwfn在載入molden文件時能正確地得知有效核電荷數，需要用文本編輯器打開molden文件并手動修改。

有兩種改法，第一種方法是把molden文件的[Atoms]字段里每個原子的第三列改為其有效核電荷數，比如
[Molden Format]
[Atoms] AU
Cl       1       7       3.212398       1.854679      23.537502
Fe       2      16       0.000000       0.000000      21.222101
Cl       3       7       0.000000       3.709358      18.906700
略...

另一種辦法是在molden文件頭部插入一個[Nval]字段，里面記錄各種元素的有效核電荷數，比如
[Molden Format]
[Nval]
Cl 7
Fe 16
[Atoms] AU
Cl 1 17 3.212398 1.854679 23.537502
Fe 2 26 0.000000 0.000000 21.222101
略...

以上兩種寫法都代表Cl的有效核電荷數是7，Fe是16。

在當前用的贗勢下，各種元素的有效核電荷數是多少可以去看基組定義，比如在CP2K的data目錄下的BASIS_MOLOPT文件中可以看到有一行Cl DZVP-MOLOPT-SR-GTH DZVP-MOLOPT-SR-GTH-q7，其中q后面的數字就是有效核電荷數（價電子數）。在目前版本Multiwfn產生的使用GTH贗勢的輸入文件中，直接看BASIS_SET后面的基組名的q后面的值立馬就知道有效核電荷數是多少。

有時候計算是在遠程服務器上，如果由于特殊原因編輯molden文件不方便的話，可以在本機寫一個文本文件，比如add.txt，里面是[Cell]和[Nval]字段，比如
[Cell]
22.55600000     0.00000000     0.00000000
-11.27800000    19.53406901     0.00000000
0.00000000     0.00000000     6.80000000
[Nval]
O 6
N 5
C 4
<---有空行

之后把add.txt上傳，然后運行cat add.txt org.molden > final.molden，就把[Cell]和[Nval]字段加入到原先的org.molden文件開頭，得到了可以用于Multiwfn做分析的final.molden。

4 關于軌道的記錄

CP2K默認是不計算空軌道的，因此默認產生的molden文件里只有占據軌道信息，這對于只依賴于占據軌道信息的分析足矣，諸如原子電荷計算、Mayer鍵級計算、IRI/IGMH分析，但沒法用于涉及到空軌道的分析，如電子激發分析、看空軌道圖形、繪制DOS（涉及費米能級以上區域時）等。

對于使用對角化做SCF的情況，在&SCF里寫ADDED_MOS N就會把最低N個空軌道也計算并輸出到molden文件里。而OT則不兼容ADDED_MOS關鍵詞。

注意Multiwfn內部要求基函數數目與軌道數目是相等的，比如你的molden文件定義了100個基函數，但只記錄了20個軌道的信息，那么載入Multiwfn后會有100個軌道，其中21~100號軌道是空的（軌道系數、占據數、能量都為0）。

用對角化時，軌道的能量會在SCF過程中計算并最終寫入到molden文件里，然而OT方式的計算則無法直接向molden文件里寫入軌道能量。molden文件里沒有能量信息的話自然沒法做牽扯到能量的分析，比如沒法繪制DOS、STM圖。用OT算大體系更快，對復雜大體系有時候不用OT還極難收斂。對于必須用OT又需要把軌道能量寫入molden文件的情況，可以先做一次OT計算，完成之后在當前目錄下會出現記錄收斂的波函數的wfn文件（注意，這和《高斯fch文件與wfn波函數文件的介紹及轉換方法》http://www.shanxitv.org/55里介紹的能被Multiwfn直接載入的那種wfn文件完全是兩碼事）。然后再用相同的設定做一次對角化的單點計算，并且要求產生molden文件、要求從wfn中讀取之前收斂的波函數當初猜，有需要的話還可以同時再加上ADDED_MOS要求計算出空軌道。此時只需要一輪或很少幾輪SCF迭代就可以完成，耗時很低，得到的molden文件里就有軌道能量信息了。舉個例子，下面對Multiwfn文件包里examples\COF_12000N2.cif體系先用OT計算，然后做對角化計算以得到記錄最低20個空軌道和軌道能量信息的molden文件，用的是CP2K 2022.1。

啟動Multiwfn，然后輸入
examples\COF_12000N2.cif
cp2k //進入產生CP2K輸入文件的界面
lycoris.inp //要產生的CP2K輸入文件名
4 //從默認的對角化切換為OT
0 //產生CP2K輸入文件
cp2k //再次進入產生CP2K輸入文件的界面
lycoris2.inp //輸出的文件名
4 //從OT切換為對角化
-2 //要求產生molden文件
-9 //其它設置
12 //設置要求的空軌道數
20 //算20個
0 //返回
0 //產生CP2K輸入文件
當前目錄下的lycoris.inp是基于OT算單點的任務的輸入文件，算完后當前目錄下就有了lycoris-RESTART.wfn。然后把lycoris2.inp里的
# WFN_RESTART_FILE_NAME lycoris2-RESTART.wfn
改為
WFN_RESTART_FILE_NAME lycoris-RESTART.wfn
并把SCF_GUESS RESTART這行開頭的#去掉以啟用之。然后運行lycoris2.inp，算完后當前目錄下就有了需要的lycoris2-MOS-1_0.molden。

如果你想檢查molden文件里有沒有軌道能量信息，把molden文件載入Multiwfn后，進入主功能0，左上角選擇Orbital info. - Show occupied orbitals，從文本窗口里顯示的能量上就能判斷，如果molden文件沒有能量信息的話所有軌道能量都為0。如果當前不方便進入圖形界面的話，也可以進入主功能6的子功能3。

5 關于做波函數分析時基組的選用

CP2K用戶通常使用CP2K開發者搞的MOLOPT系列贗勢基組做計算，此系列基組用于波函數分析完全可以，但如果你是做實空間函數的分析（對三維空間中一批點做特定函數的計算），比如AIM拓撲分析、IGMH/IRI/ELF分析等，而且體系又很大、CPU又不給力，用MOLOPT系列的話在Multiwfn中的計算耗時會很高，主要在于MOLOPT系列基組是完全廣義收縮基組，源高斯函數特別多。因此希望耗時更低的話，可以用諸如TZVP-GTH或更好的TZV2P-GTH贗勢基組，其質量做波函數分析就算不錯了。如果想用全電子基組，pob-TZVP是很好的選擇，大小適中，精度也不錯。

Multiwfn支持的分析當中有一部分在原理上是不兼容彌散函數的，如果基函數的彌散特征比較明顯，分析結果會很糟糕，比如Mulliken分析、SCPA分析、Pipek-Mezey軌道定域化等。我實測發現這些分析往往不適合用諸如DZVP-GTH、TZVP-GTH這樣的基組，因為有彌散程度較高的基函數。對這些分析可以用pob-TZVP，或者MOLOPT-SR-GTH系列基組。雖然MOLOPT-SR基組涉及的一些源高斯函數也有彌散特征（指數很小），之所以它此時能用是在于它是完全廣義收縮的，并不獨立存在彌散程度明顯的基函數。

6 總結

本文詳細介紹了如何用CP2K產生能給Multiwfn做波函數分析用的molden文件。最關鍵的就是這幾點：
(1)別忘了把盒子信息添加進去
(2)用贗勢時別忘了把有效核電荷數添加進去
(3)需要空軌道時別忘了讓CP2K求解空軌道
(4)別忘了用OT時molden文件里是沒有軌道能量信息的
(5)注意基組的選擇