广义价键计算及初始轨道的构造

上次推送广义价键波函数的简介后，很多朋友在后台询问用什么软件可以做广义价键（GVB）计算，以及如何做GVB计算？为此，我们简单地整理了一下，今天推送给大家。

GVB计算的常用量子化学软件是Gaussian、GAMESS等。Gaussian是商业软件，最大能够进行50对的GVB计算（手册上这么说），而GAMESS是开源软件，最多可处理12对大小的体系。然而，GVB计算的真正困难在于有效地构造初始轨道，它并不是自动的，不如Hartree−Fock方便（只需提供分子结构就可计算），需要手动地选择轨道并一一配对（这一点它比多组态自洽场（MCSCF）或完全活性空间自洽场（CASSCF）更麻烦，不仅需要选择轨道还需要配对轨道）。下面以Gaussian软件对CH₄分子做GVB计算为例示范传统方法如何构造初始以及其弊端。

第一步，明确研究体系或感兴趣的过程涉及几根化学键以及哪些轨道，从而确定GVB计算的对数NPair。这里为了简便，假设我们要研究CH₄分子的四根CH键，即NPair=4。但是注意，在大家的实际研究中，一个具体体系或过程涉及究竟几根化学键可能未必那么明白，需要丰富的计算经验和准确的化学直觉才能够正确确定，所以初学者使用起来并不友好。

第二步，进行Hartree−Fock（HF）计算，拿到分子轨道。传统GVB计算的初始构造是基于HF轨道，且通常是限制性HF（RHF）轨道，所以一定要先进行HF计算拿到分子轨道。比如在6-31G基组水平对平衡结构的CH₄分子进行简单的RHF计算的输入文件如：

%chk=CH4_6-31g.chk
#p HF/6-31G

RHF

0 1
C    -0.86325962    0.58011049    0.00000000
H    -0.50660520   -0.42869951    0.00000000
H    -0.50658678    1.08450868    0.87365150
H    -0.50658678    1.08450868   -0.87365150
H    -1.93325962    0.58012367    0.00000000

细心的朋友可能立马会问，对远平衡位置的结构，RHF不稳定（非限制性HF（UHF）解更低）时怎么办？我们研究发现，此时可以采用UHF自然轨道（UHF natural orbital, UNO）。具体而言，对RHF不稳定的体系（判定RHF稳定性，可以上面关键词行加上stable进行稳定性检查），先进行UHF计算，再将得到两组UHF轨道转化为一组 UNO（原理后面会有叙述）。输入文件如下：

%chk=xxx.chk
#p UHF/6-31G nosymm guess=mix stable=opt

UHF

0 1
坐标略

--link1--
%chk=xxx.chk
#p UHF chkbasis guess=(read,NaturalOrbitals,save,only) geom=allcheck nosymm

注意做完UNO之后得到的是一个含两组轨道的fchk文件，需要先转化为一组轨道的fchk文件。因为UNO实际上只有一组，高斯fchk文件此处有bug，会存两组一模一样的UNO轨道，而我们只需要一组用来局域化和后续做GVB。将UHF的fchk文件转化为RHF或ROHF的fchk文件，可以自己打开fchk文件手动删除所有Beta部分（不推荐），也可以下载fch_u2r小程序（https://gitlab.com/jxzou/mokit/-/releases）一键运行（推荐做法）。注意，此处一定要加上only，否则又只能得到UHF轨道（详见Gaussain手册）。在得到RHF轨道或UNO轨道后，将轨道局域化（参见《局域分子轨道简介》和《利用常见的程序做轨道局域化》）。

%chk=CH4_6-31g.chk
#p HF chkbasis geom=allcheck guess=(read,local,save,only) nosymm

在GVB计算中，局域化是非常关键的一步，因为它的好坏直接影响着后面识别轨道和配对。局域化较好时，轨道会很好地局域在一个化学键上，易于确定并配对轨道；不好时，其轨道离域在几个化学键上甚至整个体系，无法辨认是哪个轨道。对占据轨道，轨道较少，局域化容易进行，但对虚轨道，特别是基函数多、体系复杂时，局域化可能不收敛或者即使收敛但效果不好，所以虚轨道的局域化一直以来还是一个非常困难的问题。

最后，手动挑选轨道并配对。局域化后，我们可以借助GaussView或者其它可视化软件查看具体轨道。比如经过局域化后，个别轨道（4，5，14，16，其中4和5是占据轨道，14和16是虚轨道）如下所示：

为方便观看，上图旋转过视角。可以看到，对平衡结构的CH₄分子，局域化效果较好，每个轨道都局域在一个化学键上。另外，我们也可以看到，成键轨道（4和5，占据轨道对应成键轨道）局域在化学键上，而反键轨道（14和16，虚轨道对应反键轨道）在化学键上有一个节面，这正是成键轨道和反键轨道的典型特征。有了局域化较好的分子轨道，现在我们将它们一一配对。配对的含义是将成键轨道和它相应的反键轨道组成一对，即构成轨道对（Pair）。而Gaussian程序表示轨道对是使反键轨道采用逆成键轨道序，即最高占据轨道（HOMO）与最低非占据轨道（LUMO）是一对，HOMO-1与LUMO+1是一对，依次类推。注意GAMESS程序里不是这个轨道顺序，以后另说。

