Neural network potential energy surface and dynamical isotope effects for the N^{+}(^{3}P) + H_{2} → NH^{+} + H reaction
Received
17th May 2019
, Accepted 6th August 2019
First published on 7th August 2019
The N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) → NH^{+}(X^{2}Π) + H(^{2}S) reaction is important for initiating the chain reaction of ammonia synthesis in the universe. To study the dynamics of this reaction, a global accurate potential energy surface (PES) of the ground state NH^{+}_{2} was constructed by combining numerous highlevel ab initio energy points with the permutation invariant polynomial neural network method. Utilizing this newly constructed PES, timedependent wave packet calculations were performed on the statetostate reactions of N^{+}(^{3}P_{0}) + H_{2} (v = 0, j = 0) and N^{+}(^{3}P_{0}) + D_{2} (v = 0, j = 0) in order to study the microscopic reaction mechanisms and dynamical isotope effects. Isotope effects have a significant influence on the rovibrational state distributions and state resolved angle distributions of the product. The total differential cross sections (DCSs) present in the aforementioned reactions are dominated by both forward and backward scattering, as expected from the observable deep well along the reaction path. Meanwhile, the rovibrational stateresolved DCSs show that both reactions are not entirely statistical at the statetostate level.
1 Introduction
The nitrogen atom and nitrogenous molecules are important tracers when studying interstellar objects; nitrogen hydrides or their ions are also of significance in astrophysics due to their abundance in the universe. One of the simplest nitrogen chemical reactions of N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) → NH^{+}(X^{2}Π) + H(^{2}S) has been suggested as the first step in the synthesis of ammonia within dense interstellar clouds.^{1,2} In addition, as a prototype for understanding ion–molecule reaction dynamics, it has highly interesting characteristics for study because of its near thermoneutrality. Especially, the dynamical behaviors are substantially different when the hydrogen undergoes isotopic substitution; numerous studies^{3–12} have focused on these isotope effects, which have been consecutively explored over the past 30 years.
A variety of experimental results for the title reaction and its isotopic counterparts, such as velocityangle distributions^{13,14} and integral cross sections (ICSs),^{4,9} have been reported and provide abundant reference data for theoretical studies. Additionally, according to reported ab initio calculations,^{15} the title reaction processes a deep potential well (∼6.5 eV) path, resulting in the formation of longlived complex NH^{+}_{2} at low energies, and only a limited number of states are included due to its slight endothermicity. Thus, various statistical methods^{11,12,16,17} have been employed to treat the reactions of N^{+}(^{3}P) with H_{2} and its isotopic variants. Gerlich^{16} applied dynamically biased statistical theory to obtain the statespecific ICSs for the title reaction. The calculated results were very comparable with known experiments; however, large discrepancies occurred upon replacing H_{2} molecule with D_{2} molecule. Recently, Grozdanov et al.^{11} formulated a new statistical model to calculate the ICSs of reactions involving N^{+} with H_{2}, D_{2} and HD molecules. The fine structure states of N^{+} ion were treated equally to other internal motions, and the calculated results for the three reactions could reproduce the experimental ICSs. In later work,^{12} the authors concluded that the dynamics results were strongly dependent on the fine structure state of N^{+} ion and the initial rotation state of the reactant H_{2} or HD.
However, theoretical treatment of statistical models is limited due to the addition of certain assumptions. A better method for accurately describing the characteristics of the title reaction is based on the precise potential energy surface (PES) that covers the entire reactant and product regions. The first ground state PES of NH^{+}_{2} was constructed by Gittins and Hirst.^{18} The energy points were calculated via the configuration interaction (CI) method, and a deep ^{3}B_{1} well at C_{2v} symmetry was represented in the PES. However, a crucial issue was that their calculations of the product generated the wrong electronic state of NH^{+}(a^{4}Σ^{−}), which is actually the first excited state of NH^{+} molecule.^{19} Gonzalez et al.^{20} presented an ioninduced single surface model for the title reaction and discovered that the ^{3}A_{2} well plays a predominant role in the energy range of 1 to 4 eV. The most recent PES for the NH^{+}_{2} system was reported by Wilhelmsson et al.^{21} in 1992. In their work, a total of 752 configurations calculated by the complete active space selfconsistent (CASSCF) and multireference configuration interaction (MRCI) methods were used to fit the analytical form, and the average leastsquares deviation was 0.026 eV. Based on these PESs, dynamics simulations of the N^{+} and H_{2} reaction and its isotopic variants using classical trajectory calculations^{22} or the quasiclassical trajectory (QCT) method^{6,7,20,23} have been extensively implemented.
Although the QCT method or statistical models can describe the mechanism of the reactions of N^{+}(^{3}P) with H_{2} and its isotopic variants to some extent, the significant quantum effects are not included. The structure of the deep well and its nearly thermoneutral character imply that the quantum scattering resonances and the zero point energy (ZPE) effect are highly crucial to the reaction dynamics. Performing quantum wave packet calculations on the title reaction is extremely expensive due to the numerous quasibound states supposed by the deep well, which require a large number of grids and a long propagation time in order to achieve evolution of the wave packet. It is well known that the only quantum dynamics results for the N^{+}(^{3}P) + H_{2} reaction were calculated using the timedependent wave packet (TDWP) method proposed by Russell and Manolopoulos.^{24} They performed TDWP calculations for the reactions of N^{+}(^{3}P) + H_{2} (v = 0 and 1, j = 0 and 1) on the Wilhelmsson PES. The ICS results of the ground vibrational reactions showed different behaviors compared with QCT calculations due to the inability of the traditional QCT method to treat the vibrational ZPE effect of the product molecule. In the QCT calculations, only the ZPE of the reactant H_{2} molecule was corrected; thus, the title reaction became exothermic.
According to previous studies, traditional QCT calculations are unable to correctly describe the dynamical behaviors of the title reaction. The ZPE problem can be solved to some extent using various corrective methods; however, the Gaussian–Binning method displays excellent corrective abilities toward this shortcoming of the QCT method. The Gaussian–Binning procedure^{25,26} has been successfully incorporated into several barrierless reactions, such as H^{+} + H_{2}^{27,28} and C(^{1}D) + H_{2}^{29} reaction systems, and the calculated results were in good agreement with the quantum dynamics calculations. Another basic problem of the QCT method is its poor capability to describe quantum resonances, which are vital for deep well reactions. Therefore, using only quantum methods can give perfect dynamics information for the title reaction. Although quantum dynamics calculations for the title reaction have been reported, the inaccuracy of the results is attributed to the lack of reliable PES. The error of the Wilhelmsson PES is relatively large due to the limitations of early computational conditions. In addition, ab initio calculations^{30} and a recent experimental study^{10} have proved that the endothermicity of the reaction (33 meV) on the Wilhelmsson PES is too large. Hence, performing quantum dynamics calculations for the title reaction on a highprecision global PES is necessary to improve the accuracy of the dynamics results. Furthermore, quantum dynamics simulations of the N^{+}(^{3}P) + D_{2} reaction have not been presented to date, and it is desirable to perform these simulations to study the influence of isotope effects on the title reaction. Herein, a new accurate ground state PES of the NH^{+}_{2} system is constructed using a large number of high level ab initio energy points and a neural network (NN) model. Based on our new PES, rigorous quantum TDWP calculations for the N^{+}(^{3}P_{0}) + H_{2} (v = 0, j = 0) → NH^{+} + H and N^{+}(^{3}P_{0}) + D_{2} (v = 0, j = 0) → ND^{+} + D reactions are implemented to study the reaction mechanism and the influence of isotope effects on the dynamics results.
2 Potential energy surface
2.1
Ab initio calculations
The electronic energies of the ground state (1^{3}A′′) NH^{+}_{2} within C_{s} symmetry were calculated by the internally contracted MRCI method^{31,32} with the Davidson correction (icMRCI + Q), as implemented in the MOLPRO package.^{33} First, the Hartree–Fock wave functions were calculated as an initial guess, followed by CASSCF^{34,35} calculations for the seven equally weighted states (1^{3}A′′, 2^{3}A′′, 3^{3}A′′, 1^{3}A′, 1^{1}A′, 2^{1}A′, 1^{1}A′′). The icMRCI + Q calculations were carried out based on the optimized wave functions of the CASSCF process. In the CASSCF and icMRCI + Q calculations, six valence electrons are included in seven active orbitals (6a′ + 1a′′). The augmented correlationconsistent polarization valence quadruplezeta (AVQZ) basis set^{36} and the Dunningweighted correlation consistent polarized corevalence triplezeta (WCVTZ)^{37} basis set were used for H and N atoms, respectively, which ensured accuracy of the calculations and convergence of key energy points by testing numerous geometries. A large number of ab initio energies over a great range of configuration spaces were chosen in the Jacobi coordinates, where the N^{+}–H_{2} region was defined by 0.9 ≤ r_{HH}/a_{0} ≤ 25, 0.8 ≤ R_{N+–HH}/a_{0} ≤ 25, and 0 ≤ θ/degree ≤ 90 and the H–NH^{+} region was defined by 0.8 ≤ r^{+}_{NH}/a_{0} ≤ 28, 0 ≤ R_{H–NH+}/a_{0} ≤ 25, and 0 ≤ θ/degree ≤ 180.
2.2 Fitting the potential energy surface
Representing global PESs from ab initio energies is challenging but can be achieved by combining electronic structure calculations with advanced machine learning technologies. Artificial NNs and Gaussian process regression are the most popular machine learning methods to fit PESs by a large number of ab initio data; these methods have been widely applied to numerous molecular systems.^{38–46} In this work, the NN method was used to obtain the analytical PES of the NH^{+}_{2} system. The expression of the global surface consists of twobody and threebody terms and is written as follows: 
 (1) 
