过渡态搜索
1.确定条件
反应路线是根据经验和实验已知现象(如监测到哪些产物生成,投料比等)进行推测的,量化计算的目的主要是通过计算这条路线上的所有体系的吉布斯自由能,说明反应能垒是合理的,可以跨越的。中间体是势能面的极小值点,过渡态是鞍点。然而,实际反应可能有很多条通路通向反应终点,量化计算并不能揭示出这些路线,因此只要计算合理即可,不必纠结。
1.1优化方法
最低级别:MM2(Chem3D自带的),后面会讲。
低级别:pm7
半经验方法,用来快速优化一遍结构,建议所有体系都先用这个优化。关键词只需要加上:
- pm7
高级别:
查阅同一领域的文献,确定常用的方法和级别。选择时也要确定理论方法是适合体系研究的,可以看:
常用泛函:B3LYP, M06-2X, M06等,
基组:6-31系列,def2系列
赝势和赝势基组:SDD, Lan系列
补充:对于比较大的体系,或者含有过渡金属的体系,一种普遍的做法是在相对较低的级别(如B3LYP-D3(BJ)/6-31G(d))下做opt,同时得到热力学校正值,再在高级别(如M06-2X/6-311++G(d,p))下计算单点能等。(两数据如何处理后面会谈到)
1.2溶剂和溶剂模型
溶剂使用和实验一致的。(频率计算涉及的真空等问题后面会谈到)计算过程中务必记得加溶剂,否则可能无法得到优化的结构或过渡态。
现在一般用SMD,精度比较高,但容易出奇怪的错误,一般解决:取报错前一两帧重新计算
关键词:
- scrf=(SMD,solvent=)
也有用PCM(IEFPCM)的,体验比较好,对于过渡态搜索步骤可以先用这个模型,避免奇怪的报错。关键词默认
- scrf=(solvent=)
sob认为溶剂模型对opt的影响不大,可以在opt任务用PCM,在单点能用SMD。也有文献这么做。经验上DFT下使用SMD总会遇到无法收敛的问题,如果是过渡态结构无法收敛会很麻烦,最终还是考虑仅在单点能使用SMD。
2.中间体优化
2.1初猜与关键词
摆出合理初猜结构,获得可信的中间体(需要确保反应是基元反应)。
摆的时候可以先在Chem3D画好,选上面的MM2进行优化,保存为mol结构,用GV打开再保存为gjf(建议去掉connective关键词),分子小或者方便也可以直接GV画。
然后用pm7优化一遍。最后使用DFT。
freq是频率计算,因为这一步计算耗时相对较高,所以可以在结构比较确定了再加上。freq还有一个重要的作用是计算后面需要的吉布斯自由能校正。
关键词:
- opt freq
注意:在所有中间体,包括过渡态结构中,与反应位点无关的基团最好不要有太大变化(比如某单键连接的基团旋转到了平面上方或下方)。对于整个反应历程要有宏观的视角,考虑好后续反应可能的情况,有些可旋转单键从一开始就尽量确定好。另一些单键旋转不涉及反应位点,则通过观察给出能量较低的合理位置。当然未必能做到完善,后续可以进行结构精修,首要的还是在反应位点发生的过渡态结构。
补充:如果中间体带电荷,可以先尝试单独优化,若结果不稳定或能量较大,再向体系中加入带相反电荷的物料。
2.2中间体无法稳定
难以收敛,可以考虑使用更小的优化步长:
- opt=(maxstep=5,notrust)
涉及质子的情况,由于质子的隧穿效应,在猜测中间体结构时,可以自己放质子的位置,不需要去计算质子转移的过渡态。比如下面的情况:
3.过渡态搜索
3.1参考资料
基本概念和操作也可以看:
计算原理:
过渡态虚频物理解释
3.2初猜
三维结构的初猜有几种思路:
A.根据化学直接摆放,调整结构,减少可能会排斥的地方,化学键的形成与断裂摆放合适的角度和距离。
B.查文献,同一领域的计算文献尤其重要,反应位点附近摆的时候按照文献给的三维结构里的键长,键角,二面角等摆,可以大大增加成功率。
C.针对反应坐标,在低计算级别下进行柔性扫描,寻找鞍点作为初猜结构:
对键长(键角,二面角)进行柔性扫描(pm7)(涉及多变量要一个一个扫面,不可一起扫)选择能量高的地方进行过渡态搜索。注意扫描的范围不要太小,可以多给点。
- opt=modredundant
在末尾加上:
- Q i j s step size
Q处选填:X为原子,B为键长,A为键角,D为二面角,后面分别跟2,3,4个原子的labels(GV里view-labels可以查看)
s表示步长模式,后接步数,步长。f表示固定住。
根据实际情况,也可以考虑从远到近的扫描。
注意:
扫描完成后将鞍点及其附近的一两个点都交过渡态,有的时候不一定是在扫出来的鞍点找到。
理想的鞍点应该是一个尖峰。
D.改变反应坐标,如拉长将要断裂的化学键,限制该反应坐标进行优化和频率分析,考察在反应坐标下是否有对应的虚频。有(绝对值100以上的虚频),可以尝试提交。
在GV里看振动频率,虚频代表有过渡态,点start animation看振动方向是否和成键预期一致。另外虚频大小一般会小于-100。
3.3关键词
一般先用ts做,关键词:
- opt=(calcfc,ts,noeigen) freq pm7
calcfc:仅在首次优化中精确计算海森矩阵,有助于过渡态收敛。
noeigen:不检测特征向量的本征值,就是原先计算结构前都需要检测是否只有一个虚频,有多个就报错,现在不检测,也就不会报错了。(因为手动摆的结构不可避免可能会有多个虚频)这个关键词可能会导致收敛到其他结果,但加上不会报错。
提交任务的时候可以微调反应位点的结构,多提交几个,增加成功率。
ts尝试较多失败,或者反应复杂,可以考虑qst2。该算法提供反应物和产物,通过差值法计算若干初猜进行搜索。
打开优化的产物,通过在Open File对话框底部的Target弹出菜单中选择Add all files to active molecule group
tools-connection
Gaussian Calculation Setup对话框指定过渡态结构优化。设置Optimize to a区域为TS(QST2)。
再文本打开gjf文件,修改其他所需关键词。
opt=(calcfc,qst2) freq pm7
QST3不要使用。
3.4找不到过渡态
经验来看除了特别简单的一般提交15个左右才能出一个。
A.观察初猜结构,在成键断键方向上分子的其他部位是否有原子形成空阻。是否有能量更低,更合理的构象。
B.检查低级错误,如溶剂是否正确。
C.势能面特殊,较为平缓,难以收敛(势能面特殊与否通过柔性扫描可以看出一些,更多的只是根据试验结果猜测势能面)
gdiis是一种优化算法,只改变体系在势能面上的运动方式,不会改变势能面。添加关键词:
- opt=gdiis
D.势能面特殊,尖锐而窄,难以搜索到,考虑采用更小的优化步长。
对于优化过渡态,本身默认就是用了notrust的(始终保持不变的优化步长上限),因此关键词为:
- opt=maxstep=5
此外,低级别计算方法由于误差,可能根本没有这个鞍点,因此需要提高计算级别,如使用DFT,甚至DFT结合这个方法。
可以将步长上限(也称置信半径)调整到0.01N Bohr或弧度,默认是30。当按照优化方法判断出的下一步几何变量的变化超过这个数值,则实际的改变量会调节到恰好等于这个上限值(具体来说,是在置信半径对应的球面上做最小化寻找最佳的位置)。如果优化发生震荡,且每步位移不是很大而受力较大的话,往往表明优化过程陷在比较深、比较窄的势阱里,此时可以降低步长上限,比如降低到N=5乃至N=3,以避免每次都走过头而始终达不到势能面极小点。减小步长的做法不要在优化过程中过早地使用,否则由于优化过程走得太慢,所需的优化步数会增加。另外,如果用了calcall,减小步长上限不一定有好处,因为在精确Hessian下进行优化,每一步走的是比较准确的,人为篡改了步长有可能还阻碍收敛。
注意,maxstep设的只是最初的步长上限。随着优化的进行,步长上限其实是动态调整的。如果想一直用自己的设定,需要在opt中加上NoTrust关键词。对于结构反复微小震荡的情况,在使用maxstep减小步长上限的同时总建议同时写上此关键词。
E.有的时候会搜索到其他过渡态,或者结构完全偏离预设,可以考虑特征向量追踪法。IOp(1/130=n)
先是随便把要断掉的键拉长,频率计算, 找到对应TS过程的频率,然后IOp(1/130)选择对应的那个频率 (因为有些不好找的过渡态 可能 初猜的 frequency 的最小mod对应的不是过渡态,所以需要iop一下)。
Eigenvector number to be followed during coordinate driving jobs (IOp(1/19 = 10)).
0 Default (1).
N Follow input structure Hessian eigenvector number N.
但是往往频率最负(也就是需要寻找的TS的)的mode就是1,此方法未必用得上,实际操作中也没有靠这个成功过。
F.参考更多类似文献。
G.该反应就是不能发生,应当重新拟定路线。(不要小瞧了这条,有的时候就是因为这个原因找不到)。判断该路线能否发生,可以将该过渡态的前后中间体作为反应物和产物在Reaxy上进行搜索,看是否有相关实验证据。
H.有的时候是低级别方法找到了过渡态,但是高级别方法计算不出来,检查是否使用了SMD模型,如有,换用PCM溶剂模型。
3.5IRC计算
低级别下对过渡态做IRC计算,可以得到过渡态两端的结构,是确保找到正确过渡态的步骤。但是得到的结构不一定是预期的中间体结构,主要是看成键断键,想获得和自己优化的中间体一致的结构,需要对IRC计算optimization曲线两端的结构再做优化才能得到。
关键词:
- irc=(calcfc)
3.6高级别计算
过渡态搜索:将低级别的过渡态最后一帧导出,用高级别做同样计算。
IRC计算:一般不需要放到文章中去,如果低级别DFT已经做了IRC确认,高级别也可以不做。如果想加速,可以使用刚才的过渡态搜索产生的chk文件做初猜。关键词:
- %chk=
- irc=(rcfc)
4.数据处理
4.1计算吉布斯自由能
文本打开log文件,检索gibbs(或者下面图片中的任何一句)可以得到热力学数据。第一项是零点能,不要和单点能(电子能量)搞混了。单点能在文件最后“/HF=”后面。
如果是使用单一方法,吉布斯自由能=下面最后一项
如果是使用之前提到的相对低级别opt+高级别单点能,
吉布斯自由能=高级别单点能+相对低级别opt的thermal correction to gibbs free energy项
批量提取脚本:
将所有需要提取的log或out文件以及脚本放在一个文件夹内,终端输入:
- python GetCorr.py
补充1:两个不同级别能组合使用的原因是对opt而言,计算级别的影响比较下,对单点能则非常大。所以综合计算开销和精度,可以这么做。
补充2:溶剂中的吉布斯自由能计算比较复杂,上面提到的方式是仅对于解释反应机理而言的。实际计算溶质吉布斯自由能,参考:
4.2计算出的能垒过高
常温一般高于23kcal/mol则能垒偏高,经验上能垒每提高1kcal/mol,反应所需温度需要提高10摄氏度。如果能垒过高检查以下原因:
A.计算级别太低或不适合当前体系,经验上高计算级别会有比较显著的不同。
B.溶剂效应没考虑或考虑不周,诸如有时应当考虑显式溶剂,也就是使用溶剂分子稳定体系,降低能量
C.计算的反应过程和实际反应机理不同 ,也就是说这条路线就是不行,路线猜测是错的
D.过渡态理论的假设不满足,需要考虑隧道效应(质子是比较常见的可以发生隧道效应的物种,因此质子转移的过渡态往往不需要找,即使找到了,由于隧道效应,质子不需要翻越这个能垒就可以转移)
E.实际中有其它物质对反应有催化但是计算时没考虑,本质就是寻找其他物种进一步稳定体系,降低活化能。
4.3绘制自由能图
文章里的自由能图是,纵轴高度表示该化合物与同一路径前一个化合物的吉布斯自由能的差。化合物名称旁边标注的数字表示以int1-1与2的吉布斯自由能之和为基准0,每个化合物的吉布斯自由能。例如,ts1-2在纵轴距离int1-2为8.7个单位,表示ts1-2比int1-2能量高8.7kcal/mol。而下面的8.4表示ts1-2相对int1-1+2能量高8.4kcal/mol。
5 报错
IRC:
Maximum number of corrector steps exceded.
成因:Gaussian的默认IRC算法HPC为了使曲线平滑需做重校正,但常有校正步不收敛的问题。换用GS2算法时,其限制性优化有时也不收敛。
解决:增加IRC=LQA选项,例如IRC=(calcfc,LQA)
6 进阶:全局视角与推理思维
前面讲的都是在已经收到合作组或老师给的初猜机理以后,进行的各种操作。如果在多次尝试了上述所有方法后仍然无法顺利进行,就需要从全局考虑反应路线是否合理等。
6.1开题
必须明确,中间过程中是否分离得到稳定的中间体(是机理猜测的核心点)?副产物是什么?目标产物是什么(化学方程式写好,注意原子守恒)?
6.2路线初猜
中间体能否继续向下反应,可以在reaxys中查询,寻找依据。
过渡态初猜,可以从前一步反应物出发,也可以从后一步产物出发,其他基团摆放符合能量最低原理
整个过程中需要抓住原子之间的静电作用进行思考。