量子化学计算方法

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

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

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

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

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

6、xe，其余步骤按提示操作即可；也可将其它机 器上将已安装好的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所在目录，根据实际情况修改根据实际情况修改

7、) source $g98root/g98/bsd/g98.login(激活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.

9、编写输入文件： 用vi命令编写或在Windows下编写完毕后ftp至Linux系统；b.运行g98： g98输出文件名& 说明：1)末尾的&符号表示将作业提交到后台计算，否则 在用户退出Linux时，作业将终止； 2)若运行g98出错，请检查环境变量是否设置正确， 尤其是用户权限上的问题； c.观看计算结果： 使用vi命令，或采用tail命令跟踪计算输出： tail -f 输出文件名 d.运行过程的控制： 采用top命令观察g98运行到那个模块； 通过renice命令改变进程的优先级来调整g98的运行速度； 采用kill命令终止运算，格式为：kill 作业号 (作业号由top命

10、令可得) G98的功能和程序结构的功能和程序结构1. 主要功能：分子构型的优化基态(Ground state)激发态(Excited state)反应过渡态(Transition state)能量计算基态和激发态能量化学键的键能电子亲合能和电离能化学反应途径和势能面光谱计算IR光谱Raman光谱电子光谱NMR其它功能电荷和电荷密度偶极矩和超极矩热力学参数适用体系：气相和溶液2.程序结构：a.由主引导模块(g98.exe)和各分模块(l?.exe)组成：b.常用模块的功能：L0初始化模块；L1读入输入，根据所给关键词确定将要使用的模块；L101,102,与构型优化和反应过渡态相关的模块；L202

11、输出距离矩阵、判断化合物点群及确定新的坐标系；L301,302309与基组和赝势有关模块；L310,319计算单及双电子积分模块；L401,402SCF初始猜测模块；L502,503,508SCF模块；L601,608Mulliken布居以及自然键轨道分析模块；L701,702计算能量一阶和二阶导数模块；L8?,9?,10?,11?与Post-SCF方法有关模块；L9999进程结束模块；说明：1.根据不同的任务，某些模块需重复调用多次；2.通常耗时较多的模块有：L5,L7,L8,L9,L10,L11等，此外， L8L11这些模块的执行对内存和硬盘的需求较大；3.若L9999未能正常执行完毕，则

12、表明计算过程存在问题，需 检查之；4.可根据各个模块的功能，对g98程序进行简化，例如如果用 户通常只用g98进行能量计算，则可只保留L16和L9999模块 其它模块可以删除去。c.g98运行过程所使用的文件： 在scratch目录下有下列文件： gxx-打头的文件为临时文件，计算结束后将自动删除，其中 对于结尾为inp的文件，记录了当前g98所执行的输入文 件内容，有时可通过该文件确定当前运行作业； chk文件，该文件记录了g98运行的结果，包括分子结构、基 组、分子轨道、电荷密度以及偶极矩等，通常该文件在计 算结束后要保留要保留，以便以后作补充计算之需；说明：1).对于chk文件并不自动产

13、生，需用户自行指定，在大多数情 况下，最好给定最好给定chk文件，并在计算结束后保留之文件，并在计算结束后保留之；2).对于rwf文件，在运行过程中，g98会自行产生gxx打头的rwf 文件，当计算非正常 中断后，可通过更改该文件来续算。 但用户在编写输入文件时，最好还是指定rwf文件，在计算 结束后，再删除。 rwf文件，该文件记录了计算的中间结果，以便在计算过程非 正常中断后用于续算，该文件通常较大，当作业正常结束 后，可删除之。 Default.Rou文件，该文件设置一些系统默认参数，例如 g98运行时内存和硬盘的大小，其内容如下： -M- 256MB(内存大小) -#- MaxDisk

14、=2000MB(硬盘大小)G98输入文件的编写与使用输入文件的编写与使用1.G98输入文件的组成：L0命令部分(可无)关键词部分标题部分体系电荷和自旋多重度分子构型(1) L0命令部分： 该部分内容均以%打头，主要用于指定计算过程所需 内存，chk以及rwf文件名，其典型内容是： %mem=100mb(计算所需内存，若无，由Default.Rou指定) %chk=h2o(chk文件的名称为h2o.chk) %rwf=h2o(中间结果文件名称为h2o.rwf)(2)关键词部分： 该部分内容由一个或多个关键词组成，用于指定了计 算类型和计算输出的控制等。(3)标题部分： 由一行文本组成，该内容是必

15、需的。(4)电荷及自旋多重度： 体系所带电荷以及自旋多重度S=2s+1=成单电子数+1 例如H2O,S=1；Ti原子的基态S=3(5)分子构型描述部分： 有三种描述分子构型的方法：a.直角坐标系方法：(适用于全自由度优化情况) 格式为：元素符号 x y z 例如： O -0.464 0.177 0.0 H -0.464 1.137 0.0 H 0.441 -0.143 0.0 说明：1)元素符号大小写均可，也可直接用原子序数； 2)有时为了便于区别，可在元素符号后加一整数，如： O -0.464 0.177 0.0 H1 -0.464 1.137 0.0 H2 0.441 -0.143 0.0

16、 3)x,y,z数值必须以小数格式小数格式输入： O -0.464 0.177 0 () O -0.464 0.177 0.() 4)g98的数据输入均为自由格式，即除了用空格来分隔 数据外，也可用逗号或混合使用；b.内坐标(z-matrix)方法：(适用于构型的局部优化) 内坐标与直角坐标之间的区别在于，它侧重于从原子之间的 键连角度来描述原子间的相对位置，具体参数包括： 1)键长：( 需用两个原子描述) 即两个原子间的距离，注：该两个原子并非要具有化学直 观意义上的成键。此外，在默认情况下，键长单位为埃。 2)键角：(需用三个原子描述) 确定了二根键之间的夹角，默认单位为度，范围为-180

17、 180deg之间，但二根键不能共线，即键角不能等于不能等于0和 180deg。 3)二面角：(需用四个原子描述) 二面角加上键长和键角就确定了四个原子的位置，其默认 单位为deg，范围为-360360deg，可以为360deg。当二 面角等于0，180和360deg时四个原子共面。键长、键角和二面角数目的总和=3N-6内坐标的输入格式为： 原子1,原子2,键长,原子3,键角,原子4,二面角1234键长键角二面角例1：OHH表示一：OH,1,1.0O,1,1.2,2,104.0H,3,1.0,1,104.0,2,170.0表示二：O1H1,O1,1.0O2,O1,1.2,H1,104.0H2,

18、O2,1.0,O1,104.0,H1,170.0O1234表示三：OH,1,r1O,1,r2,2,a1H,3,r1,1,a1,2,d1 Variables:(本行内容可省)r1=1.0r2=1.2a1=104.0d1=170.0当对分子的构型进行局部优化局部优化时，需采用该表示方法例2：乙烯CC123456CC 1 1.3H 1 1.0 2 120.0H 1 1.0 2 120.0 3 180.0H 2 1.0 1 120.0 3 0.0H 2 1.0 1 120.0 3 180.0在同侧共面共面但不同侧对同一构型，内坐标的表示并不唯一对同一构型，内坐标的表示并不唯一CC241356HC,1,

19、1.0H,2,1.0,1,120.0C,2,1.3,1,120.0,3,180.0H,4,1.0,2,120.0,1,0.0H,4,1.0,2,120.0,1,180.0虚原子的使用：二根键共线时，为了避免键角为0或180，此时需引入虚原子；此外，有时为了保证所描述的构型符合一定的点群，利用虚原子便于做到这一点。虚原子的符号为X。例1：CO2COOX1234XC,1,1.0O,2,1.1,1,90.0O,2,1.1,1,90.0,3,180.0该键长值可任给COO213OC,1,1.0O,2,1.0,1,180.0例2：NH3NHHHX12345要使输入的构型满足C3v点群，需要准确提供H-N

20、-H键角以及四面体相邻两个平面间的二面角。为此，在3个H所在三角形中心引入一个虚原子X，则：XN 1 2.0H 1 1.0 2 90.0H 1 1.0 2 90.0 3 120.0H 1 1.0 2 90.0 3 -120.0注：在本例中N-X和H-X不能任给说明：1)根据需要，有时可同时用到多个虚原子； 2)在大多数场合，虚原子通常取在对称元素所处位置 或它们相交处；c.直角坐标和内坐标混合输入方法： 对于该方法，只需在采用直角坐标方法输入的原子的元素 符号后加一个整数0即可，例如： X 0 1.0 1.0 1.0 N 1 2.0 H 1 1.0 2 90.0 H 1 1.0 2 90.0