where R_{1}, R_{2} and R_{3} are the internuclear distances of NH^{+}_{a}, HH and NH^{+}_{b}, respectively. For the ion–molecule reaction, the accuracy of the longrange potential has substantial influence on the dynamics results; thus, a switch function f(R_{1}, R_{2}, R_{3}) was used to better describe the PES in the longrange region, and it is written as follows: 
 (2) 
where R_{n} represents the bond length of two atoms and R_{d} and R_{w} are the central position of f(R_{1}, R_{2}, R_{3}) and the switch strength constant, respectively. In the process of twobody term fitting, the NN structure included two hidden layers, and 6 neurons were included in each hidden layer. 97 and 92 ab initio points were generated to obtain the potential energy curves of HH and NH^{+}, and the root mean square errors (RMSE) are 0.10 and 0.12 meV, respectively. Fig. 1 shows the potential energy curves (PECs) of the H_{2}(X^{1}Σ^{+}_{g}) and NH^{+}(X^{2}Π) molecules, where the fitting of the PECs is in good agreement with the original ab initio data. The equilibrium bond length and dissociation energy of H_{2}(X^{1}Σ^{+}_{g}) on the fitted PECs are 1.402 Bohr and 4.71 eV, respectively, which compare well with the corresponding experimental data^{47} of R_{e} = 1.401 Bohr and D_{e} = 4.74 eV. In the case of NH^{+}(X^{2}Π) molecule, the equilibrium bond length and dissociation energy are 2.047 Bohr and 4.64 eV, respectively, showing good agreement with the experimental results^{48} of R_{e} = 2.043 Bohr and D_{e} = 4.65 ± 0.5 eV. Thus, the fitting results for the twobody potential are sufficiently accurate.

 Fig. 1 Comparison of ab initio data and fitting results for the potential energy curves of H_{2}(X^{1}Σ^{+}_{g}) and NH^{+}(X^{2}Π).  
