浅谈为什么优化和振动分析不需要用大基组

sobereva

浅谈为什么优化和振动分析不需要用大基组

A brief discussion on why optimization and vibration analysis do not require use of large basis sets

文/Sobereva @北京科音

First release: 2017-Aug-26   Last update: 2024-Feb-27

量化计算中有两个常识，大部分人想必都知道
(1)优化和振动分析必须在相同级别下进行
(2)在高耗时的优化+振动分析时用中小基组，而只在随后的单点能计算时才用大基组，这样做比起所有过程都用大基组耗时大为降低而精度降低很少

其中(1)的原因很好解释，因为如果振动分析用的级别（包括各种影响势能面的因素，如积分格点精度、溶剂模型等）与几何优化时不同，相当于振动分析时并不是在势能面的极小点的位置进行的（即便用的结构是之前优化得到的），此时振动分析结果显然没有意义。而(2)的原因，是因为优化和振动分析结果对基组的敏感性远低于能量计算。正好最近在QQ群里有人又提到此事，笔者遂通过一个简单的HF分子的计算数据来展现这一点。虽然不能光用一个分子就以偏概全，但很大程度上是能够以小见大的。计算使用Gaussian16完成，理论方法为B3LYP。如果缺乏基组选用的基本知识，看此文前强烈建议仔细阅读《谈谈量子化学中基组的选择》（http://sobereva.com/336）。

来看不同基组下计算的H-F键能（即E(H)+E(F)-E(HF)），优化的键长，以及谐振频率。计算分子能量时键长用的是相应基组下优化的。

首先看def2-SVP/TZVPP/QZVP这个序列的数据。按照这个顺序，基组质量是系统性地提升的，其基组尺寸类似于cc-pVDZ/TZ/QZ序列。我们通常说，优化用2-zeta带极化这个档通常就够了，从上表确实可见相对于很大基组def2-QZVP的结果，def2-SVP的键长误差才0.3%，振动频率误差才0.13%，如果要求不高的话这点误差完全可以接受，然而其键能的误差则达到26.4kJ/mol，相当于4.5%，完全不能忽略。再对比def2-TZVPP和def2-QZVP的结果，前者相对于后者的键长、振动频率的误差根本察觉不到（仅0.0003埃和0.4cm-1），但能量上还是有3kJ/mol的误差，对于高精度热力学数据计算，这种程度误差是仍是不可忽略的，而且高精度计算通常是用后HF，此时def2-TZVPP与def2-QZVPP在算能量相关问题上的差距会更大。上面的数据充分说明优化+振动分析结果对基组的敏感性远低于能量，因此前者用的基组建议比后者低一个档来节约时间。另外值得一提的是，由于做振动分析时通常还会考虑频率校正因子（详见http://sobereva.com/221），因此用大基组在精度上实际上并不会带来任何好处。虽然原理上来说在大基组优化的结构下用大基组算能量，比起在中小基组优化的结构下用大基组算能量，对于能量相关问题的计算准确度要更高，但由于中小基组下优化的误差已经不大，再加上在极小点附近势能面非常平缓（即结构的稍微变化基本不会带来什么能量变化），因此不用担心优化+振动分析用中小基组来节约时间的做法会带来能量计算精度上的损失。如果你真是想要很高精度结构和能量数据，那么优化用到def2-TZVPP就已经绝对足够了，再提升完全是白费时间，而算能量时基组则完全有必要用到def2-QZVP，特别是对于后HF计算而言。

我们再关注一下上面表中其它基组的结果。如《谈谈量子化学中基组的选择》提到的，对于一般问题没有必要给氢加极化，但是和氢有密切相关的问题则应当给氢加极化，当前就是此种情况。对于HF分子，def2-SV(P)等价于砍掉了def2-SVP上氢的p极化函数的版本，def2-TZVP等价于砍掉了def2-TZVPP上氢的d极化函数的版本。由数据可见，砍掉氢的极化对结果影响确实不小，特别是连p极化都砍掉时，哪怕是键长、振动频率的误差的增加都已经到令人恼火的程度。由于给氢加p极化对耗时增加其实不算太多（相对于给重原子增加高角动量极化函数而言），因此有必要给氢加p极化时不要吝啬，即起码用def2-SVP或6-31G**这个档。再看超恶心的连重原子的极化都没有的6-31G，误差比def2-SV(P)又明显进一步增大，键能误差已经快100kJ/mol了。对比def-TZVP和def2-TZVP的结果，虽然后者比前者昂贵很多（很大程度上是因为对F加了f极化所致），但是各方面精度都差不太多，因此def2-TZVP无福消受时用def-TZVP是很好选择。

