计算化学及其应用.

计算化学及其应用势能面的描写DescriptionofPotentialEnergySurface势能面模型分子的势能面•我们的对象是分子,因此在势能面中,有意义的坐标只有3n-6个,6个对应于3个平动和3个转动自由度•对于双原子分子,有1个有意义的自由度,即其键长R0键长R二次曲线能量势能面扫描Gaussian程序有一个关键词scan,进行势能面扫描#RHF/STO-3GscanNOSYMWaterRHFscan0,1oh,1,rh,1,r,2,ar0.8550.05:变量起点步数步长a100.0101.00:总步数:(步数+1)势能面中的化学对象1.极小点体系的稳定结构2.一级鞍点稳定结构之间的过渡态3.极小点附近形状振动频率势能面的数学描写22211322231133()()nnnnEExxxEEExExxxxEHessian矩阵能量微商,对应于力的负值,E=Fx!化学对象的数学描写1.极小点体系的稳定结构一级微商=0二级微商本征值02.一级鞍点稳定结构之间的过渡态一级微商=0二级微商本征值0有一个0(虚频)3.极小点附近形状振动频率二级微商的质量加权本征值计算化学及其应用振动频率VibrationalFrequencyCalculations双原子分子的谐振子近似khxkxnuc21)2/1v(212ˆ2222H–振动能级的能量–振动频率0键长R二次曲线能量多原子分子的谐振子近似ki,j–笛卡尔坐标下的谐振子力常数(势能面的二阶微商)–质量加权的笛卡尔坐标jijijiiiijijijinucjijijijijiinucmmkkxmkxxREkxxkxmii,,,22,22,,22,2~~212ˆ)(212ˆHH多原子分子的谐振子近似I–质量加权的笛卡尔力常数矩阵的本征值qi–简正振动模式ijijiiiijinucmMqqi/2~212ˆ,,222,2MxLLqMLkMLLkLHtttt力常数矩阵及其本征值221/21/221111311112223130000()00000000()nnnEEmmxxxmmmmEExxxF•|FI|=03n个本征值i(i=1,3n)•其中有6个等于零,对应于3个平动和3个转动自由度•频率•如果本征值是负值,那么频率就变成虚数/2iiPople,J.A.;Schlegel,H.B.;Krishnan,R.;DeFrees,D.J.;Binkley,J.S.;Frisch,M.J.;Whiteside,R.A.;Hout,R.F.;Hehre,W.J.;Molecularorbitalstudiesofvibrationalfrequencies.Int.J.Quantum.Chem.,QuantumChem.Symp.,1981,15,269-278.振动频率的校正因子•计算得到的简正频率比实验值一般高10%•这是由于谐振子近似和理论的近似而产生的方法频率零点能HF/3-21G0.90850.9409HF/6-31G(d)0.89290.9135MP2/6-31G(d)0.94340.9676B3LYP/6-31G(d)0.96130.9804振动强度•振动强度用于光谱指认•IR光谱的振动带的强度由偶极矩对简正模式的微商确定•Raman光谱的振动带强度由极化率对简正模式的微商的平方确定振动频率的计算•振动频率只与极小点的附近有关系•可以用解析方法,套用公式把二级微商直接计算出来(解析方法)–G03对HF,DFT,MP2等都可以–运算速度快•也可以先把其附近点的能量算出来,用数值微商的方法计算二级微商,(数值方法)–对所有体系通用–对无法用解析法处理的,可以用它,比如较大分子的MP2频率,没有实现解析法的高级方法等–计算时间长–G03关键词:Freq=numer,附近一般是x0.001Å计算化学及其应用几何优化GeometryOptimization几何优化的目的•寻找势能面上的极小点,确定分子的可能的稳定结构•极小点满足的条件:220,0iiEEFxx几何优化算法的必要性•势能面随着分子中原子数目的增加而迅速增加,m3n个能量值,对中等体系的势能面都无法实际执行•可以给定一个初始的结构,按照力的方向去优化,把3n维的稳定点寻找变成近似一维的寻找•几何优化得到的仅仅是势能面上的局部极小点!能量寻找极小值的算法•单变量寻找,–收敛慢–只需要能量,不需要梯度•共轭梯度法或准牛顿法–较快收敛–需要梯度(用数值方法或解析方法计算)–Fletcher-Powell,DFP,MS,BFGS,OC•牛顿法–收敛迅速–要求二阶微商能量微商•可以用解析方法直接求一阶微商的有:–Hartree-Fock–DFT–Møller-Plesset微扰理论•MP2,MP3,MP4(SDQ)–组态相互作用方法,CIS,CID,CISD–CASSCF–耦合簇方法,CCSD和QCISD•可以用解析方法直接求二阶微商的有:–Hartree-Fock–DFT–MP2–CASSCF–CIS初始猜测几何结构和Hessian矩阵计算能量及其梯度沿着当前点和前一个点的方向得到一个极小值点更新Hessian矩阵(Powell,DFP,MS,BFGS,Berny,等等)从Hessian矩阵得到位移(Newton,RFO,Eigenvectorfollowing)在梯度和位移上判断是否收敛更新几何结构是完成否几何结构收敛的判据•精确的极小点位置是不可能找到的•只能逼近,Fi→0,Dxi→0•三种判据:–能量变化很小10-8Hartree–力很小:最大力0.00045,力的均方差0.00030–前后两次的坐标位移很小:最大位移0.0018,均方差位移0.0012Gaussian中的梯度法优化•初始猜测Hessian矩阵–在冗余内坐标下,用简单的价层力场得到Hessian矩阵的经验猜测(TCA66,333,(1984)•线性寻找极小值–按照当前和前一个函数的值以及梯度拟合一个限制二次曲线–即限制其二阶微商永远是正值–在二次曲线上取得极小点,并且用插值法计算出梯度•更新Hessian矩阵和位移–使用来自前一个点的梯度信息–用BFGS方法求极小点冗余内坐标•由程序自动生成•从笛卡尔坐标开始•按照共价半径来确认成键(检验氢键和分子片成键)•构造出成键原子之间的所有角(对接近直线的角构造特殊的直线弯曲坐标)•构造出成键原子之间的所有二面角(要考虑成平面的原子组)•估算出初始Hessian矩阵的对角元(包括氢键和分子片成键)冗余内坐标Dioxetane(HF/3-21G)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.65D1D(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.00几何优化效率比较内坐标,笛卡尔坐标,混合坐标,冗余内坐标NumberNumberofNumberofOptimizationStepsMoleculeofAtomsSymmetryVariablesInternalCartesianMixedRedundant2fluorofuran9Cs157776norbornane19C2v157555bicyclo[2.2.2]octane22D3111119147bicycol[3.2.1]octane22Cs336675endohydroxybicyclopentane14C136818912exohydroxybicyclopentane14C13610201111ACTHCP16C142658172281,4,5trihydroxyanthroquinone27Cs511011178histamineH+18C148421004619TAXOL113C158Kodakdye50C130检验极小值•计算整个Hessian矩阵(在优化过程中迭代的Hessian矩阵准确度不够,而且没有包括低对称性的信息)•检验负本征值的数目:–0个对应于极小点.–1个且只有1个对应于过渡态•如果要求极小值,而Hessian矩阵有一些负值,就沿着对应的本征矢量的方向求能量更低的结构.•如果要得到过渡态,而Hessian矩阵没有负本征值,就沿着最小的本征矢量寻找鞍点几何优化失败时的策略•步数超出–检查非常容易改变的坐标和/或强烈耦合的坐标–增加几何优化的最大步数OPT=(Restart,Maxcyc=N)•最大步长超出–如果这种情况经常发生,检查非常容易改变的坐标和/或强烈耦合的坐标•在优化过程中点群改变了–检查结构和/或不使用对称性(NoSymm)坐标方面的选项•Opt=Cartesian:在笛卡尔坐标下进行优化•Opt=Z-Matrix:在Z-矩阵下进行优化•Opt=Redundant:在冗余内坐标下进行优化(默认)•Opt=ModRedundant:增加或修改冗余内坐标–N1N2[N3[N4]][value][D|F|A|R]or[Hfc]•D-数值差分•F-冻结坐标•A-活化坐标•R-删除坐标•H-使用“fc”估算力常数对角元•可以对体系进行部分优化,只需要把不优化的部分冻结即可•对于感兴趣的是一些特殊键长,比如远距离氢键,范德华力等,程序无法自动把它包含到冗余坐标中,这时就必须手工把它添加上,以保证结构的可靠性OPT的更多选项•Maxcycle=n:设定优化的最大步数•NoEigenTest:在BernyTS优化时要求检验曲率•NoFreeze:活化所有冻结的坐标(constants).•Expert:在力常数和步长上放松限制•Tight,VeryTight:提高收敛标准(力和步长)•Loose:在工作初期使用•MaxStep=m:最大步长=0.01*mGEOM的选项•Checkpoint:从Chekpoint文件中读入分子数据(通常还用到Guess=Read).•Modify:读入并修改Ceheckpoint文件中的分子数据.•NoCrowd:允许原子间距小于0.5A.•NoKeep:丢弃关于冻结的信息.•Step=n:在前一次优化失败后,把它的第n步作为起点进行再次优化计算化学及其应用过渡态寻找SearchforTransitionStructures过渡态的目的•寻找势能面上的一级鞍点,确定分子可能的稳定结构之间的变换可能经过的态•一级鞍点满