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$。其中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,但结果仍不准确。 img [图15]LEPS势能面上不同k值的PEB结果

6.3 Elber-Karplus方法

与PEB函数定义相似。第一项定义为$1/Li=1,P-1$,其中L为链子由0点到P-1点的总长,$d(i,i+1)$为R(i)与R(i+1)的距离,此项可视为所有插入点总能量除以点数,即插入点的平均能量。第二项为$γi=1,P^2$,其中代表相邻点的平均距离,是所有d(i,j)的RMS。此项相当于将弹簧自然长度设为了当前各个弹簧长度的平均值,由γ参数控制d(i,j)在平均值上下允许的波动的范围。此方法最初被用于研究蛋白质体系的构象变化。

6.4 SPW方法(Self-Penalty Walk)

在Elber-Karplus方法的基础上增加了第三项互斥项,$∑[i=0,P-1]∑[i=j+1,P-1]U(ij)$,其中$U(ij)=ρ*exp(-d(i,j)/(λ))$,定义同上。此项相当于全部点之间的“非键作用能U(ij)”之和,不再仅仅是相邻点之间才有限制势。任何点之间靠近都会造成能量升高,可以避免Elber-Karplus方法中出现的在能量极小点处结构点聚集、路径自身交错的问题,能够使路径充分地展开,确保过渡态区域有充足的采样点。式中ρ和λ都是可调参数来设定权重。此外相对与Elber-Karplus方法还考虑了笛卡儿坐标下投影掉整体运动的问题。

6.5 LUP方法(Locally Updated planes)

特点是优化过程中,只允许每个结构点R(i)在垂直于R(i-1)R(i+1)向量的超平面上运动。由于每步优化后R(i-1)与R(i+1)连线方向也会变化,故每隔一定步数重新计算这些向量,重新确定每个点允许移动的超平面。但是LUP缺点是结构点之间没有以上述弹簧势函数相连来保持间隔,容易造成结构点在路径上分布不均匀,甚至不连续,还可能逐渐收敛至两端的极小点。

6.6 NEB方法(Nudged Elastic Band)

NEB方法集合了LUP与PEB方法的优点,其函数形式基于PEB。从PEB方法的讨论可以看出,弹簧势是必须的,它平行于路径切线(R(i)-R(i-1)与R(i+1)-R(i)矢量和的方向)的分量保证结构点均匀分布在MEP上来描述它;但其垂直于路径的分量造成的弊端也很明显,它改变了这个方向的实际的势能面,优化后得到的MEP’就与真实的MEP发生了偏差,造成corner-cutting问题。解决这个问题很简单,在NEB中称为nudge过程,即每个点在平行于路径切线上的受力只等于弹簧力在这个方向分量,每个点在垂直于路径切线方向的受力只等于势能力在此方向上分量。这样弹簧力垂直于路径的分量就被投影掉了,而有用的平行于路径的分量完全保留;势能力在路径方向上的分量也不会再对结构点分布的均匀性产生影响,被保留的它在垂直于路径上的分量将会引导结构点地正确移动。这样优化收敛后结构点就能正确描述真实的MEP,矛盾得到解决。弹簧力常数的设定也比较随意,不会再对结果产生明显影响。但是当平行于路径方向能量变化较快,垂直方向回复力较小的情况,NEB得到的路径容易出现曲折,收敛也较慢,解决这一问题可以引入开关函数,即某点与两个相邻点之间形成的夹角越小,此点就引入更多的弹簧势垂直于路径的分量,使路径不易弯曲而变得光滑,但也会带来一定corner-cutting问题。也可以通过将路径切线定义为每个点指向能量更高的相邻点的方向来解决。

6.7 DNEB方法(Double Nudged Elastic Band)

弹簧势垂直于路径的分量坏处是造成corner-cutting问题,好处是避免路径卷曲。更具体来说,前者是由于它平行于势能梯度方向的那个分量造成的,若只将这个分量投影掉,就可避免corner-cutting问题,而其余分量的力F(DNEB)仍可以避免路径卷曲,这便是DNEB的主要思想。故DNEB与NEB的不同点就是DNEB保留了弹簧势垂直于路径的分量其中的垂直于势能梯度的分量。