For the threebody term, permutation invariant polynomials (PIP)^{49,50} were employed to solve the problem of the adaptation of permutation symmetry caused by the two identical H atoms. The geometries after PIP preprocessing were used in the input layer of the NN. First, the bond length between two atoms was transformed as

 (3) 

 (4) 

P_{HH} = exp(−0.2R_{HH})  (5) 
Then, the following PIPs were adopted:

 (6) 

 (7) 
Finally,
G_{1},
G_{2}, and
G_{3} were given normalized treatment as

 (9) 
where
G_{i,min} and
G_{i,max} are the minimum and maximum values of
G_{i}, respectively. The NN included two hidden layers with 12 neurons, and the Levenberg–Marquardt algorithm
^{51} was applied to train the NN. The final threebody NN expansion is presented as

 (10) 
where
ω and
b represent the corresponding weights and bias, which were optimized by training the NN. The linear and hyperbolic tangent functions were used as the transfer functions
ϕ in the output layer and hidden layers, respectively. To avoid overfitting, 30
519 energy points were divided randomly into three sets, namely the training set (90 percent), testing set (5 percent), and validation set (5 percent). The overall RMSE of the PIPNN PES is only 2.7 meV.
2.3 Features of the potential energy surface
Fig. 2 shows the potential energies at four different insertion angles in the reactant Jacobi coordinates, where r_{HH} and R_{N+–HH} are the HH bond length and the distance between N^{+} and the center of mass of the H_{2} molecule, respectively, and θ is the angle between r_{HH} and R_{N+–HH}. The energies are relative to the N^{+} + H_{2} dissociation limit, and the maximum energy was set as 5 eV in order to clearly display the structures of the wells. There are two wells at every angle, which correspond to the ^{3}B_{1} and ^{3}A_{2} wells at the C_{2v} structure, respectively. As the insertion angle decreases, the deeper well corresponding to the ^{3}B_{1} state at the C_{2v} structure becomes less accessible and the depth of the outer ^{3}A_{2} well slowly decreases. The ^{3}B_{1} well is located at r_{HH} = 3.76 Bohr at the Tshaped symmetry, indicating the involvement of a bond stretching process prior to complex NH^{+}_{2} formation. Moreover, there is a saddle point between the two wells, which can be attributed to prevention of the crossing of the two lowest electronic states.^{20} In the case of C_{2v} symmetry, the complex NH^{+}_{2} processes form at the deep ^{3}B_{1} state well via the ^{3}A_{2} → ^{3}B_{1} intersection, and the two lowest states ^{3}A_{2} and ^{3}B_{1} return to the adiabatic ground state ^{3}A′′ in the C_{s} geometry. In the present work, the ground state PES is constructed in the adiabatic representation, and the conical intersection is hampered by the limitation of the Born–Oppenheimer approximation. Thus, the origin of the repulsive regions is the interaction of the first excited state.

 Fig. 2 Potential energy surfaces of the ground state NH^{+}_{2} for insertion approaches at four θ angles (0°, 30°, 60° and 90°).  
The geometries and energies of the stationary points obtained on the new PES and original ab initio calculations are listed in Table 1 and compared with reported results. The energy values are relative to the N^{+}(^{3}P) + H_{2} dissociation limit. For the C_{2v} geometry, the equilibrium structure corresponds with the ^{3}B_{1} well, which is a highly crucial structure to the title reaction. As shown, the results calculated on the new PES are in good agreement with the original and previous ab initio calculations,^{15} and the Wilhelmsson PES underestimates the depth of the well. However, there is a large deviation between our data and previous ab initio calculations^{52} for the minimal energy of the collinear structure. This may be because the electron correlation energy cannot be treated well using the selfconsistent field method with rough CI corrections. The new PES can better reproduce the highlevel ab initio calculations at the stationary points, suggesting enhanced reliability in reaction dynamics studies.
Table 1 Stationary points of the ground state NH^{+}_{2}
C
_{2v}

r
_{HH} (Bohr) 
R
_{N+–HH} (Bohr) 
Energy (eV) 
NN PES 
3.76 
0.48 
−6.49 
Wilhelmsson PES^{21} 
3.66 
0.64 
−6.15 
icMRCI + Q 
3.76 
0.48 
−6.50 
MRDCI^{15} 
3.67 
0.49 
−6.45 
C
_{∞v}

