使用PySCF程序做CASSCF和MC-PDFT计算（更新中）

使用PySCF程序做CASSCF和MC-PDFT计算
CASSCF and MC-PDFT Calculations Using PySCF Program

Jiawei Xu
Released: 2024-08-09 / Updated: 2024-08-30

MC-PDFT (multi-configuration pair density functional theory)是最近几年逐渐兴起的CASPT2和NEVPT2的平替方法，耗时相比后两者明显降低，但精度却损失不明显，目前主要被OpenMolcas、PySCF-forge和GAMESS等程序所支持。本文讲解PySCF程序结合PySCF-forge扩展包做MC-PDFT计算的基本方法，也同时讲解PySCF做CASSCF计算的基本方法，部分流程会使用到MOKIT程序作为协助。

1. 离线安装PySCF和PySCF-forge扩展包

对于有root权限、在线安装的情况：

pip install pyscf
pip install pyscf-forge

对于最一般的、无root权限、离线安装的情况，推荐使用Anaconda Python。首先在一台能联网的机子上安装与集群相近版本的Anaconda Python，并下载PySCF及其必要依赖的安装包：

pip download pyscf
pip download pyscf-forge

上述命令将会下载得到一系列的.whl安装包，将之上传至集群。我们选择使用Anaconda Python，由于PySCF是比较常用的库，可以直接安装在base环境。首先激活base环境，然后使用pip安装：

source ~/programs/Anaconda3/bin/activate base
which pip # 确认pip路径，此时应对应于~/programs/Anaconda3/bin/pip，默认安装至~/programs/Anaconda3/lib/python3.X/site-packages
pip install pyscf-2.xxx.whl

2. 简单CASSCF计算

2.1 苯
苯分子基态的活性空间一般选为(6o,6e)，即包含全部6个π轨道和6个π电子。该计算可直接使用MOKIT中的automr程序全自动完成，输入文件如下（原子坐标略）：

%nprocshared=28
%mem=84GB
#p CASSCF/def2SVP

mokit{}

0 1
[Cartesian Coordinates]

注意：def2-SVP基组对于CASSCF计算而言太粗糙，此处仅作为演示时使用，该任务不到一分钟就能迅速算完。
苯分子比较简单，所以使用默认的设置即可。MOKIT会首先对所提供结构进行RHF和UHF计算，并根据能量判断体系为闭壳层或开壳层。对于闭壳层的情况，会首先对虚轨道进行投影，然后做定域化，并匹配成键和反键轨道。对于开壳层的情况则是基于UHF波函数得到UNO (unrestricted natural orbital)，然后匹配成键和反键轨道。基于此，调用GAMESS（或Gaussian）做广义价键计算(generalized valence bond, GVB)，并将收敛的GVB波函数作适当处理选出活性空间，传递给PySCF进行CASSCF计算。本例MOKIT自动选出的活性空间即为(6o,6e)，计算完成后得到一系列的输出文件。

观察ph_uhf_gvb15_CASSCF_NO.fch中的轨道，可见活性空间内的NO均为非整数占据，总占据数为6，其等值面图像也符合苯环π轨道的图像：

我们简单讲解一些有关轨道自动选择的计算细节。相比于Full CI，CASSCF实际上是选择一部分最具化学意义的轨道和电子作为活性空间，在活性空间范围内做Full CI。与之类似，GVB实际上也是针对一部分成反键轨道对进行有限范围内的VB计算，因此GVB轨道对也可以看做是某种意义上的活性空间。
在由UNO向GVB传递初猜时，MOKIT默认选择占据数超过10^-5的空轨道，也即ph_uhf_uno_asrot.fch中的22~36号轨道，将其以及与其一一对应的7~21号占据轨道全部纳入GVB活性空间，共15对GVB轨道对，故后续文件命名使用”gvb15″。这个标准很严格，绝大多数情况下实际上已经包含了所有的价层轨道。这个阈值可以通过关键词修改，如uno_thres=0.001则会将此标准放宽至0.001。如果采用0.001的阈值，GVB计算将仅限于π轨道，即ph_uhf_uno_asrot.fch中的19~24号轨道，对于本例来说实际上影响不大，因为苯环价层轨道中的其他C-C C-H σ轨道本身多参考性就几乎可以忽略。此外，也可以通过excludeXH关键词直接去除GVB计算中的C-H键轨道。
类似地，在由GVB向CASSCF传递初猜时，也是通过占据数决定活性空间所包含的轨道范围。MOKIT默认将占据数在0.02~1.98范围内的轨道全部纳入考虑，但这一范围没有现成的关键词可以修改（至少对于撰文时是如此）。这就对应于ph_uhf_uno_asrot2gvb15_s.fch中的19~24号轨道，恰好是所有的π轨道。也可以在GVB完成后终止计算，并通过关键词on_thresh=0.01,readno=”ph_uhf_uno_asrot2gvb15_s.fch”,ist=5以NO的形式读取GVB波函数，并取0.01作为占据数偏离0和2的阈值。

