1. 引言

由于NH的衰退在氨氮热解和氮化合物气相燃烧反应中起关键作用，H + NH的基态吸氢反应在亚胺根(NH)的衰变中起着关键作用。近年来人们对H + NH的基态氢提取反应越来越感兴趣[1] [2]。从理论上报道的文献来看，许多经典动力学和量子动力学的研究都集中在NH2体系上[3]。Adam等人基于多参考构型相互作用(MRCI)计算，构建了NH2系统(4A'') [4]和NHD系统(2A'') [5]基态的全局势能面(PESs)。然后，他们使用其全局势能面进行准经典轨线(QCT)计算[4]。几个版本的过渡态理论(TST)被用来从理论上确定较宽温度范围内的速率常数[6]。Rascual等人[7]利用QCT和变分TST方法在各自的势能面上进行了计算。同时，Poveda和Varandas [8]报道了基于MRCI/augc-cc-pVQZ水平的双多体展开(DMBE)理论的真实全局PES。Han等人[9]在DMBE PES上使用QCT和量子力学(QM)方法进行了计算。重要的是，不同研究小组[10] [11]用三个PESs上详细研究了Renner-Teller (RT)耦合对NH2体系的影响。最近有研究小组对基态用多体展开法求解二体项，用神经网络方法求解三体项，构建了一个新的高质量的PES [12]。在0~1.0 eV的碰撞能量范围内，利用新PES对H + HN反应进行了量子反应散射动力学计算。

Figure 1. The definition of the attack angle θ which is the angle between lines HO and DO

图1. HO和DO连线间攻击角的定义

尽管H + HN反应及其同位素效应引起了人们的广泛关注，但对于该体系的研究仍存在一些关键问题。首先，人们对其同位素变效应的动力学知之甚少。此外，以往的一些研究主要集中在产物极化和反应标量性质上。为了更好地理解H(D) + H(D)N反应的动力学，还需要研究化学立体动力学，它可以为碰撞过程提供有价值的信息。众所周知，反应几率取决于反应物的相互取向。然而，到目前为止，还没有报道解决H + DN反应的这个问题。相互取向用攻击原子的接近方向与双原子间轴的夹角(称为攻击角)来描述，我们在图1中定义攻角。O点是反应物DN的质心，由于N原子较重，它位于N原子一侧附近。θ为线DN与线HO之间的夹角，即攻击角。

目前的工作将试图通过量子波包计算来揭示标题反应及其同位素效应的交换反应动力学。此外，为了分析反应碰撞过程的机理，本文还在新的解析势能面上(4A'')进行了详细的QCT研究，揭示标题反应依赖攻击角的反应动力学。

2. 理论

在讨论计算结果之前，我们简要概述本工作中采用的QCT [13]和TDWP [14]方法的相关计算细节。

利用标准蒙特卡罗采样生成QCT计算的初始条件，并对Hamilton运动方程进行六阶辛积分三维积分。以0.1 eV为步长，在0 eV和1.0 eV之间连续变化碰撞能量，进行了计算。每个能量点总共运行10万条轨迹。为了比较，将反应物的初始旋转振动量子数v和j设为v = 0,1和j = 0。为了保证总能量和总角动量守恒，本文采用了0.1 fs的数值积分步长。所有的轨迹都是从A-BC核间距15 Å开始的，这个值被认为足够大，可以确保它们之间没有相互作用。我们使用公式$b=\sqrt{J\left(J+1\right)}\hslash /\mu \nu$ 来确定J依赖的反应几率，其中b为碰撞参数，J为总角动量，$\mu$ 和$\nu$ 分别为系统约化质量和相对速度。在本工作中，随碰撞能量和反应物振动激发演化的反应几率定义为P^J=0 (E_col, v) = N_R/N_T (NR表示总计算的反应轨迹数NT中的反应轨线数)。反应积分截面用式${\sigma }_R=\pi b_{\max }^2{N}_R/N_T$ 计算。定义抽样误差为$\sigma =\sqrt{\left(N_T-N_R\right)/N_T{N}_R}\times 100\%$ ，统计不确定度为${\sigma }_{rs}=\sigma {\sigma }_R$ 。

TDWP方法是一种强大的计算工具，可以提供精确的量子反应动力学，并已被许多研究小组广泛使用。本文还采用该方法计算了初始状态选择的总反应几率和总反应截面。量子计算中使用的参数如表1所示。

Table 1. Parameters for the quantum calculations

表1. 量子计算参数

对于原子-分子反应H + DN，用反应物雅可比坐标表示的哈密顿量可以写成

