D+SiD+ →D2+Si+反应量子波包动力学研究*

赵文丽 孙丰伟 张红 王永刚 高峰† 孟庆田

1)(山东农业大学信息科学与工程学院,泰安 271018)

2)(山东师范大学物理与电子科学学院,济南 250358)

本文基于最新构建的 势能面,运用切比雪夫波包方法,对初始态为ν=0,j=0的D +SiD+反应体系在 1 .0×10−3 —1.0 eV的碰撞能量范围内进行动力学研究.通过对角动量量子数 J ≤110的所有分波进行计算,运用耦合态近似和考虑科里奥利耦合效应的精确量子力学两种方法,得到该反应的反应概率、积分散射截面和速率常数.文章详细分析了碰撞能量和科里奥利耦合效应对该反应的动力学性质的影响,并与H+SiH+反应做了对比.结果表明,忽略科里奥利耦合效应会使D +SiD+反应积分散射截面和速率常数增大.与H+SiH+反应动力学相比较,D +SiD+反应的积分反应截面的值较小且随着能量的衰减非常缓慢,在低能区没有H+SiH+随能量增大而剧烈衰减的情况.计算发现,相同温度下,D +SiD+反应速率常数小于H +SiH+反应速率常数,随着温度的升高,二者的差距减小.这表明,同位素替代对反应的动力学性质有明显的影响.

1 引言

2019年,笔者团队[15]首次报道了关于(X2A1)的一个精确的全维势能面(GMS-PES).计算中应用多体展开(MBE)方法拟合所有的从头计算能量点,所有能量点都采用了包含Davidson 修正的多参考组态相互作用[MRCI(Q)][16]方法.从头计算能量点的跨度范围为24.30 eV,总方均根偏差(rmsd)为0.04 eV.前期应用该势能面的部分研究成果[17]如Si++H2的动力学计算结果和其他理论值相比,用该势能面计算的结果与实验值更接近,证明GMS-PES是精确且适用的.2021年,笔者小组[18]应用GMS-PES 研究了H(2S)+SiH+反应的量子动力学,发现在 1.0×10−3—1.0 eV 碰撞能量范围内,用包含科里奥利耦合(CC)效应的精确量子方法和耦合态近似(CS)方法计算所得的积分散射截面(ICS)差别不大,但是在300 K 左右的温度下,CC 方法计算所得的速率常数比CS 计算值大50%左右,因此H(2S)+SiH+反应中,科里奥利耦合效应非常重要,不可忽略.

有关同位素取代动力学的研究,已有一些工作进行了报道[19−21].文献[19,20]报道了分子间同位素会对反应立体动力学性质产生一定的影响,如N(4S)+H2,H+NH 反应物中的H 被D 代替,则微分反应截面会减小,因为质量因子的减小会抑制反应的发生.文献[21]研究了同位素效应对C++H2→CH++H 反应立体动力学的影响,研究发现将反应中的H2用HD和HT 代替以后,同位素效应会越来越重要.本文的主要工作是基于GMSPES,应用切比雪夫波包方法[22,23]计算 D +SiD+反应的量子动力学,并且与前期的H +SiH+反应动力学[18]进行对比,详细讨论该反应的同位素效应.计算中采用了CS 近似和考虑CC 效应的精确量子方法.本文首先介绍所采用的理论方法以及GMS-PES的结构特征,然后对计算结果进行讨论.除特别说明之外,文中的所有公式和计算都使用原子单位.

2 理论方法

2.1 切比雪夫波包计算

其中R表示H 原子与SiH+离子质心之间的距离,r表示Si 与H+之间的距离,µr=mSimH/(mSi+mH)为SiH+双体约化质量,µR=mH(mSi+mH)/(2mH+mSi)为H+SiH+体系三体约化质量.V(R,r,γ)是体系的相互作用势能,其可以从势能面中获得.分别是是双原子SiH+轨道角动量算符和转动角动量算符.轨道角动量算符的平方可以表示为

