使用CP2K結合Multiwfn對周期性體系模擬UV-Vis光譜和分析電子激發態
使用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論壇上發帖咨詢。