化学反应动力学是现代化学一个非常重要的分支，通过研究化学反应过程中的反应路径和化学结构变化等信息，人们不仅可以揭示化学反应的内在机理，还可以对化学反应进行调控，使之输出人们感兴趣的产物^[1]。化学反应动力学是化学、化工类学生非常重要的一门专业课程，要求学生在学过物理化学的基础上，从唯象规律开始，在基元反应的层次上深入学习化学反应动力学的理论、观点和方法，为以后工作和科学研究掌握正确的思维方法和手段奠定基础^[2]。对化学反应进行动力学研究时，总是从动态的观点出发，从宏观的、唯象的研究进而到在微观的分子水平上采用理论化学和计算化学的方法阐明和揭示化学反应^[3]。微观反应动力学从原子、分子的微观性质出发，分析分子间的运动及相互作用；在课程教学过程中引入了很多新的概念、理论和分析方法，具有理论性较强、抽象难懂等特点，容易使学生在学习过程中产生理解方面的困难。准经典轨线方法是一种以经典散射理论为基础的分子反应动力学计算方法^[4]，该方法将反应过程中的原子、分子的原子核视为经典力学的质点，求解这些质点在电子构成的势能场中运动，得到各个质点的运动轨线，通过对上万条乃至数十万条轨线进行统计平均，可以验证、预测实验结果。该方法计算量较小、计算结果可靠且能够提供清晰明了的反应过程图像。本文通过借助Matlab软件的图形用户界面(Graphical User Interface，简称GUI，又称图形用户接口)，结合准经典轨线方法的Fortran数值计算，将a + bc类型的基元反应过程进行仿真模拟。最终通过动画的形式将化学微观反应的具体过程展示出来，使得抽象、晦涩难懂的教学内容变得形象立体，方便直观地帮助学生理解化学反应的微观过程，深入理解反应机理。

从分子或原子的微观层次来看，化学反应可以视为具有一定量子态的反应物分子原子间的互相碰撞，进行原子重排，产生新的、处在特定量子态的产物分子的一个过程，在碰撞过程中反应体系的核构型不断发生改变。

在基元反应的碰撞过程中，分子原子间的相互作用包括了原子核与原子核之间、电子之间、电子与原子核之间的相互作用。为了简化计算，通常可采用核绝热近似(波恩-奥本海默近似)，即假定电子的运动与核的运动可以分离。这样就引入了势能面的概念^[5]，即原子核在电子组成的势场(势能面)中运动。碰撞过程中各个分子、原子的原子核在势能面上的运动遵循经典力学的Hamilton方程。通过求解Hamilton方程，可以得到各个原子核广义坐标和广义动量随时间的变化情况(称为“运动轨线”)，在此基础上可以进一步确定原子之间是否发生了重新组合，是否发生了化学反应，以及碰撞前后原子或分子能量状态是否发生了变化。通常利用Monte Carlo方法来选择初始条件，利用四阶龙格库塔加四阶预估校验(RKG4-AMH4)方法对轨线进行积分传播。每条轨线的初、末态采用量子力学进行描述，即利用振动、转动量子数来描述反应物和产物分子的内能态。对准经典轨线方法的基本原理和详细推导过程感兴趣的读者可以进一步参阅俞书勤^[4]编写的《微观化学反应》或文献^[6]。

基于准经典轨线的基本原理，人们用Fortran语言编写了很多计算程序包，用来计算三原子反应体系或多原子反应体系的动力学性质，其中应用比较广泛地有Truhlar及其合作者^[7]编写的“CLASTR”程序包，Hase等^[8]编写的“VENUS96”程序包，前者用来计算A + BC类型化学反应，后者可以计算更为复杂的多原子反应体系。基于CLASTR程序包，段志欣^[6]重新编写了新的准经典轨线计算程序，增加了很多新的计算模块，同时也能够输出每条轨线在传播过程的一些空间位置、动量、能量、时间等信息，这些信息为后续动画设计提供了基础数据。在程序中，准经典轨线的计算包含如下几个步骤^[6]：(1)读入各输入参数，如碰撞能、最大碰撞参数、碰撞壳半径、反应物分子的初始振动、转动量子数。这些参数在输入文件中可以设置、修改；(2)利用随机数，设置碰撞参数、反应物分子轴的初始空间方位等信息，得到轨线的初始条件；(3)对Hamilton方程进行积分，传播轨线；(4)对每一条轨线进行分析，确定产物，判断是否计算完毕所有轨线；(5)对计算的所有轨线进行统计分析，得到反应几率、积分/微分反应截面、反应速率常数、产物态分布等化学反应的相关信息。准经典轨线程序运行完毕后，我们可以得到反应轨线每一个时间步上三个原子的空间位置、动量等信息。每一条轨线的相关传播信息储藏在单独的文件中。