为了直观地展现为什么几何优化和振动分析对基组敏感度低（当然前提是起码已用到def2-SVP这个档），下面在不同基组下对HF的势能面进行了扫描，如下所示。图中的能量零点是相应基组下H与F原子无限分离的情况。

由图可见def2-TZVPP与def2-QZVP的曲线已经几乎精确重合，直接反映在优化结果和振动分析结果几乎相同上。def2-SVP的曲线谷底明显与def2-TZVPP的有一段距离，所以能量误差明显，但谷底横坐标位置却几乎和def2-TZVPP的没区别，肉眼难以分辨，所以几何优化结果误差几乎可以忽略。而6-31G，不仅势阱深度离谱，肉眼都可以看出谷底位置与def2-TZVPP的有一定偏离。

势能曲线谷底位置的曲率正比于振动频率。虽然def2-SVP的振动频率误差已很小，但由上图，肉眼看过去还是容易凭直觉觉得其极小点的曲率和def2-TZVPP的有明显差异。为了更容易地对比曲率，下图把def2-SVP和6-31G的曲线进行平移，使其极小点与def2-TZVPP的相重合。

由上图可见，极小点附近def2-SVP和def2-TZVPP的曲线极为吻合，因此曲率也几乎完全相同，只有到了距离平衡位置较远的区域二者才有可察觉的差异。因此也只有计算比如倍频时，由于涉及到势能面更广的范围，两种基组下的结果才可能有不可忽略的差异。而6-31G真是糟糕透顶，不仅能量计算超烂，连极小点曲率都和def2-TZVPP的大相径庭，明显低估了曲率，因此频率也严重低估。

本文的数据清楚地体现，对几何优化和振动分析不需要太大基组，比如def2-SVP、6-31G**这样2-zeta带极化函数的档次一般就够满足需要了。“不需要太大”绝对不意味着可以偷工减料到使用连非氢原子的极化函数都没有的6-31G之流。还要注意的是，不同体系对基组的敏感性是有差异的，比如优化过渡金属的配位键键长，def2-SVP和def2-TZVP的结果差异往往挺明显。再比如，18碳环在6-31G*基组下优化的结果是有明显不合理之处的，我在《我对一篇存在大量错误的J.Mol.Model.期刊上的18碳环研究文章的comment》（http://sobereva.com/584）中专门提到了，这个时候就得用更贵的如6-311G*、def-TZVP。“2-zeta带极化函数的基组对几何优化和振动分析一般够用”主要是对普通有机体系而言的。这里说的“一般”是超过至少95%的概率，如果你实测发现用更大的如def-TZVP会导致研究结论发生显著改变，那就属于小概率的情况了，此时当然就得以更大基组的结果为准。

不光是能量，诸如偶极矩、(超)极化率、NMR等诸多与电子结构直接相关的属性对基组的敏感度也是远高于几何优化和振动分析的，因此计算它们之前的优化过程也建议用中小基组。值得一提的是计算Raman和ROA计算实际上分为两个过程(1)优化和振动分析得到正则坐标 (2)计算极化率对正则坐标的导数。由于如前所述(1)对基组要求低，弥散函数更起不到什么作用，而弥散函数对(2)的结果改进明显（原因容易理解，没弥散函数时极化率算得超烂），因此Raman或ROA在Gaussian中允许两步分开做，用户可以在(1)的时候用不带弥散的基组节约时间，而只在(2)的时候用带弥散的基组。

