思想家公社的門口：量子化學·分子模擬·二次元 - CP2K

CP2K做雜化泛函計算的關鍵要點和簡單例子

sobereva — Tue, 07 Nov 2023 17:53:00 +0800

CP2K做雜化泛函計算的關鍵要點和簡單例子

First release: 2023-Nov-7 Last update: 2024-Feb-18

0 前言

對周期性體系，雜化泛函的計算耗時眾所周知遠高于純泛函。著名的第一性原理程序CP2K不光純泛函的計算超級快，對大體系做雜化泛函計算也非常快（盡管還是遠比純泛函貴），遠遠快于Quantum ESPRESSO (QE)、VASP等純粹基于平面波的程序。CP2K結合像樣的基組做雜化泛函的單點計算，在一般雙路服務器上算到上千原子都沒問題。然而，使用CP2K做高效的雜化泛函計算有很多非常關鍵性的要點必須知道，遠遠不是像Gaussian程序那樣光寫個雜化泛函名字就了事的。然而，筆者經常在網上回答CP2K問題時看到很多人完全不具備這些最基本常識，胡用瞎用CP2K做雜化泛函計算，導致耗時巨高、始終出不來結果、任務崩潰、SCF完全不收斂等各種問題。因此筆者覺得很有必要寫個小文將CP2K做雜化泛函計算的最關鍵常識和要點挨個提一下。Multiwfn程序（http://www.shanxitv.org/multiwfn）具有極其便利的創建CP2K輸入文件的功能，在本文的最后會還用Multiwfn產生C60晶體的CP2K雜化泛函計算的輸入文件作為實際例子。

本文的內容至少對CP2K 2024.1，以及2023-Nov-7及以后更新的Multiwfn是適用的。

筆者在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）的“能量的計算及相關問題”部分里講雜化泛函計算遠比本文全面、深入得多，并列舉諸多經驗技巧和例子，還詳細介紹HSE06、CAM-B3LYP、wB97X等各種范圍分離形式的雜化泛函以及revDSD-PBEP86-D3(BJ)等各種雙雜化泛函的用法和要點，使用雜化泛函結合k點計算的要點，歡迎想學習更多知識者關注和參加。

1 CP2K雜化泛函計算的特性

---------
2024-Feb-18更新：以下關于k點的說法已經過時。從2024.1版開始支持了雜化泛函計算時考慮k點，并且可以結合ADMM。此時需要用RI計算HF交換部分。

首先要知道CP2K做雜化泛函計算的一個關鍵不足是不支持考慮k點，因此如果被計算的是周期性體系，晶胞必須大到只考慮gamma點就足夠。如果你要計算的體系的晶胞原本比較小，顯然必須先擴胞到足夠大才行。如果不知道擴到多大夠用，嚴格來說需要做你感興趣的屬性的計算結果隨擴胞倍數的收斂性測試。雖然CP2K能夠對含有非常多原子的超胞做得動雜化泛函計算，但耗時會明顯高于QE等程序對原胞在考慮k點的情況下做雜化泛函計算。

最能充分體現CP2K雜化泛函計算效率優勢的情況是算一些較大分子晶體、MOF、COF等一些大晶胞體系，以及做第一性原理動力學等情況。對小晶胞體系如上所述經過擴胞后也能算，只不過不體現CP2K的長處。

CP2K顯然沒法用雜化泛函算能帶，因為不支持k點，這種目的只能用QE等程序代替。必須用CP2K算的話，應當用CP2K支持的純泛函里算帶隙相對比較好的比如HLE17。
---------

由于雜化泛函遠比普通泛函昂貴，尤其是用于幾何優化、NEB、動力學等涉及到結構變化的任務。因此如果你算的問題用純泛函就足以描述得不錯、用雜化泛函不會帶來明顯好處的話，顯然應該用純泛函。記住，非必要甭考慮雜化泛函，純泛函的耗時僅有雜化泛函的一個微不足道的零頭。對幾何優化、振動分析、NEB等高耗時任務用純泛函，算勢壘、反應能的時候用更好的雜化泛函是可以接受且常見的策略。

周期性體系的雜化泛函計算不僅耗時遠高于純泛函計算，對內存要求更是遠遠高于純泛函計算。經常用CP2K做雜化泛函任務的服務器的內存容量絕對是多多益善，選擇服務器的時候要注意這一點。也別天真地試圖拿個人計算機將就著跑CP2K的雜化泛函計算。

CP2K做雜化泛函計算的前提是編譯時必須帶著libint電子積分庫，否則沒法算雜化泛函中HF交換項中涉及的雙電子積分。CP2K的編譯方法見http://www.shanxitv.org/586。

雜化泛函用的贗勢基組和贗勢就用常見的給PBE等純泛函優化的就可以。從CP2K 9.1開始，data目錄下的BASIS_MOLOPT_UZH和POTENTIAL_UZH文件里還分別給出了專門給PBE0優化的贗勢基組和贗勢，如C的DZVP-MOLOPT-PBE0-GTH-q4和GTH-PBE0-q4，原理上來說它們用于雜化泛函計算會更理想一些（但實際未必有什么優勢）。

2 庫侖截斷

雜化泛函計算之所以比純泛函昂貴得多，在于計算雜化泛函的HF交換項時要計算大量的雙電子積分（two-electron integral, ERI），尤其是對于周期性體系來說，不光要計算晶胞內基函數之間的ERI，還要計算中心晶胞與周圍鏡像晶胞基函數之間的ERI。再加上由于只能算gamma點而晶胞必然較大，原子數、基函數通常非常多，無疑要算的ERI是巨量的。

由于周期邊界條件，原理上來說要計算的ERI是無窮的，顯然這會導致計算無法進行。因此實際計算時必須使用庫侖截斷（Coulomb truncation），也就是將ERI中描述電子庫侖相互作用的1/r算符截斷在某個距離（截斷半徑），超過截斷半徑時1/r視為0，相應地ERI就不需要計算了，這個策略使得ERI的計算數目是有限的。如果你計算的是非周期性體系，由于本身ERI數目就是有限的，因此雖然也可以用庫侖截斷，但這就不是必須的了。

雜化泛函的庫侖截斷的半徑設置顯然既影響耗時也影響結果。截斷半徑越小，要計算的ERI越少，耗時就越低，但結果越偏離泛函原本的精確結果。顯然不能為了一味地節約時間而把截斷半徑設得太小，否則結果會很垃圾。6埃通常是耗時和精度的較好權衡，設得明顯更小會令精度有明顯損失，設得更大并不會令精度有明顯提升卻會增加耗時不少。如果對精度要求很高也可以用8埃，設得更大就完全沒必要了。

要注意像6、8埃這種程度的截斷半徑其實還沒有大到令絕對能量充分收斂，但能量在相互求差時，由截斷半徑帶來的系統性誤差能很大程度抵消。可以認為用不同的截斷半徑相當于用不同的理論方法，直接對比能量的情況必須保持截斷半徑的設置相統一！

值得一提的是，不一定截斷半徑小結果就一定差，有可能截斷半徑造成的對ERI忽略的誤差和泛函本身的誤差相抵消導致恰好對某些體系某些問題的計算比用更大截斷半徑時結果還更好。

截斷半徑應當小于晶胞最短尺寸的一半以免導致得到缺乏物理意義的結果。這里所謂的最短尺寸，是(001)晶面間距、(010)晶面間距、(100)晶面間距三者中的那個最小值。也因此為了能夠使用6埃截斷半徑，晶胞最短尺寸應當大于12埃，如果晶胞不夠大則應當擴胞（即便拋開庫侖截斷問題不談，單從只能計算gamma點這點來說，晶胞也不能過小）。

Multiwfn創建的雜化泛函輸入文件中，若周期性設成了NONE，則不使用庫侖截斷，若是周期性，Multiwfn會自動加上庫侖截斷的設置。如果晶胞最短尺寸夠大，Multiwfn會自動把截斷半徑設為6埃，如果不夠大，則會根據當前晶胞最短尺寸自動設成能設的最大半徑。

以下是庫侖截斷的設置方式，加入到&DFT/&XC里面。POTENTIAL_TYPE TRUNCATED代表用庫侖截斷的1/r算符，CUTOFF_RADIUS設置截斷半徑（埃）。T_C_G_DATA t_c_g.dat是默認的可以不寫，這指定的是庫侖截斷用的gamma函數的數據文件t_c_g.dat，在CP2K的data目錄下。
&HF
&INTERACTION_POTENTIAL
    POTENTIAL_TYPE TRUNCATED
    CUTOFF_RADIUS 6.0
    T_C_G_DATA t_c_g.dat
&END INTERACTION_POTENTIAL
&END HF

3 積分屏蔽

積分屏蔽（integral screening）是指盡可能在不明顯犧牲精度的前提下減少要計算的ERI的量，這對于CP2K雜化泛函計算的耗時和內存使用量有至關重要的影響！上一節的庫侖截斷本身就是一種積分屏蔽策略，CP2K里還有另外兩種非常重要的積分屏蔽方法，也是普遍應用在Gaussian等量子化學程序中的：
(1)利用Schwarz不等式做積分屏蔽：利用Schwarz不等式關系，可以根據數目較少的二指標ERI估計出數目甚巨的三、四指標ERI中哪些數值肯定小于閾值從而可以忽略掉，這可以大幅減少ERI要算的數目（原本總共要算的ERI中占比最大的就是四個基函數序號不同的那些ERI）。閾值通過&DFT/&XC/&HF/&SCREENING中EPS_SCHWARZ設置，此值越大，被忽略的ERI數目越多、計算精度越差。默認的1E-10雖然足夠保證精度，但此時太貴，沒必要。通常設為1E-6就可以了，是精度和耗時的較好權衡，想要精度更好點用1E-8。
(2)結合密度矩陣做積分屏蔽：ERI對能量的貢獻取決于它與特定密度矩陣元的乘積。因此哪怕ERI不是很小，但與之相乘的密度矩陣元很小，那么這樣的ERI也完全不需要計算。CP2K定義了特定方法快速估計一個ERI與相應密度矩陣元乘積值的上限，如果小于上述的EPS_SCHWARZ，則此ERI就會忽略掉。這種依賴于密度矩陣的積分屏蔽默認是關閉的，如果開啟則把&DFT/&XC/&HF/&SCREENING里的SCREEN_ON_INITIAL_P設為T。由于這種積分屏蔽能令耗時巨幅降低，因此我強烈建議總是開啟，也因此Multiwfn產生的雜化泛函的輸入文件里都自動把SCREEN_ON_INITIAL_P設成T。

基于密度矩陣做積分屏蔽有一個極其關鍵的要點是初猜的密度矩陣質量不能太差！然而很多人對此一無所知！默認情況下初始的密度矩陣是CP2K自動做初猜產生的，由于它和SCF收斂時的密度矩陣相差很大、質量很糙，因此此時積分屏蔽通常做得極爛，不僅可以忽略的ERI可能沒忽略，還會把一些重要的ERI給忽略掉，由此可能導致KS矩陣質量很爛、SCF極難收斂。而且忽略哪些ERI是在SCF過程一開始決定好的，即便之后隨著SCF迭代密度矩陣質量不斷變好，也不會由此使得積分屏蔽逐步變得合理，這導致即便最后SCF收斂了，得到的能量的精度也可能很爛。因此當SCREEN_ON_INITIAL_P設T時，必須讀取一個基本合理的波函數當初猜。通常是在雜化泛函計算前先用相同的基組做一個便宜的純泛函計算，比如PBE0計算前先用PBE做一個單點計算，得到記錄了PBE收斂的波函數信息的wfn文件，然后PBE0計算時用SCF_GUESS RESTART，并把WFN_RESTART_FILE_NAME指定成那個wfn文件，以讀取其作為初猜。這不僅使得基于密度矩陣的積分屏蔽可以合理進行，還有一個額外的好處是由于初猜波函數質量較好（至少比默認自動產生的好得多），使得雜化泛函計算時達到SCF收斂所需的迭代次數會比較少，這進一步節約了時間（注意即便拋開ERI計算耗時不談，雜化泛函在構造KS矩陣時花的時間也明顯高于純泛函，因為涉及到密度矩陣元與ERI的大量的乘加操作）。

一般的量子化學程序雖然也基于密度矩陣做積分屏蔽，但不是非得像CP2K這樣得讀取一個像樣的初猜波函數，這是因為它們會根據當前最新的密度矩陣判斷這一輪SCF用的ERI哪些是能被忽略的。CP2K的做法和它們不同，主要是因為CP2K傾向于讓用戶做incore形式的SCF而非一般量化程序一般用的direct形式的SCF，見后文。

當SCREEN_ON_INITIAL_P設T時，對于雜化泛函的幾何優化、分子動力學等任務的續算，要注意光有.restart文件還不行，還必須提供上一步的.wfn文件當初猜波函數，否則也是會由于基于自動初猜的密度矩陣做積分屏蔽非常惡心而令任務無法正常跑下去。

如果你嫌麻煩而不想先做純泛函計算產生收斂的波函數再做雜化泛函計算，那就不要寫SCREEN_ON_INITIAL_P T。對于小體系、小基組特別是非周期性的計算，不基于密度矩陣做積分屏蔽時耗時也不高。而對于計算量較大的雜化泛函計算，則強烈建議不要偷這個懶而多花巨量不必要的時間。

4 ADMM

CP2K還支持其開發者獨創的ADMM (Auxiliary Density Matrix Methods)方法減少要算的ERI數目從而進一步顯著加快雜化泛函的計算，其原理見以下北京科音CP2K第一性原理計算培訓班的ppt。

ppt里提到的主基組就是指平時實際用的基組，諸如DZVP-MOLOPT-SR-GTH。對于不同的主基組，有一些專門構造的與之匹配的ADMM輔助基組，比如MOLOPT系列基組適合搭配的是CP2K的data目錄下的BASIS_ADMM_UZH文件里的那些輔助基組，從小到大依次是admm-dz、admm-dzp、admm-tzp、admm-tz2p。admm-dzp是2-zeta結合極化函數的輔助基組，通常算是保底質量。DZVP-MOLOPT-SR-GTH主基組適合搭配admm-dzp，更好的主基組如TZVP-MOLOPT-GTH結合admm-dzp也可以接受，若結合admm-tzp精度會更好些，這都比起不用ADMM能巨幅節約計算時間和內存。輔助基組的大小與主基組越接近，ADMM給雜化泛函HF交換項帶來的誤差越小，但需要計算越多ERI因而越耗時。如果輔助基組和主基組設為相同的，那么ADMM不會造成任何誤差，相應地也不會節約任何耗時（反倒還會增加一些耗時）。

data目錄下的BASIS_ADMM文件中的FIT3、pFIT3等輔助基組都比較過時了，對元素涵蓋也遠不如BASIS_ADMM_UZH里的完整，因此如今不再建議使用（它們出現在不少官方例子文件和文獻中，如今不要效仿）。目前版本的Multiwfn產生利用ADMM的雜化泛函的輸入文件時默認設的是admm-dzp。

注意前述這些輔助基組都是對GTH贗勢基組用的，千萬別做全電子計算時也用它們。要想全電子計算時用ADMM，可以用(aug-)pcseg系列基組作為主基組，在《在線基組和贗勢數據庫一覽》（http://www.shanxitv.org/309）中介紹的BSE基組數據庫中有和它匹配的(aug-)admm輔助基組，諸如pcseg-2適合搭配admm-2。pcseg是一類構造得很好的對DFT計算非常劃算的基組，《談談量子化學中基組的選擇》（http://www.shanxitv.org/336）里簡要提及了，在《北京科音基礎量子化學培訓班》（http://www.keinsci.com/workshop/KBQC_content.html）里講基組的地方做了具體介紹。

用MOLOPT系列基組做周期性體系的雜化泛函計算幾乎是一定要結合ADMM的，否則極難算得動！因為MOLOPT是完全廣義收縮基組，而CP2K利用的計算源高斯函數（PGTF）之間ERI的libint庫又沒有專門考慮廣義收縮的情況，因此造成等效的PGTF數目甚巨，ERI的數目又與PGTF數目的四次方近似成正比，數目顯然多得恐怖。如果你用片段收縮基組，比如def2-SVP、6-31G**、pcseg-1、pob-DZVP-rev2之類去做雜化泛函計算，由于需要考慮的PGTF數目遠比檔次與之相仿佛的DZVP-MOLOPT-SR-GTH要少，因此ADMM方法不是非用不可，但不用的話耗時會顯著高于DZVP-MOLOPT-SR-GTH搭配admm-dzp的情況。

不能有的原子用ADMM加速計算而有的原子不用。因此如果你用混合基組，其中有的原子用的主基組有ADMM輔助基組，有的沒有，你又想用ADMM，那么沒有ADMM輔助基組的那些原子的輔助基組就只能設定成其主基組。雖然這些原子沒法被ADMM技術加速，但至少其它那部分原子能被加速。

前面提到的各種節約耗時的策略是彼此完全兼容的，同時利用可以最大程度節約時間。即周期性體系雜化泛函計算時一般要用庫侖截斷、利用Schwarz不等式做積分屏蔽、結合密度矩陣做積分屏蔽、開啟ADMM。所有這些設定都會對能量產生影響，因此寫文章的時候最好注明，橫向對比時應統一。

啟用ADMM需要在&KIND里加上BASIS_SET AUX_FIT [輔助基組名]，同時增加一行BASIS_SET_FILE_NAME [輔助基組文件名]。并且在&DFT里加入以下內容
&AUXILIARY_DENSITY_MATRIX_METHOD
&END AUXILIARY_DENSITY_MATRIX_METHOD
如果用的是對角化而非OT，還需要在以上字段里加入ADMM_PURIFICATION_METHOD NONE要求不做純化。

5 ERI的儲存和內存使用量的控制

SCF有兩種常見形式，incore SCF是指SCF計算一開始將所有要用到的ERI全都算出來并儲存在內存里，之后SCF每一輪在構造KS矩陣時會直接從內存里讀取ERI而不再重新計算。direct SCF是指SCF每一輪都重新計算要用到的ERI，這也叫on-the-fly方式計算ERI。顯然，incore比direct總耗時低得多，但代價是機子的必須內存非常大，因為要存的ERI通常非常多。在CP2K的雜化泛函計算中，如果你給的內存足夠存得下所有ERI，那么只有SCF第一輪的時候由于要計算所有ERI所以耗時很高，而之后SCF每一輪的耗時就很低了。如果給的內存不夠儲存所有ERI，那么能在內存里存多少ERI就存多少，存的那部分在之后每一輪SCF中會直接讀取，存不下的ERI在每一輪SCF中會重新算，相當于介于incore和direct之間的情況。顯然，做大體系雜化泛函計算時應當盡量在大內存機子上做、分配給CP2K盡可能多的內存，以讓盡可能多的ERI能存到內存中、盡可能減少每輪SCF過程中要on-the-fly計算的ERI的量來節約耗時。

常有人問為什么他做雜化泛函計算時老是卡在SCF剛開始的位置不動，這是因為SCF第一輪肯定是要計算所有需要考慮的ERI，耗時顯然往往很高，對大體系、大基組的情況花一、兩個小時都是常事，看上去好像一直停在這里似的。如果你給的內存較少而導致大部分ERI都得on-the-fly方式去算，那么之后SCF每一輪都得花將近這么多時間，整個任務跑完可能得半天甚至一天時間。如果PRINT_LEVEL設的是MEDIUM，在SCF第一輪算完后就會輸出Number of sph. ERI's calculated on the fly和Number of sph. ERI's stored in-core，前者是之后每一輪SCF都需要重新算的ERI數，后者是已存到內存中而之后不需要SCF每一輪重新算的ERI數。顯然前者相對于后者越小，之后每一輪SCF比第一輪時間少得越多。

CP2K用戶通常默認用popt版本，每個MPI進程在計算HF交換項時最大內存使用量由&DFT/&XC/&HF/&MEMORY里的MAX_MEMORY控制，單位為MB，扣除掉儲存零碎數據的占用量外，其余都用來儲存ERI。MAX_MEMORY與MPI并行進程數的乘積必須小于當前機子空余物理內存量，顯然應當盡量給大以盡可能減少on-the-fly方式算的ERI量，但也不要太頂著頭分配，因為進入HF交換項計算模塊之前CP2K還會占一定量的內存，故要有適當余量。在SCF第一輪不斷往內存里寫入ERI、內存占用量不斷增大的過程中，若空余內存已用完時CP2K還在繼續往里存ERI，CP2K就會馬上崩潰，甚至還可能導致計算機暫時失去響應。所以MAX_MEMORY應當設得恰到好處。

如果空余物理內存很有限，為了能讓MAX_MEMORY能設得比較大，可以減小MPI進程數。為了不因此浪費CPU運算能力，此時可以用psmp版，此時每個MPI進程下屬的OpenMP線程都可以共享這個MPI進程儲存的ERI。

6 雜化泛函的定義方式

這一節簡要說一下雜化泛函怎么正確在輸入文件里定義。對于非周期性體系，比如用PBE0，直接在&DFT里寫上以下內容即可。此時不用庫侖截斷而用完整的1/r算符，EPS_SCHWARZ為默認的1E-10、SCREEN_ON_INITIAL_P為默認的F。
&XC
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC

顯然對周期性體系絕對不能簡單寫成上面那樣，至少也得自己去指定為使用庫侖截斷。結合前面的講解，周期性體系PBE0計算的泛函定義部分一般寫為下面這樣。PBE0泛函由75%的PBE交換項+25%的HF交換項+100%的PBE相關項構成，所以需要恰當地定義&PBE里面的參數，并在&HF里用FRACTION指定HF成份。
&XC
&XC_FUNCTIONAL
    &PBE
      SCALE_X 0.75
      SCALE_C 1.0
    &END PBE
&END XC_FUNCTIONAL
&HF
    FRACTION 0.25
    &SCREENING
      EPS_SCHWARZ 1E-6
      SCREEN_ON_INITIAL_P T
    &END SCREENING
    &INTERACTION_POTENTIAL
      POTENTIAL_TYPE TRUNCATED
      CUTOFF_RADIUS 6.0
    &END INTERACTION_POTENTIAL
    &MEMORY
      MAX_MEMORY 3000
    &END MEMORY
&END HF
&END XC

