量子化学计算方法.ppt

1、量子化学计算方法,章永凡 福州大学化学系 2003年3月,课程主要内容,一、简介； 二、Gaussian 98程序的使用及其结果处理； 三、其它相关软件的使用； 主要参考资料: G98用户手册或G98的帮助文件 相应网址: ,简 介,量子化学软件目的在于将量子化学复杂计算过程程序化，从而便于人们的使用、提高计算效率并具有较强的普适性。 绝多数量子化学程序是采用Fortran语言编写的(Fortran 77或Fortran 90)，通常有几千至上万行语句。,软件分类,计算原理,基于从头算方法(ab initio) Gaussian、ADF、Dalton、Gamess、Crystal等,基于半经验

2、或分子力学方法 MOPAC、EHMO、NNEW3等,研究对象,有限尺度体系(分子、簇合物等) Gaussian、ADF、Dalton、Gamess、MOPAC、EHMO等,无限周期重复体系(晶体、固体表面等) Crystal、ADF、NNEW3等,本研究室目前常用的量化软件： Gaussian 98: 由Pople等人编写，经过几十年的发展和完善，现 该软件已成为国际上公认的、计算结果具有较高 可靠性的量子化学软件，它包含从头算、半经验 以及分子力学等多种方法，可适用于不同尺度的 有限体系，除了部分稀土和放射性元素外，它可 处理周期表中其它元素形成的各种化合物； ADF 2000： 该软件主要

3、集中于采用各种密度泛函方法来研究 化合物的电子结构，它尤其适合于研究含有过渡 金属的化合物，例如金属团簇分子等； Crystal 98： 该软件由意大利都灵大学理论化学研究所开发， 采用从头算方法来研究固体及表面的电子结构； Nnew3： 由Hoffmann小组编写的基于EHT方法的量化软件， 主要用于处理固体及表面电子结构。,要解决的问题,当前的研究状况，包括计算水平、已解决和尚未解决的问题,计算过程,化合物构型的确定，具体途径包括：实验测定数据、利用软件进行构造等,根据现有的计算条件、模型的大 小以及所要解决的问题，选择可 行的计算方法,对计算结果进行加工和提取有用 的信息，一般包括能量、

4、轨道组 成、电荷、布居和成键分析等,计算模型和方法的选取是保证计算结果可靠性的关键，最 理想的情况是：1.所选取的计算模型与实际情形一致；2.采用 高级的计算方法。但是，由于受到计算软硬件的限制，在多数 情况下，很难同时做到上述两点要求，实际操作中，当计算模 型较大时，只能选择精确度较低的计算方法，只有对较小的模 型才能选取高级的计算方法。 因此，当确定了一种计算模型和方法后，需对其进行验证， 以保证计算结果的可靠性。假使当前的研究对象是化合物A， 可通过下列途径进行验证： 1. 与A化合物现有实验结果之间的比较； 2. 若无实验方面的报道，可对与A类似的化合物B进行研究，此 时以B的实验结果

5、作为参照； 3. 当上述方法行不通时，可以采用较大模型和较为高级的计算 方法得到的计算结果作为参照，该方法主要用于系列化合物 的研究：如对A1, A2, A3，先用大模型和基组对A1进行研究， 然后以该结果为参照，确定计算量适中的模型和方法并应用 于A2,A3。,Gaussian98程序的使用,G98的安装和运行； G98的功能和程序结构； 输入文件的编写与主要功能的使用； 补充说明；,G98程序的安装和运行,1. G98程序的安装： (1).确定运行平台：Windows或Linux？ (2).对Windows平台： 直接运行setup.exe，其余步骤按提示操作即可；也可将其它机 器上将已安

6、装好的G98直接拷贝到本机，但需设置运行环境。 对Linux平台： a.若G98是经过压缩过的(文件结尾为gz)，用gunzip命令解压： 例如： gunzip g98.linux.tar.gz b.若G98是打包的(文件结尾为tar)，用tar命令将其释放： 例如： tar xvf g98.linux.tar c.设置环境变量，以c shell为例，在用户根目录下的.cshrc文件 添加下列内容：(也可在执行g98前逐条运行) setenv g98root /home/test/src(设置g98所在目录，根据实际情况修改) source $g98root/g98/bsd/g98.login

7、(激活g98运行时所需环境变量) setenv GAUSS_EXEDIR $g98root/g98(设置g98可执行文件所处目录),注：对Linux平台，运行g98时，需注意权限问题，可用chmod 命令更改权限，将所安装的g98对所有用户开放。,2.G98程序的运行： (1).对Windows平台： a.对于刚安装好的g98，先检查环境设置情况：,左侧至上而下依次为：默认的文本编辑器；g98可执行文件所在目录； 计算中间结果存放目录；缺省的计算结果存储目录；缺省的输入文件 所在目录；PDB分子构型浏览器； 右侧至上而下依次为：设置显示属性(如背景色等)；设置文本编辑器 属性；计算过程控制属性

