H3反应势能面的构建——计算化学实验设计*

袁汝明 傅钢 韩国彬

(厦门大学化学化工学院 福建厦门 361005)

教学内容包含基本原理的讲解，Gaussian03W、Gaussview5.0、Origin等软件的使用，实际的上机操作，以及数据处理和分析，问题思考等。

1 基本原理

简单地讲，势能面就是把分子的能量表示为几何坐标的函数。值得一提的是，根据量子力学原理, 势能面并不是理所当然的, 而是建立在波恩-奥普海默(Born-Oppenheimer)近似的基础上，即将原子核运动与电子运动分开处理。对于N原子反应体系，体系的能量E是3N-6(非线形分子)或3N-5个(线形分子)内坐标的函数，可写为[2]：

E=f(R1,R2,…,R3N-6 or 3N-5)

(1)

若用几何图形表示该函数，其图像为3N-6或3N-5维空间的超曲面，即势能面[2]。对于3原子相互作用体系，变量数为3N-6=3，通常用两个核间距RAB，RBC以及其夹角θ作为变量，此时能量函数可表示为：

E=f(RAB,RBC,θ)

(2)

由式(2)得到的势能面为四维空间曲面，故难以直观地画出其图形。如果固定其中一个变量(常固定夹角θ)，此时体系只剩两个独立变量RAB和RBC，能量函数可简化为：

E=f(RAB,RBC)

(3)

以E为纵坐标，RAB、RBC为平面上的两个坐标，则得到一个三维曲面。这种曲面很像起伏的山峦, 也有山峰、山谷和洼地。分子体系的势能面可以决定分子的诸多性质。分子的稳态，如反应物和产物，对应于势能面上的极小点；而反应过渡态对应于连接极小点的山谷中的最高点，即势能面上的一级鞍点。当过渡态确定后, 根据它和稳态的能量差以及断面的形状可以计算反应速率。还可以将极小点附近的势能面简化为二维，在谐振子近似下计算体系的振动频率和零点能。

势能面函数一般无解析的形式, 但可通过量子化学方法计算出给定结构的能量，并通过相应的软件(如Gaussview或Origin等)绘制出来。常见的计算方法有半经验方法(AM1,PM3), HF方法，密度泛函方法(如B3LYP)，后自洽场方法(如CCSD，MPn系列，QCISD等)。需要注意的是，各种计算方法实际上都或多或少引入了近似，因此它们所描述的势能面也有一定的差异。

2 软硬件设备

计算软件：Gaussian03W(并行版)[3]，收费软件，厦门大学授权。

作图软件：Gaussview5.0[4]， 收费软件，厦门大学授权；Origin6.0以上版本或Excel。

计算硬件：PC服务器(本文使用的是双路四核xeon5345服务器)。

3 上机部分

3.1 计算内容

采用Gaussian03W中的QCISD(T)方法，结合6-311++G(2d,2p)基组，计算H3反应势能面。

① 共线碰撞的势能面(θ=180°)，RAB和RBC的初始长度设为0.58Å，步数为20，步长为0.09Å，共计算441个点(21×21)。输入文件(通常以gjf或com为后缀名)和说明如下[5]：

%nproc=8 //并行CPU数

%chk=h3.chk //checkpoint文件

# QCISD(T)/6-311++G(2d,2p) scan //任务类型

//空白行

H+H2 -> H2 + H PES Scan //文件说明项

//空白行

0 2 //电荷 自旋多重度

H //H3体系内坐标

X 1 1.0 //赝原子，常用于构建共线体系的内坐标

H 1 RAB 2 90.0

H 1 RBC 2 90.0 3 180.0

RAB 0.58 20 0.09 //键参数 初始键长 扫描步数 扫描步长

RBC 0.58 20 0.09

② 计算θ=180°，扫描RAB=RBC时的能量变化，初始长度设为0.88Å，步数为90，步长为0.002Å。输入文件略。

③ 改变①中的键角，分别设θ=150°、120°、90°，可和①同时进行。建议可将每个实验组分成4个小组，考察不同键角下势能面变化(输入文件略)。

