用ORCA计算旋轨耦合矩阵元
一、理论简介
旋轨耦合的理论涉及相对论量子力学，此处仅以定性的形式粗略介绍相关背景。相对论效应是指进行电子结构计算时Dirac方程与Schrödinger方程这两个理论模型之间的差别。Dirac于1928年建立了电子运动的相对论方程——Dirac方程，但是Dirac本人却认为在化学问题中，价电子受内层电子的屏蔽，其运动速度比光速小很多，相对论效应很小。但在后来的研究中，人们逐渐认识到相对论效应的重要性。自旋-轨道耦合（spin-orbit coupling, SOC），简称旋轨耦合，是一种相对论效应，指电子的自旋和轨道运动之间的相互作用。在非相对论量子力学中，自旋态改变的跃迁是禁阻的；当考虑旋轨耦合时，这样的过程才能发生，比如系间窜越（intersystem crossing, ISC）、磷光发射等过程。在《用高斯计算磷光发射能》一文中我们提到，用TD-DFT直接计算T1和S0之间的跃迁，得到的振子强度始终为0，只有当考虑旋轨耦合后，振子强度才不为0。
在常见的量子化学程序中，能在TD-DFT级别下进行旋轨耦合计算的程序有ADF、BDF、Dalton和ORCA等。此外，用Gaussian结合PySOC程序也可以实现SOC的计算。本文以ORCA 5.0.1为例介绍旋轨耦合矩阵元的计算。
二、计算实例1
体系一来自文献Chem. Sci., 2019, 10, 5031–5038一文中的mBrPhCz分子，分子结构如下：
这类材料有很强的磷光发射，并且作者提出可以直接通过S0→T1的吸收提高三重态的布居。我们便尝试计算这个过程的SOC矩阵元。
(1) 优化基态结构：
%pal nproc 24 end
%maxcore 3000
! B3LYP TZVP opt freq tightscf
*xyz 0 1
Br 5.45600000 1.42100000 0.77200000
C 6.70300000 2.00100000 2.08700000
C 10.80500000 0.75600000 1.08300000
H 10.22200000 1.41500000 0.72700000
N 9.94400000 0.90100000 3.45000000
C 8.75600000 1.65400000 3.25600000
C 8.45800000 2.73000000 4.08600000
H 9.05800000 2.97400000 4.78100000
C 7.28500000 3.44500000 3.89800000
H 7.08300000 4.17800000 4.46700000
C 10.75300000 0.39100000 2.42600000
C 6.40100000 3.09800000 2.87800000
H 5.60900000 3.60100000 2.73000000
C 9.77000000 0.41500000 5.92600000
H 9.06900000 1.03100000 6.10400000
C 10.30900000 0.28000000 4.64700000
C 11.74800000 0.11500000 0.28900000
C 7.86600000 1.26700000 2.25100000
H 8.05300000 0.52000000 1.69500000
C 11.63800000 -0.56500000 2.97500000
C 10.29200000 -0.38500000 6.92800000
C 11.34500000 -0.64200000 4.38600000
H 11.80300000 0.33900000 -0.63300000
C 12.62100000 -0.85400000 0.81000000
C 12.56500000 -1.20000000 2.15400000
H 9.94100000 -0.30600000 7.80700000
C 11.32100000 -1.30600000 6.69000000
C 11.86000000 -1.43900000 5.41100000
H 13.25300000 -1.27500000 0.24000000
H 13.15000000 -1.85900000 2.50800000
H 11.65400000 -1.84000000 7.40200000
H 12.55900000 -2.05800000 5.24000000
*(2) S0→Tn的旋轨耦合矩阵元计算：
%pal nproc 24 end
%maxcore 3000
! B3LYP def2-TZVP tightscf
%tddft
nroots 6
dosoc true
tda false
end
*xyzfile 0 1 mbpc_gsopt.xyz最后一行表示直接从上一步优化生成的xyz文件中读取结构。在%tddft区域内设置TD-DFT计算的相关参数：nroots表示计算6个激发态；dosoc设为true表示计算旋轨耦合，并且此时会自动计算三重态；ORCA中TD-DFT默认使用TDA近似，此处将tda设为false，表示使用标准的TD-DFT。如果需要获得更多的输出，还可以加上printlevel 3或更高的选项。
以下我们来分析输出文件中的结果。首先进行基态的DFT计算，随后便是TD-DFT的计算，分别输出单重态和三重态的激发能。在此之后进入旋轨耦合计算，旋轨耦合矩阵元的输出形式如下：
--------------------------------------------------------------------------------
CALCULATED SOCME BETWEENTRIPLETS AND SINGLETS
--------------------------------------------------------------------------------
Root <T|HSO|S> (Re, Im) cm-1
T S Z X Y
--------------------------------------------------------------------------------
1 0 ( 0.00 , -0.07 ) ( 0.00 , -0.16 ) ( -0.00 , -0.22 )
1 1 ( 0.00 , 0.27 ) ( 0.00 , 0.02 ) ( -0.00 , 0.05 )
1 2 ( 0.00 , 0.14 ) ( 0.00 , -0.00 ) ( -0.00 , 0.20 )
1 3 ( 0.00 , 0.21 ) ( 0.00 , 0.01 ) ( -0.00 , -0.08 )
1 4 ( 0.00 , -0.22 ) ( 0.00 , 0.14 ) ( -0.00 , 0.64 )
1 5 ( 0.00 , -0.37 ) ( 0.00 , 0.30 ) ( -0.00 , 0.57 )
1 6 ( 0.00 , -1.13 ) ( 0.00 , 1.25 ) ( -0.00 , 2.21 )
2 0 ( 0.00 , 0.21 ) ( 0.00 , 0.96 ) ( -0.00 , 0.14 )
2 1 ( 0.00 , -0.37 ) ( 0.00 , -0.02 ) ( -0.00 , 0.67 )
2 2 ( 0.00 , -0.30 ) ( 0.00 , -0.10 ) ( -0.00 , -0.03 )
2 3 ( 0.00 , -0.93 ) ( 0.00 , -0.23 ) ( -0.00 , 1.94 )旋轨耦合矩阵元通常表示为<Sm|HSO|Tn>的形式，本例中我们计算6个单重态和6个三重态，加上基态为单重态，所以共42个矩阵元。左边两列依次列出了态的序号。SOC矩阵元是一个复数，程序中会分别输出其实部和虚部，此外，SOC矩阵元还可写成X、Y、Z三个分量，程序也分别输出，而最终我们一般报道的是SOC矩阵元的模，具体算法是将6个数的平方求和再开平方。SOC矩阵元一般用cm−1为单位。
三重态事实上包含了三个子态，分别对应着磁量子数+1、0、−1。如果使用printlevel 3，程序还会输出以下内容：
--------------------------------------------------------------------------------
CALCULATED SOCME BETWEENTRIPLETS AND SINGLETS
--------------------------------------------------------------------------------
Root <T|HSO|S> (Re, Im) cm-1
T S MS= 0 -1 +1
--------------------------------------------------------------------------------
1 0 ( 0.00 , -0.07) ( -0.15 , -0.11) ( -0.15 , 0.11)
1 1 ( 0.00 , 0.27) ( 0.03 , 0.01) ( 0.03 , -0.01)
1 2 ( 0.00 , 0.14) ( 0.14 , -0.00) ( 0.14 , 0.00)
1 3 ( 0.00, 0.21) ( -0.06 , 0.01) ( -0.06 , -0.01)
1 4 ( 0.00 , -0.22) ( 0.45 , 0.10) ( 0.45 , -0.10)
1 5 ( 0.00 , -0.37) ( 0.40 , 0.21) ( 0.40 , -0.21)
1 6 ( 0.00 , -1.13) ( 1.56 , 0.88) ( 1.56 , -0.88)
2 0 ( 0.00 , 0.21) ( 0.10 , 0.68) ( 0.10 , -0.68)
2 1 ( 0.00 , -0.37) ( 0.47 , -0.01) ( 0.47 , 0.01)
2 2 ( 0.00 , -0.30) ( -0.02 , -0.07) ( -0.02 , 0.07)
2 3 ( 0.00 , -0.93) ( 1.37 , -0.16) ( 1.37 , 0.16)同样将6个数的平方和开方，得到的结果与上面的是相同的。
此时我们分析一下该体系的计算结果。原文中作者使用Dalton程序计算SOC矩阵元，计算水平为B3LYP/6-31G(d)，所得结果如下：
文中得到的S0→T1的SOC矩阵元为8.39 cm−1，这与当前的计算结果有些出入，我们得到的结果为0.27 cm−1。两个程序结果的不同是因为使用的算法、基组、积分格点等有所不同，但是我们也注意到S0与其他几个三重态之间有较大的旋轨耦合。
在SOC矩阵元之后，会输出如下内容：
- SOC stabilization of the ground state:，即SOC对基态能量的稳定化能。
- Eigenvalues of the SOC matrix:，为SOC矩阵的本征值，即混合之后的SOC态的能量。
- TD-DFT-EXCITATION SPECTRA：原始的TD-DFT激发光谱数据
- SOC CORRECTED TD-DFT/TDA-EXCITATION SPECTRA：考虑SOC后的光谱数据。
三、计算实例2
第二个体系为Ir(ppy)3，结构如下图所示
含有重金属元素Ir，在Phys. Chem. Chem. Phys., 2014,16, 14523–14530一文中研究过其磷光性质。
(1) 优化结构：
! B3LYP def2-SVP opt CPCM(CH2Cl2) D4 nopop
*xyz 0 1
C 3.393895879 -0.289566719 2.776304932
C 4.076284431 -1.374299959 2.244246939
C 3.562788518 -2.026286944 1.144617313
C 2.365194919 -1.600441421 0.563698783
C 1.656920791 -0.490717084 1.087023288
C 2.209908242 0.138316955 2.204838066
C 1.779556093 -2.262575837 -0.593808565
C 2.310998596 -3.397126504 -1.214621244
C 1.678633212 -3.946679457 -2.310344980
C 0.513725976 -3.353425417 -2.779039051
C 0.041132843 -2.238327622 -2.115943483
N 0.643735830 -1.703611386 -1.062120264
C -3.228320892 -2.843015891 2.244246939
C -1.947720074 -2.794416690 2.776304932
C -0.985168124 -1.982995155 2.204838066
C -1.253433856 -1.189576955 1.087023288
C -2.568620388 -1.248098174 0.563698783
C -3.536210228 -2.072321893 1.144617313
C -2.849226199 -0.409852865 -0.593808565
N -1.797238653 0.294314111 -1.062120264
C -1.959015004 1.083541724 -2.115943483
C -3.161014589 1.231812963 -2.779039051
C -4.257241276 0.519600723 -2.310344980
C -4.097497151 -0.302820240 -1.214621244
C -0.403486935 1.680294039 1.087023288
C -1.224740118 1.844678200 2.204838066
C -1.446175805 3.083983409 2.776304932
C -0.847963539 4.217315849 2.244246939
C -0.026578290 4.098608837 1.144617313
C 0.203425469 2.848539595 0.563698783
C 1.069670107 2.672428702 -0.593808565
C 1.786498555 3.699946745 -1.214621244
C 2.578608064 3.427078734 -2.310344980
C 2.647288613 2.121612454 -2.779039051
C 1.917882161 1.154785898 -2.115943483
N 1.153502823 1.409297275 -1.062120264
Ir 0.000000000 0.000000000 0.043448572
H 3.792342787 0.222162593 3.641507453
H 5.002872009 -1.707739869 2.689139743
H 4.102575829 -2.869480849 0.738507580
H 1.696641422 0.983426029 2.641159087
H 3.213678708 -3.847263406 -0.836341679
H 2.081671478 -4.823466100 -2.794933558
H -0.013921102 -3.747584685 -3.633027960
H -0.865398626 -1.751267403 -2.451342694
H -3.980382114 -3.478744318 2.689139743
H -1.703772945 -3.395346490 3.641507453
H 0.003351213 -1.961047587 2.641159087
H -4.536331226 -2.118194464 0.738507580
H -1.083942747 1.625090897 -2.451342694
H -3.238542989 1.885848371 -3.633027960
H -5.218079916 0.608952668 -2.794933558
H -4.938667199 -0.859495698 -0.836341679
H -1.699992635 0.977621558 2.641159087
H -2.088569843 3.173183897 3.641507453
H -1.022489895 5.186484186 2.689139743
H 0.433755396 4.987675313 0.738507580
H 1.724988490 4.706759104 -0.836341679
H 3.136408438 4.214513433 -2.794933558
H 3.252464091 1.861736314 -3.633027960
H 1.949341373 0.126176507 -2.451342694
*(2) 做SOC-TD-DFT计算：
! B3LYP ZORA ZORA-def2-TZVP SARC/J CPCM(CH2Cl2) RI-SOMF(1X) nopop
%tddft
nroots 25
dosoc true
tda false
end
%basis newgto Ir "SARC-ZORA-TZVP" end
end
*xyzfile 0 1 fac-Irppy-opt.xyzZORA (Zeroth-Order Regular Approximation to the Dirac Equation)是一种近似求解Dirac方程的方法，所用基组为ZORA-def2-TZVP，这是一种专门为ZORA方法适配的def2基组，支持的元素为H-Kr。对Ir元素使用SARC-ZORA-TZVP基组，SARC是segmented all-electron relativistically contracted的缩写，除了此处用到的SARC-ZORA外，还有SARC-DKH系列，它们分别是为ZORA和DKH2哈密顿设计的全电子基组，用于Kr之后的元素。SARC/J是一种通用的辅助基，对DKH和ZORA系列的基组都可使用。SOMF为自旋轨道平均场，是一种处理旋轨耦合积分的方法，关于1X的含义可参见ORCA手册的9.42.2 The Spin-Orbit Coupling Operator一节，RI则表示使用RI近似计算积分。对RI不了解的读者可阅读往期文章《在ORCA的SCF计算中使用RI近似》。
SOC矩阵元如下：
--------------------------------------------------------------------------------
CALCULATEDSOCME BETWEEN TRIPLETS AND SINGLETS
--------------------------------------------------------------------------------
Root <T|HSO|S> (Re, Im) cm-1
T S Z X Y
--------------------------------------------------------------------------------
1 0 ( 0.00 , 187.66 ) ( 0.00 , 9.76 ) ( -0.00 , 9.17 )
1 1 ( 0.00 , -0.22 ) ( 0.00 , 10.42 ) ( -0.00 , 10.07 )
1 2 ( 0.00 , 1.58 ) ( 0.00 , -116.41 ) ( -0.00 , 39.28 )
1 3 ( 0.00 , -0.53 ) ( 0.00 , -37.41 ) ( -0.00 , -122.18 )
1 4 ( 0.00 , 12.78 ) ( 0.00 , 357.82 ) ( -0.00 , -840.24 )
1 5 ( 0.00 , 9.20 ) ( 0.00 , 835.08 ) ( -0.00 , 356.54 )
1 6 ( 0.00 , -114.46 ) ( 0.00 , 2.75 ) ( -0.00 , 11.46 )
1 7 ( 0.00 , 0.69 ) ( 0.00 , 108.68 ) ( -0.00 , 77.58 )
1 8 ( 0.00 , 2.60 ) ( 0.00 , 80.56 ) ( -0.00 , -77.23 )
1 9 ( 0.00 , 135.13 ) ( 0.00 , 11.09 ) ( -0.00 , 1.46 )由于重元素的存在，其数值就大了很多。
在不考虑SOC时，三重态的三个子态是简并的，而考虑SOC后，能量则会发生改变，即所谓的零场分裂（zero-field splitting, ZFS）。混合后的态的能量如下：
Eigenvalues of the SOC matrix:
State: cm-1 eV
0: 0.00 0.0000
1: 20562.88 2.5495
2: 20623.57 2.5570
3: 20626.93 2.5574
4: 20758.21 2.5737
5: 20768.93 2.5750
6: 20778.83 2.5762因此其零场分裂能为State 1和State 3的能量差，即64 cm−1，与实验值85-150 cm−1相符。
四、小结
本文简单介绍了用ORCA计算SOC矩阵元的方法，希望对大家的研究有帮助。目前不少论文中对ISC过程还仅仅是拿能量接近来说事，显然是不太够的，用SOC更有说服力。最后说一点关于SOC的计算使用的结构的问题。笔者认为可以使用所研究的过程的起始状态的结构，如研究S1到三重态ISC过程则可以用S1的结构。更严格一些，也可以先尝试寻找S1与想研究的三重态的MECP点，作为计算SOC的几何结构。