8、(尤其是批作业过程)；Default.Rou文件的编 辑(该文件内容为默认情况下，计算所花费的内存及硬盘大小),需设置正确， 否则运行将出错！,b.编写或打开g98输入文件,点击RUN，并给定输出 文件名后开始运行,c.g98运行过程的控制：,最上行按钮的功能从左至右依次为： 开始运行g98；暂停进程；运行至下一模块(link)时暂停进程； 重新启动进程；清除进程(停止运算)；编辑批作业；运行完 当前任务后，暂停批作业；停止批作业的运算；观看计算结 果；打开文本编辑器；,不要随意点击！,(2). Linux平台： a.编写输入文件： 用vi命令编写或在Windows下编写完毕后ftp至Linu

9、x系统； b.运行g98： g98输出文件名 B3LYP=-1.1647796; MP2=-1.129582 MP3=-1.1339601; MP4D=-1.1355271; MP4DQ=-1.13547 MP4SDQ=-1.13547; MP4SDTQ=-1.13547; CCSD=-1.1361895 CISD=-1.1361895,(2).基组的选择： 1). 全电子基组： 关键词：sto-3g, 3-21g, 4-31g, 6-21g, 6-31g, 6-311g, d95/d95v 说明：I). 不同的基组适用范围是不同的： STO-3G(H-Xe)；3-21G(H-Xe)；6-21

10、G(H-Cl) 4-31G(H-Ne)；6-31G(H-Kr)；6-311G(H-Kr) D95(H-Cl 除了Na, Mg)；D95V(H-Ne),说明：ii)基组的大小决定了基函数的数目，即体系的原子轨道 数目，因此可从所选择的基组来推断MO数目： sto3g：为最小基组，每个原子轨道用三个高斯函数 (GF)来描述，原子轨道数即为基函数数目。 如O：1s2s2p，原子轨道数为1+1+3=5 GF数目为3*5=15 321g：为劈裂(split)基组，其含义是：内层的每个AO 用3个GF描述，价层的AO劈裂为两组，分别用 2个和1个GF描述。显然，321g的GF数与sto- 3g是相同的。

11、如O：内层为1s，AO数为1，GF数为3 价层2s的AO数为2*1=2，GF数为2+1=3 价层2p的AO数为2*3=6，GF数为3*2+3*1=9 共1+2+6=9个AO和3+3+9=15个GF 对Mg：1s 2s 2p 3s 3p 内层1s, 2s和2p共有1+1+3=5个AO和3*5=15 个GF，价层3s有2个AO和3个GF，价层3p有 6个AO和9个GF，故共5+2+6=13个AO和 15+3+9=27个GF,说明：631g：为劈裂(split)基组，其含义与3-21g类似，内层的 每个AO用6个GF描述，价层的AO劈裂为两组，分 别用3个和1个GF描述。 如O：内层为1s，AO数为

12、1，GF数为6 价层2s的AO数为2*1=2，GF数为3+1=4 价层2p的AO数为2*3=6，GF数为3*3+3*1=12 共1+2+6=9个AO和6+4+12=22个GF 对Mg：1s 2s 2p 3s 3p 内层1s, 2s和2p共有1+1+3=5个AO和6*5=30 个GF，价层3s有2个AO和4个GF，价层3p有 6个AO和12个GF，故共5+2+6=13个AO和 30+4+12=46个GF 6311g：也为劈裂基组，自是价层的AO劈裂为3组，分别 用3个、1个和1个GF描述。 对于4-31g和6-21g类似。,Iii)极化(polarization)函数的使用： 在实际计算中，有时

13、需在上述标准基组的基础上，添加一个 或多个极化函数，极化函数是指具有比原子价轨道更高角量子 数的高斯函数。例如H的价轨道为1s，则其极化函数为p型GF， 同样对C、O等价层为p轨道的原子，它们的极化函数应为d型 或f型轨道，类似地，对于过渡金属原子的极化函数为f型轨道。 极化函数的使用目的在于使原子价轨道在空间取向上变得 更“柔软”，从而使之易于与其它原子的轨道成键：,例如对于羰基基团中的C，O 原子，它们极化轨道(即d轨道 )之间形成的d轨道，使得C 原子的p轨道朝O方向极化， O的p轨道向C方向极化，从 而增强了C与O之间的作用。,由于极化轨道的使用增强价轨道的空间柔软程度，因此其使用 场

14、合主要在于环状化合物，例如有机环状化合物以及含有桥联 结构的金属簇合物等。 使用方法： 在标准基组后加“*”，“*”，(d)，(d,p)，(3df,3pd)等 例：6-31G*, 6-31G(df,pd) 说明： 当只添加一个极化函数时，可使用下述表示： “*”等价于(d)，含义是对非氢原子添加一个极化轨道； “*”等价于(d,p)，含义是对非氢原子添加一个极化轨道的同时 对H原子添加一个p极化函数； 当添加多个极化函数时，使用下述表示： (?d,?p)，表示对非H原子和H原子分别添加？个极化轨道 如：6-31G(3d,2p),(?df,?pd), 表示对非H原子添加?个极化函数和1个具有更高

