《量子化学计算方法》PPT课件

《《量子化学计算方法》PPT课件》由会员分享，可在线阅读，更多相关《《量子化学计算方法》PPT课件（213页珍藏版）》请在装配图网上搜索。

1、量子化学计算方法,章永凡 福州大学化学系 2009年2月,课程主要内容,一、有限尺度体系(分子、团簇等)电子结构计算 方法G03程序的使用 二、无限周期体系(一维链状化合物、二维层状 化合物或固体表面、三维固体体相)电子结 构计算方法VASP/CASTEP程序的使用 主要参考资料: G03用户手册或G03的帮助文件 相应网址: VASP程序用户手册,简 介,量子化学软件目的在于将量子化学复杂计算过程程序化，从而便于人们的使用、提高计算效率并具有较强的普适性。 绝多数量子化学程序是采用Fortran语言编写的(Fortran 77或Fortran 90)，通常由上万行语句组成。,软件分类,计算原

2、理,基于从头算或第一性原理方法 (ab initio/first principles) Gaussian、ADF、Dalton、Gamess、Crystal、VASP、Wien、Dmol等,基于半经验或分子力学方法 MOPAC、EHMO、NNEW3等,研究对象,有限尺度体系(分子、簇合物等) Gaussian、ADF、Dalton、Gamess、MOPAC、EHMO等,无限周期重复体系(晶体、固体表面、链状聚合物等) Crystal、NNEW3、VASP、Wien等,本研究室目前常用的量化软件： Gaussian 98/03: 由Pople等人编写，经过几十年的发展和完善， 该软件已成为国际

3、上公认的、计算结果具有较高 可靠性的量子化学软件，它包含从头算、半经验 以及分子力学等多种方法，可适用于不同尺度的 有限体系，除了部分稀土和放射性元素外，它可 处理周期表中其它元素形成的各种化合物； Crystal 98/03： 该软件由意大利都灵大学理论化学研究所开发， 采用基于原子轨道线性组合的从头算方法来研究 固体及表面的电子结构； VASP: 该软件由奥地利维也纳大学开发，采用基于平面波 基组的密度泛函理论来研究固体及表面的构型以及 动力学过程； CASTEP： MS软件模块之一，与VASP程序类似； Dmol: MS软件模块之一，主要用于有限尺度体系电子结 构研究；,采用理论方法要解

4、决的问题,当前的研究状况，包括实验和理论研究现状、已解决和尚未解决的问题,计算过程,化合物构型的确定，具体途径包括：利用实验测定结果、或者采用软件进行构造等,根据现有的计算条件、模型的大 小以及所要解决的问题，选择可 行的计算方法和相应程序,对计算结果进行加工和提取有用 的信息，一般包括构型描述、 能量分析、轨道组成、电荷和成 键分析等，并与实验结果比较,计算模型和方法的选取是保证计算结果可靠性的关键， 理想的情况是：1.所选取的计算模型与实际情形一致；2.采用 高级别的计算方法。但是，由于受到计算软硬件的限制，在多数 情况下，很难同时做到上述两点要求，实际操作中，当计算模 型较大时，只能选择

5、精确度较低的计算方法，只有对较小的模 型才能选取高级的计算方法。 因此，当确定了一种计算模型和方法后，最好对其进行验证， 以保证计算结果的可靠性。假设当前的研究对象是化合物A， 可通过下列途径进行验证： 1. 与A化合物现有实验结果之间的比较； 2. 若无实验方面的报道，可对与A类似的化合物B进行研究，此 时以B的实验结果作为参照； 3. 当上述方法行不通时，可以采用较大模型和较为高级的计算 方法得到的计算结果作为参照，该方法主要用于系列化合物 的研究：如对A1, A2, A3，先用大模型和基组对A1进行研究， 然后以该结果为参照，确定计算量适中的模型和方法并应用 于A1,A2,A3。,Gau

6、ssian03程序的使用,G03的安装和运行； G03的功能和程序结构； 输入文件的编写与主要功能的使用； 补充说明；,G03程序的安装和运行,1. G03程序的安装： (1).确定运行平台：Windows或Linux？ (2).对Windows平台： 直接运行setup.exe，其余步骤按提示操作即可；也可将其它机 器上将已安装好的G03直接拷贝到本机，但需设置运行环境。 对Linux平台： a.若G03是经过压缩过的(文件结尾为gz)，用gunzip命令解压： 例如： gunzip g03.linux.tar.gz b.若G03是打包的(文件结尾为tar)，用tar命令将其释放： 例如：

