应用案例
穆斯堡尔谱
1958 年，R. L. Mössbauer在研究γ射线共振吸收现象时，发现了 穆斯堡尔效应 。 同质异能位移（isomer shift；\(\delta^{IS}\)）是穆斯堡尔谱的重要观测参量之一，它源于具有一定尺寸的原 子核与周围电子分布之间的库仑相互作用。当原子处在不同的外界环境，导致原子核附近的库仑能发生变化， 而 \(\delta^{IS}\) 对这种变化极其敏感，因此可以用来研究原子的氧化态，自旋态，以及配位环境。 穆斯堡尔谱的另一个重要观测参量是核四级分裂（nuclear quadrupole splitting；\(\Delta E_{Q}\)）， 它来源于核自旋量子数I > 1/2的原子核产生的电四极矩（eQ；这里e是一个质子电荷，Q是核四极矩NQM；nuclear quadrupole moment） 与原子核周围电场梯度（electric field gradient；EFG）之间的相互作用。 此外，当原子置于外磁场中，由于原子核存在核磁矩，穆斯堡尔谱会进一步发生 Zeeman 分裂。
\(\delta^{IS}\) 可以表示为重元素在测试体系A和参照体系R中“有效接触密度”（effective contact density；ED） 或“接触密度”（contact density；CD）变化量的线性函数：
以上两个公式称为 校准方程 （calibration equation），其中 \(\delta^{IS}\) 是测试体系A相对于参照系R的同质异能位移实验值， ED或CD值 \(\rho_{A}\) 和 \(\rho_{R}\) 可以通过理论计算获得， \(\alpha\) 和 \(\beta\) 是待拟合的参数，其中 \(\alpha\) 也称为核标定常数（nuclear calibration constant）， \(C\) 是任意量，一般取ED或CD的整数部分。考虑到 \(\rho_{R}\) 理论值存在误差，一般用后一个公式进行拟合。
备注
CD假定原子核附近的电子密度是均匀分布的，因此可以用一个点的电子密度代替（一般取原子核中心位置的密度， 但也有程序取一系列样点的密度做加权平均）；ED考虑了电子密度的非均匀分布，原则上比前者更合理。 很多程序计算的是CD，而在BDF中计算ED，二者近似满足换算关系（参见文献 [88] 的Table S2、S3）， 换算因子可以吸收到 \(\alpha\) 中。
为了准确计算ED及其相对变化量，需要考虑以下两个因素：
原子核具有一定的尺寸分布，而默认的点电荷近似可能会导致几个数量级的误差！为此需要在
xuanyuan模块设定nuclear=1。需要考虑相对论效应。对于重元素这是显而易见的；即便对于某些轻元素也必须要考虑相对论效应， 这是因为非相对论情况下的 p 电子在原子核附近没有分布（参见文献 [88] 的Table S6）， 从而对 p 区元素的ED造成定性错误。在BDF中可以用sf-X2C哈密顿或其局域变体考虑标量相对论效应， 通过
xuanyuan模块中的heff=21（标准sf-X2C），22（sf-X2C-AXR），或23（sf-X2C-AU）指定。
铁（ \(\ce{^{57}Fe}\) ）化合物的有效接触密度
γ射线吸收和发射的几率正比于 \(\exp(-E_\gamma^2)\) ，当核激发能 \(E_\gamma\) 超过200 keV后，一般很难观测到穆斯堡尔谱， 由此只能对少数同位素探测穆斯堡尔谱。 \(\ce{^{57}Fe}\) 就属于适合实验测量的同位素之一，不过在理论计算中，一般不对这些同位素进行区分。
计算ED需要非常陡峭（也就是高斯指数非常大）的 s 型高斯基函数才能准确描述电子在铁原子核附近的分布；
对于存在 p 价电子的 p 区元素，还需要非常陡峭的 p 型高斯基函数，而基组库中的标准收缩基组通常不符合要求。
建议采用cc-pVnZ型或ANO型全电子相对论基组，并把其中的 s 函数（ p 区元素还有 p 函数）进行非收缩处理。
在以下的全电子相对论计算中，铁采用ANO-R2基组（具有三ζ精度），
并把 s 函数做非收缩处理，也就是删除 s 函数的收缩因子部分，并把收缩度从6改为0，然后存为其它文件名（如ANO-R2-ED）。
由于铁没有4 p 价电子， p 函数不需要修改。
把ANO-R2-ED复制到$BDF_WORKDIR目录下，供后面的计算调用（下载链接 ano-r2-ed.zip ）。
注意：对于非标准基组，基组名的字母必须全部大写。
我们对于铁的一系列模型体系化合物进行相对论密度泛函理论计算，泛函选取PBE0，相对论哈密顿用sf-X2C-AU。 除了铁以外的轻元素全部用def2-TZVPP基组，它在Kr元素之前属于全电子基组，虽然是非相对论的，但用于前18号元素的相对论计算是允许的。 自旋多重度和分子坐标来自文献 [89] 。以 \(\ce{FeF6^{4-}}\) 为例，输入如下：
$compass
title
FeF_6^4-
basis-block
def2-tzvpp
Fe = ANO-R2-ED
end basis
geometry # 分子直角坐标，单位：埃
Fe -0.000035 0.000012 0.000014
F 2.116808 -0.003546 0.032360
F -2.116824 0.001611 -0.030945
F -0.003602 2.164955 0.001902
F 0.001648 -2.165219 -0.003295
F 0.032586 0.003638 2.109790
F -0.030580 -0.001452 -2.109825
end geometry
MPEC+cosx # 使用MPEC+COSX加速
$end
$xuanyuan
heff # sf-X2C-AU；计算ED必须选21-23中的一个
23
nuclear # 高斯有限核模型；ED必须设为1
1
$end
$scf
charge
-4
spinmulti
5
uks
dft functional
pbe0
grid # DFT计算ED需要用精密格点
ultra fine
reled
26 # 只计算Fe的ED（对于本例，10至26的整数等价）
$end
计算完成后，在SCF布居分析信息之后可以找到ED结果：
Relativistic effective contact densities for the atoms with Za > 25
----------------------------------------------------------------
No. Iatm Za RMS (fm) Rho (a.u.)
----------------------------------------------------------------
1 1 26 3.76842 14552.65555
----------------------------------------------------------------
以此为例，完成其它铁化合物分子的ED计算（输入文件下载链接 ed-fe.zip ）。
ED结果以及 \(\delta^{IS}\) 实验值 [89] 列于下表：
分子 |
2S+1 |
\(\delta^{IS}\) (mm/s) |
ED ( \(bohr^{-3}\) ) |
|---|---|---|---|
\(\ce{FeCl4^{2-}}\)
\(\ce{Fe(CN)6^{4-}}\)
\(\ce{FeF6^{4-}}\)
\(\ce{FeCl4^-}\)
\(\ce{Fe(CN)6^{3-}}\)
\(\ce{FeF6^{3-}}\)
\(\ce{Fe(H2O)6^{3+}}\)
\(\ce{FeO4^{2-}}\)
\(\ce{Fe(CO)5}\)
|
5
1
5
6
2
6
6
3
1
|
+0.90
-0.02
+1.34
+0.19
-0.13
+0.48
+0.51
-0.87
-0.18
|
14551.76
14555.78
14552.68
14553.98
14556.08
14553.01
14554.12
14558.17
14556.37
|
用这些数据进行拟合，得到校准方程
可见拟合误差比较大，这可能是以下原因造成的：
样本太少
穆斯堡尔谱是对固态的真实体系测量的，与计算所用的气态离子模型不一致。用团簇模型、溶剂化模型 [90] 、嵌入模型 [91] 可能更合适。
铁的某些化合物存在强关联，需要测试其它泛函，或者换成适合描述强关联体系的方法
有了校准方程后，就可以对一些铁的体系预测 \(\delta^{IS}\) 。例如交错状的二环戊二烯基铁 [92] ， 通过以上密度泛函理论计算得到ED为14554.25 a.u.，代入校准方程得到 \(\delta^{IS}\) 为0.37 mm/s， 与实验值0.53 mm/s [92] 基本接近。
计算重元素化合物有效接触密度的注意事项
对于4d以上的元素，经验表明默认的高斯指数还不足以描述原子核附近的电子分布，需要额外补充一些更陡峭的高斯指数。 例如，选择cc-pVnZ型或ANO型标准基组中最陡峭的4-6个 s 型高斯指数α（ p 区重元素还要考虑 p 型高斯指数），它们近似满足以下线性关系：
通过线性拟合得到参数A、B，再通过外推（i的间隔取-0.5或-1），即可得到更陡峭的高斯指数。 一般加入2-5个更陡峭的 s 函数、1-3个更陡峭的 p 函数即可满足要求，但是要避免用10 11 以上的高斯指数， 因为这可能会造成数值不稳定。
铁（ \(\ce{^{57}Fe}\) ）化合物的EFG计算
EFG计算对相对论哈密顿的要求与有效接触密度的计算类似，但对基组的要求不同。
只有 s 电子以及少量 p 电子在原子核附近存在非零的分布，因此有效接触密度计算中只需要修改 s 、 p 基函数
原子核形变产生的电四极矩只能与轨道角动量 L > 0 电子的EFG发生相互作用，因此不必修改 s 基函数。首先把 p 函数进行非收缩处理（删除 p 函数的收缩因子部分，并把收缩度改为0），并酌情添加1-2个陡峭的 p 型高斯函数。对于存在 d 价轨道的过渡元素（镧系、锕系还有 f ），需要把 d （ f ）函数进行非收缩处理；由于 d 、 f 轨道距离原子核较远，不必添加更陡峭的 d 、 f 函数。
若同时计算有效接触密度与EFG，对基函数的修改要满足以上两条要求。
EFG计算的关键词为 relefg 。
例如，同时计算有效接触密度与EFG，以上算例的 SCF 模块输入需要改为：
$scf
charge
-4
spinmulti
5
uks
dft functional
pbe0
grid # DFT计算EFG需要用精密格点
ultra fine
relefg
26 # 只计算Fe的EFG张量
reled
26 # 只计算Fe的ED
$end
计算完成后，在SCF布居分析信息以及ED结果之后，可以找到EFG张量的结果：
Relativistic electric field gradients for the atoms with Za > 25
-----------------------------------------------------------------------------
No. Iatm Za RMS (fm) EFG tensor (a.u.)
-----------------------------------------------------------------------------
1 1 26 3.76842 -0.1061 -0.0023 0.1850
-0.0023 0.0395 -0.0018
0.1850 -0.0018 0.0666
eta Vaa Vbb Vcc
0.64736 0.0395 0.1844 -0.2239
NQCC = -8.4172 MHz with Q(ISO-057) = 160.00 mbarn
-----------------------------------------------------------------------------
在EFG 张量的9 个分量中，6 个非对角元是行列对称的；3 个对角线元之和为零。如果选择一个特殊的坐标系 \(\{\vec{a},\vec{b},\vec{c}\}\) （即EFG 张量的主轴或特征矢量）， 使得非对角元为零，而对角元（即特征值）满足 \(|V_{aa}| \le |V_{bb}| \le |V_{cc}|\) ，此时EFG 张量只需要两个非独立参数来表示就可以了， 即主值 \(V_{cc}\) 和非对称参数 \(\eta = |(V_{aa} − V_{bb})/V_{cc}|\) （\(0 \le \eta \le 1\)）。当η = 0 时，EFG 张量为轴对称。 在本例中，η = 0.64736， \(V_{cc}\) = -0.2239 a.u. 。
注意
非阿贝尔群分子的简并态在计算EFG时，单个分支的 \(V_{cc}\) 和η 一般没有意义。必须对简并态的所有分支（通过在SCF中指定占据数）分别计算EFG张量， 对它们做平均后再计算 \(V_{cc}\) 和η。
对于孤立原子， \(V_{aa} = V_{bb} = V_{cc} = 0\) ；对于线形分子（包括双原子分子）， \(V_{cc} = V_{zz}\) （分子轴为z）。 利用这一特点，BDF可以对开壳层原子、线形分子简并态的EFG结果进行校正。
核四极矩与EFG 之间的相互作用通常用核四极耦合常数（NQCC；nuclear quadrupole coupling constant） \(eQq\) 来衡量（在一些文献中也写作 \(eqQ\) ），定义为
其中 \(V_{cc}\) 仍取原子单位，核四极矩Q的单位是Barn（1 Barn = 1.0e-28 平方米）， \(eQq\) 的单位是MHz。 当同位素的Q实验值已知时，程序会打印 \(eQq\) ，在本例中是-8.4172 MHz。
穆斯堡尔谱测量的核四极分裂 \(\Delta E_{Q}\) 与NQCC满足一定的关系。例如，对于 \(\ce{^{57}Fe}\) 的I=1/2 → I=3/2核激发跃迁， γ 线激发能为14.412497 KeV（约34.85e11 MHz），有
其中，单位换算因子为 1 mm/s = 11.6248 MHz。 \(\Delta E_{Q}\) 的理论结果可以直接和穆斯堡尔谱实验值进行对比，还可以结合之前的ED结果，验证Fe的价态指认是否正确。
理论揭示DPO-TXO2的热激活延迟荧光（TADF）发光机制
热激活延迟荧光（TADF）材料是继荧光材料和贵金属磷光材料之后发展起来的第三代纯有机延迟荧光材料，其典型的特征是较小的单三重态能隙（ΔES-T）和温度正相关性。
2012年，日本九州大学的Chiahaya Adachi课题组首次报道外量子效率（EQE）超过20%的4CzIPN 分子[ ]，该材料的单线态和三线态能级差几乎为0，在室温下(298 K)的这样的热扰动下激子完全能够从三线态再回到单线态而发射荧光，因此命名为TADF（Thermally activated delayed fluorescence）。
当S1与T1的激发都是HOMO->LUMO特征，二者的能量差为2*K，K是HOMO与LUMO间的交换积分。随着HOMO与LUMO分离的增加，K会迅速减小。所以分离较大的时候，S1与T1 gap就较小，易于发生TADF需要的RISC。
为了保证高效的RISC，TADF材料需要具有较小的单三重态能隙，对应其HOMO/LUMO的有效分离，因此，TADF材料一般采用给体（D）−受体（A）、D−A−D的结构以不同的给受体作用实现HOMO/LUMO分离，同时兼顾其跃迁振子强度。
不同给受体的电子特性、三重态能级、结构刚性及扭曲程度等均均会影响材料的△EST、振子强度、态密度、激子寿命等，最终反映在材料的光物理性能和对应OLED器件的光电性能上。
本专题将以一个典型的TADF分子DPO-TXO2为例，介绍如何计算结构优化、频率、单点能、激发能、自旋轨道耦合等。同时介绍如何读取数据用于结果分析，帮助用户深入了解BDF软件的使用。
结构优化和频率计算
生成结构优化和频率输入文件
在Device Studio中导入准备的分子结构DPO-TXO2.xyz得到如图1.1-1所示界面，选中 Simulator → BDF → BDF，在弹出的界面中设置参数。
1.1-1
计算结构优化时计算类型选择Opt+Freq，方法、泛函、基组等选项用户可根据计算需要设置参数。例如Basic Settings面板设置为图1.1-2，SCF面板取消“Use MPEC+COSX”勾选（图1.1-3）、OPT 、Freq等面板的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。
1.1-2
1.1-3
生成的输入文件 bdf.inp参数部分如下：
$compass
Title
C39H28N2O4S
Geometry
C 3.86523 0.67704 0.08992
C 2.59676 1.19847 -0.21677
C 1.38236 0.46211 -0.14538
C 1.50274 -0.90633 0.05433
C 2.74673 -1.48909 0.32003
C 3.89360 -0.68925 0.41062
C 0.05129 1.23073 -0.21431
C -1.26041 0.42556 -0.14322
C -1.34326 -0.94957 0.03351
S 0.09634 -1.96093 -0.00226
C -2.49139 1.13510 -0.19404
C -3.75015 0.57230 0.07933
C -3.75167 -0.80689 0.33485
C -2.57699 -1.57032 0.24960
N 5.05789 1.50514 0.05106
N -4.95552 1.38707 0.07338
C 5.09111 2.89319 0.50297
C 6.28464 3.63010 0.39676
O 7.47953 3.08357 0.01235
C 7.47002 1.78524 -0.41733
C 6.30967 0.99832 -0.48773
C 8.72243 1.30821 -0.82591
C 8.84826 0.02519 -1.33737
C 7.70856 -0.74821 -1.50329
C 6.45512 -0.24869 -1.12362
C 4.01062 3.58921 1.07620
C 4.07062 4.96296 1.37442
C 5.24860 5.67030 1.18784
C 6.36600 4.99303 0.72541
C -6.19457 0.91553 -0.52385
C -7.33964 1.73082 -0.48834
O -7.34248 3.01488 -0.01720
C -6.17443 3.51502 0.46887
C -4.99409 2.75189 0.59422
C -6.34490 -0.31630 -1.18638
C -7.59189 -0.76699 -1.64640
C -8.71481 0.03325 -1.52666
C -8.57997 1.30489 -0.97531
C -6.24475 4.86124 0.86098
C -5.14195 5.49110 1.41274
C -3.98465 4.75621 1.61916
C -3.93157 3.39823 1.25512
O 0.11666 -2.61281 -1.29752
O 0.10373 -2.72112 1.23297
C 0.03300 2.06197 -1.51772
C 0.04308 2.16169 1.03932
H 2.54886 2.24058 -0.51595
H 2.82840 -2.56453 0.47286
H 4.82173 -1.17141 0.70878
H -2.46593 2.19212 -0.44272
H -4.67197 -1.32502 0.59460
H -2.63456 -2.65479 0.35810
H 9.59544 1.95023 -0.74373
H 9.82187 -0.35477 -1.63187
H 7.78471 -1.74349 -1.93391
H 5.60034 -0.87480 -1.35499
H 3.08415 3.09348 1.32929
H 3.19316 5.47421 1.76453
H 5.30763 6.72822 1.42899
H 7.31255 5.51704 0.61863
H -5.50297 -0.96874 -1.38412
H -7.67454 -1.75102 -2.10194
H -9.68032 -0.30389 -1.89032
H -9.43942 1.96697 -0.92291
H -7.17589 5.40700 0.73318
H -5.19606 6.53771 1.70383
H -3.11983 5.23203 2.07660
H -3.02635 2.86997 1.52459
H 0.02919 1.39736 -2.38952
H 0.89268 2.72961 -1.61468
H -0.84000 2.71525 -1.59635
H 0.04113 1.57168 1.96645
H -0.82598 2.82200 1.07532
H 0.91163 2.82447 1.08397
End Geometry
Basis
Def2-SVP
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
444
IOpt
3
Hess
final
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
PBE0
Molden
$end
$resp
Geom
$end
此时Device Studio图形界面如图1.1-4所示
1.1-4
备注
此处为保证结构优化和频率计算的条件相同，计算类型选择Opt+Freq，可以的单独做Opt计算或Freq计算。
BDF计算
在做BDF计算之前，需连接装有BDF的服务器，具体配置过程见鸿之微云操作指南。
连接好服务器，在做计算之前，用户可根据需要打开输入文件并查看文件中的参数设置是否合理，若不合理，则可选择直接在文件中编辑或重新生成，再进行BDF计算。
在图1.1-4所示的界面中，选中 bdf.inp → 右击 → Run。在弹出的界面导入相应的脚本，点击Run提交作业，如图1.1-5。
1.1-5
计算完成后点击下载按钮弹出计算结果界面如图1.1-6所示。选择.out结果文件，点击 Download下载。（提交作业操作为重复内容，在后面的计算中将不再赘述）
1.1-6
结构优化结果分析
右击下载后的out文件，选择Open with/Open containing folder即可查看结果文件。找到如下所示部分。
Force-RMS Force-Max Step-RMS Step-Max
Conv. tolerance : 0.2000E-03 0.3000E-03 0.8000E-03 0.1200E-02
Current values : 0.7369E-05 0.4013E-04 0.1843E-03 0.1041E-02
Geom. converge : Yes Yes Yes Yes
当Geom.converge的4个值均为YES时，证明结构优化收敛。上方和下方分别为收敛的分子结构笛卡尔坐标和内坐标。优化后的坐标信息可以作为初始结构用于后续计算。
检查频率，若不存在虚频证明结构已经优化到极小点。
单点能计算
生成单点能输入文件
将优化后的坐标导入Device Studio，名字改为DPO-TXO2-sp.xyz，此时图形界面如图1.2-1。
1.2-1
选中 Simulator → BDF → BDF，在弹出的界面中计算类型选择Single Point（默认值），方法、泛函、基组等选项用户可根据计算需要设置参数。例如泛函选PBE0，基组Def2-TZVP，其他参数仍为默认值，之后点击 Generate files 即可生成对应计算的输入文件。 生成的输入文件bdf.inp参数部分如下：
$compass
Title
C39H28N2O4S
Geometry
C 3.470732 -0.452949 0.333229
C 2.350276 -0.443126 -0.503378
C 1.255134 -1.275716 -0.258388
C 1.358849 -2.111496 0.851996
C 2.440432 -2.124490 1.711142
C 3.517727 -1.285828 1.451230
C -0.000048 -1.278142 -1.147435
C -1.255154 -1.275779 -0.258269
C -1.358725 -2.111574 0.852120
S 0.000118 -3.243604 1.269861
C -2.350358 -0.443230 -0.503130
C -3.470738 -0.453151 0.333573
C -3.517603 -1.286054 1.451551
C -2.440223 -2.124643 1.711370
N 4.564102 0.414026 0.042506
N -4.564206 0.413761 0.042962
C 4.451652 1.797113 0.288414
C 5.529066 2.638200 -0.032130
O 6.712474 2.137493 -0.580518
C 6.813862 0.759847 -0.795860
C 5.755871 -0.112762 -0.496962
C 7.999623 0.286590 -1.327509
C 8.161221 -1.076261 -1.582122
C 7.118160 -1.950624 -1.301513
C 5.922124 -1.471078 -0.764717
C 3.313452 2.367422 0.857787
C 3.242304 3.742953 1.084847
C 4.311909 4.564914 0.751035
C 5.460487 4.001069 0.193102
C -5.755562 -0.112971 -0.497448
C -6.813628 0.759568 -0.796285
O -6.712738 2.137080 -0.579852
C -5.529582 2.637766 -0.030885
C -4.452105 1.796731 0.289592
C -5.921333 -1.471159 -0.766141
C -7.116971 -1.950658 -1.303865
C -8.160095 -1.076369 -1.584473
C -7.998981 0.286358 -1.328883
C -5.461319 4.000541 0.194998
C -4.313011 4.564332 0.753554
C -3.243348 3.742416 1.087286
C -3.314166 2.366978 0.859540
O 0.000119 -4.563841 0.371547
O 0.000187 -3.483649 2.840945
C -0.000061 -2.561317 -2.024419
C -0.000112 -0.071391 -2.097897
H 2.353966 0.240214 -1.341805
H 2.400109 -2.768057 2.584222
H 4.382110 -1.260026 2.103052
H -2.354159 0.240153 -1.341521
H -4.381950 -1.260326 2.103422
H -2.399783 -2.768226 2.584432
H 8.781734 1.005474 -1.536628
H 9.092578 -1.440924 -1.998141
H 7.222431 -3.011204 -1.498846
H 5.108894 -2.153421 -0.550989
H 2.483350 1.726165 1.126879
H 2.346598 4.161499 1.529031
H 4.264620 5.633193 0.924336
H 6.321189 4.600814 -0.074686
H -5.108047 -2.153429 -0.552391
H -7.220889 -3.011140 -1.501914
H -9.091141 -1.440996 -2.001221
H -8.781175 1.005175 -1.537926
H -6.322045 4.600258 -0.072770
H -4.265977 5.632537 0.927382
H -2.347852 4.160920 1.531932
H -2.484014 1.725744 1.128541
H -0.000061 -3.470168 -1.414898
H 0.891657 -2.554225 -2.661972
H -0.891789 -2.554218 -2.661957
H -0.000071 0.880895 -1.555239
H -0.877870 -0.116199 -2.750591
H 0.877553 -0.116195 -2.750715
End Geometry
Basis
Def2-TZVP
Skeleton
Group
C(1)
$end
$xuanyuan
Direct
RS
0.33
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
CAM-B3LYP
MPEC+COSX
Molden
$end
BDF计算
同结构优化计算相同，连接好装有BDF的服务器后，选中 bdf.inp → 右击 → Run，检查脚本没有问题，点击Run提交作业。计算完成后点击下载按钮弹出计算结果，选择.out结果文件，点击 Download下载。
单点能结果分析
右击下载后的out文件，选择Open with/Open containing folder即可查看结果文件。找到E_tot为系统总能量，E_tot=E_ele + E_nn，本例中系统总能量为-2310.04883102 Hartree。E_ele是电子能量，E_nn是原子核排斥能，E_1e是单电子能量，E_ne 是原子核对电子的吸引能，E_kin 是电子动能，E_ee 是双电子能，E_xc 是交换相关能。
Final scf result
E_tot = -2311.25269871
E_ele = -7827.28555013
E_nn = 5516.03285142
E_1e = -14125.30142654
E_ne = -16425.97927385
E_kin = 2300.67784730
E_ee = 6514.27065120
E_xc = -216.25477479
Virial Theorem 2.004596
下方为轨道的占据情况，以及轨道能、HUMO-LOMO gap等信息。HOMO为-5.358 eV，LUMO为-1.962 eV，HOMO-LUMO gap为3.396 eV，Irrep为不可约表示，代表分子轨道对称性，本例中HOMO、LUMO不可约表示序号均为A。
[Final occupation pattern: ]
Irreps: A
detailed occupation for iden/irep: 1 1
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
...
Alpha HOMO energy: -0.24254414 au -6.59996455 eV Irrep: A
Alpha LUMO energy: -0.04116321 au -1.12010831 eV Irrep: A
HOMO-LUMO gap: 0.20138093 au 5.47985625 eV
最底部为Mulliken和Lowdin电荷布局、偶极矩信息。
[Mulliken Population Analysis]
Atomic charges:
1C -0.0009
2C -0.3029
3C 0.2227
4C -0.0143
5C -0.1228
6C -0.1890
7C 0.0046
8C 0.2227
9C -0.0150
10S 0.7787
11C -0.3023
12C -0.0022
13C -0.1888
14C -0.1223
15N -0.0121
16N -0.0121
17C 0.0563
...
[Lowdin Population Analysis]
Atomic charges:
1C -0.1574
2C -0.0592
3C -0.0682
4C -0.2154
5C -0.1050
6C -0.0869
7C -0.2270
8C -0.0682
9C -0.2154
10S 1.0012
11C -0.0591
12C -0.1574
13C -0.0869
14C -0.1050
...
[Dipole moment: Debye]
X Y Z |u|
Elec:-.3535E+01 0.8441E-03 -.1954E+01
Nucl:-.1254E-12 -.4210E-12 -.2935E-13
Totl: -3.5348 0.0008 -1.9541 4.0389
查看HOMO轨道图
为了更清楚的了解电子结构，往往需要做前线分子轨道分析。目前发布的版本BDF2022A中还无法实现数据的后处理，HOMO、LUMO轨道图可以用第三方软件Multiwfn+VMD渲染，需要用到scf.molden文件，软件的使用方法在量化论坛有专门的帖子可以学习，此文不做涉及。
HOMO轨道分布图
LUMO轨道分布图
得到的最高占据轨道(HOMO)与最低非占据轨道（LUMO）如图所示，由于两侧对称分布的吩恶嗪杂环是一个典型的给电子结构，而中心的磺酰化的四氢化萘是一个典型的吸电子的结构，因此整个分子是非常典型的D-A-D结构。可以看到HOMO轨道主要分布在两翼，LUMO轨道分布在中心，HOMO和LUMO轨道几乎没有重叠，符合TADF分子的电子结构特征。当然并不是所有HOMO和LUMO轨道分离的分子都具有TADF的光电特性，还需要满足S1和T1激发都是HOMO->LUMO轨道跃迁才行，因此我们可以进一步用BDF软件计算该分子的激发态电子结构。
激发态计算
生成激发态计算输入文件
读取优化好的结构做TDDFT计算，右键复制导入的优化后结构，命名为DPO-TXO2-td。计算类型选择TDDFT，方法、泛函、基组等选项用户可根据计算需要设置参数，前面的单点计算显示HOMO和LUMO轨道明显分离，对于这类具有明显D-A结构的分子，其激发态往往也会呈现电荷转移的特征，因此这儿我们选择最适合这类体系的范围分离泛函，如cam-B3LYP或者ω-B97xd。例如将Basic Settings面板按图1.3-1设置，TDDFT面板按图1.3-2设置，之后点击 Generate files 即可生成对应计算的输入文件。
1.3-1
1.3-2
生成的输入文件 bdf.inp参数部分如下：
$compass
Title
C39H28N2O4S
Geometry
C 3.56215000 -0.34631300 0.45361300
C 2.39970800 -0.43121500 -0.31807500
C 1.26327600 -1.11500900 0.12738900
C 1.35885600 -1.69579600 1.40258100
C 2.49771000 -1.60285400 2.19867100
C 3.61595700 -0.93278100 1.71813300
C 0.00021500 -1.24592200 -0.73874600
C -1.26297700 -1.11486500 0.12717900
C -1.35882900 -1.69562600 1.40235700
S -0.00010100 -2.61984500 2.07323100
C -2.39926700 -0.43096700 -0.31848800
C -3.56181100 -0.34590900 0.45301500
C -3.61589000 -0.93235000 1.71754500
C -2.49780200 -1.60255300 2.19826800
N 4.68577300 0.35565000 -0.05695800
N -4.68524700 0.35616500 -0.05781800
C 4.85522500 1.71734000 0.22325100
C 5.96987000 2.38879800 -0.31382300
O 6.88491700 1.74830700 -1.09915200
C 6.71947900 0.41903200 -1.36430000
C 5.62682300 -0.30753500 -0.85481400
C 7.67346300 -0.19823700 -2.15908800
C 7.56580700 -1.55645700 -2.46709500
C 6.49405000 -2.28575300 -1.96795600
C 5.53176100 -1.66610500 -1.16680600
C 3.96124200 2.44515800 1.01262100
C 4.17031100 3.80330200 1.26473400
C 5.27551600 4.45343400 0.73047600
C 6.17535900 3.73680700 -0.06194800
C -5.62705300 -0.30735400 -0.85450500
C -6.71928700 0.41938500 -1.36464300
O -6.88329900 1.74927200 -1.10167600
C -5.96897100 2.38946600 -0.31526500
C -4.85474100 1.71783400 0.22245400
C -5.53310000 -1.66639800 -1.16475900
C -6.49610200 -2.28636200 -1.96480800
C -7.56751800 -1.55693200 -2.46448300
C -7.67406700 -0.19823400 -2.15820200
C -6.17456800 3.73743100 -0.06324300
C -5.27514800 4.45388200 0.72982000
C -4.17031500 3.80359400 1.26465900
C -3.96122400 2.44545700 1.01253300
O -0.00015400 -3.96830000 1.50483700
O -0.00019500 -2.47109100 3.52665800
C 0.00020300 -2.64509100 -1.40495400
C 0.00034300 -0.20466000 -1.86117000
H 2.41118900 0.06372500 -1.28828500
H 2.48620300 -2.04935500 3.19547800
H 4.52498100 -0.84886800 2.31658900
H -2.41056900 0.06394100 -1.28871700
H -4.52499700 -0.84831800 2.31586200
H -2.48649200 -2.04903700 3.19508500
H 8.50056300 0.41098100 -2.52869800
H 8.32203900 -2.03354800 -3.09349600
H 6.39429300 -3.34933700 -2.19485600
H 4.69465500 -2.24580100 -0.77484200
H 3.09145400 1.94045700 1.43579900
H 3.45545900 4.34652000 1.88647900
H 5.44614600 5.51436800 0.92329600
H 7.05577600 4.20903800 -0.50207500
H -4.69625700 -2.24619000 -0.77237400
H -6.39717200 -3.35029700 -2.19042200
H -8.32431800 -2.03427800 -3.09000200
H -8.50081300 0.41112900 -2.52836600
H -7.05465600 4.20980200 -0.50387800
H -5.44580600 5.51480700 0.92266800
H -3.45579100 4.34667900 1.88689800
H -3.09175200 1.94062000 1.43619700
H 0.00013000 -3.45332000 -0.66309300
H 0.89243900 -2.75169300 -2.04060300
H -0.89196300 -2.75164000 -2.04071100
H 0.00033500 0.82736500 -1.47979800
H -0.87501100 -0.33812800 -2.51032400
H 0.87579000 -0.33816300 -2.51019000
End Geometry
Basis
Def2-TZVP
Skeleton
Group
C(1)
$end
$xuanyuan
Direct
RS
0.33
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
CAM-B3LYP
D3
MPEC+COSX
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
6
MPEC+COSX
Istore
1
$end
$tddft
NtoAnalyze
0
$end
$tddft
Imethod
1
Isf
1
Idiag
1
Iroot
6
MPEC+COSX
Istore
2
$end
$tddft
NtoAnalyze
0
$end
备注
Device studio中同名文件会被覆盖，输入文件默认名皆为bdf.inp。因此为避免数据被覆盖，我们每次计算需新建一个项目。
TDDFT面板Method一般建议选TDDFT，Multiplicity可选单重或三重或单重加三重。激发态数目默认计算6个，建议计算数目比实际想要的激发态数目多3个，如想计算10个态，此处可写13。
若想做NTO分析，TDDFT面板需勾选“Perform NTO Analyze”。
BDF计算
连接好装有BDF的服务器后，选中 bdf.inp → 右击 → Run，检查脚本没有问题，点击Run提交作业。计算完成后点击下载按钮弹出计算结果，选择.out结果文件，点击 Download下载。
激发态结果分析
激发能分析
右击下载后的out文件，选择Open with/Open containing folder即可查看结果文件。得到单重和三重激发能、振子强度、跃迁偶极矩等信息，isf=0为单重激发态信息，isf=1为三重激发态信息。
No. Pair ExSym ExEnergies Wavelengths f D<S^2> Domin
ant Excitations IPA Ova En-E1
1 A 2 A 3.4840 eV 355.86 nm 0.0023 0.0000 69.9% CV(0)
: A( 162 )-> A( 163 ) 5.584 0.164 0.0000
2 A 3 A 3.4902 eV 355.24 nm 0.0005 0.0000 69.3% CV(0)
: A( 161 )-> A( 163 ) 5.592 0.167 0.0061
3 A 4 A 3.8143 eV 325.05 nm 0.0003 0.0000 31.6% CV(0)
: A( 162 )-> A( 164 ) 6.182 0.482 0.3302
4 A 5 A 3.8152 eV 324.97 nm 0.0040 0.0000 31.0% CV(0)
: A( 161 )-> A( 164 ) 6.189 0.485 0.3312
5 A 6 A 4.1185 eV 301.05 nm 0.0163 0.0000 30.7% CV(0)
: A( 161 )-> A( 168 ) 6.944 0.583 0.6344
6 A 7 A 4.1229 eV 300.72 nm 0.1369 0.0000 30.8% CV(0)
: A( 162 )-> A( 168 ) 6.936 0.582 0.6388
*** Ground to excited state Transition electric dipole moments (Au) ***
State X Y Z Osc.
1 0.0003 -0.1642 0.0004 0.0023 0.0023
2 0.0579 -0.0010 0.0549 0.0005 0.0005
3 0.0019 0.0580 -0.0012 0.0003 0.0003
4 -0.1789 0.0007 0.1034 0.0040 0.0040
5 -0.0070 -0.4020 0.0039 0.0163 0.0163
6 1.0339 -0.0028 -0.5353 0.1369 0.1369
---------------------------------------------
---- End TD-DFT Calculations for isf = 0 ----
---------------------------------------------
...
No. Pair ExSym ExEnergies Wavelengths f D<S^2> Domin
ant Excitations IPA Ova En-E1
1 A 1 A 2.7522 eV 450.49 nm 0.0000 2.0000 25.3% CV(1)
: A( 162 )-> A( 167 ) 6.920 0.659 0.0000
2 A 2 A 2.7522 eV 450.49 nm 0.0000 2.0000 25.1% CV(1)
: A( 161 )-> A( 167 ) 6.928 0.659 0.0000
3 A 3 A 3.3404 eV 371.17 nm 0.0000 2.0000 33.1% CV(1)
: A( 154 )-> A( 163 ) 8.200 0.672 0.5882
4 A 4 A 3.3862 eV 366.15 nm 0.0000 2.0000 20.9% CV(1)
: A( 154 )-> A( 165 ) 8.983 0.649 0.6340
5 A 5 A 3.4620 eV 358.13 nm 0.0000 2.0000 50.3% CV(1)
: A( 162 )-> A( 163 ) 5.584 0.322 0.7098
6 A 6 A 3.4757 eV 356.72 nm 0.0000 2.0000 32.5% CV(1)
: A( 161 )-> A( 163 ) 5.592 0.466 0.7235
*** Ground to excited state Transition electric dipole moments (Au) ***
State X Y Z Osc.
1 0.0000 0.0000 0.0000 0.0000 0.0000
2 0.0000 0.0000 0.0000 0.0000 0.0000
3 0.0000 0.0000 0.0000 0.0000 0.0000
4 0.0000 0.0000 0.0000 0.0000 0.0000
5 0.0000 0.0000 0.0000 0.0000 0.0000
6 0.0000 0.0000 0.0000 0.0000 0.0000
---------------------------------------------
---- End TD-DFT Calculations for isf = 1 ----
---------------------------------------------
绘制成表格如下：
主要跃迁轨道 |
激发能/eV |
振子强度 |
贡献 |
偶极矩 |
波长/nm |
绝对重叠积分 |
|---|---|---|---|---|---|---|
A(162) -> A(163) |
3.4840 |
0.0023 |
69.9% |
0.1642 |
355.86 |
0.164 |
A(161) -> A(163) |
3.4902 |
0.0005 |
69.3% |
0.0798 |
355.24 |
0.167 |
A(162) -> A(164) |
3.8143 |
0.0003 |
31.6% |
0.0580 |
325.05 |
0.482 |
A(162) -> A(167) |
2.7522 |
0.0000 |
25.1% |
0.0000 |
450.49 |
0.659 |
A(161) -> A(167) |
2.7522 |
0.0000 |
25.3% |
0.0000 |
450.49 |
0.659 |
A(154) -> A(163) |
3.3404 |
0.0000 |
33.1% |
0.0000 |
371.17 |
0.672 |
表中依次给出激发态由低到高排序、多重度、不可约表示、占主要贡献的电子-空穴对激发、激发能、振子强度、跃迁轨道贡献占比、偶极矩、波长和绝对重叠积分。从表中我们能够看出，所研究的6个单激发态能级在2.7-4.0eV之间，分布较密集，其中前两个单重激发态波长在355nm左右，主要组分跃迁分别由HOMO→LUMO和HOMO-1→LUMO，表现出电荷转移特征。
文献报道的DPO-TXO2在溶剂环境下的能量最低吸收峰大约位于380nm左右，且随着溶剂极性的增大而红移。这主要是因为在极性越大的溶剂对极性越高的激发态稳定化程度也越高。n轨道极性最大，pi*次之，pi轨道极性最小。
计算显示DPO-TXO2分子的基态偶极矩是2.842 D，S1态的激发态偶极矩是19.4 D，显然激发态偶极矩明显大于基态偶极矩，因此激发态与溶剂环境的静电作用导致的能量降低比基态能量的降低更大，所以吸收光谱发生红移。
NTO分析
在激发态计算后，有时我们想更清楚的了解激发态跃迁的结果，此时可以做自然跃迁轨道（NTO）分析，对NTO分析的原理感兴趣的读者可以参考相关的博文（http://sobereva.com/91）。
假设我们对S1态感兴趣，可以单独对S1态做NTO分析。Basic Settings面板仍然按图1.3-1设置，TDDFT面板此时需要勾选“Perform NTO Analyze”，如图1.3-6所示。
1.3-6
备注
生成的输入文件第二个tddft模块也可手动修改为如下形式：
$tddft
NtoAnalyze
1 #对一个态NTO分析
1 #指定对第一个激发态做NTO分析
$end
计算结束后会产生nto1_1.molden格式文件，此文件中记录的已经不是scf.molden中MO轨道的信息了，而是NTO轨道信息，我们直接通过第三方软件Multiwfn主功能0并调整orbital info处理，得到的即为NTO轨道对的本征值与轨道图，软件的使用方法在科音论坛有专门的帖子可以学习，此文不做涉及。
DPO-TXO2分子的S1激发态的电子跃迁需要用两组NTO轨道才能较好地描述，下面是用VMD软件渲染出来的两组hole-particle轨道。
Hole1->particle1(73.26%)
Hole2->particle2(26.59%)
S1态NTO分析后可以看到占据轨道NTO1→非占据轨道NTO3的跃迁起主导，贡献为73.26%，占据轨道NTO2→非占据轨道NTO4贡献为26.59%。S1激发态的电子从两侧的吩恶嗪给电子基团跃迁到了中心的吸电子基团。
吸收光谱分析
对于激发态我们往往需要理论预测吸收谱，也就是将每个激发态按一定的半峰宽进行高斯展宽。在TDDFT计算正常结束后，我们需要进入终端用命令调用BDF安装路径下的plotspec.py脚本执行计算。若用户使用鸿之微云算力资源，进入命令端方式请查阅鸿之微云指南，此文不做赘述。 进入终断后，在目录下运行$BDFHOME/sbin/plotspec.py bdf.out，会产生两个文件，分别为bdf.stick.csv和bdf.spec.csv，前者包含所有激发态的吸收波长和摩尔消光系数，可以用来作棒状图，后者包含高斯展宽后的吸收谱（默认的展宽FWHM为0.5 eV），将bdf.spec.csv用第三方软件Origin作图如下：
说明位于基态的电子更容易吸收300nm波长的光发生跃迁。
激发态优化计算
生成激发态优化输入文件
导入优化好的基态结构，计算类型选择TDDFT-OPT，泛函PBE0，基组Def2-SVP，此时Basic Settings面板如图1.4-1所示，SCF面板同样消除“Use MPEC+COSX”勾选，如上图1.1-3。在优化S1态时，TDDFT面板的多重度选择Singlet，Target State为1，此时注意勾选“Calculate Dipole Moments of Target State”，如图1.4-2，OPT面板均保持默认值，点击 Generate files 即可生成对应计算的输入文件。
1.4-1
1.4-2
生成的输入文件 bdf.inp参数如下：
$compass
Title
C39H28N2O4S
Geometry
C 3.56215000 -0.34631300 0.45361300
C 2.39970800 -0.43121500 -0.31807500
C 1.26327600 -1.11500900 0.12738900
C 1.35885600 -1.69579600 1.40258100
C 2.49771000 -1.60285400 2.19867100
C 3.61595700 -0.93278100 1.71813300
C 0.00021500 -1.24592200 -0.73874600
C -1.26297700 -1.11486500 0.12717900
C -1.35882900 -1.69562600 1.40235700
S -0.00010100 -2.61984500 2.07323100
C -2.39926700 -0.43096700 -0.31848800
C -3.56181100 -0.34590900 0.45301500
C -3.61589000 -0.93235000 1.71754500
C -2.49780200 -1.60255300 2.19826800
N 4.68577300 0.35565000 -0.05695800
N -4.68524700 0.35616500 -0.05781800
C 4.85522500 1.71734000 0.22325100
C 5.96987000 2.38879800 -0.31382300
O 6.88491700 1.74830700 -1.09915200
C 6.71947900 0.41903200 -1.36430000
C 5.62682300 -0.30753500 -0.85481400
C 7.67346300 -0.19823700 -2.15908800
C 7.56580700 -1.55645700 -2.46709500
C 6.49405000 -2.28575300 -1.96795600
C 5.53176100 -1.66610500 -1.16680600
C 3.96124200 2.44515800 1.01262100
C 4.17031100 3.80330200 1.26473400
C 5.27551600 4.45343400 0.73047600
C 6.17535900 3.73680700 -0.06194800
C -5.62705300 -0.30735400 -0.85450500
C -6.71928700 0.41938500 -1.36464300
O -6.88329900 1.74927200 -1.10167600
C -5.96897100 2.38946600 -0.31526500
C -4.85474100 1.71783400 0.22245400
C -5.53310000 -1.66639800 -1.16475900
C -6.49610200 -2.28636200 -1.96480800
C -7.56751800 -1.55693200 -2.46448300
C -7.67406700 -0.19823400 -2.15820200
C -6.17456800 3.73743100 -0.06324300
C -5.27514800 4.45388200 0.72982000
C -4.17031500 3.80359400 1.26465900
C -3.96122400 2.44545700 1.01253300
O -0.00015400 -3.96830000 1.50483700
O -0.00019500 -2.47109100 3.52665800
C 0.00020300 -2.64509100 -1.40495400
C 0.00034300 -0.20466000 -1.86117000
H 2.41118900 0.06372500 -1.28828500
H 2.48620300 -2.04935500 3.19547800
H 4.52498100 -0.84886800 2.31658900
H -2.41056900 0.06394100 -1.28871700
H -4.52499700 -0.84831800 2.31586200
H -2.48649200 -2.04903700 3.19508500
H 8.50056300 0.41098100 -2.52869800
H 8.32203900 -2.03354800 -3.09349600
H 6.39429300 -3.34933700 -2.19485600
H 4.69465500 -2.24580100 -0.77484200
H 3.09145400 1.94045700 1.43579900
H 3.45545900 4.34652000 1.88647900
H 5.44614600 5.51436800 0.92329600
H 7.05577600 4.20903800 -0.50207500
H -4.69625700 -2.24619000 -0.77237400
H -6.39717200 -3.35029700 -2.19042200
H -8.32431800 -2.03427800 -3.09000200
H -8.50081300 0.41112900 -2.52836600
H -7.05465600 4.20980200 -0.50387800
H -5.44580600 5.51480700 0.92266800
H -3.45579100 4.34667900 1.88689800
H -3.09175200 1.94062000 1.43619700
H 0.00013000 -3.45332000 -0.66309300
H 0.89243900 -2.75169300 -2.04060300
H -0.89196300 -2.75164000 -2.04071100
H 0.00033500 0.82736500 -1.47979800
H -0.87501100 -0.33812800 -2.51032400
H 0.87579000 -0.33816300 -2.51019000
End Geometry
Basis
Def2-TZVP
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
444
IOpt
3
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
PBE0
D3
Molden
$end
$tddft
Imethod
1
Isf
0
Ialda
4
Idiag
1
Iroot
4
MPEC+COSX
Istore
1
$end
$resp
Geom
Method
2
Nfiles
1
Iroot
1
$end
备注
对T1态优化时，将TDDFT面板的多重度改为Triplet，其余参数同S1优化。
BDF计算
连接好装有BDF的服务器后，选中 bdf.inp → 右击 → Run，检查脚本没有问题，点击Run提交作业。计算完成后点击下载按钮弹出计算结果，选择.out结果文件，点击 Download下载。
激发态优化结果分析 右击下载后的out文件，选择Open with/Open containing folder即可查看结果文件。类似基态结构优化，当Geom.converge的4个值均为YES时，证明结构优化收敛，如上图1.1-8。将优化后的T1与S1能量相减，粗略计算ΔEST=2.425 eV。
自旋轨道耦合计算
生成自旋轨道耦合输入文件
对优化好的结构做SOC计算。计算类型选择TDDFT-SOC，哈密顿选择sf-x2c，方法、泛函可根据计算需要设置，基组选择相对论基组，例如cc-pVDZ-DK，此时Basic Settings面板如图1.5-1设置，SCF、TDDFT面板仍为默认值，之后点击 Generate files 即可生成对应计算的输入文件。
1.5-1
生成的输入文件 bdf.inp参数如下：
$compass
Title
C39H28N2O4S
Geometry
C 3.56214999 -0.34631300 0.45361300
C 2.39970799 -0.43121500 -0.31807500
C 1.26327600 -1.11500899 0.12738900
C 1.35885600 -1.69579600 1.40258100
C 2.49771000 -1.60285400 2.19867100
C 3.61595699 -0.93278100 1.71813299
C 0.00021500 -1.24592199 -0.73874600
C -1.26297700 -1.11486500 0.12717899
C -1.35882900 -1.69562600 1.40235700
S -0.00010100 -2.61984500 2.07323099
C -2.39926700 -0.43096700 -0.31848800
C -3.56181100 -0.34590900 0.45301500
C -3.61588999 -0.93235000 1.71754500
C -2.49780200 -1.60255299 2.19826800
N 4.68577300 0.35565000 -0.05695800
N -4.68524700 0.35616500 -0.05781800
C 4.85522499 1.71734000 0.22325100
C 5.96987000 2.38879800 -0.31382300
O 6.88491699 1.74830700 -1.09915199
C 6.71947899 0.41903200 -1.36430000
C 5.62682299 -0.30753500 -0.85481400
C 7.67346299 -0.19823700 -2.15908800
C 7.56580700 -1.55645700 -2.46709500
C 6.49404999 -2.28575300 -1.96795600
C 5.53176100 -1.66610499 -1.16680600
C 3.96124200 2.44515800 1.01262099
C 4.17031099 3.80330200 1.26473400
C 5.27551599 4.45343399 0.73047600
C 6.17535900 3.73680700 -0.06194800
C -5.62705300 -0.30735400 -0.85450500
C -6.71928699 0.41938500 -1.36464300
O -6.88329900 1.74927200 -1.10167600
C -5.96897099 2.38946600 -0.31526500
C -4.85474099 1.71783400 0.22245400
C -5.53310000 -1.66639800 -1.16475900
C -6.49610199 -2.28636200 -1.96480800
C -7.56751799 -1.55693200 -2.46448300
C -7.67406700 -0.19823400 -2.15820200
C -6.17456799 3.73743100 -0.06324299
C -5.27514799 4.45388200 0.72982000
C -4.17031500 3.80359399 1.26465899
C -3.96122400 2.44545700 1.01253299
O -0.00015400 -3.96830000 1.50483700
O -0.00019500 -2.47109099 3.52665799
C 0.00020300 -2.64509099 -1.40495400
C 0.00034300 -0.20466000 -1.86117000
H 2.41118899 0.06372500 -1.28828499
H 2.48620300 -2.04935499 3.19547800
H 4.52498100 -0.84886800 2.31658900
H -2.41056900 0.06394100 -1.28871699
H -4.52499699 -0.84831800 2.31586200
H -2.48649200 -2.04903700 3.19508500
H 8.50056299 0.41098100 -2.52869799
H 8.32203900 -2.03354800 -3.09349600
H 6.39429300 -3.34933699 -2.19485600
H 4.69465500 -2.24580100 -0.77484200
H 3.09145400 1.94045700 1.43579900
H 3.45545899 4.34651999 1.88647900
H 5.44614599 5.51436800 0.92329600
H 7.05577599 4.20903799 -0.50207500
H -4.69625700 -2.24618999 -0.77237400
H -6.39717200 -3.35029699 -2.19042199
H -8.32431799 -2.03427800 -3.09000200
H -8.50081300 0.41112900 -2.52836600
H -7.05465599 4.20980199 -0.50387800
H -5.44580600 5.51480700 0.92266800
H -3.45579100 4.34667899 1.88689800
H -3.09175200 1.94062000 1.43619699
H 0.00012999 -3.45332000 -0.66309300
H 0.89243900 -2.75169300 -2.04060299
H -0.89196300 -2.75164000 -2.04071099
H 0.00033500 0.82736500 -1.47979799
H -0.87501100 -0.33812800 -2.51032400
H 0.87579000 -0.33816300 -2.51019000
End Geometry
Basis
cc-pVDZ-DK
Skeleton
Group
C(1)
$end
$xuanyuan
Heff
21
Hsoc
2
Direct
RS
0.33
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
CAM-B3LYP
D3
MPEC+COSX
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
6
MPEC+COSX
Istore
1
$end
$tddft
Imethod
1
Isf
1
Idiag
1
Iroot
6
MPEC+COSX
Istore
2
$end
$tddft
Isoc
2
Nfiles
2
Imatsoc
-1
Imatrsf
-1
Imatrso
-1
$end
BDF计算
连接好装有BDF的服务器后，选中 bdf.inp → 右击 → Run，检查脚本没有问题，点击Run提交作业。计算完成后点击下载按钮弹出计算结果，选择.out结果文件，点击 Download下载。
耦合矩阵元结果分析
右击下载后的out文件，选择Open with/Open containing folder即可查看结果文件。在Print selected matrix elements of [Hsoc]部分打印耦合矩阵元信息。
SocPairNo. = 8 SOCmat = < 0 0 0 |Hso| 2 1 1 > Dim = 1 3
mi/mj ReHso(au) cm^-1 ImHso(au) cm^-1
0.0 -1.0 -0.0000018883 -0.4144393040 -0.0000012470 -0.2736747987
0.0 0.0 0.0000000000 0.0000000000 -0.0000076582 -1.6807798007
0.0 1.0 -0.0000018883 -0.4144393040 0.0000012470 0.2736747987
SocPairNo. = 9 SOCmat = < 0 0 0 |Hso| 2 1 2 > Dim = 1 3
mi/mj ReHso(au) cm^-1 ImHso(au) cm^-1
0.0 -1.0 0.0000038630 0.8478326909 -0.0000006073 -0.1332932016
0.0 0.0 0.0000000000 0.0000000000 -0.0000037537 -0.8238381363
0.0 1.0 0.0000038630 0.8478326909 0.0000006073 0.1332932016
...
绘制表格
矩阵元的模/cm^-1 |
T1 |
T2 |
|---|---|---|
S0 |
1.822 |
1.467 |
S1 |
0.522 |
0.842 |
计算得到S0态与T1态旋轨耦合1.822 cm^-1 ，如果能隙足够小，就会引起系间窜越的发生。
BDF-QM/MM案例教程一
本专题将介绍一种量子化学与分子力学结合的方法（QM/MM方法），该方法既包括量子化学的精确性，又利用分子力学的高效性，其基本思想是用量子力学处理感兴趣的中心，其余部分用经典分子力学来处理。
本章节以一个典型的没食子酸分子（Gallic Acid，GA）为例，介绍输入文件准备，QM/MM 单点计算，QM/MM结构优化，QM/MM激发态计算。其中，BDF程序主要完成量子化学计算部分，其余部分由BDF开发成员修改的pDynamo2.0程序包完成。同时介绍如何读取数据用于结果分析，帮助用户深入了解BDF软件的使用。
输入文件准备
一般来说，QM/MM计算之前，需要对目标体系进行分子动力学模拟，得到适合的初始构象。当采用PDB、MOL2或xyz文件作为输入时，pDynamo2.0程序包仅支持OPLS力场，对于小分子和非标准氨基酸力场参数不全，不推荐使用。建议优先采用Amber程序，通过拓扑文件输入力场参数。以Amber为例，从动力学模拟轨迹提取感兴趣的结构存储于.crd文件中，与对应的参数/拓扑文件.prmtop一起可以作为QM/MM计算的起始点。Python脚本如下：
from pBabel import AmberCrdFile_ToCoordinates3, AmberTopologyFile_ToSystem
# 读取输入信息
molecule = AmberTopologyFile_ToSystem(Topfile)
molecule.coordinates3 = AmberCrdFile_ToCoordinates3(CRDfile)
其中，需要提前安装好AmberTools，python2.0版本，并正确设置好 AMBERHOME 和 PDYNAMO 环境变量，关于如何将GallicAcid.pdb初始结构文件(图1，晶胞为2*1*1)生成使用AmberTools21程序相对应的坐标文件GallicAcid.crd和参数/拓扑文件GallicAcid.prmop的方法如下：
运行antechamber程序将Pdb文件转化为mol2文件：
antechamber -i GallicAcid.pdb -fi pdb -o GallicAcid.mol2 -fo mol2 -j 5 -at amber -dr no - -i 指定输入文件 - -fi 指定输入文件类型 - -o 指定输出文件 - -fo 指定输出文件类型 - -j 匹配原子类型和键类型 - -at 定义原子类型
运行parmchk2程序生成对应体系的力场参数文件
parmchk2 -i GallicAcid.mol2 -f mol2 -o GallicAcid.frcmod
运行tleap程序构建系统拓扑并为分子定义力场参数步骤如下:
使用 tleap 命令启动tleap程序
-I: Adding /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/prep to search path.
-I: Adding /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/lib to search path.
-I: Adding /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/parm to search path.
-I: Adding /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/cmd to search path.
Welcome to LEaP!
(no leaprc in search path)
>
确定并加载体系力场：source leaprc.gaff(此为GAFF力场)
> source leaprc.gaff
----- Source: /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/cmd/leaprc.gaff
----- Source of /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/cmd/leaprc.gaff done
Log file: ./leap.log
Loading parameters: /es01/jinan/hzw001/home/hzw1011/anaconda3/envs/AmberTools21/dat/leap/parm/gaff.dat
Reading title:
AMBER General Force Field for organic molecules (Version 1.81, May 2017)
>
调入配体mol2文件：GA = loadmol2 GallicAcid.mol2
> GA = loadmol2 GallicAcid.mol2
Loading Mol2 file: ./GallicAcid.mol2
Reading MOLECULE named WAT
>
检查导入的结构是否准确或缺失参数：check GA
调入体系分子的模板，并补全库文件中缺失的参数:loadamberparams GallicAcid.frcmod
准备生成的Sustiva库文件：saveoff GA GallicAcid.lib
修改生成的Sustiva库文件并调入该文件：loadoff GallicAcid.lib
> loadoff GallicAcid.lib
Loading library: ./GallicAcid.lib
保存.crd和.prmop文件：saveamberparm GA GallicAcid.prmtop GallicAcid.crd
> saveamberparm GA GallicAcid.prmtop GallicAcid.crd
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 112 improper torsions applied
Building H-Bond parameters.
Incorporating Non-Bonded adjustments.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:
res total affected
WAT 1
)
(no restraints)
>
退出tleap程序：quit
分子动力学模拟
此处采用amber软件进行分子动力学模拟，首先对体系进行能量最小化模拟，输入文件min.in如下：
Initial minimisation of GallicAcid complex
&cntrl
imin=1, maxcyc=200, ncyc=50,
cut=16, ntb=0, igb=1,
&end
- imin=1：运行能量最小化
- maxcyc=200：能量最小化的最大循环数
- ncyc=50：最初的0到ncyc循环使用最速下降算法, 此后的ncyc到maxcyc循环切换到共轭梯度算法
- cut=16：以埃为单位的非键截断距离
- ntb=0：关闭周期性边界条件
- igb=1：Born模型
使用如下命令运行能量最小化:
sander -O -i min.in -o GallicAcid_min.out -p GallicAcid.prmtop -c GallicAcid.crd -r GallicAcid_min.rst &
其中GallicAcid_min.rst为输出包含坐标和速度的重启文件
接下来利用最小化模拟得到的重启文件升温系统，从而完成分子动力学模拟，输入文件md.in如下：
Initial MD equilibration
&cntrl
imin=0, irest=0,
nstlim=1000,dt=0.001, ntc=1,
ntpr=20, ntwx=20,
cut=16, ntb=0, igb=1,
ntt=3, gamma_ln=1.0,
tempi=0.0, temp0=300.0,
&end
imin=0：进行分子动力学(MD)
irest=0：读取先前保存的重新启动文件读取坐标和速度
nstlim=1000：运行的MD步数
dt=0.001：时间步长（单位：ps）
ntc=1：不启用SHAKE约束
ntpr=20：每ntpr步输出能量信息mdout一次
ntwx=20：每ntwx步输出Amber轨迹文件mdcrd一次
ntt=3：Langevin恒温器控制温度
gamma_ln=1.0：Langevin恒温器的碰撞频率
tempi=0.0：模拟的初始温度
temp0=300.0：模拟的最终温度
使用如下命令运行分子动力学模拟:
sander -O -i md.in -o md.out -p GallicAcid.prmtop -c GallicAcid_min.rst -r GallicAcid_md.rst -x GallicAcid_md.mdcrd -inf GallicAcid_md.mdinfo
其中GallicAcid_md.mdcrd文件即为MD模拟的轨迹文件，可借助VMD软件进行可视化显示分子结构，并从动力学模拟轨迹提取感兴趣的结构存储于.crd文件中。
QM/MM 总能量计算
分子动力学模拟后提取文件为GallicAcid.prmtop， GallicAcid.crd，可对体系进行全量子化学总能量计算，python代码如下：
import glob, math, os
from pBabel import AmberCrdFile_ToCoordinates3, AmberTopologyFile_ToSystem
from pCore import logFile
from pMolecule import QCModelBDF, System
# 读取水盒子坐标和拓扑信息
molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 定义能量计算模式，此处为全体系密度泛函计算，可以定义方法和基组，分别为GB3LYP和6-31g，
model = QCModelBDF("GB3LYP:6-31g")
molecule.DefineQCModel(model)
molecule.Summary() #输出体系计算设置信息
# 计算总能量
energy = molecule.Energy()
除了可以用全量子化学QM计算体系总能量，也可对感兴趣的分子进行QM/MM计算（本例为指定第五个分子用QM方法计算），QM/MM组合能量计算python脚本如下：
import glob, math, os
from pBabel import AmberCrdFile_ToCoordinates3, AmberTopologyFile_ToSystem
from pCore import logFile, Selection
from pMolecule import NBModelORCA, QCModelBDF, System
# 定义能量计算模式
nbModel = NBModelORCA() #处理QM和MM区相互作用
qcModel = QCModelBDF("GB3LYP:6-31g")
# 读取体系坐标和拓扑信息
molecule = AmberTopologyFile_ToSystem("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 关闭体系对称性
molecule.DefineSymmetry(crystalClass = None) # QM/MM方法不支持使用周期性边界条件，故关闭周期性边界条件
# 指定QM区
qm_area = Selection.FromIterable(range (72, 90)) # 指定第五个分子用QM方法计算，其中(72, 90)指明原子列表索引值为72，73，74…..87,88,89，该值=原子序数-1
# 定义能量计算模式
molecule.DefineQCModel (qcModel, qcSelection = qm_area)
molecule.DefineNBModel (nbModel)
molecule.Summary()
# 计算总能量
energy = molecule.Energy()
QM/MM模拟的输出总结了MM部分，QM部分，QM区和MM区相互作用部分的计算细节如下：
----------------------------------- Summary for MM Model "AMBER" -----------------------------------
LJ 1-4 Scaling = 0.500 El. 1-4 Scaling = 0.833
Number of MM Atoms = 288 Number of MM Atom Types = 6
Number of Inactive MM Atoms = 18 Total MM Charge = 0.00
Harmonic Bond Terms = 288 Harmonic Bond Parameters = 7
Harmonic Bond Inactive = 18 Harmonic Angle Terms = 400
Harmonic Angle Parameters = 9 Harmonic Angle Inactive = 25
Fourier Dihedral Terms = 592 Fourier Dihedral Parameters = 5
Fourier Dihedral Inactive = 37 Fourier Improper Terms = 112
Fourier Improper Parameters = 1 Fourier Improper Inactive = 7
Exclusions = 1216 1-4 Interactions = 528
LJ Parameters Form = AMBER LJ Parameters Types = 5
1-4 Lennard-Jones Form = AMBER 1-4 Lennard-Jones Types = 5
----------------------------------------------------------------------------------------------------
------------------- Summary for QC Model "BDF:GB3LYP:STO-3g" -------------------
Number of QC Atoms = 18 Boundary Atoms = 0
Nuclear Charge = 88 Orbital Functions = 0
Fitting Functions = 0 Energy Base Line = 0.00000
--------------------------------------------------------------------------------
----------------------------- ORCA NB Model Summary ----------------------------
El. 1-4 Scaling = 0.833333 QC/MM Coupling = RC Coupling
--------------------------------------------------------------------------------
------------------------------- Sequence Summary -------------------------------
Number of Atoms = 288 Number of Components = 16
Number of Entities = 1 Number of Linear Polymers = 0
Number of Links = 0 Number of Variants = 0
--------------------------------------------------------------------------------
输出体系总能量信息以及各部分能量贡献如下：
--------------------------------- Summary of Energy Terms --------------------------------
Potential Energy = -1671893.4718 RMS Gradient = None
Harmonic Bond = 1743.3211 Harmonic Angle = 124.9878
Fourier Dihedral = 269.8417 Fourier Improper = 0.1346
MM/MM LJ = -138.0022 MM/MM 1-4 LJ = 474.4044
QC/MM LJ = -42.2271 BDF QC = -1674325.9320
------------------------------------------------------------------------------------------
QM/MM 结构优化
QM/MM几何构型优化计算的python脚本如下：
import glob, math, os.path
from pBabel import AmberCrdFile_ToCoordinates3, \
AmberTopologyFile_ToSystem , \
SystemGeometryTrajectory , \
AmberCrdFile_FromSystem , \
PDBFile_FromSystem , \
XYZFile_FromSystem
from pCore import Clone, logFile, Selection
from pMolecule import NBModelORCA, QCModelBDF, System
from pMoleculeScripts import ConjugateGradientMinimize_SystemGeometry, \
FIREMinimize_SystemGeometry , \
LBFGSMinimize_SystemGeometry , \
SteepestDescentMinimize_SystemGeometry
# 定义结构优化接口
def opt_ConjugateGradientMinimize(molecule, selection):
molecule.DefineFixedAtoms(selection) #固定原子
#定义优化方法
ConjugateGradientMinimize_SystemGeometry(
molecule,
maximumIterations = 40, # 最大优化步数
rmsGradientTolerance = 0.1, #优化收敛控制
trajectories = [(trajectory, 1)]
) # 定义轨迹保存频率
# 定义能量计算模式
nbModel = NBModelORCA()
qcModel = QCModelBDF("GB3LYP:6-31g")
# 读取体系坐标和拓扑信息
molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 关闭体系对称性
molecule.DefineSymmetry(crystalClass = None) # QM/MM方法不支持使用周期性边界条件
#. Define Atoms List
natoms = len(molecule.atoms) # 系统中总原子数
qm_list = range(72, 90) # QM 区原子
activate_list = range(126, 144) + range (144, 162) # MM区活性原子（优化中可以移动）
#定义MM区原子
mm_list = range (natoms)
for i in qm_list:
mm_list.remove(i) # MM 删除QM原子
mm_inactivate_list = mm_list[:]
for i in activate_list :
mm_inactivate_list.remove(i)
# 输入QM原子
qmmmtest_qc = Selection.FromIterable(qm_list)
# 定义各选择区
selection_qm_mm_inactivate = Selection.FromIterable(qm_list + mm_inactivate_list)
selection_mm = Selection.FromIterable(mm_list)
selection_mm_inactivate = Selection.FromIterable(mm_inactivate_list)
# . Define the energy model.
molecule.DefineQCModel(qcModel, qcSelection = qmmmtest_qc)
molecule.DefineNBModel(nbModel)
molecule.Summary()
#计算优化开始时总能量
eStart = molecule.Energy()
#定义输出文件目录名
outlabel = 'opt_watbox_bdf'
if os.path.exists(outlabel):
pass
else:
os.mkdir (outlabel)
outlabel = outlabel + '/' + outlabel
# 定义输出轨迹
trajectory = SystemGeometryTrajectory (outlabel + ".trj" , molecule, mode = "w")
# 开始第一阶段优化
# 定义优化两步
iterations = 2
# 顺次固定QM区和MM区进行优化
for i in range(iterations):
opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM区优化
opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM区优化
# 开始第二阶段优化
# QM区和MM区同时优化
opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)
#输出优化后总能量
eStop = molecule.Energy()
#保存优化坐标，可以为xyz/crd/pdb等。
XYZFile_FromSystem(outlabel + ".xyz", molecule)
AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)
PDBFile_FromSystem(outlabel + ".pdb" , molecule)
输出体系收敛信息如下（此处仅展示前20步优化收敛结果）：
----------------------------------------------------------------------------------------------------------------
Iteration Function RMS Gradient Max. |Grad.| RMS Disp. Max. |Disp.|
----------------------------------------------------------------------------------------------------------------
0 I -1696839.69778731 2.46510318 9.94250232 0.00785674 0.03168860
2 L1s -1696839.82030342 1.38615730 5.83254788 0.00043873 0.00126431
4 L1s -1696839.90971371 1.41241184 5.29242524 0.00067556 0.00172485
6 L0s -1696840.01109863 1.41344485 4.70119338 0.00090773 0.00265969
8 L1s -1696840.09635696 1.44964059 5.72496661 0.00108731 0.00328490
10 L1s -1696840.17289698 1.28607709 4.73666387 0.00108469 0.00354577
12 L1s -1696840.23841524 1.03217304 3.00441004 0.00081945 0.00267931
14 L1s -1696840.30741088 1.40349698 5.22220965 0.00162080 0.00519590
16 L1s -1696840.43546466 1.32604042 4.51175225 0.00158796 0.00455431
18 L0s -1696840.52547251 1.27123125 4.20616166 0.00158796 0.00428040
20 L0s -1696840.60265453 1.08553355 3.12355616 0.00158796 0.00470223
----------------------------------------------------------------------------------------------------------------
输出体系总能量信息如下：
--------------------------------- Summary of Energy Terms --------------------------------
Potential Energy = -1696841.6016 RMS Gradient = None
Harmonic Bond = 3.0295 Harmonic Angle = 3.6222
Fourier Dihedral = 32.0917 Fourier Improper = 0.0040
MM/MM LJ = -69.3255 MM/MM 1-4 LJ = 43.9528
QC/MM LJ = -47.2706 BDF QC = -1696807.7057
------------------------------------------------------------------------------------------
备注
QM/MM几何构型优化一般不容易收敛，在实际操作中需要的技巧较多。常见的有，固定MM区，优化QM区；然后固定QM区优化MM区。如此往复循环几次后，再同时优化QM区和MM区。优化是否收敛，和QM区的选择及QM/MM边界是否有带电较多的原子等关系很大。为了加速优化，可以在计算时固定MM区，仅选择离QM区较近的合适区域，作为活性区域，在优化中坐标可以变化。
QM/MM 激发态计算
基于上一步的QM/MM几何构型优化，继而即可将MM区活性原子添加到QM区进行QM/MM-TDDFT计算，完整的代码如下:
import glob, math, os.path
from pBabel import AmberCrdFile_ToCoordinates3, \
AmberTopologyFile_ToSystem , \
SystemGeometryTrajectory , \
AmberCrdFile_FromSystem , \
PDBFile_FromSystem , \
XYZFile_FromSystem
from pCore import Clone, logFile, Selection
from pMolecule import NBModelORCA, QCModelBDF, System
from pMoleculeScripts import ConjugateGradientMinimize_SystemGeometry, \
FIREMinimize_SystemGeometry , \
LBFGSMinimize_SystemGeometry , \
SteepestDescentMinimize_SystemGeometry
# 定义结构优化接口
def opt_ConjugateGradientMinimize(molecule, selection):
molecule.DefineFixedAtoms(selection) #固定原子
#定义优化方法
ConjugateGradientMinimize_SystemGeometry(
molecule,
maximumIterations = 40, # 最大优化步数
rmsGradientTolerance = 0.1, #优化收敛控制
trajectories = [(trajectory, 1)]
) # 定义轨迹保存频率
# 定义能量计算模式
nbModel = NBModelORCA()
qcModel = QCModelBDF("GB3LYP:6-31g")
# 读取体系坐标和拓扑信息
molecule = AmberTopologyFile_ToSystem ("GallicAcid.prmtop")
molecule.coordinates3 = AmberCrdFile_ToCoordinates3("GallicAcid.crd")
# 关闭体系对称性
molecule.DefineSymmetry(crystalClass = None) # QM/MM方法不支持使用周期性边界条件
#. Define Atoms List
natoms = len(molecule.atoms) # 系统中总原子数
qm_list = range(72, 90) # QM 区原子
activate_list = range(126, 144) + range (144, 162) # MM区活性原子（优化中可以移动）
#定义MM区原子
mm_list = range (natoms)
for i in qm_list:
mm_list.remove(i) # MM 删除QM原子
mm_inactivate_list = mm_list[:]
for i in activate_list :
mm_inactivate_list.remove(i)
# 输入QM原子
qmmmtest_qc = Selection.FromIterable(qm_list)
# 定义各选择区
selection_qm_mm_inactivate = Selection.FromIterable(qm_list + mm_inactivate_list)
selection_mm = Selection.FromIterable(mm_list)
selection_mm_inactivate = Selection.FromIterable(mm_inactivate_list)
# . Define the energy model.
molecule.DefineQCModel(qcModel, qcSelection = qmmmtest_qc)
molecule.DefineNBModel(nbModel)
molecule.Summary()
#计算优化开始时总能量
eStart = molecule.Energy()
#定义输出文件目录名
outlabel = 'opt_watbox_bdf'
if os.path.exists(outlabel):
pass
else:
os.mkdir (outlabel)
outlabel = outlabel + '/' + outlabel
# 定义输出轨迹
trajectory = SystemGeometryTrajectory (outlabel + ".trj" , molecule, mode = "w")
# 开始第一阶段优化
# 定义优化两步
iterations = 2
# 顺次固定QM区和MM区进行优化
for i in range(iterations):
opt_ConjugateGradientMinimize(molecule, selection_qm_mm_inactivate) #固定QM区优化
opt_ConjugateGradientMinimize(molecule, selection_mm) #固定MM区优化
# 开始第二阶段优化
# QM区和MM区同时优化
opt_ConjugateGradientMinimize(molecule, selection_mm_inactivate)
#输出优化后总能量
eStop = molecule.Energy()
#保存优化坐标，可以为xyz/crd/pdb等。
XYZFile_FromSystem(outlabel + ".xyz", molecule)
AmberCrdFile_FromSystem(outlabel + ".crd" , molecule)
PDBFile_FromSystem(outlabel + ".pdb" , molecule)
# TDDFT计算
qcModel = QCModelBDF_template ( )
qcModel.UseTemplate (template = 'head_bdf_nosymm.inp' )
tdtest = Selection.FromIterable ( qm_list + activate_list )
# . Define the energy model.
molecule.DefineQCModel ( qcModel, qcSelection = tdtest )
molecule.DefineNBModel ( nbModel )
molecule.Summary ( )
# . Calculate
energy = molecule.Energy ( )
输出体系总能量信息如下：
--------------------------------- Summary of Energy Terms --------------------------------
Potential Energy = -5088333.3818 RMS Gradient = None
Harmonic Bond = 0.0000 Harmonic Angle = 0.0000
Fourier Dihedral = 0.0000 Fourier Improper = 0.0000
QC/MM LJ = -112.3207 BDF QC = -5088221.0611
------------------------------------------------------------------------------------------
同时生成.log结果文件，和普通的激发态计算一样，可以看到振子强度，激发能，激发态的总能量等信息:
No. 1 w= 4.7116 eV -1937.8276358207 a.u. f= 0.0217 D<Pab>= 0.0000 Ova= 0.6704
CV(0): A( 129 )-> A( 135 ) c_i: 0.7254 Per: 52.6% IPA: 7.721 eV Oai: 0.6606
CV(0): A( 129 )-> A( 138 ) c_i: 0.2292 Per: 5.3% IPA: 9.104 eV Oai: 0.8139
CV(0): A( 132 )-> A( 135 ) c_i: 0.4722 Per: 22.3% IPA: 7.562 eV Oai: 0.6924
CV(0): A( 132 )-> A( 138 ) c_i: -0.4062 Per: 16.5% IPA: 8.946 eV Oai: 0.6542
随后还打印了跃迁偶极矩:
*** Ground to excited state Transition electric dipole moments (Au) ***
State X Y Z Osc.
1 0.0959 0.1531 0.3937 0.0217 0.0217
2 0.0632 -0.1286 0.3984 0.0207 0.0207
3 -0.0797 -0.2409 0.4272 0.0287 0.0287
4 0.0384 -0.0172 -0.0189 0.0003 0.0003
5 1.1981 0.8618 -0.1305 0.2751 0.2751
QM/MM案例教程二 二苯甲酮
Benzophenone结构准备
首先准备二苯甲酮Benzophenone的坐标文件,命名为BPH.xyz
24
C -2.54700 0.45510 0.06680
C -2.54160 -0.01810 1.38630
C -3.74290 -0.40660 1.99760
C -4.94170 -0.34290 1.28250
C -4.94480 0.12330 -0.03620
C -3.74920 0.52640 -0.64160
C -1.27680 -0.08120 2.18450
O -1.26930 0.16880 3.37250
C -0.02150 -0.46400 1.46430
C 1.18620 0.13430 1.85330
C 2.37660 -0.21530 1.21040
C 2.36490 -1.17300 0.19100
C 1.16310 -1.78220 -0.18680
C -0.03080 -1.42830 0.44700
H 1.18770 0.86620 2.66440
H -3.73280 -0.75010 3.03460
H 3.31310 0.25350 1.50860
H -5.87330 -0.64990 1.75530
H 3.29390 -1.44820 -0.30740
H -5.88040 0.17660 -0.59220
H 1.15790 -2.53420 -0.97410
H -3.75550 0.89780 -1.66500
H -0.96650 -1.90720 0.15500
H -1.61620 0.77400 -0.40440
使用Open Babe、Amber插件antechamber得到键、电荷等信息 命令行操作:
obabel BPH.xyz -O BPH_mid.mol2
#默认得到分子名为NUL,可将mol2中分子名字替换为BPH,此示例未进行该操作
antechamber -i BPH_mid.mol2 -fi mol2 -o BPH.mol2 -fo mol2 -c bcc -at gaff
使用Amber中parmchk工具得到力场参数 命令行操作:
parmchk -a Y -i BPH.mol2 -f mol2 -o BPH.frcmod
使用tleap进行溶剂化处理,并得到小分子的lib文件以及体系的溶剂化处理 准备文件tleap.in,溶剂化处理得到top、crd文件(文件名为BPH_solv.top BPH_solv.crd)
source leaprc.protein.ff14SB
source leaprc.water.tip3p
loadamberparams BPH.frcmod
BPH=loadmol2 BPH.mol2
check BPH
saveoff BPH BPH.lib
solvateoct BPH TIP3PBOX 18.0
saveamberparm BPH BPH_solv.top BPH_solv.crd
quit
命令行运行:
tleap -f tleap.in
得到初始构想的BPH_solv.top BPH_solv.crd.
动力学平衡
创建文件夹md/,在此文件夹中准备动力学模拟所需文件:最小化输入文件
01_Min.in,
升温输入文件
02_Heat.in,
平衡输入文件
03_Prod.in.
使用Amber中sander依次进行分子动力学最小化、升温、平衡;
命令行依次运行:
### Optimization
sander -O -i 01_Min.in -o 01_Min.out -p ../BPH_solv.top -c ../BPH_solv.crd -r 01_Min.rst -inf 01_Min.mdinfo
### Heat
sander -O -i 02_Heat.in -o 02_Heat.out -p ../BPH_solv.top -c 01_Min.rst -r 02_Heat.rst -x 02_Heat.mdcrd -inf 02_Heat.mdinfo
### Production
sander -O -i 03_Prod.in -o 03_Prod.out -p ../BPH_solv.top -c 02_Heat.rst -r 03_Prod.rst -x 03_Prod.mdcrd -inf 03_Prod.mdinfo
动力学结果分析
随机选取单帧结构，截取部分水的构象
1 使用cpptraj获取单帧构象(仅作示例，故随机选取一帧)
准备输入文件snap.trajin
parm ../BPH_solv.top
trajin 03_Prod.mdcrd 2976 2976 1 # read from mdcrd frames 2976 to 2976 (1 frame)
center :1 # put BPH in the center
image familiar # re-image
trajout snapshot_2976.rst rest # write the coordinates of this frame
go
命令行运行:
cpptraj -i snap.trajin
2 截取部分水的构象
删去与BPH中距离C7原子 > 20Å 的水分子,准备输入文件strip.trajin
parm ../BPH_solv.top
trajin snapshot_2976.rst # read the snapshot
reference snapshot_2976.rst rest # use it as reference (necessary for strip command)
strip @7>:20.0 # strip all waters further than 20A around atom C7
trajout strip_2976.pdb pdb # write pdb output
go
命令行运行
cpptraj -i strip.trajin
得到新的溶剂化体系
strip_2976.pdb.
QM/MM计算准备
1 top和crd文件准备
pDynamo使用Amber的top和crd文件作为输入,依据strip_2976.pdb文件,使用前面生成的力场文件得到该文件对应的Amber的top和crd文件。
新建并进入文件夹md/get_topcrd/,准备tleap的输入文件
tleap.in.
ource leaprc.protein.ff14SB
source leaprc.water.tip3p
loadamberparams ../../BPH.frcmod
loadoff ../../BPH.lib
a=loadpdb strip_2976.pdb
check a
saveamberparm a BPH_new.top BPH_new.crd
savepdb a BPH_new.pdb
quit
命令行运行,得到新的top和crd文件(
BPH_new.top、
BPH_new.crd、
BPH_new.pdb
)
2 活性区域水分子层选取 VMD 中选择距离二苯甲酮周围3Å的水作为可运动的水分子层,vmd中按照如下设置可显示二苯甲酮以及其周围3Å的水层
BPH和3Å内的水构象如下图所示
整体构象如下图所示
整个体系QM/MM区域划分如下图所示
使用VMD 中TkConsole得到MM区活性区域原子index,后续需用于QM/MM输入文件中; TkConsole控制台依次键入:
如图所示
QM/MM计算
基态优化
新建文件夹qmmm/,将BPH_new.crd,BPH_new.top文件拷贝至该目录
新建qmmm/ground_opt/文件夹,进行基态的QM/MM几何构型优化
准备QM/MM输入文件opt.py,其中定义QM区域和可活动的MM区域:
#. Define Atoms List
natoms = len ( molecule.atoms )
qm_list = range(24)
activate_list = [387,388,389,390,391,392,402,403,404,552,553,554,624,625,626,1104,1105,1106,
1203,1204,1205,1359,1360,1361,1419,1420,1421,1554,1555,1556,1572,1573,1574,
1611,1612,1613,1617,1618,1619,1845,1846,1847,1944,1945,1946,2139,2140,2141,
2262,2263,2264,2337,2338,2339,2460,2461,2462,2568,2569,2570,2736,2737,2738]
mm_list = range ( natoms )
for i in qm_list :
mm_list.remove( i )
mm_inactivate_list = mm_list[ : ]
for i in activate_list :
mm_inactivate_list.remove( i )
# . Define the selection for the first molecule.
qmmmtest_qc = Selection.FromIterable ( qm_list )
# . Define Fixed Atoms
selection_qm_mm_inactivate = Selection.FromIterable ( qm_list + mm_inactivate_list )
最小化100步
ConjugateGradientMinimize_SystemGeometry ( molecule ,
logFrequency = 2,
maximumIterations = 100,
rmsGradientTolerance = 0.1,
trajectories = [ ( trajectory, 2 ) ])
QM模型选择
qcModel = QCModelBDF ( "GB3LYP:6-31g" )
文件全文见
opt.py
QM/MM基态几何优化结果
S1态优化
新建qmmm/s1_opt/文件夹,进行s1态的QM/MM几何构型优化
准备QM/MM输入文件opt.py,其中定义QM区域和可活动的MM区域(同基态)
QM模型选择QCModelBDF_TDGRAD1类使用模板文件进行激发态的几何构型优化;
qcModel = qcModel = QCModelBDF_TDGRAD1 ( template = 'temple.inp', exgrad = 1 )
文件全文见
opt.py;
其中模板文件
opt.py
中激发态梯度设置为s1激发态的梯度。
QM/MM s1态几何优化结果
膜蛋白体系的QMMM计算
体系视图
QM/MM计算准备
1. QM/MM输入文件准备
体系的top文件和坐标文件
01_Min.in
01_Min.in
2. QM区域选取 将LEU 30设置为QM区,同时包含邻近的C=O基团 如图所示,淡黄色区域为QM区原子
其中LEU 30的index为
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
连接的C=O的index为
431 432
3. MM活性区域分子选取 MM活性区域分子选取 VMD 中选择距离LEU30周围8Å的分子作为MM活性区域
QM/MM区域划分,QM 区域包含LEU30和周围8Å原子
使用VMD 中TkConsole得到MM区活性区域原子index
TkConsole控制台依次键入
# 选择LEU30周围8Å的原子
set sel [atomselect 0 "same residue as exwithin 8 of residue 29"]
#获取LEU30周围8Å的原子的index
$sel get index
如同所示
QM/MM优化
结构优化
- 新建文件夹qmmmopt文件夹进行基态的QM/MM几何构型优化
- 准备QM/MM输入文件
opt.py
其中定义QM区域和可活动的MM区域:
#. Define Atoms List
natoms = len ( molecule.atoms )
qm_list = range (431, 452 )
activate_list = [118,119,120,121,122,123,124,125,126,127,128,129,130,131,
132,314,315,316,317,318,319,320,321,322,323,324,325,326,
327,328,329,330,331,332,333,334,335,336,337,338,339,340,
341,342,343,344,345,346,347,348,349,350,351,352,353,354,
355,356,357,358,359,360,361,362,363,364,365,366,367,368,
369,370,371,372,373,374,375,376,377,378,379,380,381,382,
383,384,385,386,387,388,389,390,391,392,393,394,395,396,
397,398,399,400,401,402,403,404,405,406,407,408,409,410,
411,412,413,414,415,416,417,418,419,420,421,422,423,424,
425,426,427,428,429,430, 452,453,454,455,456,457,
458,459,460,461,462,463,464,465,466,467,468,469,470,471,
472,473,474,475,476,477,478,479,480,481,482,483,484,485,
486,487,488,489,490,491,492,493,494,495,496,497,498,499,
500,501,502,503,504,505,506,507,508,509,510,511,512,513,
514,515,516,517,518,519,520,521,522,523,524,525,526,527,
528,529,530,531,532,533,534,535,536,537,538,539,540,541,
542,543,544,545,546,547,548,549,550,551,552,553,554,555,
1064,1065,1066,1067,1068,1069,1070,1071,1072,1073,1074,1075,
1076,1077,1078,1079,1080,1081,1082,1083,1084,1085,1086,1087,
1088,1089,1090,1091,1092,1093,1113,1114,1115,1116,1117,1118,
1119,1120,1121,1122,1123,1124,1125,1126,1127,1128,1129,1130,
1131,1132,1133,1134,1135,1136,1137,1138,1139,1140,1141,1142,
1143,1144,1145,1146,1147,1148,1149,1150,1151,1152,1153,1154,
1155,1156,1157,1158,1159,1160,1161,1162,1163,1164,1165,1166,
1167,1168,1169,1170,1171,1172,1173,1174,1175,1176,1177,1178,
1179,1180,1181,1182,1183,1184,1185,1186,1187,1188,1189,1190,
1191,1192,1193,1194,1195,1196,1197,1198,1199,1200,1201,1202,
1203,1204,1205,1206,1207,1208,1209,1210,1211,1212,1213,1214,
1215,1216,1217,1218,1219,1220,1221,1222,1223,1224,1225,1226,
1227,1228,1229,1230,1231,1232,1233,1234,1235,1236,1237,1238,
1239,1240,1241,1242,1243,1244,1245,1246,1247,1248,1249,1250,
1251,1252,1253,1254,1255,1256,1257,1258,1259,1260,1261,1262,
1263,1264,1265,1266,1267,1268,1269,1270,1271,1272,1273,1274,
1275,1276,1277,1300,1301,1302,1303,1304,1305,1306,1307,1308,
1309,1310,1311,1312,1313,1314,1315,1316,1317,1318,1319,1580,
1581,1582,1583,1584,1585,1586,1587,1588,1589,1590,1591,1592,
1593,1594,1595,1596,1597,1598,1599,1614,1615,1616,1617,1618,
1619,1620,1621,1622,1623,1624,1625,1626,1627,1635,1636,1637,
1638,1639,1640,1641,1642,1643,1644,1645,1646,1647,1648,1649,
1650,1651,1652,1653,1654,1655,1656,1657,1658,1659,1660,1661,
1662,1663,1664,1665,1666,1667,1668,1669,1670,1671,1672,1673,
1674,1675,1676,1677,1678,1679,3785,3786,3787,3788,3789,3790,
3791,3792,3793,3794,3795,3796,3797,3798,3799,3800,3801,3802,
3803,3804,3805,3806,3807,3808,3809,3810,3811,3812,3813,3814,
3815,3816,3817,3828,3829,3830,3831,3832,3833,3834,3835,3836,
3837,3838,3839,3840,3841,3842,3843,3844,3845,3846,3847,3848,
3849,3850,3851,3852,3853,3854,3855,3856,3857,3858,3859,3860,
3861,3862,3863,3864,3865,3866,3867,3868,3869,3870,3871,3872,
3873,3874,3875,3876,3877,3878,3879,3880,3881,3882,3883,3884,
3885,3886,3887,3904,3905,3906,3907,3908,3909,3910,3911,3912,
3913,3914,3915,3916,3917,3918,3919,3920,3921,3922,3923,3924,
3925,3926,3927,3928,3929,3930,3931,3932,3933,3934,3935,3936,
3937,3938,3939,3940,3941,3942,3943,3944,3945,3946,3947,3948,
3949,3950,3951,3952,3953,3954,3955,3956,3957,3968,3969,3970,
3971,3972,3973,3974,3975,3976,3977,3978,3979,3980,3981,3982,
3983,3984,3985,3986,3987,3988,3989,3990,3991,3992,3993,3994,
3995,5503,5504,5505,5506,5507,5508,5509,5510,5511,5512,5513,
5514,5515,5516,5517,5518,5519,5520,5521,5522,5523,5524,5525,
5526,5527,5528,5529,5530,5531,5532,5533,5534,5535,5536,5537,
5538,5539,5540,6044,6045,6046,6047,6048,6049,6050,6051,6052,
6053,6054,6055,6056,6057,6058,6059,6060,6061,6062,6063,6064,
6110,6111,6112,6113,6114,6115,6116,6117,6118,6119,6120,6121,
6122,6123,6124,6125,6126,6127,6128,6129,6130,6131,6132,6133,
6134,6135,6136,6137,6138,6139,6140,6141,6142,6818,6819,6820,
6821,6822,6823,6824,6825,6826,6827,6828,6829,6830,6831,6832,
6833,6834,6835,6836,6837,6838,50511,50512,50513,50538,50539,
50540,50544,50545,50546,50547,50548,50549,50559,50560,50561,
50586,50587,50588,50619,50620,50621,50652,50653,50654,50748,
50749,50750,51009,51010,51011,51030,51031,51032]
mm_list = range ( natoms )
for i in qm_list :
mm_list.remove( i )
mm_inactivate_list = mm_list[ : ]
for i in activate_list :
mm_inactivate_list.remove( i )
# . Define the selection for the first molecule.
qmmmtest_qc = Selection.FromIterable ( qm_list )
# . Define Fixed Atoms
selection_qm_mm_inactivate = Selection.FromIterable ( qm_list + mm_inactivate_list )
selection_mm = Selection.FromIterable ( mm_list )
selection_mm_inactivate = Selection.FromIterable ( mm_inactivate_list )
迭代优化MM区,QM区
# . Optimization.
# . Define the number of iterations
iterations = 2
# . QM region and MM region were fixed in turn for optimization
for i in range ( iterations ):
opt_LBFGSMinimize ( molecule, selection_qm_mm_inactivate )
opt_LBFGSMinimize ( molecule, selection_mm)
# . QM region and MM region were optimized simultaneously
opt_LBFGSMinimize ( molecule, selection_mm_inactivate)
QM模型选择 NB模型选择
# . Define the energy models.
nbModel = NBModelBDF ( )
nbModel.SetOptions (qcmmCoupling = 'Z1 Coupling')
qcModel = QCModelBDF ( "GB3LYP:6-31g" )
优化结构和原始结构构型
QM/MM边界选择算例教程
本教程示例QMMM边界的选择对几何构型优化的影响，错误的边界选择可能会导致意向不到的构象变化。
体系视图
本示例以五肽作为计算体系,以测试成建体系QM区的选择。体系
(5ala.pdb)
由五个ALA组成,N端、C端分别使用ACE、NME进行封端。
QM区域的选取
以五肽体系为例：
QM区域选择一:
选取第三个丙氨酸ALA4作为QM区,邻近的ALA作为其active MM区
qm_list = range (26, 36 )
activate_list = range ( 16, 26 ) + range ( 36, 46 )
QM/MM优化结果
QM区选择时,邻近MM区若带有较强的电荷(如本例的C=O),QM区和MM区之间静电相互作用较强, 体系无法收敛,如上动画所示,由于邻近MM区电荷较大,QM-MM之间的静电相互作用会导致QM区结构优化出现问题。
QM区域选择二:
选取第三个丙氨酸ALA4作为QM区,以及其邻近的C=O,其余邻近的ALA作为其active MM区
qm_list = range (24, 36 )
activate_list = range ( 16, 24 ) + range ( 36, 46 )
QM/MM优化结果
将C=O基团包含在QM区,体系优化收敛,体系结构在合理范围中进行优化。
蓝光HLCT分子的光物理特征研究
有机电致发光二极管（OLED）在显示和照明领域有广泛应用，其中发光染料是OLED器件中最关键的一环。当空穴和电子在发光层复合时，发光染料分子就会受激来到激发态（图1a）。根据激子统计规则，这种复合会产生25%的单重态激子和75%的三重态激子。根据选择定则，由于纯有机体系的三重态到单重态的跃迁通常是禁阻的，这75%的三重态激子会被浪费。通过合适的方法来利用三重态激子，是OLED发光染料分子设计的挑战。
马於光课题组提出的杂化局域-电荷转移激发态（HLCT），是一种可以利用三重态激子的发光机理。通过构筑具有HLCT特征的 \(\rm S_{1}\) 激发态，可以增加 \(\rm T_{n}\) - \(\rm S_{1}\) 的旋轨耦合（El-Sayed 规则），以促进 \(\rm T_{n}\) - \(\rm S_{1}\) 的反系间窜越（RISC）；除此之外，增加 \(\rm T_{n}\) - \(\rm T_{1}\) 的能级差并减小 \(\rm T_{n}\) - \(\rm S_{1}\) 的能级差，也有利于反系间窜越。把给体和芳香环相连是一种构筑HLCT分子的常用手段。
由于能级较高，蓝光HLCT分子的光稳定性通常较差。芘（Py）是一类具有高光稳定性的芳香环，而二苯基胺（TPA）是一种常见的强力给体片段。把这两者通过单键连接也许是一种构筑高稳定性蓝光HLCT分子的手段（图1b）。本文将通过量子化学计算来考察4,9-二甲基-1-N,1-N,6-N,6-N-四苯基芘-1,6-二胺（2TPA-Py）是否具有蓝光HLCT特征，验证donor-π-donor分子设计的可行性。
图 1. （a）HLCT分子的电致发光机理示意图。（b）2DPA-Py的分子结构。
基态 \(\rm S_{0}\) 优化
在对分子做精确计算前我们要得到可靠的基态结构（ \(\rm S_{0}\) ），也就是对基态做结构优化和振动分析。首先使用Gaussian对分子基态 \(\rm S_{0}\) 进行结构优化，采用泛函B3LYP、基组6-31g** ，并考虑色散矫正。选用B3LYP泛函的原因是它计算效率高、对积分格点精度的依赖性低，并加上DFT-D3(BJ)色散校正以描述可能的非共价键作用。另外，由于对于结构优化和振动分析，其计算结果对基组的敏感性较小，选用小基组可以节约时间。 gjf 输入文件如下：
%nprocshared=32
%chk=opt.chk
#P B3LYP/6-31g** em=GD3BJ opt freq
Title Card Required
0 1
C 1.23142086399057 -0.01598991339220 -0.21183114910141
C 2.72871444990389 -0.00550853390551 -0.13727964385924
C 3.41061821385248 -1.20117932087473 -0.03319288637410
C 4.79904649004982 -1.26160428059098 0.03464278981730
N 5.41364048996203 -2.51932450120475 0.15418467741142
C 4.92747565422558 -3.57250001859245 -0.62664098365207
C 4.55022212480360 -3.33536995952102 -1.94901522652513
C 4.05642159683765 -4.36315245643423 -2.72763628208729
C 3.93617240694700 -5.64193595032452 -2.20823626842045
C 4.30884346659409 -5.88288539461667 -0.89597870946301
C 4.79580537663213 -4.85982627352166 -0.10603578378833
C 6.36147678282021 -2.74991702364612 1.15337779307768
C 6.31970072442968 -2.03371232006357 2.35000679508950
C 7.27393251117782 -2.25099389609825 3.32426315813784
C 8.27590697003017 -3.18766118708190 3.13164605217437
C 8.32033728954936 -3.90430287118972 1.94668478315316
C 7.37766488670707 -3.68777915010114 0.96152123805281
C 5.55426311556430 -0.07485171406329 -0.01954737716259
C 4.87062379765459 1.16805677077445 -0.07171305294141
C 3.45350429641720 1.20058735786973 -0.14489187126324
C 2.81032456952191 2.46548965694624 -0.24263234170365
C 3.51112662376826 3.62218562784879 -0.20116625218499
C 4.92681570379474 3.62501895614023 -0.06751464021658
C 5.61185537735068 2.38194680680533 -0.05584863189100
C 7.03069025627284 2.34956069156028 -0.03333330974132
C 7.67833561978987 1.08357530517932 -0.06820637808522
C 6.97580495150633 -0.07284382674662 -0.05720475575518
C 7.75357292329654 3.55630546250584 0.00437058511537
C 7.06657662311645 4.75304774446334 0.04578195006735
C 5.67683694346796 4.81291513479720 0.02109801284253
N 5.05364278941578 6.07063412409978 0.08634601203580
C 4.04430276290905 6.30650817126718 1.02225220987505
C 4.01813497397634 5.60452422886364 2.22773071875455
C 3.00342478983393 5.82342539394774 3.13839692752355
C 2.00692660848240 6.74799998492145 2.87285054139819
C 2.02998237520816 7.45096162278722 1.67912787314289
C 3.03391132115340 7.23257419582758 0.75684192709414
C 5.58014732987927 7.11518216961340 -0.67889062460768
C 6.04635677519681 6.86001855649965 -1.96924843092115
C 6.57901675320538 7.88011253315299 -2.73218118328052
C 6.65054552833193 9.16913005476952 -2.22939545341485
C 6.19012117760983 9.42795662197789 -0.94876996878721
C 5.66385569959879 8.41283112869426 -0.17399231367271
C 9.25252702611257 3.56676709389692 0.02982018386394
H 0.85832992156275 -1.03632876280351 -0.16904659946900
H 0.88933685332284 0.43501238502809 -1.14286926646989
H 0.79789574594339 0.54473304731639 0.61552131117326
H 2.85548199101679 -2.12853502172196 0.00413662861935
H 4.64873169324299 -2.33946969900849 -2.35587693447180
H 3.76953281079583 -4.16566889381519 -3.75015021922513
H 3.55256748783067 -6.44401899690169 -2.82037387929308
H 4.20934558011874 -6.87446192152345 -0.47900026581524
H 5.06879641140089 -5.04938936946588 0.92181407561931
H 5.53591434202605 -1.30746574426274 2.50753626302272
H 7.22917194113817 -1.68928048139264 4.24582973230498
H 9.01732005764785 -3.35753750290448 3.89755908403090
H 9.10183011620465 -4.63207070344374 1.78327798868360
H 7.42600935609078 -4.23604224027954 0.03176328355712
H 1.73842332597976 2.49952836140925 -0.35594215471611
H 2.99326193531308 4.56447070267168 -0.28988821662774
H 8.75546098773894 1.04836100119705 -0.10774432599440
H 7.49820003451026 -1.01601354833581 -0.09649295442825
H 7.61798445128552 5.68140497373268 0.10650429470265
H 4.79705375266048 4.88765786929865 2.44212003924378
H 2.99594325952520 5.27247674972545 4.06748299539966
H 1.21761156723492 6.91895418257396 3.58905207574990
H 1.25350653787419 8.16906880804979 1.45894118248852
H 3.03809992945878 7.76959529262987 -0.18069135423844
H 5.98545198075639 5.85630503175884 -2.36409915773676
H 6.93473355133126 7.66838923944818 -3.72993019870071
H 7.06415340233785 9.96525107867144 -2.82966866197266
H 6.25071932363493 10.42779312892758 -0.54419248783976
H 5.32231798384009 8.61698483651345 0.83028525237055
H 9.62929685852355 3.02395636423240 0.89604599037109
H 9.62231936463369 4.58813097786218 0.07624892230056
H 9.65623043406863 3.09668205158134 -0.86638315237225
Gaussian以均方根力（Force-RMS）、最大力（Force-Max）、均方根步长（Step-RMS）、最大步长（Step-Max）四个标准来判断分子结构是否收敛。作业结束后，打开 log 输出文件，找到如下关键词
Item Value Threshold Converged?
Maximum Force 0.000012 0.000450 YES
RMS Force 0.000002 0.000300 YES
Maximum Displacement 0.001298 0.001800 YES
RMS Displacement 0.000336 0.001200 YES
Predicted change in Energy=-5.865194D-09
Optimization completed.
-- Stationary point found.
将达到收敛后的结构提取出来，作为初始结构用于后续计算。下图为部分截取：
Standard orientation:
---------------------------------------------------------------------
Center Atomic Atomic Coordinates (Angstroms)
Number Number Type X Y Z
---------------------------------------------------------------------
1 6 0 1.565649 -4.134284 -0.176953
2 6 0 1.643599 -2.628388 -0.152751
3 6 0 2.886000 -2.003080 -0.112539
4 6 0 3.020859 -0.613479 -0.088550
5 7 0 4.325420 -0.044576 -0.036499
6 6 0 5.278240 -0.466495 -0.995386
7 6 0 4.884618 -0.626003 -2.332720
当分子处于势能面的极小点时，一般不存在虚频（负频）。为了验证结构的可靠性还要检查频率计算结果。在 log 输出文件找到如下关键词，由于振动频率是从小到大排列的，观察前几个频率没有虚频说明分子处于势能面的局域极小点，其分子结构一般是可靠的。
1 2 3
A A A
Frequencies -- 12.2419 16.5051 20.1875
Red. masses -- 5.7419 5.7737 5.1174
Frc consts -- 0.0005 0.0009 0.0012
IR Inten -- 0.0197 0.0135 1.5193
Atom AN X Y Z X Y Z X Y Z
1 6 0.01 0.00 -0.04 -0.01 -0.00 -0.09 -0.02 0.03 -0.18
2 6 0.00 0.00 -0.04 -0.01 -0.01 -0.07 -0.02 0.03 -0.12
3 6 0.00 0.00 -0.03 -0.01 -0.01 -0.06 -0.02 0.03 -0.11
4 6 0.00 0.00 -0.03 -0.00 -0.01 -0.04 -0.01 0.03 -0.06
5 7 -0.00 0.01 -0.01 -0.01 -0.01 -0.00 -0.01 0.02 -0.05
吸收光谱
大多数有机分子的基态是闭壳层，通常是单重态，一般用 \(\rm S_{0}\) 表示。根据光化学第一定律（Stark-Einstein），分子需要吸收一个光子来使单电子从占据轨道跃迁到未占轨道，该光子的能量须与基态和激发态之间的能量差值一致。
那么只要一个分子受适当能量的光照射，就会吸收光子达到电子激发态么？
在量子力学中，从费米黄金规则（Fermi’s Golden Rule）可以推导出：分子从一种状态跃迁到另一种状态是需要满足多个条件的，这就是选择规则（selection rule）。具体的推导过程这儿不做赘述，建议读者直接查阅光物理的相关书籍，这里只简单列举适用于纯有机体系的部分选择规则：
满足Franck-Condon 原理：即电子跃迁过程中，起始和最终状态分子的几何结构（原子核位置）不变。
如果激发态的几何结构与基态的差别很大，那么把电子由基态的最低振动能级激发到激发态的最低振动能级时，分子的几何构型需发生变化。但是，由于原子核的质量远大于电子，且跟不上电子的运动速度，因此这种激发的概率很小。
电子自旋不变
在单电子近似下，自旋轨道正交归一。如果跃迁前后自旋不同，则自旋重叠积分必然为0，即跃迁禁阻。这就是自旋选择规则：“单重态→单重态、三重态→三重态允许；单重态→三重态、三重态→单重态禁阻”。
轨道重叠
跃迁前后分子轨道必须要有重叠，否则电子跃迁偶极矩积分为0，即跃迁禁阻。
轨道对称性改变
如果分子轨道具有对称性，除了轨道重叠，跃迁前后的轨道对称性必须不同。根据中心反演对称性把轨道分成g（中心对称）和u（中心反对称）两种，具体说法：“g→u，u→g允许，u→u，g→g禁阻”。
人们往往习惯以某占据轨道的电子向某虚轨道跃迁的形式来描述电子跃迁问题，例如自然跃迁轨道（NTO）分析：以一对占主导地位的轨道跃迁模式来阐述跃迁的本质。对于激发态电子结构分析需要做含时密度泛函（TDDFT）计算，本案例激发态选取M062x泛函，Def2SVP基组，对单重态和三重态分别计算8个态。选择M062x泛函是猜测分子的基态可能具有部分电荷转移（CT）特性。对于这一类的激发态，如果选择HF成分较低的泛函，
可能出现ghost态（不存在的态）。为了保险起见，选择高HF成分的M062x。当然，其他高HF成分的泛函，如CAM-B3LYP和ωB97XD也可以使用。 IOP(9/40=4) 关键词是为了输出更多轨道信息，以便MO变换成NTO后找到对电子激发贡献最大的一对NTO。 gjf 输入文件如下：
%nprocshared=32
%mem=6GB
%chk=tddft.chk
#P M062x/Def2SVP TD(nstates=8, 50-50) IOP(9/40=4)
Title Card Required
0 1
C 1.565649 -4.134284 -0.176953
C 1.643599 -2.628388 -0.152751
C 2.886000 -2.003080 -0.112539
C 3.020859 -0.613479 -0.088550
N 4.325420 -0.044576 -0.036499
C 5.278240 -0.466495 -0.995386
C 4.884618 -0.626003 -2.332720
C 5.799622 -1.064339 -3.285338
C 7.121820 -1.335486 -2.928150
C 7.515238 -1.172730 -1.598991
C 6.603320 -0.752040 -0.633991
C 4.698175 0.744156 1.076980
C 4.091832 0.542012 2.325047
C 4.427797 1.352319 3.406784
C 5.379893 2.362883 3.271786
C 5.988191 2.561344 2.030562
C 5.647767 1.769732 0.937881
C 1.876752 0.211899 -0.120738
C 0.588755 -0.407573 -0.126231
C 0.473733 -1.832327 -0.148612
C -0.837555 -2.410247 -0.178205
C -1.959943 -1.640762 -0.161854
C -1.876694 -0.212016 -0.120592
C -0.588696 0.407452 -0.126517
C -0.473668 1.832190 -0.149905
C 0.837623 2.410088 -0.179888
C 1.960009 1.640614 -0.162984
C -1.643530 2.628249 -0.154622
C -2.885933 2.002973 -0.113981
C -3.020804 0.613392 -0.088992
N -4.325380 0.044559 -0.036541
C -4.698230 -0.743314 1.077511
C -4.091928 -0.540291 2.325457
C -4.428022 -1.349755 3.407783
C -5.380213 -2.360326 3.273499
C -5.988485 -2.559645 2.032401
C -5.647934 -1.768883 0.939143
C -5.278162 0.465750 -0.995788
C -4.884527 0.624147 -2.333249
C -5.799513 1.061738 -3.286228
C -7.121702 1.333222 -2.929264
C -7.515137 1.171550 -1.599978
C -6.603236 0.751616 -0.634632
C -1.565577 4.134127 -0.179894
H 2.565826 -4.573123 -0.163055
H 1.054624 -4.497512 -1.075888
H 1.015461 -4.526102 0.685940
H 3.789152 -2.604760 -0.095814
H 3.860188 -0.404796 -2.611283
H 5.478783 -1.182339 -4.316028
H 7.834889 -1.669709 -3.674725
H 8.537506 -1.389762 -1.303667
H 6.910123 -0.642734 0.399820
H 3.355698 -0.245579 2.436971
H 3.946094 1.182749 4.365139
H 5.642450 2.988231 4.118776
H 6.724707 3.349483 1.905053
H 6.110094 1.939850 -0.027830
H -0.934317 -3.489037 -0.217496
H -2.937842 -2.104213 -0.185630
H 0.934386 3.488850 -0.219935
H 2.937915 2.104041 -0.187063
H -3.789080 2.604673 -0.097685
H -3.355727 0.247315 2.436831
H -3.946346 -1.179515 4.366033
H -5.642871 -2.985014 4.120945
H -6.725084 -3.347795 1.907449
H -6.110250 -1.939666 -0.026456
H -3.860104 0.402673 -2.611627
H -5.478667 1.178888 -4.317013
H -7.834757 1.666861 -3.676114
H -8.537402 1.388848 -1.304838
H -6.910046 0.643140 0.399265
H -1.015334 4.526556 0.682686
H -2.565752 4.572979 -0.166245
H -1.054607 4.496714 -1.079120
作业完成后根据激发能绘制低激发态能级图。可以看到 \(\rm S_{1}\) 态和 \(\rm T_{2}\) 、 \(\rm T_{3}\) 态之间的能级差较小，若旋轨耦合矩阵元较大则存在发生系间窜越和反系间窜越的可能。
衡量两个态之间跃迁强度可以用振子强度来衡量，它是一个无量纲量。原子单位下|i> → |j>跃迁的振子强度表达式为
其中
其中 \(\rm E_{j}\) 和 \(\rm E_{i}\) 分别为两个态的能量。基态与某个激发态之间的振子强度越大，就越容易吸收相应频率的电磁波而跃迁到那个激发态上，那么在吸收光谱中相应的吸收峰也越强。一般情况下，振子强度小于0.001可认为是跃迁禁阻。
低激发态激发能、振子强度、跃迁偶极矩如表中所示。
激发态 |
激发能/eV |
振子强度 |
跃迁偶极矩/Debye |
|---|---|---|---|
S1 |
3.1509 |
0.6012 |
19.7948 |
T1 |
2.1539 |
0.0000 |
0.0000 |
T2 |
3.2507 |
0.0000 |
0.0000 |
绘制的吸收光谱如下。
将 chk 文件转换为 fchk 文件。用Multiwfn+VMD渲染NTO轨道。
\(\rm S_{0}\) → \(\rm S_{1}\) 跃迁贡献最大的NTO对贡献值为96.40%。
\(\rm S_{0}\) → \(\rm T_{1}\) 跃迁贡献最大的NTO对贡献值为95.52%。
\(\rm S_{0}\) → \(\rm T_{2}\) 跃迁贡献最大的NTO对贡献值为86.41%。
\(\rm S_{0}\) → \(\rm T_{3}\) 跃迁贡献最大的NTO对贡献值为62.93%。
从图中可以观察到, \(\rm T_{1}\) 和 \(\rm T_{3}\) 态是典型的局域激发（LE），而 \(\rm S_{1}\) 和 \(\rm T_{2}\) 态既有电荷转移，也有局域激发的成分，属于HLCT态。
激发态 \(\rm S_{1}\) 优化
荧光是冷光现象，一般指发生在自旋单重态之间的辐射过程。根据Kasha规则，它是最低激发态到基态的发射，一般为 \(\rm S_{1}\) 态到 \(\rm S_{0}\) 态。为了模拟荧光过程，还需要对激发态 \(\rm S_{1}\) 做结构优化和频率计算，得到 log 文件和 fchk 文件，为后续MOMAP的计算做准备。泛函和基组分别为M062x和Def2SVP， gjf 文件如下：
%nprocshared=32
%mem=6GB
%chk=s1opt.chk
#P opt freq M062x/Def2SVP TD(nstates=3,root=1)
Title Card Required
0 1
C 1.565649 -4.134284 -0.176953
C 1.643599 -2.628388 -0.152751
C 2.886000 -2.003080 -0.112539
C 3.020859 -0.613479 -0.088550
N 4.325420 -0.044576 -0.036499
C 5.278240 -0.466495 -0.995386
C 4.884618 -0.626003 -2.332720
C 5.799622 -1.064339 -3.285338
C 7.121820 -1.335486 -2.928150
C 7.515238 -1.172730 -1.598991
C 6.603320 -0.752040 -0.633991
C 4.698175 0.744156 1.076980
C 4.091832 0.542012 2.325047
C 4.427797 1.352319 3.406784
C 5.379893 2.362883 3.271786
C 5.988191 2.561344 2.030562
C 5.647767 1.769732 0.937881
C 1.876752 0.211899 -0.120738
C 0.588755 -0.407573 -0.126231
C 0.473733 -1.832327 -0.148612
C -0.837555 -2.410247 -0.178205
C -1.959943 -1.640762 -0.161854
C -1.876694 -0.212016 -0.120592
C -0.588696 0.407452 -0.126517
C -0.473668 1.832190 -0.149905
C 0.837623 2.410088 -0.179888
C 1.960009 1.640614 -0.162984
C -1.643530 2.628249 -0.154622
C -2.885933 2.002973 -0.113981
C -3.020804 0.613392 -0.088992
N -4.325380 0.044559 -0.036541
C -4.698230 -0.743314 1.077511
C -4.091928 -0.540291 2.325457
C -4.428022 -1.349755 3.407783
C -5.380213 -2.360326 3.273499
C -5.988485 -2.559645 2.032401
C -5.647934 -1.768883 0.939143
C -5.278162 0.465750 -0.995788
C -4.884527 0.624147 -2.333249
C -5.799513 1.061738 -3.286228
C -7.121702 1.333222 -2.929264
C -7.515137 1.171550 -1.599978
C -6.603236 0.751616 -0.634632
C -1.565577 4.134127 -0.179894
H 2.565826 -4.573123 -0.163055
H 1.054624 -4.497512 -1.075888
H 1.015461 -4.526102 0.685940
H 3.789152 -2.604760 -0.095814
H 3.860188 -0.404796 -2.611283
H 5.478783 -1.182339 -4.316028
H 7.834889 -1.669709 -3.674725
H 8.537506 -1.389762 -1.303667
H 6.910123 -0.642734 0.399820
H 3.355698 -0.245579 2.436971
H 3.946094 1.182749 4.365139
H 5.642450 2.988231 4.118776
H 6.724707 3.349483 1.905053
H 6.110094 1.939850 -0.027830
H -0.934317 -3.489037 -0.217496
H -2.937842 -2.104213 -0.185630
H 0.934386 3.488850 -0.219935
H 2.937915 2.104041 -0.187063
H -3.789080 2.604673 -0.097685
H -3.355727 0.247315 2.436831
H -3.946346 -1.179515 4.366033
H -5.642871 -2.985014 4.120945
H -6.725084 -3.347795 1.907449
H -6.110250 -1.939666 -0.026456
H -3.860104 0.402673 -2.611627
H -5.478667 1.178888 -4.317013
H -7.834757 1.666861 -3.676114
H -8.537402 1.388848 -1.304838
H -6.910046 0.643140 0.399265
H -1.015334 4.526556 0.682686
H -2.565752 4.572979 -0.166245
H -1.054607 4.496714 -1.079120
作业完成后在 log 文件中找到最后一个Excited State 1为 \(\rm S_{1}\) 激发能，Total Energy为电子态能量。
Excited State 1: Singlet-A 2.7938 eV 443.79 nm f=0.8006 <S**2>=0.000
149 ->150 0.69410
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-DFT) = -1727.22867894
Copying the excited state density for this state as the 1-particle RhoCI density.
激发态 \(\rm T_{2}\) 、 \(\rm T_{3}\) 优化 ——————————————-——————————————————————————————————————————————————————————————
由于本案例后续会用MOMAP做 \(\rm T_{2}\) → \(\rm S_{1}\) 态和 \(\rm T_{3}\) → \(\rm S_{1}\) 态的反系间窜越速率计算，所以前期还需要量化软件做激发态 \(\rm T_{2}\) 和 \(\rm T_{3}\) 的结构优化和频率计算，得到 log 文件和 fchk 文件。泛函和基组分别为M062x和Def2SVP， T2.gjf 文件如下：
%nprocshared=32
%mem=6GB
%chk=t2.chk
#P opt freq M062x/Def2SVP TD(triplets,nstates=6,root=2)
Title Card Required
0 1
C 1.565649 -4.134284 -0.176953
C 1.643599 -2.628388 -0.152751
C 2.886000 -2.003080 -0.112539
C 3.020859 -0.613479 -0.088550
N 4.325420 -0.044576 -0.036499
C 5.278240 -0.466495 -0.995386
C 4.884618 -0.626003 -2.332720
C 5.799622 -1.064339 -3.285338
C 7.121820 -1.335486 -2.928150
C 7.515238 -1.172730 -1.598991
C 6.603320 -0.752040 -0.633991
C 4.698175 0.744156 1.076980
C 4.091832 0.542012 2.325047
C 4.427797 1.352319 3.406784
C 5.379893 2.362883 3.271786
C 5.988191 2.561344 2.030562
C 5.647767 1.769732 0.937881
C 1.876752 0.211899 -0.120738
C 0.588755 -0.407573 -0.126231
C 0.473733 -1.832327 -0.148612
C -0.837555 -2.410247 -0.178205
C -1.959943 -1.640762 -0.161854
C -1.876694 -0.212016 -0.120592
C -0.588696 0.407452 -0.126517
C -0.473668 1.832190 -0.149905
C 0.837623 2.410088 -0.179888
C 1.960009 1.640614 -0.162984
C -1.643530 2.628249 -0.154622
C -2.885933 2.002973 -0.113981
C -3.020804 0.613392 -0.088992
N -4.325380 0.044559 -0.036541
C -4.698230 -0.743314 1.077511
C -4.091928 -0.540291 2.325457
C -4.428022 -1.349755 3.407783
C -5.380213 -2.360326 3.273499
C -5.988485 -2.559645 2.032401
C -5.647934 -1.768883 0.939143
C -5.278162 0.465750 -0.995788
C -4.884527 0.624147 -2.333249
C -5.799513 1.061738 -3.286228
C -7.121702 1.333222 -2.929264
C -7.515137 1.171550 -1.599978
C -6.603236 0.751616 -0.634632
C -1.565577 4.134127 -0.179894
H 2.565826 -4.573123 -0.163055
H 1.054624 -4.497512 -1.075888
H 1.015461 -4.526102 0.685940
H 3.789152 -2.604760 -0.095814
H 3.860188 -0.404796 -2.611283
H 5.478783 -1.182339 -4.316028
H 7.834889 -1.669709 -3.674725
H 8.537506 -1.389762 -1.303667
H 6.910123 -0.642734 0.399820
H 3.355698 -0.245579 2.436971
H 3.946094 1.182749 4.365139
H 5.642450 2.988231 4.118776
H 6.724707 3.349483 1.905053
H 6.110094 1.939850 -0.027830
H -0.934317 -3.489037 -0.217496
H -2.937842 -2.104213 -0.185630
H 0.934386 3.488850 -0.219935
H 2.937915 2.104041 -0.187063
H -3.789080 2.604673 -0.097685
H -3.355727 0.247315 2.436831
H -3.946346 -1.179515 4.366033
H -5.642871 -2.985014 4.120945
H -6.725084 -3.347795 1.907449
H -6.110250 -1.939666 -0.026456
H -3.860104 0.402673 -2.611627
H -5.478667 1.178888 -4.317013
H -7.834757 1.666861 -3.676114
H -8.537402 1.388848 -1.304838
H -6.910046 0.643140 0.399265
H -1.015334 4.526556 0.682686
H -2.565752 4.572979 -0.166245
H -1.054607 4.496714 -1.079120
作业完成后在 log 文件中找到最后一个Excited State 2为 \(\rm T_{2}\) 激发能，Total Energy为电子态能量。
Excited State 2: Triplet-A 3.0388 eV 408.01 nm f=0.0000 <S**2>=2.000
138 ->150 -0.14038
148 ->150 0.61959
149 ->155 0.16846
149 ->157 -0.14448
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-DFT) = -1727.22151297
Copying the excited state density for this state as the 1-particle RhoCI density.
同样对 \(\rm T_{3}\) 优化可得 \(\rm T_{3}\) 态能量。结果显示 \(\rm T_{3}\) 电子态能量小于 \(\rm T_{2}\) 态能量，这表示 \(\rm T_{2}\) 与 \(\rm T_{3}\) 态在远离Frank-Condon区优化的过程中可能存在势能面交叉，使得最终优化的 \(\rm T_{3}\) 极小点能量小于 \(\rm T_{2}\) 。
Excited State 3: Triplet-A 2.9283 eV 423.40 nm f=0.0000 <S**2>=2.000
149 ->151 0.67322
149 ->155 0.11403
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-DFT) = -1727.22240086
Copying the excited state density for this state as the 1-particle RhoCI density.
自旋轨道耦合
自旋轨道耦合（SOC）体现的是电子自旋与电子绕核转动的相互作用。计算单重态和三重态之间的跃迁时，如果未考虑自旋轨道耦合（也就是二者耦合严格为0），那么他们的跃迁就是禁阻的；而把旋轨耦合引入哈密顿量之后，二者耦合不严格为0，此时单重态与三重态之间的跃迁就有可能发生。我们往往关心在某个特定结构下 \(\rm S_{i}\) 态与 \(\rm T_{j}\) 态之间的自旋轨道耦合。其中 \(\rm <S_{i}|SOC|T_{j}>\) 表示自旋轨道耦合矩阵元，取它的模来衡量 \(\rm S_{i}\) 、 \(\rm T_{j}\) 电子态之间的旋轨耦合作用大小。这一物理量还可以用于计算系间窜越（ISC）速率和反系间窜越（RISC）速率。
本例用BDF计算 \(\rm S_{1}\) - \(\rm T_{2}\) 和 \(\rm S_{1}-T_{3}\) 态之间的旋轨耦合矩阵元，采用M062x泛函，Def2SVP基组， inp 文件如下：
$compass
Title
C42H32N2
Geometry
C 1.586003 -4.127364 -0.277679
C 1.687147 -2.632259 -0.222611
C 2.912299 -2.005731 -0.177241
C 3.045692 -0.604559 -0.100709
N 4.333135 -0.047256 -0.013973
C 5.334680 -0.519389 -0.889609
C 5.015755 -0.759732 -2.235838
C 5.984347 -1.247032 -3.105784
C 7.281444 -1.498204 -2.653995
C 7.598599 -1.262567 -1.315545
C 6.635646 -0.783602 -0.432749
C 4.653716 0.850763 1.024080
C 3.972142 0.778096 2.248981
C 4.252686 1.694201 3.257436
C 5.214213 2.687133 3.070218
C 5.894013 2.759318 1.852004
C 5.615474 1.857669 0.831601
C 1.878302 0.232663 -0.125495
C 0.592246 -0.400870 -0.142073
C 0.488519 -1.828638 -0.190107
C -0.793696 -2.412608 -0.223962
C -1.944019 -1.641486 -0.186177
C -1.878300 -0.232665 -0.125495
C -0.592244 0.400868 -0.142073
C -0.488517 1.828637 -0.190108
C 0.793698 2.412606 -0.223963
C 1.944022 1.641485 -0.186177
C -1.687145 2.632258 -0.222613
C -2.912297 2.005730 -0.177243
C -3.045690 0.604558 -0.100710
N -4.333133 0.047256 -0.013974
C -4.653717 -0.850761 1.024079
C -3.972142 -0.778097 2.248980
C -4.252689 -1.694201 3.257435
C -5.214220 -2.687128 3.070218
C -5.894022 -2.759311 1.852005
C -5.615479 -1.857663 0.831601
C -5.334678 0.519389 -0.889611
C -5.015753 0.759734 -2.235839
C -5.984345 1.247034 -3.105785
C -7.281443 1.498203 -2.653996
C -7.598599 1.262562 -1.315546
C -6.635645 0.783598 -0.432750
C -1.586001 4.127363 -0.277682
H 2.581935 -4.588013 -0.284388
H 1.048934 -4.459802 -1.181513
H 1.027774 -4.525258 0.585788
H 3.824462 -2.607245 -0.175768
H 4.002023 -0.553508 -2.582465
H 5.725422 -1.423160 -4.150923
H 8.039250 -1.878188 -3.339953
H 8.605471 -1.467565 -0.948393
H 6.879710 -0.616710 0.617096
H 3.223855 -0.001794 2.394755
H 3.715870 1.625084 4.204829
H 5.431113 3.401551 3.864968
H 6.641732 3.536971 1.688418
H 6.130708 1.931595 -0.127046
H -0.888701 -3.496890 -0.277700
H -2.914627 -2.134719 -0.218301
H 0.888703 3.496889 -0.277701
H 2.914629 2.134717 -0.218301
H -3.824459 2.607244 -0.175771
H -3.223853 0.001791 2.394754
H -3.715872 -1.625085 4.204828
H -5.431123 -3.401546 3.864969
H -6.641744 -3.536960 1.688420
H -6.130714 -1.931587 -0.127046
H -4.002021 0.553512 -2.582466
H -5.725421 1.423163 -4.150924
H -8.039249 1.878187 -3.339954
H -8.605471 1.467558 -0.948395
H -6.879709 0.616703 0.617094
H -1.027772 4.525257 0.585785
H -2.581933 4.588012 -0.284392
H -1.048931 4.459800 -1.181516
End Geometry
Basis
Def2-SVP
Skeleton
Group
C(1)
$end
$xuanyuan
Heff
21
Hsoc
2
Direct
RS
0.33
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
M062X
MPEC+COSX
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
3
MPEC+COSX
Istore
1
$end
$tddft
Imethod
1
Isf
1
Idiag
1
Iroot
3
MPEC+COSX
Istore
2
$end
$tddft
Isoc
2
Nfiles
2
Ifgs
1
Imatsoc
2
1 1 1 2 1 2
1 1 1 2 1 3
$end
作业完成后，在 out 输出文件中找到如下关键词即为旋轨耦合矩阵元结果
[tddft_soc_matsoc]
Print selected matrix elements of [Hsoc]
SocPairNo. = 1 SOCmat = < 1 1 1 |Hso| 2 1 2 > Dim = 1 3
mi/mj ReHso(au) cm^-1 ImHso(au) cm^-1
0.0 -1.0 -0.0000031365 -0.6883766845 0.0000019744 0.4333410617
0.0 0.0 0.0000000000 0.0000000000 -0.0000000001 -0.0000289692
0.0 1.0 -0.0000031365 -0.6883766845 -0.0000019744 -0.4333410617
SocPairNo. = 2 SOCmat = < 1 1 1 |Hso| 2 1 3 > Dim = 1 3
mi/mj ReHso(au) cm^-1 ImHso(au) cm^-1
0.0 -1.0 -0.0000000002 -0.0000481328 0.0000000002 0.0000411125
0.0 0.0 0.0000000000 0.0000000000 -0.0000069588 -1.5272872617
0.0 1.0 -0.0000000002 -0.0000481328 -0.0000000002 -0.0000411125
这里的 \(\rm SOCmat=<1 1 1 |H_{SO}| 2 1 2>\) 表示矩阵元 \(\rm <S_{1}| H_{SO} |T_{2}>\) ，ReHso和ImHso分别表示实部和虚部，单位为au或 \(\rm cm^{-1}\) 。将三个mj分量的SOC矩阵元的模平方求和后再开平方得到后续MOMAP会用到的 \(\rm S_{1}\) 态与 \(\rm T_{2}\) 态旋轨耦合矩阵元，即1.15035 \(\rm cm^{-1}\) ； \(\rm S_{1}\) 态与 \(\rm T_{3}\) 态旋轨耦合矩阵元1.52729 \(\rm cm^{-1}\) 。
重整能
重整能是指当分子得失电子后，因几何结构的弛豫而导致的体系能量的变化。它既是影响电子转移速率（基于Marcus理论）的关键的物理量，也是影响光谱和辐射速率的重要因素。具体来讲，分子在跃迁初态、末态平衡位置时的能量差，就分别是基态和激发态的重整能 \(\lambda_{S0}=E_{3}-E_{1}\) ， \(\lambda_{S1}=E_{2}-E_{4}\) 。
重整能也可定义为
其中 \(\rm S_k\) 和 \(\rm ω_k\) 分别为第k个模式的黄昆因子(Huang−Rhys)和频率，D为模式位移。
黄昆因子
用MOMAP做电子振动耦合即evc计算可以得到重整能、黄昆因子等数据，输入文件需要 \(\rm S_0\) 和 \(\rm S_1\) 的结构优化和频率计算log文件和fchk文件，以及 momap.inp 文件， momap.inp 内容如下：
do_evc = 1
&evc
ffreq(1) = "s0.log"
ffreq(2) = "s1.log"
set_cart = t
/
作业完成后产生 evc.cart.dat 文件，找到下方关键词为 \(\rm S_0\) 和 \(\rm S_1\) 重整能。如下图， \(\lambda_{S0} = 1610.605 cm^{−1}\) ， \(\lambda_{S1} = 1864.085 cm^{−1}\) ，即基态和激发态重整能相差不大，说明两个态构型相差不大，属于同一个Franck-Condon区。
Total reorganization energy (cm-1): 1610.605075 1864.085048
在Device Studio中打开 evc.cart.dat 文件得到 \(\rm S_0\) 态和 \(\rm S_1\) 态重整能和黄昆因子在每个振动模式下的贡献。
分析振动模式，发现 \(\rm S_{0}\) 态的重组能主要贡献来自于1676.69 \(\rm cm^{-1}\) 的高频C-C伸缩振动和1308.32 \(\rm cm^{-1}\) 的高频弯曲振动， \(\rm S_1\) 态的重组能主要贡献来自于1683.31 \(\rm cm^{-1}\) 、1695.91 \(\rm cm^{-1}\) 和1414.86 \(\rm cm^{-1}\) 的高频弯曲振动。低频振动模式的黄昆因子显著，其中 \(\rm S_{0}\) 态的黄昆因子的最大模式为12.24 \(\rm cm^{-1}\) 的低频弯曲振动，S1态的黄昆因子最大模式为18.30 \(\rm cm^{-1}\) 的低频弯曲振动。
在分子的光物理过程中，由跃迁初态和末态的正则模式之间的相互交叠而引起的Duschinsky转动效应也会对光谱和速率产生重要影响，S0和S1极小点下分别得到的3N-6个正则坐标是不同的，它们彼此间是线性变换关系可表达为Q’’=J*Q’+ΔQ，这里Q’和Q’’分别代表两个电子态极小点下的正则模式，J称为Duschinsky矩阵。在Device Studio中打开 evc.cart.abs 文件就可以得到S0和S1态之间的Duschinsky转动矩阵的二维图。
荧光光谱
用MOMAP进行荧光辐射速率的计算需要上一步的结果 evc.cart.dat 以及新的输入文件 momap.inp ， momap.inp 内容如下：
do_spec_tvcf_ft = 1
do_spec_tvcf_spec = 1
&spec_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.09626 au
EDMA = 8.18309 debye
EDME = 9.64296 debye
FreqScale = 1.0
DSFile = "evc.cart.dat"
Emax = 0.3 au
dE = 0.00001 au
logFile = "spec.tvcf.log"
FtFile = "spec.tvcf.ft.dat"
FoFile = "spec.tvcf.fo.dat"
FoSFile = "spec.tvcf.spec.dat"
/
对于绝热激发能Ead，由于用Gaussian计算 \(\rm S_0\) 和 \(\rm S_1\) 时采用了不同的计算级别，因此我们用 \(\rm S_1\) 的结构在 \(\rm S_0\) 的计算级别下，再做一次单点计算用于修正：用此能量- \(\rm S_0\) 能量+ \(\rm S_1\) 激发能得到绝热激发能，即Ead=0.09626 au。对于吸收跃迁偶极矩EDMA，从S1.log文件读取第一个基态到激发态跃迁电偶极矩Dip.S.，开根号并单位换算得到8.18309 debye。对于发射跃迁偶极矩EDME，从S1.log文件读取最后一个基态到激发态跃迁电偶极矩Dip.S.，开根号并换单位换算得到9.64296 debye。（Ead、EDMA、EDME结果均在加了苯溶剂条件下获得(scrf(solvent=benzene,SMD))，模拟薄膜环境）。
作业完成后在 spec.tvcf.log 文件的末端能够读取辐射速率，本例 \(\rm S_1\) → \(\rm S_0\) 辐射速率为 \(\rm 1.77 \times 10^8 s^{-1}\) ，荧光寿命5.64 ns。
I^-1 = 3.05463233E+00 Hartree = 6.70414303E+05 cm-1 = 8.31208105E+01 eV
radiative rate (0): 4.28614462E-09 1.77195105E+08 /s, 5.64 ns
在Device Studio中打开 spec.tvcf.spec.dat 文件得到吸收和发射谱，如下图，吸收光在388 nm，发射光在497 nm。
反系间窜越速率
系间窜越是光化学中一种重要的无辐射过程。它是指分子受激后，由于不同自旋多重度的态的势能面之间存在交叉，导致体系经历这样的结构时以非辐射方式改变自旋多重度。在一般有机体系中，系间窜越（RISC）指的是从单重跃迁到三重态，反系间窜越（RISC）指的是从三重跃迁到单重态。反系间窜越速率，例如 \(\rm S_0 → t_2\) ，还与他们的能级差 \(\Delta E_{ST}\) 有关。这里 \(\Delta E_{ST}\) 可以用 \(\rm S_{1}\) 的激发态能量与 \(\rm T_{2}\) 激发态能量相减得到。算得 \(\rm S_1 - T_2\) 态 \(\Delta E_{ST}\) =0.05518 au， \(\rm S_1 - T_3\) 态 \(\Delta E_{ST}\) =0.05528 au。
在MOMAP程序中要想计算反系间窜越速率，首先要计算 \(\rm S_1\) 和 \(\rm T_2\) 的电子振动耦合，文件需要 \(\rm S_1\) 和 \(\rm T_2\) 的log频率计算文件以及fchk文件，还有 momap.inp 输入文件， momap.inp 内容如下：
do_evc = 1
&evc
ffreq(1) = "s1.log"
ffreq(2) = "t2.log"
set_cart = t
/
作业结束后生成的 evc.cart.dat 文件再与新的 momap.inp 文件放在同一目录下计算非辐射速率。此时输入文件 momap.inp 如下：
do_isc_tvcf_ft = 1
do_isc_tvcf_spec = 1
&isc_tvcf
DUSHIN = .t.
Temp = 298 K
tmax = 1500 fs
dt = 1 fs
Ead = 0.05518 au
Hso = 1.15035 cm-1
DSFile = "evc.cart.dat"
Emax = 0.3 au
logFile = "isc.tvcf.log"
FtFile = "isc.tvcf.ft.dat"
FoFile = "isc.tvcf.fo.dat"
/
Ead为 \(\Delta E_{ST}\) ， \(\rm H_{SO}\) 为 \(\rm S_1\) 态与 \(\rm T_2\) 态旋轨耦合矩阵元，计算得到的isc.tvcf.log文件末端为系间窜越速率和反系间窜越速率，本例 \(\rm k_{ISC} = 4.53 \times 10^4 s^{-1}\) ， \(k_{RISC} = 1.48 \times 10^2 s^{-1}`\) 。
# Intersystem crossing Ead is 0.0551800 au, rate is 4.53103856E+04 s-1, lifetime is 2.20699954E-05 s
# Reverse Intersystem crossing Ead is -0.0551800 au, rate is 1.47691362E+02 s-1, lifetime is 6.77087667E-03 s
同样地， \(\rm S_1\) 态与 \(T_3\) 态系间窜越速率 \(\rm k_{ISC} = 8.75 \times 10^7 s^{-1}\) ，反系间窜越速率 \(\rm k_{ISC} = 1.32 \times 10^7 s^{-1}\) 。
# Intersystem crossing Ead is 0.0552800 au, rate is 8.75255907E+07 s-1, lifetime is 1.14252300E-08 s
# Reverse Intersystem crossing Ead is -0.0552800 au, rate is 1.31729899E+07 s-1, lifetime is 7.59129104E-08 s
通过计算，我们发现 \(\rm T_2\) 态到 \(\rm S_1\) 态之间的反系间窜越速率很小，不满足HLCT分子的要求。 \(\rm T_3\) 态到 \(\rm S_1\) 态的反系间窜越速率很大，表明三重态激子有可能在 \(\rm T_3\) 态发生反系间窜越并转化为 \(\rm S_1\) 态。
小结
本文基于DFT和TDDFT理论，计算了2TPA-Py分子的激发态光物理过程。结果表明，2TPA-Py分子的 \(\rm S_1\) 态具有HLCT特征，其最大发射波长在497 nm为天蓝光。这个分子的 \(\rm T_3 → S_1\) 的反系间窜越速率高达 \(4.39 \times 106 s^{-1}\) ，基本满足了通过反系间窜越来利用三重态激子的要求。可见，donor-π-donor的分子设计策略有望成为构筑高稳定性蓝光HLCT分子的有效手段。
\(\rm Ir(ppy)_3\) 磷光发射机制的理论探究
自从磷光材料被应用到OLED以来，基于磷光机制的有机发光器件研究快速发展，其中代表性的有红、绿、蓝单基色磷光器件和全磷光白色磷光器件。由于自旋禁阻的原因，磷光量子产率要比荧光低得多，因此通常采用重原子效应等方法来提高磷光量子产率。当引入重原子后，自旋轨道耦合作用加强，使禁阻的三重态向单重态的跃迁变为允许，分子停滞在三线态的时间大幅缩短，极大地提高了器件的内量子效率，常用的重金属原子有Ir、Pt、Re、Os、Cu等。 磷光材料一般应具备良好的光热稳定性、较大的分子吸光截面、较高的系间窜越能力、室温下有较高的磷光量子产率和较短的三线态寿命等特征。由于有机材料中磷光的绿光材料最容易获得，因而研究最早也是从磷光绿光OLED开始，其中 \(\rm Ir(ppy)_3\) 是最具代表性的材料，被广泛地应用到磷光器件中。磷光发射速率是发光机理研究中的重要参数，本文将以 \(\rm Ir(ppy)_3\) 为例，使用BDF和MOMAP软件计算其磷光发射速率。首先需要 通过BDF量化软件进行结构优化、频率计算和旋轨耦合计算，之后基于BDF的结构优化和频率计算结果文件、旋轨耦合计算结果文件使用MOMAP软件进行磷光辐射速率的计算。
基态优化
首先使用量化软件BDF进行 \(\rm Ir(ppy)_3\) 的基态S0和第一激发三重态T1的结构优化和频率计算。
准备 \(\rm Ir(ppy)_3\) 分子结构的xyz文件如下：
61
Ir -2.606160230000 -0.262817540000 0.032585640000
C -3.837298770000 2.407777490000 0.243683290000
C -1.553000180000 2.622608110000 0.521271830000
C -3.991643180000 3.786525490000 0.422094600000
C -4.929719550000 1.476909230000 -0.003856060000
C -1.634158680000 3.988444110000 0.709067180000
H -0.591261780000 2.126459120000 0.554253960000
C -2.889293280000 4.581787990000 0.656041260000
H -4.972006580000 4.231616950000 0.376099570000
C -6.263398690000 1.888616610000 -0.076321100000
H -0.744068550000 4.570116630000 0.889344200000
H -2.999939670000 5.647195600000 0.795286530000
C -5.623129880000 -0.778047400000 -0.424041550000
C -7.264449980000 0.973962040000 -0.319493370000
H -6.527094820000 2.928025340000 0.057085040000
C -6.940359500000 -0.364034500000 -0.495248330000
H -5.397358100000 -1.824987250000 -0.570011260000
H -8.293855190000 1.298147200000 -0.375111620000
H -7.723296650000 -1.084118350000 -0.689505650000
C -2.780095460000 -2.271307610000 -0.073978550000
C -2.962145630000 -2.905938610000 1.179649240000
C -2.704720750000 -3.101798200000 -1.194651040000
C -3.053849400000 -4.297470620000 1.273150410000
C -3.052375550000 -2.037406190000 2.345262930000
C -2.797310890000 -4.477732940000 -1.095358170000
H -2.574810200000 -2.657655530000 -2.171485550000
C -2.970830350000 -5.080470100000 0.142706420000
H -3.190782260000 -4.777448690000 2.231582660000
C -3.251325250000 -2.483356810000 3.656088400000
H -2.735191030000 -5.089028800000 -1.985144810000
H -3.042943740000 -6.155896270000 0.221114210000
C -2.995026490000 0.148733040000 3.092507280000
C -3.319187740000 -1.576319800000 4.692855730000
H -3.353725050000 -3.536644150000 3.859562860000
C -3.186631940000 -0.222801710000 4.408749290000
H -2.888815760000 1.194668600000 2.833837140000
H -3.472767910000 -1.912494740000 5.707724440000
H -3.232610730000 0.522480480000 5.186904590000
N -2.927015180000 -0.717593640000 2.090156960000
C 0.125501920000 -0.309636820000 -1.075330180000
C 0.242318970000 -0.710559210000 1.197687130000
C 1.518232510000 -0.401785790000 -1.166132500000
C -0.762473880000 -0.044853840000 -2.198985150000
C 1.619464440000 -0.811598610000 1.179643800000
H -0.298309850000 -0.828666640000 2.128258480000
C 2.270161400000 -0.653514600000 -0.037593100000
H 2.007628100000 -0.277927520000 -2.118201950000
C -2.150460220000 -0.005742870000 -1.917121780000
C -0.287950380000 0.157363030000 -3.497950870000
H 2.165728000000 -1.009392350000 2.088259610000
H 3.346017750000 -0.726918300000 -0.099292420000
C -3.004998770000 0.237880450000 -2.994864780000
C -1.165335710000 0.397865570000 -4.532504760000
H 0.771521710000 0.127613130000 -3.708268180000
C -2.529020200000 0.436535570000 -4.277665910000
H -4.071982740000 0.267717660000 -2.824596670000
H -0.792761550000 0.552584980000 -5.535050450000
H -3.220629300000 0.622039860000 -5.087916480000
N -0.487689310000 -0.467248190000 0.117089220000
N -2.608631830000 1.850704260000 0.298598520000
C -4.578169480000 0.115167240000 -0.176155400000
打开Device studio，点击File-new project，命名为 phosphorescence.hpf ，将 Ir(ppy)3_s0.xyz 拖入Project中，双击 Ir(ppy)3_s0.hzw 。接下来进 \(\rm Ir(ppy)_3\) 基态S0的结构优化和频率计算。选中Simulator → BDF → BDF，在界面中设置参数。在Basic Settings界面中的Calculation Type选择Opt+Freq，方法选择PBE0泛函，基组在Basis中的All Electron类型中，选择Def2-SVP。Basic Settings界面中的其它参数使用推荐的默认值，不需要做修改。
在SCF Settings界面中，DFT Integral Grid选择Coarse，Convergence Threshold选择Tight。SCF Settings界面中的其它参数使用推荐的默认值，不需要做修改。
在OPT Settings界面中，Convergence Threshold选择Tight。OPT Settings界面中的其它参数使用推荐的默认值，不需要做修改。
Freq Settings界面中的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。选中生成的 bdf.inp 文件，右击选择open with，打开该文件，如下所示：
$compass
Title
C33H24IrN3
Geometry
Ir -2.60616023 -0.26281754 0.03258564
C -3.83729877 2.40777749 0.24368329
C -1.55300018 2.62260811 0.52127183
C -3.99164318 3.78652549 0.42209460
C -4.92971955 1.47690923 -0.00385606
C -1.63415868 3.98844411 0.70906718
H -0.59126178 2.12645912 0.55425396
C -2.88929328 4.58178799 0.65604126
H -4.97200658 4.23161695 0.37609957
C -6.26339869 1.88861661 -0.07632110
H -0.74406855 4.57011663 0.88934420
H -2.99993967 5.64719560 0.79528653
C -5.62312988 -0.77804740 -0.42404155
C -7.26444998 0.97396204 -0.31949337
H -6.52709482 2.92802534 0.05708504
C -6.94035950 -0.36403450 -0.49524833
H -5.39735810 -1.82498725 -0.57001126
H -8.29385519 1.29814720 -0.37511162
H -7.72329665 -1.08411835 -0.68950565
C -2.78009546 -2.27130761 -0.07397855
C -2.96214563 -2.90593861 1.17964924
C -2.70472075 -3.10179820 -1.19465104
C -3.05384940 -4.29747062 1.27315041
C -3.05237555 -2.03740619 2.34526293
C -2.79731089 -4.47773294 -1.09535817
H -2.57481020 -2.65765553 -2.17148555
C -2.97083035 -5.08047010 0.14270642
H -3.19078226 -4.77744869 2.23158266
C -3.25132525 -2.48335681 3.65608840
H -2.73519103 -5.08902880 -1.98514481
H -3.04294374 -6.15589627 0.22111421
C -2.99502649 0.14873304 3.09250728
C -3.31918774 -1.57631980 4.69285573
H -3.35372505 -3.53664415 3.85956286
C -3.18663194 -0.22280171 4.40874929
H -2.88881576 1.19466860 2.83383714
H -3.47276791 -1.91249474 5.70772444
H -3.23261073 0.52248048 5.18690459
N -2.92701518 -0.71759364 2.09015696
C 0.12550192 -0.30963682 -1.07533018
C 0.24231897 -0.71055921 1.19768713
C 1.51823251 -0.40178579 -1.16613250
C -0.76247388 -0.04485384 -2.19898515
C 1.61946444 -0.81159861 1.17964380
H -0.29830985 -0.82866664 2.12825848
C 2.27016140 -0.65351460 -0.03759310
H 2.00762810 -0.27792752 -2.11820195
C -2.15046022 -0.00574287 -1.91712178
C -0.28795038 0.15736303 -3.49795087
H 2.16572800 -1.00939235 2.08825961
H 3.34601775 -0.72691830 -0.09929242
C -3.00499877 0.23788045 -2.99486478
C -1.16533571 0.39786557 -4.53250476
H 0.77152171 0.12761313 -3.70826818
C -2.52902020 0.43653557 -4.27766591
H -4.07198274 0.26771766 -2.82459667
H -0.79276155 0.55258498 -5.53505045
H -3.22062930 0.62203986 -5.08791648
N -0.48768931 -0.46724819 0.11708922
N -2.60863183 1.85070426 0.29859852
C -4.57816948 0.11516724 -0.17615540
End Geometry
Basis
Def2-SVP
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
366
TolGrad
3.0E-5
TolStep
1.2E-4
IOpt
3
Hess
final
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
PBE0
D3
Grid
Coarse
ThreshConv
1.0D-9 5.0D-7
MPEC+COSX
Molden
$end
$resp
Geom
$end
选中 bdf.inp 文件，右击选择Run，弹出如下服务器提交的界面：
点击Run提交作业。任务结束后 bdf.out ， bdf.out.tmp ， bdf.scf.molden 三个结果文件会出现在Project中。
选择 bdf.out ，右击show view，在Optimization对话框中，显示结构已经达到收敛标准。
在Frequency对话框中，检查频率，若不存在虚频证明结构已经优化到极小点。
激发态优化
选择 bdf.out 文件，右击open with containing folder打开所在文件夹，在 bdf.out 文件中查找 converged in ，紧接着输出的 Molecular Cartesian Coordinates (X,Y,Z) in Angstrom : 下的结构即为优化好的 \(\rm Ir(ppy)_3\) 的S0结构。
将其保存为 Irppy3_t1.xyz 文件，将 Irppy3_t1.xyz 拖入Device Studio中进行T1激发态的结构优化和频率计算。
Irppy3_t1.xyz 内容如下：
61
Ir -0.00021963 0.00084588 0.01424181
C 2.59517396 -1.31710199 -0.58086411
C 2.23709967 0.40664133 -2.11684705
C 3.82729349 -1.60375453 -1.18851600
C 2.03843393 -2.01080680 0.57861773
C 3.44334868 0.17103124 -2.75937571
H 1.56522101 1.20579483 -2.43942631
C 4.25160770 -0.86138490 -2.27959559
H 4.44860577 -2.40719663 -0.79331056
C 2.69382363 -3.08153995 1.20802708
H 3.74085930 0.78654308 -3.60925966
H 5.21146469 -1.08097154 -2.75293386
C 0.24421139 -2.16970311 2.17811922
C 2.12763720 -3.69300459 2.31682204
H 3.65478554 -3.44259873 0.83261331
C 0.89831764 -3.22978876 2.79882363
H -0.71249803 -1.82386403 2.57651491
H 2.63779958 -4.52517522 2.80699129
H 0.44660698 -3.70582388 3.67403286
C -1.72035469 0.07933387 1.04722001
C -2.76313413 -0.76101290 0.56881686
C -2.01025266 0.87257612 2.17113445
C -4.02037491 -0.79383502 1.19368759
C -2.43582629 -1.59048558 -0.58889316
C -3.25751526 0.83538180 2.78746398
H -1.23410642 1.52839366 2.57249446
C -4.27126073 0.00210869 2.30157161
H -4.81180920 -1.44586060 0.81528272
C -3.29617560 -2.51740929 -1.19724703
H -3.44731484 1.46422538 3.66217358
H -5.24881003 -0.02406750 2.78811438
C -0.75837785 -2.14128396 -2.11906742
C -2.86048279 -3.25738595 -2.28544733
H -4.30328496 -2.65589425 -0.80436165
C -1.56140370 -3.07140479 -2.76199645
H 0.27022263 -1.95736074 -2.43892538
H -3.52729093 -3.98185076 -2.75888376
H -1.17316801 -3.63790053 -3.60937499
N -1.18202102 -1.42059810 -1.07683756
C -0.15520348 2.90563657 -0.58971088
C -1.46338895 1.72871866 -2.12679953
C -0.52211987 4.11414752 -1.20165791
C 0.72001848 2.77388462 0.57312738
C -1.86121725 2.88919512 -2.77366419
H -1.81791319 0.74611096 -2.44764729
C -1.37377275 4.10694585 -2.29536844
H -0.13934280 5.05543436 -0.80755998
C 0.92315032 1.45037015 1.05224777
C 1.31701199 3.87857243 1.20194236
H -2.54044030 2.83638038 -3.62545723
H -1.66304625 5.04659184 -2.77183207
C 1.74847348 1.30397901 2.18095829
C 2.12513377 3.69714085 2.31457434
H 1.15050564 4.89023888 0.82308669
C 2.33662182 2.40226678 2.80117156
H 1.92559465 0.30363467 2.58285826
H 2.58863363 4.55643063 2.80428368
H 2.97087444 2.25159426 3.67958431
N -0.63391888 1.73510439 -1.07938881
N 1.82331270 -0.31615001 -1.07209969
C 0.78888390 -1.52564414 1.05342887
选择Simulator → BDF → BDF，在界面中设置参数。在Basic Settings界面中的Calculation Type选择TDDFT-OPT+Freq，方法采用默认的PBE0泛函，基组在Basis中的All Electron类型中，选择Def2-SVP。Basic Settings界面中的其它参数使用推荐的默认值，不需要做修改。
SCF Settings界面中将Use MPEC+COSX Acceleraton的默认勾选去掉。SCF Settings界面中的其它参数使用推荐的默认值，不需要做修改。
TDDFT Settings面板中将Use MPEC+COSX Acceleraton的默认勾选去掉。Mutiplicity选择Triplet，Convergence Threshold选择Very Tight。TDDFT Settings界面中的其它参数使用推荐的默认值，不需要做修改。
OPT Settings和Freq Settings面板的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。选中生成的 bdf.inp 文件，右击选择open containing folder进入所在文件夹，在 bdf.inp 的 $tddft 模块中加入：
Gridtol
1E-6
bdf.inp 内容如下：
$compass
Title
C33H24IrN3
Geometry
Ir -0.00021963 0.00084588 0.01424181
C 2.59517396 -1.31710199 -0.58086411
C 2.23709967 0.40664133 -2.11684705
C 3.82729349 -1.60375453 -1.18851600
C 2.03843393 -2.01080680 0.57861773
C 3.44334868 0.17103124 -2.75937571
H 1.56522101 1.20579483 -2.43942631
C 4.25160770 -0.86138490 -2.27959559
H 4.44860577 -2.40719663 -0.79331056
C 2.69382363 -3.08153995 1.20802708
H 3.74085930 0.78654308 -3.60925966
H 5.21146469 -1.08097154 -2.75293386
C 0.24421139 -2.16970311 2.17811922
C 2.12763720 -3.69300459 2.31682204
H 3.65478554 -3.44259873 0.83261331
C 0.89831764 -3.22978876 2.79882363
H -0.71249803 -1.82386403 2.57651491
H 2.63779958 -4.52517522 2.80699129
H 0.44660698 -3.70582388 3.67403286
C -1.72035469 0.07933387 1.04722001
C -2.76313413 -0.76101290 0.56881686
C -2.01025266 0.87257612 2.17113445
C -4.02037491 -0.79383502 1.19368759
C -2.43582629 -1.59048558 -0.58889316
C -3.25751526 0.83538180 2.78746398
H -1.23410642 1.52839366 2.57249446
C -4.27126073 0.00210869 2.30157161
H -4.81180920 -1.44586060 0.81528272
C -3.29617560 -2.51740929 -1.19724703
H -3.44731484 1.46422538 3.66217358
H -5.24881003 -0.02406750 2.78811438
C -0.75837785 -2.14128396 -2.11906742
C -2.86048279 -3.25738595 -2.28544733
H -4.30328496 -2.65589425 -0.80436165
C -1.56140370 -3.07140479 -2.76199645
H 0.27022263 -1.95736074 -2.43892538
H -3.52729093 -3.98185076 -2.75888376
H -1.17316801 -3.63790053 -3.60937499
N -1.18202102 -1.42059810 -1.07683756
C -0.15520348 2.90563657 -0.58971088
C -1.46338895 1.72871866 -2.12679953
C -0.52211987 4.11414752 -1.20165791
C 0.72001848 2.77388462 0.57312738
C -1.86121725 2.88919512 -2.77366419
H -1.81791319 0.74611096 -2.44764729
C -1.37377275 4.10694585 -2.29536844
H -0.13934280 5.05543436 -0.80755998
C 0.92315032 1.45037015 1.05224777
C 1.31701199 3.87857243 1.20194236
H -2.54044030 2.83638038 -3.62545723
H -1.66304625 5.04659184 -2.77183207
C 1.74847348 1.30397901 2.18095829
C 2.12513377 3.69714085 2.31457434
H 1.15050564 4.89023888 0.82308669
C 2.33662182 2.40226678 2.80117156
H 1.92559465 0.30363467 2.58285826
H 2.58863363 4.55643063 2.80428368
H 2.97087444 2.25159426 3.67958431
N -0.63391888 1.73510439 -1.07938881
N 1.82331270 -0.31615001 -1.07209969
C 0.78888390 -1.52564414 1.05342887
End Geometry
Basis
Def2-SVP
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
366
IOpt
3
Hess
final
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
PBE0
D3
Molden
$end
$tddft
Imethod
1
Isf
1
Ialda
4
Idiag
1
crit_e
1E-9
crit_vec
1E-7
Gridtol
1E-6
Iroot
6
Istore
1
$end
$resp
Geom
Method
2
Nfiles
1
Iroot
1
$end
选中 bdf.inp 文件，右击选择Run提交作业，任务结束后 bdf.out ， bdf.out.tmp ， bdf.scf.molden 三个结果文件会出现在Project中。
选择 bdf.out，右击show view，在Optimization对话框中，显示结构已经达到收敛标准。
在Frequency对话框中，检查频率，若不存在虚频证明结构已经优化到极小点。
自旋轨道耦合
选择 bdf.out 文件，右击open with containing folder打开 bdf.out 文件，在文件中查找 converged in ，紧接着输出的 Molecular Cartesian Coordinates (X,Y,Z) in Angstrom : 下的结构即为优化好的T1激发态的结构。将其保存为 Irppy3_t1_soc.xyz 文件，如下：
61
Ir 0.00713728 0.02772384 0.06844143
C 2.49525480 -1.44901550 -0.61634342
C 2.18832036 0.30085414 -2.14613716
C 3.68634391 -1.80881598 -1.26572189
C 1.93194560 -2.11689508 0.55823360
C 3.35838993 -0.00562745 -2.82371008
H 1.54555778 1.13535499 -2.43828057
C 4.11644204 -1.08671357 -2.36826138
H 4.27131595 -2.65056635 -0.89578008
C 2.53568350 -3.23676194 1.15281696
H 3.66807720 0.59265321 -3.68133338
H 5.04582829 -1.36185464 -2.87261245
C 0.15985754 -2.20796739 2.19060975
C 1.95468524 -3.83725789 2.26057143
H 3.46642195 -3.64596143 0.75209976
C 0.76249168 -3.31842903 2.77624738
H -0.76777546 -1.81381956 2.61329026
H 2.42616559 -4.70662836 2.72403491
H 0.30108846 -3.78788395 3.64972556
C -1.72817262 0.21988877 1.05055833
C -2.80684294 -0.57231379 0.57552059
C -1.98377974 1.07446425 2.13652018
C -4.07348284 -0.50293868 1.17614116
C -2.51722058 -1.44616477 -0.55935718
C -3.24105830 1.13344573 2.72846833
H -1.17254968 1.69178400 2.52835606
C -4.29332764 0.34509124 2.25152759
H -4.89835192 -1.11318656 0.80076906
C -3.42583031 -2.33456216 -1.15446766
H -3.40666531 1.80444609 3.57586126
H -5.27935386 0.39610056 2.71819787
C -0.85701735 -2.13799807 -2.04878703
C -3.02057249 -3.12865177 -2.21547404
H -4.44525951 -2.39959512 -0.77498376
C -1.70631730 -3.03592702 -2.67708276
H 0.18295061 -2.02278722 -2.36320871
H -3.72428268 -3.82273458 -2.68079360
H -1.34337957 -3.64618311 -3.50492143
N -1.25281509 -1.36491844 -1.03498749
C 0.05749757 2.91146589 -0.57266019
C -1.32777267 1.80183369 -2.13392316
C -0.20378718 4.13789922 -1.23242993
C 0.84833732 2.74053836 0.60027468
C -1.62207961 2.97568834 -2.79963589
H -1.76529075 0.85235254 -2.45705604
C -1.02279372 4.18710974 -2.33345871
H 0.25619858 5.05119986 -0.85064151
C 0.99228869 1.37116718 1.10523883
C 1.50408647 3.78492492 1.29091761
H -2.29824567 2.96275979 -3.65460398
H -1.21968527 5.13470890 -2.83803374
C 1.79964051 1.14876808 2.23660071
C 2.27478596 3.51131149 2.40946143
H 1.40861651 4.81693356 0.94742301
C 2.43450283 2.19478112 2.89597173
H 1.90681895 0.12796182 2.60984200
H 2.77105979 4.33756352 2.92655136
H 3.04508360 2.00761950 3.78145403
N -0.50508694 1.73366277 -1.08285478
N 1.77567220 -0.40171722 -1.08777429
C 0.72548984 -1.57229627 1.07484739
基于T1激发态结构，做S0和T1之间的旋轨耦合SOC计算。将Irppy3_t1_soc.xyz拖入Device Studio中，选择Simulator → BDF → BDF，界面中设置参数。在Basic Settings界面中的Calculation Type选择TDDFT-SOC，Functional选择PBE0泛函，取消勾选Use Dispersion Correction。Hamiltonian选择sf-X2C，Basis选择All Electron类型下的x2c-SVPall基组。
在TDDFT Settings面板中的Spin-Orbit Coupling内容框中，勾选Including Ground State。
Basic Settings、SCF Settings和TDDFT Settings界面中的其它参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。选中生成的 bdf.inp 文件，右击选择open with，打开该文件，如下所示：
生成的输入文件 bdf.inp 为：
$compass
Title
C33H24IrN3
Geometry
Ir 0.00713728 0.02772384 0.06844143
C 2.49525480 -1.44901550 -0.61634342
C 2.18832036 0.30085414 -2.14613715
C 3.68634391 -1.80881598 -1.26572189
C 1.93194560 -2.11689507 0.55823360
C 3.35838993 -0.00562745 -2.82371008
H 1.54555778 1.13535498 -2.43828056
C 4.11644203 -1.08671356 -2.36826137
H 4.27131594 -2.65056635 -0.89578008
C 2.53568350 -3.23676194 1.15281696
H 3.66807719 0.59265321 -3.68133337
H 5.04582829 -1.36185464 -2.87261245
C 0.15985754 -2.20796738 2.19060975
C 1.95468524 -3.83725789 2.26057143
H 3.46642195 -3.64596142 0.75209976
C 0.76249168 -3.31842902 2.77624738
H -0.76777546 -1.81381956 2.61329025
H 2.42616559 -4.70662835 2.72403490
H 0.30108846 -3.78788394 3.64972555
C -1.72817262 0.21988877 1.05055833
C -2.80684294 -0.57231379 0.57552059
C -1.98377973 1.07446424 2.13652018
C -4.07348283 -0.50293868 1.17614116
C -2.51722058 -1.44616477 -0.55935718
C -3.24105830 1.13344573 2.72846833
H -1.17254967 1.69178400 2.52835606
C -4.29332764 0.34509124 2.25152758
H -4.89835191 -1.11318655 0.80076906
C -3.42583030 -2.33456215 -1.15446765
H -3.40666531 1.80444608 3.57586125
H -5.27935385 0.39610056 2.71819787
C -0.85701735 -2.13799807 -2.04878703
C -3.02057249 -3.12865176 -2.21547404
H -4.44525950 -2.39959511 -0.77498376
C -1.70631730 -3.03592701 -2.67708275
H 0.18295061 -2.02278722 -2.36320870
H -3.72428268 -3.82273458 -2.68079359
H -1.34337957 -3.64618310 -3.50492143
N -1.25281508 -1.36491844 -1.03498749
C 0.05749757 2.91146589 -0.57266019
C -1.32777267 1.80183369 -2.13392316
C -0.20378718 4.13789922 -1.23242992
C 0.84833732 2.74053836 0.60027467
C -1.62207960 2.97568834 -2.79963588
H -1.76529074 0.85235254 -2.45705603
C -1.02279371 4.18710974 -2.33345870
H 0.25619858 5.05119985 -0.85064151
C 0.99228869 1.37116718 1.10523883
C 1.50408647 3.78492491 1.29091760
H -2.29824567 2.96275978 -3.65460398
H -1.21968527 5.13470889 -2.83803374
C 1.79964050 1.14876807 2.23660070
C 2.27478596 3.51131149 2.40946142
H 1.40861651 4.81693355 0.94742301
C 2.43450283 2.19478112 2.89597172
H 1.90681894 0.12796182 2.60984200
H 2.77105978 4.33756351 2.92655136
H 3.04508359 2.00761950 3.78145402
N -0.50508694 1.73366277 -1.08285477
N 1.77567220 -0.40171722 -1.08777429
C 0.72548984 -1.57229627 1.07484738
End Geometry
Basis
x2c-SVPall
Skeleton
Group
C(1)
$end
$xuanyuan
Heff
21
Hsoc
2
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
PBE0
MPEC+COSX
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
6
MPEC+COSX
Istore
1
$end
$tddft
Imethod
1
Isf
1
Idiag
1
Iroot
6
MPEC+COSX
Istore
2
$end
$tddft
Isoc
2
Nfiles
2
Ifgs
1
Imatsoc
-1
Imatrsf
-1
Imatrso
-2
$end
选中 bdf.inp 文件，右击选择Run提交作业，任务结束后 bdf.out ， bdf.scf.molden 结果文件会出现在Project中。选择 bdf.out ，右击选择show view，在TDDFT面板中，选择Spinor，在Dominant Excitations中确定第2、3和4为T1态的三个分量。
点击bdf.out文件，右击选择Open Containing Folder进入所在文件夹，打开bdf.out，查找 *** List of SOC-SI results *** ，读取第2、3和4态的ExEnergies，分别为：2.1906 eV，2.1961 eV和2.2052 eV。
*** List of SOC-SI results ***
No. ExEnergies Dominant Excitations Esf dE Eex(eV) (cm^-1)
1 -0.0054 eV 99.8% Spin: |Gs,1> 0-th A 0.0000 -0.0054 0.0000 0.00
2 2.1906 eV 43.5% Spin: |S+,3> 1-th A 2.2232 -0.0326 2.1961 17712.45
3 2.1961 eV 75.0% Spin: |S+,2> 1-th A 2.2232 -0.0272 2.2015 17756.09
4 2.2052 eV 42.1% Spin: |S+,1> 1-th A 2.2232 -0.0180 2.2106 17829.67
5 2.5334 eV 49.1% Spin: |So,1> 1-th A 2.6854 -0.1520 2.5388 20477.15
6 2.5861 eV 42.4% Spin: |S+,3> 2-th A 2.6312 -0.0452 2.5915 20901.71
7 2.6064 eV 82.9% Spin: |S+,2> 2-th A 2.6312 -0.0248 2.6118 21065.69
查找 E_tot ，读取相应的能量为 -19265.29575859。T1的三个子态的能量分别为E_tot的能量与ExEnergies的能量相加，以第2个子态为例，计算方法为：-19265.29575859+2.1906/27.2114=-19265.215256 au。第3个子态的能量为-19265.215053 au。第4个子态的能量为-19265.214719 au。
Final scf result
E_tot = -19265.29575859
E_ele = -25841.98940694
E_nn = 6576.69364834
E_1e = -39510.05277256
E_ne = -66428.66809936
E_kin = 26918.61532681
E_ee = 14091.21945939
E_xc = -423.15609377
Virial Theorem 1.715687
使用上述相同的方法和基组，以S0基态的结构计算SOC。在bdf.out中查找 E_tot ，读取相应的能量为：-19265.30415493 au。
Final scf result
E_tot = -19265.30415493
E_ele = -25838.09048037
E_nn = 6572.78632544
E_1e = -39502.28526599
E_ne = -66421.04136762
E_kin = 26918.75610162
E_ee = 14087.38176801
E_xc = -423.18698239
Virial Theorem 1.715683
T1的三个子态相对于S0态的能量为上述三个子态能量减去S0态的能量，分别为：0.088899 au，0.089102 au，0.089436 au。在T1结构的SOC计算的out文件中查找 [tddft_soc_matrso]: ，输出了考虑SOC之后各激发态相对于基态的能量和跃迁偶极矩。
[tddft_soc_matrso]: Print selected matrix elements of [dpl]
No. ( I , J ) |rij|^2 E_J-E_I fosc rate(s^-1)
-------------------------------------------------------------------------------
1 1 2 0.104E-02 2.196064924 0.000056149 0.117E+05
Details of transition dipole moment with SOC (in a.u.):
<I|X|J> <I|Y|J> <I|Z|J> (also in debye)
Real= 0.279E-01 0.161E-01 -0.216E-02 0.0710 0.0409 -0.0055
Imag= 0.123E-07 -0.291E-07 -0.867E-08 0.0000 -0.0000 -0.0000
Norm= 0.279E-01 0.161E-01 0.216E-02
No. ( I , J ) |rij|^2 E_J-E_I fosc rate(s^-1)
-------------------------------------------------------------------------------
2 1 3 0.354E-03 2.201474871 0.000019090 0.401E+04
Details of transition dipole moment with SOC (in a.u.):
<I|X|J> <I|Y|J> <I|Z|J> (also in debye)
Real= 0.587E-02 0.179E-01 0.126E-03 0.0149 0.0454 0.0003
Imag= -0.108E-06 -0.357E-07 0.361E-07 -0.0000 -0.0000 0.0000
Norm= 0.587E-02 0.179E-01 0.126E-03
No. ( I , J ) |rij|^2 E_J-E_I fosc rate(s^-1)
-------------------------------------------------------------------------------
3 1 4 0.259E-01 2.210597826 0.001400915 0.297E+06
Details of transition dipole moment with SOC (in a.u.):
<I|X|J> <I|Y|J> <I|Z|J> (also in debye)
Real= 0.905E-08 -0.356E-07 -0.418E-08 0.0000 -0.0000 -0.0000
Imag= -0.535E-01 0.148E+00 0.316E-01 -0.1360 0.3771 0.0802
Norm= 0.535E-01 0.148E+00 0.316E-01
No. ( I , J ) |rij|^2 E_J-E_I fosc rate(s^-1)
-------------------------------------------------------------------------------
4 1 5 0.154E+00 2.538843563 0.009594212 0.268E+07
Details of transition dipole moment with SOC (in a.u.):
<I|X|J> <I|Y|J> <I|Z|J> (also in debye)
Real= -0.236E+00 0.158E+00 0.271E+00 -0.5998 0.4010 0.6899
Imag= -0.271E-06 0.183E-06 0.310E-06 -0.0000 0.0000 0.0000
Norm= 0.236E+00 0.158E+00 0.271E+00
No. ( I , J ) |rij|^2 E_J-E_I fosc rate(s^-1)
-------------------------------------------------------------------------------
5 1 6 0.275E-02 2.591483156 0.000174312 0.508E+05
Details of transition dipole moment with SOC (in a.u.):
<I|X|J> <I|Y|J> <I|Z|J> (also in debye)
Real= 0.339E-01 0.292E-01 0.273E-01 0.0861 0.0743 0.0693
Imag= -0.132E-07 0.447E-07 0.428E-07 -0.0000 0.0000 0.0000
Norm= 0.339E-01 0.292E-01 0.273E-01
其中， 1 2 表示第一个与第二个旋量态间的跃迁偶极矩，以此类推。这里需要第1、2和3激发态的激发能和跃迁偶极矩数据。
跃迁偶极矩的数据在 Details of transition dipole moment with SOC (in a.u.): 中列出，前三列为单位为au的偶极矩数据，后三列为单位为Debye的偶极矩数据。
将debye单位下六个数的平方之和开方，得到该态的跃迁偶极矩。另一种方法为将Norm后的三个数平方之和开方，再乘以2.5417。得到第1、2和3激发态的跃迁偶极矩为0.082058 debye, 0.047881 Debye, 0.407979 Debye。
得到的这6个参数将用于MOMAP软件的磷光发射速率的计算中。
至此，MOMAP对 \(\rm Ir(ppy)_3\) 的磷光辐射速率计算需要的BDF量化软件的结构优化频率结果文件、旋轨耦合计算结果文件和参数部分都已完成。
磷光辐射速率
接下来开始使用MOMAP进行 \(\rm Ir(ppy)_3\) 的磷光辐射速率的计算。
首先计算T1→S0的磷光辐射速率，第一步为做电子振动耦合（electron-vibration coupling, EVC）计算，该计算基于量化计算输出的分子振动频率、力常数矩阵，同时在内坐标以及直角坐标系下，计算分子跃迁发生初末态间的模式位移、黄昆因子、重整能以及Duschinsky转动矩阵。
将BDF软件计算得到的 \(\rm Ir(ppy)_3\) 的基态的优化频率计算文件改名为 irppy3_s0.out ，和T1优化频率计算文件改名为 irppy3_t1.out ，同时放在EVC计算文件夹中。
EVC计算的输入文件 momap.inp 为：
do_evc=1
&evc
ffreq(1)="irppy3_s0.out"
ffreq(2)="irppy3_t1.out"
proj_reorg=.t.
/
提交脚本文件 momap.slurm 运行任务。任务正常结束后，使用 evc.cart.dat 进行接下来的T1→S0磷光发射速率计算。对T1的三个态依次进行计算，其中第一个子态，即第2个态的磷光发射速率，输入文件 momap.inp 为:
do_spec_tvcf_ft = 1
do_spec_tvcf_spec = 1
&spec_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.088899 au
EDMA = 1 debye
EDME = 0.082058 debye
FreqScale = 1
DSFile = "evc.cart.dat"
Emax = 0.3 au
dE = 0.00001 au
logFile = "spec.tvcf.log"
FtFile = "spec.tvcf.ft.dat"
FoFile = "spec.tvcf.fo.dat"
FoSFile = "spec.tvcf.spec.dat"
/
提交脚本文件 momap.slurm 运行任务。任务结束后，确认关联函数是否收敛。
打开 spec.tvcf.log ，磷光辐射速率在第一个数和第二个数读取，单位分别为a.u.和s-1，第三个数为寿命，单位为ns。
第二个子态，即第3个态的磷光发射速率，输入文件 momap.inp 为:
do_spec_tvcf_ft = 1
do_spec_tvcf_spec = 1
&spec_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.089102 au
EDMA = 1 debye
EDME = 0.047881 debye
FreqScale = 1
DSFile = "evc.cart.dat"
Emax = 0.3 au
dE = 0.00001 au
logFile = "spec.tvcf.log"
FtFile = "spec.tvcf.ft.dat"
FoFile = "spec.tvcf.fo.dat"
FoSFile = "spec.tvcf.spec.dat"
/
提交脚本文件 momap.slurm 运行任务。
第三个子态，即第4个态的磷光发射速率，输入文件 momap.inp 为:
do_spec_tvcf_ft = 1
do_spec_tvcf_spec = 1
&spec_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.089436 au
EDMA = 1 debye
EDME = 0.407979 debye
FreqScale = 1
DSFile = "evc.cart.dat"
Emax = 0.3 au
dE = 0.00001 au
logFile = "spec.tvcf.log"
FtFile = "spec.tvcf.ft.dat"
FoFile = "spec.tvcf.fo.dat"
FoSFile = "spec.tvcf.spec.dat"
/
根据三个子态相对能量的玻尔兹曼分布（可参考http://sobereva.com/462和http://sobereva.com/165中的介绍以及相应的excel表格）对三个子态的磷光发射速率乘以相应的比例，加和，最终得到T1态的磷光发射速率。
量化理论计算探究薁(azulene)的反Kasha规则荧光机制
根据Kasha规则，由于高激发态之间存在快速的无辐射跃迁，分子的荧光或磷光的初始状态是最低的单重态或三重态。薁(azulene)是典型违反 Kasha 规则的例子。azulene的荧光来自于S2态，可以认为是由于S2→S1能隙比较大，所以降低了S2→S1内转换速率。另外，由于S1和S0之间的能隙相对较小，所以其内转换速率很大，从而降低了S1→S0的荧光量子效率，所以S1→S0的荧光很难被发现。这里以具有反常荧光现象的azulene为例，使用BDF软件和MOMAP软件，计算azulene的S1→S0的辐射速率和内转换速率，从而解释azulene第一激发态极低的量子效率导致其荧光难以被观测到的实验结果。
MOMAP对azulene的S1→S0的辐射速率和内转换速率的计算需要BDF量化软件的结构优化频率结果文件、非绝热耦合结果文件和参数部分都已完成。首先完成量化软件BDF的计算部分。
BDF计算部分
准备azulene分子结构的xyz文件如下：
18
C -0.48100000 0.74480000 0.00000000
C -0.56240000 -0.71320000 0.00020000
C -1.75790000 1.17860000 -0.00030000
C -1.96510000 -1.08880000 0.00000000
C 0.66870000 1.60890000 0.00030000
C 0.45100000 -1.58850000 0.00030000
C -2.66930000 0.05180000 -0.00010000
C 1.95140000 1.22210000 0.00020000
C 1.86730000 -1.29960000 -0.00020000
C 2.49720000 -0.11610000 -0.00040000
H -2.09080000 2.20560000 -0.00040000
H -2.35750000 -2.09240000 0.00010000
H 0.46620000 2.67860000 0.00070000
H 0.22090000 -2.65340000 0.00040000
H -3.74370000 0.14050000 -0.00030000
H 2.70720000 2.00650000 0.00050000
H 2.49320000 -2.19150000 -0.00060000
H 3.58670000 -0.13930000 -0.00080000
打开Device studio，点击File-new project，命名为 fluorescence.hpf ，将 azulene.xyz 拖入Project中，双击 azulene.hzw ，得到如图所示界面。
首先使用BDF进行azulene的结构优化和频率计算。选中Simulator → BDF → BDF，界面中设置参数。在Basic Settings界面中的Calculation Type选择Opt+Freq，方法采用默认的GB3LYP泛函，基组在Basis中的All Electron类型中，选择6-31G(d,p)。Basic Settings界面中的其它参数以及SCF Settings、OPT Settings、Freq Settings等面板的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。
选中生成的 bdf.inp 文件，右击选择open with，打开该文件，如下所示：
$compass
Title
C10H8
Geometry
C -0.48100000 0.74480000 0.00000000
C -0.56240000 -0.71320000 0.00020000
C -1.75790000 1.17860000 -0.00030000
C -1.96510000 -1.08880000 0.00000000
C 0.66870000 1.60890000 0.00030000
C 0.45100000 -1.58850000 0.00030000
C -2.66930000 0.05180000 -0.00010000
C 1.95140000 1.22210000 0.00020000
C 1.86730000 -1.29960000 -0.00020000
C 2.49720000 -0.11610000 -0.00040000
H -2.09080000 2.20560000 -0.00040000
H -2.35750000 -2.09240000 0.00010000
H 0.46620000 2.67860000 0.00070000
H 0.22090000 -2.65340000 0.00040000
H -3.74370000 0.14050000 -0.00030000
H 2.70720000 2.00650000 0.00050000
H 2.49320000 -2.19150000 -0.00060000
H 3.58670000 -0.13930000 -0.00080000
End Geometry
Basis
6-31G(d,p)
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
108
IOpt
3
Hess
final
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
GB3LYP
D3
MPEC+COSX
Molden
$end
$resp
Geom
$end
选中 bdf.inp 文件，右击选择Run，弹出如下界面：
点击Run提交作业，会自动弹出结果文件的实时输出。
任务结束后 bdf.out ， bdf.out.tmp ， bdf.scf.molden 三个结果文件会出现在Project中。
选择 bdf.out ，右击show view，在Optimization对话框中，显示结构已经达到收敛标准。
在Frequency对话框中，检查频率，若不存在虚频证明结构已经优化到极小点。
在Summary对话框中，Total Energy为-385.87807600 a.u.，为需要的基态azulene的单点能。
点击Job Manager中该计算任务，点击服务器，已经进入到了该任务所在文件夹下，输入 /data/hzwtech/DS-BDF_2022A/sbin/optgeom2xyz.py bdf.optgeom ，回车，生成 bdf.xyz 文件。点击文件传输工具，进入文件夹下，将 bdf.xyz 文件拖出，即为下一步激发态结构优化需要的输入文件。改名为 azulene_s0.xyz ，打开文件夹，将第二行描述行去掉，拖入Device Studio中。
使用BDF进行azulene的S1激发态结构优化和频率计算。选中Simulator → BDF → BDF，界面中设置参数。在Basic Settings界面中的Calculation Type选择TDDFT-OPT+Freq，方法采用默认的GB3LYP泛函，基组在Basis中的All Electron类型中，选择6-31G(d,p)。SCF Settings和TDDFT Settings面板中将Use MPEC+COSX Acceleraton的默认勾选去掉。Basic Settings、SCF Settings、TDDFT Settings界面中的其它参数以及OPT Settings、Freq Settings等面板的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。
选中生成的 bdf.inp 文件，右击选择open with，打开该文件，如下所示：
$compass
Title
C10H8
Geometry
C 0.79273796 0.49102542 -0.00003307
C -0.70229649 0.61186591 0.00000000
C 1.30022932 1.80163337 -0.00006272
C -0.99262499 1.98726800 0.00007812
C 1.54415132 -0.67887418 0.00000000
C -1.63318173 -0.42094563 -0.00002837
C 0.21877157 2.69859813 0.00000000
C 1.10656346 -2.00562676 0.00005788
C -1.41619168 -1.80093044 -0.00004814
C -0.20258112 -2.49333483 0.00000000
H 2.35092512 2.06249889 -0.00009828
H -1.98777600 2.41348149 0.00017650
H 2.62424717 -0.53731745 0.00001117
H -2.67585843 -0.10561277 -0.00001521
H 0.30641472 3.77916915 0.00002386
H 1.88966566 -2.75951313 0.00017581
H -2.31053950 -2.41870505 -0.00009019
H -0.29054446 -3.57807510 0.00000000
End Geometry
Basis
6-31G(d,p)
Skeleton
Group
C(1)
$end
$bdfopt
Solver
1
MaxCycle
108
IOpt
3
Hess
final
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
GB3LYP
D3
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
6
Istore
1
$end
$resp
Geom
Method
2
Nfiles
1
Iroot
1
$end
选中 bdf.inp 文件，右击选择Run提交作业，任务结束后 bdf.out ， bdf.out.tmp ， bdf.scf.molden 三个结果文件会出现在Project中。
选择 bdf.out ，右击show view，在Optimization对话框中，显示结构已经达到收敛标准。
在Frequency对话框中，检查频率，若不存在虚频证明结构已经优化到极小点。
选择 bdf.out.tmp ，右击open with containing folder，打开 bdf.out.tmp ，在文件开始向下查找第一个tddft计算模块 module tddft 。根据tddft模块的 Ground to excited state Transition electric dipole moments (Au) 中的State 1的跃迁电偶极矩，得到momap需要的参数EDMA，计算方法为： \(\sqrt{(0.3456)^2+(-0.1159)^2+(-0.0000)^2}\) = 0.3646 a.u.。需要将单位a.u.转换为Debye，因此EDMA= 0.3646*2.5417=0.9267 Debye。
*** Ground to excited state Transition electric dipole moments (Au) ***
State X Y Z Osc.
1 0.3456 -0.1159 0.0000 0.0079 0.0079
2 0.0538 0.1576 0.0000 0.0025 0.0025
3 -0.6407 0.2166 0.0000 0.0527 0.0527
4 0.9123 2.7142 0.0000 1.0366 1.0366
5 -0.0001 0.0002 0.0200 0.0001 0.0001
6 -1.2021 0.4024 -0.0000 0.2383 0.2383
在文件末尾向上查找第一个tddft计算模块 module tddft 。根据tddft模块的 Ground to excited state Transition electric dipole moments (Au) 中的State 1的跃迁电偶极矩，得到momap需要的参数EDME，计算方法为： \(\sqrt{(-0.2427)^2+(0.0816)^2+(-0.0000)^2}\) = 0.2560 a.u.。需要将单位a.u.转换为Debye，因此EDMA= 0.3646*2.5417=0.6507 Debye。
*** Ground to excited state Transition electric dipole moments (Au) ***
State X Y Z Osc.
1 -0.2427 0.0816 0.0000 0.0026 0.0026
2 0.0403 0.1199 0.0000 0.0013 0.0013
3 -0.2655 0.0888 -0.0000 0.0090 0.0090
4 -0.8679 -2.5831 0.0000 0.8696 0.8696
5 -1.2356 0.4150 0.0000 0.2404 0.2404
6 -0.0008 0.0003 0.0006 0.0000 0.0000
在该 tddft 模块的 Statistics for [dvdson_rpa_block]: 中读取No. 1态的能量为 -385.8030480568 a.u.，即为S1态的单点能。
------------------------------------------------------------------
Statistics for [dvdson_rpa_block]:
No. of blocks = 1
Size of blocks = 50
No. of eigens = 6
No. of HxProd = 73 Averaged = 12.167
Eigenvalues (a.u.) =
0.0593694732 0.1241240214 0.1718082072
0.1756634611 0.2122947029 0.2131964457
------------------------------------------------------------------
No. 1 w= 1.6155 eV -385.8030480568 a.u. f= 0.0026 D<Pab>= 0.0000 Ova= 0.4881
CV(0): A( 33 )-> A( 36 ) c_i: -0.1361 Per: 1.9% IPA: 5.454 eV Oai: 0.7151
CV(0): A( 34 )-> A( 35 ) c_i: -0.9847 Per: 97.0% IPA: 2.582 eV Oai: 0.4815
No. 2 w= 3.3776 eV -385.7382935086 a.u. f= 0.0013 D<Pab>= 0.0000 Ova= 0.8097
CV(0): A( 33 )-> A( 35 ) c_i: -0.6839 Per: 46.8% IPA: 4.083 eV Oai: 0.8197
CV(0): A( 34 )-> A( 36 ) c_i: 0.7092 Per: 50.3% IPA: 3.953 eV Oai: 0.8078
No. 3 w= 4.6751 eV -385.6906093228 a.u. f= 0.0090 D<Pab>= 0.0000 Ova= 0.7195
CV(0): A( 32 )-> A( 35 ) c_i: 0.5413 Per: 29.3% IPA: 5.770 eV Oai: 0.7332
CV(0): A( 33 )-> A( 36 ) c_i: 0.8170 Per: 66.7% IPA: 5.454 eV Oai: 0.7151
CV(0): A( 34 )-> A( 38 ) c_i: 0.1417 Per: 2.0% IPA: 7.164 eV Oai: 0.7494
No. 4 w= 4.7800 eV -385.6867540689 a.u. f= 0.8696 D<Pab>= 0.0000 Ova= 0.7950
CV(0): A( 32 )-> A( 36 ) c_i: 0.1707 Per: 2.9% IPA: 7.141 eV Oai: 0.8644
CV(0): A( 33 )-> A( 35 ) c_i: -0.6794 Per: 46.2% IPA: 4.083 eV Oai: 0.8197
CV(0): A( 33 )-> A( 38 ) c_i: 0.1022 Per: 1.0% IPA: 8.665 eV Oai: 0.8000
CV(0): A( 34 )-> A( 36 ) c_i: -0.6312 Per: 39.8% IPA: 3.953 eV Oai: 0.8078
No. 5 w= 5.7768 eV -385.6501228271 a.u. f= 0.2404 D<Pab>= 0.0000 Ova= 0.7166
CV(0): A( 31 )-> A( 36 ) c_i: -0.1797 Per: 3.2% IPA: 8.008 eV Oai: 0.7623
CV(0): A( 32 )-> A( 35 ) c_i: -0.7750 Per: 60.1% IPA: 5.770 eV Oai: 0.7332
CV(0): A( 33 )-> A( 36 ) c_i: 0.4518 Per: 20.4% IPA: 5.454 eV Oai: 0.7151
CV(0): A( 34 )-> A( 38 ) c_i: 0.1740 Per: 3.0% IPA: 7.164 eV Oai: 0.7494
CV(0): A( 34 )-> A( 40 ) c_i: -0.2680 Per: 7.2% IPA: 8.035 eV Oai: 0.6274
No. 6 w= 5.8014 eV -385.6492210843 a.u. f= 0.0000 D<Pab>= 0.0000 Ova= 0.4336
CV(0): A( 29 )-> A( 35 ) c_i: 0.9969 Per: 99.4% IPA: 7.064 eV Oai: 0.4335
将S1态的单点能与S0态的单点能相减，即得momap需要的两态的能量差参数Ead=0.07502 a.u.。
基于基态结构，做S0和S1之间的非绝热耦合计算。点击 azulene_s0.hzw ，右击点击copy，设置new file name为nacme，在Project中出现 nacme.hzw 。双击 nacme.hzw ，使用BDF进行azulene的非绝热耦合计算。选中Simulator → BDF → BDF，界面中设置参数。在Basic Settings界面中的Calculation Type选择TDDFT-NAC，
方法采用默认的GB3LYP泛函，基组在Basis中的All Electron类型中，选择6-31G(d,p)。SCF Settings和TDDFT Settings面板中将Use MPEC+COSX Acceleraton的默认勾选去掉。在TDDFT Settings面板中的Non-Adiabatic Coupling内容框中，在默认的Coupling between Ground and Excited-State下，点击”+”号，Basic Settings、
SCF Settings、TDDFT Settings界面中的其它参数以及OPT Settings、Freq Settings等面板的参数使用推荐的默认值，不需要做修改。之后点击 Generate files 即可生成对应计算的输入文件。
选中生成的 bdf.inp 文件，右击选择open with，打开该文件，如下所示：
$compass
Title
C10H8
Geometry
C 0.79273796 0.49102542 -0.00003306
C -0.70229648 0.61186591 0.00000000
C 1.30022931 1.80163336 -0.00006271
C -0.99262499 1.98726799 0.00007812
C 1.54415131 -0.67887417 0.00000000
C -1.63318173 -0.42094562 -0.00002837
C 0.21877157 2.69859812 0.00000000
C 1.10656346 -2.00562675 0.00005788
C -1.41619168 -1.80093044 -0.00004813
C -0.20258112 -2.49333482 0.00000000
H 2.35092512 2.06249888 -0.00009827
H -1.98777599 2.41348149 0.00017649
H 2.62424717 -0.53731744 0.00001117
H -2.67585843 -0.10561277 -0.00001520
H 0.30641472 3.77916915 0.00002385
H 1.88966565 -2.75951312 0.00017580
H -2.31053950 -2.41870504 -0.00009018
H -0.29054446 -3.57807510 0.00000000
End Geometry
Basis
6-31G(d,p)
Skeleton
Group
C(1)
$end
$xuanyuan
Direct
$end
$scf
RKS
Charge
0
SpinMulti
1
DFT
GB3LYP
D3
MPEC+COSX
Molden
$end
$tddft
Imethod
1
Isf
0
Idiag
1
Iroot
6
Istore
1
$end
$resp
Quad
FNAC
Norder
1
Method
2
Nfiles
1
Single
States
1
1 1 1
$end
选中 bdf.inp 文件，右击选择Run提交作业，任务结束后 bdf.out ， bdf.scf.molden 三个结果文件会出现在Project中。
至此，MOMAP对azulene的S1→S0的辐射速率和内转换速率的计算需要的BDF量化软件的结构优化频率结果文件、非绝热耦合结果文件和参数部分都已完成。
MOMAP计算部分
在使用量化软件BDF完成azulene的基态和激发态的结构优化、频率计算以及非绝热耦合的计算，并计算完momap输入文件中需要的参数后，接下来将使用MOMAP软件对azulene的S1→S0的辐射速率和内转换速率进行计算，通过对比辐射速率和内转换速率来解释azulene的S1→S0的荧光难以被观测到的原因。
首先计算S1→S0的荧光辐射速率，第一步为做电子振动耦合（electron-vibration coupling, EVC）计算，该计算基于量化计算输出的分子振动频率、力常数矩阵，同时在内坐标以及直角坐标系下，计算分子跃迁发生初末态间的模式位移、黄昆因子、重整能以及Duschinsky转动矩阵。
将bdf的S0优化频率计算结果改为 zulene-s0.out ，将S1的计算结果改为 azulene-s1.out ，同时放在EVC计算文件夹中。
EVC计算的输入文件 momap.inp 为：
do_evc = 1
&evc
ffreq(1) = "azulene-s0.out"
ffreq(2) = "azulene-s1.out"
proj_reorg = .t.
/
提交脚本文件 momap.slurm ，任务运行结束后，生成如下文件
其中 evc.cart.dat 为利用直角坐标系计算得到的模式位移、黄昆因子、重整能以及 Duschinsky 转动矩阵的结果；而 evc.dint.dat 为内坐标计算模式位移、黄昆因子、重整能，直角坐标计算Duschinsky转动矩阵的结果。
打开 evc.cart.dat 文件，查看总重整能的数值：
--------------------------------------------------------------------------------------------------------------------------------------------
Total reorganization energy (cm-1): 4032.869339 5126.265767
--------------------------------------------------------------------------------------------------------------------------------------------
并与 evc.dint.dat 文件中的数值进行比较：
--------------------------------------------------------------------------------------------------------------------------------------------
Total reorganization energy (cm-1): 4070.407661 5114.173064
--------------------------------------------------------------------------------------------------------------------------------------------
比较 evc.cart.dat 以及 evc.dint.dat 文件内的重整能，若重整能相差不大（< \(1000 cm^{-1}\) ），使用 evc.cart.dat 文件进行后续计算，若重整能相差较大（一般情况 evc.cart.dat 大于 evc.dint.dat ），使用 evc.dint.dat 文件进行后续计算。这里两者较为接近，且数值较为合理（< \(10000 cm^{-1}\) ），可使用 evc.cart.dat 进行下一步S1→S0的荧光辐射速率计算。
此外，还可以根据EVC计算的结果文件做更多的后处理。
evc.dx.x.com 和 evc.dx.x.xyz 为两个电子态分子叠加图，其中 evc.dx.x.com 可用GaussView打开，在View-Display Format-Molecule中选择Tube类型，显示如下：
evc.dx.x.xyz 可用Jmol软件打开，显示如下：
evc.dx.v.xyz 也是分子态叠加图，用Jmol软件打开，选择显示-矢量-0.1A，显示如下：
evc.cart.abs 为Duschinsky矩阵文件，可用来画Duschinsky矩阵二维图。可以在Device Studio中，选择Simulator-momap-analysis，打开 evc.cart.abs 文件，显示如下：
同样的方法打开 evc.dint.abs 文件，可在ColorType下拉选框选择显示颜色，显示如下：
Duschinsky矩阵都是用直角坐标计算的，二者选一即可。
将 evc.vib1.xyz 和 evc.vib2.xyz 文件与 evc.cart.dat 文件放在同一路径下，在Device Studio中，选择Simulator-momap-analysis，打开 evc.cart.dat 文件，出现基态和激发态下各个振动模式对应的重整能和黄昆因子图，显示如下：
点击Choose Color可以改变线的颜色，改变Set Width中的数值可改变线的粗细。点击图中的线，右侧结构将展示相应的振动模式。振动的快慢可通过Animation Frequency调整，振动的幅度可根据Displacement Amplitude显示。
接下来进行S1→S0的荧光光谱及荧光辐射速率的计算。此部分计算需要EVC计算所得的 evc.*.dat 文件作为输入。如前所述，将 evc.cart.dat 放在荧光光谱及荧光辐射速率的计算文件中，输入文件 momap.inp 为:
do_spec_tvcf_ft = 1
do_spec_tvcf_spec = 1
&spec_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.07502 au
EDMA = 0.9267 debye
EDME = 0.6507 debye
FreqScale = 0.70
DSFile = "evc.cart.dat"
Emax = 0.3 au
dE = 0.00001 au
logFile = "spec.tvcf.log"
FtFile = "spec.tvcf.ft.dat"
FoFile = "spec.tvcf.fo.dat"
FoSFile = "spec.tvcf.spec.dat"
/
提交脚本文件 momap.slurm ，任务计算结束后，生成文件如下所示：
spec.tvcf.ft.dat为关联函数文件，内容如下：
#F(t) file
#time(fs) abs_FC_Re abs_FC_Im emi_FC_Re emi_FC_Im
-1000.00000 -3.222659695E-06 -3.680115014E-05 -3.738401821E-06 1.830624064E-05
-999.00000 -1.314545717E-05 -4.054982998E-05 -1.645610065E-07 1.674654956E-05
-998.00000 -2.552186539E-05 -3.557857053E-05 1.690265080E-06 1.428576477E-05
-997.00000 -3.030022340E-05 -2.305671104E-05 2.187359504E-06 1.222900063E-05
-996.00000 -2.620263828E-05 -1.307399386E-05 2.137378198E-06 1.094934359E-05
-995.00000 -1.984646799E-05 -8.898814274E-06 2.126481370E-06 1.032271509E-05
-994.00000 -1.480939587E-05 -8.192889658E-06 2.544875562E-06 9.992907379E-06
-993.00000 -1.150600369E-05 -8.808752896E-06 3.437467322E-06 9.316576838E-06
-992.00000 -9.618163543E-06 -9.696518247E-06 4.080694570E-06 7.874673573E-06
-991.00000 -8.767387663E-06 -1.012196928E-05 3.827322982E-06 6.315202077E-06
-990.00000 -8.083181010E-06 -9.635894264E-06 3.049519176E-06 5.302908607E-06
-989.00000 -6.767527992E-06 -8.787573767E-06 2.230010837E-06 4.817300527E-06
-988.00000 -4.968246898E-06 -8.254624823E-06 1.540853782E-06 4.667945882E-06
-987.00000 -3.046038369E-06 -8.166135596E-06 1.028297048E-06 4.730444119E-06
-986.00000 -1.166755532E-06 -8.469415398E-06 7.338287121E-07 4.901384390E-06
计算完成后先确认关联函数是否收敛，将 spec.tvcf.ft.dat 拖入origin中，选择第一列和第二列作图：
随时间趋于0表示吸收光谱计算收敛。选择第一列和第四列作图：
随时间趋于0表示发射光谱计算收敛。
spec.tvcf.spec.dat 为光谱文件：
#spectrum
#1Energy(Hartree) 2Energy(eV) 3WaveNumber(cm-1) 4WaveLength(nm) 5FC_abs 6FC_emi 7FC_emi_intensity
1.63970997E-05 4.46187977E-04 3.59874741E+00 2.77874462E+06 2.22315237E-05 1.35978424E-12 5.71912365E-20
9.23885920E-05 2.51402258E-03 2.02769521E+01 4.93170766E+05 1.24781305E-04 2.44296355E-10 5.78932299E-17
1.68380084E-04 4.58185718E-03 3.69551568E+01 2.70598229E+05 2.26516688E-04 1.48614131E-09 6.41864375E-16
2.44371577E-04 6.64969178E-03 5.36333615E+01 1.86451114E+05 3.27482697E-04 4.56299919E-09 2.86018106E-15
3.20363069E-04 8.71752639E-03 7.03115662E+01 1.42224111E+05 4.27602118E-04 1.03315707E-08 8.48987397E-15
3.96354561E-04 1.07853610E-02 8.69897710E+01 1.14956045E+05 5.27021609E-04 1.96507539E-08 1.99781616E-14
4.72346054E-04 1.28531956E-02 1.03667976E+02 9.64618045E+04 6.25550501E-04 3.34234748E-08 4.04952733E-14
5.48337546E-04 1.49210302E-02 1.20346180E+02 8.30936218E+04 7.23426877E-04 5.25150757E-08 7.38625716E-14
6.24329038E-04 1.69888648E-02 1.37024385E+02 7.29797108E+04 8.20443555E-04 7.78979043E-08 1.24747472E-13
7.00320530E-04 1.90566994E-02 1.53702590E+02 6.50607125E+04 9.16771962E-04 1.10429551E-07 1.98369368E-13
7.76312023E-04 2.11245340E-02 1.70380794E+02 5.86920611E+04 1.01226209E-03 1.51170543E-07 3.01020401E-13
8.52303515E-04 2.31923686E-02 1.87058999E+02 5.34590693E+04 1.10713559E-03 2.00942496E-07 4.39297283E-13
9.28295007E-04 2.52602032E-02 2.03737204E+02 4.90828371E+04 1.20107277E-03 2.60943515E-07 6.21333781E-13
1.00428650E-03 2.73280378E-02 2.20415409E+02 4.53688790E+04 1.29454751E-03 3.31900017E-07 8.54982733E-13
1.08027799E-03 2.93958724E-02 2.37093613E+02 4.21774330E+04 1.38703730E-03 4.15196637E-07 1.15048722E-12
将 spec.tvcf.spec.dat 拖入origin中，选择第三列和第五列作图，得到吸收光谱：
选择第三列和第六列作图，得到发射光谱：
打开spec.tvcf.log，文件末尾输出了荧光辐射速率值，
radiative rate (0): 7.21227543E-12 2.98165371E+05 /s, 3353.84 ns
荧光辐射速率在第一个数和第二个数读取，单位分别为a.u.和 \(s^{-1}\) ，第三个数为寿命，单位为ns。这里，azulene的S1→S0的荧光辐射速率为 2.98165371E+05 /s，荧光寿命为3353.84 ns。
接下来计算azulene的内转换速率。第一步为EVC振动分析计算。
内转换速率的EVC振动分析计算需要非绝热耦合计算结果文件。将非绝热耦合计算结果文件改名为 azulene-nacme.out ，与 azulene-s0.out 和 azulene-s1.out 放在内转换速率计算文件夹下，输入文件 momap.inp 为：
do_evc = 1
&evc
ffreq(1) = "azulene-s0.log"
ffreq(2) = "azulene-s1.log"
fnacme = "azulene-nacme.log"
proj_reorg = .t.
/
提交脚本文件 momap.slurm ，任务计算结束后，生成如下文件：
与前述荧光辐射速率的evc计算相比，多出一个 evc.cart.nac 文件，该文件将和 evc.cart.dat 文件一起用于接下来的内转换速率的计算。
内转换速率输入文件 momap.inp 为：
do_ic_tvcf_ft = 1
do_ic_tvcf_spec = 1
&ic_tvcf
DUSHIN = .t.
Temp = 300 K
tmax = 1000 fs
dt = 1 fs
Ead = 0.07502 au
FreqScale = 0.40
DSFile = "evc.cart.dat"
CoulFile = "evc.cart.nac"
Emax = 0.3 au
logFile = "ic.tvcf.log"
FtFile = "ic.tvcf.ft.dat"
FoFile = "ic.tvcf.fo.dat"
/
计算结束后，生成文件如下：
其中 ic.tvcf.ft.dat 为关联函数文件，内容如下：
#time(fs) IC_ft_Re(au) IC_ft_Im(au)
-1000.00000 -0.302000787E-13 0.332374045E-13
-999.00000 -0.205081112E-13 0.364896559E-13
-998.00000 -0.114899559E-13 0.364092975E-13
-997.00000 -0.403418594E-14 0.342290784E-13
-996.00000 0.187447714E-14 0.309462504E-13
-995.00000 0.652215564E-14 0.270373111E-13
-994.00000 0.100875106E-13 0.226564272E-13
-993.00000 0.125729252E-13 0.178887681E-13
-992.00000 0.138917246E-13 0.128953745E-13
-991.00000 0.139788795E-13 0.793558174E-14
-990.00000 0.128516224E-13 0.332590373E-14
-989.00000 0.106275531E-13 -0.595694484E-15
-988.00000 0.755502051E-14 -0.349394006E-14
-987.00000 0.406056718E-14 -0.511604680E-14
-986.00000 0.703538619E-15 -0.544879345E-14
-985.00000 -0.203569213E-14 -0.479493437E-14
计算结束后首先确认关联函数是否收敛，将 ic.tvcf.ft.dat 拖入origin中，选择第一列和第二列作图：
随时间趋于0表示关联函数收敛。
其中 ic.tvcf.fo.dat 为谱函数文件，内容如下：
#1Energy(Hartree) 2Energy(eV) 3WaveNumber(cm-1) 4WaveLength(nm) 5radi-spectrum 6kic(s^{-1}) 7log(kic)
0.00000000E+00 0.00000000E+00 0.00000000E+00 Infinity 0.42121076E-05 0.17413431E+12 0.11240884E+02
0.75991492E-04 0.20678346E-02 0.16678205E+02 0.59958492E+06 0.44127214E-05 0.18242796E+12 0.11261091E+02
0.15198298E-03 0.41356692E-02 0.33356409E+02 0.29979246E+06 0.46209108E-05 0.19103480E+12 0.11281112E+02
0.22797448E-03 0.62035038E-02 0.50034614E+02 0.19986164E+06 0.48375050E-05 0.19998910E+12 0.11301006E+02
0.30396597E-03 0.82713384E-02 0.66712819E+02 0.14989623E+06 0.50633471E-05 0.20932572E+12 0.11320823E+02
0.37995746E-03 0.10339173E-01 0.83391024E+02 0.11991698E+06 0.52978153E-05 0.21901896E+12 0.11340482E+02
0.45594895E-03 0.12407008E-01 0.10006923E+03 0.99930820E+05 0.55357841E-05 0.22885692E+12 0.11359564E+02
0.53194045E-03 0.14474842E-01 0.11674743E+03 0.85654988E+05 0.57737058E-05 0.23869293E+12 0.11377840E+02
0.60793194E-03 0.16542677E-01 0.13342564E+03 0.74948115E+05 0.60096238E-05 0.24844610E+12 0.11395232E+02
为检查是否满足能隙定律，将 ic.tvcf.fo.dat 拖入origin中，选择第三列和第七列作图：
此外， ic.tvcf.fo.dat 文件中第一列和第六列表示不同Ead下的非辐射速率。
在 ic.tvcf.log 文件的末尾，读取S1→S0的非辐射速率值，
#1Energy(Hartree) 2Energy(eV) 3WaveNumber(cm-1) 4WaveLength(nm) 5radi-spectrum 6kic(s^{-1}) 7log(kic) 8time(ps)
7.50036029E-02 2.04095275E+00 1.64613880E+04 6.07482186E+02 6.54396018E-07 2.70536301E+10 1.04322255E+01 36.96361624
这里S1→S0的非辐射速率为2.70536301E+10 \(s^{-1}\) 。
对比azulene的荧光辐射速率和非辐射速率，荧光辐射速率为2.98165371E+05 /s，非辐射速率为2.70536301E+10 \(s^{-1}\) ，非辐射速率比荧光辐射速率高五个数量级，从而降低了azulene的S1→S0的荧光量子效率，所以S1→S0的荧光难以被观测到，从而表现出反kasha规则的荧光现象。
用MSSM模型研究自旋禁阻的多态反应
在大多数化学反应过程中，体系从反应物开始，经过一个或多个过渡态以及中间体，最后到达产物，整个反应都发生在基态的势能面上。 然而在有些情况下，反应物、过渡态、中间体、产物可能处于不同自旋态的势能面上，称为 自旋禁阻反应 （Spin-Forbidden Reaction）或 自旋禁阻的多态反应，以下简称为 多态反应 （Multiple-State Reaction）。这类反应的发生是由自旋轨道耦合导致的，进行理论研究具有一定的挑战性。
备注
大多数多态反应仅涉及两个自旋态，故在一些文献中也称为 两态反应 。
自旋守恒的化学反应也可能发生在自旋相同的两个（或多个）势能面上，但是如果不优化IRC，很难察觉。 这种情况用常规理论方法即可处理，故不属于上面所说的多态反应。
在化学反应中严格考虑自旋轨道耦合，需要使用二分量或四分量相对论方法（特别是二分量和四分量相对论密度泛函理论， 已经在一些量子化学程序中实现了解析梯度），但是由于使用的困难，这类方法很少得到实际应用。 如果自旋轨道耦合不强，可以用最小能量交点方法近似处理，即优化两种自旋态势能面的最低能量交点（参见 MECP ）。 尽管MECP方法得到广泛应用，但是也存在一些不足：(1) 如果MECP对应新的过渡态，由于MECP并不是真正的驻点，无法获该点的振动频率和热化学量；(2) MECP 无法考虑两个以上自旋态的同时交叉，或在小区域内的相互交叉；(3) 如果MECP导致新过渡态，即便是不太重的3d元素， MECP的位置和能量有时候也会和二分量相对论结果相差甚远 [93] ；(4) MECP不适用于锕系这类重元素体系；对于成键不饱和的5d体系，适用性存疑。
2018年，Truhlar等人提出了两态自旋混合模型（TSSM），用于模拟化学反应过程中两个自旋态之间的旋轨耦合 [80] 。 最近，我们将其推广为多态自旋混合模型（MSSM），适用于多个自旋态的情况，不仅避免了MECP方法的前三个缺点， 并且对于5d以内的元素都是适用的 [94] 。 对5d化合物反应的测试结果表明，MSSM-DFT与二分量DFT计算的反应（吸）放热、过渡态能垒高度相差在3 kcal/mol以内 （对于成键非饱和的5d化合物，需要用二分量或四分量方法做单点能校正）。
作为示例，以下我们用MSSM模型重复Truhlar等人用TSSM模型研究的钨络合物脱氢过程的两态反应 \(W(CH_3)(n-C_3H_7)(OH)_2 \rightarrow WH(CH_3)(=CHC_2H_5)(OH)_2\) [80] 。 在文献中，作者用该体系做为原子层沉积多相钨催化剂的近似，其中载体骨架被两个OH基团取代。 反应的示意图如下所示，其中蓝线是单重态反应路径，红线是三重态反应路径，黑线是两态混合后的最低能量路径，黄线是单重态在最低混合态中所占的比例。 在整个反应中，钨的氧化态中从+4变为+6。由于在反应物中钨是部分配位的，因此三重态能量最低，自旋混合基态以三重态为主； 而在产物中的钨成键饱和，因此闭壳层单重态的能量最低，自旋混合基态以单重态为主。 三重态、单重态交点与自旋混合基态的过渡态仍有一定的距离（虽然能量比较接近），其中单重态在后者中占了90%左右， 也就是说，两种自旋态在过渡态并不是按照 1:1 的比例混合的，这是MSSM与MECP结果的最大区别。
一般来说，用MSSM模型进行计算之前，首先要用常规的标量方法对各种可能的自旋态优化反应物，产物，过渡态，和中间体结构， 从中找到对反应起重要作用的自旋态以及MSSM优化计算的初始结构。 由于该反应已经在文献中进行了细致的研究 [80] ，我们跳过了这一步， 直接用MSSM优化这个两态反应在自旋混合基态势能面上的反应物，产物，和过渡态结构，并计算这三个驻点结构的振动频率。 初始分子坐标直接取自文献。 W的基组用def2-TZVP，其它原子的基组为6-31G*。SOC经验参数采用文献中的302 meV（2436 \(\rm cm^{-1}\)）。 唯一不同的是，我们把文献中采用的M06L泛函改为B3LYP泛函，因为后者具有更好的数值稳定性。计算结果列于下表， 并附上文献中的TSSM结果 [80] 、二分量相对论结果 [94] 作为对照：
方法 |
过渡态 |
产物 |
虚频 |
|---|---|---|---|
MSSM-B3LYP |
25.6 |
3.2 |
924i |
TSSM-M06L |
28.3 |
3.2 |
962i |
2c-M06L |
25.8 |
0.0 |
更多的应用案例可参考Takayanagi课题组的系列论文。 他们用TSSM模型研究了大量涉及3d、4d元素的两态反应（见文献 [87] 以及引用该文献的其它论文）。
BDF计算输入
优化自旋混合基态反应物的输入如下，其中一些复制命令的含义参见 ZnS分子自旋混合态算例 的注释。
$compass
title
W system: reactant
basis-block
6-31G*
W=def2-tzvp
end basis
# geom from PCCP 2018, 20, 4129
geometry
C -2.50599000 0.80875000 0.27831700
H -2.42978000 0.87169600 1.37614800
H -3.29329400 1.52978300 0.00369300
C -1.17932200 1.22962500 -0.35262900
H -1.30089800 1.27432600 -1.45149600
W 0.49596700 -0.04848100 0.05050100
H -0.94880200 2.26591700 -0.04284700
C -2.96148400 -0.58790000 -0.11445500
H -2.28153800 -1.36238600 0.26893400
H -3.96244500 -0.81403100 0.27115200
H -2.98974600 -0.70153100 -1.20640700
O 1.64397700 1.44403300 -0.35823300
H 1.22595200 2.31740600 -0.43846900
O 0.49410200 -0.94895200 1.73966100
H -0.19914300 -0.83307900 2.40480100
C 0.26970600 -1.43014300 -1.54486300
H -0.59187900 -2.09433400 -1.38019500
H 0.07131900 -0.88842400 -2.48170400
H 1.15657400 -2.06035000 -1.71035900
end geometry
nosymm
$end
$bdfopt
multistate
2soc 2436
solver
1
hess
final
$end
$xuanyuan
$end
%cp $BDF_WORKDIR/$BDFTASK.scforb.1 $BDF_WORKDIR/$BDFTASK.scforb 2>/dev/null || :
$scf
rks
dft
b3lyp
grid
fine
$end
%cp $BDF_WORKDIR/$BDFTASK.scforb $BDF_WORKDIR/$BDFTASK.scforb.1
$resp
geom
$end
%cp $BDF_WORKDIR/$BDFTASK.egrad1 $BDF_WORKDIR/$BDFTASK.egrad.1 2>/dev/null || :
%cp $BDF_WORKDIR/$BDFTASK.hess $BDF_WORKDIR/$BDFTASK.hess.1 2>/dev/null || :
%cp $BDF_WORKDIR/$BDFTASK.scforb.2 $BDF_WORKDIR/$BDFTASK.scforb 2>/dev/null || :
$scf
uks
dft
b3lyp
spinmulti
3
grid
fine
$end
%cp $BDF_WORKDIR/$BDFTASK.scforb $BDF_WORKDIR/$BDFTASK.scforb.2
$resp
geom
$end
%cp $BDF_WORKDIR/$BDFTASK.egrad1 $BDF_WORKDIR/$BDFTASK.egrad.2 2>/dev/null || :
%cp $BDF_WORKDIR/$BDFTASK.hess $BDF_WORKDIR/$BDFTASK.hess.2 2>/dev/null || :
优化自旋混合基态产物的输入类似，以下只给出 compass 模块的输入，其余输入同上。
$compass
title
W system: product
basis-block
6-31G*
W=def2-tzvp
end basis
# geom from PCCP 2018, 20, 4129
geometry
W -0.42387000 -0.03021600 0.04243900
O -0.54695900 -1.59597800 -1.06031400
H -1.38476400 -2.05907900 -0.89227700
O -2.07683400 -0.36779800 1.01946300
H -2.09920000 -0.45365300 1.97932500
C -1.19429700 1.74966400 -0.81144000
H -1.16507500 2.46066600 0.03099700
H -0.64596000 2.20075800 -1.64296300
H -2.24817300 1.61217900 -1.08859500
H 0.10061200 0.01381000 1.68243000
C 1.39659000 0.37442700 -0.02642700
H 1.22850000 0.24599300 -1.14781000
C 2.78803600 0.62890000 0.40356700
H 2.79144700 0.71717400 1.49877400
H 3.12713400 1.60782200 0.02677600
C 3.77244700 -0.45270200 -0.03426000
H 4.78897300 -0.22737400 0.30879300
H 3.48486900 -1.42988500 0.36906600
H 3.80168700 -0.54392400 -1.12680400
end geometry
nosymm
$end
优化自旋混合基态过渡态只需要修改反应物输入中 compass 、 bdfopt 两个模块的参数，如下：
$compass
title
W system: TS
basis-block
6-31G*
W=def2-tzvp
end basis
# geom from PCCP 2018, 20, 4129
geometry
C 2.75476500 -0.49692100 0.66265200
H 2.66846800 -0.03811600 1.65776400
H 3.19563100 -1.49505600 0.81976100
C 1.39220000 -0.59199800 0.03929900
H 1.40125200 -1.14322200 -0.93707500
W -0.48542700 -0.04574900 -0.00978200
H 0.62772500 -0.92650400 1.07995600
C 3.68608400 0.33536600 -0.20966500
H 3.29355400 1.35035000 -0.34119700
H 4.68482600 0.40829800 0.23538400
H 3.79640600 -0.10658300 -1.20747300
O -1.41339700 -1.47110300 -0.89277200
H -1.32731100 -2.36616000 -0.52320400
O -1.79388100 0.67727500 1.22320500
H -1.75031500 0.44148600 2.15765700
C -0.22357400 1.89123400 -0.83239900
H 0.11804400 2.54791300 -0.01489800
H 0.44852800 2.01363400 -1.68618700
H -1.23385200 2.22393900 -1.11944000
end geometry
nosymm
$end
$bdfopt
multistate
2soc 2436
solver
1
hess
init+final
iopt
10
$end