4.内禀反应坐标(intrinsic reaction coordinate,IRC)
MEP指的是势能面上,由一个点到达另一个点的能量最低的路径,满足最小作用原理。若质量权重坐标下的MEP连接的是反应物、过渡结构和产物,则称为IRC。所谓质权坐标在笛卡儿坐标下即r(i,x)=sqrt(m(i))*R(i,x),m(i)为i原子质量,R(i,x)为i原子原始x方向坐标,同样有r(i,y)、r(i,z)。IRC描述了原子核运动速度为无限小时,质权坐标下由过渡态沿着势能负梯度方向行进的路径(最陡下降路径),其中每一点的负梯度方向就是此处核的运动方向,在垂直于路径方向上是能量极小点。注意质量权重和非权重坐标下的路径是不一样的。
IRC可看作0K时的实际在化学反应中原子核所走的路径,温度较低时IRC也是一个很好的近似。但是当温度较高,即核动能较大时,实际反应路径将明显偏离IRC,而趋于沿最短路径变化,即便经历的是势能面上能量较高的的路径,这时就需要以动力学计算的平均轨迹来表征反应路径。
5.IRC算法
5.1 最陡下降法(Steepest descent)
最简单的获得IRC的方法就是固定步长的最陡下降法,由过渡态位置开始,每步沿着当前梯度方向行进一定距离直到反应物/产物位置,也称Euler法。由于最陡下降法及下文的IMK、GS等方法第一步需要梯度,而过渡态位置梯度为0,所以第一步移动的方向沿着虚频方向。最陡下降方法与IRC的本质相符,但是此法实际得到的路径是一条在真实IRC附近反复震荡的曲折路径,而非应有的平滑路径,对IRC描述不够精确。虽然可以通过更小的步长得以一定程度的解决,但是太花时间,对于复杂的反应机理,需要更多的点。也可以通过RK4(四阶Runga-Kutta)来走步,比上面的方法更稳定、准确,但每步要需要算四个梯度,比较费时。
5.2 IMK方法(Ishida-Morokuma-Kormornicki)
它是最陡下降法的改进,解决其震荡问题。首先计算起始点X(k)的梯度g(k),获得辅助点X'(k+1)=X(k)-g(k)*s,其中s为可调参数。然后计算此点梯度g'(k+1),在g(k)与-g'(k+1)方向的平分线上(红线所示)进行线搜索,所得能量最小点即为X(k+1),之后再将X(k+1)作为上述步骤的X(k)重复进行。整个过程类似先做最陡下降法,然后做校正。此方法仍然需要相对较小的步长,获得较精确IRC所需计算的点数较多。
[图12]IMK方法示意图
Schmidt,Gordon,Dupuis改进了IMK的三个细节,使之更有效率、更稳定。(1)将X'(k+1)的确定方式改为了X(k)-g(k)/|g(k)|*s,即每一步在负梯度方向上行进固定的s距离,与梯度大小不再有关。(2)线搜索步只需在平分线上额外计算一个点的能量即可,这个点和X'(k+1)点的能量以及g'(k+1)在此平分线上的投影三个条件作联立方程即可解出曲线方程,减少了计算量。IMK原始方法则需要在平分线上额外计算两个点的能量与X'(k+1)的能量一起拟和曲线方程。(3)第一步在过渡态位置的移动距离Δq如此确定:ΔE=k*(Δq^2)/2,k为虚频对应的力常数,ΔE为降低能量的期望值(一般为0.0005 hartree),这样可避免在虚频很大的鞍点处第一步位移使能量降低过多。
5.3 Müller-Brown方法
这是通过球形限制性优化找IRC的方法。首先将过渡态和能量极小点位置定义为P1和P2,由P1开始步进,当前步结构以Q(n)表示。每一步,在相距Q(n)为r距离的超球面上用simplex法优化获得能量极小点Q'(图中绿点),优化的起始点是Q(n-1)Q(n)与Q(n)P2方向的平分线b上距Q(n)为r距离的位置S(红点)。若Q(n)Q'与Q(n)P2的夹角较小,则Q'可当作是下一步位置Q(n+1)。如此反复,直到符合停止标准,比如下一步能量比当前更高(已走过头了)、与P2距离已很近(如小于1.2r)、或者与P2方向偏离太大(P1与P2点通过此法无法找到IRC)。最终所得到全部结构点依次相连即为近似的IRC,减小步长r值可使结果更贴近实际IRC。基于此方法也可以用于寻找过渡态,先将反应物和产物作为P1和P2,将二者距离的约2/3作为r,由其中一点在P1-P2连线上相距其r位置为初始位置进行球形优化得到O点,在O与P1、O与P2上也如此获得P1'与P2',根据P1、P1'、O、P2'、P2的能量及之间距离信息以一定规则确定其中哪两个点作为下一步的P1和P2,确定新的P1和P2后重复上述步骤,直至P1与P2十分接近,即是过渡态。此方法计算IRC可以步长可设得稍大,第一步不需要费时的Hessian矩阵确定移动方向,缺点是获得的路径曲率容易有问题,对于曲率较大的反应路径需要减小步长。
[图13]Müller-Brown方法示意图
5.4 GS(Gonzalez-Schlegel)方法
这是目前很常用,也是Gaussian使用的方法,见图14。首先计算起始点X(k)的梯度,沿其负方向行进s/2距离得到X'(k+1)点作为辅助点。在距X'(k+1)点距离为s/2的超球面上做限制性能量最小化,找到下一个点X(k+1)。因为这个点的负梯度(黑色箭头)在弧方向上分量为0,故垂直于弧,即其梯度方向在X'(k+1)到X(k+1)的直线上。这必然可以得到一段用于描述IRC的圆弧(虚线),它通过X(k)与X(K+1)点,且在此二点处圆弧的切线等于它们的梯度方向,这与IRC的特点一致,这段圆弧可以较好地(实线)。之后再将X(k+1)作为上述步骤的X(k)重复进行。
GS方法对IRC描述得比较精确,在研究反应过程等问题中,由于对中间体结构精度有要求,GS是很好的选择,而且用大步长可以得到与小步长相近的结果,优于IMK、Müller-Brown等方法。若只想得到与过渡态相连的反应物和产物结构,或者粗略验证预期的反应路径,对IRC精度要求不高,使用最陡下降法往往效率更高,尽管GS可以用更大步长,但每步更花时间。
[图14]GS方法示意图
除上述外,IRC也可以通过已提及的EF、最缓上升法、球形优化等方法得到,它们的好处是不需要事先知道过渡态的结构。赝坐标法除了简单的反应以外,只能得到近似的IRC,由于结构的较小偏差会带来能量的较大变化,容易引入滞后效应,所以这样得到的势能曲线难以说明问题。
6. chain-of-states方法
这类方法主要好处是只需要提供反应物和产物结构就能得到准确的反应路径和过渡态。首先在二者结构之间以类似LST的方式线性、均匀地插入一批新的结构(使用内坐标更为适宜),一般为5~40个,每个结构就是势能面上的一个点(称为image),并将相邻的点以某种势函数相连,这样它们在势能面上就如同组成了一条链子。对这些点在某些限制条件下优化后,在势能面上的分布描述的就是MEP,能量最高的结构就是近似的过渡态位置。
6.1 Drag method方法
这个方法最简单,并不是严格的chain-of-states方法,因为每个结构点是独立的。插入的结构所代表的点均匀分布在图8所示的短虚线上,也可以在过渡态附近位置增加点的密度。每个点都在垂直于短虚线的超平面上优化,在图中就是指平行于长虚线方向优化。这种方法一般是奏效的,但也很容易失效,图8就是一例,优化后点的分布近似于从产物和反应物用最缓上升法得到的路径(黑色粗曲线),不仅反应路径错误,而且两段不连接,与黑色小点所示的真实MEP相距甚远(黑色点是用下文的NEB方法得到的)。目前基本不使用此方法。
6.2 PEB方法(plain elastic band)
这是下述Chain-of-state方法的基本形式。也是在反应物到产物之间插入一系列结构,共插入P-1个,反应物编号为0,产编号物为P。不同的是优化不是对每个点孤立地优化,而是优化一个函数,每一步所有点一起运动。下文用∑[i=1,P]X(i)符号代表由X(1)开始加和直到X(P)。PEB函数是这样的:S(R(1),R(2)...R(P-1))=∑[i=1,P-1]V(R(i)) + ∑[i=1,P]( k/2*(R(i)-R(i-1))^2 )。其中R(i)代表第i个点的势能面上的坐标,V(R(i))是R(i)点的能量,k代表力常数。优化过程中反应物R(0)和产物R(P)结构保持不变,优化此函数相当于对一个N*(P-2)个原子的整体进行优化,N为体系原子数。
优化过程中,式中的第一项目的是让每个点尽量向着能量极小的位置移动。第二项相当于将相邻点之间用自然长度为0、力常数为k的弹簧势连了起来,目的是保持优化中相邻点之间距离均衡,避免过大。当只有第一项的时候,函数优化后结构点都会跑到作为能量极小点的反应物和产物位置上去而无法描述MEP,这时必然会有一对儿相邻结构点距离很大。当第二项出现后,由于此种情况下弹簧势能很高,在优化中不可能出现,从而避免了这个问题。drag method法在图8中失败的例子中,也有一对儿相邻结构点距离太远,所以也不会在PEB方法中出现。简单来说,PEB方法就是保持相邻结构点的间距尽量小的情况下,优化每个结构点位置。可以近似比喻成在势能面的模型上,将一串以弹簧相连的珠子,一边挂在反应物位置,另一边挂在产物位置,拉直之后松手,这串珠子受重力作用在模型上滚动,停下来后其形状可当作MEP,最高的位置近似为过渡态。
但是PEB方法的结果并不能很好描述MEP。图15描述的是常见的A、B、C三原子反应的LEPS势能面,B可与A或C成键,黑色弧线为NEB方法得到的较真实的MEP。左图中,在过渡态附近PEB的结构点没有贴近MEP,得到的过渡态能量过高,称为corner-cutting问题。这是因为每点间的弹簧势使这串珠子僵硬、不易弯曲,由图15右图可见,R(i)朝R(i-1)与R(i+1)方向都会受到弹簧拉力,其合力牵引R(i),使R(i-1)、R(i)、R(i+1)的弧度有减小趋势。如果将弹簧力常数减小以减弱其效果,就会出现图15中间的情况,虽然结构点贴近了MEP,但相邻点间距没有得到保持,过渡态附近解析度很低,错过了真实过渡态,若以能量最高点作为过渡态则能量偏低,这称为sliding-down问题。可见弹簧力常数k的设定对PEB结果有很大影响,为权衡这两个问题只能取折中的k,但结果仍不准确。