有的泛函在&XC_FUNCTIONAL里有特定的字段，可以省得自己去定義細節參數，但HF成份依然必須自己定義，否則相當于只做了殘缺的純泛函計算，結果毫無意義。比如HF成分為10%的TPSSh泛函可以寫為
&XC
&XC_FUNCTIONAL
    &HYB_MGGA_XC_TPSSH
    &END HYB_MGGA_XC_TPSSH
&END XC_FUNCTIONAL
&HF
    FRACTION 0.1
    ...同上
&END HF
&END XC

更多的泛函定義方式可以直接看Multiwfn產生的相應雜化泛函的輸入文件里的寫法，都是嚴格正確的。

7 一個常見警告：The Kohn Sham matrix is not 100% occupied

做雜化泛函計算，特別是用了ADMM時，很容易在SCF開始后看到The Kohn Sham matrix is not 100% occupied警告。具體原因在https://www.cp2k.org/faq:hfx_eps_warning有介紹，這里就不多說了。對于實際計算來說，看到這個警告時，先把&QS里的EPS_PGF_ORB設為比默認明顯更小的1E-12，如果計算正常，即便還有這個警告也可直接無視。EPS_PGF_ORB設得越小，越無法利用重疊矩陣和KS矩陣的稀疏性節約時間，但這對耗時的影響相對于雜化泛函計算的總耗時來說微不足道。目前版本的Multiwfn產生的雜化泛函的輸入文件里自動就把EPS_PGF_ORB設成了1E-12。

對于晶胞較小的情況，本身矩陣稀疏性就差，開發者建議上來就在&QS里設MIN_PAIR_LIST_RADIUS -1從而完全不利用矩陣的稀疏性節約時間，此時絕對不會再出現那個警告，但比起用EPS_PGF_ORB 1E-12時耗時會高一些。

8 簡單例子：C60分子晶體

最后，這里以C60分子晶體的單點任務作為例子演示一個標準的周期性雜化泛函的計算。如果你還不了解Multiwfn創建CP2K輸入文件的功能，建議先閱讀《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。

此體系的cif文件和下面涉及的輸入文件都可以在http://www.shanxitv.org/attach/690/file.zip中找到。此體系結構如下所示，一共240個C原子。晶胞邊長是14.07埃，對于分子晶體來說不擴胞的情況下只考慮gamma點可以接受。

啟動Multiwfn，輸入C60.cif的路徑載入之，然后輸入
cp2k
PBE.inp   //要產生的PBE計算的輸入文件名
0   //基于默認設置產生輸入文件
現在當前目錄下就有了PBE.inp，對應PBE/DZVP-MOLOPT-SR-GTH單點任務。然后接著在Multiwfn里輸入
cp2k
PBE0.inp   //要產生的PBE0計算的輸入文件名
1   //選擇理論方法
-6   //PBE0結合ADMM
0   //產生輸入文件
由PBE0.inp內容可見，Multiwfn自動指定的計算設置完全合適，對應的是PBE0/DZVP-MOLOPT-SR-GTH結合admm-dzp輔助基組、6埃庫侖截斷、EPS_SCHWARZ 1E-6、SCREEN_ON_INITIAL_P T、開啟OT。

手動把PBE0.inp中的WFN_RESTART_FILE_NAME前頭的#去掉并設為PBE-RESTART.wfn，把SCF_GUESS RESTART開頭的#去掉，從而使得PBE0計算時讀取PBE收斂的波函數當初猜。然后根據你的機子的空余物理內存量和要用的MPI進程數設置合適的MAX_MEMORY。

把PBE.inp和PBE0.inp放到同一個目錄下，先運行PBE.inp，算完后當前目錄下就出現了PBE-RESTART.wfn，然后再運行PBE0.inp。在《淘寶上購買的雙路EPYC 7R32 96核服務器的使用感受和雜談》（http://www.shanxitv.org/653）介紹的筆者的2*7R32 512GB雙路服務器上用96核并行、MAX_MEMORY設5000，用CP2K 2023.2 popt版，PBE計算花了8秒，PBE0計算花了178秒。PBE0任務的SCF迭代過程如下所示，可見第一輪耗時較高，由于當前給的內存足夠儲存所有ERI，因而ERI沒有on-the-fly計算的，所以之后每一輪SCF的耗時遠遠低于第一輪。

9 總結

本文把強大的CP2K高效地做雜化泛函計算的各方面要點進行了介紹，以避免初學者對各種關鍵信息毫不知情而憑感覺胡用瞎用。本文的內容也充分體現出用Multiwfn產生CP2K輸入文件非常容易，令研究者從復雜的輸入文件編寫中充分解放，但這絕不意味著CP2K的使用者就可以對關鍵性技術細節毫不知情，否則會各種踩坑、白浪費時間、算出無意義的結果。CP2K用戶必須掌握的知識、要領遠多于大部分其它計算程序。筆者開設的北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）可以使得學員完整系統透徹掌握全部這些知識。

使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像

sobereva — Sun, 04 Jun 2023 16:01:00 +0800

使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡（STM）圖像

文/Sobereva@北京科音 2023-Jun-4

1 前言

強大的波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）有非常方便好用的模擬隧道掃描顯微鏡（STM）圖像的功能，模擬的原理以及Multiwfn中此功能的用法我在《使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像》（http://www.shanxitv.org/549）里有專門說明，但那篇博文是以孤立體系為例。Multiwfn支持對CP2K產生的周期性體系波函數進行大量的分析，在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里我做了很多介紹，其中包括模擬周期性體系的STM圖。本文的目的是舉一個完整的例子，能令甚至之前沒用過CP2K的人都能輕易使用Multiwfn模擬周期性體系的STM圖。沒讀過http://www.shanxitv.org/549的讀者必須要先仔細閱讀此文。不了解Multiwfn的話建議看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。值得一提的是，CP2K自己也有模擬STM的功能，但遠沒有Multiwfn來得好用和方便，所以這里就不多說了。

2D Mater., 6, 015005 (2019)文中給出了幾層厚度的黑磷薄片表面的0.7 V偏壓下的常電流模式的STM圖，如下所示，本文將對黑磷計算同樣條件下的STM圖并與之對照。

本文為了節約時間只用單層黑磷，而且圖省事直接從三維晶體中截出來一層結構直接做單點計算得到波函數。嚴格來說應該用多層模型，并且固定最下層黑磷并令上層的自發弛豫。

本文的例子涉及的文件都可以在http://www.shanxitv.org/attach/671/file.rar中得到。讀者務必使用2023-Jun-4及以后更新的Multiwfn版本，否則情況與本文所述不同。

2 得到單層黑磷的molden文件

如《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）所述，用Multiwfn分析CP2K的波函數必須讓CP2K產生記錄波函數信息的molden文件。這一節我們就對單層黑磷做一個單點計算來得到此文件。

本文文件包里的Phosphorus-black.cif是黑磷的晶體結構，用GaussView打開，只保留中間的單層部分而刪除其它原子，如下所示，然后保存為Phosphorus-black.gjf。

Multiwfn有極其便利的創建CP2K輸入文件的功能，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587），這里我們用Multiwfn創建產生molden文件的單點任務文件。由于CP2K沒法產生考慮k點時的molden文件，因此必須構造超胞用gamma點計算。當前的晶胞參數是a=3.31埃、b=4.38埃，在a方向復制為7倍、b方向復制為5倍后兩個方向都有>20埃，此時只考慮gamma點是可以接受的。當前模擬STM用的偏壓為正，此時能量最低的一個或多個空軌道對STM會產生貢獻，因此必須要求CP2K求解出一些空軌道才行，幾十個就夠了。模擬STM用的偏壓越大應當計算的空軌道越多，因為越高的空軌道會被涉及。

啟動Multiwfn，載入Phosphorus-black.gjf，然后輸入
cp2k //產生CP2K輸入文件
SP.inp //將要產生的CP2K輸入文件名
-7 //設置周期性
XY //XY二維周期性
-2 //產生molden文件
-11 //進入結構編輯界面
19 //產生超胞
7 //a方向復制的倍數
5 //b方向復制的倍數
1 //c方向（垂直于表面的方向）不變
-10 //返回
-9 //其它設置
12 //計算空軌道
20 //算最低20個空軌道
0 //返回
0 //產生輸入文件

此時SP.inp就產生在當前目錄下了，在本文的文件包里提供了，計算級別是Multiwfn默認的PBE/DZVP-MOLOPT-SR-GTH。用CP2K運行SP.inp，在我的2*EPYC 7R32雙路服務器上不到10秒鐘就算完了。之后按照《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）所述，把晶胞信息和價電子信息手動加入.molden文件，即在開頭插入以下字段。本文文件包里的SP-MOS-1_0.molden是已經改好的。

[Cell]
23.17000000     0.00000000     0.00000000
0.00000000    21.90000000     0.00000000
0.00000000     0.00000000    12.11600000
[Nval]
P 5

3 用Multiwfn模擬黑磷常距離模式的STM

首先模擬一下常距離模式的STM。啟動Multiwfn，載入SP-MOS-1_0.molden，然后輸入
300 //其它功能 Part 3
4 //模擬STM
2 //設置偏壓
0.7 //0.7V，和文獻里的相同
0 //開始計算

默認計算的是Z最大值原子上方0.7埃的XY平面，X和Y方向計算范圍和當前晶胞的X和Y方向的跨度一致。瞬間就算完了，在新出現的作圖設置菜單里輸入0繪圖，圖像如下所示

可以再做一些調整，只讓最上層的P用粗體字標注，下層的用細體字顯示以作區分，并且用和文獻里相似的配色，并且調整坐標軸讓刻度整齊。為此，把圖像關閉后，接著輸入
4 //修改顯示原子標簽的距離閾值
0.71 A //在圖上顯示距離作圖平面0.71埃以內的原子的標簽
y //將距離作圖平面更遠的原子的標簽以細體字顯示
9 //修改色彩變化方式
13 //黑-橙-黃
7 //修改標簽間隔
3,3,0.0005
8 //修改色彩刻度
0,0.002 //下限和上限
1 //保存圖像

此時的圖像如下所示

此圖的效果已經很好了。和實驗STM圖像展現出的信息近似一致。

4 用Multiwfn模擬黑磷常電流模式的STM

這次模擬常電流模式的STM，這與文獻里的STM圖對應。啟動Multiwfn，載入SP-MOS-1_0.molden，然后輸入
300 //其它功能 Part 3
4 //模擬STM
2 //設置偏壓
0.7 //0.7V，和文獻里的相同
1 //切換模式為常電流模式
0   //開始計算
3   //繪制STM圖像
0.0009   //常電流值（反復嘗試，取一個圖像效果較好、和實驗圖像較接近的即可）
4 //修改顯示原子標簽的距離閾值
2.6 A //在圖上顯示距離被計算區域頂端2.6埃以內的原子的標簽（計算的區域頂端的z坐標是計算開始之前界面里7 Set range in Z direction選項后面顯示的第二個值。當前設2.6埃可以只把頂層磷原子納入進來）
y //將距離更遠的原子的標簽以細體字顯示
9 //修改色彩變化方式
13   //黑-橙-黃
7 //修改標簽間隔
3,3,0.1
1 //保存圖像
現在看到下圖，和常高模式的STM圖展現出的信息基本一致

5 其它

本文的例子充分體現出用Multiwfn+CP2K繪制周期性體系表面的STM真是又快又方便，鼓勵大家應用在實際當中。請別忘了恰當引用Multiwfn。作為練習，大家可以繪制一下單層MoS2和石墨炔的STM圖。

對于無gap的體系，如金屬表面（包括吸附小分子后），繪制STM有一點需要注意。這種情況一般都是開smearing的，此時費米能級附近的軌道占據數不為整數，導出的molden文件里的軌道占據數也是如此，此時Multiwfn沒法直接繪制STM（Multiwfn會以為是記錄自然軌道的多組態波函數）。解決方法是先進入主功能300的子功能9（計算費米能級功能），此功能會將電子按照軌道能量由低到高以整數方式填充來修改軌道占據數，相當于0 K的占據方式。此外，在此功能里輸入當前的溫度，程序會給你費米能級值，將之記下來。退出此功能后，就可以照常進入STM繪制界面，屆時再把費米能級值輸入進選項3 Set Fermi level即可。另外值得一提的是，用Multiwfn創建CP2K輸入文件時開啟smearing的話，產生的輸入文件里自動就會有計算空軌道的設置，就不用自己再單獨設置了。

CP2K中遇到SCF難收斂時的解決方法

sobereva — Sun, 14 May 2023 13:25:00 +0800

CP2K中遇到SCF難收斂時的解決方法

文/Sobereva@北京科音

First release: 2023-May-15 Last update: 2024-Feb-22

0 前言

經常有人問CP2K中遇到SCF難收斂怎么解決，這是很家常便飯的問題。這里就基于筆者講授的“北京科音CP2K第一性原理培訓班”（http://www.keinsci.com/workshop/KFP_content.html）里面相關部分的內容專門寫個博文，供CP2K用戶參考。這里假定讀者已經了解CP2K計算的基本的背景知識（如k點、CUTOFF、對角化、OT、&BS、smearing等等），如果不了解相關名詞、算法的話，十分推薦參加這個培訓全面、系統學習一遍。

筆者之前專門寫過Gaussian用戶遇到SCF不收斂的解決方法的文章《解決SCF不收斂問題的方法》（http://www.shanxitv.org/61），本文和那篇文章一樣，已經把所有可能有用的方法和要注意的問題說得盡可能全面、充分了，讀者遇到難收斂/不收斂時需要結合實際情況一個個嘗試。切勿問筆者“這文章我看了，但還是解決不了”之類的，我沒有任何可額外補充的（至多是給出當前計算的盡可能詳細的信息，比如給出輸入輸出文件，然后我能回答哪些問題需要優先考慮、哪些方法值得優先嘗試）。

下面第1節先說說遇到SCF不收斂時一般情況需要考慮的問題和解決策略，然后第2、3節再分別說使用對角化和OT時分別可以進一步考慮的。最后再說一下幾何優化、分子動力學時與SCF收斂有關的問題。在此之前先看培訓班ppt了解一下CP2K是怎么判斷SCF收斂的：

1 一般情況時考慮的

? 檢查結構合理性。SCF是否容易收斂和結構的合理性關系極大，化學意義越強的結構通常越容易收斂。因此遇到SCF很難收斂時，需要結合結構化學、固體物理常識檢查結構是否靠譜，如模型是否有意義、是否漏了原子、是否某些原子所處的化學環境和實際明顯不符（如過渡金屬本該與配體成鍵的地方卻光禿禿地懸著）等等。且要結合圖像肉眼檢查，確保不存在不該出現的結構嚴重扭曲、盒子邊緣原子和周期鏡像原子發生過近接觸等問題。若有這些問題，即便SCF收斂了，結果顯然也毫無意義。

? 確認關鍵詞、設置沒明顯不合理的。CP2K的關鍵詞很復雜，對使用者的知識水平要求很高，若用了明顯不當的關鍵詞和設置，SCF很難收斂或收斂到離譜結果是極其正常的事，然而很多要點是初學者難以注意到的。比如計算自旋極化體系的時候必須得在&DFT里寫UKS（或等價的LSD），否則即便自旋多重度設的是>1，程序也會愣當成是閉殼層算。再比如用wavelet的Poisson solver算孤立體系的時候，如果體系有明顯電子密度分布的區域跨越了盒子，結果就無意義。

CP2K屬于那種是完全手動擋還很個性的車。如果經驗不夠豐富的話，筆者建議總是用Multiwfn（http://www.shanxitv.org/multiwfn）產生其輸入文件，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。Multiwfn的CP2K輸入文件創建功能將CP2K盡最大可能包裝成黑箱，不僅用起來極其方便省事，對初學者還能大大減小關鍵詞使用不當的可能性。比如將自旋多重度設成>1時自動就會添加UKS、用wavelet做非周期性計算時會自動設置足夠大又不太浪費的盒子尺寸并且加入自動讓體系在盒子里居中的設置。

? 確保體系的凈電荷、自旋多重度的設置合理。這倆如果設置不當的話，相當于計算了無意義或激發態的電子結構，SCF極難收斂也是很可能出現的事。特別是對于磁性體系，自旋多重度的設置對SCF收斂往往影響極大。

? 如果SCF有收斂趨勢，可增大SCF迭代次數上限再試，即&SCF里的MAX_SCF，其默認的50對于很多情況明顯太小。Multiwfn產生的CP2K輸入文件里默認用128（128是為了優雅，2的7次方。和Gaussian程序一致），絕大多數情況都夠了，但也可能碰上極個別體系仍需要加大。

筆者超級反感一碰到SCF不收斂，也不觀察收斂趨勢，上來就增加迭代次數上限的做法，因為這是最典型的菜鳥思維，絕大多數情況毫無用處而白浪費時間。因為如果SCF過程都發生顯著震蕩了、沒收斂的苗頭，再怎么迭代下去都是白搭，純粹是毫無意義地浪費時間和計算資源。

? 如果發現用其它泛函或基組下計算能收斂的話，可以讀取其收斂的波函數作為初猜，這至少比自動產生的初猜波函數要好、更容易收斂。也即&SCF里用SCF_GUESS RESTART，且&DFT里WFN_RESTART_FILE_NAME指定成之前任務的.wfn文件（若考慮了k點是.kp文件）。注意兩次的基組特征必須定性一致，比如之前計算用的大核贗勢，而現在用小核贗勢或全電子基組，這時候讀之前的波函數就沒意義了。

通常較小的基組下計算比大基組下計算更容易收斂，比如直接用TZV2P檔次的MOLOPT基組不收斂的話，可以嘗試用DZVP檔次的算，若收斂了，讀它收斂的波函數為初猜。

往往HF成份越大的泛函越容易收斂，很大程度是因為其HOMO-LUMO gap更大、SCF過程中占據和空軌道混合程度相對更緩和所致。因此遇上SCF難收斂時可以嘗試用HF成份比實際要用的泛函更大的泛函，若能收斂，讀取其收斂的波函數當初猜。HF成份可以參考《不同DFT泛函的HF成份一覽》（http://www.shanxitv.org/282）。

如果原本用的是純泛函，想嘗試雜化泛函但用不動，也可以用DFT+U試試，耗時和純泛函基本一致，但此時不支持k點，此外怎么獲得+U用的Ueff參數也是個問題。如果實際計算不打算+U，也可以嘗試3或4 eV的+U看看能否令SCF收斂，如果能，則不+U時讀取其收斂的波函數當初猜。

? 做雜化泛函計算時（特別是對于周期性體系）通常會在&DFT / &XC / &HF / &SCREENING里開啟SCREEN_ON_INITIAL_P從而基于初始的密度矩陣做雙電子積分的屏蔽來巨幅節約耗時，Multiwfn產生的雜化泛函計算的輸入文件里默認也開了這個。若此時沒有讀取之前收斂的波函數當初猜（通常是便宜得多的純泛函計算的），則CP2K會基于自動初猜的波函數對應的密度矩陣做積分屏蔽，此時結果通常無意義，SCF能收斂的概率也非常低！然而這個雜化泛函計算的關鍵要點卻經常被初學者所忽略！筆者強烈建議在做雜化泛函計算之前仔細把此文認真讀一遍：《CP2K做雜化泛函計算的關鍵要點和簡單例子》（http://www.shanxitv.org/690）。

? 確保恰當考慮了k點、k點數夠多。如果晶胞邊長不夠大，顯然必須考慮k點而不能只用默認的Gamma點。而且k點設得還得夠多，否則SCF難收斂極為正常。k點如果設置不當，就算SCF能收斂了也沒意義。眾所周知，算能量問題，原則上k點應當大到足以令能量隨k點數增加而基本收斂的程度。

? 嘗試增大&DFT / &MGRID里的CUTOFF、REL_CUTOFF。這倆截斷能設置控制格點精度，不夠大的話不僅結果不準，甚至可能造成SCF難收斂。如果當前明顯已經夠大了，就沒必要再盲目嘗試再增加了。注：還有種說法是嘗試&DFT / &XC / &XC_GRID里寫USE_FINER_GRID T將計算交換-相關泛函部分的格點精度提升為CUTOFF的四倍，但根據我的經驗這并沒有什么實際用處，而耗時暴增。

? CP2K默認的產生初猜波函數的方法是孤立原子密度的疊加，在整個任務一開始時會對每個KIND對應的原子在孤立狀態下做計算，用PRINT_LEVEL MEDIUM時可以看到具體細節。計算孤立原子所用的電子組態可以通過&KIND里的MAGNETIZATION（定義alpha電子數減beta電子數）或者&BS（可以具體控制各個主量子數的各個角量子數上的alpha和beta電子數）來設置。初猜波函數中的孤立原子的組態越接近整個體系收斂的波函數中原子的實際組態則SCF越容易收斂、SCF不收斂的可能性越低。

如果體系有磁性，遇到SCF不收斂時一定別忘了嘗試用MAGNETIZATION或&BS盡可能令初猜波函數中原子組態與實際情況接近，這對SCF是否能夠順利收斂有關鍵性影響。如果當前體系在文獻（可用google搜索來找）或高通量數據庫（如Material Project）中已經被計算過，可以參考那里面報道的原子自旋磁矩或原子組態來設。有時同一種元素在體系里處于不止一種化學環境，差異很大時（比如明顯處于不同價態）可分成多個&KIND分別設置。

即便體系沒有磁性，利用&BS使初猜波函數中的原子帶電狀態使之盡可能符合實際也往往是有益的。比如對于NaCl這樣高度離子性體系計算時SCF不收斂的話，可以嘗試令初猜中Na和Cl的原子組態分別對應Na+和Cl-的狀態（默認情況下分別對應中性的Na和Cl）。曾經我用PBEsol/pob-TZVP-rev2計算10層的NaCl板時通過這種做法成功解決了SCF不收斂問題。

? 對gap較小的體系可嘗試使用能級移動，即&SCF / LEVEL_SHIFT，設比如0.3、0.4左右，默認單位是Hartree。這使得SCF過程中人為加大空軌道與占據軌道之間的能級差，從而減輕彼此間混合的劇烈程度而有時令SCF更容易收斂。

? 嘗試將&QS / EPS_DEFAULT減小幾個數量級，默認為1E-10。EPS_DEFAULT這個參數控制一堆涉及數值計算精度的參數，如果設得太大，有可能造成SCF難收斂。若這個數值當前已經很小了，比如都1E-12了，嘗試繼續減小則沒必要。