DNEB的这个设定却导致结构点不能精确收敛到MEP上。正确的MEP上的点在垂直于路径方向上受势能力一定为0,但是当用了DNEB方法后,若其中某一点处路径是弯曲的,即弹簧力在垂直于路径方向上有分量F’,而且此点势能梯度方向不垂直于此点处路径的切线,即F’不会被完全投影掉,F’力的分量F(DNEB)将继续带着这个点移动,也就是说结构点就不在正确的MEP上了。只有当结构点所处路径恰为直线,即F’为0则不会有此问题。为了解决此问题有人将开关函数加入到DNEB,称为swDNEB,当结果越接近收敛,即垂直于路径的势能力越小的时候,F(DNEB)也越小,以免它使结构点偏离正确MEP。一些研究表明DNEB和swDNEB相比NEB在收敛性(结构点受力最大值随步数降低速度)方面并没有明显提升,DNEB难以收敛到较高精度以内,容易一直震荡。

6.8 String方法

与NEB对力的投影定义一致,但点之间没有弹簧势连接,保持点的间距的方法是每步优化后使这些点在路径上平均分布。

6.9 Simplified String方法

String中计算每个点的切线并投影掉势能力平行于路径的分量的过程也去掉了,所有点之间用三次样条插值来表述路径,每一个点根据实际势能力运动后,在路径上重新均匀分布。优化方法最好结合RK4方法。NEB在点数较小的情况下比Simplified String方法能在更短时间内收敛到更高精度,但点数较多情况下则Simplified String更占优势。

6.10 寻找过渡态的chain-of-state方法 除非势能面对称且结构点数目为奇数,否则不会有结构点恰好落在过渡态。以能量最高的点作为过渡态只是近似的,为了更好地描述过渡态,可以增加结构点数,或者增加局部弹簧力常数,使过渡态附近点更密。根据已得到的点的能量,通过插值方法估算能量最高点是另一个办法。近似的过渡态也可以作为QN法的初猜寻找准确的过渡态。

6.10.1 CI-NEB方法

NEB与String等方法都可以结合Climbing Image方法,它专门考虑到了定位过渡态问题。CI-NEB与NEB的关键区别是能量最高的点受力的定义,在CI-NEB中这个点不会受到相邻点的弹簧力,避免位置被拉离过渡态,而且将此点平行于路径方向的势能力分量的符号反转,促使此点沿着路径往能量升高的方向上爬到过渡态。这个方法只需要很少的点,比如包含初、末态总共5个甚至3个点就能准确定位过渡态,是最有效率的寻找过渡态的方法之一。如果还需要精确描述MEP,可以在此过渡态上使用Stepwise descent方法、最陡下降法、RK4等方法沿势能面下坡走出MEP,整个过程比直接使用很多点的NEB方法能在更短时间内得到更准确的MEP。

6.10.2 ANEBA方法(adaptive nudged elastic band approach)

这个方法也是基于NEB,专用来快速寻找过渡态。一般想得到高精度的过渡态区域,NEB的链子上必须包含很多点,耗费计算时间。而ANEBA方法中链子两端的位置不是固定的,而是不断地将它们移动到离过渡态更近的位置,仅用很少几个点的链子就可以达到同样的精度。具体来说,设链子两端的点分别叫A点和B点(对于第一步就是反应物和产物位置),先照常做NEB,收敛至一定精度后(不需要精度太高),改变A和B的位置为链子中能量最高点相邻的两个点,然后再优化并收敛至一定精度,再如此改变A和B的位置,反复经历这一步骤,最终链子上能量最高点就是精确的过渡态。ANEBA相当于不断增加原先NEB链子的过渡态附近的点数,但实际上点数没有变。有研究表明ANEBA比CI-NEB效率更高,如果结合ANEBA与CI(称CI-ANEBA),即先用ANEBA方法经上述步骤移动几次A、B点,使之聚焦到过渡态附近,再用CI-NEB方法,效率可以进一步提高。 img [图16]ANEBA方法示意图

6.11 chain-of-states方法的一些特点

