简谈Gaussian里找过渡态的关键词opt=TS和QST2、3

sobereva

简谈Gaussian里找过渡态的关键词opt=TS和QST2、3

文/Sobereva@北京科音  2019-Feb-3

Gaussian里找过渡态的关键词有三：opt=TS、opt=QST2、opt=QST3。经常在网上看到有很多初学者居然把QST2当做找过渡态默认用的关键词，明显是被一些坑爹言教程毒害的后果。这里就抽一点时间非常简单地说说常用的opt=TS关键词，以及为什么我很不建议在一般情况下用QST2、QST3。本文并不举实际计算例子，我认为学习找过渡态必须通过大量例子才能培养出感觉，在北京科音初级量子化学培训班（http://www.keinsci.com/workshop/KEQC_content.html）我会给出涵盖各种类型反应的过渡态搜索例子，在相当于中级班的北京科音中级量子化学培训班（http://www.keinsci.com/workshop/KBQC_content.html）里我会讲得明显更深入、给出更复杂更高难的例子，并介绍诸多高阶的经验技巧。推荐参加这两个培训。关于过渡态搜索的具体算法，本文也不做介绍，感兴趣者参看各种方法的汇总《过渡态、反应路径的计算方法及相关问题》（http://sobereva.com/44）。


1 关于opt=TS

一般情况下我都建议用opt=TS找过渡态，其算法和优化极小点结构非常类似，只不过它要求优化的是势能面的一阶鞍点。具体来说一般是用opt(TS,calcfc,noeigen)关键词，calcfc代表精确计算初始结构的Hessian矩阵，noeigen代表优化过程中不对每一步做Hessian矩阵本征值数目的检测（这个要求太严而且没太大必要）。用opt=TS只需要提供一个初猜的过渡态结构就够了，很方便。只要初猜结构和实际过渡态比较接近，而且计算级别和模型合理，一般都能收敛到实际过渡态。一般建议再同时写上freq关键词使得优化结束后自动做个振动分析。

怎么验证过渡态找没找对？如果优化出来的过渡态结构在期望的反应物和产物之间，同时有且只有一个虚频，一般能有6成把握认为过渡态找对了。如果通过gview查看，发现仅有的那个虚频的振动方向对应反应坐标方向，那8成可以认为过渡态找对了；如果走IRC还能连通自己期望的反应物和产物结构，那么说明100%找对了。走IRC的相关信息看《在Gaussian中计算IRC的方法和常见问题》（http://sobereva.com/400）。

opt=TS常见问题如下：
·如果优化出现震荡、难收敛，但是用gview观看结构变化感觉收敛趋势对，可确认优化过程是在往过渡态模样的结构走，最后一帧结构有过渡态的大致模样，按照《量子化学计算中帮助几何优化收敛的常用方法》（http://sobereva.com/164）适当处理即可。
·如果是找出来的过渡态不仅存在一个对应反应坐标方向的虚频，还存在额外的虚频，按照《Gaussian中几何优化收敛后Freq时出现NO或虚频的原因和解决方法》（http://sobereva.com/278）消除多余虚频即可。
·如果优化不收敛或者报错，且最后一帧结构完全不像实际过渡态，或者虽然任务正常结束了但是根本都没有虚频，或者虽然只有一个虚频但是并不是对应当前反应方向的，那就是说明自己提供的初猜结构不合理（也有可能是计算级别不合理，导致当前级别势能面上根本没有对应实际过渡态的一阶鞍点），那就再反复尝试其它自己觉得靠谱的初猜结构，或尝试其它自己觉得合理的计算级别。

找过渡态和优化极小点一样，不需要用到太昂贵的计算级别，见《浅谈为什么优化和振动分析不需要用大基组》（http://sobereva.com/387）。找过渡态和优化极小点用的计算级别应当一致，因为并没有用不同级别的合理理由。对于溶剂下的反应，特别是极性溶剂下，优化过渡态的时候最好就带着隐式溶剂模型，因为有的时候可能有定性影响。对于某些反应，溶剂可能直接影响反应机理，此时用隐式溶剂模型时还得同时用显式溶剂模型。

