【科普】臭氧的电子结构和开壳层单重态

本文为了梳理整体思路，部分内容采用插叙和倒叙，如果暂时不理解可以尝试先往后看。
看似简单的臭氧分子其实电子结构非常复杂，这点是有些出乎意料的。对于这类电子结构复杂的体系，通常我们会有“多参考特征强”这种描述方法。那么，这究竟指的什么呢？臭氧的电子结构究竟是什么样子呢？我们一点一点来阐述。
本文计算的臭氧分子均为气态单分子，在M06-2X/def2-SVP水平下用“正确”的方法优化得到的结构。(什么叫“正确”的方法后面会说。)计算使用Gaussian软件，注意M06-2X这种Minnesota系列的泛函对DFT积分格点要求比较高，因此所有用M06-2X的计算都加了Integral=UltrafineGrid关键词。至于为什么用这个泛函，是因为这个泛函对主族描述非常好，且HF成分很高，容易找到对称性破缺态(这个概念后面会阐述)。
优化的结果为键长1.255Å，键角116.0°，属于C2v点群即两根O-O键等长，和水分子几何结构对称性相同。
我们先用RM06-2X/def2-SVP(Gaussian里写作RM062X/def2SVP)计算它的波函数，得到电子能量（这个能量绝对数值没有意义，会和后面对比做差）为-225.123589439 Hartree。这里的R代表限制型方法，即让α自旋的电子和β自旋的电子完全配对，从而只用空间部分波函数描述一对电子，让Slater行列式的维度减小了一半，显著节省了计算量。但是这个真的就是它合理的波函数么？
我们用Stable关键词对波函数稳定性进行测试，结果提示
The wavefunction has an internal instability.
。这说明我们得到的波函数不是Kohn-Sham方程的极小值解，而是得到了一个能量较高的解。换言之，我们计算的态其实某种意义上是分子在该结构下的“激发态”。那么我们怎么才能得到真正的基态呢？
我们上面用了Restricted方法，这在Gaussian等程序中对自旋多重度为1（即单线态）的体系是默认的。绝大部分分子，尤其是有机分子，只要是总共偶数个电子，而且基态还是单线态(α自旋和β自旋电子数相同，比如氧气基态α电子比β电子多两个，是三线态，不在本文讨论范围内)。但是对于个别体系（事实上，有机体系也有一些这种特征的体系，比如乙烷的碳碳键拉长到5Å，乙烯的双键两边的二面角转到90°，丁烷-1,4-二自由基等体系，但是通常都是严重偏离能量最小几何构型的结构。常见有机甚至主族体系中，能量最低结构居然还有这类自旋极化特征的单线态分子，尤其还是常见的分子中，几乎只有臭氧了。），这个假设可能是不够合理的。
我们尝试打破这种限制，改用Unrestricted方法，即UM06-2X/def2-SVP(Gaussian里写作UM062X/def2SVP)进行计算。但是即使这样，由于α电子和β电子的初猜和之前是相同的，又由于自旋部分正交，不会相互混合，因此自洽场迭代以后结果还是和刚才一样的。所以我们再做一些调整，人为得到两种自旋电子不配对的初猜来尝试。一个方法就是使用Guess=Mix关键词，这对非限制型方法在初猜的时候会把α电子和β电子各自的HOMO和LUMO以不同方式混合（波相相加，或相减，然后归一化）这样得到的初猜波函数，α电子的HOMO电子和β电子的HOMO电子的波函数就有了差别。如果最终能量最低的态真的是两种自旋HOMO和LUMO分别不相同的态(这种特征通常也只在前线轨道附近最明显)，那么经过自洽场的迭代，就有希望找到这个解。像这样的α电子和β电子空间部分不完全相同的态，我们统一叫做自旋极化态。由于α电子和β电子各自的轨道这时候已经不符合分子本身C2v的对称性，而是有类似一个单电子主要在左边，另一个主要在右边这种情况，对称性被打破，所以我们也称为对称性破缺态。
另一种构建初猜的方法是拆成几个片段，认为指定初猜α电子和β电子的单电子一个单电子主要在左边，另一个主要在右边。这叫做片段组合波函数获得对称性破缺初猜的方法。
但是对称性破缺的态未必只有一个，能量也不相同。真正想找到能量最低的态是很困难的。为此很多量化软件都有测试波函数稳定性和想尽各种办法优化到最优的波函数。这方面Gaussian是做的相当好的软件之一。用Stable=Opt关键词，并加上NoSymmetry关键词不使用对称性即可，这样不管什么初猜，原理上最后优化后的波函数结果应该没有差别。对于臭氧，当前水平下，用片段组合波函数做初猜可以得到能量最低的态，但是直接混合HOMO和LUMO得到的态有问题。
这个态大多数电子都是基本成对的(如果要配对更好，可以做双正交化)。但是α电子和β电子各自的HOMO在空间部分是不同的，自旋密度等值面图如下图所示。这才是单参考行列式能描述的臭氧最真实的结构。前面说的“正确的方法优化”，也是在这个方法构建出某个几何构型的对称性破缺波函数后，用这个波函数做初猜进行优化的。
这样得到的对称性破缺的单组态波函数，α和β电子的HOMO轨道，也就是两者自旋相反但是空间上不成对的那两个单电子的轨道如下图所示。红色为波函数正值，蓝色为负值。
α-HOMO

β-HOMO

