1.   引　言

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

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

综上可知, 研究者们对${\rm CH}_{2}^{+}$体系进行了大量的理论研究. 大部分势能面都是采用多项式拟合的方式构建的, 而神经网络拟合方法与传统的拟合方法相比具有精度高、运算快的优势. 本文的目标就是构建更高精度的势能面并在其基础上进行动力学计算.

2.   势能面

2.1   从头算计算

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

2.2   势能面拟合

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

式中, ${P}_{{\mathrm{C}\mathrm{H}}_{a}^{+}}=1/{R}_{{\mathrm{C}\mathrm{H}}_{a}^{+}}$, ${P}_{{\mathrm{C}\mathrm{H}}_{b}^{+}}=1/{R}_{{\mathrm{C}\mathrm{H}}_{b}^{+}}$, ${P}_{\mathrm{H}\mathrm{H}}=1/{R}_{\mathrm{H}\mathrm{H}}$, 其中${R}_{\mathrm{H}\mathrm{H}}$为氢分子之间的键长, ${R}_{{\mathrm{C}\mathrm{H}}_{a}^{+}}, {R}_{{\mathrm{C}\mathrm{H}}_{b}^{+}}$为CH⁺分子的键长. 在计算中, 采用了2个隐藏层, 每个隐藏层包含15个神经元. 此外, Levenberg-Marquardt方法^[39]用来训练神经网络. 神经网络的表达式为

其中, ${N}_{\mathrm{t}\mathrm{o}\mathrm{t}}$是从头算能量点的个数$;{w}_{i}^{1}$, ${w}_{j}^{2}$和${w}_{k}^{3}$是连接权重$;{b}_{1}^{3}$, ${b}_{k}^{2}$和${b}_{j}^{1}$是相位$;{f}^{1}$和${f}^{2}$是训练函数, 其数学表达式为

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

其中, N是总的能量点个数, ${y}_{i}^{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{p}\mathrm{u}\mathrm{t}}$是拟合值, ${y}_{i}^{ab}$是从头算能量. 对于本文的${\mathrm{C}\mathrm{H}}_{2}^{+}$体系, 拟合RMSE仅为1.27 meV.

2.3   势能面特征

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

表 1  CH⁺和H₂分子的光谱常数

Table 1.  Spectroscopic constants of the CH⁺ and H₂ molecules.

	re /Bohr	ωe /cm–1	ωexe /cm–1	βe /cm–1	αe /cm–1	De /eV
CH+(${X}^{1}{{\Sigma } }^{+}$)
本文结果	2.136	2860.31	59.32	14.217	0.501	4.257
理论[11]	2.136	2853.03	58.52	14.201	0.489	4.252
理论[25]	2.136	2851.0	58.1	14.199	0.489	4.244
理论[26]	2.144	2849.03	66.45	14.094	0.490
理论[17]	2.135	2861.95	59.63	14.311	0.447	4.257
实验[40]	2.137	2857.56	59.32	14.178	0.495	4.26
H2(${X}^{1}{{\Sigma } }_{\mathrm{g} }^{+}$)
本文结果	1.401	4407.63	139.43	60.85	3.012	4.751
理论[17]	1.401	4404.61	126.64	60.861	2.233	4.749
理论[27]	1.401	4403.60	126.60	60.864	2.232	4.748
理论[28]	1.403	4395.22	126.12	60.735	2.221	4.748
理论[29]	1.401	4389.66	121.56	60.826	3.162	4.711
实验[41]	1.401	4401.21	121.33	60.853	3.062	4.746

下载: 导出CSV 

| 显示表格

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

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