哪怕是有丰富经验的研究者，也没有那个敢打包票说看见一个反应就能保证一次就找到过渡态。找过渡态的难度远大于优化极小点，因为对初猜结构要求高。显然初学者不应当因为对于一个反应尝试了几个过渡态初猜都失败就简简单单放弃，或者投奔QST2。如果体系比较大，每次试初猜都比较耗时，可以先用更便宜但对于当前体系也能给出定性合理结果的计算级别，等收敛到了期望的过渡态，再把这个结构当做过渡态初猜，结合用原本打算用的更好的级别去进一步准确优化过渡态。其实对于大部分有机反应（也有很多例外），用PM6、PM7这样的极便宜的半经验方法就可以给出定性合理的过渡态结构，由于很便宜，因此建议同时加上calcall增加收敛成功的几率。另外也可以考虑尝试用现在很火的Grimme的GFN-xTB，普适性比半经验方法更好而且速度飞快，见《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）。

有时候不好确定过渡态初猜结构的时候可以借助柔性扫描，参看此文柔性扫描部分的讲解和例子：《详谈使用Gaussian做势能面扫描》（http://sobereva.com/474）。


2 关于QST2、3

网上有一堆乌七八糟的Gaussian学习资料，还有官方的exploring那本书的第3版，给初学净拿QST2举例子，对opt=TS居然反倒不怎么提，导致初学者以为Gaussian里找过渡态本来就应该用QST2，有经验的人都知道这真是赤果果地在坑初学者。相对于最常用的opt=TS，用QST2有很大弊端：
(1)使用QST2必须同时提供要研究的基元过程的反应物和产物的结构（不需要很精确也不需要非得事先优化，只要模样差不多对就可以）。这比起opt=TS只需要提供一个过渡态初猜结构显然在准备工作上更麻烦。而且有些情况你根本都不好估计其中一方是什么结构。比如双分子反应（或单分子裂解反应），你往往根本就不知道反应复合物（或产物复合物）长什么样子，得找出过渡态之后走完IRC才知道，如果你试图去猜的话，这比起你用opt=TS且直接给个合理的过渡态初猜的难度通常更大。
(2)用户必须确保提供的反应物和产物结构间原子顺序完全对应，否则任务一定失败，而保证这点往往很麻烦。虽然gview里提供了connection工具可以辅助你去把两帧结构里的原子顺序进行自动匹配，但是也往往匹配不成功，或者虽然表面上告诉你匹配成功了但其实还是有原子根本没有正确匹配，你还得手动根据两帧的结构图和原子序号通过反复交换原子使顺序一致，相当费事！
(3)QST2本身不是不需要过渡态初猜结构，只不过QST2计算过程中过渡态初猜结构是程序根据你提供的反应物和产物坐标根据插值产生出来的（优化过程中QST2算法还会根据反应物和产物结构对优化方向进行引导）。这样简单粗暴、机械化产生的过渡态结构，怎么可能比起有基本化学直觉的人摆出来的过渡态初猜结构更合理？QST2用多了你就会发现，QST2自动产生的初猜往往很差劲，虽然有时候很烂的初猜结构最终还是收敛到了实际过渡态，但是达到收敛所需的步数明显高于你用opt=TS并直接提供一个基本靠谱的过渡态初猜，因此浪费了很多计算时间。还有的时候，哪怕你提供的反应物和产物结构基本和实际相符，QST2也会给出一个一看就离谱的初猜结构，导致根本不可能收敛到实际过渡态去（即便QST2还会利用反应物和产物结构来引导优化也于事无补）。
(4)Gaussian的QST2代码本身也往往有问题，有时候明明输入文件完全正确，但是算起来却会莫名其妙报错。

