查看文章 |
几何优化 研究化学反应
2008-03-26 23:59
一、优化
前面讨论了在特定几何构型下的能量的计算,可以看出,分子几何构型的变化对能
量有很大的影响.由于分子几何构型而产生的能量的变化,被称为势能面.势能面 是连接几何构型和能量的数学关系.对于双原子分子,能量的变化与两原子间的距 离相关,这样得到势能曲线,对于大的体系,势能面是多维的,其维数取决与分子的 自由度. 3.1势能面 势能面中,包括一些重要的点,包括全局最大值,局域极大值,全局最小值,局域极小 值以及鞍点.极大值是一个区域内的能量最高点,向任何方向的几何变化都能够引起 能量的减小.在所有的局域极大值中的最大值,就是全局最大值;极小值也同样,在 所有极小之中最小的一个就是具有最稳定几何结构的一点.鞍点则是在一个方向上 具有极大值,而在其他方向上具有极小值的点. 一般的,鞍点代表连接着两个极小值的过渡态. 寻找极小值 几何优化做的工作就是寻找极小值,而这个极小值,就是分子的稳定的几何形态. 对于所有的极小值和鞍点,其能量的一阶导数,也就是梯度,都是零,这样的点被称为 稳定点.所有的成功的优化都在寻找稳定点,虽然找到的并不一定就是所预期的点. 几何优化有初始构型开始,计算能量和梯度,然后决定下一步的方向和步长,其方向 总是向能量下降最快的方向进行.大多数的优化也计算能量的二阶导数,来修正力矩 阵,从而表明在该点的曲度. 3.2 收敛标准 当一阶导数为零的时候优化结束,但实际计算上,当变化很小,小于某个量的时候, 就可以认为得到优化结构.对于Gaussian,默认的条件是 力的最大值必须小于0.00045,均方根小于0.0003 为下一步所做的取代计算为小于0.0018,其均方根小于0.0012 这四个条件必须同时满足,比如,对于非常松弛的体系,势能面很平缓,力的值已经小 于域值,但优化过程仍然有很长的路要走.对于非常松弛的体系,当力的值已经低于 域值两个数量级,尽管取代计算仍然高于域值,系统也认为找到了最优点.这条规则 用于非常大,非常松弛的体系. 3.3 几何优化的输入 Opt关键字描述了几何优化 例3.1 文件 e3_01 乙烯的优化 输入文件的设置行为 #R RHF/6-31G(d) Opt Test 表明采用RHF方法,6-31G(d)基组进行优化 3.4 输出文件 优化部分的计算包含在两行相同的 GradGradGradGradGradGradGradGradGradGrad........... 之间,这里有优化的次数,变量的变化,收敛的结果等等.注意这里面的长度单位是 波尔. 在得到每一个新的几何构型之后,都要计算单点能,然后再在此基础上继续进行优 化,直到四个条件都得到满足.而最后一个几何构型就被认为是最优构型. 注意,最终构型的能量是在最后一次优化计算之前得到的. 在得到最优构型之后,在文件中寻找 --Stationmay point found. 其下面的表格中列出的就是最后的优化结果以及分子坐标. 随后按照设置行的要求,列出分子有关性质 例3.2 文件 e3_02 氟代乙烯的优化 3.5 寻找过渡态 Gaissian使用STQN方法确定反应过渡态,关键词是Opt=QST2 例3.3 文件 e3_03 过渡态优化 例中分析的是H3CO --> H2COH 的变化,输入文件格式 #T UHF/6-31G(d) Opt=QST2 Test H3CO --> H2COH Reactants 0,2 structure for H3CO 0,2 structure for H2COH Gaussian也提供QST3方法,可以优化反应物,产物和一个由用户定义的猜测的过渡态. 3.6 难处例的优化 有一些系统的优化很难进行,采用默认的方法得不到结果,其产生的原因往往是所计 算出的力矩阵与实际的相差太远.当默认方法得不到结果时,就要采用其他的方法. Gaussian提供很多的选择,具体可以看User's Reference.下面列举一些. Opt=ReadFC 从频率分析(往往是采用低等级的计算得到的)所得到的checkpoint文件 中读取初始力矩阵,这一选项需要在设置行之前加入%Chk= filename 一句,说明文 件的名称. Opt=CalCFC 采用优化方法同样的基组来计算力矩阵的初始值. Opt=CalcAll 在优化的每一步都计算力矩阵.这是非常昂贵的计算方法,只在非常 极端的条件下使用. 有时候,优化往往只需要更多的次数就可以达到好的结果,这可以通过设置MaxCycle 来实现.如果在优化中保存了Checkpoint文件,那么使用Opt=Restart可以继续所进 行的优化. 当优化没有达到效果的时候,不要盲目的加大优化次数.这是注意观察每一步优化的 区别,寻找没有得到优化结果的原因,判断体系是否收敛,如果体系能量有越来越小 的趋势,那么增加优化次数是可能得到结果的,如果体系能量变化没有什么规律,或 者,离最小点越来越远,那么就要改变优化的方法. 也可以从输出文件的某一个中间构型开始新的优化,关键词Geom=(Check,Step=n) 表示在取得在Checkpoint文件中第n步优化的几何构型 3.7 练习 练习3.1 文件 3_01a (180), 3_01b (0) 丙烯的优化 从两种丙烯的几何异构体进行优化,一个是甲基的一个氢原子与CCH形成180度二 面角,另一个是0.优化结果表明,二者有0.003Hartree的差别,0度的要低. 练习3.2 文件 3_02a (0), 3_02b (180), 3_02c (acteald.) 乙烯醇的优化 乙 醇氧端的原子与OCC平面的二面角可以为0和180,优化得到的结果时, 0度的能量比180度的低0.003Hartree,但同时做的乙醛的优化表明,乙醛的能量 还要低,比0度异构体低0.027hartree. 练习3.3 文件 3_03 乙烯胺的优化 运行所有原子都在同一平面上的乙烯胺的优化. 比较本章的例子和练习,可以看到不同取代基对乙 ┨ 碳双键的影响. 练习3.4 文件 3_04 六羰基铬的优化 本例采用STO-3G和3-21G基组,在设置行中加入SCF=NoVarAcc对收敛有帮助. 3-21G基组的优化结果要优于STO-3G 练习3.5 文件 3_05a (C6H6), 3_05b (TMS) 苯的核磁共振 采用6-31G(d)基组,B3LYP方法优化几何构性,采用HF方法,6-311+G(2d,p)基组在优 化的几何构型基础上计算碳的化学位移.注意,核磁共振的可靠程度依赖准确的几 何结构和大的基组.输入文件如下 %Chk=NMR %Chk=NMR #T B3LYP/6-31G(d) Opt Test Opt molecule specification --Link1-- %Chk=NMR %NoSave #T RHF/6-311+G(2d,p) NMR Geom=Check Guess=Read Test NMR charg & spin 同样,还需要采用同样方法计算TMS.下面是计算结果 绝对位移 相对位移 实验值 TMS Benzene 188.7879 57.6198 131.2 130.9 练习3.6 文件 3_06a (PM3), 3_06b (STO-3G) 氧化碳60的优化 C60中有两种碳碳键,一是连接两个六元环的6-6键,另一是连接六元环和无元环的5-6键. 氧化C60就有两种异构体.本例采用PM3和HF/STO-3G方法来判断那种异构体是稳定的,以 及氧化后的C-C键的变化. 采用Opt=AddRedundant关键词可以在输出文件中打印所要求的键长,键角,这一关键词 需要在分子构型输入结束后在增加关于所要键长键角的信息,键长用两个原子的序列号 表示,键角则用三个原子表示. 计算结果显示,6-6键的氧化,碳碳键仍然存在,接近环氧化合物,而5-6键已经打开. 采用不同的方法,得到的几何结构相差不多,但在能量上有很大差异. 在采用MNDO,PM3,HF/3-21G方法得到的能量数据中,5-6键氧化的异构体的能量低,但采用 HF/STO-3G得到的结果,确实6-6键氧化的能量低. Raghavachari在其进行的上述研究中阐述 动力学因素同样是重要的; 实验上还没有发现那个是能量最低的异构体; 应该进行更精确的计算 练习3.7 文件 3_07 一个1,1消除反应的过渡态优化 分析反应 SiH4 --> SiH2 + H2, 可以采用Opt=(QST2, AddRedundant)关键词来进行过 渡态优化,同时特别关注过渡态结构中的某个键长 练习3.8 文件 3_08 优化进程比较 采用下述三种方法优化二环[2,2,2] 直接采用默认方式冗余内坐标优化 Opt; 采用迪卡尔坐标优化 Opt=Cartesian; 采用内坐标优化 Opt=Z-Matrix 结果显示,冗余内坐标优化的优化次数最短,内坐标优化的次数最多. 二 、研究化学反应和反应性 本章讨论应用电子结构理论研究化学反应.我们将从电子密度开始,然后 回顾第四章中有关反应势垒的讨论,再讨论反应研究中的更复杂的技术, 最后,通过对相应反应的计算,来研究未知体系的反应热. 本章将引入两种新的计算方法 * 势能面 * 反应路径分析 8.1 预测电子密度 将电子密度或静电势可视化是研究一个分子体系的反应性的重要的第一步. 例8.1 文件 e8_01a, e8_01b 取代苯的电子密度 在有机化学中,亲电芳香取代反应的定位效应是已经被深入研究的课题. 在这里,我们采用电子密度对这一现象进行研究. 已经知道氯苯和硝基苯的硝化是基于同样的反应机理:苯环首先受NO2+的 攻击,产生各种异构体的阳离子异构体.当硝化完成后,产物分布如下. 邻位 间位 对位 氯硝基苯 29% 1% 70% 二硝基苯 7% 88% 1% 我们在这里检验间位和对位异构体的中间体. 分子采用B3LYP/6-31G(d)进行优化,电子密度在HF/6-31G(d)等级计算.将 电子密度按照平行苯环平面的方向切片,得到不同厚度位置的电子密度图. 间位的氯硝基苯和对位的二硝基苯的电子密度分布显示,其保留了有较大 共振范围的电子结构,相反,另两个构型的电子密度分布显示其电子分布相 对局域化,并且向苯环外的方向集中. 通过电子密度的图形,可以定性的理解电子密度和反应性的关系,在得到结论 之前,检查这个体积的电子密度是必要的.关于这方面的进一步资料可以参见 Gaussian出版的白皮书Visualizing Results from Gaussian. 8.2 计算反应焓变 例8.2 文件e8_02 水解反应 现在分析水解反应 H+ + H2O --> H3O+ 目的是计算标准反应焓变dH298.其计算方法可以表示为 dH298 = dE298 + d(PV) dE298 = dEe0 + d(dEe)298 + dEv0 + d(dEv)298 + dEr298 + dEt298 其中 dEe0: 0K时产物与反应物的能量差; d(dEe)298: 0K到298K电子能量的变化.对于这个反应,这一项可以忽略; dEv0: 0K时反应物和产物的零点能之差; d(dEv)298: 0K到298K振动能量的变化; dEr298: 产物和反应物的旋转能之差; dEt298: 产物和反应物的平动能之差; d(PV): 由于有一摩尔分子消失,PV=-RT. dEe0由单点能得到,本例采用的计算方法是B3LYP/6-311+G(2df,2p).其他的各项 都要考虑内能校正,通过频率分析得到.这样,所要做的工作就是进行优化然后进 行频率分析得到所需数值.采用B3LYP/6-31G(d)就能够得到足够精确的结果. 这里注意我们不用计算H+,由于没有电子,它的电子能量显然是0;由于只有一个 原子,其振动,转动能显然也是零,这样,其只有平动能,其值为 1.5RT = 0.889kcal.mol.(详见统计热力学). 最终计算得到dH298=-163.3kcal.mol.实验值为-165.3+-1.8kcal/mol. 两者符合的相当好. 8.3 研究势能面 对于势能面的研究对反应路径分析来讲,可能产生出人意料的好的结果.本节中, 我们通过实例研究势能面的应用方法. 考虑丙烯基正离子的旋转异构体的变化, H1 H1 | | H2a C H2b H2a C H3b \ / \ / <---> \ / \ / C C C C | | | | H3a H3b H3a H2b (I) (I') 曾经认为两个异构体之间的变化是通过一个具有Cs对称性的过渡态完成的,在该 构型中,H2b-C-H3b组成的平面垂直与碳原子平面.采用HF/6-311++G(d,p)能够找 到这样的过渡态,但是进一步的采用MP2和QCISD以及同样基组的研究却没有得到 过渡态,而得到了极小值!这个新的具有Cs对称性的结构中,H1迁移到了端位的碳 原子上.这个新结构的能量比势能面中平衡结构的能量高10kcal/mol. 这样就有了另一条反应路线: 中间碳上的氢迁移到端位的碳原子上; 新形成的甲基旋转; 旋转后的甲基上的一个氢原子迁移回中间碳原子. 在这个例子中,应用了IRC计算来确定过渡态的确是连接产物与反应物的. 本章后面将对这一方法进行讨论. HF方法的研究得到了假的过渡态,原因是,由于HF方法本身的限制,其计算的亚甲 基旋转的势垒要低于氢原子迁移的势垒. 8.4 势能面扫描 势能面扫描可以研究一个区域内的势能面.一般的扫描都是由一系列的在不同结 构上的单点能计算组成的.当进行势能面扫描时,要设置分子结构的变量,设置需 要变化的结构的范围和步长. 在Gaussian中,势能面扫描是自动进行的,下面是一个进行势能面扫描的算例. #T UMP4/6-311+G(d,p) Scan Test CH PES Scan 0 2 C H 1 R R 0.5 40 0.04 该算例要求一个对于CH的势能面扫描,所用的关键词是scan,变量的设置格式是: 名称 初始值 [点数 步长] 当只有一个参数时,变量在整个扫描中是不变的,当三个参数都设定时,变量将在 一定范围内变化. 当有多个变量时,所有的可能构型都要计算. 所有等级的计算结果都在输出文件中列出,比如进行的MP2的势能面扫描也将列出 HF方法的结果. 根据得到的扫描结果,可以得到所要的势能面,通过它,可能得到极小值的可能位 置.势能面扫描过程中不进行几何优化. 8.5 反应路径分析 在第四章中我们提到,得到一个过渡态机构不能说明它就是连接产物和反应物的 结构.分析其是不是所需过渡态的一个方法是分析虚频的简正振动状态.有时,对 振动的分析也不能够确定.本节讨论更为精确的方法. IRC方法检验过渡态分子的趋势.计算从过渡态开始,根据能量降低的方向来寻找 极小值,就是说,寻找过渡态所连接的两个极小值. 反应路径是连接反应物与产物的,但是连接反应物和产物可以有不止一条路径, 通过不同的过渡态连接,通过IRC计算,我们可以寻找真正的反应路径,也就是能量 最低的反应路径. 反应路径计算可以确认得到的过渡态就是连接反应物和产物的过渡态,一旦确认, 还可以计算活化能(注意零点能校正). 运行IRC 在Gaussian中,运行IRC的关键词是IRC.需要注意的是,IRC计算是从过渡态开始的, 在两个反应方向上各进行固定步骤的计算(默认是6步). IRC计算的方法是这样的: * 优化过渡态 * 进行频率分析,确认所得到的是过渡态,计算零点能,生成进行IRC计算的力矩阵 * 运行IRC,在鞍点的能量下降方向,寻找极小值.一般的,需要增大寻找的次数, 从而尽可能的接近极小值.方法是设置MaxPoints. 确定反应势垒,一般还要进行更多的工作, * 对过渡态的高等级的能量计算 * 对反应物和产物进行优化和频率分析,得到零点能,进行高等级能量计算 8.6 势能面研究实例 我们现在用Gaussian的反应路径分析来研究甲醛的势能面.这个势能面上有很多 极小值,包括甲醛,羟基卡宾,以及H2和CO.每一组之间都可以组合成不同的反应物 产物对.这里研究两个反应 H2CO <--> CO + H2 H2CO <--> HCOH 甲醛的解离 我们要确定反应过渡态的结构,预测反应的活化能.为此,我们需要以下信息: * 甲醛,氢分子,一氧化碳分子的考虑零点能的能量. * 过渡态的几何构型和零点能校正的能量. 计算在HF/6-31G(d)水平进行,结果如下 SCF能量 零点能 总能量 H2 -1.12683 0.00968 -1.11716 CO -112.73788 0.00508 -112.7280 H2 + CO -113.84996 H2CO -113.86633 0.02668 -113.83966 计算过渡态的能量,方法是 * 过渡态几何构型优化,计算SCF能量 * 频率分析,计算零点能 * IRC计算,确认过渡态 优化过渡态 例8.3 文件 e8_03 CH2O --> H2 + CO IRC 首先考虑氧原子垂直于CHH平面的构型,同时增大OCH夹角.计算中设置 Opt=(TS,CalcFC). CalcFC一般对于过渡态的优化是有帮助的. 得到的该点几何构型与猜测的结构接近,SCF能量-113.69352 频率分析 频率分析表明其有一个虚频,零点能0.01774(校正后),总能量-113.68578 IRC计算 IRC计算需要优化好的过渡态和相应的力矩阵,得到的方法是 * 从临时文件中获得(IRC=RCFC),或 * 在IRC计算的初始进行计算(IRC=CalcFC) IRC计算在输出文件末尾对计算进行总结,列出能量和优化的变量的值.第一个 值和最后一个值是整条路径的起点和终点. 在起点上,我们得到了一个类似甲醛分子的结构,可以认定该反应路线是通向 甲醛的,在终点上,得到了一个C-H键伸长的结构,C-O键略微缩短,也表明这条 反应路线是通向解离分子的. 计算活化能 IRC计算确认了所得到的就是我们所要的过渡态,下面就可以计算活化能了. 能量 活化能(kcal/mol) 过渡态 -113.67578 反应物 -113.83966 102.8(正向) 产物 -113.84996 109.3(反向) 计算表明两个反应方向的势垒相似. 注意IRC得到的产物的能量不一定等于两个单独的体系的和,因为当IRC计算得到 分子配合物的极小值,与两个分离体系的能量和有些差别 1,2氢迁移反应 现在用同样的步骤研究第二个反应. 反式氢基卡宾的包含零点能的总能量是-113.75709,计算方法是RHF/6-31G(d). 寻找过渡态 猜测过渡态在碳原子上的一个氢原子象氧原子方向迁移,处于同时与碳原子和 氧原子作用的位置,对其进行的频率分析表明其位一阶鞍点,包含零点能的总能 量为-113.67941 反应路径分析 IRC分析得到的两个结构,一个类似于HCOH,一个类似于H2CO,说明该结构为该反 应的过渡态. 活化能预测 计算得到的活化能为100.6(正向)和48.7(反向)kcal/mol IRC的注意事项 虽然实际的反应结构如极小点,极大点,鞍点等在势能面上存在几何的和数学的 意义,但不能简单推广到物理的和化学的意义.实际的分子是有动能的,这样它就 可以不遵循反应路径.当然,计算得到的结果提供了最经济的反应途径. 8.7 等构反应(Isodesmic Reactions) 等构反应是指反应前后各种键的数量不变的反应,比如乙醛与乙烷生成丙酮和甲 烷的反应,反应前后,各种价键的数量都没有变化. 由于这一特点,对这样体系的研究可以得到相当精确的结果. 例8.5 文件e8_05 等构反应的反应焓变化 本例计算上面提到的反应的反应焓变,步骤如下: * 在HF/6-31G(d)水平优化结构 * 进行频率分析计算零点能 * 在B3LYP/6-311+G(3df,2p)水平计算单点能 结果得到反应焓变-9.95kcal/mol,实验值-9.9+-0.3kcal/mol. 例8.6 文件e8_06 通过等构反应确定二氧化碳分子生成焓 本例中利用等构反应CO2 + CH4 --> 2H2CO 来确定二氧化碳分子的生成焓. 原理如下: dHcalc = 2 E0(H2CO) - (E0(CO2) + E0(CH4)) dHf(CO2) = -(dHcalc - dHf(CH4) + 2dHf(H2CO)) 甲烷和甲醛的生成焓实验值分别为-16.0和-25.0kcal/mol(0K) 计算方法同上例,所得到的反应焓变为60.64kcal/mol,这样计算得到的 二氧化碳的生成焓为-94.64kcal/mol.实验值为-93.96kcal/mol. 等构反应的局限 等构反应的局限 等构反应对于反应体系和生成焓的研究非常重要,但其也有明显的缺点,如下 * 必须依赖良好的实验数据,实验数据不准确,得到的生成焓自然也不可信 * 该技术不能推广到活化势垒方面 * 该技术不能用于实际上发生不了的等构反应 * 不同的等构反应可以得到不同的生成焓数值. 例8.7 文件e8_07 等构反应的局限 本例通过两个不同的等构反应计算乙烷的生成焓,研究理论 MP2/6-31G(d)//HF/6-31G(d) 反应一 丙烷 + 氢 --> 乙烷 + 甲烷 反应二 乙烷 + 氢 --> 2 甲烷 另外通过6个反应计算SiF4的生成焓,研究方法 MP2/6-31G(d,p)//HF/6-31G(d) SiH4 + F2 --> 2H2 + SiF4 (i) SiH4 + 4HF --> 4H2 + SiF4 (ii) SiF3H + HF --> SiF4 (iii) SiF2H2 + 2HF --> 2H2 + SiF4 (iv) SiFH3 + 3HF --> 3H2 + SiF4 (v) SiH4 + 4F2 --> SiF4 + 4HF (vi) 结果如下 乙烷生成焓实验值 -20.0+-0.1 kcal/mol 反应一 -17.361 反应二 -22.258 SiF4 实验值 -386.0+-0.3 kcal/mol (i) -375.983 (ii) -366.588 (iii) -380.506 (iv) N/A (v) -376.112 (vi) -385.379 虽然得到的乙烷生成焓和实验值比较吻合,但两个结果之间竟然有高达5kcal/mol 的差距,对于这样的简单的烃的体系仍然有这样大的差异,在使用这一方法是就 不得不小心了. 对于SiF4,不同反应得到的数据差别也很明显,虽然也有很符合的结果,但最大的 差距竟然达到20kcal/mol.同时注意这些硅化合物的生成焓本身的可靠性不高. 8.8 练习 练习8.1 文件 8_01a~c 水和反应 计算两个反应的反应焓 Li+ + H2O --> H2OLi+ at 298.15K H2O + H2O --> (H2O)2 at 373K 两个反应的实验值分别是-34.0+-0.2kcal/mol和-3.6+-.5kcal/mol 计算方法是B3LYP/6-311+G(2df,2p)//B3LYP/6-31G(d) 由于本例中要计算298.15K和373K两个温度的水分子的热力学数据,所以,可以在 一个计算中完成,输入文件如下 %Chk=water #T RHF/6-31G(d) Freq=ReadIso ... --Link1-- %Chk=water %Nosave #T RHF/6-31G(d) Geom=Check Freq=(ReacFC, ReadIso) Guess=Read Test ... 373.0 1.0 0.9135 16 16 1 1 计算得到的第一个反应反应焓变为-34.0kcal/mol.注意对Li+不能进行频率分析, 其焓的校正为1.5RT(只有平动能) 对于第二个反应,注意由于DFT方法在处理弱作用体系上有些困难,要确认所得到 的分子是真正的极小值,而不是过渡态.本例中需要采用Opt=CalcAll来进行优化. 计算结果为反应焓-2.9kcal/mol 练习8.2 文件8_02a, 8_02b 键的解离 本例中通过势能面扫描研究键的断裂过程,研究的体系是CH键 CH UMP4/6-311+G(d,p) 键长范围 0.2-2.5A CH4 RQCISD(T)/6-311++G(d,p) 和键长范围 .75-3.15A UQCISD(T,E4T)/6-311++G(d,p) UQCISD的E4T选项进行的是在MP4(SDTQ)水平进行MP4计算,而不是默认的MP4(SDQ). 计算中还需要设置的关键词有 IOP(2/16=1)表示在扫描中忽略对称性变化 Guess=(Always,Mix)表示将HOMO和LUMO轨道混合来消除自旋对称性的影响,并且 在每一个点都重新计算新的猜测波函数. 下面是检查势能面扫描数据的方法 * 将HF, MP2, MP4, QCISD(T)不同水平计算的键能对键长作图 * 利用这些图来分析键解离: * 限制性和非限制性方法的比较 * 电子相关 * 电子相关 * 三重态对QCISD水平的贡献 对于CH,UHF曲线在其他曲线之上,相对而言,HF的结果是最差的.MP2的曲线比 MP3和MP4(SDTQ)曲线要高些,但三者相差不多. 对于CH4,UHF得到的曲线远高于其他结果,MP2的结果也比其他结果要高些,而 其他各个计算方法得到的曲线基本相似.每一次提高微扰的次数,能量曲线就 向低的方向移动. QCISD曲线与MP4曲线很接近,在对于整个曲线的描述上,QCISD(T)比MP4(SDTQ)的 结果要好很多. 限制性方法和非限制性方法得到的结果只在键分裂之后才显现出来,RHF方法 得到的曲线趋势就是错误的,因为其得到的结果是,当原子进一步远离时,能量 继续上升,实际上当两个原子距离达到一定程度后,能量将基本保持不变.而对 于QCISD(T)方法,两者的差别不大,限制性方法得到的曲线在键解离后要略高于 非限制性方法的结果.这说明,在描述键解离的过程时,非限制性方法的结果要 比限制性方法好,特别是在采取低等级计算模型时. 练习8.3 文件 8_03 H2CO的势能面 我们以前研究过H2CO的两个反应的势能面,得到了五个稳定点:其中三个极小值, 甲醛,反式羟基卡宾以及一氧化碳和氢分子;两个过渡态分别连接甲醛和另外两个 产物.现在要做的,就是研究这两个产物之间的变化途径.反应分两步: 反式羟基卡宾 <--> 顺式羟基卡宾 <--> 一氧化碳 + 氢 在HF/6-31G(d)水平研究这一问题,研究步骤是: * 寻找过渡态 * 确定得到的稳定点是过渡态,计算零点能 * 确定这一过渡态连接的两个极小值 一个可能的连接顺式和反式结构的过渡态就是HCO平面与COH平面成90度的二面 角.我们采用Opt=QST3,分别给出两个异构体和这一猜测的构型.然后进行IRC计 算,得到的是接近顺式和反式异构体的构型,说明这个过渡态结构就是所要找的 过渡态. 用同样的方法,可以得到有顺式异构体到一氧化碳和氢的反应的过渡态. 反应的活化能为 正向 反向 反式 <--> 顺式 25.8 20.5 顺式 <--> 解离 68.0 131.6 这样就得到了关于H2CO的新的势能面. 练习8.4 文件 8_04 原子电荷分析 这个练习采用非Mulliken布局分析方法计算原子电荷,考察这些结果. 原子电荷计算并不是有量子化学理论推导出来的,所有的计算原子电荷的方法都 是武断的. 练习采用的方法有 * Muliken布局分析(默认) * Natural布局分析(Pop=NPA) * 采用CHelpG方法的Breneman静电势衍生电荷 * 采用Merz-Kollman-Singh方法的静电势衍生电荷 计算方法为MP2/6-31G(d).为了节省计算时间,在计算第二个以及后面的几个算 例时可以采用在临时文件中的电子密度,Geom=CHeckpoint Density=(Checkpoint,MP2),同时注意在第一个算例中加入Density=MP2从而采 用MP2方法得到的电子密度来计算原子电荷,默认方法为HF方法. 计算的分子是丙烯阳离子. Mulliken方法得到的原子电荷,三个碳原子上都有一定的负电荷; Natural方法的结果是,中间碳原子上是负电荷,其余原子上均为正电荷; CHelpG和MKS方法得到的结果,也只有中间的碳原子上有负电荷.后面三个方法 得到的中间碳原子上负电荷的值越来越少. 练习8.5 文件 8_04 基团电荷 计算上例中的CH2和CH基团的电荷,结果如下 Mulliken NPA CHelpG MKS CH +0.18 -0.05 -0.10 +0.09 CH2 +0.41 +0.52 +0.47 +0.45 练习8.6 文件 8_06 分子中的原子电荷和键级 分子中的原子理论(theory of atoms in molecules),提供了另外的更加复杂的 计算原子电荷和相关性质的方法. 这一理论利用图论和电子密度函数来定义化学性质,如键和原子电荷.这样的理论 有清楚的量子化学意义. Gassian中的AIM关键词涉及到这一算法. 这里用AIM=BondOrders来计算丙烯正离子的原子电荷和键级. 计算的碳碳键键级为1.4,碳氢键键级为0.9. 所有原子上都是正电荷. 练习8.7 文件 8_07 Si+ 和硅烷的势能面 硅簇反应是一个新的领域,非常值得用电子结构方法研究. 这个练习检验硅正离子进攻硅烷(SiH4)的势能面.这个反应是硅簇反应的中心: Si+ + SiH4 --> Si2H2+ --> Si3H4+ --> Si4H6+ --> Si5H10+ ... 每一步反应都生成氢 我们只研究第一个反应 H Si+ | | H Si-H + Si+ --> Si + | / \ / \ H H H H H 在Krishnan Raghavachari的研究中,发现了如下的极小值 H H H H H H2 \ \ \ / \ / | Si-H...Si Si-Si Si-Si Si--Si / | / | \ / \ / H H H H H H H Si SiH4...Si+ H3Si...SiH+ H2Si...SiH2+ H2Si-Si+...H2 最右边的结构是产物配合物,结合能1kcal/mol,在我们的研究中可以认为是反应 的终点. 另外还有如下过渡态 H H-H \ / \| Si--Si / H 确定哪一个是连接过渡态和的极小值就可以确认反应势垒. 计算采用HF/6-31G(d)进行频率分析和IRC,采用MP4方法进行能量计算, IRC计算中要包含下列关键词 IRC=(RCFC, StepSize=30, MaxPoints=15) SCF=QC 这样的设置增加了两端搜索的点的数量. 具体的计算方法是 * 进行频率分析,计算零点能,准备IRC * IRC计算,确定与过渡态相连接的极小值 * UMP4/6-31G(d,p)计算能量 最终确认该过渡态连接的是H3Si-SiH+和H2Si-Si+...H2 其他的极小值当然与其他的过渡态相连.有兴趣的可以继续寻找其他的过渡态, 完成整个反应的势能面. 练习8.8 文件 8_08 等构反应 这里研究等构反应 CH3COX + H3CCH3 <--> CH3COCH3 + CH3X X = H, F, Cl 计算反应焓变,方法为B3LYP/6-311+G(3df,2p)//B3LYP/6-31G(d) 结果为-9.95, 16.71, 7.14;实验值-9.9+-0.3, 17.9+-1.3, 6.6+-0.3 采用其他方法,HF方法的结果很糟糕,MP2方法过分考虑了电子相关,MP3 方法的结果也比较好.对于AM1方法,结果比HF方法还要差. 练习8.9 文件 8_09a~b 通过等构反应计算生成焓 通过等构反应研究氟代甲烷和苯的生成焓.所用反应为 CHF3 + 2CH4 --> 3CH3F C6H6 + 6CH4 --> 3C2H6 + 3C2H4 结果为-54.27和23.89,实验值-55.9+-2.0,24.0+-0.2 练习8.10文件 8_10 一个SN2反应 研究SN2反应Cl- + H3CF --> ClCH3 + F- 方法如下 * 过渡态结构优化 * 过渡态频率分析 * 从过渡态的IRC分析 * 寻找中间体两个极小值的几何结构 * 频率分析计算零点能 得到的过渡态是Cl...CH3...F 两个极小值是F-...CH3Cl和CH3F...Cl- |
最近读者：