谈谈谐振频率校正因子

 

前言：本文对量子化学计算振动问题涉及的频率校正因子的概念做简要的介绍，这是改进谐振近似下算的振动光谱的常用做法，对于改进ZPE的计算精度也有宜。详细、深入得多的介绍、大量实例，以及振动分析方面非常全面的背景知识，在北京科音中级量子化学培训班（http://www.keinsci.com/KBQC）里笔者有专门系统的讲授，并且此培训的振动分析一节还十分深入、全面讲解非谐振计算的知识，比起本文介绍的谐振频率乘以频率校正因子的简单做法可以得到明显高得多得多的精度，欢迎参加！

 

1 频率校正因子的基本知识

对化学体系做振动分析时，通常理论计算的频率是谐振频率，由于忽视了非谐振（或称非简谐）效应，而且理论方法结合特定的基组对势能面本身的描述也有误差，所以直接基于谐振频率计算基频（振动基态到第一振动激发态的跃迁能量对应的光的频率）、零点能（ZPE）都可能有较大误差，对熵、从0K升温对焓的改变量也有一定误差。而将频率校正因子（后文简称校正因子，scale factor）乘在理论计算出的谐振频率上，此时对应的基频以及振动对热力学校正量的贡献就变得更准确了。例如一些杂化泛函算出来的较高频率的基频的平均绝对误差往往能到100多cm^-1，而校正后能降到几十cm^-1。

校正因子分为多种，通常有
(1)用于获得准确基频（fundamental）的校正因子。通常对于所有频率都是相同的校正值，但有些人为了校正得更准确，把频率划分成不同范围（比如高频和低频范围），分别给出不同的校正因子。
值得一提的是，由于忽略了非谐振效应会倾向于高估频率，导致杂化泛函的基频的校正因子都小于1.0，而且HF成份越高校正因子越接近0.9，不了解HF成份的话参看《不同DFT泛函的HF成份一览》（http://sobereva.com/282）。而纯泛函由于对势能面描述的系统性误差和谐振近似的系统性误差恰好抵消得较充分，导致此时基频校正因子很接近1，因此就没必要考虑频率校正因子了。
(2)用于获得准确谐振频率的校正因子。注意谐振频率不是实验可观测的而只是个概念，但可以通过一定手段将实验数据转化得到
(3)用于获得准确的ZPE的校正因子
(4)用于获得准确的焓变ΔH=H(T)-H(0)的校正因子。注：H(0)=U(0)=电子能量+ZPE
(5)用于获得准确的熵的校正因子
这5种校正因子各不相同，不要随意混用。ΔH和熵的数值是依赖于计算时设定的温度的，其校正因子也因此受温度影响。

由于谐振子模型下ZPE是频率的线性函数，即ZPE=∑hν/2，所以第(3)类校正因子乘到原始频率上还是乘到基于原始频率算出的ZPE上效果都一样。但是对于(4)和(5)类校正因子，要分别乘到原始频率上再算ΔH和熵，而切不要用原始频率算完了之后再对所得的ΔH和S乘上校正因子。

还有一点要说明，由于ZPE=∑hν/2，表面上看仿佛只要对振动频率ν乘上校正因子成为准确基频后ZPE就应该准确了，因此(1)和(2)类校正因子仿佛应当是相同的，这是误解！因为这种ZPE的计算公式是对于谐振模型而言的，既然校正成了实际的非谐振频率，那么当然ZPE也得按照非谐振子来算才行。严格的非谐振的ZPE必须直接做非谐振计算才能得到。ZPE的校正因子则实际上是先乘到谐振频率上将之转化成一种特殊频率，然后再按照ZPE=∑hν/2来算的时候，得到的结果就能正好和实际的ZPE很接近。

 

2 频率校正因子的获取

校正因子一般都是用先人拟合好的，通过向相应类型实验数据拟合来得到。用于拟合的分子少则几种，多则上百种。常用的校正因子可以在这些地方得到：
1 http://comp.chem.umn.edu/freqscale/version3b2.htm。这是Truhlar组拟合的一大堆校正因子，相当全面，包括前述的(1)、(2)、(3)类校正因子。不过很遗憾的是没有给出ΔH和熵的校正因子。

2 http://cccbdb.nist.gov/vibscalejust.asp。这是Computational Chemistry Comparison and Benchmark DataBase(CCCBDB)数据库中的一页，给出的是基频的校正因子。

3 J. Phys. Chem., 100, 16502-16513 (1996)中的表10，给出的是前述(1)、(3)、(4)、(5)类校正因子。注意表10的ΔH和S的校正因子是298.15K下的，其它温度下的见文中图2。此文作者后来在2007年又发了一篇文章J. Phys. Chem. A, 111, 11683-11700 (2007)，给出了更多方法和基组（主要是pople基组）下的校正因子。

