计算化学公社
标题: LookForMECP: 一个验证极小势能面交叉点（MECP）充分条件的工具 [打印本页]
作者Author: bnulk 时间: 2019-5-26 08:42
标题: LookForMECP: 一个验证极小势能面交叉点（MECP）充分条件的工具
本帖最后由 bnulk 于 2019-5-26 10:38 编辑
1 前言
对于稳定结构和过渡态，用谐振分析做最优性检验，是计算可靠性的基本保证。对于极小势能面交叉点，也需要做最优性检验，才能保证优化得到的结果是真正的极小点。这里，提供了一个既约谐振分析的程序，来做这项工作。
预印本网址：https://doi.org/10.26434/chemrxiv.8150657.v1
2 原理
主要依据为： (a) Nocedal, J.; Wright, S. J. Numerical optimization, Springer: New York, 1999。 也参考了 (b)Bertsekas, D. P. Nonlinear programming,2nd, Athena Scientific, Belmont,1999.。其中第一本书有影印本，第二本书有中文的译本。
2.1 问题描述
极小势能面交叉点（MECP），数学上是一个约束极值点。即：在“两个电子态能量相等”这个约束条件下，能量的极小值。
(, 下载次数 Times of downloads: 90)
式中E1(q) 和 E2(q) 分别是两个态核坐标q的能量函数，约束条件写为c(q)= E1(q) – E2(q)。这里需要引入一个基本假定，即：函数E1(q) 和 E2(q) 、它们的一阶导数和二阶导数都是存在的和连续的。
通过构建拉格朗日函数，把约束极小值化为无约束极小值。
(, 下载次数 Times of downloads: 100)
等式左边为拉格朗日函数，λ 是拉格朗日因子。这样，MECP最优性检验和稳定结构的最优性检验，就很类似。
2.2 必要条件
在满足线性独立约束条件（linear independence constrainted qualification，LICQ）的前提下，一阶必要条件描述如下：
(, 下载次数 Times of downloads: 106)
式中 q* 和 λ* 分别是MECP候选解的坐标和拉格朗日因子。
对于MECP来说，第一个等式即拉格朗日函数对坐标的一阶导数为零；第二个等式表示两个态的能量相等。一阶必要条件称为Karush-Kuhn-Tucker（KKT）条件；满足一阶必要条件的候选解，称为KKT点。
2.3 充分条件
在满足KKT条件的前提下，二阶充分条件描述如下：
(, 下载次数 Times of downloads: 96)
式中d 是切面M = {d : ∇c(q*) d = 0} 中的任意矢量，并且 d≠ 0。 L = ∇qql(q*,λ*) 是拉格朗日函数在KKT点处的二阶导数。
对于MECP来说，充分条件是拉格朗日函数对坐标的二阶导数，在其切面上的投影正定。
2.4 关于充分条件
对于必要条件，算法上实现非常简单。而对于充分条件，则有些复杂。情况如下图所示。
(, 下载次数 Times of downloads: 101)
图中给出了低维下充分条件的几何图像。（a）和（d）分别是极小点和极大点，类似于稳定结构的充分条件。（b）和（c）示意了MECP充分条件验证中的关键部分。在（c）图中，约束梯度和切面上任意矢量d 正交，q*是一个鞍点。否则，在（b）图中是一个真正的MECP。
我们的想法是，把约束梯度作为一个矢量，从整个分子坐标空间中投影出去，在剩下的互补空间中，拉格朗日函数正定。
3 使用方法
程序暂时只能和高斯（Gaussian 09）联用，需要高斯程序计算梯度和Hessian矩阵。
3.1 输入文件
输入文件也和高斯程序的输入文件类似。例如：
%chk=State1.chk State2.chk
%mem=60gb
%nprocshared=16
#p b3lyp/aug-cc-pvtz scf=maxcycle=300 {task=mecp method=ln hessianN=5 stepSize=0.1 Lambda=1.17 mecpFreq=liu}
MECP
0 1 3
H
N 1 B1
C 2 B2 1 A1
O 3 B3 2 A2 1 D1 0
B1 = 1.03289
B2 = 1.609249
B3 = 1.13
A1 = 103.080264
A2 = 113.938387
D1 = 179.998702
有三处改变：
（1） %Section
%chk关键词的参数，分别是两个态相应的chk文件名。
（2）Route Section
加了大括号部分。大括号内的关键词和参数，后面详述。
（3）Charge& Multipl.:
分别是体系的电荷，第一个态的自旋多重度，第二个态的自旋多重度。
3.2 关键词列表
大括号内的关键词和参数，用“=”号隔开，不区分大小写。
## task: 计算任务。
= mecp: 优化MECP。
## method: 优化方法.
= ln: Lagrange-Newton方法，对于优化MECP来说，是默认方法。
## cyc: 优化最大步数。
## stepSize: 优化每步的步长。
## hessianN: 每N步计算一次Hessian矩阵。如果没有计算Hessian矩阵，则会用猜测的方法获取Hessian矩阵。
= N: 设定N值.
## guessHessian: 猜测Hessian矩阵的方法。
= bfgs: Broyden–Fletcher–Goldfarb–Shanno算法。默认方法。
= powell: Powell算法。
## energyCon: 指定能量差的收敛限。默认值是 0.00001 hartree。
## maxCon: 指定最大力收敛限。默认值是 0.001a.u.。
## rmsCon: 指定均方根力收敛限。默认值是0.0005 a.u.。
## Lambda: Lagrange因子的初始值。理论上可以是任意值，相应的设置技巧，需要专题叙述。也可参考文章。
## mecpFreq: MECP的频率分析.
= liu: 约束频率分析。 默认方法。
= simple: 不使用约束频率分析。简单的检查充分条件。
3.3 输出文件
必要条件：
(, 下载次数 Times of downloads: 104)
“Delta Energy”表示两个态的能量差； “Maximum KKT Force”和“RMS KKT Force” 分别表示最大KKT力和均方根KKT力。
KKT点信息：
bnulk@foxmail.com-MECPResult
*********************************************
Energy= -168.6070306
Lambda= 1.4178592621661
-Lambda/(1-Lambda)= 3.39315025546209
Gradient ratio between two states:
B1 = smallGradient
B2 = 3.4
B3 = 3.47
A1 = smallGradient
A2 = 3.39
D1 = smallGradient
“Energy” 指在KKT点处，两个态的平均能量； “Lambda” 是KKT处的拉格朗日因子值；“Gradient ratio between two states”是两个态键参数的梯度比值。
“在MECP处，两个态的键参数梯度成比例“是MECP必要条件的另一种表述。详见： “Koga, N.; Morokuma, K.Chem. Phys. Lett. 1985, 119, 371.” 和“Liu, K.; Li, Y.; Su,J.; Wang, B. J. Comput.Chem .2014, 35, 703.”。
充分条件：
(, 下载次数 Times of downloads: 102)
既约振动分析。有虚频表示该KKT点并不是真正的MECP，而是一个鞍点。
3.4 最终结果
(, 下载次数 Times of downloads: 95)
左边图是在B3LYP/Aug-cc-PVTZ水平下，优化得到的两个异氰酸单三态的MECP候选点，图中给出了键长和能量（单位为埃和哈垂）。
右边图中给出了既约谐振分析（RVA）。右上图有负频率，说明KKT-1不是真正的极小点（不是MECP）；右下图没有负频率，KKT-2是真正的MECP。
既约振动分析（RVA）共有3N-7个振动，N是分子中原子个数。3N个运动中，6个平动和转动被分离出去（如果是线性MECP，则投掉5个，RVA是3N-6个振动），另一个被投影掉的是约束条件。在网址给出的文章中，提供了详细的计算细节。
需要注意的是，约束振动的模式，并不表示真正的（势能相关的）振动。但是，沿着负频率方向改变分子结构，对获取真正的MECP是有帮助的。
作者Author: pyscf 时间: 2019-5-26 14:14
你这个方法还需要hessian 实用性大打折扣
很多高端的方法撑死只有梯度
作者Author: bnulk 时间: 2019-5-26 21:52
理论上必须有Hessian的。实际上我也没找到好的技巧避开。
作者Author: pyscf 时间: 2019-5-26 22:40
MECP涉及激发态的计算 普通的dft理论做激发态研究局限性很大
你要用好点的方法 多参考之类的 它们最多只支持到解析梯度
作者Author: k64_cc 时间: 2019-5-27 11:27
都上多参考了谁还用dE搜MECP……
作者Author: Kite1990 时间: 2019-6-10 16:08
按照你的方法做了一个尝试，提交任务以后出了问题，能帮忙看一下吗？
作者Author: bnulk 时间: 2019-6-10 18:41
把输入文件发过来。或者发我邮箱：bnulk@foxmail.com
作者Author: 一颗赛艇 时间: 2019-10-19 06:19
请教楼主一个问题，
对于MECP来说，充分条件是拉格朗日函数对坐标的二阶导数，在其切面上的投影正定。
是不是指lambda函数的hessian所有的特征值全是正的？谢谢
作者Author: bnulk 时间: 2019-10-21 16:30
不是这个意思。
分三个角度解释下：
（1） 您说的情况（拉格朗日函数Hessian矩阵特征值全部是正的），确实是MECP的充分条件，但是要求太苛刻了（太充分了）。如果强行使用，就有把真正的MECP误判为非MECP的风险。
（2） 如果拉格朗日函数Hessian矩阵特征值只有一个负值（不是全部正值），且这个负值方向被约束条件约束了，那么这种一个负特征值的Hessian矩阵，也是MECP的充分条件。
（3）对于拉格朗日函数的Hessian阵，投影掉约束方向，在剩下的空间中，Hessian矩阵所有特征值全是正的，是约束优化问题中常见的充分条件表达。稍进一步，就有充要条件，可以作为MECP较为实用的判据。
按图解释下：
（1）在KKT点，拉格朗日函数的Hessian阵，如果特征值全部是正值，那么KKT点一定是MECP。1楼2.4节中的图（a）就是这种情况。 （您说的就是这种情况）
（2）在KKT点，拉格朗日函数的Hessian阵，如果在约束方向上有一个负本征值，仍然是MECP。1楼2.4节中图（b）的情况。（拉格朗日函数的Hessian阵有负本征值，但仍是MECP）
（3）在KKT点，拉格朗日函数的Hessian阵，如果在非约束方向上有一个负本征值，就不是MECP了。1楼2.4节中图（c）的情况。
（4）在KKT点，拉格朗日函数的Hessian阵，如果有多于一个的负本征值，那铁定不是MECP了。
其它：您的拉格朗日（ Lagrange）函数错写成lambda了，没有lambda函数。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |