《量子化学软件基础》习题 （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 是指两个位点处的未配对电子数目，输入文件如下：

BrokenSym计算会先进行高自旋（三重态）计算，再进行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关键词的输入文件如下：

在BDF中可以利用 FLMO (Fragment localized molecular orbital) 手动分片计算开壳层单重态。具体输入文件如下：

BDF的FLMO方法是一种基于分片的方法，将一个体系分成若干个片段，在每个片段的SCF收敛后，程序对每个片段的波函数进行局域化，再用所得的局域轨道产生总体系计算的初猜。对于正庚烷双自由基的开壳层单重态，将单电子所在的两个C原子分在不同的两个片段。对分子分片后，分别指定每个片段的自旋多重度，对于该分子，每个片段自旋多重度为2，总的自旋多重度为1，通过计算可以得到开壳层单重态。

本例的具体做法是把正庚烷双自由基分子手动分成两片，每个分子片是由中心原子、缓冲区原子（B）和链接 H 原子（L）组成。两个分子片的缓冲区原子都为与其末端C原子直接相连的C原子及其所带的氢原子组成。即第一个分子片为正戊烷自由基，第二个分子片为正丁烷自由基。每个分子片都需要进行UHF及 轨道局域化。然后，在轨道局域化模块localmo后插入Shell命令

将存储定域轨道的文件$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模块用于对分子自动分片，并产生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计算完全相同。输入文件如下：

活性轨道28，29的图像及其自然轨道占据数如下：

图3 正庚烷双自由基两个SOMO轨道

从图中可以看出，两个SOMO轨道是混合的，如果想要得到两个不混合的SOMO轨道，在ORCA中，可以在casscf 模块使用actorbs locorbs关键词，对活性轨道做局域化，从而得到以下轨道：

图4 正庚烷双自由基两个SOMO轨道（局域轨道）

使用BDF计算，同样先进行ROHF/cc-pVDZ计算，CASSCF读ROHF计算结果做初猜进行计算，两个SOMO轨道是28，29号轨道，轨道图像及其自然轨道占据数信息不再详细列出，与ORCA计算结果同。具体输入文件如下：

在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的输入文件都和三重态的计算类似：

使用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)

附录

优化后正庚烷双自由基，单位：Å