《量子化学软件基础》习题（7）
题目：
阅读文献《Full configuration interaction potential energy curves for breaking bonds to hydrogen: An assessment of single-reference correlation methods》(REF1)，重复文献中的计算结果。
解答：
在Born-Oppenheimer近似下，full configuration interaction (FCI)是选定单电子基组后多电子薛定谔方程的精确解。在REF1中，David Sherrill及其合作者考察了各种单参考态电子相关方法相对于FCI的精度。我们在重复过程中使用了与文献相同的基组，FCI结果来自于文献，iCISCF (iteration configuration interaction SCF)的计算使用BDF完成。由于很多计算涉及到Cartesian型基函数（6D, 10F），ORCA和BDF软件不支持，因此文章中FCI以外的计算均由Gaussian完成。
以BH分子的UCCSD计算为例，使用Gaussian进行键长扫描的输入文件为
%nprocshared=8
%oldchk=uhf.chk
%chk=001.chk
#P UCCSD/aug-cc-pVQZ guess=read scan nosymm
UHF scan
0 1
B
H 1 B1
B1 6.0 53 -0.1
在Gaussian中可以使用“scan”关键词进行键长扫描。在扫描时，可以先从解离曲线中H原子和B原子距离最远、BH键完全断裂的一端开始。这是因为此类情况下，体系的基态对应对称破缺的波函数，即开壳层单重态。在Gaussian中，这类计算可使用“guess=mix”关键词打破初猜波函数的对称性，确保参考态波函数不再是RHF波函数。只使用“guess=mix”有时也未必能收敛到能量最低的基态波函数，因此可以加上“stable=opt”来对收敛波函数的稳定性进行检测、优化。Gaussian会自动优化到能量最低的基态波函数。在进行UHF解离曲线的扫描时，可读取前一个构型的稳定波函数作为初猜。本文所有UHF计算均读取前一步扫面的键长的波函数作为初猜进行UHF计算，并使用“stable=opt”关键词。
【小编注：上述是常用的、较为保险的做法，比键长从短扫到长靠谱得多。但更保险的做法是不采用扫描，每个结构单独算单点、检验波函数稳定性。有的分子在不同键长下需要用不同形式的片段组合波函数构建初猜，见文末“相关阅读”】
对BH体系在aug-cc-pVQZ水平下进行RHF和基于RHF的单参考方法的计算结果的比较如图1所示：
图1 RHF及基于RHF的近似相关方法得到的BH分子键解离势能曲线图
从图中可以看出：
(1) RHF的势能曲线与FCI的解离曲线偏差较大，在键长趋于无穷大时出现了错误的渐近行为，这是因为两个原子在相距无穷远时RHF无法定性正确地描述该体系。
(2) 基于RHF的MP2方法在键长趋于无穷大时也存在错误的渐近行为。从数值角度来看，这是由于RHF参考波函数较差导致的。由于CCSD波函数的singles可以在某种程度上矫正参考态波函数，因此其结果要好于RHF-MP2。一般来说，对于单键解离曲线，RCCSD方法还是较为准确的。而RHF-CCSD(T)的结果也存在错误的渐进行为。
BH在aug-cc-pVQZ水平下进行UHF和基于UHF的单参考方法的计算结果的比较如图2所示，可以看到所有计算结果的渐近行为都是合理的，即UHF和基于UHF的单参考方法可以定性正确地描述BH的H键解离过程，但UHF和UMP2的计算得到的势能曲线的绝对能量和FCI相差较大。而UCCSD和UCCSD(T)计算得到的势能曲线与FCI结果非常接近，说明UCCSD和UCCSD(T)是可以准确描述BH分子的键解离，最终波函数的自旋污染也较小。
图2 UHF及基于UHF的近似相关方法得到的BH分子键解离势能曲线图
我们对文献中另外两个分子HF和CH4也进行了计算。对比文献中的数据，发现除个别数据点外，几乎所有的计算结果均与文献一致，这里就不再单独给出。与文献结果相差大于0.1 milli-Hartree (mEh)的数据见附录。
在HF和CH4分子的计算中，文献中所使用的6-31G**和6-31G*基组均为Cartesian型基函数，即，使用x2、y2、z2、xy、xz和yz等6个基函数来表示d轨道（6d）。对于变分计算，由于6d相对于spherical型基函数（5d）引入了额外的自由度，因此6d的能量结果一般比5d更低。注意，即使是Pople开发的6-31类型基组，BDF、ORCA等量子化学程序也使用5d基函数。
【小编注：本文研究的CH4分子仅解离一根C-H键，其余3根C-H键固定在1.086 Å】
以CH4分子为例，图3、图4分别是CH4分子使用基于RHF和UHF的相关方法在6-31G*水平下5d与6d计算结果之差。可以看出5d基组的结果略高于6d基组的结果。Hartree-Fock方法能量相差较小，MP2、CCSD和CCSD(T)方法能量相差可超过1 mEh。
图3 CH4分子基于RHF/6-31G*的键解离过程5d与6d基函数的能量差曲线图
图4 CH4分子基于UHF/6-31G*的键解离过程5d与6d基函数的能量差曲线图
由于轨道数较多，我们无法使用BDF或者ORCA重复文章中的FCI计算。因此，我们使用了BDF中的iCISCF，获得FCI的近似结果。iCI方法是刘文剑等人发展的一种选择性CI方法。在计算中，iCI方法会舍弃掉贡献较小的CSF，最终得到FCI的近似解。采用合理的截断参数，iCI和FCI方法计算出来的解离曲线只有非常微小的差别。iCIPT2是在iCI波函数的基础上进一步做ENPT2校正。iCISCF可以在iCI的波函数的基础上考虑活性空间内部轨道优化，得到更加接近FCI的结果。iCISCF(2)则是在iCISCF的基础上，在活性空间内部做ENPT2校正。以BH分子0.8 Å键长为例，使用BDF进行iCISCF计算的输入文件为
$compass
title
BH
basis
aug-cc-pVQZ
geometry
B 0.0 0.0 0.0
H 0.0 0.0 0.8
end geometry
group
C(2v)
norotate
saorb
$end
$xuanyuan
$end
$mcscf
core
1 0 0 0
active
49 16 30 30
close
0 0 0 0
actel
4
spin
1
symmetry
1
guess
hforb
nograd
ici
cmin
0.00001
actopt
2
$end为与文献中结果比较，在计算过程中要冻结B的1s轨道，并使用Jacobi rotation方法优化活性空间内的活性轨道。
iCISCF和iCIPT2的计算结果与文献中所给的FCI结果之差如图5所示。可以看到BH分子键解离过程中iCISCF计算结果（使用10-5截断波函数）与文献所给的FCI计算结果非常接近。这表明使用iCISCF或iCIPT2，即使不使用外推方法，也可以得到非常接近FCI的结果。对于HF和CH4分子，iCISCF和iCISCF(2)的结果与文献所给的FCI结果的差异，应该也主要来自于5d和6d基函数的差异。
图5 iCISCF和iCIPT2与FCI结果比较的能量差曲线图
附录
表1 BH分子相关计算中与文献结果相差较大数据 (in Hartree)
表2 HF分子相关计算中与文献结果相差较大数据 (in Hartree)
【小编注：这个差异是有原因的，这两个UHF解大家可以用Gaussian软件重复出来。前者波函数稳定，后者波函数不稳定。后者可以采用片段组合波函数为初猜得到】
表3 CH4分子相关计算中与文献结果相差较大数据 (in Hartree)