切比雪夫波包方法不存在离散或者插值误差,波包传播的递归关系为[24,25]

其中初始波包|ψ0〉表示为

式中N是归一化常量,k0为初动量,R0与δ分别表示初始波包的中心位置和宽度,|φi〉为振转本征函数,是归一化哈密顿量,H±=(Hmax±Hmin)/2,Hmax(Hmin)通过光谱的上(下)限计算获得[26].在格点边缘,需要用高斯型函数D来描述出射波的边界条件:

其中,xd表示阻尼起始点,阻尼范围由dx决定.

基于流计算方法[27]得到总反应概率:

式中rf决定了产物通道的分界面,θ ≡arccos[(E−H+)/H−]是切比雪夫角,ai(E)是初始波包的振幅[28]:

其中hli(2)(kiR)为第二类球汉克尔函数.

对于一个给定的碰撞能量Ec,将所有的总角动量量子数对应的反应概率求和即可得出给定初态的ICS:

式中,f是电子的简并因子,kB和T分别为玻尔兹曼常数和热力学温度.

2.2 势能面与MEP

势能面的特征对分子反应动力学行为有深刻的影响,下面讨论一种典型的情况.图1为平衡构型∠[H-Si-H]=118.38°时的三维势能面,相邻等势线的能量间隔为0.2 eV,图中的黑色曲线是等势线,可以看出所有的等势线都非常光滑,势能面最显著的拓扑特征是存在一个全局最小(GM)点,RHH'=4.820a0,RSiH=RSiH＇=2.806a0.显然,势能面存在完美的对称性,这是由于体系中两个氢原子的交换对称性所决定的.图中的红色曲线为该角度下的最小能量路径(MEP),这条MEP是以RSiH和RSiH′定义的函数,其中RSiH和RSiH′分别代表Si 原子和两个H 原子的核间距离.显然在这条路径上不存在任何势垒,只有一个5.82 eV的深势阱,这个深势阱与反应过程中的中间络合物的形成密切相关.

图1 势能面(∠[H-Si-H]=118.38º),图中等势线间隔为0.2 eV,红色实线为最小能量路径Fig.1.The PES of for bending angle∠[H-Si-H]=118.38º starting from–5.82 eV,the contour drawn with an increment of 0.2 eV,the solid red line is minimum energy path..

为了提供更多的关于势能面的信息,可以研究一般意义下的MEP.图2 所示为H+原子从不同角度插入时,体系的势能随着反应坐标RSiH+−RHH的变化关系.从图2 可以看出,D+SiD+→D2+Si+是放热反应,放出的热量约为1.30 eV.在30°时,存在一个2.39 eV的深势阱.随着角度的增大,势阱的深度变小.值得注意的是在30°和60° 角处存在一个双势阱结构,30° 情况下靠近反应物渐近通道的势阱深为2.40 eV,对于60°的情况靠近反应物渐近通道势阱的深度为0.72 eV,在此位置附近有一个过渡态(TS).在直线构型的情况下存在着一个鞍点,靠近产物渐近通道出口处有一个1.45 eV的深势阱紧随其后,这个鞍点对应着体系第一激发态的一个极小值.T 构型下MEP 有一个1.66 eV的深势阱,其他一些极小值一般都靠近产物渐近通道,与T 构型的情况类似.

图2 最小能量路径(∠[D-Si-D]分别为30°,60°,90°,120°,180°)Fig.2.The MEP of different approaching angles(∠[D-Si-D]=30°,60°,90°,120°,180°)for the title reaction.

2.3 计算细节

表1 波包计算中的数值参量(除特殊说明,均采用原子单位 a.u.)Table 1.Model parameters of wave packet calculation(The atomic unit is used in the calculation unless otherwise stated).

3 讨论

3.1 反应概率

