Anti-Kasha荧光的计算
Anti-Kasha荧光是违反Kasha规则的发光现象，即分子从高激发态（如S₂、S₃）直接发射荧光，而非仅从最低激发态S₁发射。此类现象在有机染料中较为罕见，但在某些特殊结构的分子中可以观察到。本文将以Azulene为例尝试对此效应进行计算。
1 理论背景
1.1 Kasha规则与Anti-Kasha现象
Kasha规则指出，分子在光激发后，无论最初被激发到哪个高能态，最终的荧光发射通常仅来源于最低激发态S₁。这是因为高激发态之间的能量差相对较小，根据能隙定律，内转换速率随能隙呈指数下降，导致S₂→S₁的内转换速率远大于S₂的辐射速率。
然而，在某些特殊情况下，分子可以违反Kasha规则，直接从S₂或更高激发态发射荧光，这被称为Anti-Kasha发射。有利于发生Anti-Kasha发射的条件：
- 大的S₂-S₁能隙：当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_{nr,S_2} + k_{IC,S_2→S_1}}\]其中：
- $k_{r,S_2}$：S₂态的辐射跃迁速率
- $k_{nr,S_2}$：S₂态的其他非辐射过程速率
- $k_{IC,S_2→S_1}$：S₂→S₁内转换速率
只有当$k_{r,S_2}$与$k_{IC,S_2→S_1}$相当或更大时，才能观察到明显的Anti-Kasha发射。
2 计算流程
Anti-Kasha荧光的理论计算需要获得以下量：
- 几何优化：基态和激发态的平衡几何结构、正则振动模式、绝热激发能、跃迁偶极矩。拥有TDDFT解析hessian的g16胜任此计算；
- 非绝热耦合：态间的NAC矩阵元。激发态间的NAC目前只有Q-Chem和BDF两家胜任，本文采用BDF。
- 激发态动力学参数：利用前述中间量计算激发态衰减速率。MOMAP拥有对BDF进行TDDFT-NAC输出的解析支持，因此本文采用MOMAP。
接下来，笔者将依次计算这些量。
3 Gaussian
平平无奇的三个opt：
1
2
3
# 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单点：
1
# td B3LYP/def2svp scrf em=gd3bj
在输出中，需要提取这些量：
- 能量
1 2 3
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 - 激发态能量 ``` 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
1
- TDM
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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
$\mu_{trans}=2.5417*\sqrt0.1778=1.0718$
| | 能量 / a. u. | TDMabs / Debye | TDMemi / Debye |
| :-: | :-----------: | :----------: | :----------: |
| S0 | -385.60401168 | - | - |
| S1 | -385.52816168 | 1.0939 | 1.0718 |
| S2 | -385.47282175 | 0.5190 | 2.8292 |
## 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 #
1
2
3
4
5
### 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
1
2
3
## 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 g16, fchk file should also be in current folder ffreq(2) = “$final_state_log_file” # »> same as ffreq(1)
1
如果要计算IC，还要补上这个：
fnacme = “$nacme_log_file” # »> set to nacme log file if calculate IC proj_reorg = .t. # I do not find this keyword in manual but device studio auto generated this when setting IC calculation. /
1
计算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”
/
1
计算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” /
1
2
3
4
5
6
可以用同一个文件算完kr与kic，模块化的好处。
## 6 结果
### 6.1 MOMAP
kr在spec.tvcf.log中会打印出来：
radiative rate (0): 2.81284588E-11 1.16286912E+06 /s, 859.94 ns Oscillator strength = 0.0000000000
1
kic则在
#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
1
2
本次计算遇到了负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版权，调试此问题较为困难，结果分析得后面补上。