Anti-Kasha的Excited State Dynamics计算流程分享

wal

Anti-Kasha荧光是违反Kasha规则的发光现象，即分子从高激发态（如S₂、S₃）直接发射荧光，而非仅从最低激发态S₁发射。此类现象在有机染料中较为罕见，但在某些特殊结构的分子中可以观察到。本文将以Azulene为例尝试对此效应进行计算。

1 理论背景

1.1 Kasha规则与Anti-Kasha现象

Kasha规则指出，分子在光激发后，无论最初被激发到哪个高能态，最终的荧光发射通常仅来源于最低激发态S₁。这是因为高激发态通常比较密集，态之间的能量差相对较小。根据energy gap law，内转换速率随能隙变大呈指数下降，而较小的能隙导致了较大的kic，导致Sn→Sn-1的内转换速率远大于Sn的辐射速率。

然而，在某些特殊情况下，分子可以违反Kasha规则，直接从S₂或更高激发态发射荧光，这被称为Anti-Kasha发射。有利于发生Anti-Kasha发射的条件：

1. 大的S₂-S₁能隙：当S₂与S₁之间的能量差足够大时，内转换速率显著降低
2. 弱的S₂-S₁振动电子耦合：非绝热耦合较弱，阻碍内转换过程
3. 激发性质差异：S₂与S₁具有不同的激发特征（如pi-pi*与n-pi*），导致跃迁禁阻。(硫光气)

1.2 理论描述

Anti-Kasha发射的量子产率可以用竞争动力学来描述：

其中：
- kr,S2：S₂态的辐射跃迁速率
- kIC,S2：S₂态到基态的内转换速率  
- kIC,S2→S1：S₂→S₁内转换速率
注：此处忽略了系间窜跃。对于系间窜跃速率不可忽视的体系，应该同时考虑在内。

当kr,S2与其余两项之和的数量级相当或更大时，才能观察到明显的Anti-Kasha发射。

2 计算流程

Anti-Kasha荧光的理论计算内容分为这几块：

- 几何优化：基态和激发态的平衡几何结构、正则振动模式、绝热激发能、跃迁偶极矩。拥有TDDFT解析hessian的g16胜任此计算。
- 非绝热耦合：态间的NAC矩阵元。激发态间的NAC目前只有Q-Chem和BDF两家胜任，由于笔者没有使用过Q-Chem，本文采用BDF。
- 激发态动力学参数：利用前述中间量计算激发态衰减速率。MOMAP拥有对BDF进行TDDFT-NAC输出的解析支持，但此程序较难获取；FCclasses3是免费开源程序，可以结合笔者的脚本进行计算，也是个不错的选择。

接下来，笔者将依次计算这些量。

3 Gaussian
平平无奇的三个opt：

1. # opt freq B3LYP/def2svp scrf em=gd3bj
   
2. # opt freq td B3LYP/def2svp scrf em=gd3bj
   
3. # opt freq td(root=2) B3LYP/def2svp scrf em=gd3bj

复制代码

以及基于opt的一个TD单点：

1. # td B3LYP/def2svp scrf em=gd3bj

复制代码

对于使用MOMAP计算ESD的情况，我们需要在输出中提取这些量：
- 基态能量：`SCF Done`后面的数值

1. SCF Done:  E(RB3LYP) =  -385.590878548     A.U. after    8 cycles
   
2.             NFock=  8  Conv=0.47D-08     -V/T= 2.0122
   
3. DoSCS=F DFT=T ScalE2(SS,OS)=  1.000000  1.000000

复制代码

- 每个激发态的能量：`E(TD-HF/TD-DFT)`后面的数值

1. Excitation energies and oscillator strengths:
   
2. 
   
3. Excited State   1:      Singlet-A      1.7066 eV  726.49 nm  f=0.0074  <S**2>=0.000
   
4.       34 -> 35        -0.70605
   
5. This state for optimization and/or second-order correction.
   
6. Total Energy, E(TD-HF/TD-DFT) =  -385.528161673   
   
7. Copying the excited state density for this state as the 1-particle RhoCI density.
   
8. 
   
9. Excited State   2:      Singlet-A      3.3980 eV  364.88 nm  f=0.0099  <S**2>=0.000
   
10.       33 -> 35        -0.46258
    
11.       34 -> 36         0.53273

复制代码

提取时应当注意到，在当前级别下S2/S1的gap足足有0.05534 a. u.(1.5 eV)。根据energy gap law，可以知道这两个态之间的内转换速率不会那么快，这为anti-kasha发射创造了条件。