? 改用其它泛函、基組。如果原先的計算級別就是極難收斂，而有其它級別性價比、精度相仿佛，且發現能收斂，那就別死磕當前的了。

? 在GPW與GAPW之間切換。如果之前用的GPW結合GTH贗勢基組不收斂，嘗試用GPAW結合一個全電子基組計算試試，反之亦然。

? 修改默認產生波函數的方式。參見手冊里&SCF / SCF_GUESS設置。但其實這沒什么用，除了默認的ATOMIC（即前述的基于孤立原子密度疊加產生的）以及RESTART（讀取之前的波函數當初猜），其它的一個能打的都沒有。

用ATOMIC產生初猜時，CP2K計算每個KIND對應的孤立原子也是需要SCF迭代的，上限是固定的50輪（沒法靠關鍵詞改），在PRINT_LEVEL MEDIUM時其迭代過程會在一開始輸出。有時候我發現到了50輪時收斂精度還很差，故有可能在某些情況下由于這個問題導致初猜質量爛，不利于實際體系的SCF收斂。我發現可以通過atom_kind_orbitals.F里atom%optimization%max_iter修改迭代次數上限，比如改成300然后重新編譯，一般孤立原子的SCF就都能收斂了。

? 如果加了外電場時SCF很難收斂，可能是外電場過大所致，可嘗試更小的數值。如果就是要用較大的外電場，可以試試讀取較小電場或不帶電場時收斂的波函數當初猜。有的體系在很強外電場下本來就不可能穩定存在，這時候SCF不收斂自然很正常，是物理上決定的，也別盲目去計算。

? 實在沒轍時干脆放寬SCF收斂限，即加大&SCF / EPS_SCF，默認是1E-5（這其實已經挺寬松了），這是走投無路時的做法。或者換程序，如免費且和CP2K同樣流行的Quantum ESPRESSO，或者嘗試CP2K里的純粹基于平面波的SIRIUS模塊。某些體系用CP2K的quickstep模塊就是很難收斂，如某些d族純金屬的slab，這和算法、基函數特征有關。

2 用對角化的情況可專門考慮的

? 對無gap（如金屬、石墨烯）或gap很小的體系（非常窄帶半導體）應當使用smearing。即加入&SCF / &SMEAR / METHOD FERMI_DIRAC，且在&SCF里用ADDED_MOS要求求解一些空軌道從而能被熱激發的電子所占據。ADDED_MOS數目夠用就行，設太大了會增加耗時而且也沒用處。Multiwfn產生CP2K輸入文件時如果開了smearing會自動做ADDED_MOS設置，并根據體系原子數估計一個通常夠用的值。用較小電子溫度（如200K、300K）時若不收斂可嘗試再提升溫度，即&SMEAR / ELECTRONIC_TEMPERATURE設的。溫度設高時ADDED_MOS往往需要相應地增加，否則熱激發的電子沒法按Fermi-Dirac分布占該占的空軌道（此時計算時會有警告）。

我老在網上看到有人不理解smearing原理是什么就胡用、亂試smearing以試圖解決SCF不收斂，所以這里強調一下。一定要注意，對于gap不是很小的體系（比如氧化鋰等絕緣體），根本沒必要嘗試smearing，本來這在原理上就不是smearing能起到幫助SCF收斂的情況，而且本文還有其它一堆方法可以嘗試用于幫助SCF收斂，完全輪不到考慮用smearing。

用smearing還會對結果產生一定影響，結果會依賴于電子溫度T的設置，用Fermi-Dirac函數做smearing時最終給出的能量相當于是T溫度下電子的自由能。如果你就是要得到基態（0K）的電子能量，原則上最好做T→0的外推。如果你不做外推，在能量橫向對比時至少也應當用相同的smearing設置。關于smearing的物理意義和實際意義，以及做外推的方法，我在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里有專門詳細的講解。

? 修改&SCF / &MIXING / METHOD設置SCF過程中將舊密度矩陣與新密度矩陣混合得到實際每一輪密度矩陣的方法。默認的DIRECT_P_MIXING往往表現極差。通常用BROYDEN_MIXING，這也是Multiwfn產生的輸入文件里默認的。亦可嘗試PULAY_MIXING，但一般不比BROYDEN_MIXING好。

? 對于使用BROYDEN_MIXING的情況，嘗試修改&MIXING中的ALPHA（新密度矩陣混入舊的的比例，一般用默認的0.4，即混入40%。難收斂時可嘗試諸如0.1、0.2、0.3）和NBROYDEN（默認的4明顯太小，改大到8通常很有助于收斂，這是Multiwfn產生的輸入文件默認的。也可以嘗試進一步改大到比如12）

? 對于使用PULAY_MIXING的情況，嘗試修改&MIXING中的NPULAY （默認的4偏小，改大往往很有助于收斂）。注：NPULAY、NBROYDEN、NBUFFER三個關鍵詞實際上等價。

? 對于使用DIRECT_P_MIXING的情況，收斂到&SCF / EPS_DIIS值后會切換到加速收斂的DIIS算法。但默認的EPS_DIIS為0.1，有點偏小，可以嘗試改大以早點啟用DIIS（但借助DIIS通常仍不如BROYDEN_MIXING）。還可以把控制最大DIIS矢量的&SCF / MAX_DIIS從默認的4改大到諸如7。

? 在對角化框架下各種方法都不行的話，嘗試改用OT。如果OT能收斂的話，而且當前由于特殊目的又必須用對角化，也可以用OT收斂的波函數當對角化計算的初猜。由于OT不支持k點，因此如果之前對角化是結合k點算的較小晶胞情況，用OT時就只能靠超胞模型了。OT也不支持用ADDED_MOS設置在SCF計算過程中求解空軌道，因此必須用smearing的導體體系就別指望OT了。如果是用GFN1-xTB，總應當用OT，這也是Multiwfn默認的。

3 用OT的情況可專門考慮的

? 如果之前沒用outer SCF則開啟之，這對解決SCF不收斂非常有效，也即加入&SCF / &OUTER_SCF字段。由于其重要性，Multiwfn產生的輸入文件里開OT后默認就是啟用outer SCF的。此時每次到達&SCF / MAX_SCF所設的inner SCF的步數上限時，就會進入outer SCF過程更新OT的preconditioner，然后進入下一次的inner SCF迭代過程。每次更新preconditioner后經常會看到比之前以更快的速度趨近收斂限。原理上來說，OT每一步都更新preconditioner是最理想、SCF最容易收斂的，但由于這樣太昂貴，因此默認情況下（即不用outer SCF時）只會在最開始計算一次preconditioner。

開啟outer SCF時inner SCF的迭代次數上限應當設得比平時更小（否則經歷很多次SCF迭代才做一次outer SCF，無法有效達到用outer SCF的目的），15至35是比較合適的范圍。outer SCF迭代次數上限通過&SCF / &OUTER_SCF / MAX_SCF控制，總SCF迭代次數上限相當于inner SCF次數上限和(outer SCF次數上限+1)的乘積。

&SCF / EPS_SCF設的是inner SCF收斂限，&SCF / &OUTER_SCF / EPS_SCF設的是outer SCF收斂限，后者應當小于或等于前者（一般設為相同）。

? 嘗試其它的preconditioner。FULL_ALL通常最好，但最貴。大體系可以嘗試明顯更便宜的FULL_SINGLE_INVERSE和FULL_KINETIC。個別時候FULL_KINETIC反倒比FULL_ALL更容易收斂。GAPW計算總應當用FULL_ALL。

? 修改OT算法，把&OT里的ALGORITHM從默認的STRICT改為IRAC，SCF經常會收斂得更穩健，非常值得嘗試！特別是我發現用pob系列全電子基組的時候，用IRAC往往容易收斂得多，因此Multiwfn產生基于pob系列基組的輸入文件時若用戶選了OT，則默認用IRAC。而且用約束性DFT（CDFT）時，我發現用IRAC也會令SCF容易收斂得多得多。

? &OT / MINIMIZER控制OT用的極小化方法，若之前用的是較快的DIIS（也是Multiwfn產生的用OT的輸入文件里默認的），可改為更穩健但耗時更高的CG（共軛梯度法），也可以嘗試BROYDEN。

? MINIMIZER用CG時，嘗試用&OT / LINESEARCH選項把線搜索方式設為比默認2PNT更貴但也更穩的3PNT。

4 幾何優化和分子動力學過程中的SCF收斂問題

在幾何優化過程中，有可能中途偶爾出現幾步SCF不收斂。此時程序還會繼續算下去。只要最終幾何優化收斂時SCF也是收斂的狀態就行了。（注：從CP2K 2024.1開始，優化中途SCF不收斂會導致任務自動停止，&SCF里寫IGNORE_CONVERGENCE_FAILURE才會在SCF不收斂的情況下也繼續優化）

有時候在幾何優化初期SCF一直不收斂，不要急于把任務停掉，可以先湊合跑跑看看情況。初始幾何結構相對于實際的極小點結構而言肯定是存在多多少少不合理性的，此時SCF也肯定比極小點結構處更難收斂。在SCF不收斂、受力也較糙的情況下湊合跑一定步數后，結構一般會變得比一開始更為合理，此時SCF也往往變得更容易收斂。

在幾何優化、分子動力學過程中，程序會自動利用之前一些結構的波函數外推出當前步的初猜波函數，這比起直接產生的初猜波函數質量通常好得多。外推方式可以通過&DFT / &QS / EXTRAPOLATION來具體設置，默認的ASPC通常是最好的選擇。由于這一點，當幾何優化、分子動力學跑一定步數之后，SCF達到收斂限所需的迭代次數明顯減小、出現SCF不收斂的可能性也大為降低。

注意當使用k點時CP2K不支持外推產生初猜波函數，因此會在幾何優化、分子動力學等涉及結構變化的任務中每一步重新產生初猜波函數，即EXTRAPOLATION USE_GUESS的情況。我一般建議改為EXTRAPOLATION USE_PREV_P，代表直接使用上一步SCF收斂的波函數當初猜。

如果想減小MD中途出現SCF不收斂的可能，可以嘗試減小步長（&MD / TIMESTEP）。如果想減小幾何優化中途出現SCF不收斂的可能，可減小優化的步長上限，即置信半徑（如通過&MOTION / &GEO_OPT / &BFGS / TRUST_RADIUS）。這樣做之后由于相鄰兩個結構間差異變小了，ASPC或USE_PREV_P給出的初始波函數對于當前結構而言就越理想，自然SCF收斂也會更容易，且有利于讓SCF能收斂的狀態保持下去。

sobNEB：產生CP2K的NEB的插點的方便的工具

sobereva — Tue, 28 Feb 2023 02:43:09 +0800

sobNEB：產生CP2K的NEB的插點的方便的工具

文/Sobereva@北京科音 2023-Feb-28

1 介紹

CP2K有很好用的NEB功能可以得到反應路徑，靠CI-NEB還可以同時得到過渡態結構，用法在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里講得很詳細。雖然CP2K的NEB功能自己就能自動線性插點，但是偶有bug（有的初始間距是其它的兩倍），而且沒法人為控制自己提供的中間結構對應哪個點，也不便于預覽初始點。

我開發的sobNEB在插點時可控性更強。雖然也是線性插點，但有以下好處
(1)可以直接控制自己提供的結構對應哪個初始點，其它點都在相鄰的提供的結構之間線性插入
(2)可直接產生對應CP2K的&BAND部分中的信息便于寫NEB輸入文件
(3)會將用戶提供的點和所有額外插的點的結構一起寫入traj.xyz文件以在VMD中觀看所有初始點構成的軌跡，以便判斷插點是否合適。

sobNEB的下載地址：http://www.shanxitv.org/soft/sobNEB_1.0.zip。有.exe后綴的是Windows版，無后綴的是Linux版。

注意sobNEB不會對輸入的結構自動消除平動轉動。

2 用法

用戶只需編輯sobNEB目錄下的sobNEB.ini，指定自行提供的各個xyz文件的路徑以及對應的點的序號（之間以冒號分隔），在啟動程序后就會開始處理。只有始、末端的結構是必須自己提供的，中間的結構可以不提供也可以提供一個或者多個。

下面是sobNEB.ini文件的例子，涉及到的三個xyz文件在sobNEB目錄下都提供了。conf1.xyz和conf2.xyz是分子的兩個構象，分別作為始、末端。TS_guess.xyz是自己擺的過渡態初猜結構，其點號位于始末、端點號的正中央。

conf1.xyz : 1
TS_guess.xyz : 11
conf2.xyz : 21

啟動sobNEB后會看到以下信息

sobNEB: Generate interpolated .xyz files and CP2K input for NEB calculation
Programmed by Tian Lu (sobereva[at]sina.com)
Version 1.0, release date: 2023-Feb-19

Loading sobNEB.ini...
3 structures are given

Generating interpolated structures between given systems 1 and 2
Generating point   2 : 2.xyz
Generating point   3 : 3.xyz
Generating point   4 : 4.xyz
Generating point   5 : 5.xyz
Generating point   6 : 6.xyz
Generating point   7 : 7.xyz
Generating point   8 : 8.xyz
Generating point   9 : 9.xyz
Generating point 10 : 10.xyz

Generating interpolated structures between given systems 2 and 3
Generating point 12 : 12.xyz
Generating point 13 : 13.xyz
Generating point 14 : 14.xyz
Generating point 15 : 15.xyz
Generating point 16 : 16.xyz
Generating point 17 : 17.xyz
Generating point 18 : 18.xyz
Generating point 19 : 19.xyz
Generating point 20 : 20.xyz

Done! .xyz file of each point has been generated in current folder
Corresponding &BAND field of NEB task of CP2K has been exported to BAND.txt in
current folder
traj.xyz is the trajectory file containing all points for previewing purpose

可見產生了2-10號、12-20號，共18個插點結構，當前目錄下出現了相應序號的xyz文件。

當前目錄下出現的BAND.txt內容如下，可以直接粘貼到&BAND字段里
    &REPLICA
      COORD_FILE_NAME conf1.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 2.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 3.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 4.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 5.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 6.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 7.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 8.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 9.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 10.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME TS_guess.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 12.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 13.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 14.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 15.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 16.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 17.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 18.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 19.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME 20.xyz
    &END REPLICA
    &REPLICA
      COORD_FILE_NAME conf2.xyz
    &END REPLICA

可以用VMD檢查當前目錄下同時產生的traj.xyz。里面21幀同時疊加顯示時如下所示，可見初始點的分布是合理的，在兩個構象間變化，可以用于NEB計算。

使CP2K計算的振動模式可以被GaussView觀看的程序：MfakeG

sobereva — Sun, 22 Jan 2023 21:01:00 +0800

使CP2K計算的振動模式可以被GaussView觀看的程序：MfakeG

文/Sobereva@北京科音 2023-Jan-22

GaussView在觀看振動模式方面非常好用，可以方便地顯示振動矢量，播放和保存振動動畫，還可以沿特定的振動模式對結構進行位移，但GaussView只支持Gaussian程序的振動分析輸出文件。之前我寫過一個程序OfakeG，見《OfakeG：使GaussView能夠可視化ORCA輸出文件的工具》（http://www.shanxitv.org/498），可以把ORCA程序振動分析的輸出文件轉化成類似Gaussian的格式，從而能通過GaussView來可視化。CP2K是非常強大且免費的第一性原理程序，為了也能借助GaussView便利地觀看其振動模式，我寫了叫MfakeG的程序，在此進行介紹。筆者講授的北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里非常詳細講授怎么用CP2K計算分子和周期性體系的紅外和拉曼光譜，其中也會充分利用MfakeG程序。

MfakeG可以在主頁http://www.shanxitv.org/soft/MfakeG免費下載，Windows和Linux版都有。

CP2K的振動分析任務的輸入文件里在&VIBRATIONAL_ANALYSIS字段中加入以下內容就可以在振動分析結束時輸出后綴為.mol的Molden文件，默認文件名為[項目名]-VIBRATIONS-1.mol。
&PRINT
&MOLDEN_VIB
&END MOLDEN_VIB
&END PRINT
如果你用Multiwfn按照《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）介紹的方式產生CP2K振動分析的輸入文件，默認也會產生.mol文件。注意這個.mol文件和常見的記錄分子結構和鍵連關系的.mol文件完全是兩碼事。

.mol文件里面有許多字段，記錄了元素、坐標、振動頻率、描述振動矢量的正則坐標、紅外強度。MfakeG干的事情就是把.mol轉換為GaussView能認的類似Gaussian振動分析輸出文件的形式。

啟動MfakeG后，輸入.mol文件的路徑，就會在相同目錄下產生與之同名但帶了-fake后綴的.out文件，之后載入到GaussView里就可以照常用Results - Vibrations界面觀看振動信息了。

CP2K一般都是用來算周期性體系的。為了能讓GaussView顯示出來晶胞邊框，對周期性體系需要自行編輯.mol文件，在里面第二行插入[Cell]字段，比如
[Cell]
19.25142628     0.00000000     0.00000000
-9.62562669    16.67336579     0.00000000
0.00000000     0.00000000    15.00000000
三行內容是晶胞的三個矢量，單位為埃。這樣MfakeG轉出來的偽Gaussian輸出文件的原子坐標部分最后會多出來三個原子信息，對應晶胞矢量。

在MfakeG的example目錄下freq.inp是CP2K對GaN超胞做振動分析的輸入文件，算完后產生了freq-VIBRATIONS-1.mol，用MfakeG轉換后就是freq-VIBRATIONS-1_fake.out。用GaussView的振動分析界面看到的是下面的效果，可見效果很好。

實際上筆者原本是打算把CP2K的輸出文件轉化成偽Gaussian格式的。但之所以后來決定轉化.mol格式，是因為其格式比CP2K輸出文件更簡單，而且這樣更有通用性，讀者可以自己寫個小程序把任意其它計算化學程序做振動分析得到的結果寫成.mol格式，之后都可以借助MfakeG用GaussView觀看。

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

sobereva — Wed, 31 Aug 2022 22:27:00 +0800

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

文/Sobereva@北京科音

First release: 2022-Aug-31 Last update: 2024-Feb-1

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文件的情況，有以下兩種做法：
? 方法1：基于OT收斂的波函數做對角化
先做一次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。

? 方法2：讓Multiwfn讀取輸出文件里的軌道能量
可以做OT時要求CP2K把軌道能量計算出來并輸出到輸出文件里，然后用Multiwfn（必須用2023-Oct-14及之后的版本）載入molden文件后，在主菜單里輸入cp，再輸入7，Multiwfn就會問你讀取哪種情況輸出的軌道能量（見下），并讓你輸入CP2K輸出文件的路徑，然后軌道能量就被讀入Multiwfn了。CP2K有以下兩種方式將軌道能量輸出到輸出文件里，都是加入到&DFT/&PRINT字段中：
(1)加入以下信息后，SCF收斂后會輸出全部求解出的軌道能量（用對角化時如果用了ADDED_MOS還會輸出空軌道能量）
&MO
ENERGIES T
OCCUPATION_NUMBERS T
&EACH
QS_SCF 0
&END EACH
&END MO
(2)加入以下信息后，SCF收斂后會顯示所有占據軌道能量和最低3個（可自己指定）空軌道的能量。即便對于SCF過程中無法求解空軌道的OT也可以用這種方法得到空軌道能量。但注意此時用OT產生的molden文件里仍然只有占據軌道的波函數，也因此沒法用Multiwfn對空軌道部分繪制PDOS（計算PDOS需要計算空軌道的軌道成份，這需要空軌道的波函數。因此此時只能是在TDOS曲線里體現空軌道）
&MO_CUBES
NHOMO 1
NLUMO 3
WRITE_CUBE F
&END MO_CUBES

如果你想檢查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)注意基組的選擇

使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線

sobereva — Thu, 10 Mar 2022 10:51:00 +0800

使用CP2K結合Multiwfn繪制密度差圖、平面平均密度差曲線和電荷位移曲線

文/Sobereva@北京科音

First release: 2022-Mar-10 Last update: 2022-Mar-18

1 前言

電子密度差（electron density difference, EDD）是分析成鍵、電子轉移非常常用的手段。波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）在繪制和分析密度差方面非常靈活、方便、強大，筆者之前有不少文章與此相關，見《使用Multiwfn作電子密度差圖》（http://www.shanxitv.org/113）、《使用Multiwfn考察分子軌道、NBO等軌道對密度差、福井函數的貢獻》（http://www.shanxitv.org/502）、《使用Multiwfn計算激發態之間的密度差》（http://www.shanxitv.org/429）、《使用Multiwfn做電子密度、ELF、靜電勢、密度差等函數的盆分析》（http://www.shanxitv.org/179）等。之前筆者寫的密度差有關的文章主要是對于孤立體系寫的，本文的目的之一是介紹怎么將非常流行、免費強大的CP2K第一性原理程序與Multiwfn相結合，對周期性體系繪制密度差圖。在此順帶特別強調，electron density difference的中文是密度差，“差分密度”是對密度差的嚴重不當稱呼，絕對不要用這個詞！

本文另一個目的是介紹怎么通過Multiwfn繪制平面平均密度差（plane-averaged EDD）曲線，它是密度差的衍生。密度差是個三維空間函數，經常通過繪制平面圖、等值面圖方式圖形化考察。雖然這兩種圖很直觀，但平面圖會受到截面選取的影響，而等值面圖會受到等值面數值（isovalue）選取的影響，因此都有一定任意性，也不容易定量討論。對于電子轉移有比較明確方向的情況，可以繪制平面平均密度差曲線，由此可以清楚地考察垂直于選取的方向上每個截面上的密度差積分值，在定量討論、對比上比較便利。例如，固體表面上吸附一個小分子，吸附導致垂直于表面方向有明顯的電子轉移，就可以在垂直于表面的方向上繪制平面平均密度差曲線來準確、細致地考察不同截面處的電子凈增、減情況。還有一種圖叫做電荷位移曲線（charge displacement curve），它相當于平面平均密度差曲線的積分曲線，對于定量討論電荷轉移量情況很有幫助，本文也會介紹怎么通過Multiwfn繪制這種圖。

