本教程软件版本AMS2020.101,以F2与CH4的反应,演示非周期性体系,也就是分子、团簇体系的过渡态搜索。
ADF搜索过渡态一般分为如下几个步骤:
反应物分子,各自进行大致的优化,参考优化分子的几何结构。这一步,不需要优化的很精确,因为接下来的计算,都不需要很高精度,只需要提高Transit State那一个步骤的精度,即可以完全弥补前面所有精度不足的问题。
反应物的进攻位点、相对位置对搜索过渡态有决定性影响。不同的进攻位点,有可能对应不同的反应。
在本例中,我们猜测F2进攻CH4的姿态如下图:
上图中,按住Shift键,选中F、C原子,右下角出现两个原子之间的距离,之所以这个距离是因为我们估计过渡态键长大约1.4~1.6埃,所以扫描范围略比这个范围大,确保过渡态键长在扫描范围内即可。因此我们将从该距离,扫描到1.3埃。如果缺乏相关经验,则应老老实实地确保扫描起点为反应物,终点为产物。
PES Scan的计算原理,参考势能面扫描PES。参数设置如下:
因为精度不需要高,所以这里我们选择计算量小的GGA,使用小基组DZP(较重的元素,可以使用TZP,特别重的元素,例如锕系可以使用TZ2P)、Frozen Core:Large。
按住Shift键,选中F、C两个原子,点击右边窗口“……(distance)”前面的加号:
从而可以设置F-C键长初始值为当前值(1.851埃),末值为1.3埃,从1.851到1.3之间扫描5个点应该就足够了:
点数越多,势能面越光滑,但计算量也越大,这里基本上大约每隔0.11埃左右就扫了一个点,已经非常光滑了。
扫描的参数可以是键长,也可以是键角、二面角,也可以多个参数配合。建议尽量只用一个“决定性”的参数即可,该参数的变化,决定了反应物到产物的变化。
保存任务,并提交作业。
计算过程中,为了节省计算时间,可以在计算过程中,在SCM - Movie中监控计算过程(如果在服务器计算,服务器没有GUI,可以将*.logifle下载到本地用SCM - Movie打开),往往并不需要算完,即可得到需要的结果,例如这里:
其中红色曲线是能量变化曲线,整个扫描过程,我们只关注收敛的点,也就是上图中圆圈所指示的前4个点。可以看到前4个点,随着F-C键长缩短,能量逐渐升高,第5个点还没有收敛,但是看到当前能量已经低于第4个点了,收敛后,能量必定更低。
因此,可以确定第4个点,就是我们想要的——过渡态结构的初始猜测。
因此基于这个结构,进行第3步计算。可以在logfile中,直接复制这个结构对应的坐标,新建Input窗口,然后ctrl v粘贴到该Input窗口,或者在Movie窗口,鼠标拖动到该帧,File - Save Geometry,保存为*.xyz文件。新建Input窗口,File - Import Coordinates导入该*.xyz文件。然后设置第三步的参数,如下。
参数设置:
这一步计算频率的作用有两个:
得到的虚频,可以用SCM LOGO > Spectra打开(在adfjobs窗口,选中任务列表中该相,变成灰色,然后点击SCM LOGO > Spectra),如下图:
下面的表格列出了振动的频率、强度,对于虚频(频率为负),点击该行,即可看到对应的振动模式动画。进而检查是否符合过渡态特征。如果最大、最强的这个虚频确实满足过渡态特征,我们就可以继续基于该分子结构真正开始搜索过渡态,也就是下一步。
这一步计算,仍然采用第三步的分子结构,参数设置:
前一步计算生成了一个*.t21文件,这里如果能够在Details > Files (Restart)读取该文件,则能大大提供过渡态搜索的效率、成功率:
计算收敛之后,*.logfile显示最大梯度(constrained gradient max)接近为0:
<Nov13-2020> <18:53:35> energy change -0.00000680 0.00007000 T <Nov13-2020> <18:53:35> constrained gradient max 0.00052115 0.00100000 T <Nov13-2020> <18:53:35> constrained gradient rms 0.00019580 0.00066667 T <Nov13-2020> <18:53:35> gradient max 0.00052115 <Nov13-2020> <18:53:35> gradient rms 0.00019580 <Nov13-2020> <18:53:35> cart. step max 0.00222910 0.01000000 T <Nov13-2020> <18:53:35> cart. step rms 0.00067977 0.00666667 T
同时检查结构,确实是一个非常不稳定的状态,因此可以基本确定这是过渡态。但完全的确认过渡态,还需要再次计算频率,来进行验证。
基于第4步最后收敛的结构,计算频率(注意整个过程最好保持相同的基组、泛函、积分精度、冻芯情况)
设置热力学性质计算参数——温度、压强,从而得到该分子结构的Gibbs自由能,用于计算活化能:
计算结束后,也如同如第3步一样观察虚频的情况,发现只有一个虚频,并且振动模式比较符合我们的期望(分别朝反应物、产物振动)。即确认这确实是一个过渡态结构。
SCM - Output - Other properties - Statistical Thermal Analysis,即可看到Gibbs自由能为-566.06kcal/mol:
Temp Transl Rotat Vibrat Total ---- ------ ----- ------ ----- 298.15 Entropy (cal/mol-K): 37.883 22.667 7.632 68.182 Nuclear Internal Energy (kcal/mol): 0.889 0.889 26.892 28.669 Constant Volume Heat Capacity (cal/mol-K): 2.981 2.981 9.094 15.055 (c) Constant Volume Heat Capacity (cal/mol-K): 2.981 2.981 8.872 14.834 Summary of energy terms hartree eV kcal/mol kJ/mol -------------------- ----------- ---------- ----------- Energy from Engine: -0.916314999066729 -24.9342 -575.00 -2405.78 Nuclear Internal Energy: 0.045687613803229 1.2432 28.67 119.95 (c) Nuclear Internal Energy: 0.045553626970109 1.2396 28.59 119.60 Internal Energy U: -0.870627385263500 -23.6910 -546.33 -2285.83 pV/n = RT: 0.000944186013486 0.0257 0.59 2.48 Enthalpy H: -0.869683199250014 -23.6653 -545.73 -2283.35 -T*S: -0.032395578008292 -0.8815 -20.33 -85.05 (c) -T*S: -0.032380708844105 -0.8811 -20.32 -85.02 Gibbs free energy: -0.902078777258306 -24.5468 -566.06 -2368.41
在Spectra窗口点击该唯一虚频,Play - Open Mode in AMSMovie可以将该振动模式的动画单独窗口显示出来。将结构从第0帧开始,手动地向右移动3帧,之后暂停在该帧,然后File-Save Geometry保存结构;以及从第0帧开始像向左移动3帧,大致到达录像的倒数第3帧)之后暂停在该帧,然后File-Save Geometry保存结构。
两个结构与“过渡态”相比,分别靠近反应物和产物方向。对这两个结构,分别进行几何结构优化(对应下载的文件中的Reactant和Product),分别得到反应物和产物的结构。
如果确实是我们所预想的反应物和产物,那么就表示这个过渡态是该反应物、产物之间的过渡态,否则意味着该反应可能不是一步完成的,中间可能还有其他单步反应。
优化结束后,使用各自优化好的分子结构,重新计算在指定温度、压强下的热力学性质,从而得到Gibbs自由能。将上面过渡态得到的Gibbs自由能与反应物的Gibbs自由能相减,即得到活化能。
Gibbs自由能的计算,参考:气相分子的频率、红外IR、零点能、转动能量、转动惯量、熵、焓、热容与Gibbs自由能,液相自由能与焓的计算
简而言之就是过渡态的“能量”-反应物的“能量”。
如果活化能是应用于厄伦尼乌斯公式,用来计算反应速率,那么“能量”一般采用计算结束之后得到的Total Bonding energy(*.logfile尾部的“Bond Energy”、*.out文件中“Total Bonding Energy”),加上零点能(SCM → Output → Other Properties → Zero-Point Energy ,计算频率默认就可以得到该值)。
如果考虑外部条件,则可以用Gibbs自由能的计算时,得到的焓(Enthalpy H)来相减。
如果用Gibbs自由能相减,这个能量差值,是不能用于厄伦尼乌斯公式的。
AMS软件提供免费试用,试用申请方式参见:AMS免费试用