使用Multiwfn結合CP2K對周期性體系做電荷分解分析（CDA）
使用Multiwfn結合CP2K對周期性體系做電荷分解分析（CDA）
Using Multiwfn in combination with CP2K to perform charge decomposition analysis (CDA) for periodic systems
文/Sobereva@北京科音 2024-Jul-1
1 前言
筆者之前寫過《使用Multiwfn做電荷分解分析(CDA)、繪制軌道相互作用圖》（http://www.shanxitv.org/166）介紹了怎么用Multiwfn程序基于Gaussian等量子化學程序產生的波函數文件對分子體系做電荷分解分析（CDA）。這種方法可以從片段軌道間相互作用的角度深入了解片段間是如何轉移電子的。由于CDA的很高價值，Multiwfn做CDA分析已經得到非常廣泛的使用。如果讀者沒看過此http://www.shanxitv.org/166的話一定要先仔細看一下，否則無法理解下文的內容。
從2024-Jun-5更新的Multiwfn開始，CDA模塊已經支持了基于CP2K產生的molden文件做周期性體系的CDA分析，下面將通過一個簡單的例子進行演示。Multiwfn可以在主頁http://www.shanxitv.org/multiwfn免費下載，不了解此程序者參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。我假定讀者已經具備了CP2K的基本使用知識，不了解者推薦通過北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）完整系統地學習。
本文的例子是計算NaCl板吸附CO分子的體系，用CDA方法考察CO與NaCl板之間電子轉移情況。這個體系優化后的結構如下，C和O分別是1和2號原子，3-110號原子是NaCl板。對這個體系，之前我在《使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線》（http://www.shanxitv.org/638）中通過密度差圖的方式定性考察了電子轉移，在《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）里通過片段電荷的方式定量考察了電子轉移。而下面用CDA分析，可以以軌道角度在明顯更深層次理解電子轉移細節。
本文例子涉及的所有文件都在http://www.shanxitv.org/attach/716/file.7z里。CP2K用的版本是2024.1，Multiwfn是2024-Jun-5更新的版本，讀者絕對不要用Multiwfn的更老版本。
2 產生輸入文件
首先我們要用CP2K計算NaCl+CO、NaCl、CO各自的molden格式的波函數文件。如果你不知道怎么用CP2K產生molden文件的話，先仔細閱讀《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651），我假定讀者已經看過本文了。這里有幾個要點：
(1)整體和片段的molden文件里的原子坐標必須對應
(2)整體和片段的計算必須使用嚴格相同的級別和設置，用的盒子也應當一致
(3)必須要求CP2K計算所有空軌道
(4)整體的molden文件里需要加入[Cell]字段定義晶胞信息，而片段的molden文件里加不加入不影響CDA結果。[Nval]字段可加入可不加入，也不影響CDA結果
本文文件包里的NaCl_CO.cif是以前我在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH下優化的NaCl板吸附CO的結構文件。用Multiwfn載入它，然后輸入
cp2k //創建CP2K輸入文件
NaCl_CO.inp //產生的輸入文件名
-2 //要求產生molden文件
-9 //其它設置
12 //設置計算的空軌道數
-1 //計算所有空軌道
0 //返回
0 //產生輸入文件
用CP2K運行NaCl_CO.inp，得到NaCl_CO-MOS-1_0.molden。然后將以下字段插入到molden文件的開頭
[Cell]
16.92168000 0.00000000 0.00000000
0.00000000 16.92168000 0.00000000
0.00000000 0.00000000 25.00000000
將NaCl_CO.inp復制為NaCl.inp，再將NaCl.inp中&COORD里的CO部分刪掉、PROJECT NaCl_CO改名為PROJECT NaCl，然后用CP2K運行之，得到NaCl-MOS-1_0.molden。將上面的[Cell]加入其中（這和CDA分析無關，加入這個的目的是為了之后在Multiwfn中能正確觀看NaCl板的周期性的晶體軌道）。
將NaCl_CO.inp復制為CO.inp，再將CO.inp中&COORD里的NaCl部分刪掉、PROJECT NaCl_CO改名為PROJECT CO，然后用CP2K運行之，得到CO-MOS-1_0.molden。這個可以不用刻意加上[Cell]，因為當前體系中CO離盒子邊界很遠，顯示軌道時是否考慮周期性對看到的圖像沒影響。
接下來就可以進行CDA分析了。
3 進行CDA分析
啟動Multiwfn，載入NaCl_CO-MOS-1_0.molden，然后依次輸入
16 //CDA分析
2 //兩個片段
CO-MOS-1_0.molden //CO的波函數文件。注意要按整個體系里原子序號順序載入片段
NaCl-MOS-1_0.molden //NaCl板的波函數文件
現在屏幕上輸出了一大堆復合物軌道的CDA結果，同時還輸出了ECDA結果：
The net electrons obtained by frag. 2 = CT( 1-> 2) - CT( 2-> 1) = 0.1580
這個ECDA結果告訴你CO（片段1）往NaCl板（片段2）轉移了0.158個電子，這和《使用Multiwfn對周期性體系計算Hirshfeld(-I)、CM5和MBIS原子電荷》（http://www.shanxitv.org/712）里通過Mulliken方法計算的CO的片段電荷嚴格相同。
當前CDA分析輸出的復合物軌道太多，一共437個（對應整個體系的占據軌道數），而絕大多數對d、b、r項的貢獻微乎其微或完全為0，沒有考察的意義。為了簡化輸出，在當前界面里輸入
-3 //設置CDA輸出時d、b、r項的閾值，它們之一的絕對值大于這個閾值的復合物軌道才會輸出
0.005
0 //顯示CDA和ECDA結果
當前結果如下
Orb. Occ. d b d - b r
272 2.000000 0.006220 -0.000012 0.006232 0.000488
275 2.000000 0.125862 -0.001874 0.127736 0.051155
297 2.000000 0.001098 0.000305 0.000792 -0.005373
318 2.000000 0.001053 0.000326 0.000727 -0.005720
341 2.000000 0.000781 0.000200 0.000581 -0.006708
352 2.000000 0.002291 0.000731 0.001560 -0.015705
360 2.000000 0.001101 0.000298 0.000803 -0.009297
414 2.000000 0.000635 0.000214 0.000421 -0.005229
427 2.000000 0.000739 0.000219 0.000520 -0.005731
-------------------------------------------------------------------
Sum: 874.000000 0.153805 0.017732 0.136074 -0.037831
Note: The "Sum" includes all terms including those not printed above
CDA的總結果是CO向NaCl轉移了d=0.154個電子，NaCl向CO反饋了b=0.018個電子，NaCl凈獲得d-b=0.136個電子。b項基本可以忽略不計，而d的主要貢獻是由275號復合物軌道造成的。這個復合物軌道的貢獻顯然值得進一步探究，看看是由CO和NaCl板的哪些片段軌道混合所造成的。為此，輸入
5 //將特定復合物軌道對CDA項的貢獻進行分解
275 //要考察的復合物軌道序號
0.005 //輸出閾值
現在看到如下結果
Occupation number of orbital 275 of the complex: 2.00000000
FragA Orb(Occ.) FragB Orb(Occ.) d b d - b r
5( 2.0000) 271( 2.0000) 0.000000 0.000000 0.000000 0.005183
5( 2.0000) 289( 2.0000) 0.000000 0.000000 0.000000 0.005642
5( 2.0000) 336( 2.0000) 0.000000 0.000000 0.000000 0.005127
5( 2.0000) 341( 2.0000) 0.000000 0.000000 0.000000 0.008960
5( 2.0000) 435( 0.0000) 0.009358 0.000000 0.009358 0.000000
5( 2.0000) 439( 0.0000) 0.011934 0.000000 0.011934 0.000000
5( 2.0000) 443( 0.0000) 0.010929 0.000000 0.010929 0.000000
5( 2.0000) 447( 0.0000) 0.008032 0.000000 0.008032 0.000000
5( 2.0000) 453( 0.0000) 0.009177 0.000000 0.009177 0.000000
5( 2.0000) 465( 0.0000) 0.010755 0.000000 0.010755 0.000000
5( 2.0000) 469( 0.0000) 0.013735 0.000000 0.013735 0.000000
5( 2.0000) 473( 0.0000) 0.011282 0.000000 0.011282 0.000000
Sum of above terms: 0.085201 0.000000 0.085201 0.024911
可以看到CO向NaCl轉移電子靠的都是CO的5號占據軌道和NaCl空軌道相互作用。NaCl板接收CO轉移來的電子分散在其大量空軌道上，比如NaCl板的443號空軌道接受了0.011個電子。在當前的輸出閾值下，上面輸出的這些項的所有的d項加和為0.085，明顯和275號復合物軌道的總d值0.153有很大差距，說明CO還向上面沒輸出的很多其它的NaCl板的空軌道轉移了電子。
下圖把275號復合物軌道、CO的5號占據軌道，以及NaCl接收CO貢獻來的電子最多的兩個軌道（439、469）展示在了一起，便于大家直觀了解情況，等值面數值用的是0.02。復合物和片段的軌道可以用Multiwfn載入相應的波函數文件后用主功能0按照《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）說的觀看。從此圖可以清晰直觀地認識到CO的孤對電子往吸附CO的那個Na位點的空軌道上轉移了電子。還有大量也具有這種特征的NaCl板的空軌道也接收了CO的孤對電子，這里沒全都畫出來。
下面再看一下上圖中275號復合物軌道是怎么由片段軌道構成的。按0返回到CDA菜單，然后進入選項2，輸入275，看到以下信息
Note: Only the fragment orbitals with contribution > 1.0 % will be shown below, the threshold can be changed by "compthresCDA" in settings.ini
Occupation number of orbital 275 of the complex: 2.00000000
Orbital 5 of fragment 1, Occ: 2.00000 Contribution: 87.42 %
Sum of values shown above: 87.42 %
這告訴你275號復合物軌道由CO的5號占據軌道貢獻了87.4%。默認情況下只有貢獻大于1%的片段軌道才會輸出，雖然NaCl板的一大堆非占據軌道都通過與CO的5號軌道混合而參與了這個復合物軌道的構成，但貢獻都小于輸出閾值，所以上面沒看到。從當前信息可以認識到275號復合物軌道主要體現CO的5號軌道特征，這是為什么上面圖中275號復合物軌道和CO的5號軌道看起來十分相似。
4 總結
此文通過一個簡單的表面吸附的例子演示了怎么用Multiwfn的CDA功能分析周期性體系的電荷轉移本質。討論電子轉移常見的套路就是繪制電子密度差、算算片段電荷，而從本文的例子看到Multiwfn可以從片段軌道相互作用對電子轉移的本質進行深入的考察，掌握了本文介紹的方法明顯可以給分析增光添彩。
使用Multiwfn做本文介紹的分析請記得在發表文章時引用Multiwfn程序啟動時提示的Multiwfn原文，以及引用介紹CDA/GCDA方法的文章（Multiwfn用的實質上是我對CDA廣義化后的GCDA形式），這在進入Multiwfn的CDA功能時屏幕上明確提示了。