4 J. Phys. Chem. A, 108, 9213-9217 (2004)。给出的是HF、MP2和B3LYP结合Dunning相关一致性基组的各类校正因子。

5 J. Comp. Chem., 33, 2380 (2012)。给出了一些泛函结合pc系列基组的高频、低频、焓、熵、ZPE的校正因子。

6 Theor. Chem. Acc., 105, 413 (2001)。给出了HF、几种泛函和MP2结合Sadlej pVTZ基组时的高频和低频校正因子

网友Liyuanhe整理了一个很好的表格，汇总了各种来源（包括上述的）的校正因子，查阅方便，值得下载留存，见http://bbs.keinsci.com/forum.php?mod=viewthread&tid=3805。

虽然这些来源拟合校正因子过程所用的测试集和方法有异，但对于同种级别，结果是接近的，因此它们是有可比性的。还有一些校正因子零零碎碎地在一些文献中报道，用得较少，可靠性难说。

 

3 关于频率计算级别的选择

注意振动频率和热力学校正量的计算绝对不是用的理论级别越高越好。CCSD(T)/cc-pVTZ可能还不如B3LYP/6-31G*算出来的频率准确，因为误差既来自于理论方法给出的力常数的误差，也来自于忽略了非谐振效应的误差。在某些低水平下，可能恰好这两部分抵消得好，导致结果不错，纯泛函算的基频几乎不需要校正就是这个原因。而如果盲目地用很高的计算级别，虽然基本消除了理论方法本身的误差，但非谐振效应缺乏较好的抵消，那么所得的基频以及热力学校正量很可能比很烂的级别算出来的还糟。另外还有关键的一点，就是在使用了校正因子后，不同方法间的差距会很大程度上拉平，原本差的方法在校正后结果可能不错。这正是为什么热力学组合方法G3，虽然目标是QCISD(T)的高精度，却采用HF/6-31G*（经过校正因子校正）这样表面看上去如此烂的级别来算ZPE。所以说，用很高级别的方法去做振动分析是⑨的行为，何况计算能量的二阶导数还那么费时，纯属白耽误功夫。半经验方法也有对应的校正因子，但是即便是经过校正后，误差还是较大，大于从头算方法很多，所以除非迫不得已不宜用半经验方法算频率。基于分子力场也可以算频率和热力学校正量，如果分子力场选得合适，那么结果可能很不错，且往往不需要做校正，例如MMFF94和MM3力场对于典型有机分子基频平均误差在60cm-1左右。

通常来说，使用B3LYP/6-31G*结合频率校正因子就能得到不错的计算结果（最起码对于有机体系和大多数无机体系而言），计算量小，而且比起其它多数方法（哪怕是更昂贵的）经过校正后的结果更好，所以完全不需要考虑其它的。

虽然频率校正因子校正后的B3LYP/6-31G*是通常推荐的做振动分析的级别，但是对有些体系，可能在B3LYP/6-31G*级别下的势能面很糟，诸如色散主导的弱相互作用显著的情况，进而可能由此导致优化出来的结构、算出的频率和热力学校正量都很差（即便已经经过频率校正因子校正）。此时就应该选择更合适的方法做结构优化和振动分析。例如对于弱相互作用对结构影响显著的时候可以用B3LYP-D3(BJ)和ωB97XD，它们都带了色散校正，色散校正的意义看《DFT-D色散校正的使用》（http://sobereva.com/210）、《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）、《谈谈量子化学研究中什么时候用B3LYP泛函优化几何结构是适当的》（http://sobereva.com/557）、《乱谈DFT-D》（http://sobereva.com/83）。

各类校正因子对于所用基组都不太敏感（除非是和很烂的基组相比），所以如果你用的基组没有对应的校正因子，可以直接用此方法在相近的基组下的校正因子。特别值得一提的是弥散函数对校正因子的影响微乎其微，所以假设你的基组用的是6-31+G*，但是找不到相应的校正因子，却能找到6-31G*下的校正因子，那么就可以放心地用6-31G*的校正因子。另外，对于同一种方法，无论是校正前还是校正后，都不代表越大的基组下结果越准确，而通常是差不多（除非基组烂到3-21G这种层次），因此一般没必要用2-zeta以上的基组、带弥散函数做振动分析。经常有人问用混合基组时怎么选择校正因子，实际上此时没有完美的做法，一种做法是索性根本不考虑校正因子，还一种也有道理的做法是使用占原子数最多的基组的频率校正因子。而如果你只对某个振动频率感兴趣，而这个振动模式涉及的原子基本上用的都是基组A，那可以使用A基组对应的频率校正因子。

要提醒的是，振动分析所用的级别和优化几何结构的级别必须严格一致，而得到的热力学校正量是可以加到其它更高水平方法算的电子能量上的。比如说，要算一个体系某温度下的气相下的自由能，那么可以用B3LYP/6-31G*去优化体系然后在考虑校正因子的情况下做振动分析，得到ZPE、ΔH和熵，再用很好的方法诸如revDSD-PBEP86-D3(BJ)/def2-QZVP去获得这个结构下的电子能量E。最后体系的自由能G=E+ZPE+ΔH-T*S。