R
_{HH} (Bohr) 
R
^{+}_{NH} (Bohr) 
Energy (eV) 
NN PES 
1.68 
2.24 
−2.20 
Wilhelmsson PES^{21} 
1.70 
2.27 
−2.05 
icMRCI + Q 
1.67 
2.24 
−2.22 
SCFCI^{52} 
1.60 
2.50 
−1.76 
The contour maps of the new PES at four different N^{+}–H–H angles (45°, 90°, 135° and 180°) are presented in Fig. 3. R_{HH} represents the bond length of the two H atoms, R^{+}_{NH} is the distance between the N^{+} and one of the H atoms, and the N^{+}–H–H angle is the included angle between R_{HH} and R^{+}_{NH}. The energy is set as zero at the N^{+} + H_{2} dissociation limit. Two valleys exist at the bottom and left, corresponding to the N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) channel and NH^{+}(X^{2}Π) + H(^{2}S) channel, respectively. The depth of the two valleys is comparable, indicating that the title reaction is approximately thermoneutral. The deep well structure connects the two valleys at all four approaching angles, highlighting the important role of quantum resonances in the reaction dynamics. At the beginning of the title reaction, the N^{+} ion and H_{2} molecule convene in the bottom valley, generating the longlived NH^{+}_{2} complex in the deep well and eventually dissociating to the NH^{+} molecule and H atom in the left valley. To more clearly describe the characteristics of the title reaction on the new PES, the minimum energy paths (MEPs) of the N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) → NH^{+}(X^{2}Π) + H(^{2}S) reaction at four N^{+}–H–H angles (45°, 90°, 135°, 180°) are shown in Fig. 4. The MEPs were determined by scanning the PES, as shown in Fig. 3, at different coordinates (R_{HH}–R^{+}_{NH}) to obtain the minimum values of the potential energy. The depth of the well decreases with increasing S–H–H angle, suggesting that the NH^{+}_{2} complex has a shorter lifetime when the reaction proceeds at a large approaching angle. It can be seen from the minimum energy paths that the title reaction is weekly endothermic, with an endothermicity of approximately 83 meV; this is consistent with the ab initio calculations. Table 2 lists the molecular constants of the reactants and products, which were obtained from the analytical PES. When the ZPEs of the reactants and products are considered, the endothermicities of the N^{+} + H_{2} and N^{+} + D_{2} reactions are approximately 8 and 23 meV, respectively, which is in good agreement with the experimental measurements^{8} of 11 ± 3 and 29 ± 3 meV.

 Fig. 3 Potential energy surfaces of the ground state NH^{+}_{2} at four N^{+}–H–H angles (45°, 90°, 135° and 180°).  

 Fig. 4 Minimum energy paths for the N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) → NH^{+}(X^{2}Π) + H(^{2}S) reaction at four N^{+}–H–H angles (45°, 90°, 135° and 180°).  
Table 2 Molecular constants of the reactants and products

ω
_{e} (cm^{−1}) 
ω
_{e}
x
_{e} (cm^{−1}) 
D
_{0} (eV) 
ZPE^{a} (eV) 
ZPE determined by ω_{e}/2 − ω_{e}x_{e}/4.

H_{2}(X^{1}Σ^{+}_{g}) 
4398.2 
121.3 
4.445 
0.269 
D_{2}(X^{1}Σ^{+}_{g}) 
3109.1 
60.2 
4.523 
0.190 
NH^{+}(X^{2}Π) 
2870.3 
66.5 
4.463 
0.176 
ND^{+}(X^{2}Π) 
2062.1 
38.8 
4.512 
0.127 
Fig. 5(a) displays the potential energy of N^{+} ion moving around the H_{2} molecule with a fixed bond length of H_{2} (1.401 Bohr); the energy is set as zero when the N^{+} ion is far from the H_{2} molecule. A 2.67 eV deep well is present at x = 0.0 Bohr and y = 2.3 Bohr. The depth of the well increases as the HH bond stretches, which can promote insertion collisions at low energies as the N^{+} ion slowly approaches the H_{2} molecule. Similar contours to Fig. 5(a) for a H atom moving around the NH^{+} molecule are presented in Fig. 5(b). The NH^{+} bond length is set to the equilibrium distance (2.043 Bohr), and the energy is set as zero when the H atom is far from the NH^{+} molecule. It can be seen that a deeper well exists around the N^{+} ion, implying that the H atom can more readily escape from the side of the H in the NH^{+} molecule in the title reaction.

 Fig. 5 (a) Contour plot of the PES for the N^{+} ion moving around the H_{2} molecule fixed at the equilibrium distance (1.401 Bohr). (b) Contour plot of the PES for a H atom moving around the NH^{+} molecule fixed at the equilibrium distance (2.043 Bohr).  
3 Dynamics calculations
The TDWP method is a powerful tool for calculating atom–diatom reactive scattering and can give accurate information about the statetostate reaction dynamics. Only an abbreviated outline is introduced here; more detailed descriptions of the method can be found in previous literature reports.^{53,54} In the body fixed representation, the Hamiltonian is expressed as follows: 
 (11) 
where R is the distance from the N^{+} ion to the center of mass of the H_{2} molecule and r is the bond length of the H_{2} molecule. μ_{R} and μ_{r} are the corresponding reduced masses associated with the R and r coordinates. Ĵ and ĵ are the total angular momentum operator of the NH^{+}_{2} system and the rotational angular momentum operator of the H_{2} molecule, respectively. is the potential energy of the NH^{+}_{2} system which is obtained on the NN PES. The statetostate Smatrix was extracted by the reactant coordinatebased method,^{55,56} and the secondorder split operator was adopted to evaluate the wave function. The statetostate reaction probability was calculated by 
 (12) 
The stateresolved ICSs were obtained by 
 (13) 