2.2 苯酚
仅考虑苯的部分，活性空间仍为(6o,6e)，但考虑到氧原子的LP电子会与π键产生一定程度的共轭，选择(7o,8e)的活性空间似乎更妥。但由于MOKIT程序是成对地构建初猜波函数，因此实际上只能计算活性轨道数和活性电子数相等，即(no,ne)的情况，因此不可能选出(7o,8e)的活性空间。我们首先仿照2.1中做法使用MOKIT程序自动做CASSCF计算，可以看到所选出的活性空间仍为(6o,6e)，包含23~28号轨道，而氧原子的LP轨道为20号轨道，不在活性空间的范围内。这是由于在GVB计算结果中，氧原子LP轨道的占据数不在默认的0.02~1.98范围内，因此被排除出了活性空间。

我们当然可以读取GVB波函数，并通过on_thresh来调整活性空间的范围，但这会同时考虑不必要的其他一些轨道，增大计算量。这时候就需要人为地调整轨道顺序。首先我们要理解PySCF的输入文件中如何定义活性空间：mcscf.CASSCF(m,(na,nb))表示活性空间包含m个空间轨道，na个α电子和nb个β电子，此时活性空间实际上是包含了前线轨道附近的这部分轨道和电子。对于单重态的苯酚，mcscf.CASSCF(6,(3,3))就对应于HOMO-2、HOMO-1、HOMO、LUMO、LUMO+1和LUMO+2及3个α电子和3个β电子。因此，额外包含氧原子的LP轨道后，活性空间的定义变为mcscf.CASSCF(7,(4,4))，同时我们还需要将20号轨道与22号轨道调换顺序，行列式波函数的性质决定了调换轨道顺序仅改变符号，而不改变波函数其他性质。使用Multiwfn的主功能6调换顺序，并将其重新保存为exchanged.fch文件：

Multiwfn phenol_uhf_gvb16_CASSCF_NO.fch
6
29
20,22
-1
100
2
7
exchanged.fch
0
q

测试表明，Multiwfn生成的这个fch文件并不能被PySCF正确读取，确切地说是不能被PySCF输入文件中所调用的MOKIT中提供的fch2py函数正确读取。无妨，我们可以用Gaussian整理其文件格式，先转换为.chk格式：

unfchk exchanged.fch exchanged.chk

编写如下Gaussian输入文件，并运行之，该任务会迅速完成：

%chk=exchanged.chk
#p ROHF/def2SVP guess=(read,only) geom=allchk nosymm int=nobasistransform

注意上述任务需要保持基组与此前的计算均一致。得到新的.chk文件后，再转换为.fch格式即可，此时我们就得到了一个可以正确被py2fch读取的.fch文件作为初猜了。由于前面已经进行过活性空间(6o,6e)的CASSCF计算，因此实际上只需要基于此前的PySCF输入文件修改活性空间定义和初猜波函数文件名即可。执行该活性空间(7o,8e)的CASSCF计算，观察所得波函数就会发现，氧原子的LP轨道相较于纳入活性空间前发生了明显的变化，而其他π轨道也由于共轭都掺入了一定LP轨道的成分。

比较两次计算得到的CASSCF能量：CASSCF(6,6) -305.406703386517 a.u.；CASSCF(7,8) -305.407351500871 a.u. 可见扩大活性空间后能量有所降低，但由于LP轨道与π轨道的共轭程度较小，能量差异并不明显。一般来说，都要求活性空间的大小相同时才可以比较能量。也可以查看CASSCF(7,8)计算的第一轮迭代的能量，即初猜波函数CASCI的能量为E0 = -305.407324890267 a.u.，这就指明了同一波函数在不同活性空间下所得的能量是不可比较的。同时，E0与CASSCF(7,8)最终收敛的能量之间的差值也可以看做是由LP轨道与π轨道共轭所带来的能量降低。

