Nonadiabatic dynamics studies of the H(²S) + RbH(X¹Σ⁺) reaction: based on new diabatic potential energy surfaces

Abstract

The global diabatic potential energy surfaces (PESs) that correspond to the ground (1²A′) and first excited states (2²A′) of the RbH₂ system PES are constructed based on 17 786 ab initio points. The neural network method is used to fit the PESs and the topographic features of the new diabatic PESs are discussed in detail. Based on the newly constructed diabatic PESs, the dynamics calculations of the H(²S) + RbH(X¹Σ⁺) → Rb(5²S) + H₂(X¹Σ_g⁺)/Rb(5²P) + H₂(X¹Σ_g⁺) reactions are performed using the time-dependent wave packet method. The dynamics properties of these two channels such as the reaction probabilities, integral cross sections, and differential cross sections (DCSs) are calculated at state-to-state level of theory. The nonadiabatic effects are discussed in detail, and the results indicate that the adiabatic results are overestimated from the dynamics values. The DCSs of these two channels are forward biased, which indicates that the abstraction mechanism plays a dominant role in the reaction.

---

1. Introduction

The Born–Oppenheimer (BO) approximation separates the motions of the electrons and nucleus of an atom and is the foundation of modern theoretical chemistry. However, the BO approximation fails at describing electronic state crossings or degeneracies (nonadiabatic effects), and these phenomena are ubiquitous in molecular reaction dynamics. Thus, to further our understanding of reaction mechanisms, accurate diabatic potential energy surfaces (PESs) are necessary. Solving for the mixing angle between different electronic states is a key step during the construction, and is often challenging. A rigorous approach to obtaining the mixing angle is to calculate the integral of the nonadiabatic coupling matrix elements (NACMEs) between different electronic states. However, some drawbacks are found when the NACME approach is applied to polyatomic systems. First of all, the NACMEs are not equal to zero at positions far from the conical intersection area. Secondly, singularities in the NACMEs and abrupt changes within certain ranges of configuration space result in discontinuities in the derivatives near the conical intersection area. Thirdly, the integration path through NACME is not unique for systems containing more than two atoms. As a result, some approximation methods have been developed to build diabatic PESs, such as the CI coefficients,^1–3 molecular properties,^4–6 and diabatization by Ansatz⁷ methods.

Reactions of X + H₂ (X = Li, Na, K, Rb, Cs) are highly endothermic, and nonadiabatic coupling occurs along the reaction path. Thus, X + H₂ reactions and their reverse reactions are ideal models for the study of nonadiabatic processes. In present work, the H + RbH reaction was chosen. Because Rb is highly endergonic, it is always excited from ground state to an excited state in experiments using a laser. For example, using a crossed beam apparatus, Cuvellier and co-worker⁸ studied the electronic-to-vibrational (E–V) energy transfer processes of the Rb(7S) + H₂(j = 1) → Rb(5D) + H₂(j = 3) reaction in 1986. The results indicated that the near-resonant process has a larger cross section than the non-resonant process at low collision velocities. In 1997, the E-V energy transfer processes between Rb(5²P_1/2,3/2) atom and H₂ molecule were investigated by Chen et al.⁹ using coherent anti-stokes Raman scattering spectroscopy. To best of our knowledge, only Fan et al.¹⁰ studied the reaction between Rb(5²D, 7²S) atom and H₂ in 1999. The relative reactivity of different excited states of Rb with H₂ followed Rb(7²S_1/2) > Rb(5²D_3/2) > Rb(5²D_5/2) and was determined by comparing the spectral intensities of RbH action spectra and Rb atomic fluorescence excitation spectra. In addition, the results also indicated that the available energy was mainly transformed into translational energy (80%).

Theoretical studies of the H + RbH reaction were mainly carried out by the group of Pascale^11,12 in 1987. The non-reaction collisions processes between Rb(5²P) and H₂, D₂ molecules were studied using close-coupling quantum mechanical methods based on a diabatic PES the authors constructed. The fine-structure transition cross sections of between Rb(5²P_1/2) and Rb(5²P_3/2) were calculated and compared with available experimental data.

As far as we know, no theoretical studies of the Rb + H₂ reaction and its reverse reactions until now. The main reason is that few PESs are available for this system. Thus, the aims of present paper are to present new global nonadiabatic PESs for the RbH₂ system and to study the nonadiabatic H + RbH (v₀ = 0, j₀ = 0) reaction process. This paper is organized as follows: the theoretical methods used in the ab initio calculations as well as the PES fitting and dynamics calculations are shown in Section 2; the properties of the diabatic PESs and the dynamical results are presented in Section 3; and the conclusions are presented in Section 4.

2. Theoretical methods

2.1 Construction of the diabatic PES

