发布时间 : 星期三 文章寻找过渡态的方法(ts qst2 qst3)更新完毕开始阅读986f47c6852458fb760b5641
过渡态简单点说就是化学反应中瞬间形成的高自由能的不稳定复合物，形象点说就是就是反应势能面上的一个鞍点，所谓的找过渡态也就是要找到这个反应势能面上的这个鞍点，这里以DA反应作为例子简单说下我用gaussian找过渡态的经验。
Gaussian里面找过渡态主要有3种方法：TS、QST2以及QST3。这里套用gaussian手册里的说明。TS是进行过渡态而不是局部最小值的优化计算。QST2是使用STQN 方法寻找过渡态结构。这个选项需要输入反应物和产物分子结构，先后通过两组连续的标题和分子说明部分定义。注意在这两个结构中的原子顺序必需一致。TS 不能和QST2 合用。QST3是使用STQN 方法寻找过渡态结构。这个选项需要输入反应物，产物和最初的过渡态结构，先后由三组连续的标题和分子说明部分定义。注意在这三个结构中的原子顺序必需一致。TS 不能和QST3 合用。接下来贴出几个例子来说明下TS、QST2和QST3的输入文件该怎么写
%chk=c2h4_c4h6_ts %nproc=1
%mem=200mb
# opt=(ts,calcfc) freq=noraman b3lyp/6-31g* Title Card Required 0 1
C -0.37960200 1.41026100 0.50972400 H -0.06409400 1.04025400 1.48028400 H -0.26589500 2.48067200 0.40089100 C -1.26020300 0.70558300 -0.28508700 H -1.84635800 1.22289600 -1.04415700 C -1.26023800 -0.70552500 -0.28508700 H -1.84641800 -1.22280900 -1.04415900 C -0.37967300 -1.41024700 0.50972300 H -0.06414500 -1.04025600 1.48028200 H -0.26601500 -2.48066300 0.40088700 C 1.45648200 0.69084700 -0.25401700 H 1.98402400 1.24692000 0.51079100 H 1.29283800 1.24312300 -1.17152100 C 1.45645300 -0.69091000 -0.25401800 H 1.29278400 -1.24317800 -1.17152200 H 1.98396500 -1.24700700 0.51079200 qst2
%chk=c2h4_c4h6_qst2 %nproc=1
%mem=200mb
# opt=(calcfc,qst2) freq=noraman b3lyp/6-31g* Title Card Required 0 1
C -3.17791906 1.13289796 -1.01187682 H -2.64475531 0.20519304 -1.01187682 H -4.24791906 1.13289796 -1.01187682 C -2.50264475 2.30787525 -1.01187682
H -3.03580850 3.23558018 -1.01187682 C -0.96264475 2.30787525 -1.01187682 H -0.42948100 3.23558018 -1.01187682 C -0.28040135 1.13693077 -1.01187682 H -0.80804382 0.20607449 -1.01187682 H 0.78957977 1.14328807 -1.01187682 C -2.59776413 0.98894279 1.10133233 H -3.14031063 0.06692559 1.12200992 H -3.12699887 1.91846375 1.07302387 C -1.24260873 0.98347968 1.11087410 H -0.70006223 1.90549689 1.09019651 H -0.71337398 0.05395873 1.13918256 Title Card Required 0 1
C 1.41432400 0.09460100 0.31826800 H 1.42370000 0.10775400 1.42899000 H 2.47445600 0.12780100 0.00402400 C 0.66866000 1.30146700 -0.16992700 H 1.26438900 2.14276400 -0.50621500 C -0.66890700 1.30134300 -0.16992700 H -1.26479300 2.14252900 -0.50621700 C -1.41434600 0.09433800 0.31826600 H -1.42372700 0.10749400 1.42898900 H -2.47448400 0.12734000 0.00402000 C 0.77040400 -1.21422100 -0.17552700 H 1.14229400 -2.05288100 0.44176200 H 1.12544800 -1.41421500 -1.20543600 C -0.77017000 -1.21437000 -0.17552800 H -1.12517200 -1.41443600 -1.20543700 H -1.14189800 -2.05310200 0.44176200 qst3
%chk=c2h4_c4h6_qst3 %nproc=1
%mem=200mb
# opt=(calcall,qst3,noeigentest) freq=noraman b3lyp/6-31g* Title Card Required 0 1
C -3.17791906 1.13289796 -1.01187682 H -2.64475531 0.20519304 -1.01187682 H -4.24791906 1.13289796 -1.01187682 C -2.50264475 2.30787525 -1.01187682 H -3.03580850 3.23558018 -1.01187682 C -0.96264475 2.30787525 -1.01187682 H -0.42948100 3.23558018 -1.01187682
C -0.28040135 1.13693077 -1.01187682 H -0.80804382 0.20607449 -1.01187682 H 0.78957977 1.14328807 -1.01187682 C -2.59776413 0.98894279 1.10133233 H -3.14031063 0.06692559 1.12200992 H -3.12699887 1.91846375 1.07302387 C -1.24260873 0.98347968 1.11087410 H -0.70006223 1.90549689 1.09019651 H -0.71337398 0.05395873 1.13918256 Title Card Required 0 1
C 1.41432400 0.09460100 0.31826800 H 1.42370000 0.10775400 1.42899000 H 2.47445600 0.12780100 0.00402400 C 0.66866000 1.30146700 -0.16992700 H 1.26438900 2.14276400 -0.50621500 C -0.66890700 1.30134300 -0.16992700 H -1.26479300 2.14252900 -0.50621700 C -1.41434600 0.09433800 0.31826600 H -1.42372700 0.10749400 1.42898900 H -2.47448400 0.12734000 0.00402000 C 0.77040400 -1.21422100 -0.17552700 H 1.14229400 -2.05288100 0.44176200 H 1.12544800 -1.41421500 -1.20543600 C -0.77017000 -1.21437000 -0.17552800 H -1.12517200 -1.41443600 -1.20543700 H -1.14189800 -2.05310200 0.44176200 Title Card Required 0 1
C 0.39335100 -1.40752200 0.50945100 H 0.07415700 -1.04074500 1.48001500 H 0.28923500 -2.47884000 0.39989000 C 1.26702700 -0.69422900 -0.28512700 H 1.85846600 -1.20567200 -1.04409200 C 1.25346800 0.71680900 -0.28478600 H 1.83499300 1.23987800 -1.04346600 C 0.36574900 1.41309100 0.50971100 H 0.05414500 1.04053600 1.48056200 H 0.24248700 2.48246500 0.40066900 C -1.45041500 -0.70335800 -0.25427100 H -1.97290400 -1.26432600 0.51044000 H -1.28115600 -1.25385900 -1.17177300 C -1.46240700 0.67835900 -0.25370700 H -1.30467900 1.23254800 -1.17112400
H -1.99537500 1.22910700 0.51125200
这里需要注意下反应物的构型以及产物的构型都必须是优化过的。QST3最后一个构型是自己猜测的过渡态的构型，所以这里需要一定的化学功底。另外使用QST2以及QST3时要注意你所假定的反应物与产物中间应该没有别的中间体并且只有一个过渡态，否则不保证你能优化出来什么东西~~~所以我比较推荐直接用TS进行过渡态优化，这当然需要涉及到自己猜过渡态的一个过程，需要对一些常见的化学反应机理有个细致的了解才能很好的猜过渡态。选定用哪一种方法后接着就是如何获取力常数，是使用rcfc，calcfc还是calcall由于这一步直接关系着整个体系的计算量所以需要慎重选择，如果你是读上一个计算chk文件进行过渡态优化，并且上一个计算做了频率分析，那毫无疑问要加上rcfc，这可以节省大量的时间。如果从头开始找过渡态的话还是calcfc的好，也就是在第一步计算下力常数，至于calcall个人推荐不到万不得已不要使用，虽然有时它也可以最快的找到过渡态，例如我曾经试过加calcall优化15步找到过渡态，而用calcfc花了47步才找到过渡态，但仍不推荐使用- -！剩下还有一种方法找过渡态——主反应坐标法，该方法通过指定体系中主要参加反应的键的原子通过对其初始Hessian 矩阵中该坐标的行和列计算数值二阶导数的方法来获得过渡态，示例如下。
%chk=c2h4_c4h6_ts_d %mem=200mb %nproc=1
# opt=(ts,modredundant) freq=noraman b3lyp/6-31g* Title Card Required 0 1
C -0.26442100 0.42403600 1.41432400 H -1.29899800 0.82842200 1.42370100 H 0.04131100 0.34416800 2.47445600 C 0.61836500 1.38088300 0.66866500 H 1.22997800 2.04931100 1.26439800 C 0.61836500 1.38088300 -0.66866500 H 1.22997800 2.04931100 -1.26439800 C -0.26442100 0.42403600 -1.41432400 H -1.29899800 0.82842200 -1.42370100 H 0.04131100 0.34416800 -2.47445600 C -0.26442100 -1.61674500 0.77039900 H -1.13801600 -2.18351800 1.14228500 H 0.62859000 -2.16741900 1.12544200 C -0.26442100 -1.61674500 -0.77039900 H 0.62859000 -2.16741900 -1.12544200 H -1.13801600 -2.18351800 -1.14228500 B 8 14 D B 1 11 D
这里有一些我个人用TS找过渡态的经验：优先推荐两种方法，第一步都是将主要参加反应的键长延长0.2-0.5埃或者延长20%-40%，然后固定该键进行构象优化，然后再用TS+calcfc找过渡态或者用主反应坐标法找过渡态，一般情况下这样都能找到相应的过渡态，其中主反应坐标法由于不用算力常数在有些条件下会稍微快一点。这样还没有找到的话可以看一下输出文件中每一步的能量以及构象，然后从中选取最可能是过渡态的那个构象缩小步长
（maxstep=7 or iop(1/8=7)）继续找过渡态。如果仍然没有找到可以对你猜测的构象进行频率分析，然后查看有没有你所要的虚频，如果有的话用GV最下面的manual displacement 调节震动，选中你需要的虚频向左调0.2-0.4或者你不要的虚频向右调0.2-0.4以微调构象，然后接着套用上一种方法的后半部分继续找过渡态。 或者从产物构象SCAN参加反应的键长，逐渐拉长键长然后看能量的变化趋势（注意这里只选择一个键），选取能量的最高点的坐标进行过渡态优化（见图示）。经过以上的方法基本上可以99%以上能找到过渡态。
如果你仍然找到过渡态那可以调整你的初始猜测构象然后重复以上的步骤。重复N次后你仍然找不到过渡态的话请重新思考你所假象的反应通道是否正确，当然如果每次找过渡态都会跑到产物或反应物构象上，则该步反应可能没有过渡态~~所有方法都用了还是找不到过渡态那可以试一下calcall，当然这也不保证你能找到过渡态~~如果你还不放心那你可以用SCAN扫描键长的方法扫描所有可能参加反应的键，然后查看势能面从中选取鞍点做过渡态优化（这会消耗大量的机时，虽然准确但绝不推荐，另本文中的势能面使用am1做的，仅供参考）……
在找过渡态的过程中刚开始的时候不要用大基组来算，用3-21甚至sto-3就可以。如果体系中全部都是前三个周期的原子那用pm6或am1来进行构象初始优化以及过渡态的寻找也是可行的。当用低的基组优化完直接读chk，ts+rcfc进一步用大基组优化过渡态即可。 最后说下过渡态优化经常碰到的两个问题：一个是由于优化步数不够的导致出错，可以在opt里面加上maxcyc=N强制优化N步；另外一个错误是Wrong number of Negative eigenvalues:…… 这个可以通过加上opt=noeigentest参数来解决。
找过渡态是件比较麻烦的事情，但也是件有意思的事情，怎么看就看你自己了。文中虽然以DA反应这么一个简单的反应作为示例，但是实际中找过渡态的过程基本如此。如有纰漏还望大家多多赐教~~