使用MOKIT调换轨道顺序
MOKIT有一内置的permute_orb函数也可用于调换轨道。编写一个python文件：

from mokit.lib.gaussian import permute_orb
permute_orb('ethanol_rhf_proj_loc_pair2gvb8_s.fch',6,13)

运行即可调换指定文件中两个轨道的顺序，轨道编号从1开始。也可以将此直接写进PySCF做CASSCF计算的输入文件。

2.3 二茂铁
在PBE0-D3(BJ)/def2-SVP下优化得到二茂铁的极小点结构，并使用MOKIT做默认的CASSCF计算，将得到以下结果。这是由于在收敛所得的GVB波函数中并没有占据数介于0.02~1.98之间的轨道，因此没有轨道被选入活性空间。一般来说，这表明当前计算的体系多参考性较弱，使用CASSCF计算的意义不大。

即便如此，我们也依然可以使用CASSCF方法计算二茂铁的电子结构，最理想的活性空间应当包含Fe的全部5个3d轨道，以及两个环戊二烯负离子的10个π轨道，总共有18个电子占据，则活性空间的大小为(15o,18e)。首先通过观察GVB轨道的图像确认这15个轨道的序号，如下图所示。如何迅速选中这些轨道，这里就本例而言谈一种策略：
① 22~24号轨道是Fe的3d⁶轨道，从配体场理论(ligand field theory, LFT)的角度而言，这三个d轨道不与配体的轨道发生混合，因此表现为LP轨道，因而也易于从图像直接辨认。
② 36/46/48和35/37/47分别是上下两个环戊二烯负离子的π轨道。
③ 根据价键理论的思想，每个成键轨道都应对应于一个反键轨道，相应地就有49/51/61和50/60/62为其反键轨道，因为在MOKIT处理过的GVB波函数中，轨道按占据数排序，因而成反键轨道的序号恰好是关于HOMO-LUMO (48-49)对称分布的。但这就会产生12个π/π*轨道，这是不符合常理的。实际上这里的6个“π轨道”已经不是单纯的π轨道，而是非占据d轨道和π轨道相互作用产生的成键轨道，反键轨道似此，只不过部分轨道仍然很大程度上保留了π/π*轨道的图像，即只有36/61和35-62两组有明显的d-π图像性质。具体可参考配体场理论以及群论的内容。

值得一提的是，GVB波函数中环戊二烯负离子的π轨道并不是传统观点中“肩并肩”的图像（即类似于下图中的Pipek-Mezey LMOs），此处更接近于Boys LMOs的双键图像。Pipek-Mezey LMOs呈现的C-C双键的结果更符合化学直觉，Boys LMOs则将此处理成了两个等价的香蕉键。这两种双键的描述形式没法说谁对谁错，从物理角度上是等价的描述，只是恰巧对于本例来说，GVB波函数的结果与Boys的结果更接近。明白这一点，对于确定需要纳入活性空间的轨道序号很重要。事实上在双占据轨道中，每个环戊二烯负离子能够找到8个图像似此的轨道，显然，其中占据数更接近于2的5个轨道对应于C-C σ轨道，占据数偏离2较多的3个轨道为π轨道。

我们需要将这15个轨道通过调换顺序使之位于前线轨道，此时需要多次调换才能完成，此处提供一个一键调换轨道顺序的小程序exmo.py。运行之，得到重新排序后的波函数文件：

python exmo.py ferrocene_uhf_uno_asrot2gvb29_s.fch "22-24 35-37 46-48 49-51 60-62"

然后使用所得的exchanged.fch作为初猜，仿照2.1和2.2中做法使用Gaussian处理文件格式问题，然后进行CASSCF计算。先产生含有分子结构、基组等基本信息的PySCF输入文件：

bas_fch2py exchanged.fch

这将会得到一个新的名为exchanged.py的PySCF输入文件。仿照2.1和2.2中做CASSCF计算的输入文件，补充开头import以及CASSCF计算的部分代码即可提交任务。此外，本例所选的(15o,18e)活性空间虽然足够完善，但计算量较大，这样尺寸的活性空间一般可以做DMRG-SCF以节省计算量。