用Gaussian寻找圆锥交叉点

一、圆锥交叉的含义

圆锥交叉（conical intersection, CI）简言之是指两个态的势能面交叉的地方，此时两个态的能量简并。在圆锥交叉区域，体系可以从激发态以无辐射形式回到基态。关于圆锥交叉的更深入的理论细节可参考Conical intersections: theory, computation and experiment一书。此外，点击文末“阅读原文”可打开一份非常不错的关于光化学计算的讲义。如果打开速度较慢，可在留言区的百度网盘链接获得。

在下图所示的势能面中，反应物经过垂直跃迁成为激发态，并弛豫到激发态势能面的能量极小点。在激发态势能面上还可能跨过能垒（若有的话），到达圆锥交叉点，随后便进入基态势能面，逐步转化为产物。

需要注意的是，S0与T1的势能面交叉称为MECP点更为合适，CI点最初只用于同一自旋多重度两个势能面之间的交叉，现在文献中把圆锥示意图泛泛地标在任两个自旋多重度间的势能面交叉，只表示交叉示意，不代表其是真的是圆锥形状的交叉，读者在读文献时务必注意。

如果仅是搜索S0与T1之间的MECP点，一般基态(U)DFT就能做，有不少程序（easyMECP, sobMECP，及二者的核心MECP程序）支持，计算量不大、操作较为方便。而对于同一自旋多重度的CI点，若用TDDFT，其激发态虽是多行列式的，但是基态仍是单行列式，既然基态、激发态简并，那么从原理上讲基态也应该考虑多个行列式的线性组合。而CASSCF是一种多组态方法，可以同时算多个根，每个根都是多行列式的，因此既可以搜索CI点，也可以搜索MECP点，若只是探究光化学、光物理过程的机理、对精度要求不十分苛刻时，是理想的选择。需要更高精度的话，后续仍需做(X)MS-CASPT2或MRCI。不过这些计算耗时都大，学习门槛也不低，本文暂只涉及用CASSCF寻找圆锥交叉点。

二、计算实例

本文以Exploring Chemistry with Electronic Structure Methods (Third Edition)中的例8.12为例，介绍相关过程的计算。这个例子计算的是苯到盆苯（benzvalene）的异构化过程中的光化学过程，盆苯是基态下能量较高的极小点。所有计算用Gaussian 16 C.01完成。

(1) 用CASSCF优化苯的基态结构。这一步建立在《用Gaussian做CASSCF计算》一文的计算基础之上。在挑选了合适的活性空间后，做结构优化及频率计算：

%oldchk=benzene_cas.chk
%chk=benzene_gs_opt.chk
#p opt freq cas(6,6)/6-31G(d) geom=allcheck guess=read

benzene_cas.chk是用DFT或其他方法优化结构，并做CASSCF单点计算生成的chk文件。

(2) 用CASSCF计算第一激发态，获得垂直激发能，同时观察轨道，确保活性空间的正确性。

%oldchk=benzene_gs_opt.chk
%chk=benzene_es.chk
#p cas(6,6,nroot=2)/6-31G(d) geom=allcheck guess=read

(3) 用CASSCF优化激发态的结构并做频率计算。

%oldchk=benzene_es.chk
%chk=benzene_es_opt.chk
#p opt freq cas(6,6,nroot=2)/6-31G(d) geom=allcheck guess=read

优化得到的激发态与基态相比，C−C键键长从1.396 Å增加到1.434 Å。

(4) 寻找圆锥交叉点。这是一项需要经验和技巧的任务。在exploring3中，使用的方法是利用柔性扫描的结果作为优化交叉点的初始结构。对C1−C2−C3−C6二面角进行扫描：

%oldchk=benzene_es_opt.chk
%chk=benzene_scan.chk
#p opt=modredundant cas(6,6,nroot=2)/6-31G(d) geom=allcheck guess=read nosymm
 
1 2 3 6 S 5 7.5

共扫描得6个结构，从输出文件中可以看到每个结构的垂直激发能在逐渐变小。这是符合预期的，因为在圆锥交叉点处基态和激发态的能量简并，垂直激发能应该趋近于0。

扫描得到的第6个结构如下图所示：

再以此结构为初始，进行圆锥交叉点优化：

%oldchk=benzene_scan.chk
%chk=benzene_ci.chk
#p opt=conical casscf(6,6,nroot=2)/6-31g(d) guess=read nosymm
 
Benzene
 