图3 画出了CS和CC 两种计算方法下D+SiD+反应在角动量量子数J=5,20,40,60,80,100 时的反应概率.显然,D+SiD+反应是无势垒的放热反应,反应概率随着碰撞能量的增大而减小,在J≤40的情况下,反应概率在高能区域有略微随能量增大而增大的趋势.随着总角动量J的增大,反应的阈能逐渐增大,这是由于离心势垒逐渐增大所引起的,可以用J挪移方法[30]或者俘获模型[31,32]解释.CS 概率呈现密集的振荡,并且振荡幅度较大,这时因为受到势能面中深势阱的影响,中间络合物会有较长的寿命.进一步观察发现,CC概率的振荡频率和幅度都比CS的小很多,这是由于反应中科里奥利耦合效应受到长程力支配的原因.当J较小(J=5)时,CC 概率与CS 概率差别很小.当J=20,40,60,80时,CC 概率在整个能量范围内都小于CS 概率.对于J较大的单一分波,如J=100时,在E<0.8 eV的区域CC 概率小于CS 概率,E>0.8 eV时,CC 概率大于CS概率.整体上来讲,在CS 概率振荡较大的区域,CC 与CS 概率差值较大,在CS 概率振荡较小的区域,CC 与CS概率值基本吻合.在CS 近似中,SiH+离子被限制在k=0(总角动量J在定体坐标z轴的分量)量子态,这将会缩短碰撞时间,增大共振进而使得反应概率值变大或者变小[33,34,35].

图3 D +SiD+ →D2+Si+反应不同分波(J=5,20,40,60,80,100)的反应概率随碰撞能量的变化,黑色实线对应CS 概率,红色虚线对应CC 概率Fig.3.The reaction probabilities of CC and CS calculations for D+SiD+ →D2+Si+ reaction at J=5,20,40,60,80,and 100,The black soid line is for CS probability and the red dashed line is for CC probability.

图4为H+SiH+与D+SiD+反应的CC 概率比较.在相同的角动量量子数下,D+SiD+的反应阈能小于H+SiH+反应阈能.在各种角动量下,H+SiH+反应概率都会随着碰撞能量的增大而减小;D +SiD+的情况却不同,在J ≤40的情况下,低能区其反应概率小于H+SiH+的反应概率,但是在高能区随着能量增大,反应概率有上扬逼近甚至超越H+SiH+反应概率的趋势.这种现象可以从两个角度分析:从势能面的角度来看,D+SiD+体系的约化质量大于H+SiH+体系的约化质量,所以D+SiD+反应物更容易被困在深势阱中处于中间络合物的状态,从而降低反应概率;从反应概率(6)式来考察,双体和三体约化质量都处于分母的位置,且D+SiD+体系的双体和三体约化质量都是H+SiH+体系对应量的两倍,能够有效降低反应概率.当J>60时,D+SiD+的反应概率却大于H+SiH+的反应概率.这是因为,当J较大时,D+SiD+体系的能量迅速增大,能够有效促进反应的发生.

图4 不同的角动量量子数下H+SiH+与D+SiD+的CC 反应概率,黑色实线对应H+SiH+反应,红色虚线对应D+SiD+反应Fig.4.The reaction probabilities of CC calculations for H +SiH+ and D +SiD+ reactions at J=10,20,40,60,80,and 90.The black solid line is for H+SiH+ reaction and the red dashed line is for D +SiD+ reaction.

3.2 积分散射截面

图5 给出了H+SiH+与D +SiD+反应在1.0×10−3—1.0 eV的碰撞能量范围下ICS 随碰撞能量的变化,计算中用到CS和CC 两种方法.观察D+SiD+反应的两种计算结果发现,忽略CC 效应使得ICS 增大,这种情况在低能区域尤其明显.两个反应的ICS 都没有阈能(零阈能),H+SiH+的ICS 在低能区域内会随着碰撞能量的增大急剧的减小,然后变化趋势趋于平缓,而 D +SiD+的ICS 反应在整个能量范围内减小的趋势都非常缓慢,这些动力学行为都符合无势垒的放热反应.此外,随着能量的增大,CC和CS 两种方法计算所得的ICS的差值越来越小,一个可能的原因就是在高能区偶发碰撞会增大,进而提高CC 碰撞概率.与H+SiH+反应比较,D +SiD+反应的ICS共振结构小很多,特别是CC 方法计算的ICS,曲线比较光滑,基本没有共振结构.这是因为ICS是反应概率加权求和的结果,在CC 计算中会用到更多的分波,导致绝大部分共振被消除.

