《量子化学软件基础》习题 （3）
习 题
1. 构造一个三重态双自由基分子，使用UHF对该双自由基分子进行结构优化。通过自旋布局(Spin Population)确定两个单占轨道(singly occupied molecular orbitals, SOMOs) 所在的原子。使用Broken Symmetry方法计算该双自由基的“开壳层单重态”。
2. 使用CASSCF(2,2)研究上述分子的三重态、开壳层单重态以及闭壳层单重态。
解答：
1. (1) 构造正庚烷双自由基分子（C7H14），去掉了正庚烷两端碳上的氢原子（C1和C19上的H），使用ORCA在高自旋（自旋多重度为3）UHF/cc-pVDZ水平下进行结构优化，优化后的坐标见附录1。
图1 C7H14双自由基分子优化后结构 (本文轨道图像均使用Multiwfn软件绘制)
从输出文件的spin population中可以看出，单电子主要分布在C1和C19（图2 中 0 C和18 C）上。即，C1和C19是两个SOMO轨道在的原子。
图2 C7H14的spin population
(2) 在ORCA中，可以使用BrokenSym或者Flipspin关键词进行Broken symmetry计算。在scf 模块中使用BrokenSym关键词，1, 1 是指两个位点处的未配对电子数目，输入文件如下：
!UHF cc-pVDZ pal2
%scf
BrokenSym 1,1
end
*xyzfile 0 3 c7h14.xyzBrokenSym计算会先进行高自旋（三重态）计算，再进行Broken Symmetry计算得到单重态。因此，这里的自旋多重度必须是3【小编注：小编在ORCA 4.2.1和5.0.3版本上都试过，BrokenSym 1,1优先级高于*xyz处的自旋多重度，因此这里的自旋多重度写1和3都不影响计算。但为了易于阅读，建议写3】。输出文件如下：
从输出文件中可以看出， 0 C和18 C的spin population绝对值均接近1，且符号是相反的。此外，输出文件中轨道占据信息显示alpha和beta轨道的占据数相等，均含有28个电子，从而可以确定计算得到的是开壳层单重态。
在ORCA中还可以使用FlipSpin关键词来进行Broken Symmetry的计算。FlipSpin关键词后为要翻转自旋的原子序号（注意ORCA是从0 开始），FinalMs为翻转后体系的总磁量子数。使用FlipSpin关键词的输入文件如下：
!UHF cc-pVDZ pal2
%scf
FlipSpin 18
FinalMs 0
end
*xyzfile 0 3 c7h14.xyz在BDF中可以利用 FLMO (Fragment localized molecular orbital) 手动分片计算开壳层单重态。具体输入文件如下：
%echo " -------------- CHECKDATA: Calculate the 1st fragment ------------------ "
$compass
title
c7h14
basis
cc-pvdz
geometry
C -0.027056 0.014499 0.007320
H 0.334747 0.551751 -0.877315
H 0.305045 0.594495 0.874341
C 0.609303 -1.374514 0.044920
H 0.262702 -1.951370 -0.819054
H 0.259278 -1.912060 0.933203
C 2.147820 -1.364760 0.037559
H 2.493149 -2.399802 -0.057560
H 2.498335 -0.839764 -0.856088
C 2.771544 -0.752032 1.258756
H 3.140041 0.266301 1.257207
H 2.687762 -1.251262 2.218097
C -1.554416 -0.023401 -0.016100 B
H -1.919914 -0.559214 0.866609 B
H -1.893486 -0.598245 -0.884716 B
H -2.193484 1.373788 -0.056811 L
end geometry
$end
$xuanyuan
$end
$scf
UHF
spin
2
$end
$localmo
FLMO
Pipek
$end
%cp $BDF_WORKDIR/$BDFTASK.flmo $BDF_TMPDIR/fragment1
%cp $BDF_WORKDIR/$BDFTASK.flmo $BDF_WORKDIR/fragment1
%ls -l $BDF_TMPDIR
%rm -rf $BDF_TMPDIR/$BDFTASK.*
%echo " ----------------- CHECKDATA: Calculate the 2nd fragment -------------------"
$compass
title
c7h14
basis
cc-pvdz
geometry
C -3.693788 1.345952 -0.080154
C -2.193484 1.373788 -0.056811
H -4.225507 1.180090 -1.010335
H -4.255203 1.224047 0.839301
C -1.554416 -0.023401 -0.016100
H -1.824925 1.905707 -0.939849
H -1.852811 1.944020 0.813356
H -1.919914 -0.559214 0.866609
H -1.893486 -0.598245 -0.884716
C -0.027056 0.014499 0.007320 B
H 0.334747 0.551751 -0.877315 B
H 0.305045 0.594495 0.874341 B
H 0.609303 -1.374514 0.044920 L
end geometry
$end
$xuanyuan
$end
$scf
UHF
spin
2
$end
$localmo
FLMO
Pipek
$end
%cp $BDF_WORKDIR/$BDFTASK.flmo $BDF_TMPDIR/fragment2
%cp $BDF_WORKDIR/$BDFTASK.flmo $BDF_WORKDIR/fragment2
%ls -l $BDF_TMPDIR
%rm -rf $BDF_TMPDIR/$BDFTASK.*
%echo " ---------- CHECKDATA: From fragment to molecular SCF calculation -------------"
$compass
title
c7h14
basis
cc-pvdz
geometry
file=c7h14.xyz
end geometry
Nfragment
2
$end
$xuanyuan
$end
$scf
UHF
spin
1
FLMO
molden
$end
&DATABASE
Fragment 1 15
8 11 12 13 14 15 16 17 18 19 20 21 5 9 10
Fragment 2 12
1 2 3 4 5 6 7 9 10 8 11 12
&ENDBDF的FLMO方法是一种基于分片的方法，将一个体系分成若干个片段，在每个片段的SCF收敛后，程序对每个片段的波函数进行局域化，再用所得的局域轨道产生总体系计算的初猜。对于正庚烷双自由基的开壳层单重态，将单电子所在的两个C原子分在不同的两个片段。对分子分片后，分别指定每个片段的自旋多重度，对于该分子，每个片段自旋多重度为2，总的自旋多重度为1，通过计算可以得到开壳层单重态。
本例的具体做法是把正庚烷双自由基分子手动分成两片，每个分子片是由中心原子、缓冲区原子（B）和链接 H 原子（L）组成。两个分子片的缓冲区原子都为与其末端C原子直接相连的C原子及其所带的氢原子组成。即第一个分子片为正戊烷自由基，第二个分子片为正丁烷自由基。每个分子片都需要进行UHF及 轨道局域化。然后，在轨道局域化模块localmo后插入Shell命令
cp $BDF_WORKDIR/$BDFTASK.flmo $BDF_TMPDIR/fragment* 将存储定域轨道的文件$BDFTASK.flmo拷贝到$BDF_TMPDIR所在的目录备用。2个分子片段计算完成后，进行整体分子的计算。在compass中，有关键词Nfragment 2，提示要读入2个分子片，分子片信息在&DATABASE中定义，其中Fragment后为分子片序号及原子数目（不含链接H原子，输入中L标记），下一行为该分子片原子序号。从最终输出文件的自旋布居信息看出计算得到的为开壳层单重态，见下图：
手动分片的输入文件过于繁琐，因为需要自己手动指认每个分子片的自旋多重度以及电荷，还要在&DATABASE模块输入每个分子片的信息。相比基于FLMO的手动分片方法来说，结合FLMO的自动分片方法来进行Broken Symmetry计算更方便。可参考
https://bdf-manual.readthedocs.io/zh_CN/latest/User%20Guide.html#id35
具体输入文件如下：
$autofrag
method
flmo
spinocc
1 +1 19 -1
$end
$compass
title
c7h14
basis
cc-pVDZ
geometry
file=c7h14.xyz
end geometry
$end
$xuanyuan
$end
$scf
UHF
spin
1
$end
$localmo
FLMO
Pipek
$endautofrag模块用于对分子自动分片，并产生FLMO计算的基本输入。autofrag模块的主要流程是先根据分子坐标自动产生原子的键连信息，然后切断某些键，产生合适的分子片段，再对分子片段增加缓冲原子，在断键处添加连接氢原子，最后生成子体系分子片段的输入文件，将子体系分子片段的计算结果组织起来，运行整体分子计算。
BDF先根据compass 模块中的分子结构与指定的autofrag的参数信息产生分子片段，以及分子片段定域化轨道计算的输入文件。然后用分子片段的定域轨道组装整体分子的 pFLMO (primitive Fragment Local Molecular Orbital) 作为全局SCF计算的初始猜测轨道，再通过全局SCF计算，在保持每一步迭代轨道都是定域轨道的前提下，得到整体分子的开壳层单重态。在输入文件中，autofrag模块可通过spinocc设定具体原子的自旋, 这里是指定第1个原子有一个未成对的alpha电子，第19个原子有一个未成对的beta电子，注意spinocc在设定时，所有开壳层原子都应该被指定。若分子带电荷，使用charge关键词可设置具体原子的电荷，格式同spinocc的设置格式。另外有关键词radbuff可以设置分子片缓冲半径（非负浮点数），默认值为2.0。ORCA和BDF的计算结果见表1。
表1 C7H14双自由基开壳层单重态能量(in Hartree)
由表1可看出，ORCA和BDF计算的双自由基开壳层单重态的能量几乎是一致的。若追求更严格的一致，可在ORCA中加入VeryTightSCF关键词，会让ORCA的结果更接近BDF和Gaussian。
另外，Gaussian中片段组合波函数方法（详细做法请参考https://gaussian.com/afc/）与BDF的FLMO手动分片方法相似，先对分子进行分片，然后指定每个分子片的自旋多重度和电荷，Gaussian生成片段组合波函数的任务开始后，会初猜整体波函数并做分析，不进行迭代，然后对每个片段初猜波函数，进行SCF迭代然后做分析，最后把得到的片段波函数组合。把最终得到的组合后的片段波函数做初猜，进行开壳层单重态的计算。使用Gaussian计算得到结果为-273.158758 Hartree，与BDF和ORCA的结果一致。
2. (1) 使用CASSCF(2, 2)计算该双自由基的三重态
先进行ROHF/cc-pVDZ计算，用该计算得到的轨道可以作为CASSCF计算的初猜轨道。使用CASSCF(2,2)计算该双自由基的三重态，仅有一个组态，两个电子自旋相同，分占两个轨道（28，29轨道），此时的CASSCF计算相当于做ROHF计算，CASSCF和ROHF计算完全相同。输入文件如下：
!cc-pVDZ
!moread
%moinp “ROHF.gbw”
%casscf
mult 3
nel 2
norb 2
nroots 1
end
*xyzfile 0 3 c7h14.xyz活性轨道28，29的图像及其自然轨道占据数如下：
图3 正庚烷双自由基两个SOMO轨道
从图中可以看出，两个SOMO轨道是混合的，如果想要得到两个不混合的SOMO轨道，在ORCA中，可以在casscf 模块使用actorbs locorbs关键词，对活性轨道做局域化，从而得到以下轨道：
图4 正庚烷双自由基两个SOMO轨道（局域轨道）
使用BDF计算，同样先进行ROHF/cc-pVDZ计算，CASSCF读ROHF计算结果做初猜进行计算，两个SOMO轨道是28，29号轨道，轨道图像及其自然轨道占据数信息不再详细列出，与ORCA计算结果同。具体输入文件如下：
$compass
title
c7h14
basis
cc-pVDZ
geometry
file=c7h14.xyz
end geometry
saorb
$end
$xuanyuan
$end
$scf
ROHF
spin
3
molden
$end
$mcscf
molden
spin
3
close
27
active
2
actel
2
guess
hforb
roots
1 1 1
$end在scf和mcscf模块使用molden关键词，计算结束可以得到molden类型文件，guess关键词可以指定CASSCF的初猜轨道。这里hforb是读取scf模块产生的轨道做CASSCF的初猜。spin为自旋多重度，active为活性轨道数目，actel是活性空间里的电子数目，roots为算的根的数目（默认值为1），第一个是态平均数，第二个是CI 中计算的态的数目，如果第三个数字是1，以相同的权重取平均值。注：在BDF中可以使用localmo模块对题中两个SOMO轨道进行局域化。
使用ORCA和BDF计算的C7H14三重态的能量见表2。从结果可见，两个软件的计算结果一致。
表2 C7H14双自由基开壳层三重态的能量(in Hartree)
(2) 使用CASSCF(2,2)计算该双自由基的开壳层单重态
读取ROHF/cc-pVDZ三重态的轨道做初猜，BDF和ORCA的输入文件都和三重态的计算类似：
!cc-pVDZ pal8
!moread
%moinp “ROHF.gbw”
%casscf
mult 1
nel 2
norb 2
end
*xyzfile 0 3 c7h14.xyz使用ORCA进行CASSCF计算时，分子结构部分的自旋多重度不会对所计算的态造成影响。计算结果见表3：
表3 C7H14双自由基开壳层单重态的能量(in Hartree)
两个SOMO轨道图像和自然轨道占据数见图5，
图5 开壳层单重态下两个SOMO轨道
注意，在判断计算得到的是否为开壳层单重态时，不能只通过占据数来判断，要结合轨道图像一起判断。ORCA做CAS默认得到的是自然活性轨道，不一定能够直接通过组态信息确定所得到的电子态，见下图：
因为CASSCF的能量对活性轨道的旋转是酉不变的，因此，如果要更直观地展现单电子，可以使用actorbs locorbs关键词得到局域化的活性轨道。使用局域活性轨道，易于通过组态信息确定电子态。
局域活性轨道如下图所示：
图6 开壳层单重态下两个SOMO轨道(actorbs locorbs)
(3) 使用CASSCF(2,2)计算该双自由基的闭壳层单重态
首先可以进行RHF/cc-pVDZ计算，读RHF结果做初猜进行CAS(2, 2)计算。由于开壳层单重态的能量要低于两个闭壳层单重态，如果仅计算基态，得不到闭壳层单重态。使用CASSCF计算多个电子态的时候往往使用态平均（state-average），即各个态具有一定的权重，通过优化平均加权CASSCF能量，得到最终的态平均分子轨道。
根据化学知识可知，两个闭壳层单重态是近简并的，其能量高于开壳层单重态。因此，我们可以计算三个单重态的态平均CASSCF，每个态的权重都相等。在ORCA中把casscf模块nroots关键词设为3，在BDF中mcscf模块roots关键词设为3 3 1即可。ORCA和BDF的计算结果见表4，结果一致。同理，如果想看某个闭壳层行列式占绝对的比重，在这个例子中可以加上actorbs locorbs关键词，然后阅读输出文件中的组态系数。
表4 C7H14双自由基的单重态能量(in Hartree)
附录
优化后正庚烷双自由基，单位：Å
C -3.693788 1.345952 -0.080154
C -2.193484 1.373788 -0.056811
H -4.225507 1.180090 -1.010335
H -4.255203 1.224047 0.839301
C -1.554416 -0.023401 -0.016100
H -1.824925 1.905707 -0.939849
H -1.852811 1.944020 0.813356
C -0.027056 0.014499 0.007320
H -1.919914 -0.559214 0.866609
H -1.893486 -0.598245 -0.884716
H 0.334747 0.551751 -0.877315
H 0.305045 0.594495 0.874341
C 0.609303 -1.374514 0.044920
H 0.262702 -1.951370 -0.819054
H 0.259278 -1.912060 0.933203
C 2.147820 -1.364760 0.037559
H 2.493149 -2.399802 -0.057560
H 2.498335 -0.839764 -0.856088
C 2.771544 -0.752032 1.258756
H 3.140041 0.266301 1.257207
H 2.687762 -1.251262 2.218097