浅谈VSCF求解分子振动问题

注：本文仅是对VSCF的简略介绍，在北京科音高级量子化学培训班（http://www.keinsci.com/workshop/KAQC_content.html）中的“分子振动问题的高级计算方法”一节专门详细深入介绍VSCF、VMP2、VCI等方法的原理并给出GAMESS-US程序做这些计算的实际例子，推荐希望专门具体学习这些内容的读者参加

文/Sobereva @北京科音  2013-Oct-2

在BO近似下，电子的运动与核运动完全分离开，而不必直接求解电子+核波函数。原子核可以当做是在绝热势能面上运动。在这样的势场下求解原子核的薛定谔方程，就能得到不同振动能级对应的核波函数（对于此文的范畴，具体称为振动波函数）和振动能量。振动基态与第一振动激发态的能量差对应的是基频（fundamental）；振动基态到更高振动激发态的吸收频率就是泛音带（overtune band），或称倍频。不同振动激发态之间跃迁的吸收频率则是热带（hot band）。

求解这种振动问题用的最简单的是normal-mode-harmonic-oscillator(NMHO)模型。它的过程是，首先对质权Hessian矩阵对角化，得到正则坐标{Q}和相应的力常数{k}，也就有了谐振势{V}。不同这样的谐振势之间显然是没有耦合的，假设是三个正则坐标，那么总的势此时可以精确写为三个谐振势的加和V(Q1,Q2,Q3)=V(Q1)+V(Q2)+V(Q3)。基于这样的势，求解N坐标的振动薛定谔方程就可以转化为单坐标的振动薛定谔方程[T_i+V(Qi)]φi=e_i*φi。这里φi就是第i个正则振动模式对应的振动波函数，e_i是振动能。解这个方程可以得到多个解，对应于Qi坐标上的振动基态和各个振动激发态波函数和能量。由于谐振势{V}之间无耦合，体系振动波函数就可以精确写为ψ(Q1,Q2,Q3...)=φ1*φ2*φ3。若{φ}都取振动基态对应的解，那么体系振动波函数ψ也就是基态的；若φ有的取的是振动激发态的，那么体系的振动波函数ψ也是对应的振动激态的。值得一提的是，由于正则振动是可区分的，所以连乘构成ψ的时候不必做反对称化构成Slater行列式。

对于实际情况，即便是在正则振动坐标下，势函数之间也是相互耦合的，不可能写成比如V(Q1,Q2)=V(Q1)+V(Q2)，只能写成V(Q1,Q2)=V(Q1)+V(Q2)+V'(Q1,Q2)，这里V'来表现两个坐标的耦合。就拿HCN来说，四个正则振动坐标中两个是键角变化，两个是键长变化。假设Q1对应键角变化，Q2对应键长变化，那么在Q1为不同值的时候，去扫描Q2的势能曲线，得到的曲线肯定是不一样的（显然不同键角下平衡键长值是不同的）。也就是说Q1与Q2是相互耦合的。

因此对于实际情况，由于V是无法直接分离成N个正则坐标上的势的加和的，只能求解[∑T_i+V(Q1,Q2,Q3...)]ψ(Q1,Q2,Q3...)=E*ψ(Q1,Q2,Q3...)这个N坐标振动薛定谔方程，故实际的振动波函数必定包含N个坐标，也就是类似post-HF多电子波函数，而没法精确写为单坐标振动波函数的乘积形式。