明显一个单电子主要分布在左边的氧原子上，另一个主要分布在右边的氧原子上。考虑到其它电子也可能未必完全配对，我们可以考察一下自旋密度。如下图。绿色和蓝色分别代表正值和负值。

我们还可以用各种原子空间划分，把自旋密度投射到某一个原子上， 这叫做自旋布居。对于这个体系虽然没有意义，但是毕竟是一种分析手段，我们也做一下，通常用模糊空间比如Hirshfeld或Becke划分效果比较好。这里我们用之前别人计算好的的原子波函数来计算的权重的Hirshfeld空间划分的自旋布居做原子填色图，蓝色为负值，红色为正值，如图。

那么能量上到底相差多少呢？做个粗略的讨论，考虑这个结构下限制型即电子全都配对的波函数和我们计算出来能量最低的波函数(-225.131722665 Hartree)做个差，为0.008133 Hartree，相当于21.4 kJ/mol，这就是本体系该结构下计算的对称性破缺使得能量降低的值。
未完待续。

这类体系，被称为自旋极化波函数/开壳层单线态/双自由基(搞量化的一般称呼的双自由基是单线态双自由基，三线态双自由基直接就叫三线态了)。自旋极化单线态的含义是体系是单线态的，但是自旋密度不处处为0，
真实的波函数是S2S2算符的本征函数。原理上只要是单线态，实验对应的S^2（自旋角动量大小的平方）应该就是s(s+1)ℏ^2，也就是0。对于这种UDFT的对称性破缺的波函数，程序给出的S2S2是0.4936。这种情况叫自旋污染，可以理解为混入了其他自旋多重度的态的波函数。因此我们用单参考行列式(Hartree-Fock,普通DFT)得到的波函数，是有问题的，尽管已经描述出了所谓的双自由基特征，能量上也是合理的，可以用于和别的分子的能量进行比较。那么，如何解决这个问题呢？
我们需要引入多参考方法，也就是用CASSCF（完全活性空间自洽场）这样的方法。这类方法以多个Hartree-Fock型行列式作为参考，在SCF（自洽场）的时候用多个行列式的线性组合来做SCF。这种方法适合描述静态相关强的体系，这类体系的多参考特征很强。它和post Hartree-Fock的区别在于后者做SCF依然只用一个行列式，而在之后再取不同的组态函数展开，做进一步线性组合，得到更好的电子密度，这样的方法可以描述动态相关效应，但是对静态相关依然描述不好。当然我们还可以在CASSCF这样的方法上进一步增加对动态相关的描述，也引入单参考的post Hartree-Fock方法，在CASSCF波函数的基础上做进一步各组态的行列式系数的变分，比如引入二阶微扰的方法之一的CASPT2。
这么说有点儿抽象，我们看一下前线轨道就明白了。
前线轨道(指HOMO和LUMO)接近简并的时候，就会出现多参考特征强的体系。
限制型波函数得到的态，这个态首先最严重的问题，就是它不是能量最低的态。

这是限制型计算的HOMO和LUMO，注意分属不同的不可约表示，计算时候不会混合。


非限制型波函数使用对称性破缺方法计算得到的是这个态。它正确描述了能量最低态的一些特征，引入了高自旋的态的特征，是一种近似描述，但是代价就是自选污染。理想的双自由基如果用非限制型单参考波函数进行对称性破缺计算，双自由基特征应该是1，S2S2也应该是1。

多参考方法描述的波函数是用这两个电子结构线性组合来描述的。

这样自旋问题也得到了解决，而且正确描述了前线轨道的简并问题。
CCSD（单激发、双激发耦合簇）下做T1诊断是最常用的衡量多参考特征强弱的方法之一，一般认为大于0.02时候多参考特征比较强。此例结果为0.026，多参考特征已经比较强了。
未完待续

做CASSCF要提供活性空间的轨道作为初猜，其中一种方式就是提供所谓的UNO轨道(非限制型自然轨道)，一种做法是上面做非限制型对称性破缺波函数得到的密度矩阵整体进行对角化，得到的本征向量就是各个UNO轨道，本征值是各个UNO轨道的占据数。原来的体系是24个电子，如果所有电子都成对就是12个占据(空间)轨道。我们看UNO的电子占据数，发现1-11号轨道都很接近2，体现了成对的电子们占的轨道。14号和更高的占据数都是0，体现了空轨道。但是12号和13号(原先的HOMO和LUMO)比较特殊，一个是1.71，另一个是0.29，都偏离预期的整数很多。这体现了这两个轨道偏离单参考很多，必须纳入活性空间去同时变分活性空间内的轨道系数和组态系数。
这是这两个UNO轨道。前面的是占据数较大的。CAS轨道和这个形状非常类似，就不画了。


按照这个占据数，可以按照一些文献的定义，计算双自由基的特征。我们只考虑前线轨道这一对的贡献，按照y=1−(2∗T)/(1+T2)y=1−(2∗T)/(1+T2)定义的双自由基特征y，其中T可以看成是双自由基特征导致的电子转移，可以看成是两个前线轨道占据数之差的一半，算得双自由基特征为0.054，越接近1特征越强，0是完全没有。可见该体系双自由基特征还不是很强。
也有按照CASSCF中，两个电子都激发上去的双激发的组态系数的平方的两倍定义双自由基特征的。这里这个组态系数是-0.4505630（基态的组态系数是0.8927446），这样计算出来的双自由基特征是0.41。
也许未完待续

走神！