- TDM：对应state的Dip. S.开根号再乘以au→debye转换因子。不仅需要提取S1结构的，还需要S0结构的。

1. Ground to excited state transition electric dipole moments (Au):
   
2.        state          X           Y           Z        Dip. S.      Osc.
   
3.          1        -0.0000      0.4217     -0.0000      0.1778      0.0074
   
4.          2         0.3447      0.0000      0.0001      0.1188      0.0099
   
5.          3         3.6573     -0.0000      0.0000     13.3755      1.3606

复制代码

而对于使用FCclasses的情况，可以不必手动提取这些信息，因为该程序自带的小工具可以自动从fchk文件中拿到比log中打印精度更高的数据。

4 BDF

BDF程序可以通过resp模块计算激发态之间的NAC，这是本程序亮点之一。

4.1 基态与激发态间NAC

计算基态与S₁态之间的NAC：

1. $resp
   
2. Quad
   
3. FNAC
   
4. Norder
   
5.   1
   
6. Method
   
7.   2  
   
8. Nfiles
   
9.   1
   
10. Single          # calculate NACME between excited state and ground state
    
11. States
    
12.   1
    
13.   1 1 1         # <istore number> <Irreducible representation order> <root order>
    
14. $end

复制代码

4.2 激发态间NAC

计算S₁与S₂态之间的NAC，这是判断Anti-Kasha发射可能性的关键：

1. $resp
   
2. Quad
   
3. FNAC
   
4. Norder
   
5.   1
   
6. Method
   
7.   2
   
8. Nfiles  
   
9.   1
   
10. Double          # calculate NACME between excited states
    
11. Pairs
    
12.   1             # Specify to calculate <n> pairs of NACME
    
13.   1 1 1 1 1 2   # S1-S2 coupling
    
14. $end

复制代码

5 MOMAP
MOMAP与BDF类似，也是模块运作的。evc模块计算振动-电子耦合(Electron-Vibration Coupling)，输入部分：

1. do_evc = 1
   
2. 
   
3. &evc
   
4. ffreq(1) = "$init_state_log_file"   # >>> set to init state log file. If use g16, fchk file with same basename should also be in current folder.
   
5. ffreq(2) = "$final_state_log_file"  # >>> the same as ffreq(1)

复制代码

如果要计算IC，还要补上这个：

1. fnacme = "$nacme_log_file"          # >>> set to nacme log file if calculate IC
   
2. proj_reorg = .t.                    # I do not find this flag in manual but Device Studio auto generated it when setting IC calculation. Prof. Zhong@ggdh told me that this flag was used to control whether to reproject the restructuring onto the intrinsic coordinates.
   
3. /

复制代码

计算kr时，接上这个：

1. do_spec_tvcf_spec = 1               # toggle fluorescence spectrum calculation
   
2. do_spec_tvcf_ft = 1                 # toggle correlation function calculation
   
3. 
   
4. &spec_tvcf                          # calculate kr
   
5. DUSHIN = .t.                        # use Duschinsky rotate
   
6. HERZ  = .f.                         # do not use Herzberg-Teller effect
   
7. EDMA = $tdm_abs debye               # >>> set to electronic dipole moment of absorption (GS)
   
8. EDME = $tdm_emi debye               # >>> set to electronic dipole moment of emission (ES)
   
9. Ead = $adia_Eex au                  # >>> set to adiabatic excitation energy
   
10. Temp = 300 K   
    
11. tmax = 1000 fs  
    
12. dt = 1 fs   
    
13. Emax = 0.3 au   
    
14. FreqScale = 1   
    
15. DSFile = "evc.dint.dat"             # input dushin file. If calculate IC rate, I recommend ALWAYS to use the dint output.
    
16. logFile = "spec.tvcf.log"   
    
17. FtFile = "spec.tvcf.ft.dat"
    
18. FoFile = "spec.tvcf.fo.dat"
    
19. FoSFile = "spec.tvcf.spec.dat"  
    
20. /   

复制代码

计算kic时，接上这个：

1. do_ic_tvcf_spec = 1                 # toggle internal conversion spectrum
   
2. do_ic_tvcf_ft = 1                   # toggle internal conversion correlation function
   
3. 
   
4. &ic_tvcf                            # calculate kic
   
5. DUSHIN = .t.   
   
6. Ead = $adia_Eex au                  # >>> set to adiabatic excitation energy
   
7. Temp = 300 K   
   
8. tmax = 1000 fs  
   
9. dt = 1 fs   
   
10. Emax = 0.3 au   
    
11. FreqScale = 1   
    
12. DSFile = "evc.cart.dat"             # input dushin file
    