21、3 120.0 H 1 1.0 2 90.0 3 -120.0d.分子构型的输入准确性是保证计算结果可靠性的前提，对 于复杂体系，在计算前均需对所输构型进行检查计算前均需对所输构型进行检查，具体包 括： 可视化处理，即采用一些分子构型软件观察所给构型是 否合理； 在g98运行到L2模块，会给出所输入分子所属点群，此 时，可检查点群是否合理。(6) g98输入文件编辑时的注意事项： 除了可采用g98所提供的输入文件编辑器来编写输入文件外， 在更多场合下，是采用其他文本编辑器来编写，此时应注意 到，在标题部分的前后各有一空行，例如： %chk=h2o# HF/6-31G(d)(此处为空行)wate

22、r energy (此处为空行)0 1O -0.464 0.177 0.0H -0.464 1.137 0.0H 0.441 -0.143 0.0练习：采用内坐标方法输入苯和甲烷的构型。2.g98主要功能的使用： g98功能的使用主要由用户所给的关键词(keyword)内容而定， 在输入关键词时注意以下事项： 1).关键词的输入是自由格式，且不区分大小写； 2).当存在多个关键词时，可用空格、逗号来隔开； 3).通常每个关键词有多个选项(option)，若要选择单个或多 个选项时，书写方式有如下几种： keyword = optionkeyword(option)keyword=(option

23、1, option2, .)keyword(option1, option2, .) 例如：opt=z-matrix，opt(z-matrix)，opt=(z-matrix)三者是 等价的。 opt(z-matrix,maxcycle=20)与 opt=(z-matrix,maxcycle=20)是等价的。 建议统一采用第四种表示方式。多个选项时4).最简单的关键词输入是#或#p，其含义是采用HF方法和 STO-3G基组计算体系的能量；a. 能量的计算： 如何计算一个体系的能量是获取分子各种性质的基础，因此 首先来看如何计算体系的能量，即进行单点能计算：(1). 计算方法的选择： g98提供的

24、常用计算方法有： 1) 半经验方法： 关键词：AM1, PM3, CNDO, INDO, MINDO 它们主要用于大的有机分子体系(由上百个原子组成)，一般 对于含金属体系不适用。这些方法只有在特殊场合适用。 2) 从头算(ab initio)方法：HF方法：即基于Hartree-Fock原理的方法 关键词：HF，RHF，UHF，ROHF说明：I)当关键词为HF时，会自动根据自旋多重度选择 RHF还是UHF；Ii)ROHF为限制性开壳层HF方法，与UHF区别在 此时除了成单电子外，其余的和电子仍配对，通常该方法得到的能量要较UHF略高。Iii)HF方法可以看作是最低级的从头算方法，该方 法除了

25、在构型优化时有使用外，不适合计算能量。密度泛函方法(DFT)：基于电荷密度自洽的方法关键词：B3LYP等根据所采用的相关和交换泛函，可以选择不同的DFT方法，具体参见g98的帮助文件。其中B3LYP方法是使用最为广泛的DFT方法，由于DFT方法考虑了电子之间的相关作用，因此得到的能量要较HF来得精确，它是目前最常用的量子化学计算方法。其缺点是，构型优化得到的键长通常偏长。MPn方法：关键词：MP2,MP3, MP4, MP5说明：I)这些方法在HF基础上，进一步根据MP微扰理论考虑电 子相关作用，微扰项截至到二阶则为MP2，截至到三阶 则为MP3，其它类推，理论上考虑的微扰项越多，得到 的能量

26、越精确，但将大大增加计算量，而且通常也无此 必要，多数场合选取MP2即可。说明：Ii)对于该类方法，硬盘和内存通常开销较大，应考虑 具体的硬件考虑之，其中对于硬盘空间的设置见文件 Default.Rou内容，另外，必须注意到由于受到操作系 统的限制，中间文件似乎不能超过2GB，此时需设置 多个中间文件，具体见g98说明；耦合簇(Coupled Cluster)方法：关键词：CCD, CCSD(T)说明：I)该类方法与MPn方法一样，也是属于较高精度的计算 方法，其中CCD方法，只考虑了双取代，CCSD则在 CCD基础上进一步考虑了单取代； Ii)与MPn方法类似，该类方法计算量较大，通常只适用