$\widehat{H}=-\frac{{\hslash }^2}{2{\mu }_R}\frac{{\partial }^2}{\partial R^2}+\frac{{\left(\widehat{J}-\widehat{j}\right)}^2}{2{\mu }_R{R}^2}+\frac{{\widehat{j}}^2}{2{\mu }_r{r}^2}+V\left(R,r,\gamma \right)+h\left(\widehat{r}\right)$ (1)

其中，R定义为H原子与DN质心之间的距离，R为DN键长，为向量R与r之间的夹角，${\mu }_R$ 和${\mu }_r$ 为系统和DN的约简质量。$\widehat{J}$ 和$\widehat{j}$ 为总角动量算符和DN旋转角动量算符，$V\left(R,r,\gamma \right)$ 为4A''状态下的相互作用势能面。反过来，双原子引用哈密顿量定义为

$h\left(\widehat{r}\right)=-\frac{{\hslash }^2}{2{\mu }_r}\frac{{\partial }^2}{\partial r^2}+V\left(\widehat{r}\right)$ (2)

式中$V\left(\widehat{r}\right)$ 为DN的双原子势函数。

本文采用分裂算子(split operator, SO)算法实现初始波包随时间在势能面上的传播。随时间变化的波函数可以用固定体(BF)平移–振动–旋转基$\left\{u_n^v\left(R\right){\varphi }_v\left(r\right)Y_{jK}^{JM\epsilon }\left(R,r\right)\right\}$ 展开为

${\psi }_{v_0{j}_0{K}_0}^{JM\epsilon }\left(R,r,t\right)={\displaystyle \sum _{n,v,j,K}{F}_{nvjK,v_0{j}_0{K}_0}^{JM\epsilon }\left(t\right)u_n^v\left(R\right){\varphi }_v\left(r\right)Y_{jK}^{JM\epsilon }\left(R,r\right)}$ (3)

其中v₀, j₀, K₀ (j₀的投影量子数)为表示初始转振态的量子数。n和v标记了平移和振动特征函数，M和K分别是总角动量J在空间固定z轴和BF z轴上的投影量子数。初始反应几率可以通过计算波函数含时部分经过与产物通道对应的足够长的时间后通过固定表面的总通量来提取。

$P_{v_0{j}_0{K}_0}^J\left(E\right)=\frac{\hslash }{{\mu }_r}\mathrm{Im}\left(\langle \psi \left(E\right)\left|\delta \left(r-r_0\right)\right|\psi \left(E\right)\rangle \right)$ (4)

一旦几率确定，总反应截面可以通过将所有部分波对反应几率的贡献加起来计算为

${\sigma }_{v_0{j}_0}\left(E\right)=\frac{\pi }{k_{v_0{j}_0}^2}{\displaystyle \sum _J\left(2J+1\right)P_{v_0{j}_0}^J\left(E\right)}$ (5)

式中$k_{v_0{j}_0}$ 为具有特定碰撞能量E的初始状态所对应的波数，为J和k₀的所有可能值的和。初始状态指定几率写为

$P_{v_0{j}_0}^J\left(E\right)=\frac{1}{2j_0+1}{\displaystyle \sum _{K_0}{P}_{v_0{j}_0{K}_0}^J\left(E\right)}$ (6)

本文中所有用TDWP方法进行的计算都采用了CS近似，并考虑了CC效应。

3. 结果与讨论

Figure 2. Comparison between the CC, CS and QCT probabilities

图2. CC，CS，QCT反应几率的比较

为了验证QCT计算的有效性，图2给出了在4A''势能面上0~1.0 eV碰撞能量范围内，总角动量J分别为0、20、30和40时，H + DN (v = 0, j = 0)→HD + N的QCT与QM结果的反应几率对比。观察到，在低碰撞能量下，J = 0时的CS几率略大于CC计算的几率，而在高能量下它们几乎保持相同的值。然而，对于更高的部分波，特别是在更高的碰撞能量下，CC和CS结果之间的差异越来越显著。这些结果表明，随着J和能量的增加，碰撞过程中的科里奥利耦合效应变得更加重要。另一个发现是，能量阈值随着J值的增加而向更高的值移动，并形成一个较长的阈值尾巴，这是由于离心势垒随着J值的增加而增加。此外，CS结果中发现了宽共振结构，并且随着J值的增加而逐渐消失，可能是由于寿命的缩短，在H + H2和H + HLi反应中也发现了这一点。这种现象的部分原因是PES的最小值随着J向更高的能量移动，与CS结果相比，CC共振结构更小，并向更高的碰撞能量移动。不同总角动量J = 0，20，30，40时，QCT计算的反应几率与CC计算结果非常吻合。这些计算结果证实NH2体系的4A''势能面非常适合量子散射和准经典计算。