13. CoulFile = "evc.cart.nac"
    
14. logFile = "ic.tvcf.log"
    
15. FtFile = "ic.tvcf.ft.dat"
    
16. FoFile = "ic.tvcf.fo.dat"
    
17. /

复制代码

可以用同一个文件算完kr与kic，这就是模块化的好处。计算完毕后，kr在spec.tvcf.log中会打印出来：

1. radiative rate     (0):     2.81284588E-11    1.16286912E+06 /s,     859.94 ns
   
2. Oscillator strength =         0.0000000000

复制代码

kic则在ic.tvcf.log中：

1. #1Energy(Hartree)       2Energy(eV) 3WaveNumber(cm-1)   4WaveLength(nm)    5radi-spectrum      6kic(s^{-1})         7log(kic)         8time(ps)
   
2.     7.58395093E-02    2.06369893E+00    1.66448483E+04    6.00786491E+02    5.05121133E-08    2.08824013E+09    9.31978044E+00      478.87212989

复制代码

本次计算遇到了负kic问题 哈哈 每个算ESD的程序逃不掉的一关XD

1. #1Energy(Hartree)       2Energy(eV) 3WaveNumber(cm-1)   4WaveLength(nm)    5radi-spectrum      6kic(s^{-1})         7log(kic)          8time(s)
   
2.     0.75991492E-01    0.20678346E+01    0.16678205E+05    0.59958492E+03   -0.20554564E-08   -0.84975391E+08               NaN   -0.11768113E-07

复制代码

由于笔者没有MOMAP的license(写邮件向帅老师乞讨一份被无视掉了哈哈哈XoX)，因此调试此问题较为困难。这一部分暂时写到这里，后续如果拿到MOMAP可能会补上。

6 FCclasses3
对于没有MOMAP license的情况，仍可以使用开源免费程序FCclasses3进行计算。对于kr，可以正常使用gen_fcc_state与gen_fcc_eldip提取数据并计算：

1. $
   
2. PROPERTY     =   EMI  ; OPA/EMI/ECD/CPL/RR/TPA/MCD/IC/NR0
   
3. MODEL        =   AH   ; AS/ASF/AH/VG/VGF/VH
   
4. DIPOLE       =   FC   ; FC/HTi/HTf
   
5. TEMP         =   298.0 ; (temperature in K)
   
6. BROADFUN     =   GAU  ; GAU/LOR/VOI
   
7. HWHM         =   0.01 ; (broadening width in eV)
   
8. METHOD       =   TD   ; TI/TD
   
9. ;VIBRATIONAL ANALYSIS
   
10. NORMALMODES  =   COMPUTE   ; COMPUTE/READ/IMPLICIT
    
11. COORDS       =   INTERNAL  ; CARTESIAN/INTERNAL
    
12. ;INPUT DATA FILES
    
13. STATE1_FILE  =   ../td_Azulene.fcc
    
14. STATE2_FILE  =   ../opt_Azulene.fcc
    
15. ELDIP_FILE   =   ../eldip_td_Azulene_fchk

复制代码

对于kic，需要找到BDF的NACME输出：

1. NACMEs including electron-translation factor:
   
2.   Gradient contribution from Final-NAC(S)-Escaled
   
3.       1       -0.0000001424       -0.0940545122       -0.0000052264
   
4.       2       -0.0645696328        0.1111467620        0.0000094841
   
5.       3        0.0666845845       -0.1257043939       -0.0000039938
   
6.       4        0.0645682184        0.1111475404        0.0000020164
   
7.       5        0.0607935966        0.1578277002       -0.0000206471
   
8.       6       -0.0666837868       -0.1257030610        0.0000070032
   
9.       7       -0.0607937180        0.1578307287        0.0000105078
   
10.       8        0.0000001690        0.0031137626       -0.0000002912
    
11.       9        0.0021204756       -0.0027583021        0.0000004147
    
12.      10       -0.0070075274       -0.0006046765       -0.0000002126
    
13.      11       -0.0021204154       -0.0027583598       -0.0000019625
    
14.      12        0.0070076098       -0.0006046934       -0.0000004938
    
15.      13        0.0965284852       -0.0503702484        0.0000173359
    
16.      14       -0.0045684190       -0.0074562651        0.0000005823
    
17.      15        0.0000009284       -0.0781063732       -0.0000022769
    
18.      16        0.0000003157        0.0048903236        0.0000025346
    
19.      17       -0.0965292017       -0.0503698835       -0.0000171479
    
20.      18        0.0045684608       -0.0074559793        0.0000023793
    
