《量子化学软件基础》习题（8）

题目：

1. 阅读文献Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3（https://doi.org/10.1063/1.2889385，下称REF1），使用ORCA重复文章中胞嘧啶（cytosine）的CASSCF及CASPT2计算（使用和文章相同的活性空间，CASPT2计算要注意shift的使用，必须使用对称性）。

2. 分别使用BDF和ORCA计算REF1中胞嘧啶所有态的FIC-NEVPT2能量，比较结果。

3. 分别使用BDF和ORCA计算REF1中胞嘧啶所有态的FIC-MRCISD能量，比较结果。

解答：

1. 首先分析胞嘧啶的对称性。REF1的SI中给出的胞嘧啶结构属于Cs点群，有不可约表示A’和A’’。REF1中胞嘧啶的CASPT2计算采用了两种不同的活性空间：a.活性空间包含10个电子和8个π轨道（在Cs点群下，8个π轨道的不可约表示都是A’’）；b.活性空间包含14个电子和10个空间轨道（除a中的8个π轨道外，又加入了两个属于点群Cs的A’不可约表示的孤对轨道n。ORCA的HF轨道中，孤对轨道的序号为25、26，在BDF中，孤对轨道是属于A’不可约表示中能量最高的两个轨道）。

小编注：此处8个π轨道指的是5个双占据π轨道和3个空π轨道，对应10个活性电子。后面的10个空间轨道指5个双占据π轨道，2个双占据轨道（作者认为是孤对电子）和3个空轨道，对应14个活性电子。这些轨道不一定是离HOMO、LUMO最近的轨道，需要肉眼观看和挑选。另外要注意，若用Gaussian程序开点群对称性计算，这些轨道的序号可能与ORCA的结果不同，下文的轨道序号不能机械照搬。

我们将以CASSCF(10,8)的计算为例说明输入文件的写法，CASSCF(14,10)的输入文件类似。若对具体输入文件的写法有疑问，可参考本文最后的GitHub链接，内有相应的输入及输出文件。使用ORCA进行CASSCF(10,8)计算的输入文件如下（其RHF结果已预先算好，波函数文件为cytosine-RHF.gbw）：

! TZVP Usesym MOread
%moinp "cytosine-RHF.gbw"
%pal
  nprocs   16
end
%scf
  rotate
   {20,25} {23,26} {30,33} {31,38}
  end
end
%casscf
 nel      10
 norb   8
 irrep  0
 mult   1
 NRoots 6
end
# cytosine CASSCF
* xyzfile 0 1 cytosine.xyz

BDF进行CASSCF计算的输入文件较为简洁，不需要过多考虑轨道旋转。

$compass
Title
 cytosine
basis
 TZVP
Geometry
file=cytosine.xyz
end Geometry
check
SAORB
$end

$xuanyuan
$end

$scf
rhf
$end

%cp cytosine-RHF.scforb cytosine-CASSCF-10-8.inporb

$MCSCF
guess
 read
spin
 1
close
 24 0
actel
 10
active
 0 8
symmetry
 1
Roots
 6 6 1
iprtmo
 1
molden
$end

使用ORCA和BDF计算胞嘧啶的两种活性空间态平均（SA，state-averaged）CASSCF计算的结果见表1。可看到，ORCA和BDF的结果在误差允许的范围内一致，故可用此结果进行CASPT2计算。

小编注：原文献对不同对称性的电子态（A',A''）使用了不同的活性空间，本文仅为重复性、练习性计算，不探究其他活性空间。

使用ORCA进行CASPT2(10,8) 的计算需改动上述输入文件中的casscf模块，此处以应用自然轨道（参数ActOrbs设置选项为NatOrbs，此为缺省设置）做CASPT2计算为例，可改为（下文对10、11两行有额外叙述）：

%casscf
 nel          10
 norb        8
 irrep       0
 mult       1
 NRoots     6
 Maxiter    200
 PTMethod   FIC_CASPT2
 PTsettings 
  CASPT2_IPEAshift    0.25
  CASPT2_rshift       0.10
 end
end