为节约时间，几何优化+振动分析不仅可以用比能量计算更小的基组，还可以用更便宜的理论方法，这种做法的合理性与上文讨论的类似。比如CCSD(T)算能量而DFT做几何优化和振动分析，或者前者用双杂化泛函后者用普通泛函，这都是非常流行、划算的搭配。虽然小基组优化+大基组算能量对内行人属于常识，但为了让外行读者也能明白此做法的合理性，论文里可以引用相关文章，可以引用比如《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）里介绍的笔者的Shermo程序的原文Comput. Theor. Chem., 1200, 113249 (2021) DOI: 10.1016/j.comptc.2021.113249，里面在计算高精度热力学量的时候就用了这种做法。

在J. Chem. Theory Comput. (2023) DOI: 10.1021/acs.jctc.3c00388还专门系统性测试了几何优化+振动分析用从小到大不同基组时对最终计算出的热力学量的影响，发现大多数情况对最终结果影响都小于1 kcal/mol，同样证明本文的观点。同时此文还指出存在低频模式时，反应自由能对优化+振动分析用的基组会更敏感，但如果使用了准RRHO模型计算热力学校正量代替Gaussian等程序默认用的RRHO，则敏感度会大幅降低。这充分体现出含有低频模式时，使用Shermo程序基于准RRHO模型计算自由能变的重要性。

另外，liyuanhe写了一篇《自由能计算对优化与单点级别的敏感性测试》（http://bbs.keinsci.com/forum.php?mod=viewthread&tid=6609），与本文的内容有一定相关性，建议阅读。

wzkchem5

一条君 发表于 2022-7-21 10:13
老师，比如氟化锂真空25度是什么状态呀，远低于其沸点，如果不是气态，热力学计算又是算的他真空25度的气 ...

当然是固体，因为不仅远低于它的沸点，也远低于它的熔点。这种问题不用查任何资料也能猜出来，但是如果猜不出来，得知道怎么查资料来解决。当然如果压力低到离谱（比如10^-100 Pa），那可能是气体，不过那是另一个问题了。
高斯默认算的是25度下、压强为1个大气压下的理想气体。换个说法就是25度下、逸度（注意不是压强）为1个大气压下的气体。但是高斯不能自动判断这个物质在25度、1个大气压下是不是气体，这是用户应该判断的内容。如果不是的话，高斯仍然能算出热力学数据，但不代表实验上存在的状态，高斯也不会告诉你实验上不存在这个状态，因为这超出了高斯的能力。事实上也超出了任何其他软件的能力，没有任何软件能只算一个LiF分子就告诉你它在你设定的温度、压力下是固体、液体还是气体，只有你明确让程序分别算固体、液体、气体，然后你自己比较哪个的自由能最低，才能知道正确的物态。

一条君

wzkchem5 发表于 2022-7-21 15:56
（1）对，因为一个化合物是不是离子化合物不是看它的解离产物，而是看它没解离的时候是离子成分为主还是 ...

老师，比如氟化锂真空25度是什么状态呀，远低于其沸点，如果不是气态，热力学计算又是算的他真空25度的气态能量吧，这块不太明白，谢谢老师

wzkchem5

一条君 发表于 2022-7-21 07:57
老师，可能基础不太好，（1）那气态LIF是离子化合物吧。另外您的意思是，离子化合物都是A+(g) + B-(g)的 ...

（1）对，因为一个化合物是不是离子化合物不是看它的解离产物，而是看它没解离的时候是离子成分为主还是共价成分为主
（2）实验已知沸点的物质，自己查沸点判断；实验不知道沸点的物质，查类似物质的沸点判断。

一条君

wzkchem5 发表于 2022-7-15 15:28
因为气相中Li+(g) + F-(g)的能量反倒高于Li(g) + F(g)，这和晶体里、水溶液里是相反的。虽然这样可能很反 ...

老师，可能基础不太好，（1）那气态LIF是离子化合物吧。另外您的意思是，离子化合物都是A+(g) + B-(g)的能量高于A(g) + B(g)，所以键解离能的定义才是解离成气态原子的能量对吗；A(g) 代表独立存在的原子，A(g) + B(g)对应您说的间隔足够远的两个原子
（2）高斯热力学计算 默认真空25度，算的是气态，那这个物质在真空温度下会是气态吗，谢谢老师