7、tar xvf g03.linux.tar ab两步合成一步方法： tar zxvf g03.linux.tar.gz,c.设置环境变量，以c shell为例，在用户根目录下的.cshrc文件 添加下列内容：(也可在执行g03前逐条运行) setenv g03root /home/$USER(设置g03所在目录，根据实际情况修改) source $g03root/g03/bsd/g03.login(激活g03运行时所需环境变量) setenv GAUSS_SCRDIR /home/$USER/g03_tmp(设置临时目录) d.运行bsd/install，自动配置并行计算环境 注：对Linux

8、平台，运行g03时，需注意权限问题，可用chmod 命令更改权限，将所安装的g03对所有用户开放。,2.G03程序的运行： (1).对Windows平台： a.对于刚安装好的g03，先检查环境设置情况：,左侧至上而下依次为：默认的文本编辑器；g03可执行文件所在目录； 计算中间结果存放目录；缺省的计算结果存储目录；缺省的输入文件 所在目录；PDB分子构型浏览器； 右侧至上而下依次为：设置显示属性(如背景色等)；设置文本编辑器 属性；计算过程控制属性(尤其是批作业过程)；Default.Rou文件的编 辑(该文件内容为默认情况下，计算所花费的内存及硬盘大小),需设置正确， 否则运行将出错！,b.

9、编写或打开g03输入文件,点击RUN，并给定输出 文件名后开始运行,c.g03运行过程的控制：,最上行按钮的功能从左至右依次为： 开始运行g03；暂停进程；运行至下一模块(link)时暂停进程； 重新启动进程；清除进程(停止运算)；编辑批作业；运行完 当前任务后，暂停批作业；停止批作业的运算；观看计算结 果；打开文本编辑器；,不要随意点击！,交换机,计 算 节 点,计 算 节 点,计 算 节 点,网 关,用户终端,基于Linux系统的计算拓扑结构,内部网(高速),外部网(普通),Window系统,Linux系统,(2). Linux平台：,说明： 网关作用类似于防火墙，用于保证内部网的安全和稳

10、定， 作为网关的计算机通常配有2个网卡，分别用于外部网和 内部网的连接。 本实验室网关机子IP地址: 219.229.140.103(非固定IP，可能发生改变) 计算作业提交过程： a. 用户登录网关通过SSH远程登录软件实现 SSH软件(SSHSecureShellClient-3.2.9.exe)可从网络上免费 下载，安装过程与通常软件安装类似。安装完毕后，设置 网关外部网的IP地址以及账号名即可使用。,点击Profiles设置IP地址及用户名,b. 从网关登录到计算节点采用telnet命令实现 例如：telnet 134.14.83.5,c. Linux常用命令: (1)ls显示文件清单

11、，相当于DOS下的dir命令：,文件属性,所属用户,大小,创建时间,注：Linux系统下字符是大小写区分的,(2)cp复制文件命令，相当于DOS下的copy命令: cp -rf,(3)mkdir创建目录，相当于DOS下的md命令：,(4)rm删除文件或目录，相当于DOS下的del命令：,删除文件,删除目录,(5)top显示当前进程和CPU以及内存使用情况,(6)kill终止某个进程，格式为：kill PID号 (PID号由top命令可得，受权限限制),(7)renice调整某个进程优先级，格式为：renice 级别 PID号 (级别为019整数，数值越大优先级越低) renice 19 79,

12、(8) cat显示文件内容，格式为：cat 文件名 (9) grep一般用于从某个或多个文件中搜索某串字符， 格式为：grep “字符串” 文件名 例：grep “F=” vasp.out (10)scp用于网关与内部网内各计算节点或外部网络之间 的文件传输 格式为： 从其它到网关： scp 文件 用户名网关IP:目录 例：scp vasp.out zyfzyf-2400:/trans 从网关到其它： scp 用户名网关IP:目录/文件名 目标目录 例：scp zyfzyf-2400:/trans/vasp.out . 该命令也可用于同一台计算机不同用户之间的文件 传输,(11) vi文本编辑

