Anti-Kasha荧光的ESD计算

Posted Aug 2, 2025 Updated Aug 4, 2025

By Bane Dysta

14 min read

Anti-Kasha荧光的ESD计算

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发射的条件：

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

1.2 理论描述

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

\[\Phi_{S_2} = \frac{k_{r,S_2}}{k_{r,S_2} + k_{IC,S_2} + k_{IC,S_2→S_1}}\]

其中：

$k_{r,S_2}$：S₂态到基态的辐射跃迁速率
$k_{IC,S_2}$：S₂态到基态内转换速率
$k_{IC,S_2→S_1}$：S₂→S₁内转换速率

只有当$k_{r,S_2}$与其余两项之和相当或更大时，才能观察到明显的Anti-Kasha发射。

2 计算流程

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

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

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

3 Gaussian

平平无奇的三个opt：

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

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

# td B3LYP/def2svp scrf em=gd3bj

对于使用MOMAP计算ESD的情况，我们需要在输出中提取这些量：

基态能量：SCF Done后面的数值

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

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

Excitation energies and oscillator strengths:
 
 Excited State   1:      Singlet-A      1.7066 eV  726.49 nm  f=0.0074  <S**2>=0.000
      34 -> 35        -0.70605
 This state for optimization and/or second-order correction.
 Total Energy, E(TD-HF/TD-DFT) =  -385.528161673    
 Copying the excited state density for this state as the 1-particle RhoCI density.
 
 Excited State   2:      Singlet-A      3.3980 eV  364.88 nm  f=0.0099  <S**2>=0.000
      33 -> 35        -0.46258
      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结构的。

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

$\mu_{trans1}=2.5417 \times \sqrt{0.1778}=1.0718$

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

4 BDF

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

4.1 基态与激发态间NAC

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

$resp
Quad
FNAC
Norder
  1
Method
  2  
Nfiles
  1
Single          # calculate NACME between excited state and ground state
States
  1
  1 1 1         # <istore number> <Irreducible representation order> <root order>
$end

类似地，S₁态与基态之间的NAC只要把States部分改为：

States
  1
  1 1 2

即可进行计算。

4.2 激发态间NAC

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

$resp
Quad
FNAC
Norder
  1
Method
  2
Nfiles  
  1
Double          # calculate NACME between excited states
Pairs
  1             # Specify to calculate <n> pairs of NACME
  1 1 1 1 1 2   # S1-S2 coupling
$end

5 MOMAP

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

do_evc = 1

&evc
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.
ffreq(2) = "$final_state_log_file"  # >>> the same as ffreq(1)

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

fnacme = "$nacme_log_file"          # >>> set to nacme log file if calculate IC
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.
/

计算kr时，接上这个：

do_spec_tvcf_spec = 1               # toggle fluorescence spectrum calculation
do_spec_tvcf_ft = 1                 # toggle correlation function calculation

&spec_tvcf                          # calculate kr
DUSHIN = .t.                        # use Duschinsky rotate
HERZ  = .f.                         # do not use Herzberg-Teller effect
EDMA = $tdm_abs debye               # >>> set to electronic dipole moment of absorption (GS)
EDME = $tdm_emi debye               # >>> set to electronic dipole moment of emission (ES)
Ead = $adia_Eex au                  # >>> set to adiabatic excitation energy
Temp = 300 K    
tmax = 1000 fs  
dt = 1 fs   
Emax = 0.3 au   
FreqScale = 1   
DSFile = "evc.dint.dat"             # input dushin file. If calculate IC rate, I recommend ALWAYS to use the dint output.
logFile = "spec.tvcf.log"   
FtFile = "spec.tvcf.ft.dat" 
FoFile = "spec.tvcf.fo.dat" 
FoSFile = "spec.tvcf.spec.dat"  
/   

计算kic时，接上这个：

do_ic_tvcf_spec = 1                 # toggle internal conversion spectrum
do_ic_tvcf_ft = 1                   # toggle internal conversion correlation function

&ic_tvcf                            # calculate kic
DUSHIN = .t.    
Ead = $adia_Eex au                  # >>> set to adiabatic excitation energy
Temp = 300 K    
tmax = 1000 fs  
dt = 1 fs   
Emax = 0.3 au   
FreqScale = 1   
DSFile = "evc.cart.dat"             # input dushin file
CoulFile = "evc.cart.nac"
logFile = "ic.tvcf.log"
FtFile = "ic.tvcf.ft.dat"
FoFile = "ic.tvcf.fo.dat"
/

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

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

kic则在ic.tvcf.log中：

 #1Energy(Hartree)       2Energy(eV) 3WaveNumber(cm-1)   4WaveLength(nm)    5radi-spectrum      6kic(s^{-1})         7log(kic)         8time(ps)
    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的程序逃不掉的一关😸

 #1Energy(Hartree)       2Energy(eV) 3WaveNumber(cm-1)   4WaveLength(nm)    5radi-spectrum      6kic(s^{-1})         7log(kic)          8time(s)
    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提取数据并计算：

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

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

 NACMEs including electron-translation factor:
  Gradient contribution from Final-NAC(S)-Escaled
     -0.0000001424       -0.0940545122       -0.0000052264
     -0.0645696328        0.1111467620        0.0000094841
      0.0666845845       -0.1257043939       -0.0000039938
      0.0645682184        0.1111475404        0.0000020164
      0.0607935966        0.1578277002       -0.0000206471
     -0.0666837868       -0.1257030610        0.0000070032
     -0.0607937180        0.1578307287        0.0000105078
      0.0000001690        0.0031137626       -0.0000002912
      0.0021204756       -0.0027583021        0.0000004147
     -0.0070075274       -0.0006046765       -0.0000002126
     -0.0021204154       -0.0027583598       -0.0000019625
      0.0070076098       -0.0006046934       -0.0000004938
      0.0965284852       -0.0503702484        0.0000173359
     -0.0045684190       -0.0074562651        0.0000005823
      0.0000009284       -0.0781063732       -0.0000022769
      0.0000003157        0.0048903236        0.0000025346
     -0.0965292017       -0.0503698835       -0.0000171479
      0.0045684608       -0.0074559793        0.0000023793
 Sum of gradient contribution from Final-NAC(S)-Escaled
              0.0000000008        0.0000100690        0.0000000063

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

bdf2fcc nacme.out nac_td_Azulene_bdf

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

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

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

 # for kr

  Radiative rate constant 
  ------------------------
 kr(s-1) =   8.8228E+05

# for kic

 ========================================================
 IC rate constant (s-1)  2.507E+09
 (for Ead =     2.064 eV)
 ========================================================

我们看到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量子产率：

\[\phi_{S1} = \frac{8.82 \times 10^5}{8.82 \times 10^5 + 2.51 \times 10^9} = 3.51 \times 10^{-4}\]

以及anti-kasha量子产率：

\[\phi_{anti-kasha} = \frac{4.86 \times 10^7}{4.86 \times 10^7 + 1.92 \times 10^2 + 6.47 \times 10^4} = 9.9 \times 10^{-1}\]

注意到由于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的研究的方面会有更多的优秀工作出现。

Quantum Chemistry, BDF

This post is licensed under CC BY 4.0 by the author.