H2+NH2-反应在全维精确势能面上的动力学模拟：验证和解释实验

on style="white-space: normal; line-height: 1.75em; box-sizing: border-box;">

研究背景

离子-分子反应在天体化学等涉及到离子的气氛中都扮演着重要的角色。实验上，由于离子的库伦相互作用，很难制备出高质量的离子束源。理论上，对于离子-分子反应体系构建全维高精度势能面比较困难，长程相互作用导致采样点数量巨大。而H₂+ NH₂^-作为一个典型的多通道离子-分子反应，实验和理论计算方面的报道都较少。

H₂+ NH₂^- ↔ H^- + NH₃ (R1)

实验上：

Bohme等人测量了297 K下R1反应正、逆向的速率常数以及平衡常数。Otto等人(Phys. Rev. Lett., 2008, 101, 063201)测量了R1反应在8-300K之间的速率常数。他们发现在8-20 K之间，速率常数随温度增加，在20-300K之间，速率常数随温度增加而减小。

理论计算：

Gianturco等人引入一个可调节的温度依赖项，采用变分过渡态方法计算了R1反应的速率常数，与实验值十分吻合。由于动力学作用可能导致统计模型失效。而准经典轨线(QCT)和量子动力学(QD)计算依赖于全维高精度势能面。于是我的导师李军就建议我去构建这个反应体系的全维高精度势能面，并进行速率计算、动力学研究等。巧合的是，中科院武汉数学与物理研究所的宋宏伟老师也在进行相应的研究，我们互相都不知情。在本工作即将完成的时候，他们发表了文章(Phys. Chem. Chem. Phys., 2021, 23, 17848)。宋老师采用CCSD(T)-F12a/AVTZ计算了大约46000个点，通过基础不变量-神经网络(FI-NN)方法拟合了势能面，并得到了速率常数，比实验值偏高。经过检查和对比发现FI-NN势能面只覆盖了能量较低的反应动力学区域TS1、TS2、反应物、产物，如Fig. 1所示。该势能面没有包括TS3和NH₄^-所在区域，并且采用了不同的计算水平。我们决定将课题继续进行下去。

▲Fig. 1 Schematic illustration of the energetics along the reaction H₂ + NH₂⁻. All energies are in kcal mol⁻¹ and relative to the H₂ + NH₂⁻ asymptote at various levels: PIP-NN PES, UCCSD(T)/AVTZ′, CCSD(T)/AVTZ′, and CCSD(T)-F12a/AVTZ, FI-NN PES from top to bottom.

02

构建全维势能面

首先确定从头算方法和基组，它们决定了数据质量，进而决定了势能面的精度和可靠性。TS3所在反应通道多参考性较强，采用多参考方法MRCI方法可能更好。但MRCI计算量太大，对于数万个高精度点的计算来说太昂贵。李军教授之前已经做过关于R1反应的产物中间体和NH₄^-的局部势能面，经过测试，发现UCCSD(T)-F12a/AVTZ对于垂直解离能描述不好(J. Chem. Phys., 2016, 144, 244311)。需要对N原子添加一组额外的弥散基组，标记为AVTZ’。我们还做了多个取向的一维扫描，发现单参考方法UCCSD(T)对于该体系可以描述的很好。最后采用UCCSD(T) /AVTZ’进行计算。其次是在构型空间采样和采用PIP-NN方法拟合势能面，详见文章和支撑材料的介绍。最终，我们选取了98154个数据点来保证反应物或产物之间的长程相互作用得到精确描述，见Fig. 2和3，拟合总误差仅为0.026 kcal/mol。

▲Fig. 2 (a)/(b) Distribution of the sampled points on the PIP-NN PES with the two reactive bond lengths up to 50 Å. The inset of (a) shows fitting errors (Efit − Etarget, in kcal mol−1) for the PIP-NN PES as a function of the ab initio energy up to 100.0 kcal mol−1. The inset of (b) displays the contour plot of the potential energies (in kcal mol-1) for the process between H2 + NH2- and NH3 + H- as a function of two reactive bond lengths, whose definitions are also shown. All other coordinates are fixed at TS1.