图5 H+SiH+与D+SiD+反应 ICS 随着碰撞能量的变化(1.0×10–3—1.0 eV)Fig.5.The ICSs of CS and CC calculations for H+SiH+ and D +SiD+ reactions over collision energy of 31.0×10–3—1.0 eV.

3.3 速率常数

基于ICS,计算了温度在100—1000 K 之间H+SiH+和D +SiD+反应的速率常数,如图6 所示.显然,忽略CC 效应使 H +SiH+反应的速率常数减小,使 D +SiD+反应的速率常数增大.这是因为,速率常数是ICS的加权积分结果,根据图5 所示的两个反应CC和CS的ICS 比较,这样的结果是可以预料的.以CC的计算结果来比较两个反应的速率常数发现,H +SiH+的速率常数基本不随温度变化,而 D +SiD+反应的速率常数随着温度的升高单调增大,特别是在常温区域附近D+SiD+反应速率常数对于温度的变化非常敏感.D+SiD+反应的速率常数整体看来都小于 H +SiH+反应的速率常数,说明质量因子对速率常数的影响是比较大的.例如,在T=200 K时,H +SiH+反应的速率常数k=1.76×10–9cm3/(molecule·s),D+SiD+的速率常数k=4.30×10–10cm3/(molecule·s).前者反应的速率常数是后者的4倍,随着温度的升高,二者的差距逐渐减小,当温度为1000 K时,前者反应的速率常数是后者的2.60 倍.

图6 H+SiH+和D +SiD+ 反应速率常数随温度的变化Fig.6.Variation of reaction rate constant of H+SiH+ and D +SiD+ with temperature.

3.4 同位素效应

动力学同位素效应(KIE)是另外一个定量衡量分子反应动力学性质的物理量,图7画出了基于CC 计算方法的同位素效应κH+SiH+(T)/κD+SiD+(T).曲线表明,KIE 值随着温度的升高而单调减小,温度越高,KIE 曲线越平缓.根据(9)式,一方面,速率常数大小与ICS的玻尔兹曼积分成正比,H+SiH+反应的ICS 在低能区域远大于 D +SiD+反应的ICS,所以,在低温区域,ICS 对KIE 起着主要作用.另一方面,速率常数大小与反应物约化质量的开方成反比,的比值大于1,为1.39,在高温区域,速率常数取决于反应物的约化质量,同时考虑到这个区域H+SiH+反应的ICS也略大于 D +SiD+反应的ICS,这两种因素叠加导致整个考察的温度范围内kH+SiH+(T)/kD+SiD+(T)值都大于2.这种结果也可以用隧穿效应来解释:与D 原子相比,H 原子质量更轻,运动速度更快,更容易发生隧穿效应与其他分子发生反应.

图7 动力学同位素效应与温度的关系Fig.7.Temperature dependence of kinetic isotope effect.

4 结论

本文用切比雪夫波包方法研究了 D +SiD+反应动力学.通过CC和CS 计算结果相比较发现忽略科里奥利耦合效应会使 D +SiD+反应的积分散射截面和速率常数增大,这归因于中间络合物机制.与H+SiH+反应相比较,在低能区域D+SiD+反应的ICS 远小于H+SiH+反应的ICS,在高能区域内略小于H+SiH+反应的ICS.在100—1000 K的温度内,D +SiD+反应的反应速率小于同温度下H+SiH+反应的反应速率.通过对同位素效应的分析发现,常温下,温度越低,κH+SiH+(T)/κD+SiD+(T)的值越大.迄今为止,还没有发现关于D+SiD+反应动力学的实验或者理论方面的报道,这些结论都期待着实验或者其他理论的验证.