15、角 动量的极化函数; 对H原子分别添加?个p极化函数和 1个d极化函数; 例如, 对H2O计算时采用6-31G(2df,3pd)基组,其含义是: 对其中的H原子,在6-31G基组上在添加3个p轨道和 1个d轨道; 对O原子, 则在6-31G基组上添加2个d轨道和 1个f轨道. 对于不同的基组,可添加的极化函数数目是不同的, 具体可参照 程序手册:,iv) 弥散函数(diffuse)的使用: 对于带有较多电荷的体系,采用标准的基组来描述是不够的, 此时需添加弥散函数,以增加价轨道在空间上的分布范围, 即: 极化函数用于改进价轨道的角度分布, 弥散函数则用于改 进价轨道的径向分布. 弥散函数是指具

16、有较小轨道指数的高斯函数, 其表示方法是 在标准基组后加上“+”或“+”, 如6-31+G, 6-31+G*等, 其中 第一个“+”表示对非H原子添加弥散函数, 第二个“+”则对H原 子添加弥散函数. 弥散函数主要用于带有电荷的体系(包括离子)以及弱作用体系.,2) 赝势基组： 关键词：lanl1mb,lanl2mb,lanl1dz,lanl2dz, cep-4g,cep-31g,cep-121g,sdd等 说明： 通常对于重元素采使用赝势基组，其中，对于过渡金属原子 一般使用lanl系列的基组，对于主族元素采用cep系列的基组， sdd基组相对较少使用； 赝势基组的适用范围要较全电子基组来得

17、广，一般除了稀土 和放射性元素外，均可使用赝势基组； 赝势基组是将内层原子的影响用势函数来代替，因此，在使 用时必须注意根据所采用的赝势基组，检查体系的总电子数 是否正确！ 对lanl1系列的基组，只考虑价层电子；如V原子为5. 对lanl2系列的基组，除了价层电子外，还考虑了次外 层的电子(对部分主族元素例外)；如V原子为5+8=13. 对其它赝势基组，请在使用时仔细核查。,对于使用较多的lanl系列基组，需注意： 对于HNe范围的原子，lanl1mb和lan2mb基组等价于 sto-3g最小基组；lanl1dz和lanl2dz基组等价于d95基组; 对于含金属金属键体系一般使用lanl2系

18、列，而不用 lanl1系列； mb结尾含义为minimal basis，dz结尾含义为double -zeta，故后者基函数数目要较前者多。,3) 混合基组的使用: 有时, 因需要, 对同一化合物中的不同原子采用不同的基组: 受基组限制: 例如, 对于Mo(CO)2+化合物, C和O原子可采用6-31+G基组, 但 Mo原子不在6-31G基组的使用范围之内. 受体系大小限制: 当体系较大时, 不可能对所有原子均采用较大的基组, 此时可 对其中的局部原子采用精度较高的基组来描述,而对其余原子 采用小基组来描述, 例如:,混合基组输入方法: 关键词: GEN和PSEUDO=READ(对赝势基组)

19、%mem=32mb %chk=hf #p b3lyp/ gen pop=full scfcyc=500 Mo(CO)2+(Line) 1,6 Mo X11.0 C1r1290.0 O1r2290.030.0 C1r1290.03180.0 O1r2290.03180.0 r1=2.228 r2=3.329 Mo 0 3-21G * C O 0 6-31+G* *,整数0为结束标志,4个*用于分隔不同基组,改行也可用数字代替,即2 3 4 5 0 (编号时不包含虚原子),%mem=32mb %chk=hf #p b3lyp/ gen pseudo=read pop=full scfcyc=500

20、 Mo(CO)2+(Line) 1,6 Mo X11.0 C1r1290.0 O1r2290.030.0 C1r1290.03180.0 O1r2290.03180.0 r1=2.228 r2=3.329 Mo 0 Lanl2mb * C O 0 6-31G* * Mo 0 Lanl2mb,使用赝势基组时,应注意不能漏!,基组定义部分,赝势定义部分,注意不能漏, 另外其后 没有*号,赝势部分和基组部分之间有一空行,(3) G98输出结果: 下面以H2O能量计算的输出结果为例:,L1输出G98版权信息,%部分内容,关键词部分,计算所需调用的模块,Title内容,电荷, S以及分子构型,各原子的直

21、角坐标,距离矩阵,需核查!,注意:计算结果是以该坐标系为准,需核查!,MO初始猜测,自洽场迭代求解部分,收敛指标,迭代次数,(通常接近2.0),(布居分析结果),MO对称性及能级(a.u.),Mulliken键级和电荷,偶极矩及多极矩,对主要计算结果进行总结,(4) 能量计算中使用到的其它常用关键词: SCF收敛指标的设置: 关键词: scf(conver=?) 含义: 当SCF前后电荷密度的差值小于10-?时, SCF收敛. 例: scf(conver=5) 表示收敛指标为10-5. 说明: 对于能量计算,缺省情况下收敛指标为10-4. 对于构型优化, 则为10-8. SCF迭代轮数的设置:

