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

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

文/Sobereva@北京科音  2018-Nov-11


摘要:本文非常简要地介绍一下对非限制性开壳层波函数的alpha和beta轨道之间做双正交化的方法,并且以三重态乙醇为例介绍如何在Multiwfn中实现,使得其alpha和beta轨道最大程度匹配,从而便于讨论轨道。本文说的Multiwfn及其手册是2018-Nov-11更新的3.6(dev)版,不要用更老的版本。


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计算给出的是完全相同的,在这点上类似于轨道定域化。双正交化后的轨道不再是HF或KS算符的本征函数,因此双正交化轨道也没有确切的轨道能量,至多以HF/KS算符期望值的方式求轨道能量。


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
从屏幕上提示的信息可见双正交化分为三步进行,比如第一步输出的信息为
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文件中。此时程序还问你是否现在就载入这个文件,如果选y将之载入的话,内存里的轨道就对应于双正交化后的轨道了,之后你可以进入主功能0去查看这些轨道的图形,也可以用Multiwfn中的各种功能去分析轨道特征,比如计算轨道成分、轨道动能等。

这里把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对应。像这样重要的信息,如果我们不做双正交化的话,是没法这么容易了解到的。

注意Multiwfn并不计算双正交化轨道的能量,因此biortho.fch里的轨道能量还是之前分子轨道的,这些能量此时已经没有意义了。另外,由于没有实际双正交化后的轨道的能量信息,此文件里的这些轨道的序号顺序并不对应于其轨道能量顺序。


2.2 丁烷双自由基

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

评论已关闭