27、 小体系。 其它后自洽场方法：组态相互作用(Configuration interaction)方法：关键词：CID, CISD, QCISD, CASSCF等。此外，g98还提供一些高精度的组合计算方法，如G1、G2等。举例：比较不同方法计算得到的H2基态能量，假使HH键长 为0.7A，采用STO-3G基组基态的电子态为：1gHF能量，单位a.u.HF= -1.117349; B3LYP=-1.1647796; MP2=-1.129582MP3=-1.1339601; MP4D=-1.1355271; MP4DQ=-1.13547MP4SDQ=-1.13547; MP4SDTQ=-1.135

28、47; 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：为最小基组，每个

29、原子轨道用三个高斯函数 (GF)来描述，原子轨道数即为基函数数目。 如O：1s2s2p，原子轨道数为1+1+3=5 GF数目为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：1s 2s 2p 3s 3p 内层1s, 2s和2p共有1+1+3=5个AO和

30、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：1s 2s 2p 3s 3p 内层1s, 2s和2p共有1+1+3=5个AO和6*5=30 个GF，价

31、层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)函数的使用： 在实际计算中，有时需在上述标准基组的基础上，添加一个或多个极化函数，极化函数是指具有比原子价轨道更高角量子极化函数是指具有比原子价轨道更高角量子数的高斯函数数的高斯函数。例如H的价轨道为1s，则其极化函数为p型GF，同样对C、O等价层为p轨道的原子，它们的极化函数应为d型或f型轨道，类似地，对于过渡金属原

32、子的极化函数为f型轨道。 极化函数的使用目的在于使原子价轨道在空间取向上变得使原子价轨道在空间取向上变得更更“柔软柔软”，从而使之易于与其它原子的轨道成键：例如对于羰基基团中的C，O原子，它们极化轨道(即d轨道)之间形成的d轨道，使得C原子的p轨道朝O方向极化，O的p轨道向C方向极化，从而增强了C与O之间的作用。由于极化轨道的使用增强价轨道的空间柔软程度，因此其使用场合主要在于环状化合物环状化合物，例如有机环状化合物以及含有桥联结构的金属簇合物等。使用方法： 在标准基组后加“*”，“*”，(d)，(d,p)，(3df,3pd)等例：6-31G*, 6-31G(df,pd)说明：当只添加一个一个

33、极化函数时，可使用下述表示：“*”等价于(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轨道; 对

34、O原子, 则在6-31G基组上添加2个d轨道和 1个f轨道.对于不同的基组,可添加的极化函数数目是不同的, 具体可参照程序手册:iv) 弥散函数(diffuse)的使用: 对于带有较多电荷带有较多电荷的体系,采用标准的基组来描述是不够的, 此时需添加弥散函数,以增加价轨道在空间上的分布范围增加价轨道在空间上的分布范围, 即: 极化函数用于改进价轨道的角度分布, 弥散函数则用于改 进价轨道的径向分布. 弥散函数是指具有较小轨道指数的高斯函数, 其表示方法是 在标准基组后加上“+”或“+”, 如6-31+G, 6-31+G*等, 其中 第一个“+”表示对非H原子添加弥散函数, 第二个“+”则对H原

35、 子添加弥散函数. 弥散函数主要用于带有电荷的体系带有电荷的体系(包括离子)以及弱作用体系弱作用体系.2) 赝势基组： 关键词：lanl1mb,lanl2mb,lanl1dz,lanl2dz, cep-4g,cep-31g,cep-121g,sdd等 说明： 通常对于重元素采使用赝势基组，其中，对于过渡金属原子 一般使用lanl系列的基组，对于主族元素采用cep系列的基组， sdd基组相对较少使用； 赝势基组的适用范围要较全电子基组来得广，一般除了稀土 和放射性元素外，均可使用赝势基组； 赝势基组是将内层原子的影响用势函数来代替，因此，在使在使 用时必须注意根据所采用的赝势基组，检查体系的总电

36、子数用时必须注意根据所采用的赝势基组，检查体系的总电子数 是否正确是否正确！ 对lanl1系列的基组，只考虑价层电子；如V原子为5. 对lanl2系列的基组，除了价层电子外，还考虑了次外 层的电子(对部分主族元素例外)；如V原子为5+8=13. 对其它赝势基组，请在使用时仔细核查。 对于使用较多的lanl系列基组，需注意： 对于HNe范围的原子，lanl1mb和lan2mb基组等价于 sto-3g最小基组；lanl1dz和lanl2dz基组等价于d95基组; 对于含金属金属键体系一般使用lanl2系列，而不用 lanl1系列； mb结尾含义为minimal basis，dz结尾含义为doubl