22、 关键词: scf(maxcyc=?) 含义: SCF最大的允许迭代次数为?. 例:scf(maxcyc=100) 表示最大迭代次数为100次. 说明: 缺省情况下, SCF迭代次数为64次, 通常对于含有过渡 金属的体系, SCF较难收敛, 迭代次数会多些.,分子轨道系数的输出: 关键词: pop=full 说明: 有时为了分析MO成分, 则需利用该关键词输出各MO的 成分, 以H2O为例:,轨道对称性,O表示占据轨道; V为空轨道,轨道能级,AO轨道,轨道系数,HOMO,LUMO,注: 输出的MO系数是未经归一化过的; 1s, 2s等符号并非真正意义上的AO, 它与计算所采用的基组 有关.

23、 对于上例, 采用的是6-31G基组, 可知, 对于价层AO是 分裂为两组的, 故对于O原子, 输出的2s和3s实际上均为价 轨道, 其它类推. 当基函数较多时, 若只考察前线附近的轨道成分, 可用关键 词: pop=regular, 此时只给出5个占据轨道和5个空轨道组成. 练习: 计算O2的基态能量, 并分析MO轨道组成,(5) 能量计算的几个问题： 基组的选取： 理论上而言，计算中所选取的基组约大，计算结果越准确，但由 于受到硬件上的限制，需根据实际情况选择基组。此外，对于环状 分子或存在作用的体系，通常需要考虑使用极化基组；对于带有 较多电荷的体系，或考察弱相互作用的体系，通常需考虑使

24、用弥 散基组； 以O2为例(O-O为1.208A)： 6-31G-150.26755a.u. 6-31G*-150.31998a.u.dE=0.05243a.u. 6-31+G-150.27678a.u.dE=0.00923a.u. 由于O2中存在作用，故极化函数的影响要比弥散函数来得大。 以O22-为例： 6-31G-149.78721a.u. 6-31G*-149.82396a.u.dE=0.03675a.u. 6-31+G-149.98117a.u.dE=0.19396a.u. 6-31+G*-150.02156a.u.dE=0.23435a.u. 此时，因带较多负电荷，故弥散函数的影响