本文使用的例子是NaCl板吸附一個CO分子，NaCl和CO將被定義為兩個片段，通過密度差考察吸附導致的電子轉移情況。本文文件包里的opt目錄下的NaCl_CO-1.restart是CP2K在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH級別下對這個體系做優化產生的restart文件，我們將基于這個結構繪制密度差。為簡化問題，在優化的時候整個NaCl板的坐標被凍結為了晶體坐標的狀態，因此不涉及表面重構。

下文涉及到的相關文件可以在http://www.shanxitv.org/attach/638/file.rar下載。cub文件由于體積太大，文件包里就不提供了。

本文的計算使用CP2K 9.1，安裝方法見《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。VMD使用的是1.9.3版，可在http://www.ks.uiuc.edu/Research/vmd/免費下載。讀者務必使用Multiwfn最新版，如果不知道是不是最新的就現在立刻去主頁http://www.shanxitv.org/multiwfn下載一個。如果對Multiwfn不了解，看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。不了解cub文件的話可以看看《Gaussian型cube文件簡介及讀、寫方法和簡單應用》（http://www.shanxitv.org/125）。

2 計算密度差格點數據

這一節我們要得到整個體系、NaCl板和CO的各自的電子密度的cub文件，然后由Multiwfn求差。

我們先創建一個計算整個體系的電子密度cub文件的CP2K輸入文件。啟動Multiwfn，然后輸入
NaCl_CO-1.restart //在本文文件包里
cp2k //創建CP2K輸入文件。在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）一文里有詳述
NaCl_CO.inp //產生的輸入文件的路徑

現在進入了創建CP2K輸入文件的界面。為了之后繪制的平面平均密度差曲線比較好觀看，最好把當前體系挪到盒子中央附近。先看一眼當前的結構是什么樣，輸入
-11 //進入幾何操作界面。此界面的詳細介紹見《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）
0 //觀看當前結構

現在出現了圖形窗口。點擊菜單上的Other settings - Toggle showing cell frame把盒子邊框顯示出來，再選Other settings - Toggle showing all boundary atoms把邊界原子顯示出來。現在看到的圖如下

可見當前整個體系在盒子底部，NaCl有三層，由于最下層的原子正好在Z=0的位置，因此這些邊界原子在盒子頂端也顯示了出來。點擊右上角的RETURN按鈕返回。為了把體系挪到盒子中央去，輸入
24 //平移體系使得選中的部分在晶胞中居中
1-110 //當前體系里所有原子。也可以直接按回車選中整個體系
0 //再次觀看體系
可見此時體系確實在盒子中央了

關閉圖形窗口，接著輸入
-10 //返回CP2K輸入文件創建界面
-7 //設置周期性
XY //XY方向二維周期性
-3 //要求產生cube文件
1 //電子密度
0 //產生CP2K輸入文件

現在當前目錄下就有了NaCl_CO.inp。讀者可根據要求修改里面的CUTOFF等設置，本例不那么講究，就用默認的。運行之，算完后當前目錄下就有了NaCl_CO-ELECTRON_DENSITY-1_0.cube，將之改名為NaCl_CO.cub。

將NaCl_CO.inp復制為NaCl.inp，用文本編輯器打開，把其中PROJECT NaCl_CO改為PROJECT NaCl，把&COORD里面C和O原子刪掉。用CP2K運行NaCl.inp，將得到的cub文件改名為NaCl.cub。

將NaCl_CO.inp復制為CO.inp，用文本編輯器打開，把PROJECT NaCl_CO改為PROJECT CO，把&COORD里面Na和Cl原子刪掉。用CP2K運行CO.inp，將得到的cub文件改名為CO.cub。

現在用Multiwfn對三個cub文件求差。啟動Multiwfn，然后輸入
NaCl_CO.cub //輸入此文件實際路徑
13 //處理格點數據
11 //格點數據運算
4 //減去另一個格點數據
NaCl.cub //輸入此文件實際路徑
11 //格點數據運算
4 //減去另一個格點數據
CO.cub //輸入此文件實際路徑

現在內存里的格點數據就已經是密度差格點數據了。想觀看一下等值面的話，可以進入選項-2 Visualize isosurface of present grid data。由于物理吸附對應的密度差的數量級往往很小，所以要用比較小的等值面數值才能清楚看到。在圖形窗口右側把等值面數值設為0.0005后，看到的圖像如下（點save picture按鈕保存的圖像比直接在圖形窗口里看到的更好）

綠色和藍色等值面分別代表格點數據數值為正和為負的部分。當前密度差是按照ρ(NaCl...CO)-ρ(NaCl)-ρ(CO)算的，顯然綠色體現吸附后電子密度增加的部分，藍色體現電子密度減少的部分。點擊右上角RETURN關閉圖形窗口。

若想得到密度差的cub文件，就選0 Export present grid data to Gaussian-type cube file (.cub)，輸入EDD.cub，然后內存里的格點數據就被導出為了當前目錄下的EDD.cub。

使用《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》（http://www.shanxitv.org/483）里提供的我寫的腳本，在VMD程序里只需要輸入一行命令cub EDD 0.0004，就可以把EDD.cub的0.0004等值面顯示出來（如果沒顯示出來，一個字一個字看http://www.shanxitv.org/483這篇文章）。在文本窗口里輸入pbc box把盒子顯示出來，再改成正交視角，經過Tachyon渲染看到的就是下面的圖。圖中藍色小球是Na，棕色是Cl，等值面的藍色和綠色和前面Multiwfn里的定義相同。

從此圖可以明顯看出，CO的碳結合到Na上后，C和O上的電子密度整體大幅減少，很大程度往Na的方向上轉移了。這也容易理解，當前體系里Na帶顯著的正電荷，CO結合上去，特別是C還帶著豐富的孤對電子，必然其電子會被往Na那邊顯著地極化。

再來看看把等值面數值設得更小的情況。在VMD文本窗口里輸入cubiso 0.0001把等值面數值調為0.0001，稍微調整視圖后看到的圖如下

可見此時的密度差圖展現了比之前更多的信息。粉色虛線圈住的地方電子密度有一定降低，由于僅在等值面數值很小的時候才能顯現出來，所以這部分降低量很少（想通過積分得到確切值的話需要用Multiwfn的盆分析功能，前面http://www.shanxitv.org/179文中介紹過）。這部分電子密度降低是因為CO向Na這邊轉移了電子，由于交換互斥作用，使得原本這部分的電子被排斥開，向更遠處轉移。

3 繪制密度差局部積分曲線和平面平均密度差曲線

這一節繪制密度差局部積分曲線和平面平均密度差曲線，首先了解一下定義。某個三維空間函數的局部積分曲線（local integral curve）可以對任意方向繪制，如果對比如Z方向繪制，則表達式為

其中p是被考察的函數。如果p取為密度差的話，那么這個曲線就是前述的密度差局部積分曲線了。

還有個概念叫積分曲線（integral curve），就是對局部積分曲線再進一步積分。例如對Z方向如下所示

其中積分的起點z_ini由自己決定。如果被考察的函數是密度差，則這個積分曲線對應的就是前述的電荷位移曲線。

平面平均曲線（plane-averaged curve）定義如下，其中A_XY是格點數據對應的盒子在XY方向的面積。這種曲線體現的是被考察的函數在特定截面上的平均值。

在Multiwfn里可以基于格點數據對任意函數繪制上述三種曲線，需要先把相應的格點數據存到內存里。比如可以通過主功能5格點數據運算功能算出來后自動存到內存里，也可以啟動時從cub等格點數據文件里直接載入，等等。

下面繪制密度差局部積分曲線。如果你之前已經關了Multiwfn，就再啟動Multiwfn，載入密度差格點數據文件EDD.cub，然后進入主功能13，再進入子功能18，這是用來繪制（局部）積分曲線和平面平均曲線的功能。如果你按照上一節操作之后還沒關Multiwfn，則內存里的格點數據正是密度差格點數據，因此直接進入子功能18就行了。

在子功能18里，Multiwfn先問你對哪個方向繪制，我們想垂直于NaCl界面繪制，因此輸入Z，因為NaCl界面平行于XY方向。之后Multiwfn問你繪制范圍的下限和上限位置是什么，此例想對整個Z范圍進行繪制，因此按照屏幕上的提示直接輸入a。之后會看到一個菜單，通過相應選項可以直接繪制積分曲線、局部積分曲線、平面平均曲線，可以保存它們的圖像，還可以把曲線數據導出為文本文件便于用Origin、gnuplot、qtgrace等第三方程序繪圖。

現在選擇屏幕上的選項2 Plot graph of local integral curve，由于當前內存里的格點數據是密度差，所以下圖就是密度差局部積分曲線

讀者如果想讓等值面圖和密度差局部積分曲線便于準確對照，可以用屏幕上相應選項讓Multiwfn導出局部積分曲線數據為當前目錄下的locintcurve.txt（每一列的含義看導出時屏幕上的提示），然后放到Origin等程序里作圖。之后在VMD里用正交視角，把盒子邊框顯示出來，把保存的等值面圖和局部積分曲線圖放到一起，把前者的白色背景設為透明色或者摳掉，再調節位置和大小讓兩張圖邊框位置正好對上，如下所示。不會用Photoshop的話用版本不很老的PowerPoint也可以完成

再改一下透明度，裁切一下邊緣去掉盒子邊框，最終得到下面的圖，效果很好。由這個圖可以比較準確地捕捉各個截面上電子的凈增減，可見在C與表層的NaCl相接的地方電子密度由于相互作用而增加明顯（但絕對的數量級并不算大）。

如果想要繪制平面平均密度差曲線，前面選2 Plot graph of local integral curve的地方改成選3 Plot graph of plane-averaged curve即可，其它不變。

4 繪制電荷位移曲線

在Multiwfn當前界面里選1 Plot graph of integral curve，由于當前內存里的格點數據是密度差，因此得到的積分曲線的就是電荷位移曲線了：

上圖橫軸左端對應于NaCl...CO體系Z=0的位置，橫軸右端對應于NaCl...CO體系Z=25埃（盒子Z尺寸）的位置。這個圖是對上一節的密度差局部積分曲線從左向右積分產生的，從肉眼可以清楚看出來對應關系。曲線最右邊數值為0，體現了密度差全空間積分為0，也即CO與NaCl結合沒有導致總電子數有凈變化。

為了便于觀看，將電荷位移曲線和結構圖+等值面圖疊加起來。我也把曲線最高的峰的橫、縱坐標位置標注上了。電荷位移曲線縱坐標是無量綱的，因為對應的是電子數。

這張圖體現出，在Z<17.1埃的區域里，電子數凈增加了0.036。如果姑且把這個位置作為CO和NaCl的分界，那么說成是CO向NaCl凈轉移了0.036個電子是可以的。對于其它吸附體系，未必在表面和被吸附物之間恰好有一個峰，讀者可以自己決定哪個Z值可作為兩部分的分界面（這有一定任意性），從而從電荷位移曲線上讀出電子凈轉移量。考察電荷位移曲線顯然不是唯一的判斷電荷轉移量的方法，更普適、更常用的是計算片段電荷，即某個片段里所有原子電荷的加和（結果明顯受原子電荷計算方法的選取所影響。原子電荷計算方法有很多，比如基于CP2K的波函數Multiwfn可以計算CM5、Hirshfeld、Mulliken等電荷）。

5 對其它函數計算（局部）積分曲線

上面介紹的Multiwfn計算（局部）積分曲線和平面平均曲線對任何三維空間函數都適用，這里順帶再舉個例子。比如想考察NaCl...CO這個體系不同Z位置的電子密度分布，就啟動Multiwfn，依次輸入
NaCl_CO.cub
13 //處理格點數據
18 //繪制（局部）積分曲線
Z //順著Z方向繪制
a //整個Z范圍
-1 //將橫坐標單位切換為埃
2 //繪制局部積分曲線
由于載入的NaCl_CO.cub里記錄的是電子密度，所以此時看到的下面的局部積分曲線體現了對應不同Z值的XY截面上的電子量。中間三個很高的峰正對應于NaCl板的每一層，18.8埃附近的小峰來自于CO的電子。

可見，不管什么函數，Multiwfn都能畫（局部）積分曲線或平面平均曲線，只要提供記錄相應函數的格點數據文件就可以。諸如自旋密度也同樣可以繪制這些曲線，比如你有個反鐵磁性體系，位于不同層上的原子自旋磁矩情況各有不同，就可以對自旋密度繪制局部積分或平面平均曲線，清楚地考察電子自旋分布情況。對于靜電勢也可以繪制平面平均曲線，特別是對于界面體系，在垂直于界面方向的這種圖是計算功函數的關鍵。用戶只要給Multiwfn提供靜電勢的格點數據文件即可，比如流行的CP2K就可以直接產生，即在Multiwfn產生CP2K輸入文件的界面里選擇導出cube文件的內容的時候選4 Hartree potential (negative of ESP)。

在http://bbs.keinsci.com/thread-29906-1-1.html一貼里有VASP用戶通過本文介紹的功能基于VASP計算的CHGCAR和LOCPOT文件，用Multiwfn分別考察了平面平均電子密度曲線和平面平均靜電勢曲線，和文獻里的結果相當一致，讀者可以看看。

6 總結&其它

本文介紹了（局部）積分曲線、平面平均密度差曲線、電荷位移曲線這些概念，并詳解介紹了如何將CP2K第一性原理程序與強大的Multiwfn波函數分析程序相結合繪制密度差圖以及這些曲線。從本文的例子可見，密度差等值面圖能很直觀地展現三維空間中各處電子密度的凈增減，而對于恰當的情況（電子轉移和極化以單一方向為主），若將之與平面平均密度差曲線、電荷位移曲線相結合，則能考察得更清楚、更便于定量討論，還便于準確地橫向對比。例如一個固體表面吸附不同分子，想考察電子轉移特征的差異，就可以把不同情況的平面平均密度差曲線同時繪制在一張圖上，差異一目了然。

本文的做法絕不僅限于CP2K用戶使用。記錄電子密度的cub文件用任何其它程序產生也都可以，cub是化學領域最常見的記錄格點數據的格式，因此被支持得很廣泛。Multiwfn支持的格點數據文件格式有很多，還包括dx、grd、vti、VASP的CHGCAR/CHG/ELFCAR/LOCPOT。對于VASP用戶，想得到密度差，可以先用Multiwfn載入CHGCAR，進主功能13用子功能0轉成cub文件，對整體和各個片段都這么做得到各自的電子密度cub文件后，再按上文方式進行操作。

Multiwfn也可以基于波函數文件直接用主功能5計算密度差格點數據。如果你是Gaussian、ORCA等主流量子化學程序用戶，一般都建議先計算出各個狀態的波函數文件（支持的格式和產生方法見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》http://www.shanxitv.org/379），之后用主功能5的自定義運算功能可輕松得到密度差格點數據（參考《使用Multiwfn作電子密度差圖》http://www.shanxitv.org/113里使用主功能5的過程），然后直接就能去主功能13里的子功能18里繪制密度差局部積分曲線、平面平均密度差和電荷位移曲線。對于CP2K也可以如此，先對涉及的各個狀態產生記錄波函數信息的molden文件，按照《使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用》（http://www.shanxitv.org/588）里的說明自行把晶胞信息作為[Cell]字段寫進去，然后用主功能5做自定義運算就可以獲得密度差格點數據。

在《全面探究18碳環獨特的分子間相互作用與pi-pi堆積特征》（http://www.shanxitv.org/572）和《一篇文章深入揭示外電場對18碳環的超強調控作用》（http://www.shanxitv.org/570）介紹的我關于18碳環體系的研究中都使用了本文提到的功能繪制了密度差的局部積分曲線，由其中的例子可見這對于考察分子間相互作用和外電場對電子分布的影響特別有好處，非常直觀。

本文的做法是非常普適的，絕不僅限于展現片段間相互作用的密度差。對于電子激發導致激發態相對于基態電子密度分布改變的密度差、外加電場導致電子分布變化的密度差、體系電離或附著電子導致的密度變化的密度差等等，都可以按本文的方法繪制相應的（局部）積分曲線來研究。

使用本文的做法通過Multiwfn產生相關數據若用于發表，記得按照程序啟動時的提示恰當引用Multiwfn。

使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和分析電子激發態

sobereva — Mon, 28 Feb 2022 22:24:00 +0800

使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和分析電子激發態

文/Sobereva@北京科音

First release: 2022-Feb-28 Last update: 2023-Feb-10

0 前言

波函數分析程序Multiwfn（http://www.shanxitv.org/multiwfn）有十分豐富的電子激發分析功能，已被量子化學研究者廣為使用，在《Multiwfn支持的電子激發分析方法一覽》（http://www.shanxitv.org/437）有全面的功能介紹。而且Multiwfn還有靈活強大的繪制電子光譜的功能，見《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）。以往都是Gaussian、GAMESS-US、ORCA等量子化學程序用戶在研究電子激發相關問題時使用Multiwfn的這些功能，而近期，Multiwfn的一些相關功能也已支持了CP2K第一性原理程序。免費、高效的CP2K算周期性大體系的TDDFT相當快。CP2K用戶使用TDDFT計算分子或周期性體系的激發態后，可以非常方便地使用Multiwfn繪制UV-Vis譜和分析電子激發特征。而且Multiwfn具有非常方便的創建CP2K的TDDFT輸入文件的功能，上手計算十分容易。本文的目的是通過一個典型的分子晶體例子，令讀者了解使用Multiwfn和CP2K通過TDDFT研究周期性體系電子激發問題的基本做法，不少關鍵性的細節也一起說明，以便讀者能舉一反三，解決自己研究中的問題。

讀者務必使用2023-Feb-10及以后更新的Multiwfn版本，否則和本文介紹的情況會有所不同。Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，不了解Multiwfn者建議看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167），使用Multiwfn發表文章別忘了按照程序啟動時的提示進行恰當引用。本文用的CP2K是9.1版，安裝方法見《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。如果讀者對TDDFT計算完全缺乏了解，建議閱讀《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）和《亂談激發態的計算方法》（http://www.shanxitv.org/265）了解相關基本知識。

本文涉及的相關文件都可以在http://www.shanxitv.org/attach/634/file.rar里得到。文件約20MB，如果下載慢請換個時段下。

本文對CP2K的電子激發計算部分僅僅做非常簡要的描述，在筆者講授的北京科音CP2K第一性原理計算培訓班有專門的一節非常詳細講授CP2K做電子激發計算（也包括X光吸收），介紹見http://www.keinsci.com/workshop/KFP_content.html，非常歡迎了解和參加。

1 體系特征和一些計算細節

本文的例子是一個Donor-π-Acceptor型分子4-(N,N-dimethyl-amino)-4′-(2,3,5,6-tetrafluorostyryl)-stilbene構成的晶體，如下所示，晶胞里有4個分子，共192個原子。我們將要繪制它的UV-Vis譜，并對其中對應強吸收峰的電子激發分析其本質。此體系X光衍射得到的晶體結構是本文文件包里的bulk.cif。

這個體系實際上是ACS Omega, 3, 10481 (2018)一文中理論研究的體系，文中深入研究了這個分子晶體的J聚集導致吸收峰相對于單體分子紅移的本質。作者用的是Quantum ESPRESSO做的TDDFT計算，而本例我們將用CP2K做TDDFT。

此晶體的原胞的三個邊長分別是7.59、5.87、41.51埃。由于前兩個晶胞尺寸較短，而CP2K的TDDFT計算時不能考慮k點，所以如果想更準確計算吸收光譜的話建議將體系復制成2*3*1的超胞。這是常見做法，J. Chem. Theory Comput., 17, 5214 (2021)在TDDFT計算分子晶體光譜時還特意對超胞尺寸做了結果的收斂性測試。搞成超胞自然會導致耗時劇增，由于本文只是示例目的，為了節約耗時，就只用原胞來算了。ACS Omega那篇文章里也是直接用的原胞來計算和分析。

由于X光衍射得到的晶體結構中的氫的位置一般都是不準的，《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）里也說過這個問題，因此電子激發計算前應當至少先對氫做優化，不過本文暫且忽略這個步驟，而自己實際研究中應當注意這個問題。

CP2K的TDDFT只支持TDA近似形式，也可以叫TDA-DFT，計算出的振子強度不如嚴格的TDDFT準確，但激發能的準確度并沒打折扣。

ACS Omega原文里是基于PBE泛函做的TDDFT，我們也用這個泛函來做。眾所周知，雖然純泛函在第一性原理程序里計算比雜化泛函快得多得多，但有明顯低估激發能的傾向（尤其是對于電荷轉移激發），故讀者實際研究時如果對激發能的準確性有要求，屆時別忘了用PBE0等適合的雜化泛函來算，這在CP2K中也是支持的，而且能算不小的體系。

CP2K 9.1 CP2K做TDDFT只支持GPW，即必須用贗勢基組而不能用全電子基組（后記：從2022.1開始支持了GAPW的TDDFT，可以用全電子基組了）。對于TDDFT目的用TZVP-GTH贗勢基組就不錯，精度足夠，而且也不很貴。再好一點還可以用TZV2P-GTH。用CP2K里常用的DZVP-MOLOPT-SR-GTH贗勢基組算TDDFT時還能比TZVP-GTH更便宜一些，但是由于此基組收縮度太高，導致之后用Multiwfn做很多分析的時候在計算格點數據時會非常昂貴，所以本文用TZVP-GTH。

2 使用CP2K做TDDFT計算

用Multiwfn創建CP2K輸入文件非常方便，在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》
（http://www.shanxitv.org/587）里有詳細說明。

啟動Multiwfn，然后輸入
bulk.cif //在本文文件包里提供了
cp2k //產生CP2K輸入文件
bulk_TDDFT.inp //將要生成的CP2K輸入文件的文件名
15 //開啟激發態計算
y //讓計算任務順帶輸出記錄分子軌道的molden文件
30 //把最低30個空軌道算出來并存到molden文件里
16 //設置激發態計算數目
50 //算最低50個激發態
2 //選擇基組
-3 //TZVP-GTH
0 //生成輸入文件

現在當前目錄下就有了bulk_TDDFT.inp。輸入文件默認的CUTOFF等設置對于當前計算來說都是適合的，所以不用做額外改動。直接用CP2K運行此文件，輸出文件是本文文件包里的bulk_TDDFT.out。任務同時產生了molden文件bulk_TDDFT-MOS-1_0.molden。