where k_{v0j0} is the momenta in the entrance channel. The stateresolved differential cross sections (DCSs) can be calculated by 
 (14) 
in which θ is the scattering angle and represents the element of the reduced Wigner rotation matrix.
The initial rovibrational states of the reactant H_{2} molecule were set as v_{0} = 0 and j_{0} = 0, and the TDWP calculations included all Coriolis couplings. All partial waves were calculated up to J = 75 for the N + H_{2} reaction and J = 81 for the N + D_{2} reaction, yielding converged ICSs and DCSs for collision energies beyond 0.7 eV, which was set as the upper limit of the collision energy in the reaction dynamics. Large grids and sufficient propagation times were used to ensure that the wave packet could reach the low energy region. The optimal numerical parameters used in the calculations are listed in Table 3 by numerous convergence tests.
Table 3 Numerical parameters used in the TDWP calculations

N^{+}(^{3}P) + H_{2} 
N^{+}(^{3}P) + D_{2} 
Grid ranges and sizes 
R (Bohr) ∈ [0.1, 22.0], N^{tot}_{R} = 389 
R (Bohr) ∈ [0.1, 22.0], N^{tot}_{R} = 389 
r (Bohr) ∈ [0.4, 16.0], N^{tot}_{r} = 299 
r (Bohr) ∈ [0.4, 16.0], N^{tot}_{r} = 299 
N
_{
j
} = 100 
N
_{
j
} = 100 
Initial wave packet 
R
_{c} = 16.0 Bohr 
R
_{c} = 16.0 Bohr 
Δ
_{
R
} = 0.20 Bohr 
Δ
_{
R
} = 0.18 Bohr 
k
_{0} = (2E_{0}μ_{R})^{1/2} with E_{0} = 0.50 eV 
E
_{0} = 0.50 eV 
Propagation time 
100000 a.u. 
100000 a.u. 
Highest J value 
74 
81 
The collision energy dependence of the total reaction probabilities for the N^{+} + H_{2} and N^{+} + D_{2} reactions with four different partial waves are presented in Fig. 6. An extremely dense oscillation feature was found in the reaction probability curves, which can be attributed to the formation of longlived complexes in the deep well. The domination of longlived resonance is similar to the typical insertion reactions of S(^{1}D) + H_{2}^{57} and C(^{1}D) + H_{2},^{58} and the title reaction has a sharper oscillation peak because of the deeper well. For J = 0, the threshold of the two reactions is consistent with the corresponding endothermicity, implying a barrierless path at low collision energies. With increasing J, the oscillation amplitude becomes less pronounced due to the increased centrifugal barrier. In the case of high J values, the threshold of the N^{+} + H_{2} reaction is larger than that of the N^{+} + D_{2} reaction due to the heavier reduced mass of ND^{+}_{2}. Because the increased collision energy diminishes the constraint of the deep well, the resonance peak becomes broader at high energy.

 Fig. 6 Total reaction probabilities for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions as a function of collision energy with four different total angular momentums (J = 0, 20, 40, 60).  
The total ICSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions obtained from the TDWP calculations on the NN PES together with the experimental measurements^{9} using a guided ion bean mass spectrometer and calculated by statistical theory^{11} are shown in Fig. 7. Previous TDWP calculations^{24} for the N^{+}(^{3}P) + H_{2} reaction on the Wilhelmsson PES are also presented. The ICS values for the N^{+}(^{3}P) + H_{2} reaction calculated by the TDWP method on the NN PES are larger than those on the Wilhelmsson PES due to the smaller reaction threshold of the NN PES. As shown, the quantum resonances remain strong on the TDWP ICS curves even though various partial waves have been summed, and the oscillatory structures calculated on the NN PES are narrower than the previous results due to the deeper well of the NN PES. The TDWP results show that the total ICS increases slightly as a function of collision energy and then nearly reaches a plateau, which is characteristic of an endothermic reaction.

 Fig. 7 Total ICSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions as a function of collision energy compared with previous TDWP calculations on the Wilhelmsson PES, the results calculated by statistical theory and available experimental values.  
