用于非限制性开壳层波函数的双正交化方法的原理与应用
用于非限制性开壳层波函数的双正交化方法的原理与应用
Principle and application of biorthogonalization method for unrestricted open-shell wavefunctions
文/Sobereva@北京科音
First release: 2018-Nov-11 Last update: 2021-Sep-10
摘要：本文非常简要地介绍一下对非限制性开壳层波函数的alpha和beta轨道之间做双正交化的方法，并且以三重态乙醇为例介绍如何在Multiwfn中实现，使得其alpha和beta轨道最大程度匹配，从而便于讨论轨道。本文说的Multiwfn及其手册是官网上最新版本的情况。
1 相关知识
众所周知，非限制性开壳层计算（如UHF、UKS，以下简称为U）的时候alpha和beta是分别求解的，因此会产生alpha和beta这两套自旋轨道。对于自旋多重度>1的体系，以及自旋多重度为1的对称破缺态，由于自旋极化，会导致alpha和beta轨道不匹配，此时虽然alpha轨道自己是正交归一的，beta轨道自己也是正交归一的，但是alpha和beta之间不满足正交归一关系。更具体来说，第i号alpha轨道和第i号beta轨道之间重叠积分不为1，这俩轨道既可能轨道形状稍有偏差，也可能完全不同。因此，对于非限制性开壳层波函数，讨论轨道的时候必须分别去考察alpha和beta轨道，这是比较麻烦的事情。而且这种情况，体系的自旋密度是由所有占据轨道所决定的，因此没法只拿某几条轨道来讨论自旋密度。（不知道什么是自旋密度的话看《谈谈自旋密度、自旋布居以及在Multiwfn中的绘制和计算》http://sobereva.com/353）。
虽然限制性开壳层(RO)计算只会产生一套轨道，没有上述U计算的麻烦，但是相对于U，RO计算耗时高、难收敛，自旋密度、体系总能量和轨道能量也都没有U那么有意义，因此用RO回避上述问题也不是好的办法。
实际上，在U计算之后，可以做双正交化(biorthogonalization)变换来解决上述问题。双正交化的具体细节和实现在Multiwfn手册3.100.12节有详述，本文就不细说了。简单来说，双正交化是对alpha和beta轨道进行酉变换，使得alpha和beta之间完美满足或基本满足正交归一关系，同时又不破坏每种自旋轨道与其自己原有的正交归一关系。双正交化具体实现方式并不唯一，Multiwfn中的双正交化方法分成三步，经过这三步变换后，每个alpha和与之序号相同的beta轨道一般可以达到近乎完美匹配，因此之后就只需要讨论一套轨道即可（类似于RO轨道的情况）。
在双正交化过程中每种自旋的占据轨道和占据轨道会发生混合、空轨道和空轨道会发生混合（具体来说是酉变换），而不同自旋轨道间不会发生混合。这种变换不会导致体系任何可观测性质发生改变，因此如果使用双正交化后的轨道去计算体系总能量、总电子密度、自旋密度等等，结果和U计算给出的是完全相同的，在这点上类似于轨道定域化。双正交化后的轨道不再是Fock算符的本征函数，因此双正交化轨道也没有确切的轨道能量，而只能以Fock算符期望值的方式求轨道能量。
2 实例
2.1 三重态乙醇
下面我们对三重态乙醇做双正交化，来看看双正交化的实际效果和价值。首先我们先看一下UB3LYP/6-31G**下直接给出的这个体系三重态的几个分子轨道等值面图，如下所示。其中14号alpha轨道是alpha轨道中的HOMO，12号beta轨道是beta轨道中的HOMO。
由图可见，只有12号alpha与12号beta轨道匹配很好，形状和相位都肉眼看不出显著差别，而其它alpha和beta轨道之间完全匹配不上。alpha比beta占据轨道多两个（13和14），但显然当前情况我们不可能认为体系的自旋密度等于13和14号alpha轨道对应的密度之和，因为由于其它alpha和beta占据轨道形状的不匹配，它们也对自旋密度有直接的影响。
在Multiwfn中做双正交化需要的输入文件必须包含基函数信息，具体看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）、《Multiwfn入门tips》（http://sobereva.com/167），对于Gaussian用户就用fch文件即可。UB3LYP/6-31G**下三重态乙醇的fch文件已经在Multiwfn自带的文件包里提供了，这里我们对它做双正交化。
启动Multiwfn，然后依次输入
examples\ethanol_triplet.fch
100
12 //产生双正交化轨道
2 //双正交化过程中考虑所有轨道（对于大体系，为节约时间，可以此处选1略过空轨道之间的双正交化，此时对于此例，在最终结果中编号在14以上的alpha和beta轨道将无意义）
0 //不计算双正交化轨道能量
从屏幕上提示的信息可见双正交化分为三步进行，比如第一步输出的信息为
Doing biorthogonalization for alpha 1 to 14, Beta 1 to 12
Singular values of orbital overlap matrix:
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9999
0.9999 0.9998 0.9995 0.9992
这说明这一步是对1~14号alpha轨道和1~12号beta轨道做双正交化。下面的数值是经过双正交化后的alpha和beta轨道间的重叠积分。当前显示了12个等于或非常接近于1.0的数值，这代表这步双正交化后1~12号alpha和1~12号beta轨道已经非常完美地一一对应了。其余的轨道的双正交化会在接下来的步骤中进行。
之后Multiwfn把双正交化后的轨道波函数自动导出为当前目录下的biortho.fch文件，也把双正交化后的轨道信息导出为当前目录下的biortho.txt文件，此文件内容如下所示。
S = Singular value, E = Energy (in eV), O= Occupancy, A=Alpha, B=Beta
Orb: 1 S= 1.0000 O(A)= 1.0 O(B)= 1.0
Orb: 2 S= 1.0000 O(A)= 1.0 O(B)= 1.0
...[略]
Orb: 11 S= 0.9995 O(A)= 1.0 O(B)= 1.0
Orb: 12 S= 0.9992 O(A)= 1.0 O(B)= 1.0
-----------------------------------------------
Orb: 13 S= 1.0000 O(A)= 1.0 O(B)= 0.0
Orb: 14 S= 1.0000 O(A)= 1.0 O(B)= 0.0
-----------------------------------------------
Orb: 15 S= 1.0000 O(A)= 0.0 O(B)= 0.0
Orb: 16 S= 1.0000 O(A)= 0.0 O(B)= 0.0
Orb: 17 S= 1.0000 O(A)= 0.0 O(B)= 0.0
...[略]
通过这个文件可以非常方便地考察双正交化后的alpha轨道与序号与之相同的beta轨道的匹配程度。由于上面列出来的轨道的S（Singular value，即相当于重叠积分）都很接近1，因此这些轨道的alpha-beta匹配程度都很好，只有alpha 12和beta 12的匹配程度稍差一点，但也很接近于1.0，因此这俩轨道的图形从肉眼上应该很难分辨。上面信息里O代表Occupancy，就是占据数的意思，为了观看方便，这个文件里把O(A)=O(B)=1.0、O(A)=1.0&O(B)=0.0和O(A)=O(B)=0.0三个部分用横线做了分隔。
此时程序还问你是否现在就载入biortho.fch文件，如果选y将之载入的话，内存里的轨道就对应于双正交化后的轨道了，之后你可以进入主功能0去查看轨道图形，也可以用主功能0的图形界面左上角的Orbital info.选项去查看这些轨道的基本信息，还可以用Multiwfn中的各种功能去分析轨道特征，比如计算轨道成分、轨道动能等。由于此例没计算双正交化轨道的能量，因此在导出的biortho.fch里“轨道能量”部分记录的是上述输出的singular value值。由于当前没有双正交化轨道的能量信息，因此当前biortho.fch里的这些轨道的序号顺序并不对应于其轨道能量顺序。
这里把11~14号双正交化后的alpha和beta轨道列出：
可见，双正交化后alpha和与之序号相同的beta在轨道形状上已经肉眼基本看不出任何区别了。此时，13和14号alpha轨道可以被冠以描述RO轨道时用的SOMO（单占据轨道）之名，这两个轨道对应的密度的加和正好就是当前体系UB3LYP/6-31G**级别计算的三重态的自旋密度。对比本文上面的两张图，也可以看到在双正交化前后，轨道特征往往差别不小。双正交化后的第13号alpha轨道基本正好对应原先的第14号alpha轨道，而双正交化后的第14号alpha轨道则不直接对应于之前任何一个alpha轨道，显然必定是两个或多个原先的alpha占据轨道发生了显著混合才产生出来的。
值得一提的是，乙醇的单重态基态在B3LYP/6-31G**下计算产生的HOMO和LUMO轨道图形如下
这个HOMO和LUMO分别非常像双正交化后的14号和13号alpha轨道。因此我们可以立刻明白，乙醇从单重态基态激发到最低三重态的过程，如果以轨道跃迁方式描述的话，可以近似说成是HOMO→LUMO的跃迁，因为这样跃迁后HOMO和LUMO就各有一个单电子了，这俩轨道又几乎恰好和双正交化后的两个SOMO对应。像这样重要的信息，如果我们不做双正交化的话，是没法这么容易了解到的。
2.2 三重态乙醇（计算双正交化轨道能量）
这个例子还是考察三重态乙醇，和上例不同的是这次我们也让Multiwfn产生双正交化轨道的能量，这样不仅可以讨论更多信息，而且Multiwfn会根据能量对轨道排序，讨论起来更方便。Multiwfn计算双正交化轨道的能量需要利用Fock矩阵，Multiwfn可以通过两种方式获得这个矩阵：
(1)直接基于分子轨道能量和系数矩阵反解出来Fock矩阵。通常建议用这个方式，最为省事。但是如果使用了大量弥散函数，量子化学程序可能会自动去除一些线性相关基函数，这个时候就没法用这个做法了，只能用下面的方法
(2)从特定文件里读取。可以从记录Fock矩阵的文本文件里读取，也可以从.47文件里读取（Gaussian可以直接产生），也可以从ORCA输出文件里读取，等等，详情见Multiwfn手册附录7。
启动Multiwfn，输入
examples\ethanol_triplet.fch
100
12 //产生双正交化轨道
2 //双正交化过程中考虑所有轨道
1 //计算双正交化轨道的能量，Fock矩阵通过分子轨道的能量和展开系数直接产生
y //根据能量对双正交化轨道排序
此时产生了biortho.fch和biortho.txt，后者内容如下
...[略]
Orb: 10 S= 1.0000 E(A)= -13.621 O(A)= 1.0 E(B)= -13.450 O(B)= 1.0
Orb: 11 S= 0.9992 E(A)= -13.057 O(A)= 1.0 E(B)= -11.821 O(B)= 1.0
Orb: 12 S= 1.0000 E(A)= -11.918 O(A)= 1.0 E(B)= -11.866 O(B)= 1.0
-------------------------------------------------------------------------------
Orb: 13 S= 1.0000 E(A)= -11.439 O(A)= 1.0 E(B)= -5.806 O(B)= 0.0
Orb: 14 S= 1.0000 E(A)= -0.293 O(A)= 1.0 E(B)= 2.817 O(B)= 0.0
-------------------------------------------------------------------------------
Orb: 15 S= 0.9992 E(A)= 13.876 O(A)= 0.0 E(B)= 15.431 O(B)= 0.0
Orb: 16 S= 1.0000 E(A)= 15.057 O(A)= 0.0 E(B)= 15.305 O(B)= 0.0
Orb: 17 S= 1.0000 E(A)= 18.508 O(A)= 0.0 E(B)= 19.028 O(B)= 0.0
...[略]
由于这回我们要求Multiwfn计算双正交化轨道的能量，因此这次的biortho.txt里显示了算出来的轨道能量。由于我们也要求按照能量进行排序，因此可以看到轨道序号顺序确实和轨道能量相同了。注意排序的时候是将序号相同的alpha和beta的轨道能量平均值作为依据来排序的。另外，排序是对每一个批次产生的双正交化轨道分别进行排序的（上面信息中的两条横线把所有轨道分割成了三个批次），因此比如第2批次的轨道的能量不一定全都高于第1批次的轨道。
注：常有Multiwfn用户问为什么他得到的序号相同的alpha和beta双正交化轨道的形状看起来非常相似，但能量却差了非常多。这是因为开壳层体系的alpha和beta电子的分布通常不对称，而且通常连alpha和beta电子数都往往不同，也因此alpha和beta对应的Fock算符不同。而由于轨道能量是相应自旋的Fock算符求期望得到的，因此哪怕形状精确相同、序号相同的alpha和beta双正交化轨道的能量也注定会不同。如果你就是死活非得要求相同序号的alpha和beta能量相同（就跟限制性开壳层RO计算的情况似的），那你姑且把相同序号的alpha和beta轨道能量取平均得了（虽说没有明确的物理意义），这时的轨道能量相当于轨道波函数相对于alpha和beta电子的平均的Fock算符的期望值。
当前biortho.fch里的“轨道能量”信息是双正交化轨道的真实能量了。如果我们选择y载入biortho.fch，用主功能0的Orbital info.选项查看轨道基本信息，会看到下面这些，确实轨道能量和biortho.txt里显示的相同。
1 E(au/eV): -19.26817 -524.3135 Occ: 1.000000 Typ: A
...[略]
11 E(au/eV): -0.47985 -13.0573 Occ: 1.000000 Typ: A
12 E(au/eV): -0.43799 -11.9184 Occ: 1.000000 Typ: A
13 E(au/eV): -0.42037 -11.4390 Occ: 1.000000 Typ: A
14 E(au/eV): -0.01078 -0.2933 Occ: 1.000000 Typ: A
76 ( 1) E(au/eV): -19.23884 -523.5155 Occ: 1.000000 Typ: B
...[略]
86 ( 11) E(au/eV): -0.43441 -11.8210 Occ: 1.000000 Typ: B
87 ( 12) E(au/eV): -0.43606 -11.8659 Occ: 1.000000 Typ: B
对于beta轨道，上面信息中括号里的序号是从第一条beta轨道开始计的，括号外的是从第一条alpha轨道开始计的。
看一下现在的第13号和第14号alpha轨道。如之前所述，这俩占据轨道几乎完全描述了当前体系的单电子分布，因为其它的占据的alpha轨道都有与之匹配程度极高（重叠积分接近1）的beta占据轨道故而对自旋密度贡献基本为0。这俩轨道的图形和能量如下所示
可见由于Multiwfn对轨道做了排序，当前轨道能量和轨道序号是一致的（而上例中没排序，它俩序号是反着的）。另外，如果你将这俩双正交化轨道和ethanol_triplet.fch里记录的MO13（-10.083 eV）、MO14（-0.276 eV）对比的话，会发现轨道形状相差不太大，能量也相差不算特别大，这也体现Multiwfn给出的双正交化轨道的能量是合理的。
2.3 丁烷双自由基
本节用到的fch文件可在此下载：http://sobereva.com/attach/448/C4H8-singlet.rar。如果对双自由基计算不了解，参看《CASSCF计算双自由基以及双自由基特征的计算》（http://sobereva.com/264）和《谈谈片段组合波函数与自旋极化单重态》（http://sobereva.com/82）。
通过非限制性开壳层以对称破缺方式算的自旋极化单重态体系，比如双自由基、反铁磁性耦合体系，双正交化也可以做，但是此时不可能使得每个alpha轨道都能与一个beta轨道完美相对应。比如对丁烷双自由基体系，双正交化时候的第一部分的输出是这样的
Doing biorthogonalization for alpha 1 to 16, Beta 1 to 16
Singular values of orbital overlap matrix:
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
0.9998 0.9998 0.9998 0.9998 0.9992 0.9990 0.9989 0.3095
可见，此时双正交化变换的效果是令alpha与beta凡是能匹配的都尽量匹配，这使得前15个轨道都有一一对应关系，重叠积分基本是1.0，但有一个alpha占据轨道就是没法与beta占据轨道相配对，重叠积分才0.3095。这样的输出信息体现在双正交化之后，我们可以就只依靠16号alpha和beta轨道来讨论体系的单电子分布以及双自由基本质特征。
顺带一提，如果想查看原本alpha和beta分子轨道之间的重叠积分，可以用Multiwfn主功能100的子功能5，具体说明详见Multiwfn手册3.100.5节。如果用这个功能里面的选项2输出对上述双自由基波函数的alpha MO和序号相同的beta MO之间的重叠积分，结果如下（以下只列出占据轨道部分）
Overlap between the 1th alpha and beta orbitals: 0.036957
Overlap between the 2th alpha and beta orbitals: -0.036727
Overlap between the 3th alpha and beta orbitals: -0.001295
Overlap between the 4th alpha and beta orbitals: 0.001066
Overlap between the 5th alpha and beta orbitals: 0.982456
Overlap between the 6th alpha and beta orbitals: -0.958735
Overlap between the 7th alpha and beta orbitals: 0.858202
Overlap between the 8th alpha and beta orbitals: -0.882210
Overlap between the 9th alpha and beta orbitals: 0.994529
Overlap between the 10th alpha and beta orbitals: -0.990360
Overlap between the 11th alpha and beta orbitals: 0.991034
Overlap between the 12th alpha and beta orbitals: -0.984323
Overlap between the 13th alpha and beta orbitals: 0.991434
Overlap between the 14th alpha and beta orbitals: -0.991413
Overlap between the 15th alpha and beta orbitals: -0.995594
Overlap between the 16th alpha and beta orbitals: 0.294899
可见，如果不做双正交化，内核轨道部分（MO 1~4）alpha-beta匹配极差，看图可知主要是因为碰巧定域在了不同原子上，这无所谓。而价层部分，7、8号轨道alpha-beta之间匹配得也不很完美，重叠积分还不到0.9，因此比如讨论自旋密度、双自由基特征时如果简单忽略它们的影响原理上并不好，而在双正交化后就仅仅需要关注第16号轨道了。请大家自行用Multiwfn基于本文提供的fch文件观看以上提及的轨道的图形。