37、e -zeta，故后者基函数数目要较前者多。3) 混合基组的使用: 有时, 因需要, 对同一化合物中的不同原子采用不同的基组:受基组限制: 例如, 对于Mo(CO)2+化合物, C和O原子可采用6-31+G基组, 但 Mo原子不在6-31G基组的使用范围之内. 受体系大小限制: 当体系较大时, 不可能对所有原子均采用较大的基组, 此时可 对其中的局部原子采用精度较高的基组来描述,而对其余原子 采用小基组来描述, 例如:混合基组输入方法: 关键词: GEN和PSEUDO=READ(对赝势基组对赝势基组)%mem=32mb%chk=hf#p b3lyp/ gen pop=full scfcyc=5

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

39、03180.0r1=2.228r2=3.329Mo 0Lanl2mb*C O 06-31G*Mo 0Lanl2mb使用赝势基组时,应注意不能漏应注意不能漏!基组定义部分赝势定义部分,注意不能漏, 另外其后没有*号赝势部分和基组部分之间有一空行(3) G98输出结果: 下面以H2O能量计算的输出结果为例:L1输出G98版权信息%部分内容关键词部分计算所需调用的模块Title内容电荷, S以及分子构型各原子的直角坐标距离矩阵需核查需核查!注意注意:计算结果是以该坐标系为准计算结果是以该坐标系为准需核查需核查!MO初始猜测自洽场迭代求解部分自洽场迭代求解部分收敛指标收敛指标迭代次数迭代次数(通常接近

40、通常接近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最大的允许迭代次数为?.

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

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

43、-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.此时，因带较多负电荷，故弥散函数的影响较极化函数来得大。体系多重度的选取：对于复杂体系，特别是含有过渡金属原子的体系，多

44、重度的选取凭经验不易得到，此时，可通过比较不同多重度时体系的能量来确定基态能量。以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

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

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

47、虑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,3O 0. 0.0 0.0O 1.208 0.0 0.0 为了便于计算结果分析，采用限制性开壳层计算方法，及将和电子强制配对；输出所有分子轨道首先确定计算结果中各原子所处的位置：21ozX or yz+12因1s为内层电子，故O与O之间成键情况可不考虑No.1No.2 与上一MO类似，也以1s成分为主，处在能量较低的位置z+z+12No

48、.3, 4分别为O2s轨道之间形成的成键(SGG)和反键(SGU)轨道：+z+-No.5为O 2pz轨道形成的成键轨道：z+zNo.6和7轨道分别为O 2px和2py形成的成键轨道：+No.8和9轨道分别为O 2px和2py形成的反键轨道：z+No.10轨道为O 2pz形成的反键轨道：z+根据上述分析结果以及各轨道的能级，可以得到O2的MO分布为：2s*2s2pz2px2py*2px*2py*2pz电荷分布： 在能量计算结果中，除了能量以及上述的MO成分外，其它有用 的信息包括原子电荷和键级，对于原子电荷：以H2O为例(B3LYP/STO-3G)：Mulliken电荷分布说明： 电荷的绝对值是

49、没有意义的，其数值受到所用方法, 尤其是 基组的影响较大： 以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 因此，计算得到的电荷只具有相对意义，即只能考察其相对数值或变化趋势。当然，当考察变化趋势时，如考察H2X(X=O,S,Se)系列化合物中，X为不同原子时，对H原子电荷的影响时，必须使用相同的计算方法和基组。键级： 键级的大小反映了某根键的强度，同样地，键级的数值也受到 计算方法和基组的影响：O-H键的Mulliken键级以

50、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 由上结果可见，键级数值与所用的基组和计算方法是有联系的，但其影响不象电荷那么显著。此外，必须注意到，弥散基组对键级的数值影响较大。 因此，在实际应用中，对于电荷和键级，我们更关心的是它们的变化情况而不是它们的绝对数值。此外，必须强调的是，上述结果是采用Mulliken布居分析的方法，若采用其它布居分析方法，结果将有所不同。自然键轨道(NBO)分析方法： 我们前面所讲分子轨道均为正则

51、分子轨道，它是未经定域化处 理，这导致计算结果与我们通常的成键概念有所不同。例如 在乙烯分子中，碳碳之间为双键，但在正则MO中，反映C与C 之间成键作用的MO可能有多个，因此根据正则MO的结果，我 们无法断定CC是单键还是双键。此时，通过对正则MO的定 域化处理，可以得到通常意义上的成键图像。 正则MO的定域化处理方法较多，其中较为常用的是NBO方法， 其使用方法是在输入文件中添加关键词： POP=NBO以乙烯分子为例：%mem=32mb#p b3lyp/3-21g pop=nbo The NBO analysis of ethylene 0,1CC,1,CCH,1,CH,2,HCCH,1,C

52、H,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)： 1s2s3s 1px2px 1py2