In general, the ICS values generated by the TDWP calculations are close to the results calculated by the statistical method and the experimental data. However, the collision energy dependence of the ICSs obtained by the TDWP calculations is in contrast with the statistical and experimental results, which present a monotonic decreasing dependence of collision energy. This arises because the statistical and experimental data include the excited spin–orbit states of N^{+} ion and the higher rotational states of the H_{2} or D_{2} molecule. In the previous experiment, the rotational temperature of H_{2} or D_{2} was 105 K, and the N^{+} ion was produced in a thermal distribution at 300 K. In this case, the rotational state populations of H_{2} at j = 0, 1, 2, and 3 were 0.243, 0.749, 0.008 and 0.0004, respectively, and the rotational state populations of D_{2} at j = 0, 1, 2, 3, and 4 were 0.472, 0.322, 0.193, 0.012 and 0.001, respectively. The spin–orbit state distributions for N^{+} ion at ^{3}P_{0}, ^{3}P_{1} and ^{3}P_{2} were 0.17, 0.39 and 0.44, respectively. In the statistical theory results, the spin–orbit states of N^{+} ion were treated equally to the vibrational and rotational states of H_{2} and D_{2}, and the endoergicities of the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions were set as 18.45 and 46.32 meV, respectively. The electronic energy of N^{+} ion is 16.2 meV for the fine structure state ^{3}P_{2}, and the rotational energy of H_{2} molecule at j = 1 is 14.4 meV, which are larger than the endoergicity of the title reaction. The averaging effects of the excited spin–orbit states and higher rotational states allow the two reactions to become exothermic. It is clear that the statistical theory results successfully reproduce the experimental data of the N^{+}(^{3}P) + H_{2} reaction even when some assumptions are included at the energy range. However, in the case of the N^{+}(^{3}P) + D_{2} reaction, the statistical theory results and experimental data have relatively large deviations. In addition, the statistical method cannot characterize the quantum resonances, which is crucial for a deep well reaction. Herein, only rigorous quantum dynamics calculations on the ground adiabatic PES for the stateresolved reactions of N^{+}(^{3}P_{0}) + H_{2} and N^{+}(^{3}P_{0}) + D_{2} were performed; thus, they displayed different dynamics behavior compared with the statistical and experimental results.
The rovibrational stateresolved ICSs of the product NH^{+} and ND^{+} molecules at four collision energies (0.1 eV, 0.3 eV, 0.5 eV and 0.7 eV) calculated by the TDWP method are shown in Fig. 8. At lower collision energies, only the ground vibrational state is populated, and the v′ = 1 and v′ = 2 channels are opened with increasing collision energy. As shown, the products are difficult to excite to higher vibrational states, implying that the energy conversion efficiency is low. The main reason for this is the presence of a deep well on the reaction path: the N^{+} ion collides with the H_{2} or D_{2} molecule, generating the longlived complex. One H or D atom is ejected only when the rotational quantum number of the continuously rotating complex is large enough to cleave the HH or DD bond. Thus, the majority of the collision energies are transformed into translational energy of the product molecules, and the transformed internal energy of the products is relatively low. In addition, a part of the collision energy is used to overcome the centrifugal barrier at the exit channel. As the collision energy increases, more rotational states are populated, and the peak of the rotational state distribution shifts to a higher j′ value; this indicates that more energy is transformed into the rotational energy of the NH^{+} or ND^{+} molecule. Compared with the N^{+}(^{3}P) + H_{2} reaction, the N^{+}(^{3}P) + D_{2} reaction more easily excites the product to high rovibrational states even though it has a larger threshold. This is because the vibrational frequency and rotational constant of ND^{+} molecule are smaller than those of NH^{+} molecule, and the title reaction is very sensitive toward this subtle difference. Moreover, the N^{+}(^{3}P) + D_{2} reaction has smaller centrifugal barriers than the N^{+}(^{3}P) + H_{2} reaction at the special partial wave, which leads to the transformation of more collision energy into internal energy. Therefore, the isotope effects have considerable influence on the product rovibrational state distributions.

 Fig. 8 Rovibrational stateresolved ICSs for the (a–d) N^{+}(^{3}P) + H_{2} and (e–h) N^{+}(^{3}P) + D_{2} reactions at four collision energies (0.1, 0.3, 0.5 and 0.7 eV).  
The DCSs can describe scattering features in greater detail, whereas the angular distributions of the NH^{+} and ND^{+} molecules formed by the two reactions have not been presented to date. The threedimensional plots of the total DCSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions at collision energies below 0.7 eV are given in Fig. 9. The angular distributions of NH^{+} and ND^{+} molecules are nearly forwardbackward symmetric in the whole selected collision energy range. The peaks appear at two extreme angles, which is attributed to the longlived complexforming mechanism. The collision energy dependence of the DCSs shows sharp, dense oscillatory structures at each scattering angle; this stems from the longlived quasibound states supported by the deep well, which lead to densely distributed quantum resonances.

 Fig. 9 Threedimensional plots of the total DCSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions as a function of collision energy.  
Fig. 10 displays the total and vibrational stateresolved DCSs of the two reactions at 0.7 eV collision energy. The total DCSs show that the forward and backward peaks are almost equal, and the angle distributions are symmetric with respect to 90°. However, the symmetries of the DCSs are broken at the vibrational stateresolved levels of the product molecules for the two reactions. The forward scattering becomes more obvious when the product is excited to higher vibrational states, which suggests that both reactions demonstrate nonstatistical behavior at the vibrational stateresolved levels of the product molecules despite their complexforming natures. Fig. 11 shows several rational states (j′ = 0, 5, 10, 15) of the resolved DCSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions at 0.7 eV collision energy, and the corresponding vibrational state is v′ = 0. It is clear that the asymmetry of the forward–backward scattering is more apparent. For example, at the rotational state j′ = 10, the forward peak of the N^{+}(^{3}P) + H_{2} reaction is three times higher than the backward peak. In the case of the N^{+}(^{3}P) + D_{2} reaction, the peak at 0° disappears in the j′ = 0 state, and the forward bias becomes dominant when the product is excited to high rotational states. This asymmetric behavior is also found in some longlived complexforming reactions, such as the H^{+} + D_{2},^{27} S(^{1}D) + H_{2}^{59} and H + O_{2}^{60} reactions. The basic reason for the asymmetric angle distributions may be the strong scattering resonances even after summation over all the partial waves. Therefore, the quenching effect of the quantum interferences between various partial waves is less efficient than intuitively expected at the stateresolved level. A comparison study^{59} of quantum and statistical methods for several typical complexforming reactions has proved that statistical theories are unable to accurately study product stateresolved reaction dynamics, and the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions have similar nonstatistical behaviors. In addition, the DCSs of the two reactions have dramatic differences at specific rotation states; thus, the isotope effects have a significant impact on the stateresolved angle distributions of the product molecules.

 Fig. 10 Total and vibrational stateresolved DCSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions at 0.7 eV collision energy.  

 Fig. 11 Several rotational stateresolved (j′ = 0, 5, 10 and 15) DCSs for the N^{+}(^{3}P) + H_{2} and N^{+}(^{3}P) + D_{2} reactions at v′ = 0 and a collision energy of 0.7 eV.  
