用于非限制性开壳层波函数的双正交化方法的原理与应用

sobereva

用于非限制性开壳层波函数的双正交化方法的原理与应用

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文件观看以上提及的轨道的图形。

granvia

Just a minor question: The biorthogonalized orbitals are not exactly doubly occupied and thus cannot constitute Slater determinants. So, how come the total energy is invariant?

sobereva

granvia 发表于 2018-11-11 14:58
Just a minor question: The biorthogonalized orbitals are not exactly doubly occupied and thus cannot ...

双正交化的轨道是整数占据的。类似于MO变换成定域化轨道，总能量是不变的。因为只不过是占据轨道之间的酉变换。

granvia

sobereva 发表于 2018-11-12 01:44
双正交化的轨道是整数占据的。类似于MO变换成定域化轨道，总能量是不变的。因为只不过是占据轨道之间的酉 ...

哦，明白了。是不是也意味着，双正交化轨道之间一般并不是严格的正交，因为变换的空间不是整个希尔伯特空间？

sobereva

granvia 发表于 2018-11-12 08:04
哦，明白了。是不是也意味着，双正交化轨道之间一般并不是严格的正交，因为变换的空间不是整个希尔伯特空 ...

双正交化之后alpha和beta也不严格正交归一

granvia

sobereva 发表于 2018-11-12 09:05
双正交化之后alpha和beta也不严格正交归一

All right.

xujian

请问老师，对于双自由基，有一个alpha占据轨道没法与beta占据轨道相配对，重叠积分的大小与双自由特征有什么关系，即重叠积分越小，说明双自由特征越明显吗？

sobereva

xujian 发表于 2018-11-14 09:45
请问老师，对于双自由基，有一个alpha占据轨道没法与beta占据轨道相配对，重叠积分的大小与双自由特征有什 ...

是

xujian

sobereva 发表于 2018-11-14 17:18
是

谢谢老师，如果在文章中用这一方法需要引用文献说明吗

sobereva

xujian 发表于 2018-11-14 20:22
谢谢老师，如果在文章中用这一方法需要引用文献说明吗

引用Multiwfn 3.6版手册3.100.12节

adia

请问老师，该方法对于过渡金属配合物可能出现不稳定的情况吗，我计算了一个五重态高自旋的Fe(III)配合物，Fe上自旋密度显示的是4.2左右，加上配位原子上的，接近4.8，应该是很标准的三价铁，但是我用该方法分析单电子轨道的时候，竟然发现有两个dz2轨道，而不是五个d轨道一样一个，很费解！

sobereva

adia 发表于 2018-12-16 19:06
请问老师，该方法对于过渡金属配合物可能出现不稳定的情况吗，我计算了一个五重态高自旋的Fe(III)配合物，F ...

本文的方法适合任何体系

那不叫自旋密度，叫自旋布居，看
量子化学中的一些常见不良写法和用词
http://sobereva.com/298

我不清楚所谓的有两个dz2指什么。双正交化只是个数学处理，不意味着双正交化得到的轨道中某些能与各个d原子轨道很好对应上。

sobereva

在今天更新的Multiwfn 3.6(dev)中，输出的biortho.fch文件里的能量部分改为了记录双正交化过程产生的singlar value，便于用户在主功能0的看轨道列表信息的界面里考察。相应地，对帖子进行了更新。

sobereva

对本文进行了重大更新，目前官网上最新的Multiwfn已经可以计算双正交化轨道的能量了。

sobereva

对程序和本文进行了更新，在最新版程序中，进入双正交化功能后，用户可以选择是只对占据轨道做双正交化还是对所有轨道都做双正交化，选前者可以略过对alpha-beta空轨道间的双正交化步骤，使得对于大体系可以节约一个数量级以上的耗时，但此时的结果里编号在N(alpha)+1以上的alpha和beta轨道将没有意义，不过这些一般也不是我们感兴趣的。