2.1.1 Ab initio calculations. The ab initio energy points used in the PES fitting were calculated using the internally contracted multi-reference configuration interaction (MRCI) method^13,14 with the aug-cc-pVQZ^15,16 basis set for the H atoms and pseudopotential basis set Def2-QZVPP¹⁷ for the Rb atom. The ab initio calculations were first performed using the Hartree–Fock calculation, and then state-averaged complete active space self-consistent field (SA-CASSCF)^18,19 calculations over three lowest electronic states (1²A′, 2²A′ and 1²A′′) with equal weight were performed. The active space of the RbH₂ system was composed of 12 active orbitals (10a′ + 2a′′). In detail, 3a′ and 1a′′ are closed and related to the 4s and 4p orbitals of Rb, and the 7a′ and 1a′′ active orbitals correspond to the 1s, 2s orbitals of H and the 5s, 5p and 4d orbitals of Rb. Subsequently, the internally contracted MRCI calculation was performed based on the reference orbits determined in the SA-CASSCF calculation. The Davidson correction was applied to the final energies to compensate for the termination errors of high-order terms. Jacobi coordinates (R_Q, R_HH, θ) were used to cover the entire configuration space, and 17 786 ab initio energy points were calculated on a non-uniform grid. In detail, R_Q is the distance between Rb atom and the midpoint of H₂ molecule, and varied between 0.4 ≤ R_Q/a₀ ≤ 25; R_HH is the bond length of H₂ molecule, and varied between 0.6 ≤ R_HH/a₀ ≤ 25; θ is the angle between R_Q and R_HH varied from 1° to 89°. The geometrical structures used in the ab initio calculations are presented in Table 1. In addition, all ab initio calculations were performed using the MOLPRO²⁰ program and using C_s symmetry.

Table 1 The geometrical structures used in the ab initio calculation

RQ (bohr)	0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 4.9, 5.3, 5.7, 6.1, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 25.0
RHH (bohr)	0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 4.9, 5.3, 5.7, 6.1, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 25.0
θ (°)	1°, 15°, 30°, 45°, 60°, 75°, 89°

---

2.1.2 Diabatization scheme. As discussed in the introduction, because it is very difficult to determine the mixing angle between the coupled electronic states from the integral of NACMEs, the molecule property method was used in the present calculations. After thorough testing, the molecule property method accurately reflected the transition characteristics between coupled states in the RbH₂ system. The diabatization scheme has been widely used in previous theoretical studies.^4–6 Herein, a brief introduction is given, and more details can be found in previous literature.^21–23 In C_s symmetry, the 1²A′ and 2²A′ states in the RbH₂ system are correlated. As the symmetry changes to C_2v, the 1²A′ and 2²A′ states change to the 1²A₁ and 1²B₂ states, respectively.

As the diabatic coupling is considered between 1²A′ and 2²A′ states, the diabatic wave function ϕ_i was obtained from the adiabatic wave function ψ_i using a unitary transformation

where γ is the mixing angle. Then the diabatic energies V_ii were obtained from adiabatic energies following

V₁₁ = cos² γE₁ + sin² γE₂ (2)

V₂₂ = cos² γE₂ + sin² γE₁ (3)

V₁₂ = (E₂ − E₁)cos γ × sin γ (4)

V₁₂ = V₂₁ (5)

where V₁₁ and V₂₂ are diabatic energies corresponding to the diabatic states, and V₁₂ and V₂₁ are the coupling potential energies between 1²A′ and 2²A′ states. In present calculation, the dipole operator transition moments are used to obtain the mixing angle. To obtain the mixing angle, the ψ₃ wave function for the 1²A′′ state was taken and was not coupled with the 1²A′ and 2²A′ states. Thus, according to eqn (1), the matrix elements 〈ψ₃|P̂|ψ₁〉 and 〈ψ₃|P̂|ψ₂〉 can be presented as

〈ψ₃|P̂|ψ₁〉 = cos γ〈ψ₃|P̂|ϕ₁〉 + sin γ〈ψ₃|P̂|ϕ₂〉 (6)

〈ψ₃|P̂|ψ₂〉 = cos γ〈ψ₃|P̂|ϕ₂〉 − sin γ〈ψ₃|P̂|ϕ₁〉 (7)

The matrix elements 〈ψ₃|P̂|ϕ₁〉 = 0 and 〈ψ₃|P̂|ϕ₂〉 = 1 at high symmetry conditions. In the present work, the approximation that 〈ψ₃|P̂|ϕ₁〉 = 0 and 〈ψ₃|P̂|ϕ₂〉 = 1 was used for all geometries, not just for the high symmetry geometries. Then, the mixing angle was obtained by the following formula

2.1.3 PES fitting. Compared with traditional fitting methods, the neural network (NN) method has several numerical advantages and is more accurate. As such, the NN method has been widely used in the construction of PESs for the F + H₂,²⁴ O⁺ + H₂,²⁵ Br + H₂ (ref. 26) reactions among others. In present work, the permutation invariant polynomial neural network (PIP-NN)^27,28 method was used to fit the PES. Herein, we only make a brief introduction to the method. For the calculations, three bond lengths in RbH₂ were transformed into P_RbHa = exp(−0.2 × r_RbHa), P_RbHb = exp(−0.2 × r_RbHb) and P_HH = exp(−0.2 × r_HH). To avoid the discontinuities in the derivatives, the permutation invariant polynomial method was introduced.