4 Conclusion
In the present work, a new PES for the ground state NH^{+}_{2} system was constructed over a large coordinate space; it can be used to obtain accurate dynamics information for the significant interstellar reaction of N^{+}(^{3}P) + H_{2}(X^{1}Σ^{+}_{g}) → NH^{+}(X^{2}Π) + H(^{2}S). The ab initio energy points are calculated using the icMRCI + Q method with high level basis sets (AVQZ for H atom and WCVQZ for N atom). A total of 30519 geometries were selected to fit the PES using the PIPNN model, and the fitting RMSE was only 2.7 meV. The structures and characteristics of the new PES are described in detail and can give more accurate information than previous PESs. TDWP calculations for the N^{+}(^{3}P_{0}) + H_{2} (v = 0, j = 0) and N^{+}(^{3}P_{0}) + D_{2} (v = 0, j = 0) reactions were carried out on the NN PES to study the reaction mechanism and dynamical isotope effects. Obvious quantum resonances are found in the reaction probabilities, ICSs, and DCS curves due to longlived complex formation in the deep well (6.49 eV). The rovibrational stateresolved ICSs reflect that the isotope effects greatly affect the state distributions of the product molecules. The total DCSs have almost forward–backward symmetric angular distributions; however, the stateresolved DCSs show nonstatistical behavior.
Herein, accurate quantum dynamics calculations for the N^{+}(^{3}P_{0}) + H_{2} (v = 0, j = 0) and N^{+}(^{3}P_{0}) + D_{2} (v = 0, j = 0) reactions were performed on the ground state adiabatic PES. A rigorous comparison of the theoretical results with experimental data is very difficult because it requires the inclusion of all accessible spin–orbit states of N^{+} ion and rotational states of H_{2} or D_{2} molecule. In addition, a set of accurate diabatic PESs correlated with the spin–orbit interaction should be used to calculate higher fine structure state reactions, which would be a very interesting but challenging task in future work.
Conflicts of interest
There are no conflicts to declare.
Acknowledgements
This work was supported by the National Natural Science Foundation of China (No. 11774043) and China Postdoctoral Science Foundation (No. 2018M631828).
References
 A. Dalgarno, E. Herbst, S. Novick and W. Klemperer, Astrophys. J., 1973, 183, L131–L133 CrossRef CAS.
 A. Dalgarno and J. H. Black, Rep. Prog. Phys., 1976, 39, 573–612 CrossRef CAS.
 N. G. Adams and D. Smith, Chem. Phys. Lett., 1985, 117, 67–70 CrossRef CAS.
 K. M. Ervin and P. B. Armentrout, J. Chem. Phys., 1987, 86, 2659–2673 CrossRef CAS.
 J. B. Marquette, C. Rebrion and B. R. Rowe, J. Chem. Phys., 1988, 89, 2041–2047 CrossRef CAS.
 M. Gonzalez, A. Aguilar and R. Sayos, Chem. Phys., 1989, 132, 137–151 CrossRef CAS.
 G. Nyman and U. Wilhelmsson, J. Chem. Phys., 1992, 96, 5198–5212 CrossRef CAS.
 P. Tosi, O. Dmitriev, D. Bassi, O. Wick and D. Gerlich, J. Chem. Phys., 1994, 100, 4300–4307 CrossRef CAS.
 L. S. Sunderlin and P. B. Armentrout, J. Chem. Phys., 1994, 100, 5639–5645 CrossRef CAS.
 I. Zymak, M. Hejduk, D. Mulin, R. Plasil, J. Glosik and D. Gerlich, Astrophys. J., 2013, 768, 86 CrossRef.
 T. P. Grozdanov and R. McCarroll, J. Phys. Chem. A, 2015, 119, 5988–5994 CrossRef CAS PubMed.
 T. P. Grozdanov, R. McCarroll and E. Roueff, Astron. Astrophys., 2016, 589, 735–745 CrossRef.
 S. G. Hansen, J. M. Farrar and B. H. Mahan, J. Chem. Phys., 1980, 73, 3750–3762 CrossRef CAS.
 B. H. Mahan and W. E. W. Ruska, J. Chem. Phys., 1976, 65, 5044–5051 CrossRef CAS.
 S. D. Peyerimhoff and R. J. Buenker, Chem. Phys., 1979, 42, 167–176 CrossRef CAS.
 D. Gerlich, J. Chem. Phys., 1989, 90, 3574–3581 CrossRef CAS.
 G. Nyman, J. Chem. Phys., 1992, 96, 3603–3612 CrossRef CAS.
 M. A. Gittins and D. M. Hirst, Chem. Phys. Lett., 1975, 35, 534–536 CrossRef CAS.
 H. P. D. Liu and G. Verhaegen, J. Chem. Phys., 1970, 53, 735–745 CrossRef CAS.
 M. Gonzalez, A. Aguilar and Y. Fernandez, Chem. Phys., 1986, 104, 57–66 CrossRef CAS.
 U. Wilhelmsson, P. E. M. Siegbahn and R. Schinke, J. Chem. Phys., 1992, 96, 8202–8211 CrossRef CAS.
 M. A. Gittins and D. M. Hirst, Chem. Phys. Lett., 1979, 65, 507–510 CrossRef CAS.
 U. Wilhelmsson and G. Nyman, J. Chem. Phys., 1992, 96, 1886–1895 CrossRef CAS.
 C. L. Russell and D. E. Manolopoulos, J. Chem. Phys., 1999, 110, 177–187 CrossRef CAS.
 L. Bonnet and J. C. Rayez, Chem. Phys. Lett., 1997, 277, 183–190 CrossRef CAS.
 L. Bonnet and J. C. Rayez, Chem. Phys. Lett., 2004, 397, 106–109 CrossRef CAS.
 P. G. Jambrina, F. J. Aoiz, N. Bulut, S. C. Smith, G. G. BalintKurti and M. Hankel, Phys. Chem. Chem. Phys., 2010, 12, 1102–1115 RSC.
 P. G. Jambrina, J. M. Alvarino, F. J. Aoiz, V. J. Herrero and V. SaezRabanos, Phys. Chem. Chem. Phys., 2010, 12, 12591–12603 RSC.
 L. Banares, F. J. Aoiz, P. Honvault, B. BusseryHonvault and J. M. Launay, J. Chem. Phys., 2003, 118, 565–568 CrossRef CAS.
 R. Tarroni, P. Palmieri, A. Mitrushenkov, P. Tosi and D. Bassi, J. Chem. Phys., 1997, 106, 10265–10272 CrossRef CAS.
 H. J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814 CrossRef CAS.
 P. J. Knowles and H. J. Werner, Chem. Phys. Lett., 1988, 145, 514–522 CrossRef CAS.
 H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schutz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253 CAS.
 H. J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053–5063 CrossRef CAS.
 P. J. Knowles and H. J. Werner, Chem. Phys. Lett., 1985, 115, 259–267 CrossRef CAS.
 R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806 CrossRef CAS.
 D. E. Woon and T. H. Dunning, J. Chem. Phys., 1995, 103, 4572–4585 CrossRef CAS.
 J. C. Yuan, D. He and M. D. Chen, Phys. Chem. Chem. Phys., 2015, 17, 11732–11739 RSC.
 J. C. Yuan, D. He, S. F. Wang, M. D. Chen and K. L. Han, Phys. Chem. Chem. Phys., 2018, 20, 6638–6647 RSC.
 Z. J. Yang, J. C. Yuan, S. F. Wang and M. D. Chen, RSC Adv., 2018, 8, 22823–22834 RSC.
 J. Cui and R. V. Krems, Phys. Rev. Lett., 2015, 115, 073202 CrossRef PubMed.
 J. Cui and R. V. Krems, J. Phys. B: At., Mol. Opt. Phys., 2016, 49, 224001 CrossRef.
 B. Kolb, P. Marshall, B. Zhao, B. Jiang and H. Guo, J. Phys. Chem. A, 2017, 121, 2552–2557 CrossRef CAS PubMed.
 Y. F. Guan, S. Yang and D. H. Zhang, Mol. Phys., 2018, 116, 823–834 CrossRef CAS.
 B. Jiang, J. Li and H. Guo, Int. Rev. Phys. Chem., 2016, 35, 479–506 Search PubMed.
 S. F. Wang, J. C. Yuan, H. X. Li and M. D. Chen, Phys. Chem. Chem. Phys., 2017, 19, 19873–19880 RSC.