21. Sum of gradient contribution from Final-NAC(S)-Escaled
    
22.               0.0000000008        0.0000100690        0.0000000063

复制代码

将之转换成FCclasses3可识别的格式：

1. bdf2fcc nacme.out nac_td_Azulene_bdf

复制代码

bdf2fcc脚本可以在此贴中找到：FCclasses3输入文件生成脚本+ORCA旋轨耦合提取脚本。转换后得到nac_td_Azulene_bdf，这个就是FCclasses期望的NAC_FILE。IC输入文件：

1. $
   
2. PROPERTY     =   IC  ; OPA/EMI/ECD/CPL/RR/TPA/MCD/IC/NR0
   
3. MODEL        =   AH   ; AS/ASF/AH/VG/VGF/VH
   
4. DIPOLE       =   FC   ; FC/HTi/HTf
   
5. TEMP         =   298.0 ; (temperature in K)
   
6. BROADFUN     =   GAU  ; GAU/LOR/VOI
   
7. HWHM         =   0.01 ; (broadening width in eV)
   
8. METHOD       =   TD   ; TI/TD
   
9. ;VIBRATIONAL ANALYSIS
   
10. NORMALMODES  =   COMPUTE   ; COMPUTE/READ/IMPLICIT
    
11. COORDS       =   INTERNAL  ; CARTESIAN/INTERNAL
    
12. ;INPUT DATA FILES
    
13. STATE1_FILE  =   ../td_Azulene.fcc
    
14. STATE2_FILE  =   ../opt_Azulene.fcc
    
15. NAC_FILE     =   ../nac_td_Azulene_bdf

复制代码

计算完成后，在输出文件末尾找到结果：

1. # for kr
   
2. 
   
3.   Radiative rate constant
   
4.   ------------------------
   
5. kr(s-1) =   8.8228E+05
   
6. 
   
7. # for kic
   
8. 
   
9. ========================================================
   
10. IC rate constant (s-1)  2.507E+09
    
11. (for Ead =     2.064 eV)
    
12. ========================================================

复制代码

我们看到FCclasses计算值(8.82E+05/2.51E+09)与MOMAP计算值(1.16E+06/2.09E+09)还是比较接近的。汇总结果：

	kr	kIC
S1-S0	8.82E+05	2.51E+09
S2-S0	4.86E+07	1.92E+02
S2-S1	-	6.47E+04

(这里笔者感觉IC计算值偏小，可能是因为没有引入溶剂效应/HT效应等引起的。不过表观上看仍然是符合观测现象的)

我们计算kasha量子产率：

以及anti-kasha量子产率：

注意到由于S2/S1、S2/S0的内转换速率较低，Azulene的anti-kasha量子产率远远高于kasha量子产率，这解释了Azulene的anti-kasha发射来源。

7 总结

使用Gaussian、BDF与MOMAP/FCclasses计算了Azulene的anti-kasha现象。笔者研究anti-kasha时间不长，如您发现有错误请在评论区斧正。过去的许多工作由于缺乏可靠的计算方法，常常只能凭激发能与实验光谱对应情况"鉴定"anti-kasha。但随着BDF的逐渐成熟，论证anti-kasha有了有力的工具，料想以后在anti-kasha的研究的方面会有更多的优秀工作出现。

wzkchem5

这个计算流程需要一个前提假设，就是ISC可以忽略，这一点对azulene成立，但对其他分子未必成立。
至于IC（尤其是S1-S0 IC）速率比实验值小，可能确实是因为溶剂效应，我试过在azulene附近加一个甲醇分子，IC速率就上升到10^11量级了，和甲醇里的实验结果符合。但隐式溶剂模型不行，因为隐式溶剂模型描述不了溶剂分子振动的贡献

wal

wzkchem5 发表于 2025-8-4 22:47
这个计算流程需要一个前提假设，就是ISC可以忽略，这一点对azulene成立，但对其他分子未必成立。
至于IC（ ...

哦对，我忘记说前提是忽略ISC了，感谢老师提醒

溶剂这块我一开始是忘记了，因为DS生成TDDFT-NAC输入那个界面溶剂是灰的，没想着加溶剂。后来想起来之后感觉这个显式溶剂考虑起来可能非常麻烦，就顺着做下来了。原因在于我最近有个体系也是要考虑溶剂下算荧光寿命，发现似乎给不同数量的显式水对速率的影响不一样，有时候给的多还能把kic干下来XD 我猜有的体系水到处形成氢键之后反而会阻碍低频振动
后面有空的话我修缮修缮，做个溶剂对比