G₂ = P_RbHb × P_RbHa (10)

G₃ = P_HH (11)

To ensure a random sample, the data for {G₁, G₂, G₃} were subjected to a normalization treatment

where G_i,min and G_i,max are the minimum and maximum values of G_i, and I_i is the input data. In this work, two hidden layers were used and 15 neurons were included in each layer. The Levenberg–Marquardt method²⁹ was used to train the neural network and the expansion of neural network can be written aswhere N_tot is the number of energy points, w_j,i^l are connection weights of the ith neuron in (l − 1)th layer and the jth neuron in lth layer, and b_j^l are biases of the jth neuron in lth layer, f¹ and f² are training function and it can be presented as

Over-fitting is a well-known problem in neural network calculations. To avoid over-fitting, the input data were randomly sampling into three portions, namely the training part (90 percent), testing part (5 percent) and validation part (5 percent).

2.2 Time-dependent wave packet method

The time-dependent wave packet (TDWP) method is a powerful theoretical tool for reaction scattering calculations, and it is widely used to calculate atom-diatom collision scattering^30–32 as well as the reactions involving polyatomic system.^33,34 Herein, we only make a brief introduction to the TDWP method used in the atom-diatom reaction scattering. In body-fixed frame, the Hamiltonian of the H(²S) + RbH(X¹Σ⁺) reaction for a given total angular momentum J can be presented aswhere R is the distance between the H atom and the center of mass of RbH molecule, and r is the vibrational coordinate of RbH. Ĵ and ĵ are angular momentum operators for the RbH₂ system and RbH molecule, respectively, and μ_R and μ_r are the reduced masses. The 2 × 2 matrix is the potential energy obtained from the diabatic PES of RbH₂ system.

Before the propagation, an initial wave packet needs to be defined. Here, at t = 0, the initial wave packet can be presented as

where u_k0(R) is a Gaussian function, and ϕ_v0j0(r) is the ro-vibrational eigenfunction of RbH.

In the calculation, the split operator scheme³⁵ was used to propagate the initial wave packet. Finally, the reactant coordinates-based method^36,37 was used to extract the state-to-state information. Using the state-to-state S-matrix, the state-to-state reaction probability was obtained as

the state-to-state integral cross section (ICS) was obtained asand the state-to-state differential cross section (DCS) was obtained aswhere is the state-to-state S-matrix, is the reduced rotation matrix element, and ϑ is the scattering angle.

The rate constant of the H(²S) + RbH(X¹Σ⁺) reaction was obtained by thermally averaging the translational energy of the corresponding ICS, and the formula is presented as

where k_b is Boltzmann's constant.

3. Results and discussion

3.1 The features of the diabatic PES

In present work, the V₁₁, V₂₂, and V₁₂ terms in the diabatic potential were individually fitted using the PIP-NN method. The root-mean-square-errors (RMSE) for V₁₁, V₂₂, and V₁₂ terms are 4.95, 6.45 and 9.17 meV, respectively. Fig. 1 shows the fitting errors for V₁₁, V₂₂, and V₁₂ as a function of the diabatic energy. Because the sign of V₁₂ in the PES abruptly changes with diabatic energy, the absolute values of V₁₂ were used in fitting process. The energy derivatives using the absolute values were discontinuity and nonzero values emerged for high symmetry geometries (C_2v, D_2h); however, these effects on the final quantum dynamics results are negligible. The relatively large fitting error was mainly due to the complex morphology of V₁₂. Due to the intrude of excited states, the fitting error for the PES of V₂₂ is larger than the PES of V₁₁.

Fig. 1 The distributions of fitting errors as a function of diabatic energy for V₁₁ (a), V₂₂ (b) and V₁₂ (c).

Fig. 2 shows a comparison of the fitting values and the ab initio calculation results as a function of R_Q at several selected angles. The potential energy curves (PECs) in the diabatic representation crossed smoothly near the conical intersection area and the PECs in the adiabatic representation overlapped in the asymptotic region. These results indicated that potential coupling term V₁₂ was zero in the asymptotic region. In addition, the position of crossing point becomes small with the increase of bond angle. The fitting values were in good agreement with the ab initio data which indicated that present PESs were well fitted.

Fig. 2 Comparisons of the fitting values and the ab initio data at several selected angles ((a), (b), (c) and (d) correspond to 1°, 30°, 60°, and 89°, respectively) at a fixed H₂ molecular bond length of 3.4 bohr.