结合QCT计算结果，用CC和CS方法对H + DN (v = 0, j = 0)→HD + N的反应截面进行了理论预测。这些结果如图3所示。很明显，在碰撞能量为0.1 eV以上，CS值总是大于CC值，并且这种趋势在一定程度上逐渐扩大。而在0.1 eV以下，CC和CS的结果表现出很好的一致性，随着碰撞能量的增加，科里奥利耦合效应更加重要。在整个碰撞能量范围内，QCT结果与CC结果的一致性较好。所有的曲线都表现出相似的阈值行为，并且随着碰撞能量的增加而单调增加，这与PES上沿最小能量路径的入口通道存在早期势垒所期望的行为一致。

Figure 3. Reaction cross sections for H + DN (v = 0, j = 0)→HD + N

图3. H + DN (v = 0, j = 0)→HD + N反应的积分截面

如前所述，依赖攻角的动力学包含了有关原子–双原子反应化学过程的丰富信息。图4为初始态(v = 0, j = 0)和图5为初始态(v = 1, j = 0)反应H + DN→HD + N在4A''面上的反应几率随碰撞能量和攻击角的变化曲线。很明显，反应物在基态和第一振动激发态的反应几率都强烈依赖于攻击角。这可以归因于在不同攻击角下，攻击角对壁垒高度的影响是不相同的。很明显，当攻击角大于50度时，反应被限制在一个有限的角度范围内，很难进行。从势能面上的特征可以看出，势垒随着攻击角的变化而变化，在小于50度的范围内，势能面上势垒非常低，H分子很容易跨过势垒发生反应。当攻击角大于50度时，势垒迅速升高，反应物没有足够的能量跨过势垒而发生反应，因此反应几乎难以进行。如上所述，本研究的碰撞能量不足以在攻击角50~180度范围内跨越壁垒并发生反应。从图5中可以看出，与图4相比，反应物振动激发对几率分布的影响是显著的。对于图4中H + DN (v = 0, j = 0)，反应倾向于主要发生在0~30度范围内，碰撞能量低于0.6 eV。随着碰撞能量的增加，这个角度范围最终转变为20~50度。而图5中H + DN (v = 1, j = 0)在碰撞能量为0.1 eV时，有效攻击角范围相对较宽，且随着碰撞能量的增加，有效攻击角范围急剧向30~50度移动，比(v = 0, j = 0)状态时更窄。从整个反应过程来看，上述特征表明H原子优先攻击反应物分子中D原子一侧，并立即生成产物HD。

Figure 4. Attack angle dependent reaction probabilities of reaction H + DN (v = 0, j = 0)→HD + N

图4. 攻击角依赖的H + DN (v = 0, j = 0)→HD + N反应的反应几率

Figure 5. Attack angle dependent reaction probabilities of reaction H + DN (v = 1, j = 0)→HD + N

图5. 攻击角依赖的H + DN (v = 1, j = 0)→HD + N反应的反应几率

图6和图7分别显示了在0~1.0 eV的碰撞能量范围内，D + HN→DH + N和D + DN→D2 + N反应在选定J值下指定初始状态(v = 0, j = 0)的QCT和QM几率。可以看出，当J值较低时，如J = 10，CC几率与CS几率之间的差异很小，但当J变大时，如J = 40，50，60，则差异变大。尤其需要注意，在所有J值情况下，在整个能量范围内，CC值总是小于CS值。在J和碰撞能量固定的情况下，D + DN的碰撞参数虽然小于D + HN，但D + HN的反应几率大于D + DN。

Figure 6. Reaction probabilities for D + HN (v = 0, j = 0)→DH + N with selected J values

图6. 不同J值下D + HN (v = 0, j = 0)→DH + N的反应几率

Figure 7. Reaction probabilities for D + DN (v = 0, j = 0)→D2 + N with selected J values

图7. 不同J值下D + DN (v = 0, j = 0)→D2 + N的反应几率

4. 结论

我们利用含时波包方法和QCT方法在一个新的势能面上对H + DN反应及其同位素变化进行了动力学研究。我们计算了碰撞能量达到1.0 eV时标题反应的反应几率、反应截面及其同位素变化。计算结果显示反应几率和截面随着碰撞能的增加而逐渐增大，这符合反应入口存在早期势垒的势能面特征。同时反应物的振动激发对反应的影响是正面的，能促进反应的进行。此外，由于不同几何构型下的反应物分子对应的势能不同，导致在不同攻击角下反应几率有很大的差异。

参考文献