这是本文档旧的修订版!
以F2与CH4的反应作为案例来阐释ADF搜索过渡态的计算过程。反应能垒=活化能=过渡态能量-反应物能量(整个过程的文件下载点击此处)。
ADF搜索过渡态一般分为如下几个步骤:
反应物分子,各自进行大致的优化,参考geoopt。这一步,不需要优化的很精确,因为接下来的计算,都不需要很高精度,只需要提高Transit State那一个步骤的精度,即可以完全弥补前面所有精度不足的问题。
反应物的进攻位点、相对位置对搜索过渡态有决定性影响。不同的进攻位点,有可能对应不同的反应。这一步,对化学直观经验和感觉的要求较高,主要的技术含量也在这里。猜测的好,计算能够少费很多力气,猜测的不好,就会导致无用的计算太多。显然,进攻位点,一般来说,要避开其他无关原子的遮挡。
在本例中,我们猜测F2进攻CH4的姿态如下图:
上图中,按住Shift键,选中F、C原子,右下角出现两个原子之间的距离,之所以这个距离是因为我们估计过渡态键长大约1.4~1.6埃,所以扫描范围略比这个范围大,确保过渡态键长在扫描范围内即可。因此我们将从该距离,扫描到1.2埃。
Linear Transit这个计算的目的,是在限定F-C键长按照我们设定的范围逐渐缩短的情况下,其他原子自由优化使得整体能量更低。参数设置如下。因为精度不需要高,所以这里我们选择计算量小的GGA,使用小基组DZP(较重的元素,可以使用TZP,特别重的元素,例如锕系可以使用TZ2P)、Frozen Core:Large:
按住Shift键,选中F、C两个原子,点击右边窗口的加号,从而可以设置F-C键长初始值为当前值(1.834埃),末值为1.2埃,从1.834到1.2之间扫描5个点应该就足够了。点数越多,势能面越光滑,但计算量也越大,这里基本上大约每隔0.15埃左右就扫了一个点,已经足够光滑了。
扫描的参数可以是键长,也可以是键角、二面角,也可以多个参数配合。但是如果多个参数同时扫描,就有一个“同步性”的问题。因此建议尽量只用一个“决定性”的参数即可,该参数的变化,决定了反应物到产物的变化。
保存任务,并提交作业。
计算过程中,为了节省计算时间,可以在计算过程中,在SCM - Movie中监控计算过程(如果在服务器计算,服务器没有GUI,可以将*.logifle下载到本地用SCM - Movie打开),往往并不需要算完,即可得到需要的结果,例如这里:
其中红色曲线是能量变化曲线,整个扫描过程,有意义的点,只是收敛的点,也就是上图中圆圈所指示的前三个点,其他点是搜索过程的点。可以看到前三个点,随着F-C键长缩短,能量逐渐升高,第四个点还没有收敛,但是看到当前能量已经低于第三个点了,收敛后,能量必定更低。
因此,可以确定第三个点,就是我们想要的——过渡态结构的初始猜测。
因此基于这个结构,进行第3步计算。可以在logfile中,直接复制这个结构对应的坐标,新建Input窗口,然后ctrl v粘贴到该Input窗口,或者在Movie窗口,鼠标拖动到该帧,File - Save Geometry,保存为*.xyz文件。新建Input窗口,File - Import Coordinates导入该*.xyz文件。然后设置第三步的参数,如下。
参数设置:
这一步计算频率的作用有两个:
得到的虚频,可以用SCM LOGO > Spectra打开(在adfjobs窗口,选中任务列表中该相,变成灰色,然后点击SCM LOGO > Spectra),如下图:
下面的表格列出了振动的频率、强度,对于虚频(频率为负),还进行了高精度数值频率矫正,得到更精确的频率、强度。如上图所示黑色方框、红色方框里面,分别是解析频率、数值频率,以及对应的强度。点击该行,即可看到对应的振动模式动画。进而检查是否符合过渡态特征。
如果频率、强度最大的这个虚频确实满足过渡态特征,我们就可以继续基于该分子结构真正开始搜索过渡态,也就是下一步。
如果遇到两个峰的频率、强度接近,如何消除不需要的那个峰呢?参考imaginaryfreq
这一步计算,仍然采用第三步的分子结构,参数设置:
前一步计算生成了一个*.t21文件,这里如果能够在Details > Files (Restart)读取该文件,则能大大提供过渡态搜索的效率、成功率:
计算收敛之后,最大梯度(constrained gradient max)接近为0:
<Dec04-2019> <17:07:07> current energy -0.91629462 Hartree <Dec04-2019> <17:07:07> energy change -0.00000626 0.00100000 T <Dec04-2019> <17:07:07> constrained gradient max 0.00064012 0.00100000 T <Dec04-2019> <17:07:07> constrained gradient rms 0.00030071 0.00066667 T <Dec04-2019> <17:07:07> gradient max 0.00064012 <Dec04-2019> <17:07:07> gradient rms 0.00030071 <Dec04-2019> <17:07:07> cart. step max 0.00501610 0.01000000 T <Dec04-2019> <17:07:07> cart. step rms 0.00210620 0.00666667 T
却是一个非常不稳定的状态,因此可以基本确定这是过渡态。但完全的确认过渡态,还需要再次计算频率,来进行验证。
基于第4步最后收敛的结构,计算频率(注意整个过程最好保持相同的基组、泛函、积分精度、冻芯情况)
如果键级不是很正常,可以在Input - Bonds -Guess bonds重新猜测键级,当然键级是否正确,并不影响计算。第一性原理计算结果并不受键级影响,根本不会读取键级信息。
如第3步一样观察虚频的情况,发现只有一个虚频,并且振动模式比较符合我们的期望(分别朝反应物、产物振动)。
点击该峰(或者该峰对应的横坐标的细的蓝色短线)左边窗口即显示该虚频对应的振动模式,通过View - Open Mode in ADFMovie可以将该振动模式的动画单独窗口显示出来。然后点击该窗口的暂停按钮,将结构从第0帧开始,手动地向右移动3帧,之后暂停在该帧,然后File-Save Geometry保存结构;以及从第0帧开始像向左移动3帧,大致到达录像的倒数第3帧)之后暂停在该帧,然后File-Save Geometry保存结构。
两个结构与“过渡态”相比,分别靠近反应物和产物方向。对这两个结构,分别进行几何结构优化(对应下载的文件中的Reactant和Product),分别得到反应物和产物的结构如下:
确实是我们所预想的反应物和产物,但值得一提的是:反应过程其实和我们设想的并不一致,我们原先以为的是,两个F原子会分别靠近H和C,形成CH3F和HF,但通过计算第6中产物的优化过程,我们发现,反应其实是这样发生的(点击下载录像)。
通过ADF的过渡态搜索,能够搜索到复杂的反应过程。也能够顺利地在搜索过程中,发现可能存在的中间反应过程。ADF搜索过渡态的成功率几乎可以达到100%。
得到反应物、产物、过渡态的能量,也就知道了反应能垒=活化能=过渡态能量-反应物能量。注意在计算能量的时候,需要用相同的参数(基组、泛函、数值精度等)。因此,一般在优化好反应物、产物、过渡态之后,可以统一的比较高精度的参数计算一下三个结构的单点能作为三者的能量。
即使初始的反应路径猜测的不是很准确,甚至是错误的,最后也能通过这个过程,得到正确的反应路径。
ADF软件提供免费试用(一般为一个月),试用申请方式参见费米科技维基百科:AMS免费试用