实际上，自行拟合频率校正因子不是难事，见《自行拟合基频校正因子与ZPE校正因子的简单方法》（http://sobereva.com/391），用自行拟合的显然比凑合借用一个更有信服力。

注意校正因子并非总能很好解决振动问题。用校正因子后只是整体统计结果变好了，并不代表对每种体系的每个频率都能得到不错的精度，比如很多级别下对于氟分子的振动频率算得都不咋样，靠频率校正因子也不灵。并且对于非谐振效应很强、很低频率的情况，用了校正因子后结果肯定也还是不好，只能直接去做非谐振计算才行，更多相关讨论见《基频频率校正因子实际效果测试》（http://sobereva.com/390）。此外，合频、倍频跃迁在谐振近似下是禁阻的，而且谐振近似下它们的跃迁频率被严重高估，计算它们也必须做非谐振计算才行。常用的VPT2形式的非谐振计算在北京科音中级量子化学培训班（http://www.keinsci.com/KBQC）有极为全面的讲解，北京科音高级量子化学培训班（http://www.keinsci.com/KAQC）里我还讲授了更高级的振动问题计算方法如VSCF、VCI的原理和实例，欢迎参加。

 

4 在Gaussian中使用频率校正因子

2020-May-17注：下文介绍的freqchk程序已经没有任何存在意义了！因为使用《使用Shermo结合量子化学程序计算分子的各种热力学数据示例》（http://sobereva.com/552）文中介绍的Shermo程序计算热力学量省事方便得多，且功能强大得多，还可以同时指定不同的校正因子，方便至极，再也不用像下文一样折腾了！

下面简单说一下在Gaussian中使用校正因子的方式。在freq任务过程中使用scale=X就代表校正因子为X，例如
#P HF/6-31G* freq scale=0.9135 temperature=323
注意输出信息中只是在计算热力学校正量的时候对所用频率乘上了校正因子，而输出的频率还是原始的频率！所以要想得到校正后的频率，得自己手动把频率乘上校正因子。Gaussian显然是不会自动用校正因子的，即默认scale=1.0。

如果已经做过一次频率计算并且保存了chk文件，想查看另一个校正因子下的热力学数据，最好的方法是用Gaussian自带的freqchk工具，它会从chk文件中读取储存的Hessian矩阵并作分析。启动它之后，输入chk文件路径，然后依次输入
N    //不输出Hyperchem文件
323   //温度。如果输入0则为默认的298.15K
1     //压力。如果输入0则为默认的1atm
0.8905   //校正因子
[回车]
[回车]
然后立刻就看到了原始频率以及等同于scale=0.8905 temperature=323时freq任务给出的热力学数据，和Gaussian自身输出的一样。

利用freqchk，就可以反复利用不同校正因子得到较准确的ZPE、ΔH和熵，并获得自由能G=H-T*S，过程在下面说一下。以下是freq任务会输出的四个量，为了省事用字母表示，使用不同校正因子时它们都会有变化。
 Zero-point correction= A
 Thermal correction to Energy= B
 Thermal correction to Enthalpy= C   注意这一项不是ΔH，而是ZPE+ΔH
 Thermal correction to Gibbs Free Energy= D
首先，在ZPE校正因子下得到ZPE（即A）。然后在熵校正因子下得到-T*S（即D-C）。之后在ΔH校正因子下得到ΔH（即C-A）。最终将这三个量加在一起，再加上高精度方法算的电子能量E，就得到了自由能G=E+ZPE+ΔH-T*S。

在作IR、Raman、VCD这些振动光谱的谱图时，要先把频率乘上校正因子后再做图。手动一一去改Gaussian输出文件里的频率值显然太费事了。Multiwfn (http://sobereva.com/multiwfn)的绘制光谱图的功能可以方便地设定校正因子。启动Multiwfn后，载入Gaussian的freq任务的输出文件，选11进入光谱绘制功能，选IR或Raman，之后选0就可以显示光谱图。如果选14 Multiply the vibrational frequencies by a factor，就可以令不同编号范围的频率乘以指定的校正因子，重新作图后就能看到效果了。用户因此不仅可以对所有频率乘上统一的校正因子，还可以根据一些文章的思想，对于不同范围的频率乘上不同的校正因子以得到更准确结果。关于Multiwfn绘制光谱的更多信息看《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（http://sobereva.com/224）。

最后顺带一提的是，除了对频率进行校正，还有人提出了直接对力常数进行校正的方法，不过这样的方法涉及到太多参数，难搞，所以没流行起来。对于NMR、UV/Vis等其它谱的计算也有人提出过不同级别下的各种校正方法，不过都远没有像频率校正因子这样被广泛使用。