0 1
 C         -0.00080300    1.41440000    0.65507800
 C          1.05274800    0.71961300   -0.07247700
 C          1.17424100   -0.73123000    0.00730800
 C          0.00097200   -1.47535000    0.15754600
 C         -1.17472000   -0.73263700    0.01547600
 C         -1.05143700    0.71680700   -0.07124800
 H         -0.00175900    2.49003400    0.65966300
 H          1.70776200    1.26435800   -0.73164000
 H          2.13962100   -1.19172400   -0.09740800
 H          0.00241700   -2.54058600    0.29322400
 H         -2.14003800   -1.19374700   -0.08641800
 H         -1.70900300    1.26006300   -0.72910300
 
0.500000000.50000000

优化圆锥交叉点的关键词为opt=conical，此时CASSCF默认为态平均CASSCF，因此需要在文件末尾加上权重。这一点在Exploring3中未提到，所给的例子中也没有写权重，实际程序运行时会报错。此外，这一步不能直接使用geom=allcheck读取scan步骤中的坐标信息，因为会将扫描的冗余内坐标沿用过来，影响圆锥交叉的优化。所以要重新保存结构，或使用geom=(allcheck,newdefinition)关键词。

优化得到的圆锥交叉点的垂直激发能为0.0012 eV，结构如下图所示：

(5) 观察圆锥交叉点的结构，可以看到，若将C2−C6键拉长，则接近苯的结构，而将C2−C6键缩短，则向盆苯的结构靠近，因此若想得到盆苯的结构，我们尝试手动将键长缩小，如至1.8 Å，以此结构为初始，进行结构优化。

%oldchk=benzene_ci.chk
%chk=benzene_pr.chk
#p opt freq casscf(6,6)/6-31g(d) guess=read
 
Benzene
 
0 1
 C     0.00046100    1.40464600    0.71242300
 C     0.89962637    0.71088527   -0.12565894
 C     1.16331800   -0.73594400    0.00659700
 C     0.00117700   -1.49102900    0.16332400
 C    -1.16241300   -0.73730500    0.01061800
 C    -0.90097037    0.70984973   -0.12239106
 H    -0.00007100    2.47952800    0.74928700
 H     1.53720837    1.27221027   -0.78782094
 H     2.14440600   -1.16412600   -0.07883700
 H     0.00209900   -2.55251600    0.32615000
 H    -2.14327100   -1.16667000   -0.07145300
 H    -1.54156837    1.27047273   -0.78223506

不幸的是，优化得到的结构有虚频，是一个过渡态的结构。观察其振动模式：

可知这可能是两个盆苯的异构过程的过渡态。因此，一个比较合适的做法就是借助IRC分析，得到靠近盆苯能量极小点的结构，再做最终的结构优化。可以先做一步粗糙的IRC分析：

%oldchk=benzene_pr.chk
%chk=benzene_irc.chk
#p irc=(rcfc,recorrect=never) casscf(6,6)/6-31g(d) guess=read geom=allcheck

结果如下图所示，IRC曲线很对称，观察两端的结构，确实印证了上述猜想。但是由于IRC所走的步数不够多，还不是很接近盆苯的能量极小点，如果直接用此时两端的结构进行优化，很有可能还是会回到过渡态的结构。

此时，可以重新做一次IRC分析，并只向一个方向走步，增大步数：

%oldchk=benzene_pr.chk
%chk=benzene_irc2.chk
#p irc=(rcfc,forward,maxpoints=50,recorrect=never) casscf(6,6)/6-31g(d) guess=read geom=allcheck

得如下结果：

最后一点的分子结构为

最后，以这一步的最终结构为初始结构，优化能量极小点：

%oldchk=benzene_irc2.chk
%chk=benzene_pr2.chk
#p opt freq casscf(6,6)/6-31g(d) guess=read
 
Benzene
 
0 1
 C    -0.30169200    0.96943500    0.72327700
 C     0.79642600    0.70088900   -0.26809600
 C     1.17355300   -0.76059400   -0.20523300
 C     0.04989100   -1.40825900    0.13313200
 C    -1.00131100   -0.33628900    0.28834500
 C    -0.69727600    0.80329900   -0.64511000
 H    -0.48556700    1.70567400    1.48406900
 H     1.49687100    1.48752000   -0.52796400
 H     2.14667400   -1.16482800   -0.42415000
 H    -0.10416600   -2.45960300    0.28151200
 H    -2.01412200   -0.53455000    0.59125500
 H    -1.27067600    1.38386100   -1.34125200

最终得到盆苯的结构如下：

在第一段提到的教程中正好有一个描述此过程的势能面，如下图所示：

总结一下便是，基态的苯在光照下得到激发态，在激发态势能面上行走到达能量极小点，并可能进一步跨越能垒（也可能没有能垒）到达圆锥交叉点，顺着圆锥交叉点进入基态势能面上的两个盆苯互变异构的过渡态，再得到盆苯的能量极小结构。同时，在图中也同时显示了在基态的势能面上也存在苯到盆苯过程的过渡态结构。但这个过渡态的能垒很高，不易发生。