Study note：自旋轨道耦合矩阵元的计算方法
自旋轨道耦合(Spin-Orbit Coupling, SOC)常简称旋轨耦合，是一个量子力学概念，描述了粒子的自旋磁矩与其轨道角动量之间相互作用。它的数学表示为：
\[H_{SOC} = \lambda \mathbf{L} \cdot \mathbf{S}\]其中：$H_{SOC}$是自旋轨道耦合的哈密顿量，$\lambda$ 是耦合常数，$\mathbf{L}$ 是轨道角动量算符，$\mathbf{S}$ 是自旋角动量算符。当某单重态$S_n$与某三重态$T_n$之间存在显著的旋轨耦合时，自旋角动量会偏离期望值，引起两个态的混合，对应波函数$\Psi_S$与$\Psi_T$之间的耦合强度由旋轨耦合矩阵元(又称旋轨耦合常数，SOCME)给出：
\[\langle \Psi_S |H_{SOC}| \Psi_T \rangle\]SOCME是可以使用量子化学方法计算得到的。有了SOCME，再加上单重态与三重态的能差$\Delta E_{ST}$，就可以计算ISC的速率了。
总之，本文记录一下SOCME的不同计算方法，以备不时之需。
1. PySOC
PySOC计算旋轨耦合常数基于Gaussian 09，需要先做一个TD(50-50,nstate=8)计算并保留swf文件，再用PySOC进行处理。由于笔者的服务器太老，无法兼容较新的glbic库，不能运行PySOC，因此跳过该法。
PS: 道听途说PySOC结合def2基组算主族重元素的SOC会很大，但这是错误的，笔者打算有机会验证一下。
参考：使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元
2. ORCA (Recommend)
ORCA 4.1开始支持计算TD-DFT级别下闭壳层体系的旋轨耦合常数，操作比较简单。首先，对体系进行结构优化，这部分可以在Gaussian中完成。随后，将Gaussian输出文件转换为xyz格式：
1
obabel coumarin.log -O coumarin.xyz
创建一个soc.inp，编辑计算输入文件：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
! WB97X-D3 def2-SVP miniprint tightSCF
%pal
nprocs 32
end
%maxcore 3000
%tddft
nroots 5
dosoc true
tda false
printlevel 3
end
* xyzfile 0 1 coumarin.xyz
//要留空行
将该文件提交计算，在输出中找到这部分：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
CALCULATED SOCME BETWEEN TRIPLETS AND SINGLETS
--------------------------------------------------------------------------------
Root <T|HSO|S> (Re, Im) cm-1
T S MS= 0 -1 +1
--------------------------------------------------------------------------------
1 0 ( 0.00 , 0.07) ( 0.06 , -0.01) ( 0.06 , 0.01)
1 1 ( 0.00 , -0.00) ( 0.04 , 0.00) ( 0.04 , -0.00)
1 2 ( 0.00 , -0.90) ( 6.95 , 3.01) ( 6.95 , -3.01)
1 3 ( 0.00 , -0.03) ( 0.16 , 0.00) ( 0.16 , -0.00)
1 4 ( 0.00 , -0.06) ( 0.09 , -0.06) ( 0.09 , 0.06)
1 5 ( 0.00 , 1.49) ( -8.36 , 7.13) ( -8.36 , -7.13)
2 0 ( 0.00 , -0.07) ( -0.01 , -0.08) ( -0.01 , 0.08)
2 1 ( 0.00 , -0.00) ( -0.02 , 0.01) ( -0.02 , -0.01)
2 2 ( 0.00 , 0.76) ( -5.88 , -2.81) ( -5.88 , 2.81)
2 3 ( 0.00 , -0.02) ( -0.01 , 0.10) ( -0.01 , -0.10)
2 4 ( 0.00 , 0.01) ( -0.04 , 0.11) ( -0.04 , -0.11)
...
这里每一行的输出是$S_n$态与$T_n$态的三个不同自旋磁量子数的子态的旋轨耦合矩阵元的模。简单来讲，想得到SOCME值，就把对应行的六个数值求平方加和再开根号。
3. Gaussian
CASSCF级别下，Gaussian支持$Z^{eff}$的方式计算旋轨耦合，但只拟合了1-17号元素。运行计算时，需要在CAS里写spinorbit，此时将默认使用stateaverage，需要输入权重：
1
2
3
4
5
6
7
8
9
10
%nprocshared=64
%mem=192GB
%oldchk=CASSCF.chk
%chk=CAS_SO.chk
#p CAS(6,6,slaterdet,spinorbit)/def2svp geom=allchk guess=read
0.500000000.50000000
1 3
CASSCF收敛后，在Spin-Orbit部分可以找到算出的旋轨耦合：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Final one electron symbolic density matrix:
1 2 3 4
1 0.200000D+01
2 0.167427D-04 0.999999D+00
3 0.236453D-04 0.133952D-06 0.100000D+01
4 0.224958D-04 -0.187947D-02 0.267663D-03 0.102263D-05
*****************************
Spin-Orbit coupling program.
*****************************
1st state is 2
2nd state is 4
Transition Spin Density Matrix
1 2 3 4
1 0.107488D-16 -0.524237D-12 0.301863D-13 -0.111019D-11
2 0.660001D-13 0.480268D-11 -0.437174D-12 -0.165754D-09
3 0.185491D-12 -0.457050D-13 -0.471827D-11 0.513046D-09
4 0.114641D-12 0.783182D-14 0.394856D-14 -0.844224D-13
Magnitude in x-direction= 0.0 cm-1
Magnitude in y-direction= -0.0 cm-1
Magnitude in z-direction= 0.0 cm-1
Total magnitude= 0.0 cm-1
MCSCF converged.