如《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）一文所述，CP2K直接產生的molden文件里不含晶胞信息，沒法用于Multiwfn做分析，因此必須手動把晶胞信息填進去，做法在此文也明確說了。現在把以下內容插入到molden文件的第二行，本文文件包里的molden文件是我已經改好的。這部分內容大家直接從CP2K輸入文件里的&CELL部分復制過去即可。

[Cell]
7.59400000     0.00000000     0.00000000
0.00000000     5.87400000     0.00000000
-3.11381184     0.00000000    41.39304623

下面說一些和當前的TDDFT計算有關的細節。

輸入文件里的&PROPERTIES - &TDDFPT是TDDFT計算的相關設置，含義在Multiwfn自動產生的注釋里都寫明了。其中的MIN_AMPLITUDE 0.01代表組態系數的絕對值大于0.01的才會輸出，這個數設得越小，之后Multiwfn做空穴-電子分析、躍遷密度分析、NTO分析等會越精確。要求較高的話建議改小到0.001，但由于輸出的組態會更多，會令輸出文件信息量更大，而且在Multiwfn的激發態分析過程中耗時也會高一些。如果只需要大概看個分析結果，用Multiwfn默認設的0.01就夠了，至少結果是定性正確的。

之所以此任務要求產生記錄了分子軌道的molden文件，是因為之后使用Multiwfn做許多電子激發分析都需要利用到分子軌道信息。這里只要求記錄最低30個空軌道而非所有空軌道，并非是為了節約耗時（相對于激發態計算過程，算出所有空軌道并記錄也不會增加多少耗時），而是因為記錄所有空軌道的話molden文件可能非常大，對大體系甚至可能達到幾GB，又占硬盤拷貝又慢。TDDFT計算中，每個激發態波函數都是由一大批組態函數線性組合來描述的，不同組態函數形式上對應于不同形式的軌道躍遷，每個組態函數對激發態波函數都各有各的貢獻量，參看《電子激發任務中軌道躍遷貢獻的計算》（http://www.shanxitv.org/230）。我們當前算出來的激發態的組態函數的系數在CP2K輸出文件里可以直接看到，如下所示。對于能量比較低的一批激發態，如果在電子激發分析時只考慮貢獻不可完全忽略的組態函數（如組態系數絕對值大于0.01的），則這些組態函數只可能涉及到占據軌道向比較低階空軌道的躍遷。對本文的例子，算最低30個空軌道基本足夠了，如果想絕對穩妥的話建議算50個空軌道。一般來說，對越高階激發態做空穴-電子分析，建議算越多的空軌道，因為越高階激發態中涉及越高空軌道的組態函數越可能無法被忽略。后文還會再專門提到怎么準確檢驗算的空軌道數目是否確實夠。

-------------------------------------------------------------------------------
-                            Excitation analysis                              -
-------------------------------------------------------------------------------
        State             Occupied              Virtual             Excitation
        number             orbital              orbital             amplitude
-------------------------------------------------------------------------------
             1   1.64779 eV
                               296                  297              -0.902183
                               294                  298               0.306244
                               295                  299               0.275238
                               293                  300              -0.124557
                               296                  302              -0.013313
                               294                  304               0.013288
                               293                  303               0.013231
                               295                  301              -0.012515
...略
            50   2.97063 eV
                               295                  307              -0.629538
                               294                  306              -0.596348
                               293                  305              -0.350078
                               296                  308              -0.348347
                               296                  310              -0.036322
                               294                  309               0.019915
                               294                  303              -0.016826
                               293                  304              -0.016232
                               296                  301               0.015761
                               295                  302               0.014885
                               284                  297               0.013140
                               283                  298               0.012456

算的激發態數目越多，計算耗時越高。當前的例子只計算了50個激發態，對于獲得較完整的UV-Vis光譜來說并不夠。從輸出文件里以下部分可見，第50激發態的激發能是2.97 eV，折合417 nm，連可見光區都沒覆蓋完整。但考慮到純泛函算的激發能整體明顯偏低，如果用諸如PBE0算，可能最高能達到比如3.6 eV （344 nm）左右，屆時可見光部分就都覆蓋了。如果要把近紫外區也都算出來，需要算的態數建議翻倍。

         State    Excitation        Transition dipole (a.u.)        Oscillator
         number   energy (eV)       x           y           z     strength (a.u.)
         ------------------------------------------------------------------------
TDDFPT|      1       1.64779 -3.4173E-08 4.5488E-09 -1.0672E-08   5.25777E-17
TDDFPT|      2       1.66084 -1.4486E-08 -8.1846E-08 9.6347E-09   2.84886E-16
TDDFPT|      3       1.66146   8.5244E-08 3.6451E-02 5.0929E-08   5.40824E-05
...略
TDDFPT|     49       2.96978   2.8362E-01 -9.3226E-07 2.7092E-01   1.11931E-02
TDDFPT|     50       2.97063 -6.2827E-06 -1.6851E-01 -3.5284E-05   2.06651E-03

從以上信息中也可以看到振子強度、躍遷偶極矩這些信息，和Gaussian等量子化學程序做TDDFT給出的信息形式非常類似。

3 繪制電子光譜

用Multiwfn繪制電子光譜功能的全面介紹、各種細節以及相關的基本原理請讀者看《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224），這里僅畫個很簡單的圖。多數UV-Vis譜橫坐標用的是nm，但在前述的ACS Omega文章里研究這個晶體的光譜時橫坐標用的是eV，如下所示。其中黑線是和我們計算模型一致的bulk狀態的吸收曲線，我們試圖重現一下。圖中M、D、C分別是作者自行基于晶體結構構造的分子單體、二聚體、鏈狀排列時的吸收曲線，這不是我們本例關心的。

啟動Multiwfn然后載入bulk_TDDFT.out，之后輸入
11 //繪制光譜
3 //UV-Vis
10 //修改橫坐標單位
1 //eV
8 //設置峰展寬用的半高全寬（FWHM）
0.07 eV //和ACS Omega文章里用的一致。順帶一提，Multiwfn默認用的FWHM比這個大得多，一般適合溶液的吸收光譜。算晶體的話可以設小一些
16 //設置顯示光譜極大極小點的標簽
1 //設置標簽狀態
1 //把曲線的極大點數值顯示出來
3 //設置標簽上的小數位數
2 //兩位小數
0 //返回
0 //作圖

現在看到的圖如下所示（注：筆者為了光譜曲線更平滑，此處先選擇-4把導出格式設為了pdf，然后選了1導出pdf文件，下圖實際上是pdf的截圖，這比直接在屏幕上顯示的光滑得多，且可以無損縮放）。讀者也可以自行再用Multiwfn界面上的選項修改橫、縱坐標范圍和標簽間隔，使得坐標軸的標簽更整齊。

由于我們算了50個激發態，比文獻里算的更多，因此上面的光譜差不多把3 eV激發能以內的電子激發對應的吸收都算出來了。上圖中紅色曲線是吸收光譜，對應左邊坐標軸，藍色標簽標注的是吸收峰的位置。由圖可見，在1.71 eV處有個很弱的吸收峰（振子強度0.04），而在1.92 eV處的吸收巨強，振子強度超過4（具體是4.13）。這和ACS Omega文章中的譜圖的特征非常吻合。那篇文章補充材料中的表S1里，正好在1.71 eV處也有個激發，振子強度僅為0.08，而在1.887 eV處有個振子強度7.96的賊強的激發。我們算的電子激發能和文獻里相符極好，但振子強度比文獻里的小了將近一倍，這主要是由于CP2K用了TDA近似所致。不過由于不同的峰的相對高度受TDA的影響遠小于絕對高度受其影響，而我們一般只關心曲線形狀，所以基于CP2K的TDDFT算的光譜是可以用的。

想搞清楚各個吸收峰對應哪個電子激發很簡單。在Multiwfn當前的繪制光譜的界面里選-1，就能把激發態序號、不同單位的激發能以及振子強度都直接顯示出來，如下所示。上面的譜圖里黑色豎線橫坐標位置是激發能，豎線高度對應右邊坐標軸，是振子強度，和以下信息一對照便知各個峰對應什么激發態。顯然，1.71 eV的極弱吸收對應的是第8激發態（S8），1.92 eV的極強吸收對應的是第13激發態（S13）。接下來我們要著重對S13這個重要的電子激發做一些分析，看看它對應的本質為何。

Index Excit.energy(eV       nm        1000 cm^-1)       Oscil.str.
   1         1.64779      752.42719       13.29032        0.00000
   2         1.66084      746.51502       13.39558        0.00000
   3         1.66146      746.23644       13.40058        0.00005
   4         1.67536      740.04512       13.51269        0.00471
   5         1.70466      727.32510       13.74901        0.00153
   6         1.70527      727.06492       13.75393        0.00000
   7         1.70540      727.00950       13.75498        0.00000
   8         1.70974      725.16406       13.78998        0.03991
   9         1.73479      714.69284       13.99203        0.00786
10         1.74660      709.86030       14.08728        0.00000
11         1.75022      708.39209       14.11648        0.00025
12         1.75461      706.61970       14.15189        0.00000
13         1.92496      644.08715       15.52585        4.13484
14         2.15741      574.69002       17.40069        0.00000
...略

4 考察主要的分子軌道躍遷

特別常見的分析電子激發的方式是看激發中哪些分子軌道的躍遷對激發有主要貢獻，然后去看軌道圖形考察。計算貢獻的方法在《電子激發任務中軌道躍遷貢獻的計算》（http://www.shanxitv.org/230）中專門說了，并且如《使用Multiwfn便利地查看所有激發態中的主要軌道躍遷貢獻》（http://www.shanxitv.org/529）所述，利用Multiwfn還可以極其便利地一鍵得到所有電子激發中的各種軌道躍遷貢獻，這里就對當前體系用一下。

啟動Multiwfn，載入bulk_TDDFT.out，之后輸入
18 //電子激發分析
15 //顯示所有激發態中的分子軌道躍遷情況
馬上看到以下內容

#   1   1.6478 eV    752.43 nm   f= 0.00000   Spin multiplicity= 1:
   H -> L 81.4%, H-2 -> L+1 9.4%, H-1 -> L+2 7.6%
#   2   1.6608 eV    746.52 nm   f= 0.00000   Spin multiplicity= 1:
   H -> L+1 75.2%, H-2 -> L 21.1%
...略
# 12   1.7546 eV    706.62 nm   f= 0.00000   Spin multiplicity= 1:
   H-3 -> L+3 61.6%, H-1 -> L+2 21.2%, H-2 -> L+1 14.5%
# 13   1.9250 eV    644.09 nm   f= 4.13484   Spin multiplicity= 1:
   H-2 -> L+2 30.3%, H -> L+3 25.6%, H-3 -> L 20.2%, H-1 -> L+1 17.6%
...略

可見基態S0到S13的激發同時由多個軌道躍遷所明顯貢獻，HOMO-2 -> LUMO+2貢獻了30.3%、HOMO -> LUMO+3貢獻了25.6%、HOMO-3 -> LUMO貢獻了20.2%、HOMO-1 -> LUMO+1貢獻了17.6%。像這種情況就極度不適合通過觀看分子軌道來分析，因為得同時看一大堆軌道，極其麻煩。這種情況最理想的分析方法就是用Multiwfn獨家的空穴-電子（hole-electron）分析，對于任何激發都能轉化為“空穴”到“電子”的躍遷，分析時只需要看這兩個函數的分布即可，在下一節我們將做這種分析。

由以上信息可見S0到S1的激發倒是有主導性的分子軌道躍遷，即HOMO -> LUMO貢獻達到81.4 %，此時看這一對軌道的特征就能大致判斷電子激發情況。這里不妨就看一下HOMO和LUMO的分布是什么樣的。關于Multiwfn顯示分子軌道在《使用Multiwfn觀看分子軌道》（http://www.shanxitv.org/269）有詳細介紹。

啟動Multiwfn，載入bulk_TDDFT-MOS-1_0.molden。注意屏幕上顯示了如下的晶胞信息：
Cell angles: Alpha= 90.0000 Beta= 94.3020 Gamma= 90.0000 degree
可見當前晶胞不是正交的，這種情況下Multiwfn在計算格點數據后直接顯示的等值面不完全正確，需要用Multiwfn導出格點數據后通過VMD或VESTA來顯示等值面。不過，由于此例的Beta角偏離90度也不太大，因此直接用Multiwfn觀看等值面也勉強可以。

進入主功能0，此時文本窗口顯示了HOMO-LUMO gap等軌道相關信息。在圖形界面里點擊右側Show Labels關閉原子標簽、點擊Show axis關閉坐標軸。選擇菜單欄Other settings - Toggle showing cell frame關閉盒子邊框。選擇菜單欄Isosur. quality - Good quality。然后在界面右下角的文本框里輸入h然后按回車，代表顯示HOMO。之后把右下角的Isovalue（等值面數值）改成一個合適的值0.026。此時看到的HOMO等值面圖如下圖左側所示。然后在右下角列表里點擊297切換到LUMO軌道，圖像如下圖右側所示。

由對比可見，HOMO和LUMO主體分布于不同分子上，因此S0-S1是分子間電荷轉移激發。關于電子激發類型介紹看《圖解電子激發的分類》（http://www.shanxitv.org/284）。

下面我們讓Multiwfn導出HOMO和LUMO的格點數據，從而之后能放到專門的可視化程序里觀看得到更好的效果。點當前圖形窗口右上角的RETURN按鈕返回主菜單，然后輸入

200 //主功能200
3 //計算并導出軌道波函數的格點文件
h //HOMO
9 //根據晶胞設置格點數據的盒子原點、邊長和格點間距
[回車] //用(0,0,0)作為原點
[回車] //格點數據的盒子的三個矢量和當前晶胞相同
[回車] //用默認的0.25 Bohr格點間距（為了降低耗時也可以設大一些，如0.4）

然后當前目錄下就出現了h.cub。將之拖入VESTA，直接就看到了下圖，效果不錯。對LUMO也這么干。

5 空穴-電子分析

下面做空穴-電子分析。此方法的原理、細節、大量用于分子體系的實例在《使用Multiwfn做空穴-電子分析全面考察電子激發特征》（http://www.shanxitv.org/434）都給出了，沒讀過的話務必一讀，也請按照此文所述恰當引用相關文章。

啟動Multiwfn，載入bulk_TDDFT-MOS-1_0.molden，然后輸入
18 //電子激發分析
1 //空穴-電子分析
bulk_TDDFT.out //CP2K輸出文件
13 //分析第13激發態(S13)，也就是那個吸收極強的激發態

注意此時屏幕上會看到提示Involved highest unoccupied orbital: 324 (HOMO+ 28)，即這個激發態涉及到的組態函數里（即由于系數絕對值>0.001因而被輸出到CP2K輸出文件里的）最高牽扯到的空軌道是第28個空軌道。前面說了，我們當前的TDDFT計算讓CP2K把最低30個空軌道信息寫入了molden文件，比28更大，所以記錄的空軌道數是夠多的。如果記錄的數目小于這里的28，則對這個激發態做的空穴-電子分析可能不準確。每次做空穴-電子分析的時候建議留意一下這里的信息，如果發現存入molden文件里的空軌道數不夠多，應該把CP2K輸入文件里ADDED_MOS后面的值改大，然后重算一次單點任務得到新的molden文件，基于它和之前的TDDFT任務的輸出文件做空穴-電子分析。

接著輸入
1 //可視化空穴、電子、躍遷密度等函數
9 //根據晶胞設置格點數據的盒子原點、邊長和格點間距
[回車] //用(0,0,0)作為原點
[回車] //格點數據的盒子的三個矢量和當前晶胞相同
[回車] //用默認的0.25 Bohr格點間距（為了降低耗時也可以設大一些，如0.4）

在筆者的8核i7-10870H筆記本上兩分鐘算完，然后出現了后處理菜單。雖然用此處的選項在Multiwfn里也能直接顯示空穴、電子等函數的等值面，但由于盒子不是正交的，所以還是先導出格點數據再用VMD、VESTA等程序來顯示。在空穴-電子分析框架里有很多定量指數和函數可以考察電子激發特征，這一節我們就考察空穴、電子以及密度差（charge density difference，當前語境下對應于電子減去空穴），因此在后處理菜單中依次選擇選項10、11、15導出它們的格點數據分別成為hole.cub、electron.cub、CDD.cub。它們在本文的文件包里也提供了。

把這三個cub文件分別拖到VESTA里顯示。等值面數值用0.001（Properties里選Isosurfaces標簽頁，點擊當前僅有的項，把Isosurface level設為0.001），在Objects - Volumetric Data里取消選擇Show Sections，即不顯示截面，然后得到的圖如下所示。

由圖可見，S0-S13激發也是個典型的電荷轉移激發，因為空穴和電子分布在明顯不同的分子上。從CDD圖上考察甚至更方便，只需一張圖就展現出了電子是怎么轉移的，青色和黃色分別是電子激發時電子減少和增加的地方。

如果用VMD來顯示等值面，還可以直接在顯示的時候把周期鏡像顯示出來，由于此時能看到完整分子，明顯更加容易觀察，如下所示

這里具體說一下用VMD怎么作上面的圖。在《在VMD里將cube文件瞬間繪制成效果極佳的等值面圖的方法》（http://www.shanxitv.org/483）中專門給了一個腳本用來十分快捷地顯示效果理想的等值面圖，先使用這個腳本把CDD.cub的0.001等值面顯示出來。然后在Graphics - Representation里把顯示分子結構的那個Representation（以下簡寫為Rep）刪掉，對于分別顯示CDD的正值和負值等值面的那兩個Rep，在Periodic標簽頁里都把-X和+Y都選上，這樣等值面就在X負方向和Y正方向各復制了一次，看上去像超胞的情況。之所以這里沒有把當前分子結構也這么操作來顯示相鄰的鏡像，是因為這么做之后跨盒子的鍵是顯示不了的，看起來略別扭。為了得到與當前等值面相對應的超胞結構文件，下面使用《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）文中介紹的做法。在Multiwfn載入bulk_TDDFT-MOS-1_0.molden后輸入
300 //主功能300
7 //對結構進行操作
19 //對體系進行平移復制
-1 //在第一個方向負方向復制一次
2 //在第二個方向正方向復制為原先兩倍
1 //在第三個方向保持不變
-2 //把當前結構導出為pdb文件
supercell.pdb
現在當前目錄下就有了supercell.pdb。將之拖到已顯示出CDD等值面的VMD圖形窗口里，把其顯示方式設為CPK，即得到上圖。

6 躍遷密度分析

前述的ACS Omega文章里通過所謂的induced charge density分布考察了分子晶體中不同分子在激發時的激子耦合對激發能的影響。實際上induced charge density就是《使用Multiwfn繪制躍遷密度矩陣和電荷轉移矩陣考察電子激發特征》（http://www.shanxitv.org/436）中筆者介紹的躍遷密度（transition density），文中的這類圖都可以用Multiwfn精確重復出來。這里我們也畫一下。

之前Multiwfn做完空穴-電子分析進入后處理菜單后，可以看到選項13是用來導出躍遷密度的cub文件的，選擇它，在當前目錄下就出現了transdens.cub，在本文的文件包里也提供了。按照上一節的方式對它繪制等值面圖，等值面數值為±0.0005的情況如下所示

讀者若仔細將上圖和ACS Omega文中圖2里的體相的induced charge density（下圖）相對比，會發現是完全一樣的，說明了計算的合理性。

至于怎么討論激子耦合本文就不詳述了，讀者可自行參看ACS Omega文中的討論，以及Multiwfn手冊4.A.9節中的計算激子耦合的介紹。

7 NTO分析

NTO（自然躍遷軌道）分析是很流行的激發態分析方法，在《使用Multiwfn做自然躍遷軌道(NTO)分析》（http://www.shanxitv.org/377）里有專門介紹和例子，本文就不再累述了。這個分析在Multiwfn中現在也可以用于周期性體系，這里我們就做一下。值得一提的是，CP2K自身也有NTO分析功能，但筆者發現很難用，一次只能對一個激發態輸出NTO軌道波函數，而且產生的非占據NTO軌道還明顯異常，因此很不推薦用CP2K自己的NTO功能，至少對目前的版本來說。

啟動Multiwfn，載入bulk_TDDFT-MOS-1_0.molden，然后輸入
18 //電子激發分析
6 //NTO分析
bulk_TDDFT.out //CP2K輸出文件
13 //考察S0-S13激發

立刻就在屏幕上給出了本征值最大的10個NTO對的本征值：

The highest 10 eigenvalues of NTO pairs:
0.311175 0.270759 0.215496 0.195281 0.000291
0.000248 0.000180 0.000177 0.000000 0.000000
Sum of all eigenvalues: 0.993608

由數據可見，NTO方法對當前這個電子激發表現極差！在本文第4節我們看到，對這個激發貢獻最大的四對MO躍遷貢獻分別是30.3%、25.6%、20.2%、17.6%，而在變換成NTO躍遷描述之后，依然還是有四對貢獻很大的躍遷，分別是31.1%、27.1%、21.5%、19.5%，比起原先的描述根本就沒什么改進。所以NTO分析在很多情況下是失敗、無用的，而空穴-電子分析則比NTO普適得多，對于任何情況都能很理想地使用。所以若無特殊情況（比如非要得到軌道相位信息等），不建議用NTO而應當用空穴-電子分析。雖然NTO對于分子體系多半情況表現得還行，但我發現對于周期性體系大多時候NTO派不上什么用處，即沒法將電子激發充分簡化描述成一對NTO躍遷所主導的情況。

雖然此例NTO分析沒用，但作為演示，還是看看NTO軌道。在Multiwfn里選擇3，再輸入S13_NTO.mwfn，在當前目錄下就出現了此文件。然后重新啟動Multiwfn，載入S13_NTO.mwfn，進入主功能0，然后在圖形窗口頂端的菜單中選擇Orbital info. - Show up to LUMO+10，此時在文本窗口可看到

