| 查看: 17667 | 回复: 42 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【sobereva个人文集】转载一篇原创精品文章：过渡态、反应路径的计算方法及相关问题已有15人参与
|
|||
|
转载声明：本文由sobereva原创，发表于分子模拟论坛www.mdbbs.org，本人读后深有受益，故而迫切与诸君共享，若转载引来非议，万望原作者sobereva见谅！ 原文地址：http://www.mdbbs.org/thread-17170-1-1.html 请转载者务必注明原作者，否则无疑学术侵权 过渡态、反应路径的计算方法及相关问题 Sobereva Department of Chemistry, University of Science and Technology Beijing, Beijing 100083, China 前言：本文主要介绍过渡态、反应路径的计算方法，并讨论相关问题。由于这类算法极多，可以互相组合，限于精力不可能面面俱到展开，所以只介绍常用，或者实用价值有限但有启发性的方法。文中图片来自相关文献，做了一定修改。由于本文作为帖子发布，文中无法插入复杂公式，故文中尽量将公式转化为文字描述并加以解释，这样必然不如公式形式严谨，而且过于复杂的公式只能略过，但我想这样做的好处是更易把握方法的梗概，有兴趣可以进一步阅读原文了解细节。对于Gaussian中可以实现的方法，文中对其在Gaussian中的使用进行了一些讨论，希望能纠正一些网上流传的误区。虽然绝大多数人不专门研究计算方法，其中很多方法也不会用到，但多了解一下对开阔思路是很有好处的。 文中指的“反应”包括构象变化、异构化、单分子反应等任何涉及到过渡态的变化过程。“反应物”与“产物”泛指这些过程的初态和末态。“优化”若未注明，包括优化至极小点和优化至过渡态。势能面是高维的，但为了直观以及表述方便，文中一般用二维势能面模型来讨论，应推广至高维情况。限于纯文本格式，向量、矩阵无法加粗表示，但容易自行判断。 目录： 1.过渡态 2.过渡态搜索算法 2.1 基于初猜结构的算法 2.1.1 牛顿-拉弗森法(Newton-Raphson,NR)与准牛顿法(quasi-Newton,QN) 2.1.2 AH方法(augmented Hessian) 2.1.2.1 RFO法(Rational Function Optimization，有理函数优化) 2.1.2.2 P-RFO法(Partitioned-RFO) 2.1.2.3 QA法(Quadratic Approximation，二次逼近) 2.1.2.4 TRIM法(trust-region image minimization，置信区域镜像最小化) 2.1.2.5 在高斯中的常见问题 2.1.3 GDIIS法(Geometry Direct Inversion in the Iterative Subspace) 2.1.4 梯度模优化(gradient norm minimization) 2.1.5 Dimer方法 2.2 基于反应物与产物结构的算法 2.2.1 同步转变方法(synchronous transit,ST) 2.2.2 STQN方法(Combined Synchronous Transit and Quasi-Newton Methods) 2.2.3 赝坐标法(pseudo reaction coordinate) 2.2.4 DHS方法(Dewar-Healy-Stewart，亦称Saddle方法)与LTP方法(Line-Then-Plane) 2.2.5 Ridge方法 2.2.6 Step-and-Slide方法 2.2.7 Müller-Brown方法 2.2.8 CI-NEB、ANEBA方法 2.3 基于反应物结构的算法 2.3.1 最缓上升法(least steep ascent,shallowest ascent) 2.3.2 本征向量/本征值跟踪法(eigenvector/eigenvalue following,EF。也称mode walking/mode following/Walking up valleys) 2.3.3 ARTn(activation-relaxation technique nouveau) 2.3.4 梯度极值法(Gradient extremal,GE) 2.3.5 约化梯度跟踪(reduced gradient following,RGF) 2.3.6 等势面搜索法(Isopotenial Searching) 2.3.7 球形优化(Sphere optimization) 2.4 全势能面扫描 3.过渡态相关问题 3.1 无过渡态的反应途径(barrierless reaction pathways) 3.2 Hammond-Leffler假设 3.2 对称性问题 3.3 溶剂效应 3.4 计算过渡态的建议流程 4.内禀反应坐标(intrinsic reaction coordinate,IRC) 5.IRC算法 5.1 最陡下降法(Steepest descent) 5.2 IMK方法(Ishida-Morokuma-Kormornicki) 5.3 Müller-Brown方法 5.4 GS(Gonzalez-Schlegel)方法 6.chain-of-states方法 6.1 Drag method方法 6.2 PEB方法(plain elastic band) 6.3 Elber-Karplus方法 6.4 SPW方法(Self-Penalty Walk) 6.5 LUP方法(Locally Updated planes) 6.6 NEB方法(Nudged Elastic Band) 6.7 DNEB方法(Double Nudged Elastic Band) 6.8 String方法 6.9 Simplified String方法 6.10 寻找过渡态的chain-of-state方法 6.10.1 CI-NEB方法 6.10.2 ANEBA方法(adaptive nudged elastic band approach) 6.11 chain-of-states方法的一些特点 6.12 高斯中opt关键字的path=M方法 6.13 CPK方法(Conjugate Peak Refinement) 1.过渡态 过渡态结构指的是势能面上反应路径上的能量最高点，它通过最小能量路径(minimum energy path,MEP)连接着反应物和产物的结构（如果是多步反应的机理，则这里所指反应物或产物包括中间体）。对于多分子之间的反应，更确切来讲过渡态结构连接的是它们由无穷远接近后因为范德华力和静电力形成的复合物结构，以及反应完毕但尚未无限远离时的复合物结构。确定过渡态有助于了解反应机理，以及通过势垒高度计算反应速率。一般来讲，势垒小于21kcal/mol就可以在室温下发生。 在势能面上，过渡态结构的能量对坐标的一阶导数为0，只有在反应坐标方向上曲率（对坐标二阶导数）为负，而其它方向上皆为正，是能量面上的一阶鞍点。过渡态结构的能量二阶导数矩阵（Hessian矩阵）的本征值仅有一个负值，这个负值也就是过渡态拥有唯一虚频的来源。若将分子振动简化成谐振子模型，这个负值便是频率公式中的力常数，开根号后即得虚数。 分子构象转变、化学反应过程中往往都有过渡态的存在，即这个过程在势能面上的运动往往都会经历满足上述条件的一点。化学反应的过渡态更确切应当成为“反应过渡态”。需要注意的是化学反应未必都经历过渡态结构。 由于过渡态结构存在时间极短，所以很难通过实验方法获得，直到飞秒脉冲激光光谱的出现才使检验反应机理为可能。计算化学方法在目前是预测过渡态的最有力武器，尽管计算上仍有一些困难，比如其附近势能面相对于平衡结构更为平坦得多、低水平方法难以准确描述、难以预测过渡态结构、缺乏绝对可靠的方法（如优化到能量极小点可用的最陡下降法）等。 搜索过渡态的算法一般结合从头算、DFT方法，在半经验、或者小基组条件下，难以像描述平衡结构一样正确描述过渡态结构，使得计算尺度受到了限制。结合分子力场可以描述构象变化的过渡态，但不适用描述反应过渡态，因为大部分分子力场的势函数不允许分子拓扑结构的改变，虽然也有一些力场如ReaxFF可以支持，有的力场还有对应的过渡态原子类型，但目前来看适用面仍然较窄，而且不够精确，尽管更为快速。 注：严格来说，“过渡结构”是指势能面上反应路径上的能量最高点，而“过渡态”是指自由能面上反应路径上的能量最高点，由于自由能变主要贡献自势能部分，所以多数情况二者结构近似一致。但随着温度升高，往往熵变的贡献导致自由能面与势能面形状发生明显偏离，从而导致过渡结构与过渡态明显偏离，两个词就不能混用了。但本文不涉及相关问题，故文中过渡态、过渡结构一律指势能面上反应路径上的能量最高点。 2.过渡态搜索算法 2.1 基于初猜结构的算法 2.1.1 牛顿-拉弗森法(Newton-Raphson,NR)与准牛顿法(quasi-Newton,QN) NR法是寻找函数一阶导数为0（驻点）位置的方法。通过对能量函数的泰勒级数的二阶近似展开，然后使用稳态条件dE/dr=0，可导出步进公式：下一步的坐标向量 = 当前坐标向量 - 能量一阶导数向量 * Hessian矩阵的逆矩阵。在势能面上以NR法最终找到的结果是与初猜位置Hessian矩阵本征值正负号一致、离初猜结构最近的驻点，由于能量极小点、过渡态和高阶鞍点的能量一阶导数皆为0，故都可以用NR法寻找。 对于纯二次形函数NR法仅需一步即可找到正确位置，而势能面远比之复杂，所以需要反复走步直至收敛。也因为势能面这个特点，为了改进优化，实际应用中NR法一般还结合线搜索步(line search)，对于优化至极小点，就是找当前点与NR法算出来的下一点的连线上的能量极小点作为实际下一步结构；若优化至过渡态，且连线方向主要指向过渡态，则找的是连线上能量极大点，若主要指向其它方向则找连线的能量极小点，若指向二者程度均等则一般不做线搜索。由于精确的线搜索很花时间，所以一般只是在连线的当前位置附近计算几个点的能量，以高阶多项式拟和后取其最小/最大点。 NR法每一步需要计算Hessian矩阵并且求其逆，所以十分昂贵。QN法与NR法的走步原理一样，但Hessian矩阵最初是用低级或经验方法猜出来的，每一步优化中通过当前及前一步的梯度和坐标对Hessian矩阵逆矩阵逐渐修正。由于只需计算一阶导数，即便Hessian矩阵不准确造成所需收敛步数增加，但一般仍比NR法速度快得多。QN法泛指基于此原理的一类方法，常用的是BFGS(Broyden Fletcher Goldfarb Shanno)，此法对Hessian的修正保持其对称性和正定性，最适合几何优化，但显然不能用于找过渡态。还有DFP(Davidon-Fletcher-Powell)，MS(Murtagh-Sargent，亦称symmetric rank 1，SR1)，PSB(Powell-symmetric-Broyden)。也有混合方法，如Bofill法是PSB和MS法对Hessian修正量的权重线性组合，比二者独立使用更优，权重系数通过位移、梯度改变量和当前Hessian计算得到，它对Hessian的修正不强制正定，很适宜搜索过渡态。 将NR步进公式放到Hessian本征向量空间下其意义更为明显（此时Hessian为对角矩阵），可看出在每个方向上的位移就是这个方向势能的负梯度除以对应的本征值，比如在i方向上的位移可写为ΔX(i)=-g(i)/h(i)，在受力越大、越平坦的方向位移越大。每一步实际位移就是这些方向上位移的矢量和。对于寻找过渡态，因为虚频方向对应Hessian本征值为负，使位移为受力相反方向，所以NR法在过渡态附近每一步都是使虚频方向能量升高，而在其它正交的方向朝着能量降低的方向位移，通过这个原理步进到过渡态。若有n个虚频，则NR法就在n个方向升高能量而其它方向降低能量找到n阶鞍点。 由于NR法的这个特点，为找到正确类型的驻点，初猜结构必须在目标结构的二次区域(quadratic region)内。所谓的二次区域，是指驻点附近保持Hessian矩阵本征值符号不变的区域，它的形状可以用多变量的二次函数近似描述，例如二维势能面情况下这样的区域可以用F(x,y)=A*x^2+B*y^2+C*x+D*y+E*x*y来近似描述。对于能量极小点，就是指初猜点在目标结构附近Hessian矩阵为正定矩阵的范围；对于找过渡态，就需要初猜点在它附近含有且仅含有一个负本征值的范围内。并且这个范围内不能有其它同类驻点比目标结构距离初猜结构更近。NR法方便之处是只需要提供一个初猜结构即可，但是由于过渡态二次区域很小（相对于能量极小点来讲），复杂反应过渡态又不容易估计，故对使用者的直觉和经验有一定要求，即便是老手，也往往需要反复尝试。 NR法对初猜结构比较敏感，离过渡态越近所需收敛步数越少，成功机率越高。模版法可以帮助给出合理的初猜，也就是如果已经知道其它机理相同的反应的过渡态结构，可以保持反应位点部分的结构不变而替换周围的原子，使之变成自己要研究的化合物反应的初猜结构。 2.1.2 AH方法(Augmented Hessian) AH方法并不是独立的寻找过渡态的方法，而是通过修改原始Hessian矩阵来调整NR法步进的长度和方向的一种方法。在NR法的步进公式中加入了一个移位参数λ，式子变为ΔX(i,λ)=-g(i)/(h(i)-λ)，NR法相当于λ等于0时的特例。λ控制着每步步进距离，它与h(i)的相对大小也控制着这个坐标上的步进方向。根据设定λ方法的不同，常见的有RFO、P-RFO和QA/TRIM。这些方法每一步也使用QN方法来快速地更新Hessian。 下面提及的置信半径R(Trust radius)是指二阶泰勒级数展开这种近似的合理的区域，可以在优化过程中固定也可以动态改变，比如下一步位置的实际能量与使用二阶泰勒级数展开预测的能量符合较好则加大R，反之减小。优化的每一步移动距离不应超过R，否则可能进入二阶泰勒级数展开近似的失效区域，NR法在势能面平坦的时候容易超过这个范围，应调整λ避免。 2.1.2.1 RFO法(Rational Function Optimization，有理函数优化) 对能量函数根据有理近似展开，而不是NR法的二阶泰勒级数近似展开，可推得与AH方法形式相同的步进公式。确定其中λ的公式是λ=∑( g(i)^2/(h(i)-λ) )，g(i)和h(i)代表此方向的梯度和本征值，加和是对所有本征向量方向加和。通过迭代方法会解出N+1个λ（N代表势能面维数），将λ按大小排列，则有λ(i)≤h(i)≤λ(i+1)。故选其中最小的λ可使各个方向位移公式的(h(i)-λ)项皆为正，保证每步位移都向着极小点。选其中大于m个Hessian本征值的λ，将会在本征值最低的m个方向上沿其上的受力反方向位移提升能量，在其余N-m个方向上降低能量，由此确保优化到m阶鞍点，若m为1即用来找过渡态。所以用了这个方法寻找指定类型驻点不再像NR法对初猜位置Hessian本征值符号有要求，而是直接通过选择λ来设定向着何种鞍点位移。如果每步步长度超过了R，则乘以一个小于1的因子来减小步长。值得一提的是，λ与势能面维数N的平方根近似成正比，随着体系尺度的增大，RFO的λ对NR法的二次近似就会趋现“校正过度”情况，产生大小不一致问题，可使用SIRFO(Size independent RFO)方法解决，即AH走步公式中的λ改为λ/N^0.5。 [ Last edited by yjcmwgk on 2010-6-22 at 20:48 ] |
回复此楼» 收录本帖的淘帖专辑推荐
» 本帖已获得的红花（最新10朵）
hookah» 猜你喜欢
厦门大学国家杰青侯旭教授课题组招聘人工智能、机器学习相关交叉方向博士研究生！
已经有3人回复
-
LSV曲线
已经有0人回复
-
物理化学论文润色/翻译怎么收费?
已经有131人回复
-
美国西北大学 生物医学实验室 招聘优秀博士后！
已经有3人回复
-
北京化工大学普通博士招生---生物方向
已经有2人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- vasp计算过渡态问题
已经有13人回复
- 关于过渡态的计算的能量问题。
已经有7人回复
- 【求助】自由基耦合反应过渡态有过渡态吗？怎样计算其过渡态？
已经有4人回复
- NEB计算过渡态出现问题（急！！！！）
已经有9人回复
- CH4在Ir（111）表面的解离吸附过渡态计算：Henkelman的PRL文章
已经有10人回复
- 有机非基元反应的过渡态计算如何确定中间体的个数和结构？
已经有4人回复
- 求助如何对应输入反应物和产物进行过渡态计算
已经有13人回复
- vasp 过渡态计算 线性插值方法。
已经有8人回复
- 【整理】今天整理硬盘上的资料，偶然发现的好玩的小东西就直接发上来吧
已经有105人回复
- 【求助】过渡态计算的QST3方法
已经有12人回复
- 【第一性原理活动贴】VASP 官方手册学习活动
已经有56人回复
- 【分享】【资源】计算模拟区区刊第一期正式发布
已经有6人回复
- 计算模拟区区刊RCS第一期正式发布
已经有119人回复
jackwangee
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 1168
- 散金: 3
- 红花: 2
- 帖子: 124
- 在线: 185.6小时
- 虫号: 517335
- 注册: 2008-03-03
- 性别: GG
- 专业: 理论和计算化学
41楼2014-11-18 10:45:51
2