▲Fig. 3 Potential energies along the MEPs for the transformations between (a) NH2-‧‧‧H2 and NH3‧‧‧H-, and (b) NH3‧‧‧H- and NH3‧‧‧H-, as a function of the reaction coordinate s (amu1/2bohr). One-dimensional cuts for the potential energy between NH2- and H2 (c) and between NH3 and H- (d) along the R distance. The symbols represent the ab initio energies calculated at the UCCSD(T)/AVTZ′ level.

03

速率常数

采用准经典轨线(QCT)方法在50K、100K、200K和300K共五个温度下分别计算了2-25万条轨线，统计误差均小于5%。由于生成NH₄^-的反应能垒较高，在计算中并没有发现生成NH₄^-的轨线，故只考虑两个夺氢反应通道。如Fig. 4所示，QCT计算结果比实验值高，这主要是因为QCT方法不能考虑量子效应，如隧穿、零点振动能等。其中零点能泄露导致计算反应几率偏高，速率常数偏高。为了部分矫正该效应，采用hard-ZPE方法，即如果产物NH₃的振动能小于其零点振动能，则不予统计该轨线。可以看出hard-ZPE的结果更接近实验值，尤其是同位素效应。但计算的速率常数还是偏高，因为难以判断轨线在传播过程中是否发生了零点能泄露。

▲Fig. 4 Theoretical and experimental thermal rate coefficients for the H2 + NH2− → NH3 + H− reaction (a) and the D₂ + NH₂⁻ → NH₂D + D⁻ reaction (b). The KIEs are plotted in (c). Available results from ref. 1, 4, 7 and 11 were included for comparison.

04

反应动态学

在Otto等人的实验中，H₂(ortho)和H₂(para)的比例是3：1。我们就直接来检查H₂(ortho)和H₂(para)对反应的影响。分别在1.0、2.0、5.0、8.0、10.0、15.0、20.0 kcal/mol 7个碰撞能，以及H₂(j=0, 1, 2)三种转动态下计算了1-70万条轨线。所有的统计误差均小于5%。如Fig. 5(a)所示，QCT的结果显示，H₂的转动激发明显抑制了反应的发生，该现象与以往离子-中性反应不同。在之前的H₂+OH⁺和H₂+H₂O⁺两个反应体系中，H₂的转动激发对反应起到促进作用。由于QCT方法不能考虑量子效应，中科院武汉物理与数学研究所的宋宏伟老师进行QD计算。在低碰撞能下，QD计算结果远小于QCT的计算结果，由于QCT方法无法考虑零点能效应。对于H₂(j=1)，QD结果与QCT结果一致。对于H₂(j=2)，QD计算过于昂贵，预期H₂(j=2)的QD结果跟QCT结果一致。H₂的转动激发抑制反应的一部分原因来自于H₂和NH₂^-之间的相互作用势能面。如Fig. 6所示，只有取向g和h为吸引作用，其余均为排斥作用。而且取向h的吸引作用较强，利于反应前体的形成。因此，H₂的转动激发使得H₂和NH₂^-在相互靠近过程中更容易形成排斥作用的取向，从而抑制了反应的发生。

▲Fig.5 (a) Integral cross-sections (ICSs) for the reaction NH₂⁻ + H₂ (j = 0, 1, 2) → NH₃ + H⁻ on the new PIP-NN PES as a function of the collision energy (in kcal mol−1). (b) K-specific ICSs calculated by QD.

▲Fig. 6 Asymptotic potential energy curves between NH₂⁻ and H₂ at different orientations. Both NH₂⁻ and H₂ are fixed at equilibrium. The lines represent PIP-NN PES and symbols for ab initio energies.

05

结论与展望

采用准经典轨线方法在新构建的全维高精度势能面上进行了速率常数计算及动态学计算。通过QCT计算发现H₂的转动激发抑制了反应的发生，并通过QD计算进行了验证。这种特殊的反应机理丰富了我们对于离子-中性反应体系微观反应动力学机制的认识。除此之外，我们在势能面上进行了H₂/D₂+NH₂⁻反应速率常数的计算，计算结果与温度呈负相关，与实验相比，变化趋势相同但比实验值高。这主要是由于QCT方法不能考虑量子效应。未来希望可以找到合适的计算方法重现实验上低温(<20K)速率常数随温度的变化趋势。