wzkchem5

一条君 发表于 2022-7-14 04:25
老师，共价键和离子键键能算的方法一样吗，比如HF和LiF，都是E(H)+E(F)-E(HF)或E(Li)+E(F)-E(LiF)吗，为什 ...

因为气相中Li+(g) + F-(g)的能量反倒高于Li(g) + F(g)，这和晶体里、水溶液里是相反的。虽然这样可能很反直觉，但是事实就是这样的，F原子看起来氧化性那么强，Li原子看起来还原性那么强，但是间隔足够远的F原子并不能氧化Li原子，之所以LiF在晶体或溶液里是离子化合物，完全是由于晶格能或溶剂化能的作用。

sobereva

一条君 发表于 2022-7-14 11:25
老师，共价键和离子键键能算的方法一样吗，比如HF和LiF，都是E(H)+E(F)-E(HF)或E(Li)+E(F)-E(LiF)吗，为什 ...

片段的状态选取有任意性，说键能的时候总要说清楚均裂还是异裂

chands

一条君 发表于 2022-7-14 11:25
老师，共价键和离子键键能算的方法一样吗，比如HF和LiF，都是E(H)+E(F)-E(HF)或E(Li)+E(F)-E(LiF)吗，为什 ...

键解离能的定义如此，这里LiF的键解离能是指LiF气态分子：LiF(g) = Li (g) + F (g)的焓变。
如果要衡量LiF(s)的离子键，用晶格能LiF(s) = Li+(g) + F-(g)

一条君

老师，共价键和离子键键能算的方法一样吗，比如HF和LiF，都是E(H)+E(F)-E(HF)或E(Li)+E(F)-E(LiF)吗，为什么氟化锂不是E(Li+)+E(F-)-E(LiF)呢，谢谢~

sobereva

granvia 发表于 2020-1-19 12:51
感觉测试体系覆盖面不是很广泛，说服力未必可靠。非共价体系，柔性构象，过渡金属配合物，过渡态等等都应该 ...

这种测试就是专门paper的事了，不是闲暇时间写个帖子能覆盖的

根据我阅读的大量已有的研究文章和长年实际研究经验，本文对各种类型体系结论都是一致的，足够以小见大。只不过算弱相互作用的时候应当用中档3-zeta进行优化和振动分析（如6-311G**、def-TZVP），而2-zeta偏小，这在我其它文章里说过。

granvia

感觉测试体系覆盖面不是很广泛，说服力未必可靠。非共价体系，柔性构象，过渡金属配合物，过渡态等等都应该系统测试一下吧

去离子水

算红外吸收强度的时候，需要算偶极矩对正则坐标导数。偶极矩的计算对弥散函数的要求高。那么算红外光谱的时候是不是最好加上弥散函数呢？

sobereva

1130240115 发表于 2017-8-27 19:54
我也遇到过类似的问题，使用M06-2X优化O2，发现优化得到的O2的键长比实验值小0.2埃左右，是否是实验值有问 ...

看自旋多重度是否设成了应有的3，而且也看基组是否合理
双原子分子的键长可以通过微波谱很准确地得到，而且这又是非常常见的分子，如果你是从NIST上查的双原子光谱常数数据，那不会有问题，除非你是取错了该取的数。

1130240115

我也遇到过类似的问题，使用M06-2X优化O2，发现优化得到的O2的键长比实验值小0.2埃左右，是否是实验值有问题？

sobereva

Daniel_Arndt 发表于 2017-8-27 13:12
我记得一篇JOC上提到过，对于[1.1.1]propellane的那根与三次旋转轴重合的碳-碳键的键长，M06-2X给出的结 ...

我对于实验值其实比较怀疑，几种泛函，包括双杂化，以及泛函横测文章中对有机分子结构优化表现的最好的PBE0在def2-TZVP下那个键都比1.6埃小0.3埃左右或者更多
用CCSD/def2-TZVP优化结果是1.5527