计算化学公社

标题: Anti-Kasha的Excited State Dynamics计算流程分享 [打印本页]
作者
Author: wal    时间: 2025-8-4 16:59
标题: Anti-Kasha的Excited State Dynamics计算流程分享
本帖最后由 wal 于 2025-8-5 01:18 编辑

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发射的量子产率可以用竞争动力学来描述：

(, 下载次数 Times of downloads: 0)

其中：
- 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:

  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.

  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
复制代码

(, 下载次数 Times of downloads: 0)

而对于使用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

  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

  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

  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

  3.   Radiative rate constant
  4.   ------------------------
  5. kr(s-1) =   8.8228E+05

  7. # for kic

  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量子产率：

(, 下载次数 Times of downloads: 0)


以及anti-kasha量子产率：

(, 下载次数 Times of downloads: 0)

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




作者
Author: wzkchem5    时间: 2025-8-4 22:47
这个计算流程需要一个前提假设，就是ISC可以忽略，这一点对azulene成立，但对其他分子未必成立。
至于IC（尤其是S1-S0 IC）速率比实验值小，可能确实是因为溶剂效应，我试过在azulene附近加一个甲醇分子，IC速率就上升到10^11量级了，和甲醇里的实验结果符合。但隐式溶剂模型不行，因为隐式溶剂模型描述不了溶剂分子振动的贡献
作者
Author: wal    时间: 2025-8-4 23:13
本帖最后由 wal 于 2025-8-4 23:22 编辑
wzkchem5 发表于 2025-8-4 22:47
这个计算流程需要一个前提假设，就是ISC可以忽略，这一点对azulene成立，但对其他分子未必成立。
至于IC（ ...

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

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




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3