The accuracy of long-range interaction potential will affect the dynamics results for the H + RbH reaction, especially in the low collision energy region. Therefore, long-range interaction potential of H close to RbH in collinear structures is shown in Fig. 3, when the RbH bond length is fixed at equilibrium. To test the accuracy of the fitting results for the long-range interaction potential, the ab initio results are also shown in Fig. 3. As shown in the figure, the fitting results are in good agreement with the ab initio values. However, when zoom in the long-range region, the differences between fitting results and ab initio values are more obvious. The fitting error is about 1 meV for the long-range region.

Fig. 3 The long-range interaction potential as H atom approach to the RbH molecule in collinear structures.

The potential energies for an H atom moving around an RbH molecule on the V₂₂ and V₁₁ PESs are displayed in Fig. 4(a) and (b), respectively, for a fixed RbH bond length of 4.78 bohr. R_Q′ is the distance between the H atom and the center of the RbH molecule. As seen in the figure, a deep well was observed near the H atom in both Fig. 4(a) and (b). This well suggests that at low collision energies, if an H atom approaches the RbH molecule from the Rb side, it will move around the RbH molecule to form an H₂ molecule. Thus, the sideward and backward scattering signals may be more obvious at low collision energies. Meanwhile, at high collision energies, the incoming H atom will be repelled and the H atom in RbH does not form chemical bond with the incoming H atom. If the H atom approach the RbH molecule from H atom side, the incoming H atom will be trapped in the deep well and an H₂ molecule will form.

Fig. 4 Color plot of an H atom moving around an RbH molecule. The bond length of RbH is fixed at the equilibrium distance (4.78 bohr). (a) Denotes the V₂₂ PES (b) denotes the V₁₁ PES.

The potential energies of the coupling term V₁₂ are shown in Fig. 5(a) and (b) for R_Q equal to 4.4 and 5.4 bohr, respectively. As seen in the figure, the largest value for the V₁₂ is about 0.6 eV in the collinear structure. In addition, V₁₂ is almost equal to zero beyond the cross region. It is very clear that the largest values of V₁₂ are mainly focused in the region 0 ≤ θ (°) ≤ 30, 1.8 ≤ R_HH (bohr) ≤ 2.7 when R_Q = 4.4 bohr and in the region 0 ≤ θ (°) ≤ 45, 1.8 ≤ R_HH (bohr) ≤ 2.5 when R_Q = 5.4 bohr. The regular shape of V₁₂ indicates that the present diabatic PESs are a reasonable description of cross region. Fig. 6 presents the mixing angles at several selected angles for R_Q = 5.4 bohr and correspond to the data in Fig. 5(b). The mixing angle is very small in the region beyond the cross area and increases in the vicinity of the cross region where it reaches a maximum value of 45°, which corresponds to the cross point. At the conical intersection at 90°, the mixing angle is equal to zero, because the coupling term V₁₂ = 0.

Fig. 5 Color plot for the coupling term V₁₂ as a function of internal coordinates when the R_Q equal to 4.4 bohr (a) and 5.4 bohr (b), respectively.

Fig. 6 The mixing angles at different selected angles when R_Q = 5.4 bohr.

The Rb + H₂ reaction studied here is very similar to the previous studied Li + H₂ reaction as both are highly exothermic reactions between an alkali metal and hydrogen. The minimum energy paths from the RbH₂ diabatic PESs calculated here and that from the LiH₂ diabatic PESs reported by He et al.³⁸ are compared in Fig. 7. As seen in the figure, the energy of the H(²S) + RbH(X¹Σ⁺) channel was set as zero. For convenient comparison, the energy of the H(²S) + LiH(X¹Σ⁺) channel was also set as zero. For the RbH₂ system, it is very clear that the H(²S) + RbH(X¹Σ⁺) reaction is highly exothermic and the energies are about −1.65 and −3.0 eV for the Rb(5²P) + H₂(X¹Σ_g⁺) and Rb(5²S) + H₂(X¹Σ_g⁺) channels, respectively. As the zero-point energy considered, the exoergic energies were −1.435 and −2.775 eV for these two channels, respectively. The H(²S) + LiH(X¹Σ⁺) reaction is also an exothermic reaction and the exoergic energies were about −1.2 and −2.4 eV for the Li(2²P) + H₂(X¹Σ_g⁺) and Li(2²S) + H₂(X¹Σ_g⁺) channels, respectively. This comparison clearly shows that the exoergic energies of the H(²S) + RbH(X¹Σ⁺) reaction are larger than that of H(²S) + LiH(X¹Σ⁺) reaction. The Rb + H₂ and Li + H₂ reactions had other similar features including a well on the V₂₂ PES and a well and a barrier on the V₁₁ PES. Additional features of these stationary points on the RbH₂ diabatic PESs, including the geometrical structure, energy, and harmonic vibrational frequency are presented in Table 2. As seen in the table, these stationary points belong to C_2v symmetry.

Fig. 7 Comparison of the global minimum energy paths of RbH₂ diabatic PESs and that of the LiH₂ diabatic PESs obtained from He et al.³⁸ as a function of the reaction coordinate.

Table 2 The spectroscopic constants of the stationary points^a

