用GAMESS中的Spin-flip TD-DFT找S0/S1交叉点
用GAMESS中的Spin-flip TD-DFT找S0/S1交叉点
寻找势能面交叉点是激发态的研究中经常遇到的问题。不同自旋多重度的势能面交叉点相关的介绍可以参考本公众号之前所发关于MECP系列文章。自旋多重度相同的势能面的交叉点常称为圆锥交叉（conical intersection, CI），我们也曾介绍过如何用CASSCF方法寻找CI点。然而CASSCF方法涉及活性空间的选择等问题，在使用上不是特别方便，对稍大一些的体系，其计算量往往也难以承受。TD-DFT是当前激发态计算中最常用的方法，不少程序支持使用TD-DFT来寻找CI点，如GAMESS、ORCA等。然而，对于S0和S1势能面的交叉点，则需要特别注意。虽然上述两个程序的TD-DFT都支持寻找S0/S1交叉点，而且碰巧的是，这两个程序官方给出的算例都是寻找S0/S1交叉点，但实际上TD-DFT在描述参考态（S0）与激发态的交叉点时是有缺陷的，原理上无法描述S0/Sn交叉点。这点在ORCA 5.0.2版的手册8.3.12节中已经指出，也有不少文献中提及此点，如J. Phys. Chem. A, 2009, 113, 12749.等文章。
若想在TD-DFT级别找S0/Sn交叉点，可以尝试使用Spin-flip框架下的TDDFT方法（简写作SF-TDDFT或SFDFT）。在SF-TDDFT中，常以三重态为参考态，通过翻转一个alpha电子为beta电子，这样S0和S1均是由多个行列式表示，原理上就可以描述S0/S1交叉点了。在SF框架下，不少测试表明Hartree-Fock成份为50%的BHandHLYP泛函在SF-TDDFT中有不错的结果。在常见的程序中，Q-Chem、GAMESS、ORCA均支持SF-TDDFT来寻找CI点。本文我们介绍如何用GAMESS来找S0/S1交叉点（GAMESS程序的安装见《GAMESS简易编译教程》一文）。以后我们再介绍如何在ORCA中做Spin-flip计算。
本文我们尝试用SF-TDDFT方法来寻找J. Phys. Chem. A 2021, 125, 559中的一个体系的S0/S1交叉点，其分子结构如下：
左图为俯视图，右图为侧视图
首先需要做一个SF单点计算，以挑选所研究的态。输入文件如下：
$CONTRL SCFTYP=ROHF MULT=3 ICHARG=0 RUNTYP=ENERGY
DFTTYP=BHHLYP ISPHER=0 MAXIT=160 NPRINT=-5 TDDFT=SPNFLP $END
$BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 NPFUNC=1 DIFFSP=.F. DIFFS=.F. $END
$SYSTEM MWORDS=2000 MEMDDI=2000 $END
$SCF DIRSCF=.T. DIIS=.T. $END
$TDDFT NSTATE=5 $END
$DATA
SF Single Point
C1
C 6.0 -1.94899995 0.80467730 -0.67386604
C 6.0 -1.94899868 0.80467686 0.67386507
C 6.0 -1.07616035 0.15778413 1.66978959
C 6.0 -1.07616095 0.15778444 -1.66978980
C 6.0 0.31004511 -0.06035182 1.50729154
C 6.0 0.31004405 -0.06035342 -1.50728982
H 1.0 -2.80935315 1.30771907 -1.12742642
H 1.0 -2.80935213 1.30771812 1.12742688
Si 14.0 1.22293498 0.57944668 0.00000050
C 6.0 2.98615545 -0.05428175 -0.00000025
H 1.0 3.52977672 0.31101873 -0.88562038
H 1.0 3.52977835 0.31102088 0.88561795
H 1.0 3.03733612 -1.15418221 0.00000097
C 6.0 1.25579636 2.45305625 -0.00000149
H 1.0 1.78241688 2.82955629 -0.89164657
H 1.0 0.24176141 2.87838801 -0.00000291
H 1.0 1.78241581 2.82955906 0.89164301
C 6.0 -1.69749080 -0.24896792 -2.86203126
H 1.0 -2.76500279 -0.05752403 -2.99315546
C 6.0 -0.99080690 -0.90211528 -3.86157863
H 1.0 -1.50117893 -1.21715270 -4.77367866
C 6.0 0.36885242 -1.14909276 -3.69396105
H 1.0 0.93797286 -1.65974548 -4.47296814
C 6.0 1.00150241 -0.72363461 -2.53088211
H 1.0 2.07260433 -0.90611006 -2.42219277
C 6.0 -1.69749083 -0.24897008 2.86203000
H 1.0 -2.76500325 -0.05752760 2.99315273
C 6.0 1.00150367 -0.72363246 2.53088399
H 1.0 2.07260598 -0.90610710 2.42219663
C 6.0 -0.99080715 -0.90211688 3.86157778
H 1.0 -1.50118004 -1.21715584 4.77367680
C 6.0 0.36885278 -1.14909243 3.69396183
H 1.0 0.93797200 -1.65974468 4.47297014
$END在GAMESS中支持ROHF和UHF两种形式的参考态，此处我们选用ROHF参考态，自旋多重度设为3。选择BHandHLYP泛函（在GAMESS中写作BHHLYP）和6-31G(d,p)基组。在CONTRL中使用TDDFT=SPNFLP关键词来实现SF-TDDFT的计算。TDDFT中的NSTATE=5关键词表示计算5个激发态。
单点计算的Spin-flip结果如下：
----------------------------
SUMMARY OF SPIN-FLIP RESULTS
----------------------------
STATE ENERGY EXCITATION S-SQUARED
1 NEGATIVE ROOT(S) FOUND.
1 A -908.4788531272 -3.392 0.0110
0 A -908.3541900392 0.000 (REFERENCE STATE)
2 A -908.3450078844 0.250 1.9898
3 A -908.2975242437 1.542 0.1496
4 A -908.2852747509 1.875 1.0330
5 A -908.2747854242 2.161 1.0291结果按照由低到高的顺序给出各个态的能量。此处出现了一个激发能为负的态，其S2期望值为0.0110，为单重态，这个态实际上是该分子的真正基态S0。而此处第二行的参考态由于我们将其多重度设为3，它就是分子的实际T1态。态2的S2期望值为1.9898，接近2，因此它也是一个三重态，即T2。态3的S2期望值为0.1496，可指认为S1态。而下面两个态则有较大的自旋污染。这是由于GAMESS中实现的不是自旋匹配的SF-TDDFT方法（ORCA中亦如此），所以会出现自旋污染。总之，通过单点计算，我们找出要研究的态为1和3。
接下来便可优化交叉点的结构，输入文件如下：
$CONTRL SCFTYP=ROHF MULT=3 ICHARG=0 RUNTYP=CONICAL
DFTTYP=BHHLYP ISPHER=0 MAXIT=160 NPRINT=-5 TDDFT=SPNFLP $END
$BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 NPFUNC=1 DIFFSP=.F. DIFFS=.F. $END
$SYSTEM MWORDS=2000 MEMDDI=2000 $END
$SCF DIRSCF=.T. DIIS=.T. $END
$STATPT NSTEP=500 $END
$CONICL OPTTYP=PENALTY IXROOT(1)=1,3 SIGMA=8.0 $END
$TDDFT NSTATE=5 TAMMD=.T. $END
[此处略去分子坐标]此时CONTRL中的RUNTYP设为CONICAL，CONICL中的OPTTPY=PENALTY表示使用penalty-constrained优化算法，也可使用GAMESS中的默认算法BPUPD；IXROOT(1)=1,3表示寻找第1个态和第3个态的交叉点；SIGMA=8.0是PENALTY算法中的一个参数，其默认值为3.5，若优化出的交叉点的能量差较大，可以尝试增大SIGMA的值来重新优化，以获得能量差更小的结构。
随着优化的进行，可以看到态1和3的能量会逐渐接近，可想而知，这两个态的序号按道理会变成两个连续的数字，这也是激发态结构优化中经常到的势能面交叉问题。此例中，随着优化的进行，SF的结果变成如下所示：
----------------------------
SUMMARY OF SPIN-FLIP RESULTS
----------------------------
STATE ENERGY EXCITATION S-SQUARED
0 NEGATIVE ROOT(S) FOUND.
0 A -908.2457867532 0.000 (REFERENCE STATE)
1 A -908.2335571615 0.333 1.7029
2 A -908.2273163680 0.503 0.2076
3 A -908.2259141432 0.541 0.2526
4 A -908.1568193434 2.421 1.0312
5 A -908.1386487788 2.915 1.1457各态的能量顺序已经发生了变化。此时，我们要将任务终止，取优化的最后一步的结构，将IXROOT(1)设为2,3重新优化。
最终成功优化出了S0/S1交叉点结构。最后一步的SF结果如下：
----------------------------
SUMMARY OF SPIN-FLIP RESULTS
----------------------------
STATE ENERGY EXCITATION S-SQUARED
0 NEGATIVE ROOT(S) FOUND.
0 A -908.3443444265 0.000 (REFERENCE STATE)
1 A -908.3313910094 0.352 1.1440
2 A -908.3193995875 0.679 0.8748
3 A -908.3186507400 0.699 0.0811
4 A -908.2290731307 3.137 1.1484
5 A -908.2171758204 3.460 1.2347最终第2个态的自旋污染还是较大的，S2期望值已经接近1，可以看作是有较多三重态的成分的混合。当然，第1个态S2期望值也接近1，此例可能还需要对结果作更深入的研究，也欢迎留言讨论。
最后，我们将优化出的结果，与原文在TDDFT和XMS-CASPT2级别下找到的交叉点对比如下：
红色为SF-TDDFT，白色为TD-DFT，蓝色为XMS-CASPT2
结构非常接近。其中关键参数C3-C2-C4-C1二面角（该二面角反映了分子的扭曲程度）如下：
- SF-TDDFT：−103.3
- TDDFT：−96.3
- XMS-CASPT2：−100.9
小 结
本文简单介绍了如何用GAMESS程序中的SF-TDDFT方法寻找S0/S1交叉点。交叉点的寻找不是一件容易的事，在用SF-TDDFT优化结构过程中，最烦琐的问题就是能量顺序的变化和自旋污染问题，因此在优化过程中要随时查看计算结果，并做出相应的调整。此外，交叉点的优化也很依赖于初始结构，文献和经验的积累也非常重要。
致 谢
感谢jxzou在成文中的建议。