注意，REF1中说其设置level shift=0.10（其SI第18页），这一说法实际并不明确，因为ORCA和Molcas均可设置三种不同的shift：IPEAshift、rshift和ishift（Molcas的相应参数分别名为IPEAshift、IMAGinary和SHIFt）。level shift具体指代不明。经测试，REF1中的shift设置实际为IPEAshift=0.25、rshift=0.10。之所以存在该问题，猜测是原作者使用MOLCAS计算，其默认设定IPEAshift=0.25，故再设置的shift参数实际上是rshift，并需使用正则轨道（即在ORCA的casscf模块输入参数ActOrbs并指定选项CanonOrbs）。至于为何使用IPEA shift并设置为0.25，合理性可参考文献Assessment of n‑Electron Valence State Perturbation Theory for Vertical Excitation Energies（下载地址为https://pubs.acs.org/doi/pdf/10.1021/ct400136y），其3.6节对此进行了一定讨论，结论为对有机闭壳层分子设置IPEAshift=0.25会显著改善CASPT2的结果。

REF1的结果为MS-CASPT2的结果，而最新的ORCA 5.0版本尚无此方法。故此处采用FIC-CASPT2计算，但是FIC-CASPT2和MS-CASPT2的结果会略有差别。要验证其数据的真实性，需使用Molcas重复上述计算，计算的输入文件为

&GATEWAY
Coord
13
Acroleincoordinates in  Angstrom
H -2.114860 -1.429678 0.000000
H -0.173973 -2.806186 0.000000
H  2.073228 -1.658021 0.000000
H  3.175240 0.564335 0.000000
H  2.235202 2.033636 0.000000
C -0.060783 -1.726152 0.000000
C  1.144884 -1.099470 0.000000
C  1.107049 0.338190 0.000000
C -1.227573 0.430359 0.000000
O -2.315109 0.998271 0.000000
N  0.000000 1.058130 0.000000
N -1.201178 -0.989148 0.000000
N  2.278974 1.024187 0.000000
basis=TZVP
&SEWARD
&SCF
&RASSCF
nactel= 10  0  0
inactive=24 0
ras2=  0 8
Symmetry=1
spin=1
ciroot=  6 6 1
&CASPT2
FROZEN=0 0 0 0
multistate=6 1 2  3  4  5 6
maxiter=40
IPEAshift=0.0

基于(10,8)、(14,10)两种不同的活性空间、在计算中使用正则轨道或自然轨道，可分为4种情况。我们使用ORCA进行4种FIC-CASPT2计算（基于态平均CASSCF轨道），能量结果见表2。可以看到，本例中不同的CASSCF轨道处理方法会略微影响计算结果。这是由于shift的引入破坏了FIC-CASPT2的酉不变性。

注：如不设置任何shift，FIC-CASPT2能量对活性轨道的旋转是酉不变的。

2. 分别使用BDF和ORCA计算文献1中胞嘧啶所有态的FIC-NEVPT2能量，比较结果。

仍以处理活性空间CASSCF(10,8) 为例说明输入文件写法。使用ORCA计算胞嘧啶所有态的FIC-NEVPT2能量只需改动上述输入文件中的casscf模块，

%casscf
nel      10
norb   8
irrep  0
mult   1
NRoots 6
PTMethod  FIC_NEVPT2
end

而使用BDF则需要在CASSCF计算后再加入TraInt模块和Xi’an-CI模块，加入内容为

$TRAINT
Orbital
MCorb
MRPT2
$end
$XIANCI
nroot
6
spin
1
symmetry
1
NEVPT2
$end

使用ORCA和BDF计算胞嘧啶所有态的FIC-NEVPT2能量（均未应用冻芯近似）的结果见表3。需要注意，不论是ORCA还是BDF，对于不同对称性的态，NEVPT2是分块计算的。

注意：

① BDF需要记录的是SS-NEVPT2而非MS-NEVPT2的结果，此处SS即specific state，MS即multi-state。

② 为使ORCA和BDF的NEVPT2计算结果能够对照，注意正确设置冻芯近似：ORCA是默认设定冻芯近似而BDF不默认设置，若都想设置冻芯近似，可在后者输入文件的TraInt模块里加入关键字Frozen并设置参数，如在NEVPT2(10,8) 的计算中此选项需设置为“8 0”（即冻结属于不可约表示A’的8个轨道，它们都是碳、氮、氧的1s2电子）；若都不想设置冻芯近似，则ORCA输入文件需设置关键字NoFrozenCore，而BDF无需为此设置参数。

③ 使用BDF和ORCA计算得到的单点能结果略有差别。这主要是因为ORCA中FIC-NEVPT2是分别对每个态，重新构造每个态的Fock矩阵，再将其对角化得到轨道能和相应正则轨道（即canonstep采取选项1）。而BDF中的NEVPT2算法，使用的是SA-CASSCF的轨道和轨道能量（即其PTSettings的参数canonstep采取选项2）。目前ORCA中的FIC-NEVPT2尚不支持这种canonstep的设置。

3. 分别使用BDF和ORCA中的FIC-MRCI方法重复上述计算，比较结果。

FIC-MRCI（Fully-internally-contracted Multi-reference Configuration Interaction）是一种MRCISD方法。在此简要叙述不收缩（uncontracted，简记为uc）、内收缩（internally-contracted，简记为ic）的概念。设由MCSCF计算得到的体系波函数|CAS〉为

其中，cI是组态|I〉在态|CAS〉中的相应系数。

uc-MRCI即将所有参与构成的|CAS〉组态|I〉及其单双激发态共同张成的行列式空间作为变分空间。设体系的占据轨道数为no，非占轨道数为nv，组态|I〉的总数为nI，则uc-MRCI计算涉及的组态数的数量级约为no2nv2nI。

ic-MRCI，即将|CAS〉及其单双激发态共同张成的线性空间作为变分空间。由于活性空间的轨道数目一般远小于占据轨道数，故ic-MRCI计算涉及的单双激发空间维度和单参考态CISD方法接近，约为no2nv2。单双激发按照电子激发的性质可以分为15种类型，若不区分双激发中的单三激发类型，则为8种类型。FIC-MRCI将所有类型都进行内收缩，所以其名称中有Fully字样。不同内收缩类型的总结可以参考（The Journal of Chemical Physics 154 (21), 214111）

使用BDF进行FIC-MRCI计算时，参考组态空间分成两部分，第一部分是参考组态，第二部分是参考组态的正交空间。默认情况下，这两部分都会进入FIC-MRCI的变分空间。若舍弃第二部分，不考虑参考态的弛豫，此即为BDF中的VVCI关键词对应的方法。使用VVCI后，BDF和ORCA的FIC-MRCI方法在处理MCSCF的单根时结果完全相同，但同时处理两根或更多根时结果仍然会略有区别。这是因为两款软件的底层实现略有不同，ORCA的FIC-MRCI会将每一个根分别构造波函数并对角化，不计入这些根之间的相互作用；但BDF的VVCI会将相同对称性的根一起处理，考虑了这些根之间的相互作用。对此可见表11-4（下文BDF的结果使用VVCI），在活性空间CASSCF(14,10)中，ORCA和BDF计算得到基态的结果一致，但激发态稍有区别。 仍以活性空间CASSCF(10,8) 为例说明输入文件写法。使用ORCA计算胞嘧啶所有态的FIC-MRCI能量需改动上述输入文件中的casscf模块并增加autoCI模块，可改为

%casscf
nel    10
norb   8
mult   1
irrep  0
Nroots 6
end
%autoCI
CItype FICMRCI
nel      10
norb   8
irrep  0
mult   1
NRoots 6
end

而使用BDF则需要在CASSCF计算后再加入TraInt模块和Xi’an-CI模块，加入内容为

$TRAINT
Frozen
 8 0
Orbital
 MCorb
$end
$XIANCI
nroot
 6
spin
 1
symmetry
 1
VVCI
$end

ORCA和BDF的FIC-MRCI计算结果见表4（均使用冻芯近似）。可看到除11A’->61A’过程外，两者的激发能结果基本一致，11A’->61A’过程的激发能结果相差稍大。

此外，我们还用BDF比较了冻芯近似对其单点能的影响，结果见表5。可看到，不使用冻芯近似降低了电子态的绝对能量，但是对激发能的影响非常小。

附录：上述计算的所涉及文件均可见链接

https://github.com/liu-chunzhang/an_introductory_tutorial_of_quantum_chemistry/tree/main/exe11

此仓库里还有许多个人的其他计算及相应说明文件可供读者学习和参考。

山东大学前沿交叉科学青岛研究院 刘纯彰

于2022年6月19日

附录：计算采用的胞嘧啶的分子坐标为

H -2.114860           -1.429678           0.000000
H -0.173973          -2.806186           0.000000
H 2.073228            -1.658021           0.000000
H 3.175240            0.564335             0.000000
H 2.235202            2.033636             0.000000
C -0.060783          -1.726152           0.000000
C 1.144884            -1.099470           0.000000
C 1.107049            0.338190             0.000000
C -1.227573          0.430359             0.000000
O -2.315109          0.998271             0.000000
N 0.000000            1.058130             0.000000
N -1.201178           -0.989148           0.000000
N 2.278974            1.024187             0.000000