RRbHa	RRbHb	RHH	Angle (°)	En (eV)	ω1 (cm^−1)	ω2 (cm^−1)	ω3 (cm^−1)
a The unit is bohr for bond length. The angle is [HRbH].b Relative to the energy of Rb(5^2P) + H2(X^1Σg^+) asymptote.c Relative to the energy of H + RbH(A^1Σ^+) channel.
Global minimum (V22 PES)
5.611	5.611	1.435	14.7	−0.145^b	286.47	601.94	4014.17
Local minimum (V11 PES)
5.242	5.242	8.720	112.56	−1.316^c	171.40	572.42	810.67
Transition state (V11 PES)
5.761	5.761	4.128	42.0	−0.393^c	312.77	421.25	1847.46i

---

3.2 Dynamics results and discussion

The optimal numerical parameters obtained from numerous numerical simulations of the reaction probabilities at total angular momentum J = 0 are shown in Table 3 and were used in all J > 0 nonadiabatic TDWP calculation.

Table 3 The optimal numerical parameters used in the nonadiabatic TDWP calculation (the atomic unit is used unless otherwise stated)^a

a N^totR denotes the total grids in R direction, N^intR denotes the grids in the interaction region of R direction, N^totr denotes the total grids in r direction, N^asyr denotes the grids in the interaction region of r direction.
Grids range and size	R ∈ [2.5, 22.0], N^totR = 179, N^intR = 159
	r ∈ [2.5, 22.0], N^totr = 179, N^asyr = 79
	jmin = 0 ∼ jmax = 259, Nj = 259
Initial wave packet	exp[−(R − Rc)^2/2ΔR^2]cos(k0R), Rc = 16.0
Matching plane	Rv = 16.0
K-block	K = 9
Total propagation time	45 000
Time step	10
Highest J value	90

---

The total reaction probabilities of the H(²S) + RbH(X¹Σ⁺) → Rb(5²S) + H₂(X¹Σ_g⁺)/Rb(5²P) + H₂(X¹Σ_g⁺) reaction for several selected total angular momentum of J = 0, 30, 60 and 90 are displayed in Fig. 8 over a collision energy range from 0.001 to 1.0 eV. As shown in the figure, the reaction probabilities of these two channels are similar at J = 0 and both decrease with an increase in the collision energy. This is a typical phenomenon seen in exothermic reaction without barrier. In fact, it is unusual that the reaction probabilities for these two reactions are almost equal at same collision energy. Before propagation, the initial wave packet is located in the H(²S) + RbH(X¹Σ⁺) reactant channel. After propagation, the wave packets are distributed between the Rb(5²S) + H₂(X¹Σ_g⁺) and Rb(5²P) + H₂(X¹Σ_g⁺) channels according to the adiabatic-to-diabatic transformation matrix which is composed of the mixing angle. If the wave packet is evenly distributed between the two channels, the reaction probability of Rb(5²P) + H₂(X¹Σ_g⁺) channel will be larger than that of Rb(5²S) + H₂(X¹Σ_g⁺) channel. The exoergic energy of Rb(5²P) + H₂(X¹Σ_g⁺) channel is less than that of Rb(5²S) + H₂(X¹Σ_g⁺) channel, so the velocity of H in the corresponding Rb(5²P) + H₂(X¹Σ_g⁺) channel is relatively low. Therefore, the reactivity of the Rb(5²P) + H₂(X¹Σ_g⁺) channel will be higher than that of the Rb(5²S) + H₂(X¹Σ_g⁺) channel. Here the results are inconsistent with the hypothesis. Therefore, we suppose that a large proportion of wave packet is distributed in the Rb(5²S) + H₂(X¹Σ_g⁺) channel, which is the main channel of the H(²S) + RbH(X¹Σ⁺) reaction. Moreover, the difference between the reaction probabilities of Rb(5²S) + H₂(X¹Σ_g⁺) and Rb(5²P) + H₂(X¹Σ_g⁺) channels increase with increasing J value. This shows that with the increase of J, the proportion of wave packet that goes into the Rb(5²P) + H₂(X¹Σ_g⁺) channel becomes smaller and the nonadiabatic effect becomes less obvious. This conclusion is consistent with the conclusion of studies of the H + LiH reaction.³⁹

Fig. 8 Total reaction probabilities of the two channels of the H(²S) + RbH(X¹Σ⁺) reactions at total angular momentums of J = 0, 30, 60 and 90.

To further understand the nonadiabatic effects in the reaction, we also performed the adiabatic dynamics calculations based on the adiabatic PES which obtained by diagonalizing the diabatic PESs. The results are presented in Fig. 9 for J = 0 and 10. As seen in the figure, the reaction probabilities from the adiabatic results are larger, similar to the results seen for the H + LiH reaction.³⁹ However, the sum of two channels in the nonadiabatic results is even larger than the adiabatic results, especially in the high collision energy range. This is very different from the results for the H + LiH reaction, in which the sum of two channels in the nonadiabatic values is smaller than the adiabatic values. We suppose such difference can be attributed to the different exothermic energies of these two reactions. However, it is important to note that the conclusions based on the adiabatic results are unreliable because the nonadiabatic effects were not considered for both H + RbH and H + LiH reactions.

