本帖最后由 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），数学上是一个约束极值点。即：在“两个电子态能量相等”这个约束条件下，能量的极小值。
式中E1(q) 和 E2(q) 分别是两个态核坐标q的能量函数，约束条件写为c(q)= E1(q) – E2(q)。这里需要引入一个基本假定，即：函数E1(q) 和 E2(q) 、它们的一阶导数和二阶导数都是存在的和连续的。
通过构建拉格朗日函数，把约束极小值化为无约束极小值。
等式左边为拉格朗日函数，λ 是拉格朗日因子。这样，MECP最优性检验和稳定结构的最优性检验，就很类似。
2.2 必要条件 在满足线性独立约束条件（linear independence constrainted qualification，LICQ）的前提下，一阶必要条件描述如下：
式中 q* 和 λ* 分别是MECP候选解的坐标和拉格朗日因子。 对于MECP来说，第一个等式即拉格朗日函数对坐标的一阶导数为零；第二个等式表示两个态的能量相等。一阶必要条件称为Karush-Kuhn-Tucker（KKT）条件；满足一阶必要条件的候选解，称为KKT点。
2.3 充分条件
在满足KKT条件的前提下，二阶充分条件描述如下：
式中d 是切面M = {d : ∇c(q*) d = 0} 中的任意矢量，并且 d≠ 0。 L = ∇qql(q*,λ*) 是拉格朗日函数在KKT点处的二阶导数。 对于MECP来说，充分条件是拉格朗日函数对坐标的二阶导数，在其切面上的投影正定。
2.4 关于充分条件
对于必要条件，算法上实现非常简单。而对于充分条件，则有些复杂。情况如下图所示。
图中给出了低维下充分条件的几何图像。（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 输出文件 必要条件：
“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.”。
充分条件：
既约振动分析。有虚频表示该KKT点并不是真正的MECP，而是一个鞍点。
3.4 最终结果
左边图是在B3LYP/Aug-cc-PVTZ水平下，优化得到的两个异氰酸单三态的MECP候选点，图中给出了键长和能量（单位为埃和哈垂）。
右边图中给出了既约谐振分析（RVA）。右上图有负频率，说明KKT-1不是真正的极小点（不是MECP）；右下图没有负频率，KKT-2是真正的MECP。
既约振动分析（RVA）共有3N-7个振动，N是分子中原子个数。3N个运动中，6个平动和转动被分离出去（如果是线性MECP，则投掉5个，RVA是3N-6个振动），另一个被投影掉的是约束条件。在网址给出的文章中，提供了详细的计算细节。
需要注意的是，约束振动的模式，并不表示真正的（势能相关的）振动。但是，沿着负频率方向改变分子结构，对获取真正的MECP是有帮助的。
| 
发表于 Post on 2019-5-26 08:42:30
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author
