C++H2 反应的动力学研究:基于一个新构建的势能面*

2022-05-26李文涛袁美玲王杰敏

物理学报 2022年9期

李文涛袁美玲王杰敏

1)(洛阳师范学院物理与电子信息学院,洛阳 471934)

2)(潍坊科技学院通识学院,寿光 262700)

3)(枣庄学院人工智能学院,枣庄 277160)

采用一种明显依赖于电子-电子距离的多参考组态相互作用方法,结合ACVQZ 基组进行了体系的从头算计算.采用置换不变多项式神经网络方法对从头算得到的18222 个单点能进行了拟合.此外,将势能面的特征与可得到的实验结果和理论结果进行了比较.结果表明新构建的势能面比以往的势能面更加精确.基于新构建的势能面,采用准经典轨线方法进行了C++H2 反应的动力学计算,报道了积分截面和微分截面等动力学信息,并与之前的理论结果进行了比较.动力学结果表明在反应过程中插入反应机理占据主导地位.

1 引言

自1941 年Douglas 和Herzberg[1]在星际介质中发现CH+离子的高丰度含量后,CH+离子受到了大量的实验[2−6]和理论[7−29]研究.CH+离子是高活性的,在宇宙环境中很容易与电子、氢原子、氢分子等发生反应,从而被迅速破坏掉.因此,CH+离子在宇宙中的高丰度就变得难以解释,这也是现代天体物理学中的一个难题.为解释其在宇宙中高丰度的现象,在接下来的几十年中,人们对CH+离子的形成和破坏机制进行了广泛的研究.在实验方面,1997 年,Hierl 等[2]采用高温流动余热装置测量了C+(2P)和H2分子及其同位素D2分子反应在400—1300 K 内的热速率常数.实验结果表明,反应物分子的转动能和平动能对反应的影响类似,而反应物分子的振动激发对反应的影响远大于平动能和转动能.与反应物的振动基态v=0比较,v=1 的振动激发,使H2分子速率常数增大约1000 倍,使D2分子速率常数增大约6000 倍.2011 年,Plasil 等[3]采用温度可变的22 极离子阱和冷喷氢原子束,研究了CH+离子与H 原子的碰撞反应.他们报道了该反应在10—1000 K 温度范围内的速率常数,发现当温度低于100 K 时,速率常数会出现一个迅速的下降.这说明非转动的CH+离子在低温情况下很难与H 原子发生碰撞反应.

在理论方面,大量的理论研究工作者采用构建高精度势能面并进行动力学计算的方式来获得指定电子态的速率常数.2005 年,Stoecklin 和Halvick[7]采用多项式拟合的方法构建了体系的基态势能面,并基于该势能面采用RIOSA-NIP 方法计算了H+CH+反应的反应概率、积分截面和速率常数.2007 年,Halvick 等[8]基于该势能面采用准经典轨线方法(quasiclassical trajectory method,QCT)和相空间理论重新研究了H+CH+反应,并报道了该反应的速率常数.2013 年,Zanchet 等[9]基于该势能面采用含时量子波包方法研究了C++H2反应,并报道了该反应量子态分辨的速率常数.2011 年,Warmbier 和Schneider[10]报道了体系的一个新势能面.他们采用AVTZ和AVQZ 基组和多参考组态相互作用方法(multireference configuration interaction method,MRCI)进行了从头算计算.此外,基于该势能面还采用了QCT 方法和ABC 非含时的量子方法对H+CH+反应进行了动力学计算,并报道了积分截面和速率常数等动力学信息.2015 年,Li 等[11]采用MRCI方法和AV6Z 基组进行了高精度的从头算计算,并采用多项式拟合的方法构建了体系高精度的势能面.该势能面得到广泛的应用.Guo 等[12],Sundaram 等[13−15]和 Wu 等[16]采用该势能面对H+CH+反应及其逆反应进行了动力学计算并报道了积分截面、速率常数等动力学信息.最新的势能面是由Guo 等[17]构建的,在从头算计算中,他们把基组外推到了完全基组极限,同样采用多项式拟合的方法构建了新的势能面.

2 势能面

2.1 从头算计算