④ 计算反应物(产物)能量，对比不同计算方法预测的反应能垒Eb。

3.2 结果输出

用文本编辑器打开Gaussian输出程序(通常以log或out为后缀名)，就可得到势能面扫描的结果(θ=180°)，截取其中的部分数据列于表1。

表1 势能面扫描的数据结果(部分)*

*其中体系能量单位为kJ·mol-1。

值得一提的是，采用QCISD(T)方法计算时，Gaussian03W还可以同时输出HF(即SCF)，MP2，MP3和QCISD的能量。因此，我们还可以方便地比较不同计算方法描述的势能面之间的差异。

3.3 数据处理

① 在GaussView5.0中，用File→Open命令打开输出文件，用Result→Scan命令打开一个绘有3D势能面的新窗口(图1)。在GaussView5.0下，按鼠标左键可以旋转图形，可以从各个角度观察H3势能面；按Shift+鼠标左键可平移图形；滚动鼠标滚轮可以缩放图形，方便观察特定区域。

图1 H3体系(H+H2→H2+H(θ=180°))的势能面

图2 H3体系(θ=180°)的势能等高图

② 为了得到更平滑的势能面，可用Origin程序中非线性内插的方法。将数据增加至50×50，绘制不同H3体系的势能等高图(图2)。由图2可见，反应物和产物分别位于图中的右下角和左上角；连接反应物, 过渡态和产物的山谷就是反应路径(虚线所示)；反应路径中的最高点为反应的过渡态。

③ 根据势能面的对称性，对过渡态附近的区域进行重点构筑，见计算内容②。根据定义，反应路径上的最高点能量和反应物能量之差，即反应能垒Eb，列入表2，并对比不同θ下能垒的高低。

表2 不同θ下反应能垒

由表2可见，θ越小，活化络合物中RAB(BC)的键长越长，H3反应所需克服的反应能垒也越高。因此，H原子进攻的最有利方式是共线碰撞。

3.4 问题思考

① 势能面扫描得到的活化能垒是否可以直接和实验活化能进行比较？(提示：还需考虑零点校正和温度校正。)

② 不同理论方法(如HF，MP2，QCISD和QCISD(T)等)计算得到的活化能垒是否相同？(提示：不同。各种计算方法均引入了近似，通常认为后自洽场方法所描述的势能面较为合理。)

③ 研究反应的势能面是计算化学的核心问题。能否结合近年来的研究进展，举一个或几个例子谈谈势能面在预测反应机理中的作用，以小论文的形式提交。

4 时间安排

本实验是为大三下学期或大四上学期学生所设计的物理化学探索性实验，需用时两天。具体的实验安排见表3。在教学过程中，我们对于各种计算方法的细节不作过多的讲解，而注重培养学生的兴趣，即如何应用成熟的计算化学软件去探索一个具体的反应过程。

表3 实验安排

通过上机练习和数据分析，可以得到不同θ角的H3反应势能面，并可通过Gausview程序进行360°全方位观察，有利于学生对过渡态理论的理解和掌握。此外，学生还可以熟悉各种计算相关的应用程序，如量子化学程序包(Gaussian03)，图形界面分子构建工具(GaussView5.0)以及数据处理程序(Origin)等。这种密切结合物理化学教学的计算机实验设计，不仅可以将前沿的科研手段融入教学中，提高教学内容的创新性，而且有助于学生将物理化学的基本概念运用到科学研究中。

参 考 文 献

[1] 傅献彩，沈文霞，姚天扬.物理化学.第4版.北京：高等教育出版社，1990

[2] 夏少武.活化能及其计算.北京：高等教育出版社，1993

[3] Frisch M J，Trucks G W，Schlegel H B，etal.Gaussian03W,Rev.E01.Wallingford CT：Gaussian Inc，2004

[4] Gaussview 5.0.Wallingford CT：Gaussian Inc，2008

[5] Frisch ☞，Frisch M J，Trucks G W.Gaussian03 User′s Reference.Wallingford CT：Gaussian Inc，2003