...略
Orb:   292 Ene(au/eV):     0.000291       0.0079 Occ: 2.000000 Type:A+B (?   )
Orb:   293 Ene(au/eV):     0.195281       5.3139 Occ: 2.000000 Type:A+B (?   )
Orb:   294 Ene(au/eV):     0.215496       5.8639 Occ: 2.000000 Type:A+B (?   )
Orb:   295 Ene(au/eV):     0.270759       7.3677 Occ: 2.000000 Type:A+B (?   )
Orb:   296 Ene(au/eV):     0.311175       8.4675 Occ: 2.000000 Type:A+B (?   )
Orb:   297 Ene(au/eV):     0.311175       8.4675 Occ: 0.000000 Type:A+B (?   )
Orb:   298 Ene(au/eV):     0.270759       7.3677 Occ: 0.000000 Type:A+B (?   )
Orb:   299 Ene(au/eV):     0.215496       5.8639 Occ: 0.000000 Type:A+B (?   )
Orb:   300 Ene(au/eV):     0.195281       5.3139 Occ: 0.000000 Type:A+B (?   )
Orb:   301 Ene(au/eV):     0.000291       0.0079 Occ: 0.000000 Type:A+B (?   )
...略

這里以a.u.為單位的能量對應的是NTO本征值。可見占據和非占據NTO是配對的。比如296和297號軌道構成一個對電子激發貢獻31.1%的NTO對，其中296和297分別是占據和非占據的NTO。軌道圖形請讀者效仿前面第四節自行觀看，筆者就不多說了。

8 總結

本文介紹了如何將Multiwfn與CP2K相結合，繪制電子光譜，并對電子激發特征進行分析。可見二者結合使用，使得研究周期性體系的電子激發相關問題相當容易和便利，比起常見的Gaussian+Multiwfn的組合研究分子體系的電子激發并沒有額外的復雜性。希望讀者在本文的基礎上舉一反三，使用CP2K+Multiwfn研究更多類型的周期性體系的電子激發。

Multiwfn的電子激發分析分析功能非常多，遠遠不限于以上所述的幾種。但這些方法中目前只有部分明確支持周期性體系，有些分析在原理或程序上不支持周期性，而有些可能能支持但尚未做測試。明確支持周期性的分析在Multiwfn手冊2.9.2.2節說了。以后Multiwfn在周期性分析方面還會不斷做擴展和完善，使得更多方法完美兼容周期性體系。如有疑問歡迎在程序啟動時提示的Multiwfn論壇上發帖咨詢。

使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用

sobereva — Sat, 27 Feb 2021 14:14:00 +0800

使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用

文/Sobereva@北京科音

First release: 2021-Feb-27 Last update: 2022-Mar-8

1 前言

2010年提出的Noncovalent interaction (NCI)分析，也即Reduced density gradient (RDG)分析，是如今使用極為普遍的圖形化考察化學體系中弱相互作用的方法，筆者寫過諸多相關文章介紹原理和實現，還錄了操作演示視頻，如下所示，故本文對基本原理不再累述，并假定讀者已經具備了背景知識
? 《使用Multiwfn圖形化研究弱相互作用》（http://www.shanxitv.org/68）
? 《用Multiwfn+VMD做RDG分析時的一些要點和常見問題》（http://www.shanxitv.org/291）
? 《通過鍵級曲線和ELF/LOL/RDG等值面動畫研究化學反應過程》（http://www.shanxitv.org/200）
? 《使用Multiwfn做NCI分析展現分子內和分子間弱相互作用》（演示視頻。https://www.bilibili.com/video/av71561024）
? Multiwfn手冊3.23.1節
? 《使用Multiwfn研究分子動力學中的弱相互作用》（http://www.shanxitv.org/186）
筆者開發的Multiwfn是做NCI分析最流行、最好用的程序。但以前的版本主要面向分子體系，基于Gaussian、ORCA等量子化學程序產生的波函數進行分析。對于周期性體系往往也可以通過團簇模型結合Multiwfn做NCI分析，比如《使用量子化學程序基于簇模型計算金屬表面吸附問題》（http://www.shanxitv.org/540）里對銀團簇表面吸附苯做了NCI分析、在Multiwfn原文J. Comput. Chem., 33, 580 (2012)里截取了尿素團簇做了NCI分析，但是還是有很多情況不便于或者幾乎不適合用簇模型描述。為了讓周期性體系的NCI分析可以非常容易地實現、使這個方法能考察更廣泛的體系，從2021年2月更新的Multiwfn開始，Multiwfn也支持了結合CP2K程序產生的波函數做NCI分析，這使得圖形化考察固體內部和表面的弱相互作用極其容易，將在本文進行詳細示例。看了第3節的例子后你會發現計算過程真是巨簡單。本文也會提及很多讀者應了解的注意事項，以便在實際研究中遇到問題時恰當變通。

另外，Independent gradient model (IGM)也是非常重要的圖形化研究弱相互作用的方法，相對于NCI的關鍵好處是可以很靈活地自定義片段，能只顯示片段間的弱相互作用而不會受到片段內的干擾。此外，Multiwfn做IGM分析時還可以給出諸多定量指標便于討論。Multiwfn做IGM分析在此文有詳述：《通過獨立梯度模型(IGM)考察分子間弱相互作用》（http://www.shanxitv.org/407）。如今Multiwfn也將IGM分析擴展到了周期性體系，將在下文示例。本文還會順帶提及筆者提出的IGM方法的改良，即IGM based on Hirshfeld partition (IGMH)，比IGM圖像效果明顯好得多，詳見《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621）。

值得一提的是，Multiwfn還支持筆者提出的IRI方法，可以同時把弱相互作用和化學鍵作用區域直觀地展現，很有實用價值，在《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）中有專門的介紹，并且給了計算孤立體系和計算周期性體系的例子，強烈建議一看！

對一般情況，筆者強烈建議用后來提出的IRI代替NCI，能明顯展現更多有價值的信息；也強烈建議用后來提出的IGMH代替IGM，圖像效果好得多，且物理意義更嚴格。IRI和NCI的分析操作幾乎完全相同，而IGMH和IGM的分析操作幾乎完全相同。在《一篇最全面介紹各種弱相互作用可視化分析方法的文章已發表！》（http://www.shanxitv.org/667）中介紹的筆者的綜述文章中對上述可視化分析方法有非常全面深入的介紹，里面也展示了用于周期性體系的例子，推薦閱讀和引用。

Multiwfn可以在官網http://www.shanxitv.org/multiwfn免費下載，不熟悉Multiwfn的話注意參看《Multiwfn FAQ》（http://www.shanxitv.org/452）和《Multiwfn入門tips》（http://www.shanxitv.org/167）。請讀者務必使用2021-Feb-25日及之后更新的Multiwfn以確保情況和下文內容相符。本文例子中周期性體系的DFT計算部分使用CP2K 8.1版實現，安裝方法見《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。CP2K結合Multiwfn可以非常容易地使用，見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）。本文繪圖使用VMD 1.9.3版，可以在http://www.ks.uiuc.edu/Research/vmd/免費下載。

下文涉及的所有文件都可以在http://www.shanxitv.org/attach/588/file.rar下載。cub文件比較大所以就沒有提供了。

2 輸入文件

在給出具體例子之前先介紹一下Multiwfn做NCI、IGM分析都需要用什么輸入文件。

后記：關于這部分知識，更多、更詳細的介紹見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）。

2.1 普通NCI分析用的輸入文件

普通形式的NCI分析需要基于電子波函數計算電子密度及其導數。CP2K產生的molden文件記錄了Gaussian型基函數描述的電子波函數，Multiwfn可以基于它做各種波函數分析，也包括NCI分析。CP2K產生molden文件的方法在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）的3.1節已經示例過。如果你習慣自己寫CP2K輸入文件的話，是在&DFT字段里加上以下內容
&PRINT
&MO_MOLDEN
NDIGITS 8
GTO_KIND SPHERICAL
&END MO_MOLDEN
&END PRINT
這樣CP2K計算完畢后就會在當前目錄下產生一個.molden文件。

CP2K產生的molden文件還不能直接用于做周期性體系的NCI分析，因為標準的molden格式里沒有晶胞信息。為了能讓molden文件給Multiwfn提供晶胞信息，我給molden格式做了擴展，定義了[Cell]字段。用文本編輯器打開molden文件后，需手動在第二行插入[Cell]并在下面三行定義晶胞的三個平移矢量（埃），下面是個例子
[Molden Format]
[Cell]
7.13358000 0.00000000 0.00000000
0.00000000 7.13358000 0.00000000
0.00000000 0.00000000 7.13358000
[Atoms] AU
C 1 4 0.000000 0.000000 0.000000
C 2 4 1.685064 1.685064 1.685064
C 3 4 0.000000 3.370128 3.370128
...略
讓Multiwfn載入這樣修改過的molden文件后，晶胞信息就會被讀入，文件載入完畢后屏幕上就可以看到晶胞信息的提示了：
Cell information:
Translation vector 1   X=   13.480512 Y=    0.000000 Z=    0.000000 Bohr
Translation vector 2   X=    0.000000 Y=   13.480512 Z=    0.000000 Bohr
Translation vector 3   X=    0.000000 Y=    0.000000 Z=   13.480512 Bohr
Cell volume:       2449.7350 Bohr^3    (        363.0134 Angstrom^3 )

有一點需要注意的是，如今CP2K計算普遍使用MOLOPT基組，這是個完全廣義收縮的基組，而Multiwfn的計算代碼并未對廣義收縮的基組進行專門的優化，因此使用MOLOPT基組得到的波函數在Multiwfn里計算耗時會比較高。如果你有個像樣的服務器那無所謂，但如果只有個個人PC，算NCI的格點數據的時候對稍大的體系會比較吃力。我一般建議用比如6-31G*、6-311G**之類基本像樣的而且是片段收縮的全電子基組在CP2K里做個GAPW單點計算來得到NCI分析用的molden文件，比起用DZVP-MOLOPT-SR-GTH基組做GPW計算得到的molden文件在NCI分析耗時上能降低好幾倍。另外，用全電子基組比起用MOLOPT這樣的贗勢基組還有個好處是內核電子被明確表達了。NCI分析對基組不敏感，一般用6-31G*就足夠了，用更大基組也不會發現結果有什么可查覺的改變。

如果由于特殊情況，你就是想用諸如MOLOPT贗勢基組在Multiwfn里做NCI等波函數分析，那么應當再手動改一下molden文件，使得Multiwfn明確知道各個原子的價電子數是多少（指的是使用此贗勢時原子在孤立存在時的價電子數。和原子在當前體系里的實際價電子數沒任何關系），這有兩個必要性：
? molden文件里并沒有明確記錄原子的價電子數，這導致Multiwfn算原子電荷的時候不知道實際上當前原子的核電荷數是多少，因此給出的原子電荷是不合理的。
? 讓Multiwfn知道被贗勢贗化的電子數的話，Multiwfn會通過自帶的EDF內核電子密度數據庫提供內核電子密度，此時做純粹基于電子密度的各種分析時結果會較接近于全電子基組的情況，即更為真實。關于這點詳見《在贗勢下做波函數分析的一些說明》（http://www.shanxitv.org/156）。
讓Multiwfn知道各個原子價電子數是多少，可以通過以下兩種修改molden文件的做法實現
(1)把元素名后頭的元素序號改成價電子數。比如molden文件的[Atoms]字段某一行原本是
C 3 6 0.000000 0.000000 0.000000
如果贗勢贗化了碳的1s的兩個電子，即價電子數為4，則改為
C 3 4 0.000000 0.000000 0.000000
(2)更方便的做法是直接在[Atoms]字段之前加入[Nval]字段，這是我對molden格式的擴展。例如插入以下內容
[Nval]
C 4
O 6
就告訴Multiwfn所有的碳都是4個價電子，氧是6個價電子。加入了[Nval]和[Cell]字段的一個例子見Multiwfn自帶的examples\PBC\CP2K_diamond_2x2_DZVP-MOLOPT.inp。

注意CP2K在考慮k點的時候不能產生molden文件，Multiwfn本身也不支持考慮k點，所以對于小晶胞體系必須用足夠大的超胞來計算。

實際上Multiwfn還可以基于Gaussian程序的PBC計算產生的fch文件做周期性體系的波函數分析。但對于二維、三維周期性體系，Gaussian的這個功能幾乎什么都算不動，毫無實用性，這里就不提了。Multiwfn也可以基于含有晶胞平移矢量信息的mwfn格式（https://doi.org/10.26434/chemrxiv.11872524）做周期性體系的波函數分析，但CP2K目前版本產生不了。

2.2 準分子近似的NCI分析以及IGM分析用的輸入文件

還有一種NCI分析是基于準分子近似（promolecular approximation）實現的，即使用原子自由狀態電子密度簡單疊加產生的近似的分子電子密度（準分子密度，promolecular density）進行計算。對于周期性體系的準分子近似的NCI分析，只需要給Multiwfn提供含有原子坐標和晶胞信息的文件作為輸入文件即可，比如cif文件、Gaussian的PBC任務的gjf文件、含有CRYST1字段的pdb文件等，詳見《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）的第2節。

準分子近似的NCI計算耗時只是普通NCI分析的一個零頭。而且由于只需要結構信息，用任何第一性原理程序優化坐標/晶胞都可以，也可以直接用實驗的cif文件。但其結果合理性比普通的NCI略差，盡管通常來說定性一致。一般來說能算得動普通的NCI分析的時候就不建議用準分子近似。

IGM分析需要的輸入文件和準分子近似的NCI分析完全一樣，因為IGM分析也是完全依賴于幾何結構和Multiwfn內置的自由狀態原子密度實現的。

3 周期性體系的NCI分析示例

3.1 金剛石中的位阻作用

3.1.1 普通NCI分析

此例用NCI分析考察金剛石中的位阻作用。首先基于晶體結構，用Multiwfn產生一個做PBE/6-31G*單點計算并輸出molden文件的CP2K輸入文件，使用2*2*2的超胞。為此，用Multiwfn載入本文文件包里的金剛石的conventional cell的文件diamond.cif，然后輸入
cp2k //創建CP2K輸入文件
diamond_222.inp //要產生的CP2K輸入文件名
2 //選擇基組
10 //6-31G*
-2 //要求產生molden文件
-11 //幾何操作
19 //構建超胞
2 //第1個方向重復的次數
2 //第2個方向重復的次數
2 //第3個方向重復的次數
-10 //返回
0 //產生輸入文件

現在diamond_222.inp文件被產生在了當前目錄下，CP2K算完后就有了diamond_222-MOS-1_0.molden。用文本編輯器打開它，然后把diamond_222.inp里的晶胞信息復制進去，此時此文件開頭部分內容為
[Molden Format]
[Cell]
7.13358000     0.00000000     0.00000000
0.00000000     7.13358000     0.00000000
0.00000000     0.00000000     7.13358000
[Atoms] AU
C        1       6       0.000000       0.000000       0.000000
C        2       6       0.000000       3.370128       3.370128
...略

啟動Multiwfn，載入diamond_222-MOS-1_0.molden，然后輸入
20 //圖形化分析弱相互作用
1 //NCI分析
9 //借助晶胞信息定義格點
[按回車] //使用0,0,0作為盒子原點。這里說的盒子是指被計算的格點均勻分布的區域，后同
[按回車] //使用晶胞的三個平移矢量的長度作為盒子三個邊的長度
0.2 //格點間距

現在Multiwfn就開始計算NCI分析用的RDG和sign(λ2)rho這兩個函數的格點數據了。在筆者的2*E5-2696v3 36核服務器上僅僅11秒就算完了。之后選3把這兩個函數的格點數據分別導出為當前目錄下的func2.cub和func1.cub。將這兩個cub文件以及Multiwfn自帶的examples目錄下的RDGfill.vmd都拷到VMD目錄下，然后啟動VMD，在VMD的文本窗口輸入source RDGfill2.vmd執行此繪圖腳本，再輸入pbc box來繪制出盒子。之后在Graphics - Representation里選擇正處于CPK風格的那個representation，通過界面上的Sphere Scale把原子球尺寸調成0.7。此時的sign(λ2)rho著色的RDG等值面圖，也即一般說的NCI圖，如下所示

可見圖像效果非常理想。紅色的等值面在金剛石中呈網狀分布，體現出金剛石里面有顯著的位阻作用。讀者有興趣的話可以通過相同的方法去對硅也繪制這種圖，你會發現硅當中也有網狀的等值面，但顏色是綠色的，體現出由于硅的孔洞相對較大，孔洞中的位阻互斥就沒那么強了。

通常做周期性體系NCI分析時格點都按照此例方式來設即可，這樣格點會均勻分布在整個盒子里，從而把體系所有地方的弱相互作用都圖形化展現出來。格點間距應根據需要恰當設置，想要更光滑的等值面就把格點間距設得更小。注意在盒子范圍固定的情況下，計算耗時、cub文件的尺寸都和格點間距的三次方呈反比。此例如果把格點間距設成更小的0.13，你會發現得到的等值面上完全看不出任何鋸齒，但耗時會增加為原先的(0.2/0.13)^3=3.6倍。

NCI分析中經常涉及到繪制RDG與sign(λ2)rho之間的散點圖，做法在第1節里提到的筆者的RDG博文中，以及《繪制有填色效果的用于弱相互作用分析的RDG散點圖的方法》（http://www.shanxitv.org/399）中都做了介紹。周期性體系的NCI分析也可以這么繪制散點圖，這里就不再演示了，這和分析孤立體系的操作沒有任何區別。

3.1.2 準分子近似的NCI分析

現在再演示一下如何對金剛石超胞做準分子近似的NCI分析。用上一節的diamond_222.inp或者diamond_222-MOS-1_0.molden當做Multiwfn的輸入文件都可以，因為它們都能給Multiwfn提供這種分析所需要的原子坐標和晶胞信息。這里用原胞diamond.cif作為輸入文件進行演示，請大家舉一反三。

啟動Multiwfn，載入diamond.cif，然后輸入
300 //其它功能(Part 3)
7 //對體系進行幾何操作
19 //構建超胞
2
2
2
-10 //返回
0 //返回主菜單。現在內存里記錄的金剛石的原胞已經成了2*2*2的超胞了
20 //圖形化分析弱相互作用
2 //準分子近似的NCI分析
9 //設置格點。同前，不再解釋
[按回車]
[按回車]
0.2 //格點間距（Bohr）

你會充分體會到準分子近似的NCI分析極快，一眨眼的功夫格點數據就算完了。之后選3導出cub文件，將導出的func1.cub、func2.cub以及Multiwfn自帶的examples目錄下的RDGfill_pro.vmd都復制到VMD目錄下。啟動VMD，在文本窗口里輸入source RDGfill_pro.vmd，然后輸入pbc box顯示盒子。之后在Graphics - Representation界面里改一下CPK原子球尺寸，然后選Isosurface那個representation，把Isovalue改成0.2并按回車。此時圖像如下，可見和普通NCI分析得到的圖像特征定性一致。

3.2 共價有機框架化合物（COF）中的弱相互作用

此例將對下面這個COF晶胞中的弱相互作用通過NCI方法進行圖形化展現。

此體系的cif文件是本文文件包里的COF_12000N2.cif。這個晶胞已經很大了，顯然沒必要再弄成超胞。如《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》（http://www.shanxitv.org/569）一文所述，X光衍射得到的晶體結構中氫的位置通常是不可靠的，所以本例我們在做NCI分析前用CP2K先優化一下氫原子的位置。

啟動Multiwfn并載入COF_12000N2.cif，然后輸入
COF_12000N2.cif
cp2k
COF_opt.inp
-2 //要求產生molden文件
-1 //設置任務類型
3 //優化結構但不優化晶胞
9 //設置原子凍結
optH //只優化氫原子
2 //選擇基組
10 //6-31G*
0 //產生輸入文件

這樣就得到了PBE/6-31G*只優化COF中氫原子位置的CP2K輸入文件COF_opt.inp。運行之，在36核服務器上十來分鐘就算完了。之后可以看到一批以COF_opt-MOS-1_為開頭的molden文件，它們記錄了優化過程中各幀的波函數。其中序號最大的是COF_opt-MOS-1_4.molden，這是最后優化完的結構下的波函數文件。將COF_opt.inp里的晶胞信息手動加入到COF_opt-MOS-1_4.molden文件里作為[Cell]字段。

用Multiwfn載入COF_opt-MOS-1_4.molden，然后按照與3.1節完全相同的操作步驟讓Multiwfn做普通的NCI分析。使用0.2 Bohr格點間距時，在筆者的36核服務器上花了不到一分鐘算完。然后用前例說的RDGfill2.vmd腳本在VMD里對Multiwfn導出的func1.cub和func2.cub作圖，結果如下

上圖非常清晰直觀地將COF晶胞中的兩層結構間的大面積pi-pi堆積展現了出來。苯環中央的紅色錐形等值面展現了環中的位阻作用。當前體系中強烈的N-H...O內氫鍵被非常藍的圓片型等值面表現了出來。

上圖一些等值面還不算特平滑，有的等值面邊緣鋸齒有點明顯，將格點間距降到0.15會明顯更好看，但耗時會多一倍多，每個cub文件增大到90兆。另外，上圖左側和右側的N-H...O內氫鍵的圓片型等值面中間凹陷了一些，像缺了一塊，這是因為這個內氫鍵很強，在作用區域有的地方電子密度超過了0.05 a.u.，而Multiwfn默認會把電子密度超過0.05 a.u.區域的RDG設為100來避免顯示出強作用區域（如化學鍵）的等值面。為解決這個問題，可以把Multiwfn目錄下的settings.ini中的RDG_maxrho從默認0.05設為改大到0.06 a.u.，這樣重新計算和繪制后就不會看到這個強內氫鍵的RDG等值面中間凹陷（被屏蔽掉）一塊了。

如果你想把上圖中下方的相鄰晶胞（-Z晶胞）的結構和RDG等值面也顯示出來，就進入Graphics - Representation界面，對每個Representation都進入Periodic標簽頁并勾選上-Z，即可看到下圖