Fig. 9 The total reaction probabilities based on the adiabatic and diabatic PESs at J = 0 and 10, respectively.

The resonance structures were observed on the reaction probability curves. Only a shallow well was found on the V₂₂ PES (seen in Fig. 7). The vibrational state-resolved reaction probabilities at J = 0 of these two channels used to determine the reaction mechanism are shown in Fig. 10(a) and (b), respectively. Comparing the results in Fig. 10(a) with those in Fig. 10(b) shows that the Rb(5²S) + H₂(X¹Σ_g⁺) channel had more vibrational states than the Rb(5²P) + H₂(X¹Σ_g⁺) channel due to the larger exoergic energy of Rb(5²S) + H₂(X¹Σ_g⁺) channel than that of the Rb(5²P) + H₂(X¹Σ_g⁺) channel. As shown in Fig. 7, no well was observed in the Rb(5²S) + H₂(X¹Σ_g⁺) channel, but some resonance structures were observed. We believe a potential well will be formed on the reaction path when the product is in the vibrational excited state, and the potential well deepens as the excited vibrational state energy level increases. Thus, the number of bound and quasi-bound states will be trapped in the well, and the resonance structures formed due to the coupling between the bound states in the well and the continuum states on the different vibrational states. We believe that these resonance structures can be classified as “Feschbach resonance”. As shown in Fig. 10(a), the amplitude of the reaction probability first increases with an increase in the number of vibrational states, then reaches a maximum value around v′ = 4, 5, and finally decreases with a further increase in the number of vibrational states. In addition, as the number of vibrational states increases, the resonance structure becomes more obvious as we speculated. For the Rb(5²P) + H₂(X¹Σ_g⁺) channel, a shallow well was observed on the minimum energy path, however, as shown in Fig. 10(b), the situation was very similar to Fig. 10(a). The “Feschbach resonance” mechanism also plays an important role in the reaction.

Fig. 10 The vibrational state-resolved reaction probabilities for the Rb(5²S) + H₂(X¹Σ_g⁺) (a) and Rb(5²P) + H₂(X¹Σ_g⁺) (b) channels at J = 0.

The total ICSs of the two channels of the H(²S) + RbH(X¹Σ⁺) reaction are plotted in Fig. 11 as a function of collision energy. The ICSs of these two channels are very similar and both decrease with an increase of collision energy. In addition, the ICSs of the Rb(5²S) + H₂(X¹Σ_g⁺) channel are larger than that of Rb(5²P) + H₂(X¹Σ_g⁺) channel over the collision energies range studied. This is a predictable result. Due to the influence of nonadiabatic effect, the gap between each J partial wave constituting the ICSs of the Rb(5²S) + H₂(X¹Σ_g⁺) channel and each J partial wave constituting the ICSs of the Rb(5²P) + H₂(X¹Σ_g⁺) channel gradually increased with an increase in the J value.

Fig. 11 Total integral cross sections of the two channels of the H(²S) + RbH(X¹Σ⁺) reaction.

The vibrational state-resolved ICSs of the H(²S) + RbH (X¹Σ⁺) reaction are displayed in the left (a) and right (b) panels of Fig. 12 and correspond to the Rb(5²S) + H₂(X¹Σ_g⁺) and Rb(5²P) + H₂(X¹Σ_g⁺) reactions, respectively. For the Rb(5²S) + H₂(X¹Σ_g⁺) channel, the vibrational states below 7 are opened due to the large exoergic energy (seen in Fig. 7 approx. 2.775 eV) and the higher vibrational exited states (7, 8 and 9) are gradually opened as the collision energy further increases. This trend reflects that more energy is available below 7 and transformed into internal energy. The amplitude of the vibrational state-resolved ICS was very similar to the vibrational state distribution of the reaction probabilities (Fig. 10). The amplitude first increases, reaches the maximum value around a number of vibrational states of 4, and then decreases as the number of vibrational states further increases. The amplitude of the Rb(5²P) + H₂(X¹Σ_g⁺) channel is very similar to the Rb(5²S) + H₂(X¹Σ_g⁺) channel, but the number of vibrational states is relatively small, which can be attributed to the relatively small exoergic energy (approx. 1.435 eV).

Fig. 12 Vibrational state-resolved integral cross sections for the two channels of the H(²S) + RbH (X¹Σ⁺) reaction (panel (a) denotes the Rb(5²P) + H₂(X¹Σ_g⁺) channel and panel (b) denotes the Rb(5²S) + H₂(X¹Σ_g⁺) channel).