K. P. Huber and G. Herzberf, Constants of Diatomic Molecules, Springer, 1979 Search PubMed.
 R. Colin and A. E. Douglas, Can. J. Phys., 1968, 46, 61 CrossRef CAS.
 B. J. Braams and J. M. Bowman, Int. Rev. Phys. Chem., 2009, 28, 577–606 Search PubMed.
 B. Jiang and H. Guo, J. Chem. Phys., 2013, 139, 054112 CrossRef PubMed.
 M. T. Hagan and M. B. Menhaj, IEEE Trans. Neural Netw. Learn. Syst., 1994, 5, 989–993 CrossRef CAS PubMed.
 D. M. Hirst, Mol. Phys., 1978, 35, 1559–1568 CrossRef CAS.
 Z. G. Sun, S. Y. Lee, H. Guo and D. H. Zhang, J. Chem. Phys., 2009, 130, 174102 CrossRef PubMed.
 Z. G. Sun, H. Guo and D. H. Zhang, J. Chem. Phys., 2010, 132, 084112 CrossRef PubMed.
 Z. G. Sun, X. Lin, S. Y. Lee and D. H. Zhang, J. Phys. Chem. A, 2009, 113, 4145–4154 CrossRef CAS PubMed.
 S. GomezCarrasco and O. Roncero, J. Chem. Phys., 2006, 125, 054102 CrossRef PubMed.
 J. C. Yuan, D. He and M. D. Chen, Sci. Rep., 2015, 5, 14594 CrossRef CAS PubMed.
 Z. P. Sun, C. F. Zhang, S. Y. Lin, Y. J. Zheng, Q. T. Meng and W. S. Bian, J. Chem. Phys., 2013, 139, 014306 CrossRef PubMed.
 P. Larregaray and L. Bonnet, J. Chem. Phys., 2015, 143, 144113 CrossRef PubMed.
 Z. G. Sun, D. H. Zhang, C. X. Xu, S. L. Zhou, D. Q. Xie, G. Lendvay, S. Y. Lee, Y. Lin and H. Guo, J. Am. Chem. Soc., 2008, 130, 14962 CrossRef CAS PubMed.

This journal is © the Owner Societies 2019 