25、较极化函数来得大。,体系多重度的选取： 对于复杂体系，特别是含有过渡金属原子的体系，多重度的选取 评经验不易得到，此时，可通过比较不同多重度时体系的能量来 确定基态能量。 以O2为例(基组为6-31G*)： S=1-150.25734a.u. S=3-150.31998a.u. S=5-149.79200a.u. S=7-148.80784a.u. 由上结果可见，O2的基态应为三重态。 键能的计算： 以O2为例(O-O键长为1.208A，计算OO键能： 步骤如下： 首先确定O2基态的能量： 假设基组为6-31G*，则总能量为： -150.31998a.u.(基态为：3g),其次，计算O原子基态

26、能量，结果为：-75.06061a.u. 则O2中OO键的键能为： 2-75.06061+150.31998 =0.19876a.u.=521.8KJ/mol 该结果略高于实验值（494KJ/mol）。 在有关键能的计算中，除了计算方法以及基组对最终的结果 有较大影响外，若要得到精确的键能，还需考虑基组叠加误差 即BSSE。 BSSE是指源于能量计算中基组带来的误差，它可以粗略理解 为化学实验中的背景误差，要消除BSSE，可采用massage关键词 来实现，以O2为例具体使用过程如下：,%mem=32mb #p b3lyp/6-31g* massage BSSE of O2 0,3 O 0.

27、0.0 0.0 O 1.208 0.0 0.0 2 Nuc 0.0,将其中一个O原子的核电荷定义为零，但其 基组在计算过程中仍保留，该原子可称之为 鬼原子(ghost)，注意：鬼原子与前面所提到 的虚原子(dummy)是不同的，后者不存在基 组的问题；,基函数数目与计算O2时相同，但电子数目只为8个,由此得到考虑了BSSE时，O原子的能量为-75.06282a.u.，根据 该数值，可得OO键能为：510.4KJ/mol。此外，进一步的修 正是在计算O2总能量时，需考虑零点能(ZPE)校正。,将上述思想推广到一般情况，若要考察某个体系中A与B之间 的结合能力：,A,B,A,B,A,B,首先将A和

28、B一起考虑，得到体系 总能量E(AB),保留A部分，将B部分的原子均设 为鬼原子，得到考虑B部分基组效 应的A部分能量E(A),保留B部分，将A部分的原子均设 为鬼原子，得到考虑A部分基组效 应的B部分能量E(B),最终，A与B之间的结合能为：E(AB)-E(A)+E(B),分子轨道成分的分析： 以O2为例： %mem=32mb #p ROB3LYP/STO-3G POP=FULL The single point calc. Of O2 0,3 O 0. 0.0 0.0 O 1.208 0.0 0.0,为了便于计算结果分析，采用限制性 开壳层计算方法，及将和电子强制 配对；,输出所有分子轨道

29、,首先确定计算结果中各原子所处的位置：,2,1,o,z,X or y,z,+,1,2,因1s为内层电子，故O与O之间成键情况 可不考虑,No.1,No.2 与上一MO类似，也以1s成分为主，处在能量较低的位置,z,+,+,z,+,1,2,No.3, 4分别为O2s轨道之间形成的成键和反键轨道：,+,z,+,-,No.5为O 2pz轨道形成的成键轨道：,z,+,z,No.6和7轨道分别为O 2px和2py形成的成键轨道：,+,No.8和9轨道分别为O 2px和2py形成的反键轨道：,z,+,No.10轨道为O 2pz形成的反键轨道：,z,+,根据上述分析结果以及各轨道的能级，可以得到O2的MO

30、分布为：,2s,*2s,2pz,2px,2py,*2px,*2py,*2pz,电荷分布： 在能量计算结果中，除了能量以及上述的MO成分外，其它有用 的信息包括原子电荷和键级，对于原子电荷： 以H2O为例(B3LYP/STO-3G)：,Mulliken电荷分布,说明： 电荷的绝对值是没有意义的，其数值受到所用方法, 尤其是 基组的影响较大： 以H2O为例： 方法/基组 O H B3LYP/STO-3G-0.3290.165 B3LYP/3-21G-0.6370.318 B3LYP/6-31G-0.7050.353 HF/6-31G-0.7950.398 因此，计算得到的电荷只具有相对意义，即只能

31、考察其相对 数值或变化趋势。当然，当考察变化趋势时，如考察H2X (X=O,S,Se)系列化合物中，X为不同原子时，对H原子电荷的影 响时，必须使用相同的计算方法和基组。,键级： 键级的大小反映了某根键的强度，同样地，键级的数值也受到 计算方法和基组的影响：,O-H键的Mulliken键级,以H2O为例： 方法/基组 O-H键级 B3LYP/STO-3G 0.243 B3LYP/3-21G 0.239 B3LYP/6-31G 0.231 B3LYP/6-31+G 0.214 HF/6-31G 0.249 由上结果可见，键级数值与所用的基组和计算方法是有联系 的，但其影响不象电荷那么显著。此外，

32、必须注意到，弥散 基组对键级的数值影响较大。 因此，在实际应用中，对于电荷和键级，我们更关心的是 它们的变化情况而不是它们的绝对数值。此外，必须强调的 是，上述结果是采用Mulliken布居分析的方法，若采用其它 布居分析方法，结果将有所不同。,自然键轨道(NBO)分析方法： 我们前面所讲分子轨道均为正则分子轨道，它是未经定域化处 理，这导致计算结果与我们通常的成键概念有所不同。例如 在乙烯分子中，碳碳之间为双键，但在正则MO中，反映C与C 之间成键作用的MO可能有多个，因此根据正则MO的结果，我 们无法断定CC是单键还是双键。此时，通过对正则MO的定 域化处理，可以得到通常意义上的成键图像。

33、 正则MO的定域化处理方法较多，其中较为常用的是NBO方法， 其使用方法是在输入文件中添加关键词： POP=NBO,以乙烯分子为例： %mem=32mb #p b3lyp/3-21g pop=nbo The NBO analysis of ethylene 0,1 C C,1,CC H,1,CH,2,HCC H,1,CH,2,HCC,3,180.,0 H,2,CH,1,HCC,3,180.,0 H,2,CH,1,HCC,4,180.,0 CC=1.31477 CH=1.07363 HCC=121.8867,进行NBO成键分析,NBO计算结果的主要内容：,自然布居分析结果,自然电荷,自然电子组态

34、,C的电子组态为： 2s1.042p3.38,自然键轨道的组成,自然键轨道输出格式说明,该键轨道的电子填充数,轨道类型： BD:成键轨道; CR: 内层轨道; LP: 孤对电子; RY: Rydberg弥散轨道； BD*:反键轨道,键轨道系数,轨道组成,对于C原子，其轨道分别是(与所用基组有关，这里为3-21G)： 1s2s3s 1px2px 1py2py 1pz2pz 故由上述轨道组成可知，该轨道为C 2px之间形成的键,对于2号轨道，为C-C 键，并可见这里C原子采用的sp2杂化 方式，对36号轨道为C-H 键.,7,8为C的内层轨道，920均为弥散轨道，电子填充数较少,21为C-C *轨

35、道，22为C-C *轨道，其余为C-H *轨道,根据各键轨道的能量，可 得到乙烯的成键图象为：,4(C-H),4*(C-H),1(C-C),1(C-C),1*(C-C),1*(C-C),由上结果可见，乙烯中 CC键为双键。,上机练习： 1. 计算CO中，CO键能； a. 首先计算CO分子的总能量: %mem=32mb #p b3lyp/6-31g* The energy of CO 0,1 O C,1,1.138 计算结果为: -113.30945a.u. b. 分别计算单个C原子的能量, 此时需将O原子造成的BSSE 考虑在内:,%mem=32mb #p b3lyp/6-31g* massa

36、ge The atomic energy of C in CO 0,3 O C,1,1.138 1 nuc 0.0 计算结果为-37.84664a.u.,c. 类似的, 计算CO中O原子的能量:,%mem=32mb #p b3lyp/6-31g* massage The atomic energy of O in CO 0,3 O C,1,1.138 2 nuc 0.0 计算结果为: -75.06280a.u. 则CO中C-O键能为: -37.84664-75.06280+113.30945) =0.40001a.u.=1050.2kJ/mol (exp. 1070.3kJ/mol),2. 分

37、别从正则和定域两种角度分析乙炔的成键情况;,H,C,C,H,X,1,2,3,4,5,%mem=32mb #p b3lyp/6-31g* pop(full,nbo) C2H2 0,1 X C 1 1.0 H 2 r1 1 90.0 C 2 r2 1 90.0 3 180.0 H 2 r3 1 90.0 3 180.0 r1=1.0666 r2=1.205 r3=2.2716,r1,r2,r1,正则MO:,NBO:,(C-H),(C-C),1(C-C),2(C-C),2(C-H),2*(C-H),1(C-C),2(C-C),2*(C-C),1*(C-C),b. 构型优化,化合物构型优化是G98的常

38、用功能之一,构型优化过程是建立 在能量计算基础之上的. 对于通常意义上的稳定构型是指具有最 低能量的构型,即是在势能面上的能量最低点(极小点). 此外G98 也能对非基态构型进行优化, 例如激发态构型以及过渡态构型等.,初始构型,能量计算,计算能量梯度,各原子受力和位移 大小是否满足标准,根据受力情况 得到新的构型,否,输出最终结果,是,因此与单点能计算相比,构型 优化只多调用了计算能量梯 度的模块L701, L702, L703等,(1) 构型优化过程说明:,(2) 构型优化的输入: 在单点能计算的输入基础上,添加关键词 OPT即可.,%mem=32mb %chk=H2O #P B3LYP/

39、6-311G OPT POP=FULL The geometry optimization of water 0,1 O H,1,1.0 H,1,1.0,2,105.4,说明: 在默认情况下, 构型优化过程中将化合物的自由度全部放开, 若 要对化合物的局部构型进行优化, 此时分子构型需采用内坐标 方式输入, 同时关键词OPT需选取z-matrix选项, 例如:,%mem=32mb %chk=H2O #P B3LYP/6-311G OPT=z-matrix POP=FULL The geometry optimization of water 0,1 O H,1,1.0 H,1,1.0,2,a

40、a=105.4,采用内坐标方式进行构型优化,要优化的参数用变量来表示,固定的构型参数直接给出,变量的初始值,本例仅对H2O中的H-O-H 键角进行优化, 而O-H键长 固定不便,另一种等价表示方法是: %mem=32mb %chk=H2O #P B3LYP/6-311G OPT=z-matrix POP=FULL The geometry optimization of water 0,1 O H,1,r H,1,r,2,a Variables a=105.4 Constants r=1.0,表明该部分参数在构型优化过程中为变量,表明该部分参数在构型优化过程中为常量,OPT关键词的其它辅助选项

41、: 该关键词选项较多,具体请参考G98的帮助文件,其中较为常用 的参数有: MaxCycle= 定义构型优化的最大轮数, 例如: opt(maxcycle=30) TS 优化对应于过渡态的构型, 例如opt(ts) QST2 根据用户给出的反应物和生成物的构型, 进行反应过渡 态构型的优化 QST3 与QST2类似, 只是用户需另外单独给出反应过渡态的 初始构型,(3) 构型优化结果的输出: 以H2O为例: (B3LYP/6-311G),与单点能计算相比,需多调用L7模块,给出了要优化变量的初始值,构型优化最大轮数可 用maxcycle选项设定,接下来的L2L6模块与单点能计算相同,不再说明,

42、L7模块输出:,各原子受力情况,最大值,均方根值,根据受力和所采用的优化 方法,确定新的构型参数,四个构型收敛指标,当前值,标准阈值(可用loose参数修改),重复上述过程对新构型计算能量及其梯度,直至满足收敛条件:,(4) 构型优化说明: 要缩短构型优化时间,需尽可能给出较为准确的初始构型, 例如 采用X衍射实验结果等等; 对于较大体系的构型优化, 为了缩短机时, 可采用分步优化的方 法, 即首先采用较小的基组进行优化, 再将该步骤得到的构型作 为初始构型, 采用较大的基组作进一步的构型优化. 该方法尤其 适合于初始构型不太确定的情形; 由于构型优化涉及到多变量的优化过程, 其最终的结果受到

43、初始 构型的影响较大, 即不能保证所得的构型对应于能量极小点. 这 也是由于化合物势能面的复杂性引起的, 例如:,local minimum,为了保证得到的构型对应于能量极小点, 通常需在构型优化的 基础上进一步进行频率计算, 若计算结果存在明显虚频, 则得到 的构型并非对应于能量极小点.,除了初始构型对优化结果有较大影响外, 体系对称性的限制也 将影响到最终得到的构型. 默认情况下, 在构型优化过程中, 体 系的对称性是保持不变的, 即分子所属的点群是不变的. 例如, 以NH3为例, 如果初始构型将NH3定义为平面型(D3h), 则最终 是得不到属于C3V群的NH3. 为了克服该困难, 可用

44、nosymm关 键词来取消对称性上的限制, 例如:,%mem=32mb %chk=H2O #P B3LYP/6-311G OPT NOSYMM The geometry optimization of water 0,1 O H,1,1.0 H,1,1.0,2,105.4,此时, H2O对称性将放开, 不再为C2v群, 而是C1群,c. 振动光谱(IR)计算: 在构型优化基础上, 通过进一步计算能量的二阶导数, 可求得 力常数，进而得到化合物的振荡光谱。与单点能和构型优化 相比，IR计算需调用L10和11模块(包括L1002, 1101, 1102, 1110等)。 (1) IR计算的输入：

45、关键词为FREQ，例如：,%mem=32mb %chk=H2O #P B3LYP/6-311G OPT FREQ Geometry optimization and the frequency calculation of H2O 0,1 O H,1,1.0 H,1,1.0,2,105.4,在构型优化的基础上作进一步的 频率计算,说明： 必须注意，只有对稳定构型进行IR计算才有意义，即所要进 行IR分析的构型是必须事先经过优化，而且构型优化与IR 计算必须使用相同的计算方法和基组，否则计算结果也是没 有意义的。 以H2O为例，要对其进行IR分析须经过下述步骤： a. 首先对H2O的构型进行优化

46、，此时输入文件为：,%mem=32mb %chk=H2O #P B3LYP/6-311G OPT Geometry optimization of H2O 0,1 O H,1,1.0 H,1,1.0,2,105.4,b. 在上一步计算结果的基础上，进一步进行IR计算:,%mem=32mb %chk=H2O #P B3LYP/6-311G geom=check guess=check freq Frequency calculation of H2O 0,1,i) 在该步骤所采用的计算方法和基组与前一步骤必须相同； ii) geom=check含义是从checkpoint文件(h2o.chk)中

47、读取构型， 注意的是该chk文件是在上一步计算过程中产生的，在opt 结束后，该文件保存了H2O经优化后的构型。由于已确定 从chk文件中读取构型，因此在输入文件的后面不要再给 出H2O的构型，否则计算将出错； iii) guess=check含义是从checkpoint文件中读取初始波函数， 由于该波函数已经SCF迭代优化过，因此可直接作为初始 波函数，这样可以减少FREQ计算中SCF所耗费的时间；,上述分步计算过程较为复杂，其等价输入文件为：,%mem=32mb %chk=H2O #P B3LYP/6-311G OPT FREQ Geometry optimization and the

48、frequency calculation of H2O 0,1 O H,1,1.0 H,1,1.0,2,105.4,FREQ的其它常用选项有： raman 在IR计算的同时，求解raman光谱强度，对于DFT和 MPn方法，在缺省情况下没有计算raman强度，此时 需使用该选项，即freq=raman，相应地，计算时间将 增加；,freq计算需要消耗较多的内存和硬盘，而且比较费时，一般 只适用于中小体系；,(2) IR计算结果分析：,偶极矩，极化率，超极化率,振动频率,依次为： 频率大小；约化质量；力常数；IR强度； Raman活性(强度)；极化率,简正坐标，通过对简正坐标的分析， 可对振动

49、类型进行归属,简正坐标分析说明：,O,H,H,z,y,1,2,3,对应于O-H键的剪式振动,O,H,H,z,y,1,2,3,对应于O-H键的伸缩振动,O,H,H,z,y,1,2,3,此外,IR计算也给出了体系的一些热力学参数, 包括零点能(ZPE), 焓以及吉布斯自由能等:,说明: 通常IR计算得到的频率要校实验结果来得大些, 若要得到较为 准确的数据, 需用校正因子校正, 该校正因子的数值与所用的计 算方法和基组均有关, 具体参考G98手册; 现有多种程序可用于将IR计算结果图示化, 例如Gaussian View 等, 从而便于简正坐标的分析; 当所优化的构型并非对应于能量极小点时, 将得

50、到较为明显的为 负值的频率, 即虚频, 此时可以通过对简正坐标的分析, 推测稳定 构型, 从而消除虚频的出现.,有一数值较大的虚频, 表面平面性NH3并非对应于势能面最低点,例, 对平面型NH3的freq计算结果为:,N,H,H,H,x,y,练习: 分别对重迭和交错式乙烷的构型进行优化, 并进行IR计算, 确定它们构型的相对稳定性.,量化计算的一般过程: 1. 首先根据计算目的选取关键词, 对于本例, 由于涉及到构型 优化, 则应选取opt关键词, 同时因要进行IR计算, 故也将使 用到freq关键词; 2. 其次是选择合适的计算方法和基组, 对于本例我们将采用 密度泛函理论的B3LYP方法和

51、6-31G基组, 对于该方法和基 组的标准表示方法是: B3LYP/6-31G; 3. 确定体系的构型(包括电荷和自旋多重度), 构型描述可用直 角坐标和内坐标(z-matrix), 具体根据那种表示方法方便来选 取, 对于本例采用内坐标表示更为简便; 4. 编写G98输入文件, 投入运行.,对重迭式乙烷：,1,2,3,4,5,6,7,8,1(2),3(7),4(6),5(8),方法一：构型优化和IR计算一步完成,%mem=64mb %chk=c2h6 #p b3lyp/6-31g opt freq Opt of C2H6 0,1 C C 1 r1 H 1 r2 2 a1 H 1 r2 2 a

52、1 3 120.0 H 1 r2 2 a1 3 240.0 H 2 r2 1 a1 4 0.0 H 2 r2 1 a1 3 0.0 H 2 r2 1 a1 5 0.0 r1=1.55 r2=1.0 a1=111.0,1 2 3 4 5 6 7 8,计算方法,所用基组,1,2,3,5,6,7,8,4,对同一化合物, 内坐标的表示方式可有多种, 但 为了能够更好的反映出体系的对称性, 需要选取 一种较好的表达方式, 通常内坐标表示中的前三 个原子的选取较为重要, 一般原则是: 1. 1号和2号原子尽量处在体系的对称轴上; 2. 若体系有对称面, 则该3个原子最好处在该对称 面上. 因此, 在书写内

53、坐标前, 最好先对体系的对称元素 分布进行考察. 在左图中, 按照其编号顺序来描述乙烷的构型, 较难反映出其D3h点群, 除非要给出非常精确的键 角和二面角.,方法二：构型优化和IR计算分步进行：,%mem=64mb %chk=c2h6 #p b3lyp/6-31g opt Opt of C2H6 0,1 C C 1 r1 H 1 r2 2 a1 H 1 r2 2 a1 3 120.0 H 1 r2 2 a1 3 240.0 H 2 r2 1 a1 4 0.0 H 2 r2 1 a1 3 0.0 H 2 r2 1 a1 5 0.0 r1=1.55 r2=1.0 a1=111.0,先构型优化,在

54、构型优化的基础上进行IR计算,%mem=64mb %chk=c2h6 #p b3lyp/6-31g freq geom=check guess=check IR calc. of C2H6 0,1,从c2h6.chk中读取构型,从c2h6.chk中读取波函数 (该法对于大尺度体系可 缩短SCF时间),文件名应相同,说明：chk文件存放经优化过的构型和波函数 等信息，在作进一步的计算时，可充分 利用该文件的内容；,构型优化结果 (本例为自由度全放开 的优化, 因此变量较 多, 为3N-6个, 但因受 到对称性的限制, 有 些变量的数值是相 同的),若只对部分变量进行优化, 需采用opt=z-ma

55、trix选项 %mem=64mb %chk=c2h6 #p b3lyp/6-31g opt=z-matrix Opt of C2H6 0,1 C C 1 1.55 H 1 r2 2 a1 H 1 r2 2 a1 3 120.0 H 1 r2 2 a1 3 240.0 H 2 r2 1 a1 4 0.0 H 2 r2 1 a1 3 0.0 H 2 r2 1 a1 5 0.0 r2=1.0 a1=111.0,此时出现在z-matrix中的数值均认为是 常量, 在构型优化过程中是不变的, 因此 采用该选项时, 需注意那些参数要固定, 那些参数需放开优化,IR计算结果:,对构型优化的补充说明: 除了采

56、用opt关键词外, 有时候也可使用scan关键词来扫描化合 物的势能面, 以C2H6为例, 考察旋转角的影响:,1(2),3,4,5,(6),(7),(8),%mem=64mb #p b3lyp/6-31g scan nosymm Scan of C2H6 0,1 C C 1 r1 H 1 r2 2 a1 H 1 r2 2 a1 3 120.0 H 1 r2 2 a1 3 240.0 H 2 r2 1 a1 4 d1 H 2 r2 1 a1 3 d1 H 2 r2 1 a1 5 d1 r1=1.5481 r2=1.0956 a1=111.7037 d1 0.0 10 12.0,因在扫描过程中

57、分子的点群将改 变, 故需取消对称 性限制,依次为初始值, 扫描点数, 步长,由图可见, 乙烷交错式和重迭式之间存在约11kJ/mol的势垒,Scan的方法只适用于维数较低的情形, 有时为了寻找构型流变的 主要影响因素时可用该方法, 此外也适用于单变量体系.,过渡态(TS)构型的优化 过渡态构型简单言之, 即为具有一个明显虚频的构型, 从能量 角度上看, 过渡态可以看作是势能面上的一阶鞍点, 它在自由度 为N的能量空间中, 只在其中一个自由度方向能量为极大值, 而 在其它N-1个自由度方向上为极小值. 根据已知构型情况, G98提供了三种TS构型优化方法, 它们的 关键词分别为: OPT=TS根据用户提供的初始构型优化TS构型, 此时用户只 需提供一个初始构型; OPT=QST2 根据用户给出的反应物和生成物的构型来优化TS 构型, 此时用户需预先给出两个构型, 分别对应于 反应物和生成物; OPT=QST3 在OPT=QST2的基础