本文采用Molpro 程序包进行从头算计算,并且对所有的构型都采用了Cs对称性.在从头算计算中,采用了一种明显依赖于电子-电子距离的内收缩多参考组态相互作用方法(MRCI-F12)[30].此外,对C+离子和H 原子采用了ACVQZ 相关一致基组[31,32].对于所有构型,首先执行包含4 个等权重电子态(3 个2A',1 个2A'')的态平均完全活化空间自洽场(state average-complete active space self-consistent field,SA-CASSCF)[33,34]计算,然后利用SA-CASSCF 计算提供的参考波函数进行了MRCI-F12 计算.为了补偿高阶截断误差,计算中考虑了Davidson 修正.在计算中采用雅克比坐标来排布空间格点(RQ,RHH,α).其中,RQ是C+离子到H2分子的中心距离,在0.6—25 Bohr 的范围内放置55 个非均匀格点;RHH是H2分子的键长,在0.6—25 Bohr 的范围内放置56 个非均匀格点;α是RQ和RHH之间的夹角,其变化范围是1°到89°,间隔为15°.在拟合过程中,删掉了能量高于基态最低能量10 eV 的从头算能量点.最终,在拟合过程中采用了18222 个从头算能量点.

2.2 势能面拟合

置换不变多项式神经网络(permutive invariant polynomial-neuronal network,PIP-NN)[35,36]方法已经广泛地应用于势能面构建,例如H +HBr[37]和K+H2反应[38]等.在此,简要地介绍一下PIP-NN 方法.众所周知,神经网络方法由3 个部分组成,即输入端、隐藏层、输出端.对于输入端,为了避免在边界处一阶导数不连续,采用PIP方法对雅克比坐标进行如下变换:

其中,Ntot是从头算能量点的个数;是连接权重;是相位;f1和f2是训练函数,其数学表达式为

过拟合问题一直是神经网络方法中的一个棘手问题,为了避免过拟合的发生,将输入数据随机分为3 个部分,即训练部分(90%)、测试部分(5%)和验证部分(5%).采用均方根误差(root mean square error,RMSE)来测试势能面的精确度,它的表达式为

其中,N是总的能量点个数,是拟合值,是从头算能量.对于本文的体系,拟合RMSE仅为1.27 meV.

2.3 势能面特征

测试势能面的精确程度的一个有效途径就是将计算得到的光谱常数与实验进行比对.当第3 个原子远离另外两个原子中心距离为30 Bohr 时,改变另外两个原子之间的间距,就可以在势能面上得到CH+和H2分子的势能曲线.将该曲线用最小二乘法拟合得到光谱常数,与实验结果和其他理论结果一并列入到了表1 中.

如表1 所列,本文拟合得到的CH+和H2分子的光谱数据与实验值十分吻合.对于CH+离子,平衡核间距re、谐振频率ωe、离解能De与实验值的差距仅为0.001 Bohr,2.75 cm–1,3 meV.对于H2分子,平衡核间距re与实验结果一致,谐振频率ωe和离解能De与实验值的差距为6.42 cm–1和5 meV.实验与理论的微小差距说明新构建的势能面很好地描述了当第3 个原子远离另外两个原子时的势能曲线.

表1 CH+和H2 分子的光谱常数Table 1. Spectroscopic constants of the CH+ and H2 molecules.

图1 和图2 展示了这个新构建的 CH+2势能面的主要地形特征.由图1 可知,在研究涉及的空间范围内等势线正确并光滑,其等势线对应反应过程中一些主要的驻点.图1(a)描述了键角∠HC+H=136.68°时,RCH+键伸缩的等势线图.很明显,图1(a)中有一个很深的势阱,对应整个势能面的最低能量点,对应的结构为RCH+=2.08 Bohr,RHH=3.866 Bohr,相比C++H2通道的能量,势阱的深度约为4.30 eV.图1(b)描述了C+离子垂直接近H2分子中心点过程的等势线.可看到两个驻点分别对应于局域最小值和整体最小值.其中局域最小值的构型为RCH+=2.602 Bohr,RHH=1.695 Bohr.图2(a)给出了当CH+离子处于平衡核间距的构型下,H原子围绕CH+离子运动的等势线图.在C+离子一侧有一个明显的深势阱.当H 原子在∠HQC+=0°—60°范围接近CH+离子时,很容易被该势阱捕获,从而形成长寿命的中间产物.图2(b)给出了当H2分子处于平衡键长时,C+离子围绕H2分子的等势线图.在90°附近有一个明显的深势阱.可以预见:当C+离子接近H2分子时,很可能被该势阱吸引,从而进入T 型结构;对于C++H2→CH++H 反应,插入反应机理对该反应起着重要作用.

图1 (a)键角在136.68°时,化学键伸缩的等势线图(等势线的起点为–9.1 eV,间隔为0.4 eV);(b) C+离子以C2v 对称性接近H2 分子中心的等势线图(等势线的起点为–9.1 eV,间隔为0.37 eV)Fig.1.(a) Contour plot for chemical bond stretching,in which the angle is fixed at 136.68° (Contours starting at–9.1 eV and equally spaced by 0.4 eV);(b) contour plot for the C+ ion approach to the midpoint of H2 molecule in the C2v symmetry (Contours starting at–9.1 eV and equally spaced by 0.37 eV).

图2 (a) 当 RCH+=2.136 Bohr 时,H 原子绕CH+离子运动的等势线;(b) 当 RHH=1.401 Bohr 时,C+离子围绕H2 分子运动的等势线.Fig.2.(a) Contour plot for the H atom moves around CH+ ion at the bond distance RCH+=2.136 Bohr;(b) contour plot for C+ion moves around H2 molecule at its equilibrium geometry RHH=1.401 Bohr.

图3 势能面不同角度的最小能量路径以及来自文献[11]的理论结果Fig.3.The minimum energy paths of the potential energy surface at different angles along with the results obtained from Ref.[11].

3 动力学

为了进一步验证势能面的准确性,基于新构建的势能面,采用QCT[42]对C++H2反应以0.05 eV的步长,在0.4—1.0 eV 的碰撞能范围内进行了动力学计算.当入射C+原子与H2分子的质心的距离小于8 Å时,轨线开始计算,并且对于每一个能量随机采样1×105条轨线.在轨线的计算中,为了确保数值的稳定性,采用的积分步长为0.02 fs.在正式开始计算之前,对于每一个能量点,需要对最大碰撞参数bmax进行调节.在本文的计算中,采取测试的轨线数目为5×103条,在测试过程中逐渐增大bmax直至没有反应轨线出现.在积分截面和微分截面的计算中,碰撞参数b在0 和最大碰撞参数bmax之间随机取值.当产物分子之间的距离大于8 Å时,轨线终止计算.积分截面和微分截面可以通过下列公式得到:

其中,Rtra是反应轨线数目,Ntra是总轨线数目,Δθ是角度步长,Nθ是角方向格点数目.

图4 给出了C++H2→H+CH+反应的积分截面.为了便于与之前的理论结果进行比较,Guo等[12]和Herráez-Aguilar 等[18]的结果也列入到了图中.如图4 所示,积分截面的大小随着碰撞能的增加而增大,Guo 等[12]和Herráez-Aguilar 等[18]的结果与本文结果具有相同的趋势.尽管如此,本文的结果与Guo 等[12]和Herráez-Aguilar 等[18]的结果之间还存在着明显的不同.Herráez-Aguilar等[18]的结果是基于Warmbier 和Schneider[10]的势能面,采用QCT 得到的.与本文的结果比较,在研究的碰撞能范围内,本文的结果明显大于Herráez-Aguilar 等[18]的结果,这是因为在计算中采用了不同的势能面.Guo 等[12]的结果是基于Li 等[11]的势能面,采用含时量子波包方法计算得到的.然而在他们的计算中,对角动量J在z轴方向上的投影进行了截断(Kblock=20).对于固定的J值,完全的科里奥利耦合计算在z轴方向上的投影应为J+1,这表明Guo 等[12]的结果在J≥20 时就可能得不到收敛的动力学结果.因此,Guo 等[12]的结果很可能在低能处是收敛的,而对于较高的碰撞能,他们的结果很可能不收敛.本文采用的QCT方法不能考虑低能处的量子效应,因此,在低能处,本文的结果与Guo 等[12]的结果有一些差距.在高能处,量子效应变得不那么明显,量子结果应与QCT 结果取得类似的结果.然而,Guo 等[12]的结果与本文的QCT 结果依然存在着较大的差距,这可能是因为Guo 等[12]截断了J在z轴方向上的投影,从而导致量子结果不收敛.

图4 C++H2 反应碰撞在0.4—1.0 eV 范围内的积分截面以及文献[12,18]的结果Fig.4.Integral cross section of C++H2 reaction in the collision energy range from 0.4 to 1.0 eV along with the values obtained from Ref.[12,18].

图5 给出了C++H2→ H+CH+反应在不同碰撞能下的微分截面.从图5 不难看出,对于所有能量点,散射信号都集中在0°和180°附近,并且前后向对称.这说明,在反应过程中插入反应机理占据主导地位.反应路径上的深势阱(图3)支持长寿命的束缚态和准束缚态,从而使得产物可以到达任何可接触到的通道.