老有一些初学者净以为QST2是万灵药，好像给了反应物和产物结构就一定能找出对应的过渡态，明显这是极度荒唐的，QST2可能哪有这么聪明？甚至于，实验上发现反应物A最终转化到产物B，而且一看就不是一个基元过程，居然有的初学者也拿A和B的结构直接输到Gaussian里跑QST2。这么做，整个复合反应中哪个过渡态的结构都优化不出来！PS：目前倒是也有一些自动找过渡态的算法，比如GSM、GRRM等等，思想、流程都远比QST2要复杂得多，可参考《寻找势能面上全部极小点、过渡态的算法GRRM简介》（http://sobereva.com/319）、WIREs Comput. Mol. Sci., e1354 (2017)、J. Chem. Theory Comput., 13, 5780 (2017)。

还有初学者用QST2的理由是“opt=TS算不出来”。我特别讨厌“算不出来”这种含糊不清的描述，因为等于什么都没说。要知道，只要恰当地使用opt=TS，什么反应的过渡态都能算出来！原理上，opt=TS找不出来的过渡态，用QST2也绝对不可能找到。如果你已经基于基本合理的过渡态初猜结构、在合理的设定下正确地做opt=TS任务，最终也没找出来过渡态，那么指望用QST2就能找出来你想要的过渡态，概率等于连续三次单抽都出UR。

QST2也不是一无是处，在恰当的时候使用也能帮助你更方便地找过渡态，但前提必须是你已经能够熟练使用opt=TS找过渡态，而且已经充分了解QST2的思想和优缺点。可以用QST2的情况有下：
(1)虽然你的头脑里知道要找的过渡态大概长什么样，自己在gview里慢慢摆也能摆出来其合适的初猜，但是涉及的几何变量太多，摆起来太花时间或者确实就是不好操作，与此同时反应物和产物的大致结构比较好确定、也容易确保原子顺序对应，那么可以先尝试用QST2，没准直接能优化到过渡态，从而节约了你反复尝试摆合理过渡态初猜的时间。
(2)对于一个基元反应，限于你的化学直觉，实在不知道过渡态应该长什么样，你可以用QST2碰运气，但是绝对不要抱太大期望。

值得一提的是，QST2计算的时候不需要像opt=TS那样需要由用户提供初始的Hessian，即哪怕calcfc都可以不写，因此对于Gaussian里一些没二阶解析导数的方法如CCSD也可以用（当然了，在CCSD级别下找过渡态一般是愚蠢的做法，贼贵而且没必要），也可以避免一些后HF方法在计算解析Hessian时的高昂代价。不过，QST2关键词里也写上calcfc（或者更贵也更好的recalc=x、calcall）有助于增加收敛到过渡态的几率。
PS：其实opt=TS也不是必须得明确提供初始的Hessian，比如用opt(TS,modRedundant,noeigen)或者opt(TS,gediis,noeigen)的话也可以不用calcfc，但当然收敛成功的几率要明显低于opt(TS,calcfc,noeigen)。

至于QST3，不仅像QST2一样需要提供反应物和产物结构，还得再提供过渡态初猜结构，麻烦得要死，而收敛到过渡态的成功几率并不比直接用opt=TS强多少，基本没有实际价值，因此无需考虑。

opt=TS、QST2、QST3找出来的过渡态结构在原理上是没有差别的。因此，不管用的哪种方法去找，只要最终的结构通过前面提到的方式验证可以确定找对了，那就绝对没有必要再用另外的方法再去计算一次。

zsu007

社长，学习了！

hakuna

好敬业，虽不怎么用G，还是认真学一下！

chemzhh

谢谢sob的分享，平时很多任务都是寻找过渡态，sob的这篇文章也明白了很多时候的困惑应该如何解决。个人的体会关键是正确的化学直觉，能够很好的给出初始结构，理论上初始结构给得好，很快就会收敛了，这样也会极大的节省计算资源。但这一点有些时候太难了，我有些时候通过柔性扫描解决，有些时候通过固定相似过渡态的核心结构，先优化外围结构，然后用作初始结构。另外感觉很多时候用gaussview画初始结构的时候，画不出自己想要的结构，不知sob是否有这方面的文章或者有兴趣写一下这方面的内容，谢谢。新年快乐！