图 1  (a)键角在136.68°时, 化学键伸缩的等势线图(等势线的起点为–9.1 eV, 间隔为0.4 eV); (b) C⁺离子以C_2v对称性接近H₂分子中心的等势线图(等势线的起点为–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 H₂ molecule in the C_2v symmetry (Contours starting at –9.1 eV and equally spaced by 0.37 eV).

下载: 全尺寸图片 幻灯片

图 2  (a) 当${R}_{{\mathrm{C}\mathrm{H}}^{+}}=2.136$ Bohr时, H原子绕CH⁺离子运动的等势线; (b) 当${R}_{\mathrm{H}\mathrm{H}}=1.401$ Bohr时, C⁺离子围绕H₂分子运动的等势线.

Fig. 2.  (a) Contour plot for the H atom moves around CH⁺ ion at the bond distance ${R}_{{\mathrm{C}\mathrm{H}}^{+}}=2.136$ Bohr; (b) contour plot for C⁺ ion moves around H₂ molecule at its equilibrium geometry ${R}_{\mathrm{H}\mathrm{H}}=1.401$ Bohr.

下载: 全尺寸图片 幻灯片

图3给出了CH${}_2^+$体系基态势能面在不同角度下的最小能量路径, 为了与之前的理论结果进行比较, Li等^[11]的结果也列入了图中. 如图3所示, H + CH⁺通道的能量设置为0. 很明显, C⁺ + H₂ → CH⁺ + H反应是一个吸热反应, 在不考虑零点能的情况下, 其吸热能约为0.5 eV, 这与Li等^[11]报道的结果十分接近. 此外, 在小角度时, 能量路径上存在两个势阱, 分别位于势垒的两侧. 然而随着角度的增加, 势垒消失, 势阱变为一个, 并且势阱的深度持续变浅. 角度为180°时, 在能量路径上出现了一个微小势垒, 该势垒的高度几乎与H + CH⁺通道齐平. 此外, Li等^[11]的最小能量路径与本文的最小能量路径基本一致, 一些细微的差别可能是因为使用不同的拟合方法导致的.

图 3  CH${}_2^+$势能面不同角度的最小能量路径以及来自文献[11]的理论结果

Fig. 3.  The minimum energy paths of the CH${}_2^+$ potential energy surface at different angles along with the results obtained from Ref.[11]

下载: 全尺寸图片 幻灯片

3.   动力学

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

其中, R_tra是反应轨线数目, N_tra是总轨线数目, Δ_θ是角度步长, N_θ是角方向格点数目.

图4给出了C⁺ + H₂ → 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轴方向上的投影进行了截断(K_block = 20). 对于固定的J值, 完全的科里奥利耦合计算在z轴方向上的投影应为J + 1, 这表明Guo等^[12]的结果在J ≥ 20时就可能得不到收敛的动力学结果. 因此, Guo等^[12]的结果很可能在低能处是收敛的, 而对于较高的碰撞能, 他们的结果很可能不收敛. 本文采用的QCT方法不能考虑低能处的量子效应, 因此, 在低能处, 本文的结果与Guo等^[12]的结果有一些差距. 在高能处, 量子效应变得不那么明显, 量子结果应与QCT结果取得类似的结果. 然而, Guo等^[12]的结果与本文的QCT结果依然存在着较大的差距, 这可能是因为Guo等^[12]截断了J在z轴方向上的投影, 从而导致量子结果不收敛.

图 4  C⁺ + H₂反应碰撞在0.4—1.0 eV范围内的积分截面以及文献[12, 18]的结果

Fig. 4.  Integral cross section of C⁺ + H₂ reaction in the collision energy range from 0.4 to 1.0 eV along with the values obtained from Ref. [12, 18].

下载: 全尺寸图片 幻灯片

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

图 5  C⁺ + H₂反应若干碰撞能下的微分截面

Fig. 5.  Differential cross sections of the C⁺ + H₂ reaction at several collision energies.

下载: 全尺寸图片 幻灯片

4.   结　论

本文采用PIP-NN构建了${\rm CH}_2^+$体系基态势能面, 使用MRCI-F12方法结合ACVQZ基组, 对${\rm CH}_2^+$体系进行了从头算计算, 详细描述了势能面的特征, 并与之前的理论报道进行了详细的比较. 结果表明, 新构建的势能面足够精确, 可以用于${\rm CH}_2^+$体系一切类型的动力学计算. 为了进一步验证新构建势能面的准确性, 基于新构建的势能面, 采用QCT在0.4—1.0 eV的碰撞能范围内对C⁺ + H₂ → H + CH⁺反应进行了计算. 报道了总积分截面和微分截面, 并与之前的理论结果进行了比较. 对于积分截面, 本文的结果大体上与之前的理论结果保持一致. 微分截面的结果表明插入反应机理在反应过程中占据主导地位.