使用Multiwfn計算激發態間的躍遷偶極矩和各個激發態的偶極矩
文/Sobereva @北京科音
First release: 2014-Mar-28 Last update: 2022-Apr-13
1 前言
計算激發態之間的躍遷偶極矩有一些用處，比如用于sum-over-states (SOS)方式計算(超)極化率（詳見http://www.shanxitv.org/232）、計算激發態到激發態的吸收光譜（用瞬態吸收光譜技術可以測定）。Gaussian的CIS、TDHF、TDA-DFT、TDDFT計算可以直接給出基態到激發態的躍遷偶極矩。CIS和TDA任務也提供了alltransitiondensities關鍵詞可以把激發態之間的躍遷偶極矩輸出出來，但是對于TDHF/TDDFT任務則沒有這個關鍵詞（雖然Gaussian有個density=transition=(n,m)關鍵詞可以產生n到m激發態的躍遷密度，然而這個關鍵詞在目前的版本中貌似功能不正常，而且光是有躍遷密度也沒用，還得利用偶極矩積分才能得到躍遷偶極矩，但Gaussian也并不給出這樣的信息）。
強大的波函數分析程序Multiwfn的電子激發分析模塊中的子功能5能夠計算所有態（包括基態和激發態）之間的躍遷電和磁偶極矩，可以用于上述研究。此外，這個功能還可以直接給出各個激發態的電偶極矩，通過與基態的偶極矩求差就可以考察電子激發過程中偶極矩的變化。在本文就簡要演示一下這個功能的使用。請讀者務必使用官網上最新版Multiwfn，可以在http://www.shanxitv.org/multiwfn免費下載。如果對Multiwfn一無所知的話，建議先閱讀《Multiwfn入門tips》（http://www.shanxitv.org/167）和《Multiwfn FAQ》（http://www.shanxitv.org/452）。
對于Gaussian用戶，使用Multiwfn的上述功能需要兩類文件：(1)Gaussian的CIS、TDHF、TDDFT或TDA-DFT的輸出文件 (2)相應任務產生的.fch/fchk文件。Multiwfn的這個功能不限于Gaussian用戶使用，也支持ORCA、GAMESS-US等程序，詳見Multiwfn手冊3.21.A節的說明。
Gaussian的輸入文件對于CIS和TDDFT，分別寫成類似這樣
# B3LYP/6-31+G(d) TD(nstates=10) IOp(9/40=4)
# CIS(nstates=10)/6-31+G(d) IOp(9/40=4)
這里假設算10個激發態。Multiwfn計算時需要利用Gaussian輸出的激發態的組態系數，默認情況下只有絕對值大于0.1的系數會被輸出出來，而較小的都不輸出，這樣的話Multiwfn算出的結果將會不太準確。IOp(9/40=x)的含義是將系數大于10^-x的組態都輸出出來，因此IOp(9/40=4)會把系數絕對值大于0.0001的組態都輸出。x設的越大，輸出的組態越多，結果越準確，但是x太大的話Multiwfn的耗時也會非常大，通常x=3或4就夠了，精度足夠，計算速度也比較快。
在計算之后得到了.out/.log文件，還同時得到了.chk文件。用formchk將.chk轉換為.fch/fchk，這個文件里記錄了基函數的定義以及基態軌道信息，這是Multiwfn要利用的。
2 實例
下面例子涉及的文件都可以在http://www.shanxitv.org/attach/227/file.rar下載。
下面以苯酚為例進行說明怎么利用Multiwfn計算激發態之間的躍遷偶極矩，輸入文件如下，計算后得到phenol.out以及phenol.fch。
%chk=C:\phenol.chk
# b3lyp/6-311G* TD(nstates=5) IOp(9/40=4)
[空行]
b3lyp/6-311G** opted
[空行]
0 1
C 0.01810200 -1.86802400 0.00000000
C 1.23175000 -1.16732400 0.00000000
C 1.23175000 0.23407600 0.00000000
C 0.01810200 0.93477600 0.00000000
C -1.19554600 0.23407600 0.00000000
C -1.19554600 -1.16732400 0.00000000
H 0.01810200 -2.93802400 0.00000000
H 2.15839700 -1.70232400 0.00000000
H 2.15839700 0.76907600 0.00000000
H -2.12219300 0.76907600 0.00000000
H -2.12219300 -1.70232400 0.00000000
O 0.01810200 2.36477600 0.00000000
H -0.88699500 2.68477600 0.00000000
啟動Multiwfn，依次輸入以下命令，//后面的是注釋。
C:\phenol.fch //先載入電子激發任務產生的fch文件
18 //電子激發分析功能，包含多個子功能。這些功能無與倫比的強大，建議參看《Multiwfn支持的電子激發分析方法一覽》（http://www.shanxitv.org/437）
5 //計算激發態間的躍遷偶極矩
C:\phenol.out //電子激發任務的Gaussian輸出文件
此時屏幕上首先輸出了激發態的匯總信息
Exc.state# Exc.energy(eV) Multi. MO pairs Normalization
1 5.11440 1 1478 0.500002
2 5.93650 1 1697 0.500000
3 6.14320 1 919 0.499999
4 6.67690 1 828 0.499994
5 6.90240 1 1593 0.499994
其中N_pairs代表這個激發態在輸出文件中通過多少組態來表示，IOp(9/40=x)的x越大顯然N_pairs也就越大。Multi.是激發態的自旋多重度。Sum coeff.^2越接近理想值說明結果精度越高，對于閉殼層和開殼層情況理想值分別是0.5和1.0。如果偏離理想值比較大，則應該加大x來保證結果精度。
然后選1，激發態間的躍遷電偶極矩就輸出到了屏幕上。也可以選2輸出到當前目錄下的transdipmom.txt里。程序先給出了基態的偶極矩，然后給出基態到各個激發態的躍遷偶極矩、激發能和振子強度：
Ground state dipole moment in X,Y,Z: -0.554241 -0.157977 0.000000 a.u.
Transition dipole moment between ground state (0) and excited states (a.u.)
i j X Y Z Diff.(eV) Oscil.str
0 1 -0.4634428 -0.0387809 0.0000000 5.11440 0.02710
0 2 -0.0410220 0.4644091 0.0000000 5.93650 0.03161
0 3 0.0000000 0.0000000 0.0024911 6.14320 0.00000
0 4 0.0000000 0.0000000 -0.0776262 6.67690 0.00099
0 5 -1.3520332 -0.2422995 0.0000000 6.90240 0.31905
接下來，輸出的是激發態之間的躍遷偶極矩<i|-r|j>、能量差和振子強度。對于i=j的情況來說，其數值<i|-r|i>對應的是第i激發態的偶極矩的由電子貢獻的部分（注意這不等于這個態的偶極矩，因為還有另一部分，即原子核電荷對偶極矩的貢獻）。
Note: In below output the case of i=j corresponds to contribution of electron t
o dipole moment of excited state i
Transition electric dipole moment between excited states (a.u.):
i j X Y Z Diff.(eV) Oscil.str
1 1 -0.5515307 0.6638909 0.0000000 0.00000 0.00000
1 2 0.1289203 0.0072212 0.0000000 0.82210 0.00034
1 3 0.0000000 0.0000000 -0.0575947 1.02880 0.00008
1 4 0.0000000 0.0000000 0.0104765 1.56250 0.00000
1 5 0.0412045 1.1173067 0.0000000 1.78800 0.05476
2 2 -0.5370763 0.7195150 0.0000000 0.00000 0.00000
2 3 0.0000000 0.0000000 -0.0271832 0.20670 0.00000
2 4 0.0000000 0.0000000 0.0723899 0.74040 0.00010
2 5 0.7076947 -0.0614123 0.0000000 0.96590 0.01194
3 3 2.0119437 -3.9241936 0.0000000 0.00000 0.00000
3 4 0.3712230 1.1487737 0.0000000 0.53370 0.01906
3 5 0.0000000 0.0000000 -0.0330295 0.75920 0.00002
4 4 -0.2932424 2.0434979 0.0000000 0.00000 0.00000
4 5 0.0000000 0.0000000 0.0048938 0.22550 0.00000
5 5 -0.4634208 1.0586426 0.0000000 0.00000 0.00000
值得一提的是，如果在Gaussian中使用下面這樣的關鍵詞，在輸出文件末尾就會由L601模塊輸出第3個激發態的偶極矩
# b3lyp/6-311G* TD(root=3,nstates=5) density=rhoci
給出的結果為（debye）
X= 5.1144 Y= -9.9743 Z= 0.0000
換算成a.u.單位后，結果為X=2.0119 Y=-3.9238 Z=0.0000，可見和Multiwfn給出的2.0119437 -3.9241936 0.0000000很相符。注意這里用了rhoci關鍵詞，如果只寫density的話，那么Gaussian傳遞給L601模塊的激發態的密度是弛豫的密度。而當用rhoci時，傳遞的是非弛豫的密度，這才是和通過組態系數和基態軌道直接產生的密度直接對應的，也正是與Multiwfn輸出的是對應的。
如果你想得到各個態之間的躍遷磁偶極矩，進入功能18的子功能5后先選0，選Magnetic，然后再按上面做法操作即可。注意躍遷磁偶極矩定義為-i<i|r×▽|j>，而Multiwfn輸出的值對應的是<i|r×▽|j>部分，這和Gaussian輸出文件里給出的躍遷磁偶極矩用的形式相同。
Multiwfn也可以直接給出各個激發態的包含了原子核貢獻的偶極矩。也就是進入功能18的子功能5后選擇4，結果就會輸出到當前目錄下的dipmom.txt中，內容為：
Note: The electric dipole moments shown below include both nuclear charge and electronic contributions
Ground state electric dipole moment in X,Y,Z: -0.554241 -0.157977 0.000000 a.u.
Excited state electric dipole moments (a.u.):
State X Y Z exc.(eV) exc.(nm)
1 -0.551525 0.663891 0.000000 5.1144 242.42
2 -0.537071 0.719515 0.000000 5.9365 208.85
3 2.011949 -3.924194 0.000000 6.1432 201.82
4 -0.293237 2.043498 0.000000 6.6769 185.69
5 -0.463415 1.058643 0.000000 6.9024 179.62
PS：之所以這里把原子核對電子態的貢獻也考慮后，偶極矩數值和之前只考慮電子貢獻的時候一樣，那是因為Gaussian默認會把體系原點放置到核電荷中心位置，此時原子核對偶極矩貢獻恰好為0。