13、命令 该命令常用但较为复杂，它有2种模式：命令模式和插入 模式，二者之间关系为： i Esc command mode insert mode command mode 在命令模式下，可实现以下功能及其对应按键： delete a character: x delete a line: dd search a string : /(向后) ?(向前) save the change: :w save the change and quit: :wq quite without saving : :q! page down: Ctrl+d page up: Ctrl+u go to file e

14、nd: shift+g go to n line: :n,(12) tar文件打包命令(适用文件扩展名为tar) 该命令用于多个文件/目录的打包或解包，常用格式有： 文件打包：tar cvf 要打包成的文件名 要打包的文件 例：tar cvf model.tar * 生成 model.tar文件 解包： tar xvf 要解包的文件名 例：tar xvf model.tar (13) gzip/gunzip文件压缩/解压命令(文件扩展名为gz) 例：gzip model.tar 生成 model.tar.gz 文件 gunzip model.tar.gz 生产model.tar文件 (14)

15、du察看当前目录所占硬盘空间大小(类似命令df) 例：du h df h (察看硬盘各分区大小),(15) rsh或ssh用于从某个节点登录到其它节点 例：rsh c0102 登录到c0102节点上，为当前用户名 rsh zyfc0102 与上相同，但用户名为zyf ssh c0102 ssh zyfc0102 根据rsh或ssh服务的具体设置来确定是否需要提供密码 (16) su从当前用户转变为超级用户或其它用户 例：su 转变为超级用户 su zyf 将用户转变为zyf用户 (17) ifconfig察看网络设置 (18) dmesg察看系统日志 (19) adduser ,passwd,

16、d. 运行g03过程： (1)编写输入文件： 用vi命令编写或在Windows下编写完毕后ftp至Linux系统； vi test.gjf (2)运行g03： g03输出文件名 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.1361895a.u.,(2).基组的选择,全电子基组,赝势基组,高斯型函数: g p ( p, r ) = Np exp ( - p r2 )

17、 原子轨道可表示为高斯型函数的线性组合： 1 S = p = 1 , L (d p g p( p , r) ),系数,指数,最终体系分子轨道为这些原子轨道的线性组合,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)基组的大小决定了基函数的数目，即体

18、系的原子轨道 数目，因此可从所选择的基组来推断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是相同的。 如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

19、 对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数为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

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

21、化函数应为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

22、,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轨道

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

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

25、基组，只考虑价层电子；如V原子为5. 对lanl2系列的基组，除了价层电子外，还考虑了次外 层的电子(对部分主族元素例外)；如V原子为5+8=13. 对其它赝势基组，请在使用时仔细核查。,对于使用较多的lanl系列基组，需注意： 对于HNe范围的原子，lanl1mb和lan2mb基组等价于 sto-3g最小基组；lanl1dz和lanl2dz基组等价于d95基组; 对于含金属金属键体系一般使用lanl2系列，而不用 lanl1系列； mb结尾含义为minimal basis，dz结尾含义为double -zeta，故后者基函数数目要较前者多。,3) 混合基组的使用: 有时, 因需要, 对同一化

26、合物中的不同原子采用不同的基组: 受基组限制: 例如, 对于Mo(CO)2+化合物, C和O原子可采用6-31+G基组, 但 Mo原子不在6-31G基组的使用范围之内. 受体系大小限制: 当体系较大时, 不可能对所有原子均采用较大的基组, 此时可 对其中的局部原子采用精度较高的基组来描述,而对其余原子 采用小基组来描述, 例如:,混合基组输入方法: 关键词: GEN和PSEUDO=READ(对赝势基组) %mem=32mb %chk=hf #p b3lyp/ gen pop=full scfcyc=500 Mo(CO)2+(Line) 1,6 Mo X11.0 C1r1290.0 O1r229

27、0.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 Mo(CO)2+(Line) 1,6 Mo X11.0 C1r1290.0 O1r2290.030.0 C1r1290.03180.0 O1r2290.03180.0 r1=2.228 r2

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

29、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次,

30、 通常对于含有过渡 金属的体系, SCF较难收敛, 迭代次数会多些.,分子轨道系数的输出: 关键词: pop=full 说明: 有时为了分析MO成分, 则需利用该关键词输出各MO的 成分, 以H2O为例:,轨道对称性,O表示占据轨道; V为空轨道,轨道能级,AO轨道,轨道系数,HOMO,LUMO,注: 输出的MO系数是未经归一化过的; 1s, 2s等符号并非真正意义上的AO, 它与计算所采用的基组 有关. 对于上例, 采用的是6-31G基组, 可知, 对于价层AO是 分裂为两组的, 故对于O原子, 输出的2s和3s实际上均为价 轨道, 其它类推. 当基函数较多时, 若只考察前线附近的轨道成分,

31、 可用关键 词: pop=regular, 此时只给出5个占据轨道和5个空轨道组成. 练习: 计算O2的基态能量, 并分析MO轨道组成 假设基组为6-31G*，O-O键长为 1.2145 A,(5) 能量计算的几个问题： 基组的选取： 理论上而言，计算中所选取的基组越大，计算结果越准确，但由 于受到硬件上的限制，需根据实际情况选择基组。此外，对于环状 分子或存在作用的体系，通常需要考虑使用极化基组；对于带有 较多电荷的体系，或考察弱相互作用的体系，通常需考虑使用弥 散基组； 以O2为例(O-O为1.208A)： 6-31G-150.26755a.u. 6-31G*-150.31998a.u.d

32、E=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. 此时，因带较多负电荷，故弥散函数的影响较极化函数来得大。,体系多重度的选取： 对于复杂体系，特别是含有过渡金属原子的体系，多重度的选取 凭经验不易得到，此时，可通过比较不同多

33、重度时体系的能量来 确定基态能量。 以O2为例(基组为6-31G*)： S=1-150.25734a.u. S=3-150.31998a.u. S=5-149.79200a.u. S=7-148.80784a.u. 由上结果可见，O2的基态应为三重态。 自旋污染判断： 对开壳层体系的计算，需判断自旋污染情况，以便保证计算结果的可靠性：,=s*(s+1)，例如对O2基态，有2个单电子，s=2*1/2=1，则 =1*(1+1)=2，对上例子s=3/2,体系电子态的选取： 多数情况下，体系的电子态由体系的点群和电子占据情况由程序自动确定，但在某些情况，特别是在团簇体系，采用这种默认方法得到的电子态可

34、能出现问题，其中包括：无法确定轨道的对称性或者电子填充顺序出现问题。此时可人为对电子的填充轨道进行调整。所采用的方法为Guess关键词中的Alter选项。此外，采用该方法也可实现对特定激发态进行研究。,#p b3lyp 6-31G opt NH2 0,2 N H 1 r H 1 r 2 a r 1.03 a 120.0,NH2自由基的基态,NH2基态的电子填充情况为：a12a12b22a12b1，即相应的电子态为2B1, 相应能量为-55.85291a.u.,如果我们想考察NH2的2A1激发态，则可对电子填充顺序进行人为调整。具体的办法是将态中的4号和5号轨道进行对调，使得最终电子填充顺序为：

35、a12a12b22b12a1，即体系电子态为2A1。,#p b3lyp 6-31G opt Guess=Alter NH2 0,2 N H 1 r H 1 r 2 a r 1.03 a 120.0 4 5,该行前面有一空行，对应于态电子填充的改变情况,最终计算结果：,键能的计算： 以O2为例，计算OO键能： 步骤如下： 首先确定O2基态的能量： 假设基组为6-31G*，则总能量为： -150.32005a.u.(基态为：3g) (O-O: 1.2145 A) 其次，计算O原子基态能量，结果为：-75.06062a.u. 则O2中OO键的键能为： 2-75.06062+ 150.32005=0.

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

37、 0.0 0.0 2,表示体系由2个碎片组成,整个体系的电荷和自旋多重度 碎片1的电荷和自旋多重度 碎片2的电荷和自旋多重度,该原子属于那个碎片,由此得到考虑了BSSE时的键能为： 0.19881-0.00436=0.19445a.u.=510.4KJ/mol。此外，进一步的修 正是在计算O2总能量时，需考虑零点能(ZPE)校正。,分子轨道成分的分析： 以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,为了便于计算结果分析，采用限制性

38、 开壳层计算方法，及将和电子强制 配对；,输出所有分子轨道,首先确定计算结果中各原子所处的位置：,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,

39、+,根据上述分析结果以及各轨道的能级，可以得到O2的MO 分布为：,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.7

40、950.398 因此，计算得到的电荷只具有相对意义，即只能考察其相对 数值或变化趋势。当然，当考察变化趋势时，如考察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 由上结果可见，键级数值与所用的基

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

42、MO的定域化处理，可以得到通常意义上的成键图像。 正则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计算结果的主要内

43、容：,自然布居分析结果,自然电荷,自然电子组态,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均

44、为弥散轨道，电子填充数较少,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键能(CO键长为1.138 )； 2. 分别从正则和定域两种角度分析乙炔的成键情况; (CC和CH键长分别为1.205和1.0666),正则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

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

46、32mb %chk=H2O #P B3LYP/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

47、 O H,1,1.0 H,1,1.0,2,a 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,表明该部分参数在构型优化过程中为变量,表明该部分参数在构

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

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

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

51、属于C3V群的NH3. 为了克服该困难, 可用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

52、, 1110等)。 (1) IR计算的输入： 关键词为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分析须经

53、过下述步骤： a. 首先对H2O的构型进行优化，此时输入文件为：,%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含义是

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

55、y optimization and the 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活性(强度)；极化

56、率,简正坐标，通过对简正坐标的分析， 可对振动类型进行归属,简正坐标分析说明：,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计算得到的频率要校实验结果来得大些, 若要得到较为 准确的数据, 需用校正因子校正, 该校正因子的数值与所用的计 算方法和基组均有关, 具体参考G03手册; 现有多种程序可用于将IR计算结果图示化, 例如Gaussian View 等, 从而便于简正坐标的分析;

57、 当所优化的构型并非对应于能量极小点时, 将得到较为明显的为 负值的频率, 即虚频, 此时可以通过对简正坐标的分析, 推测稳定 构型, 从而消除虚频的出现.,有一数值较大的虚频, 表面平面性NH3并非对应于势能面最低点,例, 对平面型NH3的freq计算结果为:,N,H,H,H,x,y,练习: 1. 分别对重迭和交错式乙烷的构型进行优化, 并进行IR计算, 确定它们构型的相对稳定性。 2. 确定NH3由三角锥型转化为平面形所需克服的能垒。,对重迭式乙烷：,1,2,3,4,5,6,7,8,1(2),3(7),4(6),5(8),方法一：构型优化和IR计算一步完成,%mem=64mb %chk=c

58、2h6 #p b3lyp/6-31g opt freq 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,1 2 3 4 5 6 7 8,计算方法,所用基组,1,2,3,5,6,7,8,4,对同一化合物, 内坐标的表示方式可有多种, 但 为了能够更好的反映出体系的对称性, 需要选取 一种较好的表达方式, 通常内坐标表示中的前三 个原子的选

59、取较为重要, 一般原则是: 1. 1号和2号原子尽量处在体系的对称轴上; 2. 若体系有对称面, 则该3个原子最好处在该对称 面上. 因此, 在书写内坐标前, 最好先对体系的对称元素 分布进行考察. 在左图中, 按照其编号顺序来描述乙烷的构型, 较难反映出其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 a

60、1 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,先构型优化,在构型优化的基础上进行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文件存放经优化过的构型和波函数 等信息，在作进一步的计算时，可充分 利用该文件的内容；,构型优化结果 (本例为自由度全放开

61、的优化, 因此变量较多, 为3N-6个, 但因受到对称性的限制, 有些变量的数值是相同的),若只对部分变量进行优化, 需采用opt=z-matrix选项 %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中的数值均认为是 常量,

62、在构型优化过程中是不变的, 因此 采用该选项时, 需注意那些参数要固定, 那些参数需放开优化,IR计算结果:,对构型优化的补充说明: 除了采用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

63、2 r2 1 a1 5 d1 r1=1.5481 r2=1.0956 a1=111.7037 d1 0.0 10 12.0,因在扫描过程中 分子的点群将改 变, 故需取消对称 性限制,依次为初始值, 扫描点数, 步长,由图可见, 乙烷交错式和重迭式之间存在约11kJ/mol的势垒,11kJ/mol,Scan的方法只适用于维数较低的情形, 有时为了寻找构型改变的主要影响因素时可用该方法, 此外也适用于单变量体系. 采用scan方法的缺点是在扫描某个或多个变量的同时，其它部分构型是固定的，未进行优化。因此，更为合理的方法是在扫描的同时对剩余部分的构型也进行优化，此时可以采用限制性优化方法，在G03，可通过opt关键词的选项ModRedundant来实现。,%mem=64mb #p b3lyp/6-31g opt=ModRedundant nosymm Scan of C2H6 0,1 C C 1 r1 H 1 r2 2