在之前不顯示相鄰鏡像的NCI圖中，我們看到在當前晶胞的底端也出現了大面積綠色等值面，這是這部分結構與-Z鏡像晶胞的上層部分的相互作用所導致的，看起來很礙眼，怎么避免顯示出來呢？在Multiwfn中有兩種做法可以實現：
(1)讓盒子起點位置的Z坐標稍微大于0，并讓盒子Z方向的尺寸適度小于晶胞的Z尺寸。具體來說，選擇做NCI分析，進入設置格點界面后，還是選擇模式9，然后輸入以下內容
0,0,1 //盒子原點的X、Y坐標都為0，Z坐標從1 Bohr開始，使得格點不會在晶胞最底層出現
0,0,10.85 //讓盒子的前兩個邊長等同于晶胞的前兩個平移矢量長度，而令盒子的Z尺寸比晶胞Z尺寸小2 Bohr，使得格點不會在晶胞頂端分布（注：晶胞尺寸在Multiwfn載入文件后屏幕上直接顯示了，Z尺寸是12.85 Bohr，減去2就是這里寫的10.85 Bohr）
0.2 //格點間距，同前
用這樣得到的cub文件在VMD里作圖如下所示，可見很理想，只保留了晶胞中間兩層結構間的層間和層內相互作用。

(2)另一種避免顯示與鏡像晶胞間的RDG等值面的方法是在Multiwfn計算時要求不考慮相應方向的周期性，我最推薦這種做法。比如對此例不想考慮Z方向的周期性，就把Multiwfn的settings.ini里的ifdoPBCxyz從原先的1,1,1改成1,1,0，這里0就代表第三個晶胞方向不考慮周期性（由于當前的晶胞第三個方向就是Z方向，所以等于不在Z方向考慮周期性）。這么改過之后照常計算和繪圖即可，格點還是按原先方式設置，結果和上圖完全相同。這種做法不僅方便，省得特意考慮格點分布范圍，而且由于忽略了一個方向的周期性，計算耗時還下降了百分之幾十，一舉兩得。

下面是對當前體系按照3.1.2節的過程繪制的準分子近似的NCI圖的結果，忽略了Z方向周期性。計算耗時極低，只花了1秒就算完了，圖像效果和普通NCI定性一致。

雖然準分子近似版NCI多數情況下和普通NCI定性一致，但是如果你對準確度要求比較高，比如考察氫鍵時想根據等值面藍色的深淺程度判斷強度，那么不要用準分子近似的版本，否則可能得到一些誤導性結論，畢竟用準分子近似相當于沒有考慮原子間相互作用引起的電子轉移、極化效應對NCI分析結果的影響。

3.3 氧化石墨烯與水的相互作用

在《使用Multiwfn非常便利地創建CP2K程序的輸入文件》（http://www.shanxitv.org/587）文中曾使用PBE-D3(BJ)/TZVP-MOLOPT-GTH對氧化石墨烯+一個水分子的表面吸附體系做了幾何優化，計算后得到了opt-1.restart文件，這其實是個CP2K的輸入文件，里面記錄了優化的最后一幀的信息。現在我們對這個結構做個NCI分析。首先基于這個文件創建一個PBE/6-311G**級別下的單點任務文件，用于得到molden文件。

啟動Multiwfn后輸入
opt-1.restart
cp2k
SP.inp
-2 //要求產生molden文件
2 //修改基組
11 //6-311G**
0 //產生輸入文件

用CP2K計算產生的SP.inp，算完后就有了SP-MOS-1_0.molden，然后把SP.inp里的晶胞信息手動寫入到這個molden文件作為[Cell]字段。

啟動Multiwfn，載入SP-MOS-1_0.molden。在計算前我們先進入主功能0看一下這個體系的結構特征，并且選圖形窗口菜單欄的Other settings里的Toggle showing cell frame把晶胞邊框顯示出來，
此時看到的圖像如下所示

做這個體系的NCI分析時要注意，由于石墨烯板緊貼著盒子最下方，因此如果直接將0,0,0作為盒子原點位置，屆時你會看到石墨烯里六元環中的位阻區域的RDG等值面被截斷，不好看。為避免這個問題，一個做法是用Multiwfn的主功能300的主功能7里的選項1先把整個體系往Z的正方向挪比如2 Bohr（約1埃），另一個辦法是定義格點的時候，讓盒子原點Z位置適當小于0，比如-2.0 Bohr。另外，當前體系水分子上方有大量的真空區，為了節約時間，可以讓Z方向盒子尺寸明顯小于晶胞的Z尺寸（設成其一半就足夠包圍感興趣的弱相互作用區域，可以節約一倍的計算時間）。下面計算時會考慮這兩點。

在當前圖形窗口的右上角點擊RETURN關閉窗口，然后輸入
20 //弱相互作用圖形化分析
1 //NCI分析
9 //借助晶胞信息定義格點
0,0,-2 //盒子原點位置（Bohr）
0,0,9 //當前晶胞Z方向尺寸是18.897 Bohr，令盒子Z尺寸取其一半，其它方向和晶胞尺寸相同
0.15 //格點間距
在筆者36核機子上13秒算完。算完后用選項3導出cub文件，和3.1.1節一樣在VMD里結合RDGfill2.vmd腳本繪制NCI圖，如下所示。為了讓大家了解當前盒子范圍，即格點數據計算的區域，在Graphics - Representation界面里我把對應Isosurface的那個representation里的Show下拉框從默認的Isosurface改為了Box+Isosurface。

從上圖可見，石墨烯每個環中央的位阻區域被清楚地通過紅色等值面展現了出來。水與氧橋之間有一塊淡藍色圓片型等值面，體現出氫鍵的存在，但顏色不算很深，說明不算很強的氫鍵。在水分子與碳的接觸區域還有一塊綠色的等值面，體現出水和氧化石墨烯之間還有明顯的范德華作用區域。上圖白色方框展現的是盒子范圍，可見盒子取得比較恰當，既沒有截斷任何等值面，也沒有明顯的浪費。由于當前體系與Z方向的鏡像距離比較遠，不會與其相互作用導致出現額外的RDG等值面，此例就沒有刻意把settings.ini里的ifdoPBCxyz像前例一樣改成1,1,0，但如果這么改一下的話倒是可以通過忽略Z方向周期性節約一些耗時。

3.4 C60晶體

此例考察一下C60富勒烯分子晶體中的弱相互作用，本文文件包里C60.cif是此體系的實驗結構。這個體系的晶胞里原子數很多，有240個，若用普通NCI分析耗時會較高，而且此體系主要是pi-pi堆積作用，這種作用靠準分子近似的NCI就可以合理地展現出來，所以此例就只用準分子近似的NCI了。

啟動Multiwfn，然后輸入
C60.cif
20
2 //準分子近似的NCI
9
[回車]
[回車]
0.15 //讓圖像較平滑，用較小的格點間距

計算完后導出cub文件并用RDGfill_pro.vmd在VMD里作圖，如下所示

這個圖像還不是很令人滿意，主要是紅色的表現位阻區域的等值面太多了，弄得圖像很凌亂，不便于觀察富勒烯之間的相互作用。為解決此問題，一方面可以用后文說的IGM，另一方面可以把sign(λ2)rho數值明顯大于0的區域手動屏蔽掉。具體做法是在Multiwfn做NCI分析的后處理菜單里選-1先把RDG vs sign(λ2)rho的散點圖繪制出來，如下所示

由圖可見在sign(λ2)rho>=0.05的區域有很多戳到底的spike，它們是導致分子內位阻作用的等值面出現的原因。為了屏蔽掉它們，在Multiwfn后處理菜單中選-2 Set RDG value where value of sign(λ2)rho is within a certain range，輸入0.02,999，然后輸入100。此時sign(λ2)rho數值大于0.02的區域的格點的RDG值就都被設為一個很大的數值100了，如果再重新繪制散點圖會看到下圖，可見sign(λ2)rho>0.02區域的spike都已經沒了。

現在重新選擇導出cub文件，然后再次用RDGfill_pro.vmd在VMD里作圖。并且在VMD里把+X的鏡像也顯示出來便于考察。然后在VMD文本窗口里輸入
color Name C tan
color change rgb tan 0.700000 0.560000 0.360000
這會把碳的顏色改成黃褐色，因為原本的青色與等值面顏色的對比度不夠強烈。再在Graphics - Representation里把分子結構顯示方式改為Licorice，并把Bond Radius改為0.1，之后在VMD Main窗口里選Display - Orthographic改為正交視角。此時的圖像如下

可見當前的等值面清楚地把每個富勒烯與相鄰富勒烯之間的pi-pi堆積展現了出來。肯定有人注意到，上圖中央的富勒烯中間有些C-C鍵是斷開的，而且圖像中央的RDG等值面中間有個黑色的細縫。這是因為VMD沒法識別和顯示當前晶胞里的原子與相鄰鏡像晶胞里的原子間的鍵，而且等值面在跨盒子的地方不可避免地會顯示出個縫隙。如果想避免這倆問題讓圖像更完美，就只能在Multiwfn做這個準分子近似的NCI分析前，先用主功能300的子功能7的選項19把體系構造成2*1*1的超胞，然后再照常做NCI分析。此時計算耗時會明顯增高，cub文件大小也會增加一倍。

3.5 Ih晶型的冰

此例用普通NCI方法考察Ih晶型的冰中的弱相互作用。Ih是常壓下不是特別低溫的情況時最穩定的冰的晶型，也是日常生活中我們接觸到的冰的晶型。大家有興趣也可以類似地考察其它晶型的冰中的弱相互作用。

先產生優化冰晶胞中氫位置的CP2K輸入文件。啟動Multiwfn，載入本文文件包里的H2O-Ice-Ih.cif，然后輸入
cp2k
opt.inp
-2 //產生molden文件
2 //選擇基組
11 //6-311G**
3 //開啟色散校正
2 //DFT-D3(BJ)
-1 //設置任務類型
3 //優化結構但不優化晶胞
9 //設置凍結
optH //只優化氫
0 產生輸入文件

現在用CP2K進行計算，很快就算完了。把晶胞信息手動寫入對應優化好的結構的opt-MOS-1_5.molden文件里作為[Cell]字段。由于冰當中的氫鍵作用相對較強，為了避免氫鍵對應的RDG等值面被截斷，把RDG_maxrho設為0.07。然后照常做NCI分析。由于當前晶胞不大，用很小格點間距要算的格點數也不很多，因此為了圖像效果盡量精甚細膩，此例用了0.13 Bohr這樣很小的格點間距。

照常使用RDGfill2.vmd腳本在VMD里對Multiwfn做NCI分析后導出的cub文件作圖。為了讓冰當中的相互作用顯得比較完整，在Graphics - Representation界面里把顯示分子結構和顯示等值面的representation都設成+X、+Y、-Y、+Z、-Z鏡像都顯示，并且設置霧化效果使得遠處的物體被一定程度遮掩而避免擾亂視線。具體來說就是確保VMD Main的Display里的Depth cueing已經打開了，然后在Display - Display settings里把Cue Mode設成Linear，并把Cue Start和Cue End分別設為1.5和3.0。此時看到的圖像如下。

上圖非常清楚地展現出在冰當中，每個水與周圍四個水形成了氫鍵，而且RDG等值面非常藍，體現出氫鍵強度較大。零星的綠色等值面體現出冰當中還有一些主要算是范德華作用的區域，但可以忽略不計。

3.6 NaCl表面吸附CO

本例繪制NaCl表面吸附一氧化碳的NCI圖。根據諸多文章的研究，最穩定的吸附結構應當是一氧化碳的碳沖著Na+，此例也優化的是這種結構。用的是兩層厚度3*3的氯化鈉板，底層原子凍結，上層和CO都在優化過程中自由弛豫。先通過PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH優化結構，然后PBE/6-31G*算單點得到波函數文件。此模型可以用GaussView很容易地搭建，保存成gjf并載入后通過Multiwfn很容易地就可以創建相應任務的CP2K輸入文件，相關過程這里就不累述了。單點計算后得到的SP-MOS-1_0.molden在本文文件包里提供了。之后對此體系做普通NCI分析，用0.15 Bohr格點間距，圖像如下所示。此體系計算耗時很高，為了節約時間，應恰當設置盒子Z方向尺寸，避免讓格點分布在根本沒原子出現的真空區。

為了看得清楚，俯視圖和側視圖也在下面給出

圖中青色圓球是Cl-，藍球是Na+，褐球是碳，紅球是氧，藍色方框展現的是計算的格點數據的分布范圍。從此圖可以清楚地看到一氧化碳的碳與Na+有明顯的相互作用，與旁邊的四個Cl-也有明顯作用，并且Na+和Cl-之間也有顯著相互作用。Na+帶顯著正電，它與碳的鍵軸末端的孤對電子明顯會形成顯著的靜電主導的吸引作用（文獻中的吸附能約為20 kJ/mol，不算小了），且Na+和Cl-之間毫無疑問是挺強的離子鍵，為何RDG等值面全都是綠色而不是藍色？這是因為鈉的半徑較大，而且價層孤立狀態下就只有一個電子，在當前體系中價電子就更少了，故在它與其它原子形成相互作用的區域必定電子密度（rho）比較低，所以sign(λ2)rho數量級很小，故著色是綠色。記住切勿盲目用RDG等值面顏色衡量作用強度，要正確理解NCI方法的本質和局限性。

上面的圖有個明顯的不足之處就是Na+和Cl-之間作用的等值面嚴重妨礙了觀看氯化鈉和一氧化碳之間的相互作用。既可以使用后文所述的IGM方法避免此問題，也可以使用Multiwfn的格點數據處理功能（主功能13）中的子功能14解決，它可以用來將NaCl板與一氧化碳交疊區域以外的區域的格點的RDG設成一個特定值，比如100，從而達到屏蔽NaCl內部相互作用的目的。下面就來這么屏蔽一下。NCI分析后導出的func2.cub是RDG的格點數據文件，將之載入Multiwfn，然后輸入
13 //主功能13
14 //設置兩個片段交疊區域以外區域格點的數值
1.3 //原子球半徑設為范德華半徑乘以1.3（片段占據的區域是片段內所有原子球的疊加）
100 //交疊區域以外的格點數值設為100
2 //手動輸入兩個片段里的原子序號
1-72 //片段1的原子序號，即NaCl板中的原子
73,74 //片段2的原子序號，即CO的原子
0 //將當前的格點數據導出為cub文件
[按回車] //格點數據導出為當前目錄下的func2.cub

現在，用新得到的func2.cub替換原先的func2.cub，再用RDGfill2.vmd腳本作圖，得到如下圖像，可見非常理想，確實NaCl內部的等值面都被屏蔽掉了。為了效果更好，作圖時筆者還設了景深霧化，并且把RDG等值面數值從默認的0.5稍微設大到了0.6，這樣等值面更膨脹一些，邊緣鋸齒看起來也更小一點。

順帶一提，還可以在VMD中使用焦距虛化效果，如下所示，這樣層次更清楚。也可以通過把原子球調大來提升層次感。

像上面這樣只需要對一小部分區域作NCI圖的情況，在Multiwfn里設置計算的格點范圍的時候可以只讓盒子框住感興趣的區域（即感興趣的等值面可能出現的區域），耗時能節約非常多。這通過恰當設置盒子原點和邊長即可實現，也可以在格點設置界面里通過選項10 Set box of grid data visually using a GUI window在圖形窗口里直觀設置盒子的位置和尺寸，正如此文2.6節所示的情況：《使用量子化學程序基于簇模型計算金屬表面吸附問題》（http://www.shanxitv.org/540）。

4 周期性體系的IGM分析

IGM分析只需要原子坐標。如果體系雖然是周期性的，但是被考察的區域與晶胞邊界有一定距離，那么完全可以按照《通過獨立梯度模型(IGM)考察分子間弱相互作用》（http://www.shanxitv.org/407）所述的當成孤立體系考察。比如Mater. Today Commun., 26, 102028 (2021)一文中，在筆者的建議下作者對他們研究的MOF吸附芳香化合物的問題做了IGM分析，如下所示，可見很好地展現出了MOF對分子的吸附時主要的作用區域。像這種情況就當成孤立體系計算做IGM計算就可以，考慮周期性的話耗時會明顯更高。

但也有一些情況IGM等值面要涉及到體系邊緣，這種情況就必須當成周期性做IGM分析了，否則盒子邊緣處的等值面就不正常了，比如3.2節的COF體系就是如此。這里我們就用IGM方法對COF這個體系做一下分析，如果沒讀過《通過獨立梯度模型(IGM)考察分子間弱相互作用》一文的話請仔細閱讀，基本知識本文不再累述。為了不顯示出中心晶胞和相鄰晶胞的相互作用對應的IGM等值面，還是把settings.ini里的ifdoPBCxyz設成1,1,0忽略掉Z方向的周期性。

啟動Multiwfn，然后輸入
COF_opt-1.restart //本文文件包里優化COF中氫原子位置任務產生的記錄最終結構的restart文件
20 //圖形化分析弱相互作用
10 //IGM
2 //定義兩個片段
1-72 //第一個片段的原子序號，即兩層COF其中一層中的原子
73-144 //第二個片段的原子序號，即另一層COF中的原子
9 //借助晶胞信息設置格點
[按回車] //用0,0,0作為盒子原點
[按回車] //用晶胞邊長作為盒子邊長
0.4 //格點間距。IGM等值面不容易像RDG的那么容易有鋸齒，所以可以用明顯大得多的格點間距節約時間

本身IGM分析就很快，而且格點間距又大，在普通4核機子上不到半分鐘就算完了。之后選3導出格點數據，在當前目錄下就有了dg_inter.cub、dg_intra.cub、dg.cub和sl2r.cub，將它們都挪到VMD目錄下。并且把Multiwfn目錄下的examples文件夾中的IGM繪圖腳本IGM_inter.vmd和IGM_intra.vmd都復制到VMD目錄下。

啟動VMD，然后在命令行窗口輸入source IGM_inter.vmd，這會繪制出IGM方法中定義的sign(λ2)rho著色的δg_inter函數的等值面。但是直接顯示出來的等值面非常窄，不能充分展現出兩層MOF間的pi-pi堆積，這是因為默認的δg_inter等值面數值偏大。因此進入Graphics - Representation界面，把Isovalue改成更小的0.005，此時圖像如下所示。可以看到此圖非常清楚地展現出了層間pi-pi堆積作用，而且片段內的相互作用不像NCI那樣被顯示出來，同時等值面很平滑，可見Multiwfn做IGM分析又簡單又快效果又非常好。IGM分析還有一個好處是可以把各個原子的貢獻給出來，并可以由此給原子球著色凸顯較重要的原子，詳見前述的《通過獨立梯度模型(IGM)考察分子間弱相互作用》。

如果在VMD里輸入source IGM_intra.vmd，可以把COF層內的相互作用，包括化學鍵作用和內氫鍵都展現出來（展現后者需要調小isovalue）。由于這不是此例我們感興趣的，就不說了。

下面對本文3.6節考察過的NaCl吸附CO的體系做IGM分析，將NaCl和CO分別作為兩個片段，操作方式同上。下圖左側和右側對應0.005和0.003的δg_inter等值面。

δg_inter數值取得越小，其等值面就可以把越弱的相互作用展現出來。如上所示，在0.005等值面的時候，C和Na+的相互作用區域很明顯，但C和周圍Cl-的范德華作用對應的等值面只顯示出了CO稍微歪過去方向的那一小塊。而將等值面數值降到0.003后，更多相互作用區域被等值面展現了出來，并且此時C和Na+接觸的那部分等值面中央顏色已經很藍了。為什么此時顏色這么藍而不像RDG等值面那樣是綠色的？這是因為此時的δg_inter等值面特別肥厚，已經延展到了比較靠近C的部分，此處電子密度已經不是很小了，sign(λ2)rho可以達到負零點零幾了，因此對應的著色明顯發藍。

IGM等值面比較肥厚是IGM方法的十分明顯的弊端，而且上面的顏色往往還會誤導研究者，在筆者2022年提出的IGMH方法中被完美解決，而且IGMH比IGM原理更嚴格。讀者務必看《使用Multiwfn做IGMH分析非常清晰直觀地展現化學體系中的相互作用》（http://www.shanxitv.org/621），里面有IGMH非常全面的介紹和與IGM的對比。在Multiwfn中的使用上，IGMH和IGM的區別僅僅在于主功能20里應選擇11而非10，并且需要載入和NCI分析一樣的含有波函數信息的輸入文件。IGMH由于是基于波函數的分析，故耗時顯著高于IGM，而且和NCI分析一樣格點間距不宜太大以避免鋸齒。對當前體系，為節約耗時，建議設置格點時讓盒子僅框住感興趣的作用區域，并且把ifdoPBCxyz設為0,0,0當孤立體系處理（因為CO與NaCl作用的位置離盒子邊緣很遠，這么做不會有不良影響）。基于IGMH方法繪制的δg_inter為0.002的等值面圖如下所示，可見等值面和RDG的等值面很相似，都比較薄，而且和IGM一樣都完全避免了NaCl板內的等值面妨礙考察CO的吸附，圖像效果極度理想。我建議凡是能得到波函數的情況，都應當用IGMH代替IGM。

5 總結

本文詳細介紹了如何使用Multiwfn基于CP2K第一性原理程序產生的波函數做NCI分析，還介紹了怎么對周期性體系做IGM分析，還對筆者提出的比IGM強得多的IGMH分析也做了簡要示例。Multiwfn靈活方便的創建CP2K輸入文件的功能，結合Multiwfn非常簡單易用的NCI/IRI/IGM/IGMH分析功能，真是使得固體和表面的弱相互作用圖形化分析異乎尋常地簡單。研究周期性體系中的弱相互作用不要再只是算算結構、算算能量、只靠數據討論了，結合本文介紹的圖形化分析方法可以令文章大為增光添彩、明顯更吸引眼球，令枯燥抽象的分析討論變得直觀生動傳神。

Multiwfn對周期性體系的分析絕不限于本文所介紹的。目前的版本已可以對周期性體系計算原子電荷、計算鍵級、做AIM拓撲分析、計算和繪制上百種實空間函數（如ELF、LOL、IRI、電子密度拉普拉斯、能量密度、自旋密度等等）、做范德華勢分析（《談談范德華勢以及在Multiwfn中的計算、分析和繪制》http://www.shanxitv.org/551），等等。Multiwfn對周期性體系的支持詳情看Multiwfn手冊2.9節的說明。筆者在未來會陸續寫更多的Multiwfn對周期性體系波函數分析的相關教程。讀者也別忘了看《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598），其中有IRI研究周期性體系的介紹，比NCI分析能展現更豐富的信息，建議代替NCI。

