用Gaussian做CASSCF计算

一、CASSCF方法简介

完全活性空间自洽场（complete active space self-consistent field, CASSCF）是一种组态相互作用（configuration interaction, CI）方法。以下我们用尽可能简单的语言阐述CI方法的核心原理，不涉及复杂的数学公式。在处理氢分子时，如果基态组态σ1s2和一定量的激发组态σ*1s2混合，会使体系能量降低，也就意味着改进了基态波函数，这种方法称为CI方法。在CI中，可以有多个激发组态参加。full CI就是在所选基组范围内，考虑所有可能的激发组态。一般情况下，full CI是不现实的，需要对组态进行截断，如只考虑单激发组态则称为CIS方法，再加入双激发组态则为CISD方法。此类CI方法中，只将各组态的组合系数进行变分优化。若在计算时，同时优化组态系数和各组态基函数的系数（也即分子轨道系数），则称为多组态自洽场（multi-configuration self-consistent field, MCSCF）方法。在MCSCF方法中，自然也是考虑越多的组态则结果越精确，但同时计算量也越大。CASSCF就是其中一种方法，将分子轨道分为三组。能量低的占据轨道和能量高的虚轨道不参与CI计算，剩下的为活性空间，通常由HF计算得到的若干前线轨道或在所研究问题中变化显著的轨道组成。参加CI计算的组态包括电子在活性空间中全部可能的组态。合理地选择活性轨道空间是CASSCF计算中的关键。CASSCF活性空间选取实例见《激发态计算之如何选取特定的分子轨道作为活性空间》一文。

二、计算实例

在Gaussian中，CASSCF方法可简写为CAS，并用两个数字表示活性空间，即CAS(n,m)，n表示活性空间中的电子数目，m表示活性空间中的轨道数目。例如CAS(4,6)，则表示活性空间中有6个轨道，其中两个为占据轨道，填充4个电子，其他4个轨道为空轨道。

通常对一个体系先做一个小基组的HF计算，用来观察轨道的特征，并选取活性轨道。通常不用太大的基组，因为那样会导致轨道比较弥散，不便于指认。程序里活性轨道的选取是从HOMO的位置往下数，LUMO的位置往上数，例如CAS(4,6)就是HOMO−1、HOMO、LUMO、LUMO+1、LUMO+2、LUMO+3的位置，一共6个轨道。若想要的轨道不在这6个位置上，需要用guess=alter关键词来调整轨道顺序，将想要的轨道调换进来。以下以苯分子的π→π*激发态计算为例，说明整个过程。

(1) 优化分子结构，例如可用M06-2X/def2-SVP水平进行：

%chk=benzene_opt.chk
#p opt freq m062x/def2svp
 
Benzene
 
0 1
 C        0.00000000    1.38980400    0.00000000
 C        1.20360500    0.69490200    0.00000000
 C        1.20360500   -0.69490200    0.00000000
 C        0.00000000   -1.38980400    0.00000000
 C       -1.20360500   -0.69490200    0.00000000
 C       -1.20360500    0.69490200    0.00000000
 H        0.00000000    2.47523400    0.00000000
 H        2.14361600    1.23761700    0.00000000
 H        2.14361600   -1.23761700    0.00000000
 H        0.00000000   -2.47523400    0.00000000
 H       -2.14361600   -1.23761700    0.00000000
 H       -2.14361600    1.23761700    0.00000000
 

(2) 做HF计算，观察轨道，选取活性空间：

%oldchk=benzene_opt.chk
%chk=benzene_hf.chk
#p hf/sto-3g geom=allcheck
 

计算完成后将chk文件转为fchk文件，并用Multiwfn程序打开，其中一些轨道如图所示：

图中O表示占据（occupied）轨道，V表示虚（virtual）轨道。当我们想研究π→π*激发态时，需要选出π特征的轨道。图中轨道17、20、21、22、23、24为π轨道。因此需要将17和19号轨道交换，构成CAS(6,6)的活性空间。

(3) 使用相同的基组做CAS计算，用以检查轨道交换是否成功：

%oldchk=benzene_hf.chk
%chk=benzene_cas.chk
#p cas(6,6)/sto-3g guess=(read,alter) geom=allcheck
 
19,17
 

将需要交换的轨道的编号写在输入文件末尾，有多组需要交换时，每一组写一行。计算完成后，观看轨道可以看到19号轨道为π轨道（此处略去该图）。

(4) 若要用CAS计算基态能量，输入文件为：

%chk=benzene_cas.chk
#p cas(6,6)/6-311G** guess=read geom=allcheck
 

最终能量为EIGENVALUE后面的数值。注意，在Exploring第三版的第384页，说读取结果是在TOTAL一行，是错误的。

(5) 若要计算激发态的能量：

%chk=benzene_cas.chk
#p cas(6,6,nroot=2)/6-311G** guess=read geom=allcheck
 

其中nroot=j表示计算CI的第j个根，其中第1个为基态，因此要算第一激发态，则需要将nroot设为2。这一步最后会给出两个EIGENVALUE，分别对应基态和激发态的能量，同时在第2个EIGENVALUE后面还会给出以eV为单位的垂直激发能。在这个计算中，轨道是在第一激发态下优化，对第一激发态是最优的，但不是对基态最优的，因此它的基态能量是不准的，激发能的误差偏大。

另一种获得垂直激发能的方法是将第(5)步的第一激发态的能量与第(4)步的能量相减。这种算法相当于是做了两个特定态计算，一个是基态（轨道在基态下优化，是对基态最优的），另一个是第一激发态（轨道在第一激发态下优化，是对第一激发态最优的），所得的激发能更准确一些。

若想获得更准确的激发能，还可以使用态平均（state-average, SA）CASSCF方法。在SA-CASSCF中，需要设定各个态的权重，一般将各个态的权重设为相等。此外，为得到更精确的结果，nroot的值可比感兴趣的态多2或3，类似TDDFT中nstates的设置。权重以F10.8（Fortran中的标准输入输出格式，表示整个数占10个字符，其中小数部分占8个字符）的格式依次写在输入文件最后，例如本例设nroot=4，每个态的权重相等，均为0.25，输入文件写成：

%chk=benzene_cas.chk
#p cas(6,6,nroot=4,stateaverage)/6-311G** guess=read geom=allcheck
 
0.250000000.250000000.250000000.25000000
 

同样在输出文件中会给出各个态的能量及相应的垂直激发能（位于EIGENVALUE后面）。

致谢

感谢jxzou在成文中的建议与帮助。