1、1实验四 应用量子化学计算方法进行分子结构优化以及异构化反应研究Experiment 4. Study on Molecular Structure Optimization and Isomerization Reaction by Using Quantum Chemistry Method4.1 目的要求 Purpose（1）了解量子化学计算的原理和用途以及几种常用的量子化学计算方法。（2）熟悉常用量子化学计算软件 Gaussian 03 的基本使用方法和操作步骤。（3）掌握如何使用 Gaussian 03 软件进行分子结构优化和异构化反应过渡态计算。（4）本实验 4 学时。4.2 背景
2、介绍 Background Information量子化学（quantum chemistry）以量子力学为理论基础，以计算机为工具，主要通过计算来阐述物质（化合物、晶体、离子、过渡态、反应中间体等）的结构、性质、反应性能及反应机理，研究物质的微观结构与宏观性质的关系，揭示物质和化学反应所具有的特性的内在本质及其规律性 1-4。随着量子化学计算方法不断发展，计算量以及计算速度不断提高，所计算的体系越来越复杂，现在可以计算有机分子甚至较大分子量的生物分子。目前常用的量子化学计算软件有 Gaussian（）、GAMESS（www.msg.ameslab.gov/GAMESS）、Spartan（ ）
3、和Molpro（）等。Gaussian 软件是使用最为广泛的量子化学计算软件，支持几乎所有的量子化学计算方法，可以计算得到分子的几乎一切性质，如稳定结构、能量、振动频率、红外和拉曼光谱、NMR 化学位移、轨道能级、静电势、极化率、电离能、电子亲和力、电子密度分布、过渡态和反应途径等。可以模拟在气相和溶液中的体系，模拟基态和激发态等问题。它最早的版本是 1970 年的 Gaussian 70，最新的版本是 Gaussian 09。本实验使用的版本为 Gaussian 03。4.3 实验原理 Experimental Principles4.3.1 量子化学计算方法和特点多体理论是量子化学的核心问
4、题。n个粒子构成的量子体系的性质原则上可通过求解n粒子体系的薛定谔（Schrdinger）方程得到体系的波函数来描述。（式4-1 ）22,11pqppi jiP iZZEmRr 量子化学是通过求解Schrdinger方程（式4-1） ，得到分子中电子的运动状态。原则上，对Schrdinger方程的求解可以获得对分子这样的多电子体系中电子结构和相互作用的全部描述。但是由于该方程包含核与电子的双重坐标，并且二者很难进行变量分离，一般性地2精确求解Schrdinger方程是相当困难的。因此，在实践中总希望发展和运用量子力学的近似方法, 从而无需进行很繁杂的计算就可以说明复杂原子体系的主要特性, 这就
5、必须在原始量子化学方程中引进一些重要的简化, 以便得到一定程度的近似解。量子化学发展到现在, 根据为解薛定谔方程而引入近似程度的不同, 大致可分为以下几种方法:4.3.1.1 密度泛函理论（Density Functional Theory，DFT ）1964年，Kohn 提出了密度泛函理论 6-14，其以Hohenberg-Kohn定理为基础，指出电子密度决定分子的一切性质，体系的能量是电子密度的泛函。各种密度泛函理论差别在于选择交换泛函和相关泛函的不同，例如：纯密度泛函包含一个相关泛函和一个交换泛函，如BDW、 BLYP 等；而杂密度泛函则包含一个相关泛函和多个交换泛函，例如B3LYP、B
6、HandHLYP 等。由于密度泛函计算结果精确，计算速度快，DFT以无可比拟的优越性成为当前国际研究的主流方向，与分子动力学结合的分子模拟，更是当前理论化学研究反映动态过程的有利工具，成为计算材料科学的重要基础和核心技术。但DFT也并不是适合所有的体系，研究表明，其对共价键体系计算结果精确，氢键体系次之，Van der Waals键体系再次之。4.3.1.2 从头算方法（Ab Initio Method）在量子化学中，从头算法是指基于Born-Oppenheimer、独立电子（Hartree）和非相对论三大近似，利用电子质量、普朗克常数和电量三个基本物理量及原子系数，对分子的全部积分严格进行计
7、算，不借助任何经验参数来求解薛定谔方程的计算方法 5, 15-19。不同的从头算法考虑了不同的相关能项，如：HF 方法只考虑了同电子自旋的相关（交换相关）问题，而没有考虑相反自旋的电子相关问题和瞬时电子相关的问题；MPn 方法给体系考虑了微扰项，而更为精确的计算应包含更多的相关能相，如组态相互作用方法（CIS，CISD）和耦合簇方法（ CASSCF）等。由于从头算法在理论上的严格性和计算结果的精确性、可靠性，从小分子体系到大分子体系，从静态性质到动态性质，各方面都有从头算法的应用。对过渡金属配合物，金属原子簇合物等大分子化合物的研究也迅速增加。但基于计算精度和计算资源的矛盾的考虑，从头算法主要
8、应用于小分子体系的高精度计算、对中等大小分子体系进行定量计算、对大分子体系的定性计算三个方面。4.3.1.3 半经验方法（Semi-empirical Method ）从头算法虽然有严谨的理论支持，能得到较好的计算结果，但是当遇到诸如酶、聚合物、蛋白质等大分子体系时，计算很耗时，其计算代价无法承受。为了在计算时间和计算精度上找到一个平衡点，科学家们以从头算法为基础，忽略一些计算量极大，但是对结果影响极小的积分，或者引用一些来自实验的参数，从而近似求解薛定谔方程，就诞生了半经验算法。其最核心的近似方法是忽略双原子轨道微分重叠的NDDO近似（neglect of diatomic differen
9、tial overlap）。在此基础上又有了所谓全略微分重叠近似的CNDO（ complete neglect of differential overlap）、间略微分重叠近似的INDO（intermediate neglect of differential overlap）和改进的间略微分重叠近似 MINDO（modified intermediate neglect of differential overlap）等。以后的半经验方法，是在 MINDO方法的基础上进行改进。目前常用的改进后的方法有ZINDO（Zemer INDO）、AM1和PM3等 20, 21。半经验方法理论上没有从
10、头算法那么严谨，因而在处理复杂体系的中间体、过渡态时3会遇到一定的困难，其计算的结果只带有定性和半定量的特性。主要用于非常大的体系的计算或处理大体系的第一步，或为了得到一些分子的初步研究结果。4.3.2 常用基组量子化学中的基组是用于描述体系波函数的若干具有一定性质的函数，是从头计算的基础，有着非常重要的意义。在量子化学计算中，根据体系的不同，需要选择不同的基组，构成基组的函数越多，基组便越大，对计算的限制就越小，计算的精度也越高，同时计算量也会随基组的增大而剧增。Slater基组：斯莱特型基组就是原子轨道基组，由体系中各个原子中的原子轨道波函数组成。斯莱特型基组是最原始的基组，函数形式有明确
11、的物理意义，但数学性质并不好，随着量子化学理论的发展，Slater型基组很快就被淘汰了。Gauss基组：Gauss基组用Gauss函数替代了原来的Slater函数。Gauss 函数在计算中有较好的性质，可以将三中心和四中心的双电子积分轻易转化为二中心的双电子积分，因而可以在相当程度上简化计算，但也会使得量子化学计算的精度下降。最小基组：最小基组又叫STO-3G基组，STO是Slater型原子轨道的缩写， 3G表示每个Slater型原子轨道是由三个Gauss函数线性组合获得。STO-3G基组用三个Gauss函数的线性组合来描述一个原子轨道，对原子轨道列出HF方程进行自洽场计算，以获得Gauss函
12、数的指数和组合系数。STO-3G基组规模小，计算精度相对差，但是计算量最小，适合较大分子体系的计算。比STO-KG更好的基组是“分裂价基” （split valence basis sets） 。分裂价基是通过每个Slater原子内层轨道用一个STO-KG逼近，而每个价层轨道用两个新的具有不同指数STO-KG来逼近，从而改善计算结果与实验结果的偏差。例如 4-31G，其中每个内层轨道用一个STO-4G逼近而每个价层轨道用一个STO-3G和一个STO-1G 逼近。对于元素周期表中第三周期以下的元素（主族及过渡金属元素） ，随着d轨道的出现，其化学性质及结构特点都将发生变化，此时含有极化函数的Ga
13、uss基组对于计算结果处理变得更合理。分裂价基组允许轨道改变其大小，但不能改变形状。极化基组则取消了这样的限制，它是通过增加描述每个原子基态需要的轨道角动量来完成的。例如，极化基组在C 原子上增加d轨道的成分，在过渡金属上增加f轨道成分，在 H原子上增加p轨道成分。例如 6-31G(d,p)，它是在6-31G 基组的基础上，增加d轨道到重原子上，增加p轨道到H 原子上。对于电子相对离原子核比较远的体系，如含孤对电子的分子、阴离子以及其它带有明显负电荷的体系、处于激发态并且具有较低的离子化能的体系，弥散函数具有重要的作用。弥散函数允许轨道占据更大的空间。例如，6-31+G(d)基组表示的是 6-
14、31G(d)基组在重原子上加上弥散函数，而6-31+G(d)基组表示对于氢原子也加上弥散函数。4.3.3 振动频率的计算振动频率对于分子势能面的表征具有十分重要的意义。首先，通过区分全部为正频率的局域极小点（local minima）和存在一个虚频的鞍点（saddle point），可以确定分子势能面上的稳定点；其次，振动频率可以被用来确认稳定的但又有高的反应活性或者寿命短的分子；最后，计算得到的正则振动频率按统计力学方法给出稳定分子的热力学性质，例如被实验化学家们广泛使用的熵、焓、平衡态同位素效应以及零点振动能估测等。分子振动涉及到由化学键相连接的原子间相对位置的移动。由于分子内化学键的作用
15、，4使得各原子核处于能量最低的平衡构型并在其平衡位置附近振动。我们可以把振动与运动幅度相对较大的平动和转动分离开来，对于一个由N 个原子组成的分子，忽略势能高次项，在其平衡态附近原子核的振动总能量可近似表述为：（式4-2）22331,1()NNeqeqiji ijijVETV式中 ，M i为原子质量，x i,eq为核的平衡位置坐标，x i代表偏离平衡位置12,()iiieqqx的坐标，V eq为平衡位置的势能，可取为势能零点。按照Lagrange方程：，i=1, 2, 3, , 3N0idTVtq代入T和V的表示式则有微分方程：，j=1, 2, 3, , 3N （式4-3）3311Nijfq式
16、中 为力常数矩阵F的矩阵元，f ij 可由势能一阶导数的数值微商或解析的2()ijeqijVf二次微商得到。最后可得到久期方程：（式4-4）31()0NijijjfC其中 ij = 1（i = j时）或 ij = 0（i j 时）。当久期行列式|FI|0时C j才有非零解，式中I为单位矩阵。解此本征方程可求出本征值和相应的本征矢量。各原子以相同的频率和初相位绕其平衡位置作简谐振动并同时通过其平衡位置，这种振动叫做正则振动。式4-3利用标准方法求得3N个正则模式下的频率模式，其中 6个（对于线性分子为5个）频率值趋于零，其物理意义是扣除了平动和转动自由度。4.3.4 过渡态理论过渡态理论是193
17、5年由Eyring和Polany等人在统计热力学和量子力学的基础上提出来的。他们认为由反应物分子变成生成物分子,中间一定要经过一个过渡态,而形成这个过渡态必须吸取一定的活化能，这个过渡态就称为活化络合物，所以又称为活化络合物理论。用该理论，只要知道分子的振动频率、质量、核间距等基本物性，就能计算反应的速率系数，所以又称为绝对反应速率理论（absolute rate theory） 。该理论认为反应物分子间相互作用的势能是分子间相对位置的函数，活化络合物与反应物之间建立化学平衡，总反应速率由活化络合物转化成产物的速率决定。传统过渡态理论建立在四条基本假定的基础之上：(1) 假定反应的分子体系是核
18、运动绝热的，即分子体系的核间相对运动可以和其它运5动（主要是电子运动）分离，这即是Born-Oppenheimer 近似。(2) 假定反应体系服从Boltzmann统计平衡分布规律。(3) 不返回假定。认为穿过分割面的轨线由反应物区进入产物区以后不再返回到反应物区。因此穿过分割面的轨线都是反应性轨线。(4) 运动分离假定。认为反应分子体系沿反应坐标变化时其变化速度很小，可作经典处理，从而可把沿反应坐标方向的运动与其它运动坐标分离，把其它所有运动自由度都认为是绝热的。根据传统过渡态理论的基本假定可以推导出反应速率表达式： ()()exp()RQTKUsh其中， 表示从反应物到产物等价的反应通道数
19、； ，k B是Boltzmann 常数；Q 1(T) 表示过渡态的单位体积配分函数；Q R (T)表示反应物单位体积配分函数； U(s )是鞍点的势垒高度；s 是过渡态点(s=0)。由此可见，传统过渡态理论是把分割面选在 s=0的过渡态处，用统计力学方法计算反应速率常数。由于这种理论基于轨线不返回假定，其结果必然会导致过高地估计反应速率常数，即传统过渡态理论提供了速率常数的经典上限值。4.4 实验操作 Experimental Procedures4.4.1 Gaussian 03 的基本操作4.4.1.1 启动程序点击 g03w.exe 执行文件，屏幕上显示 G03W 的主窗口。工具条上各按
20、钮的功能如下：开始作业暂停作业当前作业正在执行的 Link 完成后暂时停止恢复暂停作业终止作业编辑批处理文件（或建立新文件）当前作业完成后终止当前的批处理终止作业和批处理任务打开 g03w.ini 指定的外部编辑器（ external editor） 浏览或编辑当前作业的输出文件4.4.1.2 准备计算作业的输入文件点击G03W “File ”下拉菜单中“New”选项，屏幕上跳出输入文件编辑对话框。框内工具条上各按钮的功能为：启动作业返回主窗口将对话框编辑内容存成输入文件（*.gjf）6放弃对话框编辑内容终止正在进行的作业运行作业通道校核工具（Check route utility）Gauss
21、ian 作业的输入文件名的后缀为*.gjf。它由6段构成：Link 0命令段（section）：检查点文件（chk *.chk，此文件将记录过程信息和结果信息）及其它临时文件（rwf*.rwf 、int*.int 、d2e*. d2e）的名称及位置。可包含多个输入行，每行的行首均用%号开始。作业控制段（route section）：规定使用的计算方法、基组以及所需的计算项目。可包含多个输入行，每行的行首均用#号开始。标题行（title section）：作业计算内容与目的的简要描述，以便于输出文件的阅读。电荷与自旋多重度（charge & multipl.）：给定分子所带的净电荷数（零或正负整
22、数）及自旋多重度（2s1，s为0或1/2 ）。如中性分子输入“0 1”，带正一价阳离子分子输入“1 2”。分子坐标输入部分（molecule specification）：给定分子中各原子的坐标。笛卡儿坐标或分子内坐标均可使用。可选的附加部分：根据某些作业的特殊需要而输入。下面以H 2O分子（见图4-1）在（从头算）HF/STO-3G水平上的几何构型优化计算为例说明各输入段的内容与格式。chk=h2o.chk 检查点文件名为H2O.chk# HF/STO-3G opt=z-matrix 指定计算方法和基组，opt=z-matrix为内坐标优化Water STO-3G structure 关于计
23、算作业的说明0 1 闭壳层中性分子，自旋多重度为1OH 1 r1H 1 r1 2 a1 分子内坐标，r1为键长，a1为键角r1=0.96a1=104.5Figure 4-1. Structure of H2O4.4.1.3 运行作业将所有内容输入后点击 按钮，屏幕上跳出保存输出文件的窗口。程序默认的输出文件与输入文件同名（后缀改为out ）。直接点击“保存”，作业立即启动。7当G03W主窗口“Run Progress”显示“Processing Complete.”时，作业正常结束。点击主窗口右上角的 按钮或者使用记事本查看结果文件*.out；也可使用绘图软件GaussView查看结果文件。G
24、aussView是一个专门设计的与Gaussian配套使用的软件，其主要用途有两个：构建Gaussian软件的输入文件；以三维结构图的形式显示 Gaussian计算的结果。4.4.2 HF/6-31G(d)水平下 CO2 和 CH4 的结构优化和频率分析构建输入文件 CO2.gjf 和 CH4.gjf，按照 4.4.1 介绍的基本操作方法运行作业，得到结果文件 CO2.out 和 CH4.out。使用记事本以及 GaussView 软件打开 CO2.out 和 CH4.out 文件，查看三维结构、几何参数、能量以及频率信息。4.4.2.1 CO2 (D2h)分子（见图 4-2）的结构优化和频率
25、分析构建输入文件 CO2.gjf：chk=co2.chk # HF/6-31G(d) opt=z-matrix freq freq为振动频率计算OPT and FREQ on CO2 0 1 Figure 4-2. Structure of CO2C O 1 rX 1 1.0 2 90.0 键角要求 0.0a180.0，所以线性分子要用虚原子 XO 1 r 3 a 2 180.0 r为键长，a为键角r=1.13a=90. 4.4.2.2 CH4 (Td) 分子（见图 4-3）的结构优化和频率分析构建输入文件CH4.gjf：chk= ch4.chk # HF/6-31G(d) opt=z-mat
26、rix freq OPT and FREQ on CH4 0 1 C Figure 4-3. Structure of CH4H 1 r H 1 r 2 a H 1 r 2 a 3 bH 1 r 2 a 3 -b r为键长，a为键角， b为二面角8r=1.09a=109.47122 b=120.4.4.3 甲醛的稳定构象以及异构化反应的理论研究Figure 4-4. Isomerization pathway of formaldehyde在 B3LYP/6-31G(d)水平下对甲醛的酮式和烯醇式两种异构体进行几何构型优化，找出能量最低构象；计算两种异构体的异构化反应过渡态，并进行内禀反应坐标
27、（IRC）计算验证过渡态连接的反应物和产物分子（见图 4-4）。4.4.3.1 HCHO 的结构优化和频率分析构建输入文件 hcho.gjf，运行作业，得到结果文件 hcho.out。使用记事本以及GaussView 软件 打开 hcho.out，查看三维结构、几何参数、能量以及频率信息。输入文件 hcho.gjf：%chk=hcho.chk # B3LYP/6-31G(d) opt=z-matrix freq OPT and FREQ on HCHO 0 1 CO 1 r1H 1 r2 2 a1H 1 r2 2 a1 3 b1r1 1.37r2 1.07a1 120.b1 180.4.4.3
28、.2 HCOH 的结构优化和频率分析9构建输入文件 hcoh.gjf，运行作业，得到结果文件 hcoh.out。使用记事本以及GaussView 软件 打开 hcoh.out，查看三维结构、几何参数、能量以及频率信息。输入文件 hcoh.gjf：%chk=hcoh.chk # B3LYP/6-31G(d) opt=z-matrix freq OPT and FREQ on HCOH 0 1 CO 1 r1H 2 r2 1 a1H 1 r3 2 a2 3 b1r1 1.37r2 1.07r3 1.10a1 120.a2 170.b1 180.4.4.3.3 HCHO 和 HCOH 的异构化反应过
29、渡态计算构建输入文件 ts.gjf，运行作业，得到结果文件 ts.out。使用记事本以及 GaussView 软件打开 ts.out，查看三维结构、几何参数、能量，检查是否有一振动模式为虚频（频率为负值）以及虚频模式对应的反应坐标方向。 输入文件 ts.gjf：%chk=ts.chk # B3LYP/6-31G(d) opt=(ts, calcfc, maxcyc=100) opt=ts 为过渡态计算，calcfc 第一点计算力Freq iop(1/11=1) 常数，maxcyc=100 最大迭代次数为 100iop(1/11=1)避免因负本征值不对而引起的TS on isomerizatio
30、n of formaldehyde 计算终止问题0 1 CO 1 r1H 2 r2 1 a1H 1 r3 2 a2 3 b110r1 1.37r2 1.40r3 1.06a1 70.a2 160.b1 180.4.4.3.4 IRC 计算确定过渡态连接何种反应物和产物拷贝做过渡态计算的 ts.chk 文件两份，存成 irc1.chk 和 irc2.chk，用于 IRC 计算。分别构建正向（Forward ）IRC 计算输入文件 irc1.gjf 和反向（ Reverse）IRC 计算输入文件irc2.gjf，运行作业，得到结果文件 irc1.out 和 irc2.out。输入文件 irc1.g
31、jf：%chk=irc1.chk # B3LYP/6-31G(d) IRC= (Forward, calcfc, Stepsize 为步长，Maxpoint 为计算的点数，Stepsize=5, Maxpoints=30) Geom=check Geom=check 和 Guess=check 分别为读取Guess=check irc1.chk 文件中的过渡态坐标和分子轨道Forward IRC on isomerization 0 1 输入文件 irc2.gjf：%chk=irc2.chk # B3LYP/6-31G(d) IRC=( Reverse, calcfc, Stepsize=5,
32、 Maxpoints=30) Geom=check Guess=check Reverse IRC on isomerization 0 1 分别读取结果文件 irc1.out 和 irc2.out 中最后一个点（第 30 个点）的坐标，用GaussView 软件打开，验证是否为反应物（HCHO）和产物（HCOH ）结构。4.5 实验报告要求 Important Information About The Report按附录 3 格式，在摘要部分要求用文字简单概括实验过程及结果；在实验目的部分简单说明实验的背景、实验方案及意义；实验过程中重要实验参数只需列出计算项目所用到的关键词及其定义；实验结果要求记录各分子的优化的三维结构图和几何参数，包括键长（单位：） 、键角和二面角（单位：度） ，键长保留 3 位小数，键角和二面角保留 1 位小数；经过零点能校正的分子的总能量（单位：hartree） ，保留 5 位小数；正向 IRC 计算和反