这些输出文件包含有Hamilton方程求解过程中每一个时间步的12个变量：(Q₁、Q₂、Q₃)为原子c相对于b为原点的坐标，(Q₄、Q₅、Q₆)是原子a相对于bc分子质心的坐标；(P₁、P₂、P₃)和(P₄、P₅、P₆)表示相应的共轭动量。利用这些信息把动画演示出来，还需要进一步进行信息提取和坐标变换，我们简要地描述一下坐标变换的基本原理，这些变换的具体实现将在Matlab中进行。首先取b原子位置为坐标原点，则各个原子的坐标矩阵满足：

(1)$B = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ 0 \end{array}} \right], \begin{array}{*{20}{c}} {}&{} \end{array}C = \left[ {\begin{array}{*{20}{c}} {{Q_1}} \\ {{Q_2}} \\ {{Q_3}} \end{array}} \right], \begin{array}{*{20}{c}} {}&{} \end{array}A = \left[ {\begin{array}{*{20}{c}} {{Q_4}} \\ {{Q_5}} \\ {{Q_6}} \end{array}} \right] + \frac{{{M_{\rm{a}}}}}{{{M_{\rm{b}}} + {M_{\rm{c}}}}}\left[ {\begin{array}{*{20}{c}} {{Q_1}} \\ {{Q_2}} \\ {{Q_3}} \end{array}} \right]$

三原子体系质心的坐标为：

(2)${\rm{O}} = \frac{{{M_{\rm{a}}}A + {M_{\rm{b}}}B + {M_{\rm{c}}}C}}{{{M_{\rm{a}}} + {M_{\rm{b}}} + {M_{\rm{c}}}}} = \frac{1}{{{M_{\rm{a}}} + {M_{\rm{b}}} + {M_{\rm{c}}}}}\left[ {\begin{array}{*{20}{c}} {{M_{\rm{c}}}{Q_1} + {M_{\rm{a}}}{Q_4} + \frac{{{M_{\rm{a}}}{M_{\rm{c}}}}}{{{M_{\rm{b}}} + {M_{\rm{c}}}}}{Q_1}} \\ {{M_{\rm{c}}}{Q_2} + {M_{\rm{a}}}{Q_5} + \frac{{{M_{\rm{a}}}{M_{\rm{c}}}}}{{{M_{\rm{b}}} + {M_{\rm{c}}}}}{Q_2}} \\ {{M_{\rm{c}}}{Q_3} + {M_{\rm{a}}}{Q_6} + \frac{{{M_{\rm{a}}}{M_{\rm{c}}}}}{{{M_{\rm{b}}} + {M_{\rm{c}}}}}{Q_3}} \end{array}} \right]$

动画演示时采用三原子体系的质心O作为坐标原点，因此还需要再进行一次坐标变换：

(3) $OA = A - O;OB = A - B;OC = C - O$

其中OA，OB，OC是以三原子体系质心为坐标原点的坐标系下，a、b、c三个原子的空间位置矩阵。

图形用户界面的制作主要包括界面设计和程序实现两个主要的步骤^[9]，首先需要构思出图形用户界面草图，制作静态的界面；之后要编写图形用户界面的动态程序功能。

Matlab的GUI功能可以开发出功能友好的图形界面，利用图形用户界面开发环境GUIDE可以通过简单的鼠标拖曳等操作设计出自己的GUI程序界面^[9]。图1显示了微观化学反应动画主界面的整体布局。主界面包含三个主要部分：左上侧的“选择数据”部分，左下侧“控制按钮”部分，以及右侧的“动画演示”部分。弹出式菜单“选择文件”，用于选择Fortran程序包计算得到的数据文件。a、b、c三个原子的质量和半径可以在主界面的6个编辑框内直接输入，也可以在M格式的文件中进行修改(修改后点击“刷新”按钮即可显示数据)。为了更好地观察动画演示过程，我们还设置了碰撞原子的颜色选择按钮、时间步数编辑框和观察视角选择框。时间步数默认20个时间步显示一帧动画，当碰撞过程进行比较快时，可以调节时间步的大小；观察视角框提供了6个不同的方位视角对碰撞过程进行观察。在左下侧，四个控制按钮“Play、Pause、Continue、Stop”可以对运行中的动画进行控制。

完成图形用户主界面的静态设置和制作后，保存运行该界面的文件，MatlabGUI会自动生成界面对应的初始化M格式文件。要实现图形用户界面的动态功能程序，需要进一步在M文件中编写程序的功能代码，即通过设置主界面上交互控件的回调函数来完成不同交互事件下后台程序完成的功能。回调函数编写的主要代码如下：