The rotational state-resolved DCSs of the H(²S) + RbH(X¹Σ⁺) reaction at several selected collision energies of 0.1, 0.3, 0.5 and 0.7 eV are presented in Fig. 13 and 14 for the Rb(5²S) + H₂(X¹Σ_g⁺) and Rb(5²P) + H₂(X¹Σ_g⁺) channels, respectively. As seen in the figure, the shapes of each vibrational state-resolved ICSs are similar. The vibrational state-resolved ICSs first increase with j′, reach a maximum value, and then decrease as j′ further increases. In all cases, the rotational quantum number is the broadest in the ground state and narrows as v′ increases. At a given collision energy, the total energy is fixed. Thus, the internal energy is transformed from rotational to vibrational energy as v′ increases. In addition, the rational quantum number increases with increasing collision energy, although the increases are not very apparent. In all cases, the amplitude of each vibrational state-resolved ICS increases with increasing v′ before reaching a maximum value and then decreases as v′ further increases.

Fig. 13 The rotational state-resolved integral cross sections at several selected collision energies for the H(²S) + RbH(X¹Σ⁺) → Rb(5²S) + H₂(X¹Σ_g⁺) reaction.

Fig. 14 The rotational state-resolved integral cross sections at several selected collision energies for the H(²S) + RbH(X¹Σ⁺) → Rb(5²P) + H₂(X¹Σ_g⁺) reaction.

Fig. 15 shows the DCSs of the two channels at collision energies of 0.1, 0.3, 0.5 and 0.7 eV, respectively. As seen in the figure, the forward scattering signals are very apparent at all collision energies while the sideward and backward scattering signals are not obvious. For most scattering angles, the scattering intensity decreases as the collision energy increases which is consistent with trends in the ICSs shown in Fig. 11. In addition, the decrements in the scattering intensity are anisotropic. At high collision energies, the H₂ molecule retain a higher forward scattering signal than sideward and backward scattering signals. For example, at a scattering angle of 1 degree, the DCSs are 25.1, 22.4 and 27.7 Ų sr⁻¹ for collision energies of 0.3, 0.5 and 0.7 eV, respectively. Whereas, at a scattering angle of 90°, the DCSs are 0.77, 0.50 and 0.37 Ų sr⁻¹ for 0.3, 0.5 and 0.7 eV, respectively. The DCSs are bias toward forward scattering, which indicates that the direct abstraction mechanism plays a dominant role in the reaction.

Fig. 15 Total differential cross sections of the two channels of the H(²S) + RbH(X¹Σ⁺) reaction at several selected collision energies.

The vibrational state-resolved DCSs of the two channels of the title reaction are shown in Fig. 16(a) and (b), respectively, at a collision energy of 0.5 eV. Similar with the results in Fig. 10 and 12, more vibrational states were seen in the Rb(5²S) + H₂(X¹Σ_g⁺) channel than in Rb(5²P) + H₂(X¹Σ_g⁺) channel, due to the relatively larger exoergic energy of the Rb(5²S) + H₂(X¹Σ_g⁺) channel. In addition, the shapes of the ICS for these vibrational states are very similar to the shape of the total ICS and are both biased towards forward scattering. This result also indicates that the direct abstraction mechanism dominates the reaction.

Fig. 16 The vibrational state-resolved differential cross sections for the Rb(5²S) + H₂(X¹Σ_g⁺) (a) and Rb(5²P) + H₂(X¹Σ_g⁺) (b) channels when the collision energy equal to 0.5 eV.

The specific initial state rate constants of the two channels of the H(²S) + RbH(X¹Σ⁺) reaction are displayed in Fig. 17 over a temperature range from 200 to 1000 K. As seen in the figure, the rate constants for the initial state (v = 0, j = 0) of the two channels are not very sensitive to the temperature. The rate constants are around 1.6 × 10⁻⁹ and 9.4 × 10⁻¹⁰ cm³ s⁻¹, for the Rb(5²S) + H₂(X¹Σ_g⁺) and Rb(5²P) + H₂(X¹Σ_g⁺) channels, respectively.

Fig. 17 Specific state rate constants of the H(²S) + RbH(X¹Σ⁺) reaction over a temperature range from 200 to 1000 K.

4. Conclusion

In present work, new global diabatic PESs for the ground and first excited states of the RbH₂ were constructed by fitting 17 786 ab initio points. The SA-CASSCF and MRCI methods with aug-cc-pVQZ and Def2-QZVPP basis sets were used in the ab initio calculations. The mixing angle was obtained from the molecular properties of the dipole transition moment. The PIP-NN method was used to fit the diabatization terms, V₁₁ and V₂₂, as well as the coupling term V₁₂ individually. These three terms were well fitted and the RMSEs were 4.95, 6.45 and 9.17 meV for V₁₁, V₂₂ and V₁₂, respectively. The features of the new diabatic PESs were discussed in detail and dynamics studies of the H(²S) + RbH(X¹Σ⁺) reaction were carried out based on the new diabatic PESs. The dynamics properties, including the reaction probabilities, ICSs, DCSs, rovibrational state distributions of the product, and rate constants are reported. The forward scattering signals indicate that the abstraction mechanism plays an important role in the reaction.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work is supported by the 13th Five-Year Science and Technology Key Project of Education Department of Jilin Province (Grant No. JJKH20200482KJ).