Vibrational self-consistent field (VSCF)方法用于求解出最好的单坐标振动波函数。虽然V是无法分离的，但是VSCF把核坐标之间的耦合效应转化成正则坐标下的等效势V`来描述，由此就可以求解出单坐标的振动波函数。VSCF还是先对Hessian矩阵对角化获得正则振动坐标{Q}，然后比如对于第i个正则振动模式，求解[T_i+V(Qi)+V`(Qi)]φi=e_i*φi。其中V是在Qi坐标上扫描出的势能函数（显然是非谐振势）。由于获得等效势V`又要借助于其它的φ，所以这个过程需要通过SCF方式求解，直到{e}变化极小，可见和Hartree-Fock过程极像。求解的时候也可以通过振动基态和激发态的{φ}构成的集合作为基来像HFR方程那样求解。通过VSCF过程得到的是最佳的正则振动波函数{φ}，将之相乘就构建了零阶的N坐标的体系振动波函数。然后可以基于此通过组态相互作用（VSCF-CI）、微扰（cc-VSCF）、MCSCF等思路来获得更精确的实际振动波函数。例如做CI，那么就是每个正则坐标对应的基态/激发态的φ以全部可能的方式相乘（或者只考虑最高到第二振动激发态，类似CISD），构成展开体系振动波函数的N坐标振动基函数，然后变分求解即可（十分类似于Slater行列式展开精确N电子波函数）。VSCF-CI的结果与VSCF的结果相差不是太大，而VSCF比NMHO有很大改进。

可见，最简单的求解振动问题的模型，即NMHO模型，有两个弊端
1：用了谐振势，和实际的非谐振势有差距，特别是描述解离问题时。
2：没考虑振动模式的耦合（原因一方面也是因为用了谐振势）。
VSCF的目的在于较准确获得多原子分子的吸收频率，既考虑非谐振效应，即彻底解决了问题1；也把正则坐标之间的耦合效应等效地考虑进去，即很大程度等效解决了问题2，但只有通过CI/微扰等相关方法才能彻底解决之。

如果不考虑耦合问题，而单纯考虑非谐振效应，那么直接在正则坐标对应的势能曲线（非谐振势）下求解核薛定谔方程，就直接得到了该正则模式的振动能级和振动波函数。这叫做uncoupled-anharmonic-oscillator (UAO)模型。UAO比起VSCF有较大差距。

GAMESS-US、NWChem、Molpro等支持VSCF。只要所用的理论方法能给出能量，即能给出正则坐标上的势能曲线，那么就可以在这样的级别下做VSCF。所以VSCF绝不仅限于SCF型电子结构方法（HF/DFT），VSCF不需要所用的方法能给出高阶导数。虽然一开始需要Hessian矩阵来得到正则坐标，但可以是在低级别下的方法得到的。假设势能曲线（特别是解析的）能方便地得到，那么VSCF可以求解巨大体系的振动问题。

提供正则坐标上的势能曲线既可以做势能面扫描，或者拟合成解析曲线，也可以在平衡位置求二阶以上的导数，这样基于泰勒展开，其它点的势能就能近似地知道。因此，在平衡位置上能提供的导数阶数越高，非谐振校正也就越准确。

对于NMHO或VSCF模型，由于采用了单振动坐标近似，体系的振动能量表示为每个正则振动模式的振动能的加和，也就是说体系的振动能量由各个正则模式的振动量子数来决定。而对于真实情况，体系的振动能量只取决于一个数，即体系振动量子数。光的吸收实际上是不同体系振动能级之间的跃迁，但我们之所以也可以通过正则模式的振动能级来讨论红外吸收，是因为单振动坐标近似往往是一个比较好的近似，特别是在VSCF情况下。

更多的关于VSCF的讨论、计算公式见Acc. Chem. Res. 1986, 19, 202-208和http://homepages.abdn.ac.uk/m.m.law/pages/research/rovib99/bk/node50.html

合频（或称组频，Combination Bands）峰是指的两个基频频率相加的峰，出现它的原因在于吸收前后体系从振动基态跃迁到了振动激发态，而这个过程又可以很好地近似描述为那两个正则振动模式的跃迁。

分子的核振动是复杂的，看似是无规律的，若将它描述成多个振动模式的叠加，每个振动模式有特定的频率，那么这些振动模式就被称作正则模式。正则模式使得分子的振动变得易于考察。正则模式只能说是一种等效的描述振动的方式。我们说体系有3N-5（线性分子）或3N-6（非线性分子）个正则模式，这也假定了分子的振动波函数可以近似写为正则振动波函数的乘积。

在实际获得正则模式的时候，都是基于谐振模型，即对角化质权Hessian矩阵得到正则坐标（本征矢）。之后，是使用谐振势、UAO还是VSCF方法，那都是随意的。只不过，考虑了非谐振势、以VSCF方式等效考虑耦合，可以让各个正则模式对应的振动能级、振动波函数更合理。

Fermi共振是指一个振动的基频与某个合频或某一个振动的倍频接近时，使得二者的吸收频率都发生了变化。这是由于靠正则振动波函数乘积所得的相应两个N电子振动波函数间的相互作用所致，可认为是VSCF-CI矩阵的相应的非对角元此时不可忽略。