(1) Fortran计算结果文件的读取

function pb1_Callback(hObject, eventdata, handles)

global Y;

[yDate, yPath]=uigetfile({'*.*'});

yPath2 = strcat(yPath, yDate);

Y=load(yPath2);

set(handles.tx1, 'string', yPath2);

(2)原子质量、半径、颜色的选择

原子的质量、半径可以在M文件中直接写入，点击“刷新”按钮，即可载入文本框内。

function pushbutton12_Callback(hObject, eventdata, handles)

global Ma Mb Mc;

global Na Nb Nc;

global Nstep;

set(handles.edit1, 'string', Ma);

set(handles.edit2, 'string', Mb);

set(handles.edit3, 'string', Mc);

set(handles.edit4, 'string', Na);

set(handles.edit5, 'string', Nb);

set(handles.edit6, 'string', Nc);

set(handles.edit7, 'string', Nstep);

也可以在文本框内用键盘直接输入(以A原子的Ma为例)：

function edit1_Callback(hObject, eventdata, handles)

global Ma;

Ma=str2num(get(hObject, 'string'));

function edit1_CreateFcn(hObject, eventdata, handles)

global Ma;

set(hObject, 'string', Ma);

end

后面的“Mb、Mc、Na、Nb、Nc、Nstep”以及表示颜色的“yAco、yBco、yCco”可以设置类似回调函数。

(3)坐标变换、球体设置以及动画运行：

B=[0];

C=Y(:, 1:3);

M1=Mc/(Mb+Mc);

Ms=Ma+Mb+Mc;

for i=1:3

A(:, i)=M1*Y(:, i)+Y(:, i+3);

end

O=(Ma*A+Mb*B+Mc*C)/Ms;

Oa=A-O;

Ob=B-O;

Oc=C-O;     %坐标变换结束

[xa, ya, za]=sphere(30); % A原子的球体设置

xa=Na*xa;

ya=Na*ya;

za=Na*za;

B、C原子球体表示省略；

for i=1:3

T(:, i)=Y(:, 13); %时间提取

end

yPau=1;

ySt=1;     %动画演示控制变量

videoif=0;

videoif=get(handles.videout, 'value');

Nmax=length(oa);

VideoY=moviein(floor(Nmax/Nstep));

axes(handles.axes1);

videoi=0;

for i=1:Nstep:Nmax

videoi=videoi+1;

if ySt==0

break;

end

if yPau==0

while yPau==0

pause(0.2);

view(yView);

end

end

surf(za+oa(i, 3), xa+oa(i, 1), ya+oa(i, 2), 'EdgeColor', 'none', 'FaceColor', yACo);

hold on;

surf(zb+ob(i, 3), xb+ob(i, 1), yb+ob(i, 2), 'EdgeColor', 'none', 'FaceColor', yBCo);

surf(zc+oc(i, 3), xc+oc(i, 1), yc+oc(i, 2), 'EdgeColor', 'none', 'FaceColor', yCCo);

axis equal

light

axis([-8, 15, -5, 5, -5, 5]);

view(yView);     %设置yview数值可以从不同视角演示动画

axis off

title(T(i, 1));

if(videoif==1)

VideoY(:, videoi)=getframe;

end

hold off

pause(0.001)

end

fps=str2num(get(handles.fps, 'string'));

picqu=str2num(get(handles.picqu, 'string'));

movie(VideoY, 1);

if(videoif==1)

movie2avi(VideoY, 'out.avi', 'fps', fps, 'quality', picqu);

end

(4)动画运行的控制

Ⅰ.动画展示步数的选择，按积分时间步的倍数进行动画展示，步数越小，动画帧数越多。

function edit7_CreateFcn(hObject, eventdata, handles)

global Nstep;

set(hObject, 'string', Nstep);

end

Ⅱ.动画的运行、暂停、复位、停止利用两个时间参数来控制。yst为全局控制变量，yst = 1时，动画开始运行，yst = 0时动画stop；ypau为循环内变量，ypau = 0时动画暂停，ypau = 1时继续运行。

高放热反应F + HCl → Cl + HF属于典型的重+轻-重质量组合反应体系，它在理解卤化氢化学激光的发光机理方面扮演着非常重要的角色^[6]。我们利用改进后的CLASTR程序包对F + HCl反应体系进行了准经典轨线计算，轨线传播过程中的三个原子的空间位置、动量、传播时间等信息储存在单个文件中。为了更好地说明不同碰撞能下碰撞类型和反应机理，我们计算了83.64 kJ·mol⁻¹、209.1 kJ·mol⁻¹、1254.6 kJ·mol⁻¹三个碰撞能下的动力学信息，每个碰撞能下计算了10000条轨线。轨线开始时，F原子到HCl分子质心的距离设定为100 nm，时间步长设定为0.1 fs。