順帶一提，如果你之前不會用CP2K的話，也完全不用擔心什么。如本文所示，靠Multiwfn可以非常簡單地創建本來特別繁瑣復雜的CP2K輸入文件，而且如前文提到的我的CP2K安裝教程所寫的，用CP2K預編譯版都完全不需要自己編譯，下載就能用。再加上CP2K又快又免費，結合Multiwfn做周期性體系的波函數分析簡直不能更好用。所以即便你以前用的是Quantum ESPRESSO、CASTEP等其它第一性原理程序，也可以毫無壓力地過度到Multiwfn+CP2K的分析組合上。

使用Multiwfn非常便利地創建CP2K程序的輸入文件

sobereva — Mon, 22 Feb 2021 22:15:00 +0800

使用Multiwfn非常便利地創建CP2K程序的輸入文件

文/Sobereva@北京科音

First release: 2021-Feb-22 Last update: 2023-Nov-7

1 前言

CP2K（https://www.cp2k.org）是速度超級快，功能又非常全面的第一性原理程序，受到越來越多的重視，有很多無可替代的重要價值。CP2K的一個問題是輸入文件寫起來極其繁瑣，一個簡單的任務的輸入文件就要動輒上百行。而且有時用戶還需要了解很多算法和數值實現細節、查閱繁復的手冊，甚至還要google半天、參考不少自帶測試文件才能最終寫出能用的輸入文件。這些問題使得CP2K這個優秀的程序不容易普及開來，初學者往往看到其極度復雜的輸入文件就打退堂鼓了。

筆者的Multiwfn程序（主頁&下載：http://www.shanxitv.org/multiwfn）的主功能100的子功能2能產生大量量子化學程序的輸入文件，諸如其中產生ORCA輸入文件的功能在此文專門介紹過《詳談Multiwfn產生ORCA量子化學程序的輸入文件的功能》（http://www.shanxitv.org/490），給很多人帶來了極大便利。為了讓CP2K這個利器對于常見計算問題用起來盡可能方便、盡可能把手動擋變得像自動擋、令使用門檻盡可能低，使之能更充分地發揮其巨大價值，筆者在Multiwfn里加入了產生CP2K輸入文件的功能，在此文將進行簡要介紹。這個功能雖然注定不可能讓完全零基礎的CP2K用戶毫無壓力地使用CP2K，但至少已經把使用難度和復雜度降到極限了，用戶只需要選擇干什么、怎么干，即可產生一個完成度>95%的輸入文件，之后一些必須根據實際情況調整的地方，比如平面波截斷能大小等，再自行改一下即可，這使得CP2K用戶得到了極大的解放！

如果Multiwfn的這個功能對你的研究產生了實際幫助，請在寫文章時順帶提及輸入文件是借助Multiwfn產生的（比如寫the CP2K input files were generated with help of Multiwfn program），并按程序啟動時顯示的要求恰當引用Multiwfn。

Multiwfn產生CP2K輸入文件的流程是：
1 啟動Multiwfn，載入一個Multiwfn可以識別的至少含有結構信息的文件
2 進入主功能100的子功能2，選擇產生CP2K輸入文件（即選項25）。也可以直接在Multiwfn主菜單里輸入cp2k進入此功能
3 輸入要產生的CP2K輸入文件的路徑
4 通過各種選項設置如何進行計算
5 選擇0，得到CP2K輸入文件

下面第2節先說載入到Multiwfn的輸入文件的要求。第3節給出一些例子演示產生CP2K輸入文件的過程。Multiwfn的這個功能始終不斷打磨并且迎合CP2K最新版本的改變，本文內容對應Multiwfn官網上最新版本的Multiwfn（不是今天剛下載的就不要覺得是最新的），一定不要用老版本！Multiwfn產生CP2K輸入文件的功能在未來還會不斷擴展，盡量做到更強大更完美，并且盡量迎合未來版本CP2K發生的各種變化，屆時本文也會相應地修改。

一定要注意，盡管靠Multiwfn創建CP2K輸入文件非常方便，也是絕對絕對不能把CP2K當黑箱用的！！！因為真正正確、準確、高效率地使用CP2K做計算需要掌握的理論背景知識真是超級多。“北京科音CP2K第一性原理計算培訓班”是從頭一次性系統、完整學習第一性原理計算和CP2K程序的最好途徑，是CP2K用戶絕對不應錯過的培訓，介紹見http://www.keinsci.com/workshop/KFP_content.html。此培訓深入淺出，能非常全面學習到第一性原理計算各種問題的背景知識、深入透徹掌握CP2K，比起自己在網上找零七八碎的往往還是有錯誤或過時的資料學習能少走無數彎路。培訓內容涵蓋面極廣，ppt多達兩千多頁，對不同主題都給了精心準備的大量例子使得學員能舉一反三。其中講解利用CP2K做各類問題的研究時都會使用Multiwfn創建輸入文件，學員從中會深切體會到Multiwfn對于CP2K的使用起到的重要價值。

2 給Multiwfn用的輸入文件

Multiwfn支持什么輸入文件這里都詳細說了：《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。對于產生CP2K的輸入文件，顯然載入到Multiwfn的文件中必須含有體系的結構才行。對于此類目的，Multiwfn支持一大波格式，如pdb/gro/xyz/mol/mol2/mwfn/wfn/wfx/molden/fch/cub/gjf/gro等等等等。

如果載入到Multiwfn的文件不包含晶胞信息，那么Multiwfn會默認將此體系視為孤立體系來設置CP2K里的計算參數，CP2K輸入文件里的晶胞信息是在分子周圍做一定距離的延展來得到的。

多數情況CP2K都是拿來算周期性體系的。若想讓Multiwfn產生的CP2K輸入文件里有實際的晶胞信息，載入到Multiwfn的文件里就必須有晶胞信息。用以下格式作為輸入文件時，Multiwfn會既讀取原子坐標也讀取晶胞信息。
(1)cif文件
(2)含有CRYST1字段的pdb或pqr文件
(3)CP2K的輸入文件（必須以.inp或.restart為后綴）
(4)GROMACS的gro文件
(5)含有@CRYSIN字段的mol2文件
(6)含有Ndim且其數值>0的mwfn文件（mwfn是Multiwfn私有的波函數文件格式，專門留出了字段記錄晶胞信息，見mwfn介紹文章https://doi.org/10.26434/chemrxiv.11872524）
(7)帶有Tv（平移矢量）信息的Gaussian輸入文件，即Gaussian做PBC任務的輸入文件，可以通過GaussView直接創建
(8)Gaussian做PBC任務產生的fch/fchk文件
(9)含有晶胞信息的xyz文件。這是我自己定義的，如果xyz文件的注釋行（第二行）包含諸如以下形式的內容（以埃為單位的平移矢量），就會被Multiwfn讀取
Tv_1: 7.426 0.0 0.0 Tv_2: 3.6 6 6.40 0.0 Tv_3: 0.0 0.0 10.0
也可以按照extended xyz格式的方式在xyz文件的第二行記錄晶胞信息，例如：
Lattice="7.426 0.0 0.0 -3.66 6.40 0.0 0.0 0.0 10.0"
(10)含有[Cell]字段的molden文件。此字段是我對molden格式做的擴展，這使得記錄波函數的molden文件也能記錄晶胞信息。大家可以直接手動往已有的molden文件里插入[Cell]信息，詳見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）
(11)VASP的POSCAR、CHGCAR、CHG、ELFCAR、LOCPOT文件。文件名里必須包含相應字樣才能識別成相應格式，例如要作為POSCAR載入，文件名可以是POSCAR_miku或miku.POSCAR
(12)Quantum ESPRESSO的輸入文件

3 產生CP2K輸入文件的一些例子

下面就通過不同體系、不同類型計算演示一下Multiwfn產生CP2K輸入文件的功能，這些例子也可以作為CP2K入門例子。這些例子只用到了一部分Multiwfn提供的選項，其它的請自行嘗試，選項的含義都非常易于理解所以就不一一說了。下文用到的文件都在http://www.shanxitv.org/attach/587/file.zip中提供了。筆者假定讀者已經以正當方式安裝了CP2K，安裝方法見《CP2K第一性原理程序在CentOS中的簡易安裝方法》（http://www.shanxitv.org/586）。由于本文只是舉例，為了省計算時間，有些任務的k點取得比實際研究時應取的明顯要少。

3.1 做COF（共價有機框架）化合物的單點任務

本文文件包里COF_12000N2.cif是個COF的晶體結構文件，本例用Multiwfn對它產生一個最簡單的任務的輸入文件，即PBE泛函結合DZVP-MOLOPT-SR-GTH基組做單點計算。

啟動Multiwfn，然后依次輸入
COF_12000N2.cif //寫實際的路徑。也可以把此文件拖到Multiwfn文本窗口里自動產生路徑，省得手動輸入了
cp2k //進入產生CP2K輸入文件的功能
[按回車] //新產生的CP2K輸入文件將為當前目錄下的COF_12000N2.inp。也可以自己輸入路徑
0 //直接生成輸入文件

現在COF_12000N2.inp就在當前目錄下出現了，可以直接用CP2K跑了。這樣在默認設置下產生的輸入文件對應PBE/DZVP-MOLOPT-SR-GTH算單點，只考慮gamma點。可以打開.inp看一下有哪些自己需要改的，一般來說自己根據實際需要改的也就是影響精度的設置，比如影響平面波部分計算精度的&MGRID里的CUTOFF和REL_CUTOFF、控制密度矩陣收斂限的&SCF里的EPS_SCF，這些問題本文就不多說了。

Multiwfn產生的CP2K輸入文件里面有很多注釋，即#后面的內容，這有助于對CP2K還不是特別熟的人了解選項的含義、明白怎么設置。輸入文件里有些選項的值是根據我的使用經驗和習慣設的，也有一些設置雖然是默認的，卻也出現在了輸入文件里，這個考慮是便于大家之后自己手動修改（要不然想改一個默認設置的話還得照著手冊寫一堆字段實在太費事了）。

下面再用COF舉個例子，使用B97M-rV泛函結合6-31G*基組做GAPW形式的單點計算，并且要求算完后產生記錄此體系波函數的.molden文件，可用于之后在Multiwfn中做波函數分析，見比如《使用IRI方法圖形化考察化學體系中的化學鍵和弱相互作用》（http://www.shanxitv.org/598）和《使用Multiwfn結合CP2K通過NCI和IGM方法圖形化考察固體和表面的弱相互作用》（http://www.shanxitv.org/588）。關于molden文件更多信息見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》（http://www.shanxitv.org/651）。本例用較小的基組僅作為演示目的，對精度有要求的計算至少用6-311G**。值得一提的是，在純泛函中，B97M-rV是對有機體系能量相關問題計算得最好之一，并且這個泛函自帶rVV10色散校正。

啟動Multiwfn后依次輸入
COF_12000N2.cif
cp2k
COF_GAPW.inp
1 //設置理論方法
11 //B97M-rV（通過libxc庫實現）
2 //設置基組
10 //6-31G*
-2 //要求產生.molden文件
0 //產生輸入文件
現在當前目錄下就有了COF_GAPW.inp，直接就可以跑了。跑完之后當前目錄下會出現COF_GAPW-MOS-1_0.molden。

順帶一提，在Multiwfn的主功能0里可以直接看體系的結構和晶胞，由此可以確認一下結構和晶胞信息有沒有被Multiwfn從輸入的文件里正常載入。啟動Multiwfn并載入COF_12000N2.cif后，選0，就會進入圖形界面。關閉標簽和坐標軸，并且點擊菜單欄的Other settings里的Toggle showing cell frame，當前窗口將為下圖所示的狀態。可見原子位置和晶胞都正確。點擊Save picture按鈕可以在當前目錄下得到線條明顯更平滑的大尺寸圖像文件。

也可以用Multiwfn載入CP2K的.inp輸入文件或者CP2K計算過程中產生的.restart文件，在主功能0里可以看其中的結構是什么樣。還可以載入它們后，進入主功能100的子功能2，選擇導出gjf文件，這樣得到的gjf文件里會有Tv描述的晶胞信息，可以放到GaussView里進行觀看和編輯修改。還可以用主功能100的子功能2里的一堆選項轉換成其它一大堆常見格式、一大堆計算程序的輸入文件（包括Quantum ESPRESSO、VASP的POSCAR），賊靈活。

3.2 對SiC超胞做結構優化

這一節演示對SiC的3*3*3超胞用PBE/DZVP-MOLOPT-SR-GTH做晶胞參數不變的結構優化，并且開OT。

啟動Multiwfn后輸入
SiC.cif //在本文文件包里
cp2k
[按回車] //產生的輸入文件名為默認的SiC.inp
-1 //選擇任務
3 //結構優化（不變胞）
-11 //進入幾何操作界面。這個界面有豐富的功能，見《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）
19 //構建超胞
3 //方向1為原先的3倍
3 //方向2為原先的3倍
3 //方向3為原先的3倍
如果你想看一下當前結構的話，可以選0進入圖形界面，會看到下圖，可見結構沒問題

接著輸入
-10 //返回到之前創建CP2K輸入文件的界面
4 //開啟OT
0 //產生輸入文件

現在當前目錄下的SiC.inp就可以直接跑了。

實際上，還有一種等價的創建SiC 2*2*2超胞輸入文件的做法。進入Multiwfn的CP2K輸入文件創建界面后，輸入
4    //開啟OT
-1   //選擇任務
3    //結構優化（不變胞）
-9 //其它設置
3 //設置三個方向晶胞重復次數
3,3,3
0 //返回
0 //產生輸入文件

在得到的輸入文件中可以看到坐標部分只記錄了原胞，但有兩處MULTIPLE_UNIT_CELL 3 3 3，體現出是對3*3*3超胞進行計算。

3.3 優化Cu的晶胞

這一節創建的輸入文件對應于在PBEsol/DZVP-MOLOPT-SR-GTH級別下對Cu晶體做原子坐標和晶胞的優化。由于是導體，所以用了smearing。由于用的是原胞算的，所以此例設了k點。

啟動Multiwfn后依次輸入
Cu.cif //在本文文件包里，是銅晶體的cif文件
cp2k
Cu_cellopt.inp //在當前目錄下產生Cu_cellopt.inp
1 //設置理論方法
-3 //PBEsol，很適合算固體晶格常數的泛函
8 //設置k點
3,3,3
-1 //設置任務
4 //優化原子坐標和晶胞
6 //使用占據數smearing
0 //產生輸入文件

現在得到的Cu_cellopt.inp就可以直接算了。

3.4 氧化石墨烯+水分子的結構優化

這一節我們做一個氧化石墨烯+水分子的結構優化。具體來說，是在單層石墨烯的一個C-C鍵上加一個氧橋，然后放一個水分子到合適位置，與氧橋形成氫鍵。構建這個初始結構可以利用GaussView之類的程序。由于文字+圖片方式表達起來太費事，我把GaussView里的操作錄了個短小的視頻：http://www.shanxitv.org/attach/587/gview_build.mp4，其中旋轉水分子的時候按住alt鍵，平移水分子的時候按住alt+shift鍵。按照視頻建模完畢后，保存成gjf文件，就是此文文件包里的graphite_oxide_H2O.gjf，可見里面有Tv（平移矢量）信息。

此例的計算將使用PBE-D3(BJ)/TZVP-MOLOPT-GTH完成。由于當前體系有弱相互作用，所以加了DFT-D3(BJ)校正，相關知識見《談談“計算時是否需要加DFT-D3色散校正？”》（http://www.shanxitv.org/413）。

啟動Multiwfn，然后依次輸入
graphite_oxide_H2O.gjf
cp2k
opt.inp //新產生的文件名
-1 //設置任務
3 //優化結構（晶胞尺寸不變）
3 //設置色散校正
2 //DFT-D3(BJ)
2 //設置基組
3 //TZVP-MOLOPT-GTH
-7 //設置周期性
XY //當前體系是XY周期性，Z方向不考慮周期性
0 //產生輸入文件

之后就可以直接跑了。

順帶一提，在CP2K輸入文件生成界面里選擇優化任務后，可以通過選項9設置對哪些原子坐標進行凍結，可以通過諸如1,5,9-12,14-18形式方便地輸入。這種選擇語句還可以在GaussView里方便地生成，即把某些原子選中成黃色后，在Tools - Atom selection界面里就可以直接拷出來粘貼到Multiwfn窗口里。另外，Multiwfn讓你輸入被凍結的原子序號時若輸入optH還可以直接實現凍結重原子而只優化氫原子的目的。

3.5 用GFN1-xTB方法跑含有216個水的盒子的動力學

GROMACS程序自帶了一個做NPT動力學預平衡好的含有216個水的結構文件spc216.gro，在本文文件包里也給了。此例基于這個文件，創建一個用GFN1-xTB級別做NVT系綜的從頭算動力學的輸入文件。GFN1-xTB是半經驗層面的DFT，簡要介紹見《盤點Grimme迄今對理論化學的貢獻》（http://www.shanxitv.org/388）和《將Gaussian與Grimme的xtb程序聯用搜索過渡態、產生IRC、做振動分析》（http://www.shanxitv.org/421）。

啟動Multiwfn然后輸入
spc216.gro
cp2k
spc216_AIMD.inp
-1 //設置任務
6 //分子動力學
10 //設置熱浴
2 //CSVR熱浴，亦即GROMACS里的V-rescale
1 //設置理論方法
30 //GFN1-xTB
0 //創建輸入文件

在得到的spc216_AIMD.inp中，應結合注釋根據實際需要做恰當修改，比如控制模擬步數的&MD里的STEPS、控制初始溫度和控溫溫度的TEMPERATURE。諸如此類直接就可以手動很容易地改的參數，在Multiwfn的CP2K輸入文件創建界面里一般就沒有留出相應的選項了，因為完全沒必要。

注：實際上 GFN1-xTB模擬純水效果很差，密度明顯偏高。此例只不過是個例子而已。

值得一提的是，在CP2K輸入文件創建界面里選擇動力學任務后，會看到還有幾個相關選項：
選項12：設置壓浴
選項-5：設置動力學軌跡格式。不用控壓時默認是xyz，用控壓時默認是尺寸明顯更小而且記錄盒子信息的二進制格式dcd。也可以用pdb格式，記錄盒子信息，文件尺寸很大
選項-6：設置每多少步把坐標寫入一次軌跡文件

3.6 對HCN氫轉移找過渡態、做振動分析

此例對HCN->NCH氫轉移用dimer方法優化過渡態，之后再做個振動分析考察虛頻情況。此例只是做個演示，對這種分子體系找過渡態的問題CP2K遠遠不如Gaussian，又不好用效率又極低。

首先用GaussView創建個看上去像這個異構化反應的過渡態的結構，保存為H2CO_TS.gjf。然后啟動Multiwfn，輸入
H2CO_TS.gjf
cp2k
dimer.inp
-1 //選擇任務
7 //搜索過渡態
0 //產生輸入文件

然后用CP2K運行產生的dimer.inp。算完后會看到有一個dimer-1.restart文件，這是CP2K輸入文件，包含的是dimer任務最后一步的結構。基于這個結構，我們創建一個振動分析的輸入文件。

啟動Multiwfn，然后輸入
dimer-1.restart
cp2k
freq.inp
-1 //選擇任務
5 //振動分析
0 //產生輸入文件

然后用CP2K運行產生的freq.inp。跑完之后可以看到有且只有一個虛頻，目錄下還出現了freq-VIBRATIONS-1.mol，里面用molden格式記錄了振動矢量，可以用Molden程序觀看振動模式。

3.7 計算氮化硼二維材料的REPEAT原子電荷

做涉及固體表面的經典力場的分子動力學模擬，需要有固體部分的原子電荷，REPEAT電荷對于這個目的是不錯的選擇，對于擬合固體表面的原子電荷比RESP更理想。此例創建一個CP2K算氮化硼板的REPEAT電荷計算的輸入文件。

本文文件包里有個氮化硼晶胞文件BN.cif，用GaussView打開，刪掉兩層BN中的任意一層，并把盒子Z方向尺寸設為10埃（真空區不夠大的話會令REPEAT電荷不合理），之后保存為BN.gjf。

啟動Multiwfn，然后輸入
BN.gjf
cp2k
BN_REPEAT.inp
-4 //輸出原子電荷
7 //REPEAT
8 //設置k點
6,6,6
0 //產生輸入文件

為了讓結果更好，手動打開.inp文件，把ELEMENT N那行下面的基組改為TZV2P-MOLOPT-GTH。這個基組對B沒有定義，所以對B還是用原先默認的DZVP-MOLOPT-SR-GTH。

然后用CP2K運行得到的BN_REPEAT.inp即可，算完后會得到一個.resp后綴的文件，里面就是各個原子的電荷值。B為0.478，很合理。當前目錄下還有個.xyz文件，拖入VMD里可以看擬合點的分布。Multiwfn產生的CP2K算REPEAT電荷的擬合點分布設置是根據我的經驗已調整到了比較合理的情況，一般不需要再做修改。

3.8 用TDDFT算激發態

請看專門的文章：《使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和考察電子激發態》（http://www.shanxitv.org/634）。

3.9 做雜化泛函計算

請看專門的文章：《CP2K做雜化泛函計算的關鍵要點和簡單例子》（http://www.shanxitv.org/690）。

4 總結

Multiwfn創建CP2K輸入文件的功能使得的CP2K程序使用變得相當簡單，大大降低了使用門檻，不用再照著手冊和例子一個輸入文件寫半天（而且還很容易寫錯、寫得不合理），而只需要在Multiwfn里敲幾下鍵盤就可以產生出一個基本能用的輸入文件，而且借助此功能還可以通過GaussView方便地為CP2K計算來建模。這個功能也很有助于初學者學習CP2K的輸入文件的編寫。

Multiwfn創建CP2K輸入文件的界面在未來還會加入更豐富的選項、讓CP2K的更多的功能用起來變得更簡單。在發表文章時提及輸入文件是借助Multiwfn創建的并恰當引用Multiwfn程序原文，是對Multiwfn這個功能繼續開發最好的鼓勵與支持。

Multiwfn創建CP2K輸入文件能做的任務遠遠不止上面演示的那些，還有NEB、NMR、路徑積分分子動力學、X吸收光譜，支持的方法還有MP2、雙雜化、GW、分子力場等等，在北京科音CP2K第一性原理計算培訓班（http://www.keinsci.com/workshop/KFP_content.html）里都有非常全面深入講解。