计算化学公社
标题: 使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱 [打印本页]
作者Author: sobereva 时间: 2019-2-10 10:17
标题: 使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱
使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱
Using ORCA to calculate spin-orbit coupling matrix elements and plot spin-orbit coupling corrected UV-Vis spectra under TDDFT
文/Sobereva@北京科音 First release: 2019-Feb-10 Last update: 2021-Jul-10
1 前言
看本文前，强烈建议先看《使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元》（http://sobereva.com/411）了解一些旋轨耦合矩阵元计算的基本知识，以及看《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）了解TDDFT计算的相关知识。TDDFT和旋轨耦合(Spin-orbit coupling, SOC)方面的理论性的内容在本文就不多提了，笔者假设读者已经看了上面两篇具备了相关常识。如果你不会装ORCA，看《量子化学程序ORCA的安装方法》（http://sobereva.com/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://sobereva.com/attach/462/file.rar。本文用的是ORCA 4.1.1。
2 旋轨耦合矩阵元计算实例1：甲醛
首先用一个非常简单的体系甲醛来说明怎么在ORCA里做考虑了旋轨耦合的TDDFT计算、怎么去理解输出。
2021-Jul-10注：从ORCA 5.0开始，不要再像下文例子中那样自己写grid和gridx关键词，用默认的即可。详见《谈谈ORCA 5.0的新特性和改变》（http://sobereva.com/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://sobereva.com/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://sobereva.com/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://sobereva.com/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://sobereva.com/214）一文有介绍。
3 旋轨耦合矩阵元计算实例2：溶剂中的Os配合物
在前述的《使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元》一文中用Gaussian+PySOC计算过下面这个Os配合物体系，基态是单重态
(, 下载次数 Times of downloads: 578)
由于有重金属，因此我们一般都会想到用赝势基组，比同等质量的全电子基组又省时间又能顺带体现标量相对论效应。如果用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://sobereva.com/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光谱。
(, 下载次数 Times of downloads: 599)
这里使用下面的关键词算对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://sobereva.com/373）里提过的，配体的地位没有处于体系中央的过渡金属高，再加上普通泛函对高角动量函数要求比较低，因此用-f的版本不至于令误差比不带-f的版本增加太多。顺带一提，ORCA的SCF收敛做得不如Gaussian好，当前这样的关键词算某些过渡金属配合物容易遇到SCF不收敛，碰见这种情况可以尝试加上slowconv或veryslowconv关键词，此时会用对一般情况来说更慢但是更稳妥、收敛成功几率更高的收敛设定和算法，实测挺奏效。当前例子我们要算的并不是S与T不同子态间的旋轨耦合矩阵元，因此printlevel 3其实可以去掉，此时由于可以少输出许多矩阵，输出文件会小很多。当前这个任务用2*E5-2696 v3的机子在36核并行时4h13m算完。
Multiwfn是绘制各种类型电子和振动光谱的又方便又强大的程序，可在http://sobereva.com/multiwfn下载。相关信息看《Multiwfn入门tips》（http://sobereva.com/167）。在《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（http://sobereva.com/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，看到的图像如下
(, 下载次数 Times of downloads: 607)
我们也可以将考虑SOC和不考虑时候的图绘制到一起，便于考察SOC对光谱的影响。这用Multiwfn非常容易实现，怎么用Multiwfn绘制多个情况下的光谱在《使用Multiwfn绘制构象权重平均的光谱》（http://sobereva.com/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的情况）。绘制的结果如下：
(, 下载次数 Times of downloads: 627)
由图可见，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以上的曲线才能保证真实、完整。
(, 下载次数 Times of downloads: 600)
ORCA考虑SOC对光谱的影响靠的是让不考虑SOC时得到的单重态和三重态混合来实现的，只有当算的态数比较多时，也相当于令“基”完备性较高时，考虑SOC后的电子态才能被较好描述，这样得到的考虑SOC后的激发态的振子强度和激发能才可能比较合理。在绘制SOC下的UV-Vis光谱时，算的态数应当比起不考虑SOC时更多。下图是计算50、100、150个态的时候在考虑SOC后的对比图
(, 下载次数 Times of downloads: 622)
我们就关注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://sobereva.com/314）中说过自发辐射速率怎么算。对于T1的一个子态，只要把其激发能、振子强度都代进去，就可以得到从这个子态发磷光的速率。要算总的磷光速率，图省事可以直接把三个子态的发射速率取平均，但原理上更严格的方法是将三个子态的发射速率取权重平均，权重可以通过Boltzmann分布获得，见《根据Boltzmann分布计算分子不同构象所占比例》（http://sobereva.com/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分布比率，在表格最下面给出算术平均和权重平均后的发射速率。对当前情况，结果如下所示
(, 下载次数 Times of downloads: 604)
在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实现磷光速率的计算，在北京科音高级量子化学培训班里会更深入地讲。
作者Author: 小苹果 时间: 2019-2-20 08:40
谢谢，正在尝试用orca计算磷光性质，期待大神关于Dalton磷光计算好文!
作者Author: pyscf 时间: 2019-3-24 07:19
测试了除甲醛外的好几个体系 不知道为什么td-dft总是报错 但是切换到tda就没事
不知道其他人有没有遇到类似情况？
作者Author: abin 时间: 2019-3-24 08:40
http://bbs.keinsci.com/thread-2455-1-1.html
作者Author: sobereva 时间: 2019-3-24 13:58
文中的配合物也是TDDFT算的，如果文中的例子都不能运行的话，那应该就是ORCA和当前运行环境兼容的事了
作者Author: pyscf 时间: 2019-3-26 00:44
比如这个例子 用tddft报错 Negative eigenvalue in MV matrix found以及Rebuilding of the RPA subspace failed because there are no old eigenvectors但是tda就没问题
(, 下载次数 Times of downloads: 13)
作者Author: Novice 时间: 2019-4-11 08:43
Sob大大，我看ORCA发表的磷光计算的原文（https://pubs.acs.org.ccindex.cn/doi/10.1021/acs.jctc.8b00841）计算磷光速率基本上数量级很对啊，是不是您用的计算磷光速率的表格和原文的方法不同所致？ 我觉得原文“2.2. Path Integral Approach for the Rates.”和“4.2. Effect of Different Schemes for SOC Integrals.”节应该包含这部分内容，可是我是小白，搞不明白
如果能改进希望您改进一下，这样估计用ORCA也能算出合理的磷光速率了。
作者Author: sobereva 时间: 2019-4-11 15:35
帖子不需要改进
之前这篇文章在论坛里其他人也提过。光是算一些轻元素，结果尚且有的并不理想，结论没有什么普适性。而且文章只用微扰方式考虑T-S跃迁偶极矩，跟Dalton的响应函数处理比起来差得远。而且文章的计算过程不仅麻烦耗时也高，不适合一般用户效仿
作者Author: 13104468322 时间: 2019-9-12 14:45
求助：sob大大，仔细看了您关于Pysoc和ORCA计算旋轨耦合常数的帖子，发现的输入文件的分子坐标均采用了基态结构，而利用BDF计算SOC通常输入的都是激发态结构坐标。能解释下同样算SOC为什么用不同结构坐标吗
作者Author: sobereva 时间: 2019-9-13 07:02
用什么结构完全取决于你算SOC的目的是什么
作者Author: kittycq520 时间: 2019-12-12 15:29
SOB大大，关注咱们思想家公社很久了，在我之前的文章中也用过了几次咱们Multiwfn，非常好用。我最近一项工作是讨论TADF材料的SOC和反向系间窜越速率krisc之间的关系。根据您这篇博文中的工作，我进行了计算，计算结果和我的实验数据吻合的非常好，现在在写文章。不过有两个小问题想请问下SOB大大：1.按照本文中的方法，我在讨论SOC和T1→S1的krisc速率的关系时，用的SOC值是基于基态结构进行计算的，这样可不可以呢（命令行详见后文）？2.我下载了全部TADF材料和SOC相关的文章，但是并没有用该方法计算SOC的，有没有相关的文献可以让我写文章时候引用一下嘞？非常感谢大大，期待SOB大大的回复。
计算SOC用的命令行：
! B3LYP/G 6-31G** miniprint grid4 gridx4 tightSCF
%tddft
nroots 5
dosoc true
tda false
printlevel 3
end
* xyz 0 1
作者Author: sobereva 时间: 2019-12-15 10:33
1 讨论这个问题，用T1和S1之间的MECP更有意义
2 方法和PySOC那篇差不多，引那个也可以。实际上只引用ORCA手册相应部分就够了。
作者Author: kittycq520 时间: 2019-12-15 13:43
懂啦！谢谢SOB大大！
作者Author: linqiaosong 时间: 2020-2-3 22:42
不知道可不可以在算磷光发射能的时候和算磷光寿命的时候用不同的理论水平？
比如用PBE0/def2-TZVP算磷光发射能
然后用B3LYP/6-31G(d)算磷光寿命
主要是考虑到计算二次响应的时间成本比算发射能高不少，如果都用def2-TZVP有点吃不消......
作者Author: sobereva 时间: 2020-2-4 10:04
原理上可以，但小心被审稿人吐槽。所以最好先找个体系做个测试，到时候有审稿人质疑这种做法时就搬出来
作者Author: linqiaosong 时间: 2020-2-4 12:27
是个不错的主意，感谢社长
作者Author: linqiaosong 时间: 2020-2-5 22:16
- Phosphorescence transition rate from excited state no. 1
- (Triplet->singlet transition, high-temperature limit)
- Transition energy: 1.073 eV
- or 1155.362 nm
- Length gauge:
- Partial rates (H_SO): X-polarization 1.35221E-02 Transition moment : 1.015E-04
- Length gauge:
- Partial rates (H_SO): Y-polarization 2.79046E-03 Transition moment : 4.609E-05
- Length gauge:
- Partial rates (H_SO): Z-polarization 2.93512E-04 Transition moment : 1.495E-05
- Phosphorescence - length gauge:
- Oscillator strength (/2PI) (H_SO) 5.289099E-11
- Dipole strength [a.u.] (H_SO) 1.264024E-08
- Dipole strength E-40 [esu**2 cm**2] 8.166196E-04
- Total transition rate (H_SO) 5.535363E-03 s-1
- Total phosphorescence lifetime (H_SO) 1.806566E+02 s
通过振子强度和两个态之间的能量差就可以根据爱因斯坦公式计算激发态寿命，对于荧光和磷光发射都适用，往往和实验对应得不错。如果和实验值差异很大，有可能发生了显著的内转换、外转换等无辐过程。计算公式为：寿命τ=1.5/(f*ΔE^2)。这里寿命的单位为秒，ΔE以cm-1为单位。寿命的倒数就是自发辐射速率。
——"Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法"(http://sobereva.com/462)
还有一个小问题，就是Dalton计算磷光速率的时候，它的输出结果磷光寿命和按照社长帖子里的公式算出来的磷光寿命刚好差个系数3，而在公社论坛http://bbs.keinsci.com/thread-2455-1-1.html这篇帖子里，beefly说这个系数3是程序的bug，不知道社长对此怎么看？
(, 下载次数 Times of downloads: 194)
作者Author: sobereva 时间: 2020-2-6 02:59
我不觉得有除的必要
你可以看看dalton作者写的磷光计算综述，把输出数据和作者的数据进行对照
Chem. Rev. 2017, 117, 6500−6537
作者Author: qazwer 时间: 2020-2-25 18:51
您好，我不太理解文章开头说的闭壳层体系下的单重态和三重态的SOCME计算，涉及到了三重态不应该是开壳层吗？另外，像http://bbs.keinsci.com/thread-16387-1-1.html中这种情况有什么处理办法吗？
作者Author: sobereva 时间: 2020-2-25 19:02
闭壳层说的是参考态
作者Author: yinhang 时间: 2020-12-22 16:25
请问老师（b）各个单重态与各个三重态的三个子态间的总旋轨耦合矩阵元与（c）各个单重态与各个三重态间的LxSx、LySy、LzSz算符（即简化的旋轨耦合算符）对应的旋轨耦合矩阵元有什么区别？为什么我的结果中有些不同。此外，（c）中的I和J是否分别代表T态和S态
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.08) ( -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) ( -16.03 , -14.78) ( -16.03 , 14.78)
1 3 ( 0.00 , -0.05) ( -0.00 , 0.00) ( -0.00 , -0.00)
2 0 ( 0.00 , -0.19) ( 0.00 , -0.00) ( 0.00 , 0.00)
2 1 ( 0.00 , 0.02) ( -0.00 , -0.00) ( -0.00 , 0.00)
2 2 ( 0.00 , -0.00) ( 9.60 , 9.87) ( 9.60 , -9.87)
2 3 ( 0.00 , 0.02) ( 0.00 , 0.00) ( 0.00 , -0.00)
3 0 ( 0.00 , -0.00) ( -7.97 , -3.85) ( -7.97 , 3.85)
3 1 ( 0.00 , 0.00) ( 14.54 , 14.16) ( 14.54 , -14.16)
3 2 ( 0.00 , -0.16) ( 0.00 , 0.00) ( 0.00 , -0.00)
3 3 ( 0.00 , -0.00) ( 0.01 , -0.04) ( 0.01 , 0.04)
--------------------------------------------------------------
CALCULATED REDUCED SOCME BETWEEN TRIPLETS AND SINGLETS
--------------------------------------------------------------
Root <I|Lx Sx|J> <I|Ly Sy|J> <I|Lz Sz|J>
I J cm-1 cm-1 cm-1
--------------------------------------------------------------
0 0 -0.00 0.00 -0.08
0 1 -0.00 0.00 0.00
0 2 -20.91 22.67 -0.00
0 3 0.00 0.00 -0.05
1 0 -0.00 -0.00 -0.19
1 1 -0.00 0.00 0.02
1 2 13.96 -13.58 -0.00
1 3 0.00 -0.00 0.02
2 0 -5.45 11.28 -0.00
2 1 20.02 -20.56 0.00
2 2 0.00 -0.00 -0.16
2 3 -0.06 -0.01 -0.00
作者Author: sobereva 时间: 2020-12-23 07:42
不用管Lx Sx，这是简化的SOC算符部分，和实际研究涉及的量没直接关系
是
作者Author: shiyz 时间: 2021-4-29 19:57
sob大大，您给的这个例子中，计算的激发态是垂直激发的还是绝热激发？如何区分？谢谢
! B3LYP/G TZVP miniprint tightSCF grid4 pal4
%tddft
nroots 5
作者Author: wzkchem5 时间: 2021-4-29 21:20
在基态结构上计算激发态的是垂直激发，在激发态结构上计算激发态能量并且减去基态结构上的基态能量是绝热激发
作者Author: rendong 时间: 2021-6-27 23:26
尊敬的卢老师您好，我想请教三个小问题，
1，在Gaussian中激发态的几何结构得到S1态的结构，然后在ORCA中计算旋轨耦合矩阵元，是不是只有S1和T1、T2、T3... 之间的旋轨耦合矩阵元值是有意义的?
2，如果想计算S2和T1、T2、T3... 之间的旋轨耦合矩阵元，是不是必需要先优化得到S2态的几何结构再用ORCA计算?
3，如果计算磷光发射，想得到T1与S0之间的旋轨耦合矩阵元，是不是必须用优化好的T1态几何结构，再用ORCA进行计算?
作者Author: sobereva 时间: 2021-6-28 07:02
1 这些是有意义的。但“只有”不知道你指什么
2 任何结构都可以计算各个态之间的SOC矩阵元，只不过根据实际问题，某些结构下的SOC矩阵元才有实际意义
3 光靠T1和S0的SOC矩阵元是不足以得到磷光发射强度的。是
作者Author: rendong 时间: 2021-6-28 22:42
非常感谢卢老师答疑，根据您的解答，我重新理解了一下，您看这样对不对：
1，想获得S1与T1、T2、T3... 之间的SOC矩阵元就必须用S1态几何结构，想获得S2与T1、T2、T3... 之间的SOC矩阵元就必须用S2态几何结构。
2，然后我想请教老师，如果用S1态几何结构在ORCA中计算SOC矩阵元，得到的列表中S2与T1、T2、T3... 之间的SOC矩阵元的意义是什么
3，另外通过学习您的博文了解到用S0几何结构作为初猜，用Gaussian优化S1几何结构，最后一步优化的第一激发能是该分子的荧光发射能，想请教老师，最后一步优化的第二激发能的意义是什么
谢谢老师
作者Author: sobereva 时间: 2021-6-29 00:14
1、2 我已经说了，任何几何结构下都可以算Sn与Tn之间的SOC，用什么结构取决于具体目的
3 没有什么实际用处
作者Author: wzkchem5 时间: 2021-6-29 01:30
当你需要计算的是S1到Tn的ISC的时候，应该用S1结构，如果是Tn到S1的ISC，应该用Tn的结构。（这里假设你研究的都是放热的ISC过程。至于吸热的ISC过程应该用初态还是末态的结构比较好，我也不是太清楚。但是不管用哪个结构，都是近似的，严格一些的处理，比如考虑Herzberg-Teller效应的话，会需要不止一个结构下的SOC矩阵元。）所以不能笼统地说算SOC矩阵元一定要用哪个结构。
S1结构下S2和Tn之间的SOC矩阵元一般用不到，但是一些高阶的理论可能还是会用到它的，比如有的考虑Herzberg-Teller效应的方法，算S1到Tn的ISC，会涉及S1-S2的NAC矩阵元，和S2-Tn的SOC矩阵元之积。
S1结构下的第二激发能本身没有什么意义，但是此时第二激发能和第一激发能之差，对应于实验的瞬态吸收。比如你做超快光谱，先用一个激光脉冲把分子激发到第一激发态，过几纳秒测瞬态吸收光谱，可以看到第一激发态的吸收。如果你看到的峰恰好对应第一激发态到第二激发态的激发，那么这个激发能可以用S1结构下的第二激发能减去第一激发能得到。
作者Author: rendong 时间: 2021-6-29 10:00
谢谢王老师，非常感谢您详细的解释，让我茅塞顿开
作者Author: rendong 时间: 2021-6-29 10:03
谢谢卢老师，我没表达清楚，我想获得S1激发态和S2激发态分子磷光发射过程中隙间穿越有关的信息
作者Author: 幸运兔 时间: 2021-7-1 19:40
请问7.png后面的百分比是怎么算出来的，如果是Weight，T1（1）是不是：0.45282+0.19397+0.19397=0.84076，T1（2）=0.40524+0.40524=0.81 T1（3）0.39122+0.23609+0.23609=0.8634，然后我又通过帖子中给出的j计算Bolzman的excel，三个分别是20.66%，20.56%，20.46%。跟7.png对不起来，应该是我哪个地方理解错了，请问%unnormalized是怎么算的？
作者Author: sobereva 时间: 2021-7-10 10:44
今日更新的Multiwfn已兼容了ORCA 5.0版的SOC-TDDFT的输出文件绘制光谱图。
作者Author: Novice 时间: 2021-10-23 16:14
本帖最后由 Novice 于 2021-10-23 16:17 编辑
社长，我最近用ORCA5.0.1算Ir配合物的旋轨耦合时，发现考虑旋轨耦合，光谱比不考虑悬轨耦合是红移的。然后我又参考了ORCA 5最新的tutorials介绍（https://www.orcasoftware.de/tutorials_orca/spec/SOC.html），发现5.0版的
“SOC stabilization of the ground state: XXX cm-1
Eigenvalues of the SOC matrix:”后面给出的各态能量，和最后面输出的“SPIN ORBIT CORRECTED ABSORPTION SPECTRUM ”各态的能量是一致的，而不是您帖子里写的需要再加上SOC对基态的稳定能。
而且该tutorials还专门举了个Ir配合物算例，分析的也是考虑旋轨耦合后光谱红移，算的更准了（和我计算的趋势一致）。
所以是不是4.X版的orca有bug，您按照该结果给出的分析也需要重新考虑一下？
作者Author: sobereva 时间: 2021-10-24 05:06
建议去ORCA论坛问问开发者，也不排除是ORCA 5.0的bug。如果开发者能确认是4.x有问题，没把SOC stabilization of the ground state考虑到最后输出的SOC CORRECTED ABSORPTION SPECTRUM给出的激发能里，我再改帖子。
作者Author: wzkchem5 时间: 2021-10-25 16:56
我测了一下，发现从4.2.1开始就改过来了。我觉得这种问题无所谓是不是bug一说，因为不管加不加基态SOC稳定化能，只要软件行为和说明书一致，就不是bug，顶多只是不方便而已。
作者Author: rendong 时间: 2021-12-25 17:08
老师好，我想请教一下用ORCA在TDDFT下计算T1->S0的SOC矩阵元（已经优化好T1结构），inp文件中自旋多重度写1还是3是对的？
作者Author: sobereva 时间: 2021-12-25 18:35
1
作者Author: rendong 时间: 2021-12-25 18:50
好嘞，谢谢老师
作者Author: 喵星大佬 时间: 2022-3-18 22:32
很多TADF分子的S1，T1，MECP和S0的结构几乎是重合的，这种时候感觉随便用哪个都差不多
作者Author: waker 时间: 2022-7-28 17:23
不好意思老师，我还是有些疑问
（1）如果计算S1到T1的跃迁时候，分子结构采用S1的；计算T1到S1的跃迁时候，分子结构采用T1，那这两个值都写作<S1|Hsoc|T1> 吗？
（2）那用S0结构计算出来的旋轨耦合矩阵元代表什么意义？
（3）如果我想计算T1到S1的跃迁，除了更改结构以外，输入文件需要有哪些修改呢？
作者Author: sobereva 时间: 2022-7-28 19:19
是。说清楚用的什么结构
用什么结构算SOC矩阵元，公社论坛首页google框搜索，早就有过讨论
作者Author: 天天121 时间: 2022-12-4 10:38
老师如果有的分子在基态并没有稳定构型，但他在S1和T1上有稳定构型，这样的SOC能算吗ORCA？（例如A分子在基态上不能发生反应，但在激发态上却可以自发的发生反应并形成异构体B）
作者Author: sobereva 时间: 2022-12-5 02:30
任何几何结构下都能算SOC，只不过结果有没有实际价值的问题
作者Author: Hilbrac 时间: 2023-10-26 21:41
sob老师你好，我想问一下旋轨耦合矩阵元和最后的能量修正有什么关系吗？
比如我想知道旋轨耦合效应最终对基态能量的修正，是否就是把所有包含基态的旋轨耦合矩阵元的模加和起来就能得到了呢？
作者Author: wzkchem5 时间: 2023-10-27 03:17
不是，需要对角化一个QDPT Hamiltonian。参见https://www.tandfonline.com/doi/full/10.1080/00268976.2013.785611
作者Author: sobereva 时间: 2023-10-27 07:34
完全不是
常见做法是用考虑不同态之间的旋轨耦合的哈密顿重新对角化得到本征值（电子态的能量）
作者Author: Hilbrac 时间: 2023-10-27 16:44
好的，谢谢大佬。
作者Author: Hilbrac 时间: 2023-10-27 16:50
谢谢sob老师，最后还有一点点疑惑。
那就是这个例子中，如果要概括一下的话，是否就相当于在以TDDFT计算的这些态（这个例子中就是那5个单重态和5个三重态）为基矢张开的Hilbert空间中，重新计算考虑了相对论效应后的哈密顿算符的本征值及其本征态呢？如果是的话，那原则上来说，是否我用TDDFT计算越多态作为基矢，最后算出来的各个本征值及其本征态就会越接近实际呢？
作者Author: wzkchem5 时间: 2023-10-27 16:59
第一个问题，是的
第二个问题，不是，原因参见我之前发的那篇文献
作者Author: Hilbrac 时间: 2023-11-16 15:24
再次感谢大佬！最近终于有空通读了一下您发的那篇文献。
但是我看里面的例子里，确实是SOC相互作用空间里含的态数目越多，结果越趋近收敛，虽然最后收敛的值会跟实验值有所偏差。
但我感觉里面还有些没有说清楚的点，比如这个相互作用空间里到底是否应该包含参考态，里面只提到说这是一个敏感的问题，需要谨慎考虑，并没有说清楚在相互作用空间中包含参考态究竟有什么影响，这是我目前比较迷惑的一个点。
另外，按照这里面的说法，如果我要考察SOC对某个能级的影响，那我选的SOC相互作用空间内就需要尽可能包含所有空间或能量上与之相近的态（轨道），是吗？
作者Author: wzkchem5 时间: 2023-11-16 16:57
文章说了包含参考态可能导致低估参考态的能量，而且说了原因。第9页“Yet, as seen from Table 2, the
energy of 1S0 becomes increasingly lower as the interaction space is enlarged, to a much larger extent than the excited spinor states. While the latter is due to the lack of spin–orbital relaxations, the former does arise from such effects as the interaction matrix elements (line 6 in Table 1) do involve rotations between the occupied and virtual orbitals.”
对
作者Author: Hilbrac 时间: 2023-11-17 15:00
感谢大佬指正。
①这里英文的大致意思是否是，随着相互作用空间的扩大，算出来的激发态和参考态的能量都逐渐降低，且参考态的能量降得会更快，而能量降低的原因，前者是缺乏自旋轨道弛豫，后者是因为相互作用矩阵元素涉及到了占据轨道和虚轨道间的旋转？
②这里的spin–orbital relaxations指的是否就是（含自旋的）轨道弛豫呢？
③“相互作用矩阵元涉及到了占据轨道和虚轨道间的旋转”有些太过抽象，它具体指的是什么意思呢？
③为避免低估参考态能量，我需要如何先验地知道是否应该将参考态包含进相互作用空间呢？
作者Author: wzkchem5 时间: 2023-11-17 21:39
（1）基本正确，但“缺乏自旋轨道弛豫”是前者能量降低得慢的原因，不是前者能量降低的原因
（2）可以认为是SOC导致的轨道弛豫。这个轨道是否含自旋不重要，因为要做SOC一般都是用RHF/RKS或ROHF/ROKS的，本来alpha、beta自旋轨道就是一样的
（3）类比一下当参考波函数不是收敛的HF波函数的时候（也就是Brillouin条件不成立的时候），你在参考波函数以及它的所有单激发行列式组成的基上对角化CIS矩阵，会得到什么结果。原理和这个是相通的
（4）不容易知道。这是SOC校正的TDDFT方法的一个根本缺陷，没有办法不依赖经验或测试就知道是否需要包括参考态，以及如果包括参考态的话激发态算多少个。对于足够轻的元素，考虑不考虑的影响都不大（但SOC对光谱本身也没有可见的影响）；对于足够重的元素，本来就要改用two-component方法；对于不是太轻也不是太重的元素，可能需要做一些测试，或参考文献类似体系，或改用NEVPT2等方法
作者Author: Hilbrac 时间: 2023-11-18 12:11
重新看了一下，“Yet, as seen from Table 2, the energy of 1S0 becomes increasingly lower as the interaction space is enlarged, to a much larger extent than the excited spinor states. While the latter is due to the lack of spin–orbital relaxations, the former does arise from such effects as the interaction matrix elements (line 6 in Table 1) do involve rotations between the occupied and virtual orbitals.”这句话意思大致是：这种方法会考虑参考态的自旋轨道弛豫，但并不会考虑到激发态的自旋轨道弛豫，所以后者能量随相互作用空间增大而降得更慢一些？
如果是这样的话，那是否可以认为是自旋轨道弛豫导致了参考态能量的低估呢？
作者Author: wzkchem5 时间: 2023-11-18 16:55
对
可以这么认为
作者Author: Hilbrac 时间: 2023-11-20 14:47
感谢大佬！
那我还有一点疑问。首先，对电子SOC的处理可以近似分为两阶，一阶为电子自旋与自身角动量的耦合，二阶为电子自旋与其它电子轨道角动量的耦合。（①这里应该没说错吧？）
这篇文献里也有提到说，这个sf-X2C-S-TD-DFT方法是用了有效单电子算符将SOC处理到了一阶。（见下图一）
再来看一下Sob老师的博文使用Gaussian+PySOC在TDDFT下计算旋轨耦合矩阵元 - 思想家公社的门口：量子化学·分子模拟·二次元 (sobereva.com)，里面有这么一段话。（见下图二）
对于Sob老师这段话，我的理解就是，这种用TDDFT来计算SOC矫正的方法确实只考虑到了电子SOC处理的第一阶，仅处理电子自旋与自身轨道角动量的耦合。（只不过Sob老师这里把电子自身轨道角动量划分为了绕核角动量与绕其它电子的角动量）
然而，这篇文献后面又说这个sf-X2C-S-TD-DFT方法里面会包含有spin-other-orbit interactions（②某电子自旋角动量与其它电子轨道角动量的相互作用？），也就是说这个方法是涵盖了些许SOC二阶项的？这不是与前面提到的“用有效单电子算符将SOC处理到了一阶”矛盾了吗？（见下图三）
还是说，这篇文献里提到的Hso算符与Sob老师提到的Hso算符划分方法不一样，导致前者还是涵盖了一些SOC二阶项的？
作者Author: wzkchem5 时间: 2023-11-20 17:10
不是的，SOC的贡献有很多种分解方法，分解成“电子自旋与自身角动量的耦合”和“电子自旋与其它电子轨道角动量的耦合”是一种做法，分解成单电子项和多电子项是一种做法，分解成SOC算符的Taylor级数展开又是一种做法。sob老师的博文是按第二种分解的，sf-X2C-S-TD-DFT-SOC的文章是按第三种分解的，第一种分解在理论讨论中用得比较多，但是实际计算用得不多，因为两项都计算不比只计算第一项难。
只有第三种分解方法才能把各项叫做“阶”（因为Taylor展开有“几阶”的概念），前两种分解方法至多只能说是第一项和第二项，不能叫一阶和二阶
作者Author: Hilbrac 时间: 2023-11-21 10:35
那这三种分解方法之间是否有对应关系呢？
①sob老师用的第二种分解方法，即单电子和多电子项，二者加和是否就对应于第一种分解方法中的“电子自旋与自身轨道角动量的耦合”呢？
②这篇文献中用到的第三种方法，Taylor展开分解方法，与前两种分解方法有何关联吗？比如它的一阶和二阶分别代表什么呢？
作者Author: wzkchem5 时间: 2023-11-21 18:02
（1）单电子项+多电子项=“电子自旋与自身角动量的耦合”+“电子自旋与其它电子轨道角动量的耦合”。进一步地，多电子项可以分成库伦项和交换项，单电子项和库伦项里只包含“电子自旋与自身角动量的耦合”的贡献，交换项既包含“电子自旋与自身角动量的耦合”的贡献又包含“电子自旋与其它电子轨道角动量的耦合”的贡献。
（2）我之前说的不准确，不是Taylor展开，是DKH展开，但是如果对DKH不熟悉，可以类比成Taylor展开来理解。文章后面提到了so-DKHn，看那里的引用文献就知道了。每一阶都既包含了一部分单电子项，又包含了一部分多电子项；而且既包含了一部分“电子自旋与自身角动量的耦合”，又包含了一部分“电子自旋与其它电子轨道角动量的耦合”。
作者Author: Hilbrac 时间: 2023-11-21 20:34
①那sob老师这里的双电子旋轨耦合算符H_SO2和这篇文献里提到的two-electron integrals虽然不同，但是否都属于双电子项的一部分，也即多电子项的一部分？（双电子项∈多电子项？）
②我之所以认为上述二者不同，是因为在我看来，前者指仅包含了电子自旋与自身轨道角动量的耦合，后者却涉及到了电子自旋与其它电子轨道角动量的耦合，所以二者理应不同。这是否是因为前者仅涉及库仑项，而后者库仑项和交换项都涉及到了？
③文献里提到的Gaunt interaction又是何方神圣？它与交换项有何关系呢？
作者Author: wzkchem5 时间: 2023-11-21 21:40
（1）我又看了一下，我可能理解错了sob老师的记号体系。像我之前说的，SOC算符既可以分解成单电子算符的贡献和双电子算符的贡献，又可以分解成spin-same-orbit和spin-other-orbit的贡献。而sob老师把H_SO1叫做“单电子旋轨耦合算符”，H_SO2叫做“双电子旋轨耦合算符”，紧接着又把前者定义成spin-same-orbit，后者定义成spin-other-orbit，这种叫法我感觉比较奇怪，因为电子-电子相互作用对spin-same-orbit项也有很大贡献（也就是电子之间的库伦排斥加上一部分的交换项贡献），也就是电子A对【电子B的自旋和电子B的轨道角动量的耦合】的贡献，这个贡献应该叫单电子还是叫双电子？很难说。我之前是按照“单电子旋轨耦合算符”、“双电子旋轨耦合算符”的名字，武断地认为前者是单电子积分的贡献，后者是双电子积分的贡献，但这个可能不是sob老师的本意。我也不确定这个叫法是不是广泛认可的，但是spin-same-orbit、spin-other-orbit是广泛认可的叫法，也不容易有歧义。至于”多电子项“是我说的不规范，应该叫双电子项更合理，因为起码在Dirac方程的框架下，没有三电子或更多电子的相互作用的贡献，如果讨论量子电动力学那另说。
如果sob老师说的H_SO2就是spin-other-orbit项，那么H_SO2是双电子项的一部分，而且起码当采用SOMF的时候，它是交换项的一部分。
（2）不是的。起码在SOMF级别下，spin-same-orbit = 单电子项+库伦项+交换项的一部分，spin-other-orbit = 交换项的另一部分。所以spin-same-orbit也包含一部分交换项的贡献，但是包括得不全
（3）看https://pubs.aip.org/aip/jcp/art ... braic-two-component
作者Author: Hilbrac 时间: 2023-11-22 16:06
感谢大佬耐心解答！实在是很详细，讲得也很清楚！
另外，对于博文中的“双电子旋轨耦合算符H_SO2”，sob老师的定义是，“电子自旋与这个电子绕其它电子运动的耦合”，所以我想sob老师这里的H_SO2指的可能是“spin-same-orbit = 单电子项+库伦项+交换项的一部分”中的“库伦项+交换项的一部分”？所以个人还是认为sob老师只考虑了spin-same-orbit。
然后我还有两点小疑惑：
①为什么这里的角动量耦合都只考虑“自旋-轨道耦合”，而没有考虑“自旋-自旋耦合”和“轨道-轨道耦合”呢？
②大佬这里说的SOMF是“spin-orbital mean field”吗？如果是的话，我目前只知道这是一种类似于Hartree平均场近似的做法，但进一步的知识我就不得而知了，如果可以的话，希望大佬能再详细讲讲这个SOMF。
作者Author: wzkchem5 时间: 2023-11-22 16:50
本帖最后由 wzkchem5 于 2023-11-22 10:02 编辑
（1）自旋-自旋耦合（SSC）比较小，但是想算也是可以算的，例如orca算零场分裂（ZFS）就考虑了SSC，orca说明书相关章节有详细讲解。（另外还有核自旋-核自旋的SSC，也就是NMR的耦合常数。）轨道-轨道耦合在非相对论情况下已经考虑过了，而且不是作为单独的一项计算的，把这一项从通常的能量表达式里拆出来反而麻烦（而且未必有严格定义），所以一般不讨论。
（2）对。参见Neese, F., Efficient and accurate approximations to the molecular spin-orbit coupling operator and their use in molecular g-tensor calculations. The Journal of chemical physics 2005, 122 (3), 034107.
作者Author: Hilbrac 时间: 2023-11-24 14:33
感谢大佬的又一篇文献，辛苦了！
①对于这篇文献中的sf-X2C-S-TD-DFT-SOC方法所考虑到的DKH展开，这种展开的第一阶里面虽然含有一部分“电子自旋与其它电子轨道角动量的耦合”，也就是一部分spin-other-orbital，但这一部分是不是只占第一阶很小的比例呢？（因为文献后面提到了一下这一项带来的能量修正确实很小），如果是的话，这个spin-other-orbital是不是在展开的第二阶里面才会很显著呢？
②还有关于前面提到的把参考态放进相互作用空间里，扩大这个空间可能会导致参考态能量的低估这件事，我还想确认一下，这个根本缺陷是来源于SI方法吗？还是另有原因，仅sf-X2C-S-TD-DFT-SOC这个方法有这个缺陷而已？
作者Author: wzkchem5 时间: 2023-11-24 19:25
（1）起码对于足够轻的元素而言，DKH展开里面的每一阶里的每一项，应该都比前一阶对应的项小。所以起码对于足够轻的元素而言，DKH第一阶的spin-same-orbit（注意是orbit不是orbital）>DKH第一阶的spin-other-orbit>DKH第二阶的spin-same-orbit>DKH第二阶的spin-other-orbit。但正如Neese 2005文章所说，第一个大于号大概大一个数量级，10%能不能叫很小的比例就见仁见智了。
推而广之，一般量化里面，很少有一阶下很小但不是0，二阶下很大的量。不论是Taylor展开还是DKH展开，还是任何能展开到无穷阶的展开，一般都是只有在这个展开收敛的时候，取前面几阶的结果才比较靠谱。那么如果一个量在二阶比一阶大，三阶是不是很可能比二阶还大？那取二阶很可能就是不靠谱的了，也就是这个方法你经常听到的阶数一定是大于2的。但是你看DKH，一般最常见的就是DKH2，虽然DKH3、DKH4也有人用，但是明显比用DKH2的少，这就说明对于大部分问题来说二阶项远大于三阶项，那么二阶项很难有理由远大于一阶项。一阶项可以忽略的情况，往往是因为一阶项有特殊的数学性质，导致一阶项恰好为0（例如因为Brillouin条件，CIS基态等于HF基态，导致CISD的D贡献远大于S贡献），但是很少有数学性质能让只取一阶项的时候一阶项特别小却不是0。
（2）来源于在SI方法里只考虑单激发。一方面，如果做FCI-SOC，那么不会有参考态低估的问题；另一方面，如果做变分的二分量TDDFT，也不会有参考态低估的问题
作者Author: Hilbrac 时间: 2023-11-27 15:20
再次感谢解答！
另外，对于SOC校正的计算，我目前的理解是这样的：像SI方法这种，是先用TD-DFT算出考虑了电子关联的波函数，再以它们作为参考态来算SOC矩阵，以此得到考虑了SOC后的各个激发态的激发能。但这是把电子关联和SOC分开处理的，如果对于SOC效应很强的情况，比如锕系元素这种，那么SOC校正会显著改变能级差，从而影响到spin-other-orbit项进而又改变SOC项自身，也就是说，这个时候我们不能把电子子关联和SOC像SI方法这样分开处理，而必须用四分量或二分量方法，或者便宜些的DFT+SOC赝势方法。不知我这样的理解是否正确呢？
如果我上述的理解是正确的，我最近又看到了这样一篇文献：https://doi.org/10.1021/jp9044594，里面讲了有关铀的计算，尤其是这里面提到了SI方法，比如Method部分的最后一段，我的疑惑就是，如果铀原子的SOC效应很大，如我上面所说那样的话，那这里为什么还能用SI方法呢？
作者Author: wzkchem5 时间: 2024-1-5 17:35
可能因为需要的精度低。不能抛开精度要求而说某某方法在某某情况下一定不准
作者Author: 1730231594 时间: 2025-3-14 20:17
sob老师我想问一下您在计算过渡金属配合物的SOC时说DKH-def2-TZVP对前4周期是全电子基组，SARC-DKH-TZVP对Xe之后的元素是全电子基组，我想问一下有没有针对第5周期的全电子基组
作者Author: sobereva 时间: 2025-3-15 00:17
2020年提出了Rb到Xe的SARC基组
作者Author: 陈AG 时间: 2025-6-30 14:31
本帖最后由 陈AG 于 2025-6-30 15:32 编辑
老师，我要画出磷光光谱的话，我先拿高斯UDFT优化完T1的文件创建inp文件,下面这个输入需要怎么改才是磷光光谱的呢
! 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=5 TDA false dosoc true printlevel 3 end
# BP86/SDD/TZVP opted
* xyz 0 1
作者Author: sobereva 时间: 2025-6-30 18:18
在multiwfn里按下文绘制光谱的时候，用选项20 Modify oscillator strengths修改振子强度，把发射态以外的态的振子强度都设为0，再绘制的光谱就对应于磷光发射谱
使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图
http://sobereva.com/224
按照Kasha规则，发射态是T1态，因此把T1态考虑SOC后分裂出的三个子态的振子强度保留即可（严格来说应该考虑子态间的Boltzmann分布，即手动修改三个子态的振子强度乘上Boltzmann分布比例）。
作者Author: 陈AG 时间: 2025-7-30 08:52
1."T1态考虑SOC后分裂出的三个子态的振子强度保留",老师，请问这个怎么操作呀，附件是orca输出文件。
2.我选了20之后，2-400的振子强度设为0，图片是图1这样子的，怎么才可以获得像实验一样的光谱（图2）呢
作者Author: sobereva 时间: 2025-7-31 13:16
Multiwfn绘制光谱的界面里可以手动设置各个激发态的振子强度，下文明确说了
使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图
http://sobereva.com/224
先把坐标轴调一致了再说。而且我根本不知道你这是什么光谱
认真看
在网上求助计算化学问题的时候必须把问题描述得详细、具体、准确、清楚、完整
http://sobereva.com/620（http://bbs.keinsci.com/thread-25787-1-1.html）
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |