高斯软件入门

量子化学计算方法,章永凡福州大学化学系2003年3月,课程主要内容,一、简介；二、Gaussian98程序的使用及其结果处理；三、其它相关软件的使用；主要参考资料:G98用户手册或G98的帮助文件相应网址:,简介,量子化学软件目的在于将量子化学复杂计算过程程序化，从而便于人们的使用、提高计算效率并具有较强的普适性。绝多数量子化学程序是采用Fortran语言编写的(Fortran77或Fortran90)，通常有几千至上万行语句。,软件分类,计算原理,基于从头算方法(abinitio)Gaussian、ADF、Dalton、Gamess、Crystal等,基于半经验或分子力学方法MOPAC、EHMO、NNEW3等,研究对象,有限尺度体系(分子、簇合物等)Gaussian、ADF、Dalton、Gamess、MOPAC、EHMO等,无限周期重复体系(晶体、固体表面等)Crystal、ADF、NNEW3等,本研究室目前常用的量化软件：Gaussian98:由Pople等人编写，经过几十年的发展和完善，现该软件已成为国际上公认的、计算结果具有较高可靠性的量子化学软件，它包含从头算、半经验以及分子力学等多种方法，可适用于不同尺度的有限体系，除了部分稀土和放射性元素外，它可处理周期表中其它元素形成的各种化合物；ADF2000：该软件主要集中于采用各种密度泛函方法来研究化合物的电子结构，它尤其适合于研究含有过渡金属的化合物，例如金属团簇分子等；Crystal98：该软件由意大利都灵大学理论化学研究所开发，采用从头算方法来研究固体及表面的电子结构；Nnew3：由Hoffmann小组编写的基于EHT方法的量化软件，主要用于处理固体及表面电子结构。,要解决的问题,当前的研究状况，包括计算水平、已解决和尚未解决的问题,计算过程,化合物构型的确定，具体途径包括：实验测定数据、利用软件进行构造等,根据现有的计算条件、模型的大小以及所要解决的问题，选择可行的计算方法,对计算结果进行加工和提取有用的信息，一般包括能量、轨道组成、电荷、布居和成键分析等,计算模型和方法的选取是保证计算结果可靠性的关键，最理想的情况是：1.所选取的计算模型与实际情形一致；2.采用高级的计算方法。但是，由于受到计算软硬件的限制，在多数情况下，很难同时做到上述两点要求，实际操作中，当计算模型较大时，只能选择精确度较低的计算方法，只有对较小的模型才能选取高级的计算方法。因此，当确定了一种计算模型和方法后，需对其进行验证，以保证计算结果的可靠性。假使当前的研究对象是化合物A，可通过下列途径进行验证：1.与A化合物现有实验结果之间的比较；2.若无实验方面的报道，可对与A类似的化合物B进行研究，此时以B的实验结果作为参照；3.当上述方法行不通时，可以采用较大模型和较为高级的计算方法得到的计算结果作为参照，该方法主要用于系列化合物的研究：如对A1,A2,A3，先用大模型和基组对A1进行研究，然后以该结果为参照，确定计算量适中的模型和方法并应用于A2,A3。,Gaussian98程序的使用,G98的安装和运行；G98的功能和程序结构；输入文件的编写与主要功能的使用；补充说明；,G98程序的安装和运行,1.G98程序的安装：(1).确定运行平台：Windows或Linux？(2).对Windows平台：直接运行setup.exe，其余步骤按提示操作即可；也可将其它机器上将已安装好的G98直接拷贝到本机，但需设置运行环境。对Linux平台：a.若G98是经过压缩过的(文件结尾为gz)，用gunzip命令解压：例如：gunzipg98.linux.tar.gzb.若G98是打包的(文件结尾为tar)，用tar命令将其释放：例如：tarxvfg98.linux.tarc.设置环境变量，以cshell为例，在用户根目录下的.cshrc文件添加下列内容：(也可在执行g98前逐条运行)setenvg98root/home/test/src(设置g98所在目录，根据实际情况修改)source$g98root/g98/bsd/g98.login(激活g98运行时所需环境变量)setenvGAUSS_EXEDIR$g98root/g98(设置g98可执行文件所处目录),注：对Linux平台，运行g98时，需注意权限问题，可用chmod命令更改权限，将所安装的g98对所有用户开放。,2.G98程序的运行：(1).对Windows平台：a.对于刚安装好的g98，先检查环境设置情况：,左侧至上而下依次为：默认的文本编辑器；g98可执行文件所在目录；计算中间结果存放目录；缺省的计算结果存储目录；缺省的输入文件所在目录；PDB分子构型浏览器；右侧至上而下依次为：设置显示属性(如背景色等)；设置文本编辑器属性；计算过程控制属性(尤其是批作业过程)；Default.Rou文件的编辑(该文件内容为默认情况下，计算所花费的内存及硬盘大小),需设置正确，否则运行将出错！,b.编写或打开g98输入文件,点击RUN，并给定输出文件名后开始运行,c.g98运行过程的控制：,最上行按钮的功能从左至右依次为：开始运行g98；暂停进程；运行至下一模块(link)时暂停进程；重新启动进程；清除进程(停止运算)；编辑批作业；运行完当前任务后，暂停批作业；停止批作业的运算；观看计算结果；打开文本编辑器；,不要随意点击！,(2).Linux平台：a.编写输入文件：用vi命令编写或在Windows下编写完毕后ftp至Linux系统；b.运行g98：g98输出文件名B3LYP=-1.1647796;MP2=-1.129582MP3=-1.1339601;MP4D=-1.1355271;MP4DQ=-1.13547MP4SDQ=-1.13547;MP4SDTQ=-1.13547;CCSD=-1.1361895CISD=-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-21G(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=5GF数目为3*5=15321g：为劈裂(split)基组，其含义是：内层的每个AO用3个GF描述，价层的AO劈裂为两组，分别用2个和1个GF描述。显然，321g的GF数与sto-3g是相同的。如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：1s2s2p3s3p内层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数为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：1s2s2p3s3p内层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个GF6311g：也为劈裂基组，自是价层的AO劈裂为3组，分别用3个、1个和1个GF描述。对于4-31g和6-21g类似。,Iii)极化(polarization)函数的使用：在实际计算中，有时需在上述标准基组的基础上，添加一个或多个极化函数，极化函数是指具有比原子价轨道更高角量子数的高斯函数。例如H的价轨道为1s，则其极化函数为p型GF，同样对C、O等价层为p轨道的原子，它们的极化函数应为d型或f型轨道，类似地，对于过渡金属原子的极化函数为f型轨道。极化函数的使用目的在于使原子价轨道在空间取向上变得更“柔软”，从而使之易于与其它原子的轨道成键：,例如对于羰基基团中的C，O原子，它们极化轨道(即d轨道)之间形成的d轨道，使得C原子的p轨道朝O方向极化，O的p轨道向C方向极化，从而增强了C与O之间的作用。,由于极化轨道的使用增强价轨道的空间柔软程度，因此其使用场合主要在于环状化合物，例如有机环状化合物以及含有桥联结构的金属簇合物等。使用方法：在标准基组后加“*”，“*”，(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个具有更高角动量的极化函数;对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)的使用:对于带有较多电荷的体系,采用标准的基组来描述是不够的,此时需添加弥散函数,以增加价轨道在空间上的分布范围,即:极化函数用于改进价轨道的角度分布,弥散函数则用于改进价轨道的径向分布.弥散函数是指具有较小轨道指数的高斯函数,其表示方法是在标准基组后加上“+”或“+”,如6-31+G,6-31+G*等,其中第一个“+”表示对非H原子添加弥散函数,第二个“+”则对H原子添加弥散函数.弥散函数主要用于带有电荷的体系(包括离子)以及弱作用体系.,2)赝势基组：关键词：lanl1mb,lanl2mb,lanl1dz,lanl2dz,cep-4g,cep-31g,cep-121g,sdd等说明：通常对于重元素采使用赝势基组，其中，对于过渡金属原子一般使用lanl系列的基组，对于主族元素采用cep系列的基组，sdd基组相对较少使用；赝势基组的适用范围要较全电子基组来得广，一般除了稀土和放射性元素外，均可使用赝势基组；赝势基组是将内层原子的影响用势函数来代替，因此，在使用时必须注意根据所采用的赝势基组，检查体系的总电子数是否正确！对lanl1系列的基组，只考虑价层电子；如V原子为5.对lanl2系列的基组，除了价层电子外，还考虑了次外层的电子(对部分主族元素例外)；如V原子为5+8=13.对其它赝势基组，请在使用时仔细核查。,对于使用较多的lanl系列基组，需注意：对于HNe范围的原子，lanl1mb和lan2mb基组等价于sto-3g最小基组；lanl1dz和lanl2dz基组等价于d95基组;对于含金属金属键体系一般使用lanl2系列，而不用lanl1系列；mb结尾含义为minimalbasis，dz结尾含义为double-zeta，故后者基函数数目要较前者多。,3)混合基组的使用:有时,因需要,对同一化合物中的不同原子采用不同的基组:受基组限制:例如,对于Mo(CO)2+化合物,C和O原子可采用6-31+G基组,但Mo原子不在6-31G基组的使用范围之内.受体系大小限制:当体系较大时,不可能对所有原子均采用较大的基组,此时可对其中的局部原子采用精度较高的基组来描述,而对其余原子采用小基组来描述,例如:,混合基组输入方法:关键词:GEN和PSEUDO=READ(对赝势基组)%mem=32mb%chk=hf#pb3lyp/genpop=fullscfcyc=500Mo(CO)2+(Line)1,6MoX11.0C1r1290.0O1r2290.030.0C1r1290.03180.0O1r2290.03180.0r1=2.228r2=3.329Mo03-21G*CO06-31+G*,整数0为结束标志,4个*用于分隔不同基组,改行也可用数字代替,即23450(编号时不包含虚原子),%mem=32mb%chk=hf#pb3lyp/genpseudo=readpop=fullscfcyc=500Mo(CO)2+(Line)1,6MoX11.0C1r1290.0O1r2290.030.0C1r1290.03180.0O1r2290.03180.0r1=2.228r2=3.329Mo0Lanl2mb*CO06-31G*Mo0Lanl2mb,使用赝势基组时,应注意不能漏!,基组定义部分,赝势定义部分,注意不能漏,另外其后没有*号,赝势部分和基组部分之间有一空行,(3)G98输出结果:下面以H2O能量计算的输出结果为例:,L1输出G98版权信息,%部分内容,关键词部分,计算所需调用的模块,Title内容,电荷,S以及分子构型,各原子的直角坐标,距离矩阵,需核查!,注意:计算结果是以该坐标系为准,需核查!,MO初始猜测,自洽场迭代求解部分,收敛指标,迭代次数,(通常接近2.0),(布居分析结果),MO对称性及能级(a.u.),Mulliken键级和电荷,偶极矩及多极矩,对主要计算结果进行总结,(4)能量计算中使用到的其它常用关键词:SCF收敛指标的设置:关键词:scf(conver=?)含义:当SCF前后电荷密度的差值小于10-?时,SCF收敛.例:scf(conver=5)表示收敛指标为10-5.说明:对于能量计算,缺省情况下收敛指标为10-4.对于构型优化,则为10-8.SCF迭代轮数的设置:关键词:scf(maxcyc=?)含义:SCF最大的允许迭代次数为?.例:scf(maxcyc=100)表示最大迭代次数为100次.说明:缺省情况下,SCF迭代次数为64次,通常对于含有过渡金属的体系,SCF较难收敛,迭代次数会多些.,分子轨道系数的输出:关键词:pop=full说明:有时为了分析MO成分,则需利用该关键词输出各MO的成分,以H2O为例:,轨道对称性,O表示占据轨道;V为空轨道,轨道能级,AO轨道,轨道系数,HOMO,LUMO,注:输出的MO系数是未经归一化过的;1s,2s等符号并非真正意义上的AO,它与计算所采用的基组有关.对于上例,采用的是6-31G基组,可知,对于价层AO是分裂为两组的,故对于O原子,输出的2s和3s实际上均为价轨道,其它类推.当基函数较多时,若只考察前线附近的轨道成分,可用关键词:pop=regular,此时只给出5个占据轨道和5个空轨道组成.练习:计算O2的基态能量,并分析MO轨道组成,(5)能量计算的几个问题：基组的选取：理论上而言，计算中所选取的基组越大，计算结果越准确，但由于受到硬件上的限制，需根据实际情况选择基组。此外，对于环状分子或存在作用的体系，通常需要考虑使用极化基组；对于带有较多电荷的体系，或考察弱相互作用的体系，通常需考虑使用弥散基组；以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.此时，因带较多负电荷，故弥散函数的影响较极化函数来得大。,体系多重度的选取：对于复杂体系，特别是含有过渡金属原子的体系，多重度的选取凭经验不易得到，此时，可通过比较不同多重度时体系的能量来确定基态能量。以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原子基态能量，结果为：-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#pb3lyp/6-31g*massageBSSEofO20,3O0.0.00.0O1.2080.00.02Nuc0.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和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#pROB3LYP/STO-3GPOP=FULLThesinglepointcalc.OfO20,3O0.0.00.0O1.2080.00.0,为了便于计算结果分析，采用限制性开壳层计算方法，及将和电子强制配对；,输出所有分子轨道,首先确定计算结果中各原子所处的位置：,2,1,o,z,Xory,z,+,1,2,因1s为内层电子，故O与O之间成键情况可不考虑,No.1,No.2与上一MO类似，也以1s成分为主，处在能量较低的位置,z,+,+,z,+,1,2,No.3,4分别为O2s轨道之间形成的成键(SGG)和反键(SGU)轨道：,+,z,+,-,No.5为O2pz轨道形成的成键轨道：,z,+,z,No.6和7轨道分别为O2px和2py形成的成键轨道：,+,No.8和9轨道分别为O2px和2py形成的反键轨道：,z,+,No.10轨道为O2pz形成的反键轨道：,z,+,根据上述分析结果以及各轨道的能级，可以得到O2的MO分布为：,2s,*2s,2pz,2px,2py,*2px,*2py,*2pz,电荷分布：在能量计算结果中，除了能量以及上述的MO成分外，其它有用的信息包括原子电荷和键级，对于原子电荷：以H2O为例(B3LYP/STO-3G)：,Mulliken电荷分布,说明：电荷的绝对值是没有意义的，其数值受到所用方法,尤其是基组的影响较大：以H2O为例：方法/基组OHB3LYP/STO-3G-0.3290.165B3LYP/3-21G-0.6370.318B3LYP/6-31G-0.7050.353HF/6-31G-0.7950.398因此，计算得到的电荷只具有相对意义，即只能考察其相对数值或变化趋势。当然，当考察变化趋势时，如考察H2X(X=O,S,Se)系列化合物中，X为不同原子时，对H原子电荷的影响时，必须使用相同的计算方法和基组。,键级：键级的大小反映了某根键的强度，同样地，键级的数值也受到计算方法和基组的影响：,O-H键的Mulliken键级,以H2O为例：方法/基组O-H键级B3LYP/STO-3G0.243B3LYP/3-21G0.239B3LYP/6-31G0.231B3LYP/6-31+G0.214HF/6-31G0.249由上结果可见，键级数值与所用的基组和计算方法是有联系的，但其影响不象电荷那么显著。此外，必须注意到，弥散基组对键级的数值影响较大。因此，在实际应用中，对于电荷和键级，我们更关心的是它们的变化情况而不是它们的绝对数值。此外，必须强调的是，上述结果是采用Mulliken布居分析的方法，若采用其它布居分析方法，结果将有所不同。,自然键轨道(NBO)分析方法：我们前面所讲分子轨道均为正则分子轨道，它是未经定域化处理，这导致计算结果与我们通常的成键概念有所不同。例如在乙烯分子中，碳碳之间为双键，但在正则MO中，反映C与C之间成键作用的MO可能有多个，因此根据正则MO的结果，我们无法断定CC是单键还是双键。此时，通过对正则MO的定域化处理，可以得到通常意义上的成键图像。正则MO的定域化处理方法较多，其中较为常用的是NBO方法，其使用方法是在输入文件中添加关键词：POP=NBO,以乙烯分子为例：%mem=32mb#pb3lyp/3-21gpop=nboTheNBOanalysisofethylene0,1CC,1,CCH,1,CH,2,HCCH,1,CH,2,HCC,3,180.,0H,2,CH,1,HCC,3,180.,0H,2,CH,1,HCC,4,180.,0CC=1.31477CH=1.07363HCC=121.8867,进行NBO成键分析,NBO计算结果的主要内容：,自然布居分析结果,自然电荷,自然电子组态,C的电子组态为：2s1.042p3.38,自然键轨道的组成,自然键轨道输出格式说明,该键轨道的电子填充数,轨道类型：BD:成键轨道;CR:内层轨道;LP:孤对电子;RY:Rydberg弥散轨道；BD*:反键轨道,键轨道系数,轨道组成,对于C原子，其轨道分别是(与所用基组有关，这里为3-21G)：1s2s3s1px2px1py2py1pz2pz故由上述轨道组成可知，该轨道为C2px之间形成的键,对于2号轨道，为C-C键，并可见这里C原子采用的sp2杂化方式，对36号轨道为C-H键.,7,8为C的内层轨道，920均为弥散轨道，电子填充数较少,21为C-C*轨道，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#pb3lyp/6-31g*TheenergyofCO0,1OC,1,1.138计算结果为:-113.30945a.u.b.分别计算单个C原子的能量,此时需将O原子造成的BSSE考虑在内:,%mem=32mb#pb3lyp/6-31g*massageTheatomicenergyofCinCO0,3OC,1,1.1381nuc0.0计算结果为-37.84664a.u.,c.类似的,计算CO中O原子的能量:,%mem=32mb#pb3lyp/6-31g*massageTheatomicenergyofOinCO0,3OC,1,1.1382nuc0.0计算结果为:-75.06280a.u.则CO中C-O键能为:-37.84664-75.06280+113.30945)=0.40001a.u.=1050.2kJ/mol(exp.1070.3kJ/mol),2.分别从正则和定域两种角度分析乙炔的成键情况;,H,C,C,H,X,1,2,3,4,5,%mem=32mb#pb3lyp/6-31g*pop(full,nbo)C2H20,1XC11.0H2r1190.0C2r2190.03180.0H2r3190.03180.0r1=1.0666r2=1.205r3=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的常用功能之一,构型优化过程是建立在能量计算基础之上的.对于通常意义上的稳定构型是指具有最低能量的构型,即是在势能面上的能量最低点(极小点).此外G98也能对非基态构型进行优化,例如激发态构型以及过渡态构型等.,初始构型,能量计算,计算能量梯度,各原子受力和位移大小是否满足标准,根据受力情况得到新的构型,否,输出最终结果,是,因此与单点能计算相比,构型优化只多调用了计算能量梯度的模块L701,L702,L703等,(1)构型优化过程说明:,(2)构型优化的输入:在单点能计算的输入基础上,添加关键词OPT即可.,%mem=32mb%chk=H2O#PB3LYP/6-311GOPTPOP=FULLThegeometryoptimizationofwater0,1OH,1,1.0H,1,1.0,2,105.4,说明:在默认情况下,构型优化过程中将化合物的自由度全部放开,若要对化合物的局部构型进行优化,此时分子构型需采用内坐标方式输入,同时关键词OPT需选取z-matrix选项,例如:,%mem=32mb%chk=H2O#PB3LYP/6-311GOPT=z-matrixPOP=FULLThegeometryoptimizationofwater0,1OH,1,1.0H,1,1.0,2,aa=105.4,采用内坐标方式进行构型优化,要优化的参数用变量来表示,固定的构型参数直接给出,变量的初始值,本例仅对H2O中的H-O-H键角进行优化,而O-H键长固定不便,另一种等价表示方法是:%mem=32mb%chk=H2O#PB3LYP/6-311GOPT=z-matrixPOP=FULLThegeometryoptimizationofwater0,1OH,1,rH,1,r,2,aVariablesa=105.4Constantsr=1.0,表明该部分参数在构型优化过程中为变量,表明该部分参数在构型优化过程中为常量,OPT关键词的其它辅助选项:该关键词选项较多,具体请参考G98的帮助文件,其中较为常用的参数有:MaxCycle=定义构型优化的最大轮数,例如:opt(maxcycle=30)TS优化对应于过渡态的构型,例如opt(ts)QST2根据用户给出的反应物和生成物的构型,进行反应过渡态构型的优化QST3与QST2类似,只是用户需另外单独给出反应过渡态的初始构型,(3)构型优化结果的输出:以H2O为例:(B3LYP/6-311G),与单点能计算相比,需多调用L7模块,给出了要优化变量的初始值,构型优化最大轮数可用maxcycle选项设定,接下来的L2L6模块与单点能计算相同,不再说明,L7模块输出:,各原子受力情况,最大值,均方根值,根据受力和所采用的优化方法,确定新的构型参数,四个构型收敛指标,当前值,标准阈值(可用loose参数修改),重复上述过程对新构型计算能量及其梯度,直至满足收敛条件:,(4)构型优化说明:要缩短构型优化时间,需尽可能给出较为准确的初始构型,例如采用X衍射实验结果等等;对于较大体系的构型优化,为了缩短机时,可采用分步优化的方法,即首先采用较小的基组进行优化,再将该步骤得到的构型作为初始构型,采用较大的基组作进一步的构型优化.该方法尤其适合于初始构型不太确定的情形;由于构型优化涉及到多变量的优化过程,其最终的结果受到初始构型的影响较大,即不能保证所得的构型对应于能量极小点.这也是由于化合物势能面的复杂性引起的,例如:,localminimum,为了保证得到的构型对应于能量极小点,通常需在构型优化的基础上进一步进行频率计算,若计算结果存在明显虚频,则得到的构型并非对应于能量极小点.,除了初始构型对优化结果有较大影响外,体系对称性的限制也将影响到最终得到的构型.默认情况下,在构型优化过程中,体系的对称性是保持不变的,即分子所属的点群是不变的.例如,以NH3为例,如果初始构型将NH3定义为平面型(D3h),则最终是得不到属于C3V群的NH3.为了克服该困难,可用nosymm关键词来取消对称性上的限制,例如:,%mem=32mb%chk=H2O#PB3LYP/6-311GOPTNOSYMMThegeometryoptimizationofwater0,1OH,1,1.0H,1,1.0,2,105.4,此时,H2O对称性将放开,不再为C2v群,而是C1群,c.振动光谱(IR)计算:在构型优化基础上,通过进一步计算能量的二阶导数,可求得力常数，进而得到化合物的振荡光谱。与单点能和构型优化相比，IR计算需调用L10和11模块(包括L1002,1101,1102,1110等)。(1)IR计算的输入：关键词为FREQ，例如：,%mem=32mb%chk=H2O#PB3LYP/6-311GOPTFREQGeometryoptimizationandthefrequencycalculationofH2O0,1OH,1,1.0H,1,1.0,2,105.4,在构型优化的基础上作进一步的频率计算,说明：必须注意，只有对稳定构型进行IR计算才有意义，即所要进行IR分析的构型是必须事先经过优化，而且构型优化与IR计算必须使用相同的计算方法和基组，否则计算结果也是没有意义的。以H2O为例，要对其进行IR分析须经过下述步骤：a.首先对H2O的构型进行优化，此时输入文件为：,%mem=32mb%chk=H2O#PB3LYP/6-311GOPTGeometryoptimizationofH2O0,1OH,1,1.0H,1,1.0,2,105.4,b.在上一步计算结果的基础上，进一步进行IR计算:,%mem=32mb%chk=H2O#PB3LYP/6-311Ggeom=checkguess=checkfreqFrequencycalculationofH2O0,1,i)在该步骤所采用的计算方法和基组与前一步骤必须相同；ii)geom=check含义是从checkpoint文件(h2o.chk)中读取构型，注意的是该chk文件是在上一步计算过程中产生的，在opt结束后，该文件保存了H2O经优化后的构型。由于已确定从chk文件中读取构型，因此在输入文件的后面不要再给出H2O的构型，否则计算将出错；iii)guess=check含义是从checkpoint文件中读取初始波函数，由于该波函数已经SCF迭代优化过，因此可直接作为初始波函数，这样可以减少FREQ计算中SCF所耗费的时间；,上述分步计算过程较为复杂，其等价输入文件为：,%mem=32mb%chk=H2O#PB3LYP/6-311GOPTFREQGeometryoptimizationandthefrequencycalculationofH2O0,1OH,1,1.0H,1,1.0,2,105.4,FREQ的其它常用选项有：raman在IR计算的同时，求解raman光谱强度，对于DFT和MPn方法，在缺省情况下没有计算raman强度，此时需使用该选项，即freq=raman，相应地，计算时间将增加；,freq计算需要消耗较多的内存和硬盘，而且比较费时，一般只适用于中小体系；,(2)IR计算结果分析：,偶极矩，极化率，超极化率,振动频率,依次为：频率大小；约化质量；力常数；IR强度；Raman活性(强度)；极化率,简正坐标，通过对简正坐标的分析，可对振动类型进行归属,简正坐标分析说明：,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计算结果图示化,例如GaussianView等,从而便于简正坐标的分析;当所优化的构型并非对应于能量极小点时,将得到较为明显的为负值的频率,即虚频,此时可以通过对简正坐标的分析,推测稳定构型,从而消除虚频的出现.,有一数值较大的虚频,表面平面性NH3并非对应于势能面最低点,例,对平面型NH3的freq计算结果为:,N,H,H,H,x,y,练习:分别对重迭和交错式乙烷的构型进行优化,并进行IR计算,确定它们构型的相对稳定性.,量化计算的一般过程:1.首先根据计算目的选取关键词,对于本例,由于涉及到构型优化,则应选取opt关键词,同时因要进行IR计算,故也将使用到freq关键词;2.其次是选择合适的计算方法和基组,对于本例我们将采用密度泛函理论的B3LYP方法和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#pb3lyp/6-31goptfreqOptofC2H60,1CC1r1H1r22a1H1r22a13120.0H1r22a13240.0H2r21a140.0H2r21a130.0H2r21a150.0r1=1.55r2=1.0a1=111.0,12345678,计算方法,所用基组,1,2,3,5,6,7,8,4,对同一化合物,内坐标的表示方式可有多种,但为了能够更好的反映出体系的对称性,需要选取一种较好的表达方式,通常内坐标表示中的前三个原子的选取较为重要,一般原则是:1.1号和2号原子尽量处在体系的对称轴上;2.若体系有对称面,则该3个原子最好处在该对称面上.因此,在书写内坐标前,最好先对体系的对称元素分布进行考察.在左图中,按照其编号顺序来描述乙烷的构型,较难反映出其D3h点群,除非要给出非常精确的键角和二面角.,方法二：构型优化和IR计算分步进行：,%mem=64mb%chk=c2h6#pb3lyp/6-31goptOptofC2H60,1CC1r1H1r22a1H1r22a13120.0H1r22a13240.0H2r21a140.0H2r21a130.0H2r21a150.0r1=1.55r2=1.0a1=111.0,先构型优化,在构型优化的基础上进行IR计算,%mem=64mb%chk=c2h6#pb3lyp/6-31gfreqgeom=checkguess=checkIRcalc.ofC2H60,1,从c2h6.chk中读取构型,从c2h6.chk中读取波函数(该法对于大尺度体系可缩短SCF时间),文件名应相同,说明：chk文件存放经优化过的构型和波函数等信息，在作进一步的计算时，可充分利用该文件的内容；,构型优化结果(本例为自由度全放开的优化,因此变量较多,为3N-6个,但因受到对称性的限制,有些变量的数值是相同的),若只对部分变量进行优化,需采用opt=z-matrix选项%mem=64mb%chk=c2h6#pb3lyp/6-31gopt=z-matrixOptofC2H60,1CC11.55H1r22a1H1r22a13120.0H1r22a13240.0H2r21a140.0H2r21a130.0H2r21a150.0r2=1.0a1=111.0,此时出现在z-matrix中的数值均认为是常量,在构型优化过程中是不变的,因此采用该选项时,需注意那些参数要固定,那些参数需放开优化,IR计算结果:,对构型优化的补充说明:除了采用opt关键词外,有时候也可使用scan关键词来扫描化合物的势能面,以C2H6为例,考察旋转角的影响:,1(2),3,4,5,(6),(7),(8),%mem=64mb#pb3lyp/6-31gscannosymmScanofC2H60,1CC1r1H1r22a1H1r22a13120.0H1r22a13240.0H2r21a14d1H2r21a13d1H2r21a15d1r1=1.5481r2=1.0956a1=111.7037d10.01012.0,因在扫描过程中分子的点群将改变,故需取消对称性限制,依次为初始值,扫描点数,步长,由图可见,乙烷交错式和重迭式之间存在约11kJ/mol的势垒,Scan的方法只适用于维数较低的情形,有时为了寻找构型流变的主要影响因素时可用该方法,此外也适用于单变量体系.,过渡态(TS)构型的优化过渡态构型简单言之,即为具有一个明显虚频的构型,从能量角度上看,过渡态可以看作是势能面上的一阶鞍点,它在自由度为N的能量空间中,只在其中一个自由度方向能量为极大值,而在其它N-1个自由度方向上为极小值.根据已知构型情况,G98提供了三种TS构型优化方法,它们的关键词分别为:OPT=TS根据用户提供的初始构型优化TS构型,此时用户只需提供一个初始构型;OPT=QST2根据用户给出的反应物和生成物的构型来优化TS构型,此时用户需预先给出两个构型,分别对应于反应物和生成物;OPT=QST3在OPT=QST2的基础上,用户还需提供TS的初始构型;其中对于QST2和QST3,用户给出的多个构型描述中,原子的次序应相同.,%mem=64mb#pRHF/6-31G(d)Opt=QST2SiH2+H2-SiH4Reactants0,1SiX11.0H11.48255.0H11.48255.03180.0H1R2A1390.0H1R5A22180.0R=2.0A1=80.0A2=22.0SiH2+H2-SiH4Products0,1SiX11.0H11.48255.0H11.48255.03180.0H1R2A1390.0H1R5A22180.0R=1.48A1=125.2A2=109.5,例:,反应物构型描述,生成物构型描述,有两行Title(Title前后均有一空行),说明:为了保证得到的构型的确对应于过渡态,通常必须作进一步的IR计算,以考察是否只有一个虚频,故可将FREQ与OPT一起使用:,%mem=64mb#pRHF/6-31G(d)Opt=QST2FREQSiH2+H2-SiH4Reactants0,1SiX11.0H11.48255.0H11.48255.03180.0H1R2A1390.0H1R5A22180.0R=2.0A1=80.0A2=22.0SiH2+H2-SiH4Products0,1SiX11.0H11.48255.0H11.48255.03180.0H1R2A1390.0H1R5A22180.0R=1.48A1=125.2A2=109.5,计算结果:,反应物,TS,产物,在得到TS构型后,有时需进一步对反应途径进行考察.虽然连接反应物-TS-和生成物的反应途径可能不只一种,但我们通常最关心的是连接反应物和生产物的能量最低的反应途径(MEP),因其具有最大的统计概率,故也是实际上最有可能发生的过程.在G98中,考察反应途径可用IRC关键词,当使用该关键词时,必须已经获得了TS构型以及对应的力常数,对于力常数可在FREQ计算时获得,也可加上irc=calcfc选项来计算.,%mem=64mb%chk=c2h6#pb3lyp/6-31girc=calcfcnosymmOptofC2H60,1CC1r1H1r22a1H1r22a13120.0H1r22a13240.0H2r21a14d1H2r21a13d1H2r21a15d1r1=1.5481r2=1.0956a1=111.7037d1=0.0,因未知道力常数,故需重新计算力常数,如果在FREQ计算中,保留了相应的check文件,则可用:irc=rcfc来读取力常数.,以重迭式乙烷为例,我们已经知道该构型对应于TS(因其只有一个虚频),则IRC计算如下:,计算结果:在默认情况下,将搜索MEP位于TS左右两测各6个点的构型,各点之间的间隔(搜索步长)为0.1a.u.(相当于反应坐标的间隔),因此实际上,在默认情况下,要优化12个构型,这些优化过的构型处在MEP上.若用户需修改构型数和步长,可用maxpoints和stepsize选项:maxpoints=N求解MEP上2*N各构型stepsize=N搜索步长为N*0.1对于上例,计算结果为:,反应坐标,对应于TS,处在MEP上各构型的参数,D.电子光谱(UV)的计算G98提供了三种电子光谱的计算方法,包括CIS(single-excitationCI),TD(time-dependent)和ZINDO方法,其中ZINDO为半经验方法,下面主要介绍CIS方法的使用.(1)CIS方法的使用:以H2为例:,%mem=64mb#pcis/6-31+g*TheCIScalc.ofH20,1HH10.76,CIS计算结果分析:CIS计算需调用L8和L9模块,其中UV计算结果由L914模块给出:,默认情况下求解三个激发态,主要信息:,该激发态对应的电子态为1u,电子跃迁所需能量,对应的UV波长,跃迁强度,该激发态的主要来源,反映了各分子轨道的贡献