使用ORCA在TDDFT下計算旋軌耦合矩陣元和繪制旋軌耦合校正的UV-Vis光譜
使用ORCA在TDDFT下計算旋軌耦合矩陣元和繪制旋軌耦合校正的UV-Vis光譜
文/Sobereva@北京科音
First release: 2019-Feb-10 Last update: 2021-Jul-10
1 前言
看本文前，強烈建議先看《使用Gaussian+PySOC在TDDFT下計算旋軌耦合矩陣元》（http://www.shanxitv.org/411）了解一些旋軌耦合矩陣元計算的基本知識，以及看《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）了解TDDFT計算的相關知識。TDDFT和旋軌耦合(Spin-orbit coupling, SOC)方面的理論性的內容在本文就不多提了，筆者假設讀者已經看了上面兩篇具備了相關常識。如果你不會裝ORCA，看《量子化學程序ORCA的安裝方法》（http://www.shanxitv.org/451）。
從ORCA 4.1開始，ORCA支持了TDDFT下閉殼層體系的旋軌耦合的計算，旋軌耦合算符對應的是Breit-Pauli哈密頓，可以做的相關的事包括：
·計算旋軌耦合矩陣元
·計算旋軌耦合對基態和激發態能量的影響
·計算考慮旋軌耦合時的振子強度和轉子強度，結合Multiwfn可以繪制考慮旋軌耦合后的光譜
本文就對ORCA的TDDFT的旋軌耦合計算進行基本介紹，其中最主要是講旋軌耦合矩陣元的計算和繪制旋軌耦合校正的UV-Vis光譜，最后也簡單談一下用ORCA的這種SOC-TDDFT計算有無可能給出靠譜的磷光發射速率的問題。
使用ORCA在TDDFT下計算旋軌耦合矩陣元相對于使用之前筆者介紹的Gaussian+PySOC的組合有下列好處：
(1)ORCA開RI的時候做TDDFT計算本身遠比Gaussian快得多，除很小體系外至少快一倍，而氣相的結果可以與Gaussian很好相符
(2)ORCA對學術用戶完全免費
(3)PySOC只能用較low的有效核電荷(Zeff)方式考慮旋軌耦合。而ORCA不僅支持Zeff，還支持更精確的旋軌耦合平均場(SOMF)方式考慮
(4)PySOC大部分是Fortran寫的，有一小部分莫名其妙地偏要用Python編寫，導致在Windows下運行不方便，而很多操作系統下由于Python版本原因導致PySOC運行不成功。而ORCA在所有操作系統下都能順利運行
(5)用ORCA計算比用Gaussian+PySOC的組合方便不少，而且還省得讓Gaussian往硬盤里保存很占地方的rwf文件
總的來說，ORCA算是目前在TDDFT下計算旋軌耦合矩陣元最理想的程序。
下面結合實例介紹ORCA的SOC-TDDFT計算。涉及的輸入輸出文件都可以在此處下載：http://www.shanxitv.org/attach/462/file.rar。本文用的是ORCA 4.1.1。
2 旋軌耦合矩陣元計算實例1：甲醛
首先用一個非常簡單的體系甲醛來說明怎么在ORCA里做考慮了旋軌耦合的TDDFT計算、怎么去理解輸出。
2021-Jul-10注：從ORCA 5.0開始，不要再像下文例子中那樣自己寫grid和gridx關鍵詞，用默認的即可。詳見《談談ORCA 5.0的新特性和改變》（http://www.shanxitv.org/604）。
下面這個輸入文件對甲醛在B3LYP/def-TZVP下做TDDFT計算。幾何結構已經在PBE0/def-TZVP級別下優化過。
! B3LYP/G TZVP miniprint tightSCF grid4 pal4
%tddft
nroots 5
dosoc true
tda false
printlevel 3
end
* xyz 0 1
C 0.00000000 0.00000000 -0.52513500
H 0.00000000 0.93987900 -1.11261300
H 0.00000000 -0.93987900 -1.11261300
O 0.00000000 0.00000000 0.67200400
*
之所以此例在B3LYP后面寫了個/G，是因為ORCA里默認的B3LYP與Gaussian的B3LYP不同，加/G代表我們這里想使用與Gaussian相同的B3LYP定義。TZVP代表用def-TZVP基組，和Gaussian里寫TZVP時用的基組相同。miniprint代表不輸出一堆亂七八糟沒太大用的信息（尤其是布居分析，又占地方又不好讀，而用筆者開發的Multiwfn隨時都可以基于ORCA產生的.molden文件做更豐富的布居分析，而且分析起來方便得多）。tightSCF代表用比默認更嚴的SCF收斂限，對于涉及到后HF、TDDFT等多組態方法的時候都建議加上。grid4代表用比默認更好的DFT積分格點，這有助于降低數值誤差。pal4代表用4核并行。%tddft...end段落代表做TDDFT計算，dosoc true代表TDDFT計算過程考慮旋軌耦合效應。ORCA默認做TDDFT時候是用TDA近似的，原理上不如完整形式的TDDFT精確（詳見《亂談激發態的計算方法》http://www.shanxitv.org/265），因此用tda false要求ORCA做完整形式的TDDFT計算。nroots 5代表計算5個激發態，當使用dosoc true的時候會自動用triplets true選項，因此會同時計算5個單重態激發態和5個三重態激發態。printlevel 3代表在SOC-TDDFT計算過程中比默認時輸出更多信息，其中有很多是很重要的，而在默認下居然沒有輸出！（這點在手冊只字未提，筆者好不容易才意外地試出來這個關鍵性的隱藏選項，簡直像發現游戲秘技一樣）。
計算開始后，程序按照以下流程計算和輸出
(1)做常規B3LYP/TZVP基態計算以得到參考態波函數
(2)做常規TD-B3LYP/TZVP計算，先計算5個單重態激發態，再計算5個三重態激發態
(3)計算旋軌耦合積分，構建所有態之間的旋軌耦合矩陣元，然后對此矩陣對角化得到本征值和本征矢
(4)輸出完整的旋軌耦合矩陣，實部和虛部分別輸出，每個三重態當做不同的三個子態考慮
(5)輸出單-三重態間旋軌耦合矩陣元
對于當前例子，其實算到這里就夠了，后面幾步的意義在之后的例子再細說。
(6)輸出旋軌耦合對基態的穩定化能，即SOC stabilization of the ground state:后面的
(7)輸出考慮了SOC后哈密頓的本征值，開頭是Eigenvalues of the SOC matrix:
(8)輸出考慮了SOC后哈密頓的本征矢，開頭是Eigenvectors of the SOC matrix:
(9)輸出不考慮SOC時的激發能、振子強度、轉子強度、躍遷偶極矩等躍遷信息，標題為TD-DFT-EXCITATION SPECTRA
(10)輸出考慮了SOC時的躍遷信息，標題為SOC CORRECTED TD-DFT/TDA-EXCITATION SPECTRA
此后還有些信息，就不是本文關注的了。
如果你用了printlevel 4，還會再額外輸出基函數間的SOC積分、MO之間的SOC積分等信息，一般沒必要輸出這些，除非自己想基于這些數據寫額外的程序。
具體來說，上述第(5)步依次輸出以下四個矩陣
(a)各個單重態與各個三重態間的SOC算符的x,y,z三個分量對應的旋軌耦合矩陣元
(b)各個單重態與各個三重態的三個子態間的總旋軌耦合矩陣元
(c)各個單重態與各個三重態間的LxSx、LySy、LzSz算符（即簡化的旋軌耦合算符）對應的旋軌耦合矩陣元
(d)各個三重態間的簡化的旋軌耦合算符對應的旋軌耦合矩陣元
我們一般最關心的信息是上面(b)部分，也正是PySOC程序最終輸出的，此例結果如下所示（如果不用printlevel 3，這個矩陣是不輸出的）
CALCULATED SOCME BETWEEN TRIPLETS AND SINGLETS
--------------------------------------------------------------------------------
Root <T|HSO|S> (Re, Im) cm-1
T S MS= 0 -1 +1
--------------------------------------------------------------------------------
1 0 ( 0.00 , -63.78) ( -0.00 , 0.00) ( -0.00 , -0.00)
1 1 ( 0.00 , 0.00) ( -0.00 , -0.00) ( -0.00 , 0.00)
1 2 ( 0.00 , 0.00) ( 5.62 , 0.00) ( 5.62 , -0.00)
1 3 ( 0.00 , -0.00) ( -3.80 , 0.00) ( -3.80 , -0.00)
1 4 ( 0.00 , -0.00) ( -0.00 , -38.41) ( -0.00 , 38.41)
1 5 ( 0.00 , 45.07) ( -0.00 , -0.00) ( -0.00 , 0.00)
2 0 ( 0.00 , -0.00) ( 0.00 , -0.00) ( 0.00 , 0.00)
2 1 ( 0.00 , 48.56) ( 0.00 , 0.00) ( 0.00 , -0.00)
2 2 ( 0.00 , 0.00) ( 0.00 , -0.30) ( 0.00 , 0.30)
2 3 ( 0.00 , 0.00) ( 0.00 , 0.19) ( 0.00 , -0.19)
...[略]
S=0對應基態，S=1~5和T=1~5是算出來的5個單重態激發態和5個三重態激發態。MS=0、-1、+1對應三重態的自旋磁量子數不同的三個子態。可見ORCA對實部Re和虛部Im是分別輸出的，而PySOC在輸出的時候輸出的是模，即實部和虛部平方和開根號。如果你想得到某個單重態和某個三重態之間總的旋軌耦合矩陣元，那就把這個單重態與三個子態的旋軌耦合矩陣元的模求平方再開根號，換句話說，就相當于把相應那一行的六個值都求平方加和后再開根號（例如PCCP,16,26184的式5就是用這個方式算的）。根據上面的輸出，我們可知諸如S1-T2之間的總旋軌耦合矩陣元是48.56 cm-1。（其實用(a)部分的對應的行的六個數值求平方和得到的結果與此也相同，此時可以不用寫printlevel 3）
默認情況下，ORCA是以SOMF方式考慮旋軌耦合效應的，在ORCA里也被叫做SOMF(1X)。SOMF全稱為Spin-orbit mean-field，它將SOC算符中難算的雙電子部分用平均場方式來近似描述，精度比較理想，雖然名義上比起更粗糙的Zeff處理耗時高，但根據筆者測試，即便對于大體系的SOC-TDDFT計算，SOMF也不算耗時有多高，而且總耗時的瓶頸也完全不在計算旋軌耦合積分上。因此若無特殊情況，一律用默認的SOMF即可。
如果你就是想在SOC-TDDFT計算中用Zeff方式考慮旋軌耦合的話，在上面的輸入文件中額外插入一行%rel SOCType 1 end即可（別寫到%tddft...end里頭去），這里%rel是設置旋軌耦合算符具體類型的段落。此時，在旋軌耦合部分計算開始處會看到
Operator type ... Effective nuclear charge
Effective nuclear charges used:
Atomtype H -> Zeff 1.0000
Atomtype C -> Zeff 3.6000
Atomtype O -> Zeff 5.6000
這里用的ORCA內置的有效核電荷來自Koseki的文章，和筆者之前在http://www.shanxitv.org/411中提到的一致。此時輸出的旋軌耦合矩陣元為
Root <T|HSO|S> (Re, Im) cm-1
T S MS= 0 -1 +1
--------------------------------------------------------------------------------
...[略]
2 0 ( 0.00 , -0.00) ( 0.00 , -0.00) ( 0.00 , 0.00)
2 1 ( 0.00 , 53.39) ( 0.00 , 0.00) ( 0.00 , -0.00)
2 2 ( 0.00 , 0.00) ( 0.00 , -0.33) ( 0.00 , 0.33)
...[略]
此時<S1|H_SO|T2>為53.39 cm-1，和SOMF的結果48.56 cm-1稍有差異。如http://www.shanxitv.org/411所示，之前用Gaussian+PySOC在相同級別下、結合Koseki有效核電荷時得到的這個值為49.81 cm-1，和此處的53.39 cm-1間的一點差異在于整個計算中涉及的各種零零碎碎的數值細節有所不同，而且PySOC計算時用的是笛卡爾型高斯函數，而當前是球諧型（ORCA里也只能用球諧型高斯函數）。
在ORCA里如果利用RI近似的話，對于大體系的基態和TDDFT計算都可以顯著加快。由于甲醛很小，RI的效果不明顯，甚至反倒令耗時更高，所以前例就沒開RI。而對大體系，強烈建議把前例的關鍵詞改為下面的關鍵詞來節約時間
! B3LYP/G TZVP RIJCOSX def2/J miniprint tightSCF grid4 gridx4 pal4
RIJCOSX代表用RI-J加速庫侖部分的計算，用COSX方法快速計算交換部分，gridx4代表用比默認更好的做COSX部分的積分格點以保證數值精度損失小。def2/J代表在RI-J部分使用def2/J輔助基組（def2/J輔助基組給def-系列基組用完全沒問題）。對當前體系，開不開RIJCOSX對算出來的旋軌耦合矩陣元的絕對值影響不超過0.1 cm-1，所以SOC-TDDFT計算中可以放心用RIJCOSX來加速。更多RI相關的信息在《大體系弱相互作用計算的解決之道》（http://www.shanxitv.org/214）一文有介紹。
3 旋軌耦合矩陣元計算實例2：溶劑中的Os配合物
在前述的《使用Gaussian+PySOC在TDDFT下計算旋軌耦合矩陣元》一文中用Gaussian+PySOC計算過下面這個Os配合物體系，基態是單重態
由于有重金屬，因此我們一般都會想到用贗勢基組，比同等質量的全電子基組又省時間又能順帶體現標量相對論效應。如果用Zeff方式考慮旋軌耦合的話，是可以使用贗勢的。雖然如筆者的PySOC那篇博文所述，Koseki提出過適合Os等重元素的有效核電荷，結果也不差，但是ORCA 4.1.1里沒有內置重元素的Koseki有效核電荷，雖然手冊里說可以由用戶自行提供，但筆者發現ORCA 4.1.1存在bug，按照手冊里的方式自行提供了也會被程序無視。這點估計以后版本會修正，筆者已經在ORCA論壇反饋過了(https://orcaforum.kofo.mpg.de/viewtopic.php?f=8&t=4478)。因此，對此體系，時下我們只能用SOMF方式考慮旋軌耦合（ORCA還支持其幾種變體，大同小異），此時只能用全電子基組，用贗勢基組+贗勢的話雖然不報錯，但結果是明顯不對的。
這次我們用以下關鍵詞通過DKH2哈密頓來做全電子標量相對論計算，旋軌耦合只體現在TDDFT部分。完整的輸入文件見本文的文件包。此任務用的幾何結構和之前PySOC博文里相同。
! B3LYP/G DKH2 DKH-def2-TZVP SARC/J RIJCOSX cpcm(CH2Cl2) tightSCF miniprint grid4 gridx4
%basis
NewGTO Os "SARC-DKH-TZVP" end
end
%pal nprocs 36 end
%maxcore 2500
%tddft nroots=10 TDA false dosoc true printlevel 3 end
* xyz 0 1
...[坐標]
end
其中，DKH2代表用流行的DKH2哈密頓做標量相對論計算，基組使用DKH-def2-TZVP，是專為DKH計算設計的重收縮版本的def2-TZVP。由于def2-TZVP對于Os是贗勢基組，所以DKH-def2-TZVP對Os必然是沒有定義的，故此例對Os元素使用SARC-DKH-TZVP，此基組是對Xe之后元素有定義的適合DKH2計算的全電子基組，與DKH-def2-TZVP搭配很適合。雖然SARC/J名字看上去好像只是給SARC-DKH-系列基組用的輔助基組，但在ORCA里它也同時具有def2/J DecontractAuxJ（對def2/J輔助基組去收縮）的含義，因此無論是SARC-DKH系列（也包括SARC2-DKH）還是DKH-def2-系列基組都可以通過寫SARC/J來提供適合的輔助基組。當前計算使用CPCM模型表現二氯甲烷CH2Cl2溶劑環境。由于當前計算較耗時，故用36個核計算，每個核內存分配上限約為2500MB。
當前這個任務，在http://bbs.keinsci.com/thread-6310-1-1.html一文提到的筆者買的E5-2696 v3雙路36核機子下只花了12分鐘就跑完了。旋軌耦合矩陣元部分輸出如下
Root <T|HSO|S> (Re, Im) cm-1
T S MS= 0 -1 +1
--------------------------------------------------------------------------------
1 0 ( 0.00 , 0.01) ( -59.62 , 7.55) ( -59.62 , -7.55)
1 1 ( 0.00 , 89.52) ( 0.03 , 0.00) ( 0.03 , -0.00)
...[略]
1 9 ( 0.00 , 3.58) ( 0.04 , 0.05) ( 0.04 , -0.05)
1 10 ( 0.00 , 0.03) ( 1.48 , -40.59) ( 1.48 , 40.59)
2 0 ( 0.00 , -198.63) ( -0.86 , 0.23) ( -0.86 , -0.23)
2 1 ( 0.00 , -1.30) (-544.99 , 16.42) (-544.99 , -16.42)
2 2 ( 0.00 , 0.94) ( -70.51 , 665.08) ( -70.51 , -665.08)
...[略]
比如<S1|H_SO|T2>的總大小為sqrt((-1.30)**2+(-544.99)**2+16.42**2+(-544.99)**2+(-16.42)**2)=771.1 cm-1。之前PySOC那篇博文里，用Gaussian+PySOC通過TD-B3LYP/def2-TZVP結合IEFPCM溶劑模型得到的結果為<S1|H_SO|T2_0>=0.69 cm-1、<S1|H_SO|T2_+1>=<S1|H_SO|T2_-1>=419.63 cm-1，總大小593.4 cm-1，和當前算的結果數量級差不多。由于當前計算明顯更為嚴格，又是用全電子標量相對論又是用更好的SOMF方式考慮旋軌耦合，特別是Zeff方法對于越重的元素往往誤差越大，因此應當認為本文的結果更合理。
2021-Jul-10注：以下關于溶劑部分的內容對ORCA 5.0及以后版本不適用，程序不再輸出下面這些內容，而且溶劑下計算結果也沒有4.x時候存在的以下提及的問題。
關于ORCA在溶劑模型下輸出的TDDFT激發能有一些特殊問題值得提及。在輸出完三重態激發能之后，有這么一段（這里State 0是指的第1個三重態激發態）
State Shift(Eh) Shift(eV) Shift(cm**-1) Shift(nm) ERPA(eV) ERPA+SHIFT(eV)
-------------------------------------------------------------------
0: -0.0006565 -0.018 -144.1 3.5 2.534 2.516
1: -0.0120876 -0.329 -2652.9 28.4 3.954 3.626
2: -0.0100690 -0.274 -2209.9 22.6 4.012 3.738
3: -0.0011292 -0.031 -247.8 2.1 4.287 4.256
4: -0.0024908 -0.068 -546.7 4.6 4.330 4.262
5: -0.0021236 -0.058 -466.1 3.8 4.365 4.307
6: -0.0297932 -0.811 -6538.8 61.7 4.461 3.650
7: -0.0083835 -0.228 -1840.0 14.9 4.475 4.247
8: -0.0007398 -0.020 -162.4 1.2 4.486 4.466
9: -0.0003368 -0.009 -73.9 0.5 4.569 4.560
我在《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）中說過，溶劑對電子激發的響應有快部分和慢部分之分，兩部分都會影響激發能。上面的ERPA下面的數值對應的只是只考慮了慢部分的情況，相當于是用溶劑模型下經過SCF得到的基態軌道以常規方式算的TDDFT激發能，ORCA是根據這個激發能來對激發態進行排序的。Shift是溶劑的快部分對激發能的影響，體現的是對電子激發的瞬時響應，二者的加和ERPA+SHIFT也正是在這個表格之前輸出的各個激發態的能量，它和ERPA的順序可能存在不同。而Gaussian則是按照最終的（快、慢部分都考慮后的）激發能對激發態來排序的。因此由于這個差異，Gaussian與ORCA給出的激發態順序可能不同。
另外，Gaussian與ORCA雖然在氣相下算的TDDFT激發能相符很好，但在溶劑模型下則有很大出入（對于電子激發，ORCA只能用線性響應溶劑模型，而這里作為對比的對象是Gaussian也用線性響應溶劑模型的情況）。這倆程序都用CPCM時，由于用的默認原子半徑不同等細節因素存在差異，結果肯定是對應不好的，這暫且不提；哪怕這兩個程序都用SMD溶劑模型，讓溶劑方面的細節盡可能一致（仍不完全相同，因為ORCA的SMD的極性部分不是用的SMD原文里的IEFPCM而是更low的CPCM形式），這倆程序給出的某些態的TDDFT激發能也是存在明顯不同的（已經通過組態系數可以確認對比的是相同的態），相差往往達到零點幾eV的程度。根據筆者的感覺，ORCA算出來的溶劑的快部分是不太可靠的，即SHIFT數值往往太大，如上可見有的居然能達到-0.8 eV的程度，更何況當前用的二氯甲烷極性很低，不可能會產生這么大影響才對。而由于Gaussian在溶劑模型方面有資歷深厚的專家Barone等人在做，發表過很多關于溶劑效應對電子激發問題的理論文章，因此我認為Gaussian的結果明顯更可信。而且，如果不考慮快部分，即直接用ORCA輸出的ERPA的結果與Gaussian給出的結果相比，則能基本相符，進一步體現ORCA給出的快部分應當是有嚴重問題的。
根據以上討論，對于溶劑模型下ORCA給出的SOC-TDDFT旋軌耦合矩陣元，就用直接輸出的態的編號討論就行了，而不建議根據實際輸出的激發能（也等同于上表中的ERPA+SHIFT）再對激發態手動重新排序。如上可見，重新排序的話有的態的序號會有明顯改變，比如從上面的表里看到ORCA給出的第7個三重態激發態的最終的激發能是3.650 eV，如果按照最終激發能來重新排序的話，它就成了T3而非T7了。
4 旋軌耦合校正后的UV-Vis光譜的繪制實例：Ir(bzq)2(acac)配合物
普通TDDFT計算時，由于自旋禁阻，單重態基態與三重態激發態之間的振子強度精確為0，而考慮旋軌耦合后，就不再如此。ORCA的SOC-TDDFT任務會輸出在TDDFT過程中考慮旋軌耦合后基態到各個激發態的激發能和振子強度。具體是這么做的：先構建基態以及所有激發態之間的旋軌耦合矩陣元并加到原本的哈密頓矩陣上，每個三重態當做不同的三個子態來考慮。然后對這個矩陣做對角化，這等同于以原先的單、三重態電子態波函數作為基在考慮了旋軌耦合算符的哈密頓下做變分。這樣得到的每個本征矢就相當于是考慮了旋軌耦合后的每個電子態的波函數，形式上是原先各個單重態電子態和三重態電子態的線性組合，組合系數就是這個本征矢的相應元素，而本征值就對應于考慮旋軌耦合后各個電子態的能量。此時，基態與激發態之間就不再自旋禁阻了，因為基態和激發態現在都同時具有單、三重態的成份，不再是S^2算符的本征函數，因此去計算基態與激發態之間的躍遷偶極矩時，其中不僅包含原本因為自旋禁阻而精確為0的<S|-r|T>項，還包括可以不為零的<S|-r|S>和<T|-r|T>項。另外，新產生的電子態的能量顯然和之前電子態的能量不同，無論對于基態還是激發態都是如此。新的基態與原先基態之間的能量變化就是ORCA輸出的SOC stabilization of the ground state。PS：由于ORCA的SOC-TDDFT在撰文時還沒有相關文章發表，所以上敘述只是來自筆者個人的判斷。
看過不少磷光速率計算的讀者應該也知道，還有一種常見的得到單-三重態間振子強度的方法是使用一階微擾理論，也就是將旋軌耦合算符作為微擾項來修正原有的波函數，經由S-T間的SOC矩陣元使三重態摻入單重態，也讓單重態摻入三重態，這樣最終得到的S-T混合態之間的躍遷偶極矩就可以不為零了，因此振子強度也就可以不為零了。從原理上說，ORCA用的變分方式的處理應當更為嚴格，因為相當于可以考慮更高階的混合效應，而且可以順帶直接給出考慮SOC后的電子態的能量，但需要額外付出的代價是還得計算三重態激發態之間的旋軌耦合矩陣元，這占了整個旋軌耦合矩陣的大部分。
基于考慮旋軌耦合后產生的新的態的激發能和振子強度、轉子強度，原理上就可以照常繪制UV-Vis光譜和ECD光譜了。但遺憾的是，目前的ORCA 4.1.1在SOC-TDDFT計算過程中給出的考慮旋軌耦合后的激發態的轉子強度是錯的。在考慮SOC前輸出的有正有負，是合理的，但是考慮SOC后全都成正的了，明顯不對，因此繪制出來的ECD光譜也無意義。這個問題我已經在ORCA論壇進行了反饋（https://orcaforum.kofo.mpg.de/viewtopic.php?f=8&t=4479），Neese說這應該是bug，但愿在未來版本中會修正。
下面筆者以下圖這個Ir(bzq)2(acac)體系作為例子演示繪制旋軌耦合校正后的UV-Vis光譜。
這里使用下面的關鍵詞算對Ir(bzq)2(acac)進行計算，體系基態是單重態。
! B3LYP/G DKH2 DKH-def2-TZVP(-f) SARC/J RIJCOSX tightSCF miniprint grid4 gridx4
%pal nprocs 36 end
%maxcore 2500
%basis
NewGTO Ir "SARC-DKH-TZVP" end
end
%tddft nroots=100 TDA false dosoc true printlevel 3 end
# BP86/SDD/TZVP opted
* xyz 0 1
...[坐標]
end
其中#開頭的行是注釋行，后面的內容用于提示自己之前這個結構是在BP86結合SDD和def-TZVP基組優化后的。DKH-def2-TZVP(-f)代表把DKH-def2-TZVP基組中的f極化函數砍掉以節約時間。如《談談贗勢基組的選用》（http://www.shanxitv.org/373）里提過的，配體的地位沒有處于體系中央的過渡金屬高，再加上普通泛函對高角動量函數要求比較低，因此用-f的版本不至于令誤差比不帶-f的版本增加太多。順帶一提，ORCA的SCF收斂做得不如Gaussian好，當前這樣的關鍵詞算某些過渡金屬配合物容易遇到SCF不收斂，碰見這種情況可以嘗試加上slowconv或veryslowconv關鍵詞，此時會用對一般情況來說更慢但是更穩妥、收斂成功幾率更高的收斂設定和算法，實測挺奏效。當前例子我們要算的并不是S與T不同子態間的旋軌耦合矩陣元，因此printlevel 3其實可以去掉，此時由于可以少輸出許多矩陣，輸出文件會小很多。當前這個任務用2*E5-2696 v3的機子在36核并行時4h13m算完。
Multiwfn是繪制各種類型電子和振動光譜的又方便又強大的程序，可在http://www.shanxitv.org/multiwfn下載。相關信息看《Multiwfn入門tips》（http://www.shanxitv.org/167）。在《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》（http://www.shanxitv.org/224）中對各種光譜的計算和繪制原理以及在Multiwfn中的操作都有詳細說明。2019年1月22日之后更新的Multiwfn支持基于ORCA 4.x的SOC-TDDFT計算的輸出文件繪制UV-Vis和ECD譜，2021年7月10日及之后更新的Multiwfn開始兼容ORCA 5.0版的輸出。
啟動Multiwfn，輸入上面計算產生的輸出文件的路徑，然后依次輸入
11 //繪制光譜
3 //UV-Vis
y //載入考慮了SOC之后的電子激發數據。如果選n則載入的是考慮SOC之前的
0 //繪制光譜，顯示到屏幕上
此時屏幕會立刻看到圖像。為了美觀，我們調整一下坐標軸范圍和標簽間隔。關閉窗口，繼續輸入
3 //設置橫坐標
150,550,50 //范圍為150~550nm，間隔50nm
4 //設置左側坐標軸
0,60000,5000
y //讓右側坐標軸同步調節
此時再選0，看到的圖像如下
我們也可以將考慮SOC和不考慮時候的圖繪制到一起，便于考察SOC對光譜的影響。這用Multiwfn非常容易實現，怎么用Multiwfn繪制多個情況下的光譜在《使用Multiwfn繪制構象權重平均的光譜》（http://www.shanxitv.org/383）一文末尾明確說過。創建一個名為multiple.txt的文件，內容如下（.out路徑寫的是實際路徑），第一部分是文件名，之后是要顯示在圖上的圖例。
E:\Ir_bzq2_acac.out with SOC
E:\Ir_bzq2_acac.out without SOC
然后啟動Multiwfn，載入multiple.txt，重復一遍上面繪制UV-Vis光譜的過程。與上例唯一不同的是載入文件的時候會載入兩次Ir_bzq2_acac.out，第一次問你是否載入SOC校正后的光譜信息時輸入y，第二次問你的時候輸入n（注意，由于考慮SOC后每個三重態作為三個不同的態考慮，態數比不考慮SOC時要多，而Multiwfn里同時繪制多個光譜時要求態數最多的情況作為第一項出現，因此multiple.txt里必須第一條對應的是SOC的情況）。繪制的結果如下：
由圖可見，SOC效應使得UV-Vis光譜整體藍移幾十nm。按理說，肯定是考慮SOC之后的光譜更合理，但對于這點，筆者持保留態度。因為對于過渡金屬配合物的計算，只要計算級別合理，其實不考慮SOC就已經可以算得很不錯了，而ORCA這樣考慮SOC后，對結果影響卻如此之大，所以可能結果變得更差。筆者手頭沒有此體系的實驗光譜，而且即便有，也有可能因為與其它因素帶來的誤差的僥幸抵消而妨礙判斷，所以考慮SOC后光譜是否更合理這點筆者在此文就暫不深究了。讀者有興趣的話可以考慮和免費的專做相對論計算的Dirac程序做的二分量TDDFT計算得到的光譜對照一下，以看看ORCA這種SOC的考慮是否對改進TDDFT光譜確實有益。
這里再提一下當前體系考慮SOC后的激發能具體怎么來的。ORCA計算三重態時輸出的原始的T1激發能是
STATE 1: E= 0.082507 au 2.245 eV 18108.3 cm**-1
155a -> 159a : 0.047168
156a -> 159a : 0.030710
157a -> 158a : 0.873592
在SOC矩陣構造完畢并對角化考慮了SOC的哈密頓矩陣后，輸出了如下信息
SOC stabilization of the ground state: -1692.8499 cm-1
Eigenvalues of the SOC matrix:
State: cm-1 eV
0: 0.00 0.0000
1: 19166.96 2.3764
2: 19201.28 2.3807
3: 19269.86 2.3892
4: 19499.20 2.4176
...略
可見此處輸出的State 1,2,3的激發能比較接近，這三個態實際上就是T1態分裂出的三個子態。原先T1的激發能是18108.3 cm-1左右，考慮SOC后分裂出的三個態激發能是19200 cm-1左右，明顯高了不少。但是這還不是最終的激發能，因為前面提示了，SOC還導致基態能量下降了-1692.8499 cm-1，因此T1最低的那個子態最終的激發能是19166.96+1692.85=20859.8 cm-1，這和我們在靠后部分看到的輸出考慮SOC后的激發能、振子強度和躍遷偶極矩的段落中的信息是一致的
-----------------------------------------------------------------------------
SOC CORRECTED ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS*
-----------------------------------------------------------------------------
State Energy Wavelength fosc T2 TX TY TZ
(cm-1) (nm) (au**2) (au) (au) (au)
-----------------------------------------------------------------------------
1 20859.8 479.4 0.000051143 0.00081 0.02788 0.00000 0.00548
...略
正是因為此例中SOC既提升了激發態能量又降低了基態能量，因此使得光譜整體出現了明顯藍移。
關于算的態數這里多說幾句。由于當前體系較大，對于繪制電子光譜目的，計算的態數不能太少，像之前計算Os配合物那樣只算10個態是絕對不夠的，最少也得算50個態。下圖是Multiwfn繪制的不考慮旋軌耦合時計算50、100、150個態的對比圖（為節約耗時，配體用的是DKH-def2-SV(P)），可見至少算50個態的時候300nm以上的曲線才能保證真實、完整。
ORCA考慮SOC對光譜的影響靠的是讓不考慮SOC時得到的單重態和三重態混合來實現的，只有當算的態數比較多時，也相當于令“基”完備性較高時，考慮SOC后的電子態才能被較好描述，這樣得到的考慮SOC后的激發態的振子強度和激發能才可能比較合理。在繪制SOC下的UV-Vis光譜時，算的態數應當比起不考慮SOC時更多。下圖是計算50、100、150個態的時候在考慮SOC后的對比圖
我們就關注300nm右側的區域就夠了。由圖可見，算50個態和算100個態的時候光譜曲線是有明顯可察覺的差異的，而計算100個態和計算150個態的光譜則看不出太大差異。如果不考慮SOC，如之前的圖可見，算50、100還是150個態在300nm以上區域是肉眼完全看不出差異的。這個現象說明，在考慮SOC時，只計算50個態還不足以令光譜達到收斂，計算100個態的時候則基本達到收斂了。像此例類似的體系，做SOC-TDDFT計算光譜時，我都建議計算不少于100個態以求穩妥。
隨著計算的態數增加，SOC-TDDFT耗時也會迅速增加，所以也不能盲目算得太多。當算的態數很多時，最最耗時的步驟就是構建完整的SOC矩陣。對本例的體系，在B3LYP下對配體用DKH-def2-SV(P)、Ir用SARC-DKH-TZVP，在前述36核機子下以SOC-TDDFT方式算50、100、150個態的耗時分別為27m、79m、169m。值得一提的是，雖然ORCA默認用的SOMF(1X)方式考慮旋軌耦合比Zeff要貴，但是對于計算的態數較多時，其實總耗時都差不多，因為瓶頸完全不在于計算旋軌耦合積分的部分。
5 基于ORCA的SOC-TDDFT輸出信息計算磷光發射速率靠譜么？
磷光發射對應T1->S0，因此，如果找出T1主要對應的在考慮旋軌耦合后產生的三個態的振子強度和激發能，我們還可以計算磷光發射速率。在《Gaussian中用TDDFT計算激發態和吸收、熒光、磷光光譜的方法》（http://www.shanxitv.org/314）中說過自發輻射速率怎么算。對于T1的一個子態，只要把其激發能、振子強度都代進去，就可以得到從這個子態發磷光的速率。要算總的磷光速率，圖省事可以直接把三個子態的發射速率取平均，但原理上更嚴格的方法是將三個子態的發射速率取權重平均，權重可以通過Boltzmann分布獲得，見《根據Boltzmann分布計算分子不同構象所占比例》（http://www.shanxitv.org/165），只不過此時公式里的各個構象相對于最穩定構象的自由能差應該改為各個子態相對于能量最低子態的激發能差。
基于ORCA的SOC-TDDFT輸出的信息按上述方式算磷光發射速率的結果靠譜么？我們還用上一節Ir(bzq)2(acac)那個例子，看看算出的結果和實驗值相符如何。我們從輸出文件中找到下面部分
SOC CORRECTED ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS*
-----------------------------------------------------------------------------
State Energy Wavelength fosc T2 TX TY TZ
(cm-1) (nm) (au**2) (au) (au) (au)
-----------------------------------------------------------------------------
1 20859.8 479.4 0.000051143 0.00081 0.02788 0.00000 0.00548
2 20894.1 478.6 0.000053714 0.00085 0.00007 0.02909 0.00002
3 20962.7 477.0 0.002733492 0.04293 0.20278 0.00003 0.04251
4 21192.1 471.9 0.000006626 0.00010 0.00001 0.01015 0.00000
5 21260.4 470.4 0.002483201 0.03845 0.19479 0.00011 0.02259
6 21273.4 470.1 0.000053372 0.00083 0.00093 0.02872 0.00012
...略
由于T1的能量低于其它激發態，而且從以上信息看，前三個態的激發能很接近（旋軌耦合造成的能量分裂程度一般不會很大），因此可以初步認為前三個態就是T1分裂出的子態。但由于S1可能與T1能量較為接近，所以光從能量上判斷哪三個態對應原始的T1也不是絕對穩妥的。我們可以考察下考慮SOC后的哈密頓矩陣對角化產生的本征矢的構成來進一步判斷。從輸出文件中可以找到下面信息
Eigenvectors of the SOC matrix:
E(cm-1) Weight Real Imag : Root Spin Ms
STATE 0: 0.00
0.95663 -0.97807 -0.00121 : 0 0 0
STATE 1: 19166.96
0.45282 0.67290 0.00545 : 1 1 0
0.01878 -0.13705 -0.00111 : 4 1 0
0.19397 0.44040 0.00356 : 1 1 -1
0.03062 -0.00138 0.17499 : 3 1 -1
0.01131 0.00079 -0.10635 : 5 1 -1
0.19397 -0.44040 -0.00357 : 1 1 1
0.03062 -0.00145 0.17499 : 3 1 1
0.01131 0.00094 -0.10635 : 5 1 1
STATE 2: 19201.28
0.03764 0.00211 0.19399 : 3 1 0
0.40524 0.63655 -0.00725 : 1 1 -1
0.01383 -0.00127 -0.11758 : 2 1 -1
0.02012 0.00149 0.14183 : 3 1 -1
0.40524 0.63655 -0.00661 : 1 1 1
0.01383 0.00129 0.11758 : 2 1 1
0.02012 -0.00159 -0.14183 : 3 1 1
STATE 3: 19269.86
0.01172 0.00207 -0.10824 : 2 0 0
0.03248 0.00345 -0.18018 : 3 0 0
0.39122 0.62536 0.01197 : 1 1 0
0.23609 -0.48581 -0.00888 : 1 1 -1
0.01004 0.00202 -0.10017 : 2 1 -1
0.01122 0.10591 0.00201 : 4 1 -1
0.23609 0.48579 0.00972 : 1 1 1
0.01004 0.00182 -0.10018 : 2 1 1
0.01122 -0.10590 -0.00205 : 4 1 1
...略（注意給出的能量是沒有加上SOC導致的基態穩定化能的）
上面的信息展現了考慮SOC后新產生的電子態是怎么由原先的各個單重態和三重態混合產生的，Weight下面是貢獻，總和為1（這里只輸出了>0.01的）。Root下面是不考慮SOC時電子態的序號，0是基態，Spin為0和1分別對應單重態和三重態。
由以上信息可見，考慮SOC后能量最低的三個激發態幾乎完全由三重態構成，而且原先的T1（Root=1,Spin=1的項）起到的貢獻是最主要的。毫無疑問考慮SOC后得到的激發能最低的三個態一起對應了T1，用這三個態的信息，從原理上就可以計算磷光發射速率。
筆者做了一個Excel表格，把T1的三個子態的激發能和振子強度輸進去，就可以計算各個子態的發射速率、壽命以及Boltzmann分布比率，在表格最下面給出算術平均和權重平均后的發射速率。對當前情況，結果如下所示
在J. Phys. Chem. C, 117, 25714 (2013)中給出了此體系的實驗磷光速率，是0.6E5 /s，明顯比我們算出來的209869.5 /s低得多，連數量級都不對應！雖說實驗值是在2-MeTHF下測的，但溶劑效應也不會導致這么大差異，比如在CH2Cl2下算的此體系的權重平均的磷光速率是182851.2 /s，也沒比當前氣相算出來的小太多。當前雖然我們用的結構是基態優化的，而計算磷光理應使用對T1態優化的結構，但結構的這點差異也不會導致算出來的磷光速率產生數量級的差異（而且JPCC這篇文章提到有人發現用基態優化的結構算磷光反倒更準，且認為由于T1極小點的勢阱很淺，因此實際磷光發射的結構介于S0和T1結構之間，故這篇文章都用S0結構算的）。總之，基于ORCA的SOC-TDDFT輸出的數據，用常規計算磷光速率的方式，結果是很不可靠的（雖然我也發現對于個別配合物體系如Ir(ppy)3的結果基本合理，但終究穩健程度低，不敢用）。而且為了令SOC效應考慮充分，還不得不算很多態，起碼100個（最好還做個結果隨態數增加的收斂性測試），此時耗時也較高，所以大家就別指望利用ORCA以上述方式算磷光速率了。
順帶一提，Dalton算是目前最理想的算磷光速率的程序，而且免費，它算磷光速率是基于響應函數理論算的，等同于算無窮數目激發態，因此原理上更高級。之前筆者對此體系，用Dalton2016，以Zeff方式考慮旋軌耦合，在TD-B3LYP/6-31G*/SDD下算的結果是6.2E4 /s，和實驗相符奇好（當然很大程度上是巧合），用2*E5-2696v3 36核機子只花了一個小時。以后有時間我會專門撰文簡要介紹怎么用Dalton實現磷光速率的計算，在北京科音高級量子化學培訓班里會更深入地講。