References

1. H. J. Werner, B. Follmeg and M. H. Alexander, J. Chem. Phys., 1988, 89, 3139–3151.
2. B. Jiang, C. J. Xie and D. Q. Xie, J. Chem. Phys., 2011, 135, 164311.
3. B. Jiang and D. Q. Xie, Chin. J. Chem. Phys., 2009, 22, 601–604.
4. D. He, J. Yuan, H. Li and M. Chen, Sci. Rep., 2016, 6, 25083.
5. J. Yuan, D. He, S. Wang, M. Chen and K. Han, Phys. Chem. Chem. Phys., 2018, 20, 6638–6647.
6. Z. J. Yang, J. C. Yuan, S. F. Wang and M. D. Chen, RSC Adv., 2018, 8, 22823–22834.
7. T. Lenzen and U. Manthe, J. Chem. Phys., 2017, 147, 084105.
8. J. Cuvellier, L. Petitjean, J. M. Mestdagh, D. Paillard, P. de Pujo and J. Berlande, J. Chem. Phys., 1986, 84, 1451–1458.
9. M. L. Chen, W. C. Lin and W. T. Luh, J. Chem. Phys., 1997, 106, 5972–5978.
10. L. H. Fan, J. J. Chen, Y. Y. Lin and W. T. Luh, J. Phys. Chem. A, 1999, 103, 1300–1305.
11. J. Pascale, F. Rossi and W. E. Baylis, Europhys. Lett., 1987, 3, 183–188.
12. J. Pascale, F. Rossi and W. E. Baylis, Phys. Rev. A, 1987, 36, 4219–4235.
13. H. J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
14. P. J. Knowles and H. J. Werner, Chem. Phys. Lett., 1988, 145, 514–522.
15. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
16. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
17. T. Leininger, A. Nicklass, W. Küchle, H. Stoll, M. Dolg and A. Bergner, Chem. Phys. Lett., 1996, 255, 274–280.
18. H. J. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053–5063.
19. P. J. Knowles and H. J. Werner, Chem. Phys. Lett., 1985, 115, 259–267.
20. H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schutz, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 242–253.
21. D. W. Schwenke, S. L. Mielke, G. J. Tawa, R. S. Friedman, P. Halvick and D. G. Truhlar, Chem. Phys. Lett., 1993, 203, 565–572.
22. P. Halvick and D. G. Truhlar, J. Chem. Phys., 1992, 96, 2895–2909.
23. M. Peric, R. J. Buenker and S. D. Peyerimhoff, Mol. Phys., 1990, 71, 673–691.
24. J. Chen, Z. G. Sun and D. H. Zhang, J. Chem. Phys., 2015, 142, 024303.
25. W. T. Li, J. C. Yuan, M. L. Yuan, Y. Zhang, M. H. Yao and Z. G. Sun, Phys. Chem. Chem. Phys., 2018, 20, 1039–1050.
26. W. T. Li, D. He and Z. G. Sun, J. Chem. Phys., 2019, 151, 185102.
27. B. Jiang, J. Li and H. Guo, Int. Rev. Phys. Chem., 2016, 35, 479–506.
28. B. Jiang and H. Guo, J. Chem. Phys., 2013, 139, 054112.
29. M. T. Hagan and M. B. Menhaj, IEEE Transact. Neural Networks Learn. Syst., 1994, 5, 989–993.
30. Z. G. Sun, D. Q. Yu, W. B. Xie, J. Y. Hou, R. Dawes and H. Guo, J. Chem. Phys., 2015, 142, 174312.
31. M. L. Yuan, W. T. Li, J. C. Yuan and M. D. Chen, Int. J. Quantum Chem., 2018, 118, e25493.
32. W. B. Xie, L. Liu, Z. G. Sun, H. Guo and R. Dawes, J. Chem. Phys., 2015, 142, 064308.
33. S. Liu, X. Xu and D. H. Zhang, J. Chem. Phys., 2012, 136, 144302.
34. B. Zhao, Z. G. Sun and H. Guo, J. Chem. Phys., 2016, 144, 064104.
35. M. D. Feit, J. A. Fleck and A. Steiger, J. Comput. Phys., 1982, 47, 412–433.
36. S. Gómez-Carrasco and O. Roncero, J. Chem. Phys., 2006, 125, 054102.
37. Z. Sun, X. Lin, S.-Y. Lee and D. H. Zhang, J. Phys. Chem. A, 2009, 113, 4145–4154.
38. D. He, J. C. Yuan, H. X. Li and M. D. Chen, Sci. Rep., 2016, 6, 25083.
39. W. T. Li, J. X. Sun and D. He, Phys. Chem. Chem. Phys., 2020, 22, 17587–17596.

---