53、py 1pz2pz故由上述轨道组成可知，该轨道为C 2px之间形成的键对于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#p b3lyp/6-31g* The energy of CO 0,1

54、OC,1,1.138计算结果为: -113.30945a.u.b. 分别计算单个C原子的能量, 此时需将O原子造成的BSSE 考虑在内:%mem=32mb#p b3lyp/6-31g* massage The atomic energy of C in CO 0,3OC,1,1.1381 nuc 0.0 计算结果为-37.84664a.u.c. 类似的, 计算CO中O原子的能量:%mem=32mb#p b3lyp/6-31g* massage The atomic energy of O in CO 0,3OC,1,1.1382 nuc 0.0 计算结果为: -75.06280a.u.则CO中

55、C-O键能为: -37.84664-75.06280+113.30945) =0.40001a.u.=1050.2kJ/mol (exp. 1070.3kJ/mol)2. 分别从正则和定域两种角度分析乙炔的成键情况;HCCHX12345%mem=32mb#p b3lyp/6-31g* pop(full,nbo) C2H2 0,1XC 1 1.0H 2 r1 1 90.0C 2 r2 1 90.0 3 180.0H 2 r3 1 90.0 3 180.0r1=1.0666r2=1.205r3=2.2716 r1r2r1正则MO:NBO: (C-H)(C-C)1(C-C)2(C-C)2(C-H)2

56、*(C-H)1(C-C)2(C-C)2*(C-C) 1*(C-C)b. 构型优化构型优化 化合物构型优化是G98的常用功能之一,构型优化过程是建立在能量计算基础之上的. 对于通常意义上的稳定构型是指具有最低能量的构型,即是在势能面上的能量最低点(极小点). 此外G98也能对非基态构型进行优化, 例如激发态构型以及过渡态构型等. 246810-138-136-134-132-130-128-126-124-122-120-118-116246810Energy初始构型能量计算计算能量梯度各原子受力和位移大小是否满足标准根据受力情况得到新的构型否输出最终结果是因此与单点能计算相比,构型优化只多调用

57、了计算能量梯度的模块L701, L702, L703等(1) 构型优化过程说明:(2) 构型优化的输入: 在单点能计算的输入基础上,添加关键词 OPT即可.%mem=32mb%chk=H2O#P B3LYP/6-311G OPT POP=FULL The geometry optimization of water 0,1OH,1,1.0H,1,1.0,2,105.4 说明:在默认情况下, 构型优化过程中将化合物的自由度全部放开, 若 要对化合物的局部构型进行优化, 此时分子构型需采用内坐标 方式输入, 同时关键词OPT需选取z-matrix选项, 例如:%mem=32mb%chk=H2O#P

58、 B3LYP/6-311G OPT=z-matrix POP=FULL The geometry optimization of water 0,1OH,1,1.0H,1,1.0,2,aa=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,1OH,1,rH,1,r,2,a Varia

59、bles a=105.4 Constants r=1.0表明该部分参数在构型优化过程中为变量表明该部分参数在构型优化过程中为常量OPT关键词的其它辅助选项: 该关键词选项较多,具体请参考G98的帮助文件,其中较为常用 的参数有: MaxCycle= 定义构型优化的最大轮数, 例如: opt(maxcycle=30) TS 优化对应于过渡态的构型, 例如opt(ts) QST2 根据用户给出的反应物和生成物的构型, 进行反应过渡 态构型的优化 QST3 与QST2类似, 只是用户需另外单独给出反应过渡态的 初始构型(3) 构型优化结果的输出: 以H2O为例: (B3LYP/6-311G)与单点能

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