chemzhh

另外也可以考虑尝试用现在很火的Grimme的GFN-xTB，普适性比半经验方法更好而且速度飞快，见《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/42）。
这篇博文的链接应该是http://sobereva.com/421

sobereva

chemzhh 发表于 2019-2-5 16:19
另外也可以考虑尝试用现在很火的Grimme的GFN-xTB，普适性比半经验方法更好而且速度飞快，见《将Gaussian与G ...

已改

sobereva

chemzhh 发表于 2019-2-5 16:16
谢谢sob的分享，平时很多任务都是寻找过渡态，sob的这篇文章也明白了很多时候的困惑应该如何解决。个人的体 ...

这些内容文字不好表述，得举一大堆例子，我都是在科音的量化培训班里通过实例详细讲

kekexili_08

QST2和QST3都是不错的方法， 关键是会不会用。 和opt=ts一样， 对结构的直觉最重要。 曾经有人说做QST2必须优化反应物产物， 我提出不是必须， 被人狠批。 事实是我用qst2做TS， 绝大多数给的结构都是没优化过的。 甚至有些时候给优化过的反应物和产物， 反而是不可行的。
再比如， IOP（5/13=1） 我也觉得是好用的。 现在G16里面不能用， 很多结构优化不容易做了。 而用g09，用iop（5/13=1) 轻松搞定。

sobereva

kekexili_08 发表于 2019-3-5 09:37
QST2和QST3都是不错的方法， 关键是会不会用。 和opt=ts一样， 对结构的直觉最重要。 曾经有人说做QST2必须 ...

一个典型的就是HCN异构化，如果结构是精确优化成直线的，直接就失败，反倒还非得人为先把结构扭曲一下
本来QST2的初始结构就是靠插值猜的，不要求多高精度，所以自然也就不需要两端的结构非得先做优化，非要说必须优化反应物和产物结构的人没正确理解QST2的思想。

XZDZ

老师，请问对于一个反应，要怎么才能确定其是否有过渡态存在呢？目前在计算的一个结构，一直找不到过渡态，也不知道是自己的方法不对，还是真的这个反应就没有过渡态

sobereva

XZDZ 发表于 2019-3-11 10:24
老师，请问对于一个反应，要怎么才能确定其是否有过渡态存在呢？目前在计算的一个结构，一直找不到过渡态， ...

欢乐多

老师，您好，下图2是A+B复合物通过构象搜索、优化（B3LYP/6-31G* em=GD3BJ opt int=fine）、单点能计算得到的玻尔兹曼稳定构象。为了让构象搜索得到的玻尔兹曼稳定构象能够朝发生反应的方向去，在两个反应位点上进行了柔性扫描，图1是与搜索过渡态、跑IRC、优化（opt=calcfc  b3lyp/6-31g(d) em=gd3bj int=fine）反应物得到的A+B的构象基本一致，那么在计算能垒的时候，以图2构象体系的能量为参考呢，还是以图1构象体系的能量为参考呢？

sobereva

欢乐多 发表于 2020-6-16 10:46
老师，您好，下图2是A+B复合物通过构象搜索、优化（B3LYP/6-31G* em=GD3BJ opt int=fine）、单点能计算得到 ...

从TS跑IRC，直接连通哪个就把哪个当这一步的反应物

欢乐多

sobereva 发表于 2020-6-17 11:16
从TS跑IRC，直接连通哪个就把哪个当这一步的反应物

老师，第一步基元反应从TS跑IRC得到的产物（已优化）与第二步基元反应从TS跑IRC得到的反应物（相同基组优化）构象不一样，这个中间体的能量哪个为准？

欢乐多

sobereva 发表于 2020-6-17 11:16
从TS跑IRC，直接连通哪个就把哪个当这一步的反应物

老师您好，在新加进去G后，TS3的能垒是与反应物的G+F能量的能量差，相对能量差怎么与A+B联系起来呢？
谢谢