NEB方法的设定只是决定了每一步结构点实际感受到的势能面是怎么构成的,并没有指定优化方法。NEB可以结合一些常见的优化方法,比如最陡下降法、共轭梯度法、quick-min、FIRE、L-BFGS法等(没有线搜索步的全局L-BFGS法效率一般最高),但只能像前述寻找IRC方法一样得到一条路径。实际上很多情况反应的路径不止一条,尤其是势能面复杂的大分子构象转变过程。当NEB结合构象搜索方法,比如分子动力学、蒙特卡罗等方法时,就可以用于寻找多条反应路径。例如有几条反应路径,彼此间都有一定高度的势垒分隔,如果初始给出的路径在第i条附近,优化后只能收敛到第i条路径,若对每个点使用分子动力学方法,设定一定温度,则这些点有机会凭借动能越过势垒到达另外一条路径k附近,随后逐渐降温减小动能,相当于对它们进行最陡下降法优化,就找到了第k条路径,若如此反复多次,有可能找到更多路径。

这类chain-of-states方法的优点还在于易于实现,算法简单,只有能量和其一阶导数是必须要算的,随着体系尺度增大计算量的增加远比基于Hessian矩阵的方法要小。对于大体系储存Hessian并求逆亦是困难的,在某些情况下Hessian矩阵受计算能力制约只能在低水平方法下得到或者无法获得,chain-of-states方法避免了这个问题,很适合用于分子力学研究生物大分子的结构变化路径以及平面波基组下的DFT方法研究固体表面化学反应。此方法也容易并行化,例如可以每个节点负责优化其中一个或几个点,只有计算弹簧力时才需要从另外节点传入相邻结构点坐标,数据通信量小,并行效率高。

6.12 高斯中opt关键字的path=M方法

与chain-of-states方法有一定类似之处,可以在一次计算中获得优化后的过渡态、产物、反应物以及用于描述IRC的中间点结构,总共M个点。此方法须结合QST2或QST3关键字。结合QST2时,除反应物和产物以外剩下的M-2个点在二者冗余内坐标下线性插值产生,结合QST3则是剩下的M-3个点在反应物与过渡态、过渡态与产物之间插值产生。之后迭代的每一步主要分为以下几个步骤:(1)初始输入的反应物、产物通过RFO法向最优构型优化。(2)能量最高的点q(k)(此点在第一步确定)通过EF法向过渡态优化,并设一段圆弧通过q(k-1)、q(k)、q(k+1),此圆弧在q(i)处的切线作为EF方法选择所跟踪的本征向量的引导,类似于STQN步。(3)其余的点执行微迭代步骤(迭代内的迭代),其中包含类似于GS法的球面优化步骤以及调整间距步骤。可参考图14,优化其中任意点q(i)前,首先获得经过q(i-1)、q(i)并与q(i-1)的梯度相切的圆弧或曲线,将其在q(i)处的切线定义为T(i),然后定义一个在q(i)处法线与T(i)平行、经过q(i-1)与q(i)的球面,使q(i)限制在此球面上优化。然后在这些点依次相连的路径上调整这些点的间距至平均,之后重复微迭代直至每一步力和位移都已收敛,或者有任何点位移超过了置信半径。(4)检查力和位移是否都已收敛至标准。这个方法看上去比单独优化反应物、产物、过渡态并计算IRC省时间,但是实际用起来并不好用,容易出错。

6.13 CPK方法(Conjugate Peak Refinement)

在某种意义上称为动态的chain-of-states方法。每条链子只含一个可动点,链子数由最初的一条开始不断增加,对MEP的描述也越来越精确。CPK中的第一步类似LST,在连接反应物和产物的直线中找到能量最高点(称为Peak),然后沿着共轭方向优化得到中间点,对中间点与反应物、中间点与产物分别再做上述步骤,先找到最大点再共轭优化,如此反复直到收敛。最后将反应物、产物以及执行CPK过程中所有优化后的点相连,就得到了近似的反应路径。CPK方法所得的反应路径可以经过很多过渡态,很适合寻找一些涉及到复杂结构重排、包含甚至上百个过渡态的构象变化路径,如蛋白质局部折叠/去折叠过程。CPK方法缺点是实现起来相对复杂,定位过渡态较为费时。 img [图17]CPK方法示意图