而局域化后，便不再有轨道能量的概念，轨道的排序也会比较任意。此时需要肉眼观察轨道之间的成键和反键关系，并按照上述规则重新排序。如上面CH₄分子的轨道，5号轨道的反键轨道应该在6号位置，但它却是14号轨道。Gaussian程序中使用guess=alter关键词调换轨道顺序。输入文件可以如下书写： 

%chk=CH4_6-31g.chk
#p HF chkbasis geom=allcheck guess=(read,alter,only,save) nosymm

6,14
7,16
8,15
9,13

加only关键词的原因同上。这样出来的结果，轨道顺序配对正确，就算是正确地为GVB计算构造了一个初始。接着再优化分子轨道即可顺利地得到收敛的GVB结果。 

%chk=CH4_6-31g.chk
#p GVB(4)/6-31G guess=read geom=allcheck nosymm

2 2 2 2

这后面4个2是说明这个计算有4个Pair，每个Pair中有2个轨道（这些信息有些多余，但Gaussian软件要求这样）。

可以看到，整个构造初始的操作是非常繁琐的，既要让程序停下来，也要使用其它软件看轨道，还要调换轨道顺序，若局域化效果不好，就无法正确地配对轨道。即使是效果好，对几对或十几对大小体系，尚且可以手动完成，但若有几十对、上百对大小的体系，轨道配对就变得非常麻烦甚至完成不了。所以，自动地为GVB计算构造初始是非常必要的。

自动构造初始

鉴于GVB计算构造初始如此麻烦，我们课题组2019提出了一个自动构造GVB初始的方法（J. Chem. Theory Comput. 2019, 15, 141−153），该方法能使上百对的GVB计算像HF计算一样自动进行。现将这个方法简单介绍如下：

第一步，确定GVB对数。简单起见，可取价轨道和最小基下虚轨道数目较小者，因为一个体系的内核轨道通常是非活性的，而描述一个体系的基组水平至少是最小基，即：

其中是价电子占据轨道的数目，是最小基下虚轨道的数目。例如对于CH₄而言n=i=b=4，即要做GVB(4)，与四根C-H键数目相等。

第二步，确定活性的占据轨道与虚轨道。这里分两种情形。对RHF稳定的体系，活性的占据轨道可由价轨道构成，而活性的虚轨道能通过奇异值分解(SVD)得到体系在其最小基上的投影虚轨道，因为该空间是描述体系最小有效虚轨道空间。

其中、 分别是某给定基组（例如cc-pVTZ）和最小基（如STO-6G）下的虚轨道。而对RHF不稳定的体系，Pulay教授在1988年曾提出过一个比较简单地确定活性轨道的方法，将两组UHF轨道转化成一组UNO，其中分数占据的部分就被认为成活性轨道。假设自旋向上、自旋向下的UHF占据轨道分别为和 ，

对应轨道（或）是自旋向上占据轨道（或自旋向下占据轨道）的酉变换，

那么两个UNO（成键和反键）可以表示如下：

其中和分别对应着和的轨道占据数。不过，这个方法在不同几何结构下，得到的活性轨道数目可能不同。这很好理解，因为不同结构下<S**2>值是不一样的。利用我们在RHF稳定情形时的做法，先确定活性空间大小，再做占据数截断，可以避免这个问题。

第三步，分别将活性的占据轨道和虚轨道局域化，这里采用Boys或Pipek−Mezey方法都可以。

最后，利用前面推送的Kuhn−Munkres算法将局域的占据轨道和虚轨道一一配对构成GVB初始轨道。只是，在这里权重要用两个轨道之间空间邻近指标来表示，比如

其中表示标准的基于分子轨道的偶极积分。

整个构造GVB初始的过程可以下面这个流程表示，

Python实现的代码已上传至GitHub（https://github.com/cingchun/GVB），欢迎大家下载使用。若使用请记得引用文章J. Chem. Theory Comput. 2019, 15, 141−153。利用这种自动构造的初始，再结合Gaussian或GAMESS软件，曾经麻烦的GVB计算可以像HF计算一样简单方便。比如分子，考虑全部价电子共有120对，其中π电子对有30对。在Cartesian型cc-pVDZ基组水平，体系的基函数共有900个。像这样一个复杂的体系，GVB计算能够轻松顺利地进行。部分GVB轨道如下所示：

此外，值得说明的是，这个自动构造初始的方法对其它多参考计算也有帮助。比如利用GVB波函数确定最小有效活性空间（见我们的原始文章）或者利用GVB波函数挑选活性轨道（马海波教授和黎书华教授课题组最近合作的一个工作（DOI：10.1021/acs.jpca.0c05216）就测试了这种做法的优越性）。

参考资料

[1] Shen, J.; Kou, Z.; Xu, E.; Li, S. The coupled cluster singles, doubles, and a hybrid treatment of connected triples based on the split virtual orbitals. J. Chem. Phys. 2012, 136, 044101.

[2] Pulay, P.; Hamilton, T. P. UHF natural orbitals for defining and starting MC-SCF calculations. J. Chem. Phys. 1988, 88, 4926−4933.

致谢

本文部分图表和算例由jxzou完成，在此表示感谢！