電子激發任務中軌道躍遷貢獻的計算
電子激發任務中軌道躍遷貢獻的計算
文/Sobereva @北京科音
First release: 2014-Apr-15 2020-Jan-28
做ZINDO、CIS、TDHF、TDDFT這類電子激發計算時，電子激發通常表示為大量分子軌道(MO)間的躍遷，Gaussian會輸出對應的躍遷系數，比如 Excited State 1: Singlet-A" 4.6966 eV 263.99 nm f=0.0003 <S**2>=0.000
14 -> 16 0.56407
14 -> 20 -0.21920
14 -> 25 -0.24192
14 -> 29 -0.17712
14 -> 33 -0.13568
這些數值確切來說是軌道躍遷對應的單激發組態函數的組態系數。對于CIS、TDHF、TDDFT、ZINDO方法，激發態波函數都是通過各種單激發組態函數線性組合來描述的。如果想了解細節，大家應當找量子化學書看看關于CIS方法的介紹，就自然不難理解其物理意義了。為了分析電子躍遷的本質特征，經常會考察每對MO躍遷對電子激發的貢獻。其實這是個極為簡單的問題，但是總看到一些人錯誤的討論，故在這里澄清一下計算方法。這里只是討論Gaussian程序的情況，對其它程序討論是類似的。
如果對Gaussian做TDDFT計算不了解，務必參看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）。
1 ZINDO和CIS計算
對于基態（確切來說是參考態）是閉殼層的計算，所有躍遷系數平方和相加等于0.5，即50%。之所以加和不為理應的1.0，是因為alpha軌道和beta軌道此時完全一致，alpha電子的躍遷和beta電子的躍遷的情況因此也是完全一致的，因此只輸出了alpha部分。例如上例的14號MO到16號MO的躍遷（14->16），其實分為兩種情況，即14A->16A和14B->16B，由于二者的躍遷系數是完全相同的，故只輸出了其一。因此，14->16對這個激發的貢獻就是0.56407^2/0.5*100%=63.6%，14->29的貢獻就是(-0.17712)^2/0.5*100%=6.3%。然而，上面的例子的5個躍遷系數的平方和相加結果為0.47453，離0.5還差得遠，這顯然不是有效位數的問題。于是有些人在計算貢獻時，比如計算14->16的貢獻時，居然這么算：0.56407^2/0.47453*100%=67.1%，這是明顯錯誤的！之所以加和不為0.5，是因為為了輸出信息簡潔，Gaussian默認只把系數絕對值大于0.1的系數輸出了出來。如果想讓更多的系數都輸出出來，可以用IOp(9/40=x)關鍵詞，它說明把絕對值大于10^-x的都輸出出來。如果你用比如IOp(9/40=3)把絕對值大于0.001的系數都輸出，那么平方和就非常接近于0.5了，但是此時屏幕上的信息量會略大，可能會有好幾十、幾百個系數。
對于基態是開殼層的計算，alpha軌道和beta軌道不再一致，故alpha和beta電子激發的情況在計算中就都明確考慮了，此時所有躍遷系數平方和為1.0。因此對于一對兒MO的躍遷，將其系數平方乘上100%就得到了它的貢獻了。
老有人問某些組態系數為負的時候是什么意思。根本無需關注正負號，算貢獻時求了平方之后都一樣，而且正負號也沒有可以直觀去理解的物理意義（分子軌道展開系數里有正有負還比較容易通過繪制軌道圖形來說明，然而組態函數本身是多電子基函數，其系數的正負號是相當抽象的，既沒法靠畫圖展現，也沒法給予明確的物理解釋）。
2 TDHF和TDDFT計算
TDHF和TDDFT計算和上一節的情況最大的不同在于不僅有激發，即“->”所表示的，還有去激發，即“<-”所表示的。“去激發”并沒有物理意義，這種稱呼僅僅是根據TDHF、TDDFT波函數的歸一化特點而給了這么個稱呼而已，不要去深究和過度解釋。這種TD形式的計算出現去激發是必然的，但通常去激發所占成份很小。例如Excited State 3: Singlet-B2 4.7665 eV 260.12 nm f=0.0044 <S**2>=0.000
35 -> 40 0.02865
36 -> 39 0.27686
36 -> 44 0.06122
36 -> 57 -0.01098
36 -> 62 -0.01239
37 -> 40 0.63762
38 -> 39 -0.10663
38 -> 50 -0.01114
36 <- 39 0.02135
36 <- 44 0.01165
37 <- 40 -0.01239
對于TDHF/TDDFT計算，對于基態為閉殼層和開殼層時還是分別歸一化為0.5和1.0，但注意，所有激發的系數平方和減去所有去激發的系數平方和得到的才是這個結果，也就是說去激發起到的是負貢獻。比如上例37->40的貢獻為(0.63762)^2/0.5*100%=81.31%，而36<-44的貢獻為-(0.01165)^2/0.5*100%=-0.027%，37<-40的貢獻為-(-0.01239)^2/0.5*100%=-0.031%。
值得一提的是，有一種TDA(Tamm-Dancoff approximation)形式的TDDFT，也叫TDA-DFT，它和一般的TDDFT的關系形同于CIS與TDHF的關系，也不存在去激發，這在Gaussian D.01及之后版本中可以用TDA關鍵詞來使用。TDA近似可以令TDDFT計算有一些加快，對激發能和振子強度有一定影響，精度往往比TDDFT還好一些，但是振子強度比TDDFT整體低估，算的轉子強度也不怎么好。更多信息看《亂談激發態的計算方法》（http://www.shanxitv.org/265）。
3 使用Multiwfn輸出貢獻值
如上所示，計算MO躍遷的貢獻很簡單，用計算器點幾下就出來了。如果懶得手算，也可以用著名的波函數分析程序Multiwfn來輸出貢獻值，Multiwfn中有一大批功能專門用來做電子激發分析。Multiwfn可以在其主頁http://www.shanxitv.org/multiwfn免費下載。如果不熟悉Multiwfn建議先參看《Multiwfn入門tips》（http://www.shanxitv.org/167）。下文對應的是3.6版及之后版本，更老版本與下面的操作過程不同。
對于ZINDO/CIS/TDHF/TDDFT任務，計算之后將chk文件轉化為fch文件，假設叫a.fch，相應的Gaussian輸出文件假設叫a.out，都放在了C:\下。啟動Multiwfn，依次輸入
C:\a.fch
18 //電子激發分析模塊
-1 //檢查、修改和導出組態系數
C:\a.out
3 //假設分析的是第3個電子激發態
然后屏幕上輸出了一些統計信息，包括激發系數的平方和以及去激發系數平方和的負值。之后馬上看到以下結果，是對電子激發貢獻最大的10個軌道躍遷信息。#后面那個數是這個軌道躍遷在所有輸出的軌道躍遷中的序號，不用管。
Some MO transitions sorted by absolute contributions:
# 1272 48 -> 49 Coeff.: 0.70422 Contri.: 99.1852%
# 2489 48 <- 49 Coeff.: -0.07085 Contri.: -1.0039%
# 1202 46 -> 49 Coeff.: -0.05577 Contri.: 0.6221%
# 748 33 -> 49 Coeff.: -0.03030 Contri.: 0.1836%
# 1274 48 -> 53 Coeff.: 0.02771 Contri.: 0.1536%
# 1273 48 -> 50 Coeff.: 0.02375 Contri.: 0.1128%
# 1203 46 -> 50 Coeff.: -0.02166 Contri.: 0.0938%
# 1224 47 -> 68 Coeff.: 0.02038 Contri.: 0.0831%
# 1183 45 -> 51 Coeff.: 0.01880 Contri.: 0.0707%
# 1218 47 -> 54 Coeff.: -0.01552 Contri.: 0.0482%
如果想把數據從屏幕上拷貝出來，見Multiwfn手冊5.4節的說明。如果只想把貢獻超過某個閾值的軌道躍遷輸出出來，可以在當前界面里用選項-2，然后輸入閾值。
使用Multiwfn還可以以更直觀的方式非常方便快捷地對一大批激發態輸出主要軌道躍遷貢獻，而且只需要提供Gaussian輸出文件就夠了而不需要fch，輸出信息例如
# 1 3.907 eV Multi= 1: H-4 -> L 81.9%, H-4 -> L+2 12.1%
# 2 4.062 eV Multi= 1: H -> L 86.0%, H-3 -> L 5.3%
# 3 4.417 eV Multi= 1: H-6 -> L 85.3%, H-6 -> L+2 11.9%
# 4 4.791 eV Multi= 1: H-2 -> L 54.5%, H -> L+1 27.6%, H-3 -> L+1 6.4%
# 5 4.887 eV Multi= 1: H -> L+3 57.3%, H-2 -> L 17.0%, H-1 -> L+2 8.8%, H
-1 -> L 8.0%
操作見《使用Multiwfn便利地查看所有激發態中的主要軌道躍遷貢獻》（http://www.shanxitv.org/529）。