计算化学10-应用示例-分子几何结构优化

分子几何结构优化势能面的方程221()()()2nieeeNeiVVrErm分子的完全SchrÖdinger方程：Born-Oppenheimer近似后方程分解为核运动和电子运动两个方程：2221()(,)(,)22nieeeNNNiVVVRrERrm2221()()()22nieeeNNNiVVVRERm几何结构优化问题的数学描述0iiiVXFX势能面的势能函数:22120nieeeNNNelNNielNNiiiVVVVEVmEVVXXX势能函数求极值问题:分子几何结构优化的数学过程早期优化方法逐点优化法基于能量本身，计算量大，收敛慢，不利于程序化现代优化方法能量梯度法基于能量的一阶,二阶导数，更准确快速，易于程序化一维优化和多维优化问题22(,)2fxyxy22(,)2fxyxy梯度的概念方向导数的定义：函数z=f(x,y)在一点P沿某一方向的变化率0(,)(,)(,)limcossinfxyfxxyyfxylffflxy梯度的定义(,)coscossinsincossin(,)iffffgradfxyijxyxyjffffflxyxyiffijxyjgradfxylgrad22((2)2s4,)cogfxradxyxiyyj一维优化方法已知函数的解析形式和极小化条件可以使用Lagrange乘因子法无法知道函数的解析形式可以有如下两种方法使用二次函数拟合线性搜索可变尺度法22(,)2fxyxy22(2)24gradxyxiyj例:从(9,9)出发，使用Lagrange乘因子法求的梯度和梯度方向的极值点.步骤1(9,9)的梯度(18,36)2负梯度方向的下一点为(-9,-27)，该方向可以用函数表达y＝2x-93使用Largrange乘因子法可以确定该方向的极值点为(4,-1)线性搜索可变尺度法划界搜索法Newton(可变尺度)法a求得函数的近似导数：b沿着梯度方向寻找极小点：()()2iiiiifxxfxxx1iiixxs多维优化方法一阶导数法最陡下降法共轭梯度(方向)法二阶导数法Newton-Raphson方法准Newton方法最陡下降法基本原理从指定点出发，循梯度的负方向搜寻到极值点后作为新的起点，进行下一步搜寻例：函数从(9,9)出发，在(-18,-36)方向找到极值点(4,-1)后，a求该点的负梯度方向(-8,4),得到下一点为(-4,3)b得到该方向的方程y=-0.5x+1，c继续使用Lagrange乘因子法，求得该方向极值点(2/3,2/3)重复上述步骤，得到(0.296,-0.074)…………22(,)2fxyxy优点对于远离驻点的结构，优化效率非常高，能很快释放分子内的力缺点每一步都要进行直角转向，收敛慢，校正过度，振荡共轭梯度(方向)法基本原理做完一次线性搜索后，后一次优化的方向取该点的梯度方向与前一次优化的方向的组合22(,)2fxyxy222281880/9(8)(4)43620/9(18)(36)kv例：使用共轭梯度法求的极值点1起始点(9,9)的负梯度为(-18,-36),此方向极值点(4,-1)2点(4,-1)处的负梯度为(-8,4),搜索方向表达式为y=-1/4x,可以找到极值点(0,0)优点：对于有M个变量的函数，可以通过M步优化找到极值两种共轭梯度法纯的二次函数Fletcher-Reeves算法非纯二次函数Polak-Ribiere算法111()kkkkkkgggggNewton-Raphson方法将势能函数展开成Taylor级数****()()()()0()/()kkkkkkxxxxxxxxxxVVVVV()()kxxVV2()()()()()()/2kkkkkxxxxxxxxVVVV如果势能函数是纯二次函数，那么存在条件在势能极小点处，对函数求导，可以得到：22(,)2fxyxy1820(9,9)(9,9)3604vv例：求的极值点一阶导数Hessian矩阵*91/20180901/4360x优点对于纯二次函数，可以一步找到极值点缺点要求Hessian必须正定，否则将得到能量更高的坐标对于非纯二次函数，需要多步计算Hessian矩阵的计算量和存储量都非常大主要适用于小分子体系准Newton方法基本原理初猜一个Hessian矩阵，开始优化后，每步更新一次Hessian矩阵，每次更新Hessian矩阵都只使用上一步的Hessian矩阵和当点的一阶导数优化算法的选择算法的选择由多种因素决定1大分子体系多使用最陡下降法或者共轭梯度法2小分子多用Newton-Raphson法3对于远离驻点的结构，结合最陡下降法和Newton-Raphson法NMeNOHNMeNOHNH2NH2+NHONHNH2+H2N收敛判据1力最大力，均方根力2位移最大位移，均方根位移(36)TggNPeterPulay,AnalyticalDerivativeMethodsinQuantumChemistry,Adv.Chem.Phys.69,241(1987)(AbInitioMethodsinQuantumChemistryII,ed.K.P.Lawley(Wiley,1987))H.B.Schlegel,OptimizationofEquilibriumGeometriesandTransitionStructuresAdv.Chem.Phys.,67,249(1987).(AbInitioMethodsinQuantumChemistryI,ed.K.P.Lawley(Wiley1987))H.B.Schlegel,GeometryoptimizationonPotentialEnergySurfacesinModernElectronicStructureTheory,ed.D.R.Yarkony(WorldScientificPress,1995)H.B.Schlegel,“GeometryOptimization”inEncyclopediaofComputationalChemistry,ed.PvRSchleyer,NLAllinger,TClark,JGasteiger,PKollman,HFSchaeferPRSchreiner,(Wiley,Chichester,1998)关于分子几何结构优化方面的重要文献初猜分子几何结构和Hessian矩阵计算分子能量和梯度根据Hessian矩阵和梯度进行优化更新Hessian矩阵根据Hessian矩阵继续优化根据力和位移判断收敛性更新几何结构是完成否分子结构优化过程Gaussian中的分子结构优化过程与输出文件示例乙烯分子优化#PRHF/6-31G(d)OptTestEthyleneGeometryOptimization01CC1CCH1CH2HCCH1CH2HCC3180.H2CH1HCC3180.H2CH1HCC4180.Variables:CC=1.31CH=1.07HCC=121.5输入文件优化作业执行过程1/18=20,38=1/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;L101L103L202L301L302L303L401L502L601L701L702L703L716L9999优化开始的标志第一次做L103初猜Hessian矩阵优化结束的标志冗余内坐标R1R(2,1)1.5351R2R(3,1)1.4858R3R(4,2)1.4858R4R(4,3)1.4968R5R(5,1)1.0765R6R(6,1)1.0765R7R(7,2)1.0765R8R(8,2)1.0765A1A(2,1,3)89.26A2A(1,2,4)89.26A3A(1,3,4)90.74A4A(2,3,4)90.74A5A(2,1,5)115.76A6A(3,1,5)111.18A7A(2,1,6)115.76A8A(3,1,6)11.18A9A(5,1,6)111.65A10A(1,2,7)115.76A11A(4,2,7)111.18A12A(1,2,8)115.76A13A(4,2,8)111.18A14A(7,2,8)111.65程序自动生成根据共价半径确认成键（包括氢键和分子间键）所有键原子之间的键角所有键原子间的二面角D1D(4,2,1,3)0.00D2D(4,2,1,5)113.27D3D(4,2,1,6)-113.26D4D(7,2,1,3)113.27D5D(7,2,1,5)-133.45D6D(7,2,1,6)0.00D7D(8,2,1,3)-113.26D8D(8,2,1,5)0.00D9D(8,2,1,6)133.47D10D(4,3,1,2)0.00D11D(4,3,1,5)-117.47D12D(4,3,1,6)117.45D13D(3,4,2,1)0.00D14D(3,4,2,7)-117.47D15D(3,4,2,8)117.45D16D(2,4,3,1)0.00Gaussiane3_08优化时的坐标选项：OPT=RedundantOPT=AddRedundantOPT=Z-matrixOPT=Cartesian部分冻结优化：#HF/6-31G*Opt=ModRedundantTestPartialoptimization11lCH1R1H1R12A1O1R22A23120.0H4R33A32180.0A1=120.0...R3=1.1451.3F5432F优化收敛问题a默认的优化次数不够：增加优化次数opt=(RestartMaxcycle=N)b初猜力常数与实际不符:解决办法：1从振动分析计算的chk文件中读入力常数，Opt=ReadFc2使用给定的方法和基组计算力常数，Opt=CalcFc3优化中的每一步都计算力常数，Opt=CalcAllc从能量最低的最优结构重起优化作业Opt=CalcFc，Geom=(Check,Step=n)其他优化选项NoEigenTest:在使用Berny优化计算过渡态时，不做曲率测试NoFreeze:激活所有冻结的变量(Geom=Check)Tight,VeryTight,LooseIOP(1/8=m):改变步长0.01*m优化结果的确认计算完整的Hessian矩阵检查Hessian矩阵的负本征值极小点0个负本征值过渡态有且仅有1个负本征值寻找极小点如果有负本征值,循负本征值方向得到能量更低的结构寻找过渡态如果没有负本征值,沿最小本征值“上山”通过振动分析计算来确认优化结果复合作业：例：HF/6-31goptfreq振动分析与优化计算必须使用同样的方法和基组振动分析计算内坐标与简正坐标下Hessian矩阵的转换2,002,,022,1(