请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 1078|回复: 8

[算法与编程] LookForMECP: 一个验证极小势能面交叉点(MECP)充分条件的工具

[复制链接]

19

帖子

1

威望

435

eV
积分
474

Level 3 能力者

发表于 2019-5-26 08:42:30 | 显示全部楼层 |阅读模式
本帖最后由 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),数学上是一个约束极值点。即:在“两个电子态能量相等”这个约束条件下,能量的极小值。



约束优化问题.png



式中E1(q) 和 E2(q) 分别是两个态核坐标q的能量函数,约束条件写为c(q)= E1(q) – E2(q)。这里需要引入一个基本假定,即:函数E1(q) 和 E2(q) 、它们的一阶导数和二阶导数都是存在的和连续的。


通过构建拉格朗日函数,把约束极小值化为无约束极小值。


拉格朗日函数.png


等式左边为拉格朗日函数,λ 是拉格朗日因子。这样,MECP最优性检验和稳定结构的最优性检验,就很类似。




2.2  必要条件

在满足线性独立约束条件(linear independence constrainted qualification,LICQ)的前提下,一阶必要条件描述如下:

理论必要条件.png


式中 q*λ* 分别是MECP候选解的坐标和拉格朗日因子。

对于MECP来说,第一个等式即拉格朗日函数对坐标的一阶导数为零;第二个等式表示两个态的能量相等。一阶必要条件称为Karush-Kuhn-Tucker(KKT)条件;满足一阶必要条件的候选解,称为KKT点。


2.3  充分条件


在满足KKT条件的前提下,二阶充分条件描述如下:

      

    理论充分条件.png


式中d 是切面M = {d : c(q*) d = 0} 中的任意矢量,并且 d0L = qql(q*,λ*) 是拉格朗日函数在KKT点处的二阶导数。

对于MECP来说,充分条件是拉格朗日函数对坐标的二阶导数,在其切面上的投影正定。


2.4  关于充分条件

      对于必要条件,算法上实现非常简单。而对于充分条件,则有些复杂。情况如下图所示。


充分条件图示.png

   

图中给出了低维下充分条件的几何图像。(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  输出文件

必要条件:

必要条件.png

“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.”。



充分条件:
充分条件.png


既约振动分析。有虚频表示该KKT点并不是真正的MECP,而是一个鞍点。

3.4  最终结果

计算结果.png
    左边图是在B3LYP/Aug-cc-PVTZ水平下,优化得到的两个异氰酸单三态的MECP候选点,图中给出了键长和能量(单位为埃和哈垂)。
    右边图中给出了既约谐振分析(RVA)。右上图有负频率,说明KKT-1不是真正的极小点(不是MECP);右下图没有负频率,KKT-2是真正的MECP。
    既约振动分析(RVA)共有3N-7个振动,N是分子中原子个数。3N个运动中,6个平动和转动被分离出去(如果是线性MECP,则投掉5个,RVA是3N-6个振动),另一个被投影掉的是约束条件。在网址给出的文章中,提供了详细的计算细节。
    需要注意的是,约束振动的模式,并不表示真正的(势能相关的)振动。但是,沿着负频率方向改变分子结构,对获取真正的MECP是有帮助的。











评分

参与人数 6威望 +1 eV +25 收起 理由
ggdh + 5 精品内容
snljty + 5 牛!
一颗赛艇 + 5 GJ!
ABetaCarw + 5 大佬orz
zsu007 + 5 牛!
sobereva + 1

查看全部评分

114

帖子

0

威望

384

eV
积分
498

Level 3 能力者

发表于 2019-5-26 14:14:13 | 显示全部楼层
你这个方法还需要hessian 实用性大打折扣
很多高端的方法撑死只有梯度

19

帖子

1

威望

435

eV
积分
474

Level 3 能力者

 楼主| 发表于 2019-5-26 21:52:49 | 显示全部楼层
pyscf 发表于 2019-5-26 14:14
你这个方法还需要hessian 实用性大打折扣
很多高端的方法撑死只有梯度

理论上必须有Hessian的。实际上我也没找到好的技巧避开。

114

帖子

0

威望

384

eV
积分
498

Level 3 能力者

发表于 2019-5-26 22:40:34 | 显示全部楼层
bnulk 发表于 2019-5-26 21:52
理论上必须有Hessian的。实际上我也没找到好的技巧避开。

MECP涉及激发态的计算 普通的dft理论做激发态研究局限性很大
你要用好点的方法 多参考之类的 它们最多只支持到解析梯度

135

帖子

0

威望

768

eV
积分
903

Level 4 (黑子)

发表于 2019-5-27 11:27:25 | 显示全部楼层
pyscf 发表于 2019-5-26 22:40
MECP涉及激发态的计算 普通的dft理论做激发态研究局限性很大
你要用好点的方法 多参考之类的 它们最多 ...

都上多参考了谁还用dE搜MECP……

3

帖子

0

威望

775

eV
积分
778

Level 4 (黑子)

发表于 2019-6-10 16:08:26 | 显示全部楼层
按照你的方法做了一个尝试,提交任务以后出了问题,能帮忙看一下吗?
TIM图片20190610160357.png

19

帖子

1

威望

435

eV
积分
474

Level 3 能力者

 楼主| 发表于 2019-6-10 18:41:30 | 显示全部楼层
把输入文件发过来。或者发我邮箱:bnulk@foxmail.com

181

帖子

0

威望

1475

eV
积分
1656

Level 5 (御坂)

发表于 2019-10-19 06:19:32 | 显示全部楼层
请教楼主一个问题,
对于MECP来说,充分条件是拉格朗日函数对坐标的二阶导数,在其切面上的投影正定。

是不是指lambda函数的hessian所有的特征值全是正的?谢谢

19

帖子

1

威望

435

eV
积分
474

Level 3 能力者

 楼主| 发表于 2019-10-21 16:30:37 | 显示全部楼层
一颗赛艇 发表于 2019-10-19 06:19
请教楼主一个问题,

是不是指lambda函数的hessian所有的特征值全是正的?谢谢

    不是这个意思。

    分三个角度解释下:
    (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函数。

评分

参与人数 1eV +5 收起 理由
一颗赛艇 + 5 谢谢

查看全部评分

您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2019-11-17 19:11 , Processed in 0.172267 second(s), 28 queries .

快速回复 返回顶部 返回列表