《计算化学10-应用示例-分子几何结构优化ppt课件.ppt》由会员分享，可在线阅读，更多相关《计算化学10-应用示例-分子几何结构优化ppt课件.ppt（53页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、分子几何结构优化分子几何结构优化势能面的方程势能面的方程分子的完全分子的完全 Schrdinger方程：方程：Born-Oppenheimer近似后方程分解为核运动和电子运动两近似后方程分解为核运动和电子运动两个方程：个方程：几何结构优化问题的数学描述几何结构优化问题的数学描述势能面的势能函数势能面的势能函数:势能函数求极值问题势能函数求极值问题:分子几何结构优化的数学过程分子几何结构优化的数学过程 早期优化方法早期优化方法 逐点优化法逐点优化法 基于能量本身，计算量大，基于能量本身，计算量大，收敛慢，不利于程序化收敛慢，不利于程序化 现代优化方法现代优化方法 能量梯度法能量梯度法 基于能量的
2、一阶基于能量的一阶,二阶二阶导数，导数，更准确快速，易于程序化更准确快速，易于程序化一维优化和多维优化问题一维优化和多维优化问题梯度的概念梯度的概念 方向导数的定义：方向导数的定义：函数函数 z=f(x,y)在一点在一点 P 沿某一方向的变化率沿某一方向的变化率 梯度的定义梯度的定义一维优化方法一维优化方法 已知函数的解析形式和极小化条件已知函数的解析形式和极小化条件可以使用可以使用Lagrange乘因子法乘因子法 无法知道函数的解析形式无法知道函数的解析形式可以有如下两种方法可以有如下两种方法 使用二次函数拟合使用二次函数拟合 线性搜索可变尺度法线性搜索可变尺度法例例:从从(9,9)出发，使
3、用出发，使用Lagrange乘因子法求乘因子法求 的梯度和梯度方向的极值点的梯度和梯度方向的极值点.步骤步骤 1 (9,9)的梯度的梯度 (18,36)2 负梯度方向的下一点为负梯度方向的下一点为(-9,-27)，该方向可以用函数表达该方向可以用函数表达 y2x-9 3 使用使用Largrange乘因子法可以确定该方向的极值点为乘因子法可以确定该方向的极值点为 (4,-1)线性搜索可变尺度法线性搜索可变尺度法 划界搜索法划界搜索法 Newton(可变尺度可变尺度)法法a 求得函数的近似导数：求得函数的近似导数：b 沿着梯度方向寻找极小点：沿着梯度方向寻找极小点：多维优化方法多维优化方法一阶导数
4、法一阶导数法最陡下降法最陡下降法共轭梯度共轭梯度(方向方向)法法二阶导数法二阶导数法Newton-Raphson方法方法准准Newton方法方法最陡下降法最陡下降法基本原理基本原理从指定点出发，循梯度的负方向搜寻到极值点后作从指定点出发，循梯度的负方向搜寻到极值点后作为新的起点，进行下一步搜寻为新的起点，进行下一步搜寻例：函数例：函数 从从(9,9)出发，在出发，在(-18,-36)方向找到极值点方向找到极值点(4,-1)后，后，a 求该点的负梯度方向求该点的负梯度方向(-8,4),得到下一点为得到下一点为(-4,3)b 得到该方向的方程得到该方向的方程 y=-0.5x+1，c 继续使用继续使
5、用Lagrange乘因子法，求得该方向极值点乘因子法，求得该方向极值点(2/3,2/3)重复上述步骤，得到重复上述步骤，得到(0.296,-0.074)优点优点对于远离驻点的结构，优化效率非常高，对于远离驻点的结构，优化效率非常高，能很快释放分子内的力能很快释放分子内的力缺点缺点每一步都要进行直角转向，收敛慢，校正每一步都要进行直角转向，收敛慢，校正过度，振荡过度，振荡共轭梯度共轭梯度(方向方向)法法基本原理基本原理做完一次线性搜索后，后一次优化的方向取该点做完一次线性搜索后，后一次优化的方向取该点的梯度方向与前一次优化的方向的组合的梯度方向与前一次优化的方向的组合例：使用共轭梯度法求例：使用
6、共轭梯度法求 的极的极值点值点1 起始点起始点(9,9)的负梯度为的负梯度为(-18,-36),此方向极值点此方向极值点(4,-1)2 点点(4,-1)处的处的负负梯度为梯度为(-8,4),搜索方向表达式为搜索方向表达式为y=-1/4x,可以找到极值点可以找到极值点(0,0)优点：优点：对于有对于有M个变量的函数，可以通过个变量的函数，可以通过M步优化找到极值步优化找到极值两种共轭梯度法两种共轭梯度法纯的二次函数纯的二次函数 Fletcher-Reeves算法算法非纯二次函数非纯二次函数 Polak-Ribiere算法算法Newton-Raphson方法方法将势能函数展开成将势能函数展开成Ta
7、ylor级数级数如果势能函数是纯二次函数，那么存在条件如果势能函数是纯二次函数，那么存在条件在势能极小点处，对函数求导，可以得到：在势能极小点处，对函数求导，可以得到：例：求例：求 的极值点的极值点一阶导数一阶导数 Hessian矩阵矩阵优点优点对于纯二次函数，可以一步找到极值点对于纯二次函数，可以一步找到极值点缺点缺点要求要求Hessian必须正定，否则将得到能量更高的坐标必须正定，否则将得到能量更高的坐标对于非纯二次函数，需要多步计算对于非纯二次函数，需要多步计算Hessian矩阵的计算量和存储量都非常大矩阵的计算量和存储量都非常大主要适用于小分子体系主要适用于小分子体系准准Newton方
8、法方法 基本原理基本原理 初猜一个初猜一个Hessian矩阵，开始优化后，每步更新一矩阵，开始优化后，每步更新一 次次Hessian矩阵，每次更新矩阵，每次更新Hessian矩阵都只使用上矩阵都只使用上 一步的一步的Hessian矩阵和当点的一阶导数矩阵和当点的一阶导数优化算法的选择优化算法的选择算法的选择由多种因素决定算法的选择由多种因素决定1 大分子体系多使用最陡下降法或者共轭梯度法大分子体系多使用最陡下降法或者共轭梯度法2 小分子多用小分子多用Newton-Raphson法法3 对于远离驻对于远离驻点的结构，结合最陡下降法和点的结构，结合最陡下降法和Newton-Raphson法法收敛判
9、据收敛判据1 力力 最大力，均方根力最大力，均方根力2 位移位移 最大位移，均方根位移最大位移，均方根位移lPeter Pulay,Analytical Derivative Methods in Quantum Chemistry,Adv.Chem.Phys.69,241(1987)(Ab Initio Methods in Quantum Chemistry II,ed.K.P.Lawley(Wiley,1987)lH.B.Schlegel,Optimization of Equilibrium Geometries and Transition Structures Adv.Chem.P
10、hys.,67,249(1987).(Ab Initio Methods in Quantum Chemistry I,ed.K.P.Lawley(Wiley 1987)lH.B.Schlegel,Geometry optimization on Potential Energy Surfaces in Modern Electronic Structure Theory,ed.D.R.Yarkony(World Scientific Press,1995)lH.B.Schlegel,“Geometry Optimization”in Encyclopedia of Computational
11、 Chemistry,ed.PvR Schleyer,NL Allinger,T Clark,J Gasteiger,P Kollman,HF Schaefer PR Schreiner,(Wiley,Chichester,1998)关于分子几何结构优化方面的重要文献关于分子几何结构优化方面的重要文献初猜分子几何结构和初猜分子几何结构和Hessian矩阵矩阵计算分子能量和梯度计算分子能量和梯度根据根据Hessian矩阵矩阵和梯度进行和梯度进行优化优化更新更新Hessian矩阵矩阵根据根据Hessian矩阵继续优化矩阵继续优化根据力和位移判断收敛性根据力和位移判断收敛性 更新几何结构更新几何结构
12、是是完成完成否否分子结构优化过程分子结构优化过程Gaussian中的分子结构优化过程中的分子结构优化过程与输出文件示例与输出文件示例乙烯分子优化乙烯分子优化乙烯分子优化乙烯分子优化#P RHF/6-31G(d)Opt Test Ethylene Geometry Optimization 0 1 C C 1 CC H 1 CH 2 HCC H 1 CH 2 HCC 3 180.H 2 CH 1 HCC 3 180.H 2 CH 1 HCC 4 180.Variables:CC=1.31 CH=1.07 HCC=121.5输入文件输入文件优化作业执行过程优化作业执行过程 1/18=20,38=1
13、/1,3;2/9=110,17=6,18=5,40=1/2;3/5=1,6=6,7=1,11=1,16=1,25=1,30=1/1,2,3;4/7=1/1;5/5=2,38=5/2;6/7=2,8=2,9=2,10=2,28=1/1;7/1,2,3,16;1/18=20/3(1);99/99;2/9=110/2;3/5=1,6=6,7=1,11=1,16=1,25=1,30=1/1,2,3;4/5=5,7=1,16=3/1;5/5=2,38=5/2;7/1,2,3,16;1/18=20/3;2/9=110/2;6/7=2,8=2,9=2,10=2,19=2,28=1/1;99/9=1/99;L
14、101L103L202L301L302L303L401L502L601L701L702L703L716L9999优化开始的标志优化开始的标志第一次做第一次做 L103初猜初猜Hessian矩阵矩阵优化结束的标志优化结束的标志冗余内坐标冗余内坐标程序自动生成程序自动生成根据共价半径确认成键（包括氢键和分子间键）根据共价半径确认成键（包括氢键和分子间键）所有键原子之间的键角所有键原子之间的键角所有键原子间的二面角所有键原子间的二面角Gaussian e 3_08Gaussian e 3_08优化时的坐标选项：优化时的坐标选项：OPT=Redundant OPT=AddRedundant OPT=Z
15、-matrix OPT=Cartesian部分冻结优化：部分冻结优化：#HF/6-31G*Opt=ModRedundant TestPartial optimization1 1 lCH 1 R1H 1 R1 2 A1O 1 R2 2 A2 3 120.0H 4 R3 3 A3 2 180.0A1=120.0.R3=1.14 5 1.3 F5 4 3 2 F优化收敛问题优化收敛问题a 默认的优化次数不够：默认的优化次数不够：增加优化次数增加优化次数 opt=(Restart Maxcycle=N)b 初猜力常数与实际不符初猜力常数与实际不符:解决办法：解决办法：1 从振动分析计算的从振动分析计
16、算的chk文件中读入力常数，文件中读入力常数，Opt=ReadFc2 使用给定的方法和基组计算力常数，使用给定的方法和基组计算力常数，Opt=CalcFc3 优化中的每一步都计算力常数，优化中的每一步都计算力常数，Opt=CalcAllc 从能量最低的最优结构重起优化作业从能量最低的最优结构重起优化作业 Opt=CalcFc，Geom=(Check,Step=n)其他优化选项其他优化选项NoEigenTest:在使用在使用Berny优化计算过渡态时，优化计算过渡态时，不做曲不做曲率测试率测试NoFreeze:激活所有冻结的变量激活所有冻结的变量(Geom=Check)Tight,VeryTig
17、ht,LooseIOP(1/8=m):改变步长改变步长 0.01*m优化结果的确认优化结果的确认计算完整的计算完整的Hessian矩阵矩阵检查检查Hessian矩阵的负本征值矩阵的负本征值极小点极小点 0 个负本征值个负本征值 过渡态过渡态 有且仅有有且仅有 1 个负本征值个负本征值 寻找极小点寻找极小点如果有负本征值如果有负本征值,循负本征值方向得到能量更低的结构循负本征值方向得到能量更低的结构 寻找过渡态寻找过渡态 如果没有负本征值如果没有负本征值,沿最小本征值沿最小本征值“上山上山”通过振动分析计算来确认优化结果通过振动分析计算来确认优化结果复合作业：复合作业：例：例：HF/6-31g
18、opt freq振动分析与优化计算必须使用同样的方法和基组振动分析与优化计算必须使用同样的方法和基组振动分析计算振动分析计算内坐标与简正坐标下内坐标与简正坐标下Hessian矩阵的转换矩阵的转换内坐标下的势能表示内坐标下的势能表示平衡点的势能规定为零平衡点的势能规定为零核运动的哈密顿核运动的哈密顿交叉项存在，无法分离变量求解交叉项存在，无法分离变量求解简正坐标和简正振动模式简正坐标和简正振动模式核运动的哈密顿核运动的哈密顿 质权坐标质权坐标再找到一个再找到一个U阵将阵将K对角化，可以得到核运动的哈密顿为对角化，可以得到核运动的哈密顿为通解为通解为通解为通解为据此能求算相关的热化学性质据此能求算
19、相关的热化学性质据此能求算相关的热化学性质据此能求算相关的热化学性质例例例例 1595(1826)cm-1 3756(4188)cm-1 3652(4070)cm-1 1.145 1.115 1.114 HF/6-31G*(实验实验)优化收敛问题优化收敛问题 (2)问题问题引起原因引起原因解决方案解决方案力过大力过大(a)初始构型太差初始构型太差(b)分子坐标选择得不好分子坐标选择得不好(a)改变分子构型重新输入改变分子构型重新输入(b)重新构造分子坐标重新构造分子坐标优化结构的优化结构的Hessian矩阵出现负本征值矩阵出现负本征值(a)结构没有优化到极小点结构没有优化到极小点(b)Hess
20、ian矩阵更新中数矩阵更新中数值计算错误值计算错误(a)沿负本征值方向优化得到能量更沿负本征值方向优化得到能量更低的结构低的结构(b)改变方法初猜改变方法初猜Hessian矩阵重新矩阵重新优化优化(c)如果如果(b)仍不行，冻结负本征值仍不行，冻结负本征值所在的坐标矢量优化收敛后再放开所在的坐标矢量优化收敛后再放开优化优化Hessian矩阵的本征矩阵的本征值太大值太大(少见少见)(a)Hessian矩阵初猜不好矩阵初猜不好(b)Hessian矩阵更新不好矩阵更新不好(c)不同坐标之间偶合过强不同坐标之间偶合过强(a)(b)重新初猜重新初猜Hessian矩阵矩阵(c)重新构造分子坐标重新构造分子
21、坐标优化步数太多优化步数太多(a)优化比较困难，需要的优化比较困难，需要的步数确实很多步数确实很多(b)输入内坐标中存在冗余输入内坐标中存在冗余(c)存在很松弛的坐标存在很松弛的坐标(d)不同坐标之间偶合过强不同坐标之间偶合过强(e)初猜不好初猜不好(a)增加优化次数增加优化次数 maxcyc=n(b)重新输入内坐标重新输入内坐标(c)冻结松弛的坐标优化冻结松弛的坐标优化收敛后再放开优化收敛后再放开优化(a)重新构造分子坐标重新构造分子坐标优化中的对称性问题优化中的对称性问题优化到高对称性优化到高对称性一般可以做到，但是收敛非常慢，而且受算法的限制，一般可以做到，但是收敛非常慢，而且受算法的限
22、制，只能得到近似的结果只能得到近似的结果*解决办法：提高到相应对称性直接优化然后作比较解决办法：提高到相应对称性直接优化然后作比较优化到低对称性优化到低对称性优化中计算得到的梯度矢量属于分子点群的对称表示，优化中计算得到的梯度矢量属于分子点群的对称表示，优化过程中分子对称性不会降低优化过程中分子对称性不会降低*解决办法：检查解决办法：检查Hessian矩阵，依负本征值所在矢量方向矩阵，依负本征值所在矢量方向降低对称性降低对称性1 2 31 O 7.839428 0.266565 0.2665652 H 0.266565 0.591542-0.0443853 H 0.266565-0.04438
23、5 0.591542Total atomic charges:11 O-0.3725582 H 0.186279（Mulliken 布居分析给出的原子净电荷）3 H 0.186279Sum of Mulliken charges=0.00000Atomic charges with hydrogens summed into heavy atoms:11 O 0.0000002 H 0.0000003 H 0.000000Sum of Mulliken charges=0.00000Electronic spatial extent(au):=17.8523Charge=0.0000 elec
24、tronsDipole moment(Debye):X=0.0000 Y=0.0000 Z=-1.6887 Tot=1.6887（永久偶极矩及其分量）Quadrupole moment(Debye-Ang):XX=-6.0922 YY=-4.1956 ZZ=-5.4610（四极矩）XY=0.0000 XZ=0.0000 YZ=0.0000Octapole moment(Debye-Ang*2):XXX=0.0000 YYY=0.0000 ZZZ=-0.1750 XYY=0.0000XXY=0.0000 XXZ=-0.0058 XZZ=0.0000 YZZ=0.0000（八极矩）YYZ=-0.5
25、891 XYZ=0.0000Hexadecapole moment(Debye-Ang*3):XXXX=-3.2295 YYYY=-6.5600 ZZZZ=-4.7118 XXXY=0.0000XXXZ=0.0000 YYYX=0.0000 YYYZ=0.0000 ZZZX=0.0000（十二极矩）15 ZZZY=0.0000 XXYY=-1.8140 XXZZ=-1.3520 YYZZ=-1.7004XXYZ=0.0000 YYXZ=0.0000 ZZXY=0.0000N-N=9.157240288768D+00 E-N=-1.969268663071D+02 KE=7.4584226988
26、61D+01Symmetry A1 KE=6.662704811097D+01Symmetry A2 KE=0.000000000000D+00Symmetry B1 KE=5.057462452019D+00Symmetry B2 KE=2.899716425620D+001|1|UNPC-UNK|SP|RHF|STO-3G|H2O1|PCUSER|06-Dec-2001|0|#RHF/STO-3G|Water Energy|0,1|O,0,0.,0.,0.|H,0,0.784,-0.554,0.|H,0,-0.784,-0.554,0.|Version=x86-Win32-G98RevA.9|State=1-A1|HF=-74.9606932|RMSD=3.693e-004（光谱项和总能量）|Dipole=0.,-0.664381,0.|PG=C02V C2(O1),SGV(H2)|Job cpu time:0 days 0 hours 0 minutes 30.0 seconds.File lengths(MBytes):RWF=10 Int=0 D2E=0 Chk=5 Scr=1Normal termination of Gaussian 98