电子激发任务中轨道跃迁贡献的计算
电子激发任务中轨道跃迁贡献的计算
Calculation of orbital transition contributions in electron excitations
文/Sobereva @北京科音
First release: 2014-Apr-15 2020-Jan-28
做ZINDO、CIS、TDHF、TDDFT这类电子激发计算时，电子激发通常表示为大量分子轨道(MO)间的跃迁，Gaussian会输出对应的跃迁系数，比如Excited State 1: Singlet-A" 4.6966 eV 263.99 nm f=0.0003
14 -> 16 0.56407
14 -> 20 -0.21920
14 -> 25 -0.24192
14 -> 29 -0.17712
14 -> 33 -0.13568
这些数值确切来说是轨道跃迁对应的单激发组态函数的组态系数。对于CIS、TDHF、TDDFT、ZINDO方法，激发态波函数都是通过各种单激发组态函数线性组合来描述的。如果想了解细节，大家应当找 书看看关于CIS方法的介绍，就自然不难理解其物理意义了。为了分析电子跃迁的本质特征，经常会考察每对MO跃迁对电子激发的贡献。其实这是个极为简单的问题，但是总看到一些人错误的讨论，故在这里澄清一下计算方法。这里只是讨论Gaussian程序的情况，对其它程序讨论是类似的。
如果对Gaussian做TDDFT计算不了解，务必参看《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（//www.umsyar.com/314）。
1 ZINDO和CIS计算
对于基态（确切来说是参考态）是闭壳层的计算，所有跃迁系数平方和相加等于0.5，即50%。之所以加和不为理应的1.0，是因为alpha轨道和beta轨道此时完全一致，alpha电子的跃迁和beta电子的跃迁的情况因此也是完全一致的，因此只输出了alpha部分。例如上例的14号MO到16号MO的跃迁（14->16），其实分为两种情况，即14A->16A和14B->16B，由于二者的跃迁系数是完全相同的，故只输出了其一。因此，14->16对这个激发的贡献就是0.56407^2/0.5*100%=63.6%，14->29的贡献就是(-0.17712)^2/0.5*100%=6.3%。然而，上面的例子的5个跃迁系数的平方和相加结果为0.47453，离0.5还差得远，这显然不是有效位数的问题。于是有些人在计算贡献时，比如计算14->16的贡献时，居然这么算：0.56407^2/0.47453*100%=67.1%，这是明显错误的！之所以加和不为0.5，是因为为了输出信息简洁，Gaussian默认只把系数绝对值大于0.1的系数输出了出来。如果想让更多的系数都输出出来，可以用IOp(9/40=x)关键词，它说明把绝对值大于10^-x的都输出出来。如果你用比如IOp(9/40=3)把绝对值大于0.001的系数都输出，那么平方和就非常接近于0.5了，但是此时屏幕上的信息量会略大，可能会有好几十、几百个系数。
对于基态是开壳层的计算，alpha轨道和beta轨道不再一致，故alpha和beta电子激发的情况在计算中就都明确考虑了，此时所有跃迁系数平方和为1.0。因此对于一对儿MO的跃迁，将其系数平方乘上100%就得到了它的贡献了。
老有人问某些组态系数为负的时候是什么意思。根本无需关注正负号，算贡献时求了平方之后都一样，而且正负号也没有可以直观去理解的物理意义（分子轨道展开系数里有正有负还比较容易通过绘制轨道图形来说明，然而组态函数本身是多电子基函数，其系数的正负号是相当抽象的，既没法靠画图展现，也没法给予明确的物理解释）。
2 TDHF和TDDFT计算
TDHF和TDDFT计算和上一节的情况最大的不同在于不仅有激发，即“->”所表示的，还有去激发，即“<-”所表示的。“去激发”并没有物理意义，这种称呼仅仅是根据TDHF、TDDFT波函数的归一化特点而给了这么个称呼而已，不要去深究和过度解释。这种TD形式的计算出现去激发是必然的，但通常去激发所占成份很小。例如Excited State 3: Singlet-B2 4.7665 eV 260.12 nm f=0.0044
35 -> 40 0.02865
36 -> 39 0.27686
36 -> 44 0.06122
36 -> 57 -0.01098
36 -> 62 -0.01239
37 -> 40 0.63762
38 -> 39 -0.10663
38 -> 50 -0.01114
36 <- 39 0.02135
36 <- 44 0.01165
37 <- 40 -0.01239
对于TDHF/TDDFT计算，对于基态为闭壳层和开壳层时还是分别归一化为0.5和1.0，但注意，所有激发的系数平方和减去所有去激发的系数平方和得到的才是这个结果，也就是说去激发起到的是负贡献。比如上例37->40的贡献为(0.63762)^2/0.5*100%=81.31%，而36<-44的贡献为-(0.01165)^2/0.5*100%=-0.027%，37<-40的贡献为-(-0.01239)^2/0.5*100%=-0.031%。
值得一提的是，有一种TDA(Tamm-Dancoff approximation)形式的TDDFT，也叫TDA-DFT，它和一般的TDDFT的关系形同于CIS与TDHF的关系，也不存在去激发，这在Gaussian D.01及之后版本中可以用TDA关键词来使用。TDA近似可以令TDDFT计算有一些加快，对激发能和振子强度有一定影响，精度往往比TDDFT还好一些，但是振子强度比TDDFT整体低估，算的转子强度也不怎么好。更多信息看《乱谈激发态的计算方法》（ //www.umsyar.com/265）。
3 使用Multiwfn输出贡献值
如上所示，计算MO跃迁的贡献很简单，用计算器点几下就出来了。如果懒得手算，也可以用著名的波函数分析程序Multiwfn来输出贡献值，Multiwfn中有一大批功能专门用来做电子激发分析。Multiwfn可以在其主页//www.umsyar.com/multiwfn免费下载。如果不熟悉Multiwfn建议先参看《Multiwfn入门tips》（//www.umsyar.com/167）。下文对应的是3.6版及之后版本，更老版本与下面的操作过程不同。
对于ZINDO/CIS/TDHF/TDDFT任务，计算之后将chk文件转化为fch文件，假设叫a.fch，相应的Gaussian输出文件假设叫a.out，都放在了C:\下。启动Multiwfn，依次输入
C:\a.fch
18 //电子激发分析模块
-1 //检查、修改和导出组态系数
C:\a.out
3 //假设分析的是第3个电子激发态
然后屏幕上输出了一些统计信息，包括激发系数的平方和以及去激发系数平方和的负值。之后马上看到以下结果，是对电子激发贡献最大的10个轨道跃迁信息。#后面那个数是这个轨道跃迁在所有输出的轨道跃迁中的序号，不用管。
Some MO transitions sorted by absolute contributions:
# 1272 48 -> 49 Coeff.: 0.70422 Contri.: 99.1852%
# 2489 48 <- 49 Coeff.: -0.07085 Contri.: -1.0039%
# 1202 46 -> 49 Coeff.: -0.05577 Contri.: 0.6221%
# 748 33 -> 49 Coeff.: -0.03030 Contri.: 0.1836%
# 1274 48 -> 53 Coeff.: 0.02771 Contri.: 0.1536%
# 1273 48 -> 50 Coeff.: 0.02375 Contri.: 0.1128%
# 1203 46 -> 50 Coeff.: -0.02166 Contri.: 0.0938%
# 1224 47 -> 68 Coeff.: 0.02038 Contri.: 0.0831%
# 1183 45 -> 51 Coeff.: 0.01880 Contri.: 0.0707%
# 1218 47 -> 54 Coeff.: -0.01552 Contri.: 0.0482%
如果想把数据从屏幕上拷贝出来，见Multiwfn手册5.4节的说明。如果只想把贡献超过某个阈值的轨道跃迁输出出来，可以在当前界面里用选项-2，然后输入阈值。
使用Multiwfn还可以以更直观的方式非常方便快捷地对一大批激发态输出主要轨道跃迁贡献，而且只需要提供Gaussian输出文件就够了而不需要fch，输出信息例如
# 1 3.907 eV Multi= 1: H-4 -> L 81.9%, H-4 -> L+2 12.1%
# 2 4.062 eV Multi= 1: H -> L 86.0%, H-3 -> L 5.3%
# 3 4.417 eV Multi= 1: H-6 -> L 85.3%, H-6 -> L+2 11.9%
# 4 4.791 eV Multi= 1: H-2 -> L 54.5%, H -> L+1 27.6%, H-3 -> L+1 6.4%
# 5 4.887 eV Multi= 1: H -> L+3 57.3%, H-2 -> L 17.0%, H-1 -> L+2 8.8%, H
-1 -> L 8.0%
操作见《使用Multiwfn便利地查看所有激发态中的主要轨道跃迁贡献》（//www.umsyar.com/529）。