在Matlab平台上运行M格式文件，点击运行(Run)后出现动画主界面。点击“选择文件”按钮，选择F + HCl反应体系准经典轨线计算的输出结果文件，在a、b、c原子对应的空白框内，填入各个原子的相对原子质量和半径。查询“元素周期表”和“原子半径周期表”可得F、H、Cl原子的相对质量和半径，数据见表1。

为各原子选择相应的颜色(F原子为红色、H原子为绿色，Cl原子为蓝色)，随后点击“Play”按钮，就可以在右侧显示区观察反应过程中各原子、分子状态随时间演变的动画。运行界面如图2所示。

利用输出的视频文件，可以对化学反应的微观碰撞过程进行详细的分析，也可以揭示不同的反应机理。下面以不同碰撞能下的典型轨线为例，从动画运行过程截图来说明不同的碰撞过程和反应机理。

F + HCl碰撞可以生成新的产物分子HF，发生化学反应；也可以通过碰撞将能量传递给反应物分子HCl，改变HCl分子的振动、转动状态，但不生成新的分子；在碰撞能较高时(如1254.6 kJ·mol⁻¹)，还可以将HCl分子解离成原子状态。这三种不同的碰撞过程分别如图3、4、5所示。

图3描述了碰撞能为209.1 kJ·mol⁻¹时F + HCl → Cl + HF化学反应的一种微观过程，图3(a)到图3(g)按照时间顺序进行排列。从图3中可以看到，F原子飞向HCl分子，从H原子端“进攻”并“夺取”H原子，直接生成新的产物HCl分子，这种反应过程是A + BC类型非常典型的反应机制，被称为“夺取机制(stripping mechanism)”^[10]。图4显示的是碰撞能为83.64 kJ·mol⁻¹情况下的一种非反应性碰撞过程，碰撞中并没有生成新的产物分子，但存在能量的传递。如图4(d–h)所示,在碰撞前HCl分子处在振动、转动基态(v = 0, j = 0)，HCl分子没有转动运动，但在碰撞过程中，部分碰撞能配分给了HCl分子的内能，HCl分子被激发到了较高的振动、转动态，其振动、转动量子数变为(v’ ≥ 0, j’ ≥ 0)。在碰撞能比较高时，反应物HCl分子可以直接解离成原子状态，图5给出了在1254.6 kJ·mol⁻¹碰撞能下一个典型的碰撞解离过程，反应物HCl分子直接解离成H原子和Cl原子，碰撞过程可表示为F + HCl → F + H + Cl。

反应机理的探究是化学反应动力学研究的重要内容，利用碰撞过程的微观动画，可以比较方便地分析、研究a + bc类型的微观反应机理。F + HCl → Cl + HF反应存在着非常丰富的反应机理^[6]，它的很多动力学行为也可以利用反应机理来进行解释。上文图3中展示了生成HF分子的“夺取机制”，而图6展现了另外一种机制，即“碰撞迁移”机理。

图6(a‒c)显示，F原子首先靠近Cl原子端，HCl分子只有振动运动；但在F―Cl间距达到最小值时，此时进入势能面的强排斥区域，电子对HCl分子产生较强的力矩作用，使得H原子绕Cl原子发生快速转动，如图6(c‒f)所示。这表明F原子先和Cl原子发生碰撞，H原子快速迁移过来，最终形成HF分子。我们进一步的分析表明，在低碰撞能83.64 kJ·mol⁻¹时，“迁移机制”占据主导地位；而在高碰撞能209.1 kJ·mol⁻¹时，“夺取机制”占据主导地位。F + HCl反应的很多动力学结果，如反应物HCl的转动激发促进反应的发生，产物HF的振、转态分布，微分散射截面，乃至产物HF分子空间转动取向极化行为都可以利用两种机制的竞争关系来进行解释。

基于准经典轨线方法和Matlab GUI工具，我们开发设计了将a + bc三原子微观碰撞过程可视化的动力学模拟程序。该程序能够以动画的形式动态地展现微观原子、分子之间的碰撞过程和碰撞结果，可以将微观的化学反应变成直观的物理图像，有助于把抽象问题具体化和形象化。模拟程序运行后的图形界面具有较好的交互性，学生可以通过控件和菜单来观察和控制各种反应类型的微观碰撞过程，从而对微观化学反应过程、反应机理进行深入的分析和研究。我们期望通过微观化学反应可视化模拟程序的开发设计，能够丰富化学反应动力学课程的教学方法和手段，提升教学的趣味性，进一步提高化学反应动力学课程的教学效果。