Global diabatic potential energy surfaces and quantum dynamical studies for the Li(2p) + H2(X1Σ+g) → LiH(X1Σ+) + H reaction

The global diabatic potential energy surfaces which are correlated with the ground state 1A′ and the excited state 2A′ of the Li(2p) + H2 reaction are presented in this study. The multi-reference configuration interaction method and large basis sets (aug-cc-pVQZ for H atom and cc-pwCVQZ for Li atom) were employed in the ab initio single-point energy calculations. The diabatic potential energies were generated by the diabatization scheme based on transition dipole moment operators. The neural network method was utilized to fit the matrix elements of the diabatic energy surfaces, and the root mean square errors were extremely small (3.69 meV for , 5.34 meV for and 5.06 meV for ). The topographical features of the diabatic potential energy surfaces were characterized and the surfaces were found to be sufficiently smooth for the dynamical calculation. The crossing seam of the conical intersections between the and surfaces were pinpointed. Based on this new analytical diabatic potential energy surfaces, time-dependent wave packet calculation were conducted to investigate the mechanism of the title reaction. At low collision energies, the product LiH molecule tends to forward scattering, while at high collision energies, the forward and backward scatterings exist simultaneously.

Although, many experimental studies have focused on the Li(2p) + H 2 (X 1 Σ + g ) → LiH(X 1 Σ + ) + H reaction, few theoretical studies of this reaction have been reported 22,23 . The reactants Li(2p) + H 2 (X 1 Σ + g ) and the products LiH(X 1 Σ + ) + H correlate with two different adiabatic states (2 2 A′ , 1 2 A′ ). Thus, the Born-Oppenheimer (BO) approximation is not valid in the vicinity of the conical intersection 24 . To appropriately deal with the breakdown of the BO approximation, a proper diabatization scheme should be chosen to adequately treat the coupling between the two states. In 1999, Lee et al. 22 investigated the Li(2p) + H 2 collisions. The PES used in their research is limited to three different point group symmetries: C ∞v , C 2v and C s where the molecular axis of the H 2 molecule makes an angle of π /4 with respect to the velocity vector of the lithium atom. Since this PES is not global, it cannot accurately describe all the features of the intermolecular interaction potential of the Li(2p) + H 2 reaction. In 2011, two adiabatic global PESs of LiH 2 (1 2 A′ , 2 2 A′ ) were constructed by Hsiao et al. 23 , and the QCT calculations were carried out based on their PESs. In their QCT calculations, the trajectories initially evolve on the excited state surface 2 2 A′ and hop to the ground state 1 2 A′ at the exit channel. The transition probability between the two states (1 2 A′ , 2 2 A′ ) was assumed to be unity, which is an approximate method for treating the diabatic process of the reaction. A valid method for describe the effect of the diabatic process for the molecular reaction accurately, is to construct a set of diabatic PESs for the coupled electronic states. However, for the title reaction, there are no such a set of diabatic PESs that can be employed for dynamical calculations. To meet these requirements above mentioned, a set of global diabatic PESs was constructed for the lowest two 2 A′ states of LiH 2 system using the neural network (NN) method 25 . Based on the diabatic PESs, the reaction dynamics of the title reaction were investigated using the TDWP method.

Results
PESs topographical attributes. The accuracy of the ab initio computations can be assessed by comparing the calculated spectroscopic constants with the experimental values. The spectroscopic constants of H 2 and LiH obtained from the diabatic PESs and the corresponding experimental data are listed in Table 1. The H 2 and LiH spectroscopic constants were calculated in the supermolecule approach with the other atom (Li and H, respectively) at a distance of 30 a 0 . Obviously, all of the theoretical results have a good agreement with the corresponding experimental results 26,27 . The calculated equilibrium bond lengths and harmonic vibrational frequencies of H 2 and LiH are nearly identical with the experimental values. Both of the dissociation energies of H 2 and LiH are slightly smaller than the corresponding experiment data. The small error of anharmonic constants for H 2 and LiH are 0.8 cm −1 and 3.4 cm −1 , respectively. Figure 1 shows the adiabatic and diabatic potential curves for a fixed internuclear distance of HH molecule (r HH = 2.9 a 0 ) as function of R Li-HH for various values of θ, which is the angle between R Li-HH and r HH . It is noted that the cross point shifts to the smaller R Li-HH with the increase of θ. At large θ, the adiabatic potentials strongly avoid each other in the neighborhood of the crossing point; on the contrary, the diabatic potentials cross over with each other smoothly. Furthermore, in all cases, the adiabatic and diabatic energies become identical when the molecular geometries are far from the crossing point. Figure 2 shows the corresponding mixing angles which were used to constructed diabatic potential curves shown in Fig. 1. It can be seen from this figure that in the vicinity of the crossing point, the change rate of the mixing angle becomes steeper with the increase of θ. The mixing angle reaches an asymptotic value gradually with the increase of R Li-HH . Otherwise, according to Eqs (2) and (3), the cross point appears when the mixing angle θ equals to 45° and this conclusion is in agreements with the result of Fig. 1. Figure 3 shows the coupling potentials obtained from the corresponding mixing angle and the adiabatic potentials. The value of the coupling potential decreases more rapidly with the increase of θ and the coupling potential gradually decrease to zero, while the R Li-HH increases to a large value. According to Eq. (4), it is can be concluded that the smaller value of the difference between E a 1 and E a 2 is, the smaller value of coupling potential is. This conclusion is consistent with the result obtained from Fig. 3. Figure 4 shows the adiabatic and diagonal diabatic potential energy curves for θ = 90° as a function of R LiH for various HH bond lengths r HH . It clearly shows that the cross point shifts to larger R Li-HH with the increase of HH bond length. The three-dimensional diabatic PESs and the corresponding contour plots of the differences 11 of the two diabatic states are presented in Fig. 5. From the upper panel of Fig. 5, it can be seen that for the top state, the reactant region was dominated by V d 22 and the product region was dominated by V d 11 , whereas the bottom state shows a reverse case. It can be seen that the transition of the title reaction occurs in the vicinity of interaction region. The lower panel of Fig. 5 shows the exact position of the crossing seam. The red lines in these contour plots represent the isopotential line 22 11 , i.e., the conical intersection seam of the two surfaces. The Li(2p) + H 2 → LiH(X 1 Σ + ) + H reaction processes on the V d 22 surface. The minimum energy paths of the title reaction as function of R HH − R LiH coordinate for various values of θ are shown in Fig. 6. From this plot it can be seen that there is a deep well (about 0.8 eV) along the reaction path. This well has a great impact on the reaction. The reactant zero point energy (ZPE) is 0.273 eV and the product ZPE is 0.086 eV, so that the endothermicity of the title reaction is 0.225 eV, taking into consideration the ZPEs of reactant and product.
Dynamical calculations. Figure 7 shows the initial state-specified (v = 0, j = 0) total reaction probabilities at five total angular momentum quantum number J values (0, 20, 30, 40, 50) as a function of the collision energy. For the reaction probabilities of J = 0, there is a threshold about 0.225 eV which corresponds to the endothermicity of the reaction. With the increase of the J value, the threshold increases due to the emergence of a centrifugal barrier. In the curves of the reaction probabilities, there are a mass of peaks, especially at low collision energies, which is a typical feature of quantum resonance. The feature probably corresponds to long-lived resonances associated with the well (about 0.8 eV relative to the Li(2p) + H 2 (X 1 Σ + g ) asymptote) on the V d 22 . The amplitude of the  oscillations decreases with the collision energy increasing, because of the decrease in lifetime of the intermediate complex.
In TDWP calculation, the maximum total angular momentum quantum number is 65, which is large enough to calculate the integral cross sections (ICSs) and differential cross sections (DCSs) at the collision energies under 1.0 eV. Figure 8 shows the total and vibrationally resolved ICSs of Li(2p) + H 2 (X 1 Σ + g ) → LiH(X 1 Σ + ) + H reaction calculated by the S-matrix method. There are a small amount of tiny oscillations on the ICSs curves, which is different from the intense and sharp peaks on the curves of probability. It is because the resonance structure is erased by summing over all the partial waves. For the total ICSs, the threshold is 0.225 eV which is consistent with the probability for J = 0. With the collision energy increasing, the total ICSs increases and the four vibrational excitation states of product appear in the order of v′ = 1 to v′ = 4. With the collision energy increasing, the ICS of ground vibrational state increases until reaching the maximum value at collision energy of 0.505 eV. When the collision energy increases from 0.505 eV, the ICSs for the vibrational excitation states still increases, but it decreases for the ground vibrational state. Therefore, the increase of the total ICSs comes from vibrational excitation states when the collision energy is higher than 0.505 eV. Figure 9 presents the DCSs of the  Li(2p) + H 2 (X 1 Σ + g ) → LiH(X 1 Σ + ) + H reaction at three collision energies to study the angular distribution of the product LiH. As shown in the figure, the product molecule LiH tends to be forward scattering at low collision energy (0.3 eV). With the collision energy increasing, the peak at 180° appears, which implies that the trend of the backward scattering of LiH molecule becomes more and more obvious. However, the peak at 0° is still higher than that at 180°, thus the forward scattering is dominant in the title reaction.

Discussion
A set of diabatic potential energy surfaces for the ground and first excited state of LiH 2 system was constructed using a proper diabatization scheme and high-quality ab initio energy data. The diabatization scheme is based on ) are showed in (c,d), respectively. The red lines represent the position of the diabatic crossing. transition dipole moments which could reflect the correct transition properties through the conical intersection.
The ab initio calculations were performed by internally contracted multi-reference configuration interaction (MRCI) 28,29 method. The augmented correlation-consistent polarization valence quadruple-ζ (AVQZ) 30 atomic basis was employed for H atom and the Dunning-weighted correlation consistent polarized core-valence quadruple-ζ (WCVQZ) 30 atomic basis was employed for Li atom. The NN methods are used to fit the diabatic potential energy surfaces and three high-accuracy potential surfaces were constructed with extremely small root   . The endothermicity of Li(2p) + H 2 (X 1 Σ g + ) → LiH(X 1 Σ + ) + H reaction is 0.225 eV. There is a deep potential well about 0.8 eV exists in the reaction pathway. On the new PESs, the TDWP method was used to study the reaction of Li(2p) + H 2 (X 1 Σ + g ) → LiH(X 1 Σ + ) + H. The reaction probabilities, ICSs and DCSs of the title reaction were calculated. The threshold is about 0.225 eV which corresponds to the endothermicity of the reaction. A mass of peaks are found on the curves of reaction probabilities especially at low collision energies, and this feature probably corresponds to the long-lived complex formed in the well. However, there is only a small amount of tiny oscillations on the total ICS curve, because the resonance is eased by summing over all the partial waves. The DCS results show the product LiH tends to be forward scattering at low collision energy (0.3 eV). With the collision energy increasing, the forward and backward scatterings both exist in the reaction. To our knowledge, there are no experimental results that can directly compare with theoretical results presented here. We are eagerly looking forward to a molecular beam experiment for the title reaction in the future, so that a comparison between theory and experiment can be fulfilled.

Methods
Potential Energy Surface. Diabatizaton scheme. The Li(2p) + H 2 → LiH(X 1 Σ + ) + H reaction is correlated with 1 2 A′ and 2 2 A′ of LiH 2 molecule in C s symmetry. For LiH 2 system, the 1 2 A′ state corresponds to 1 2 A 1 state in C 2v symmetry, while the 2 2 A′ state corresponds to 1 2 B 2 state in C 2v symmetry. There are many methods of diabatization for the construction of the diabatic states [31][32][33][34][35][36][37][38][39][40] . In this work, the molecule property which could reflect the character of the transition between those coupled states is used to construct the diabatic PESs. The diabatization scheme used in this work is briefly described here and more detail presentation can be found in previous literatures [33][34][35] . Considering two coupling states, the diabatic wave functions φ i d can be constructed from adiabatic wave functions ψ i a by the following unitary transformation where α is the mixing angle. The diabatic energies V ii d can be obtained in terms of adiabatic energies E ii a by  12 21 The V d 11 and V d 22 are the corresponding diabatic potential energies in the diabatic states, and the V d 12 and V d 21 are the coupling potential energies between the two states. The mixing angle α can be obtained by molecular properties [35][36][37]39,40 . In our work, the transition dipole moments are used to obtain the mixing angle. According the Eq. (1), the matrix elements ψ ψ P a a 3 1 and ψ ψ P a a 3 2 can be written as follows: where ψ a 3 is a third state which not involved in the coupling. The operator P is taken to be dipole moment P z which is parallel to the z axis. An approximate treatment is that the ψ φ P a a 3 1 is set as zero and the ψ φ P a a 3 2 was set as one not just for the high symmetry geometries (C 2v , D 2h ), but for all geometries. Then the mixing angle α can be derived from: region are about 0.2 Bohr, 2.0 Bohr and 0.5 Bohr, respectively. Moreover, the ab initio energy points above 20.0 eV relative to the energy of the Li-H-H dissociation limit are excluded from the fitting process. All of these ab initio energy points are calculated at the internally contracted MRCI level and using the complete active space self-consistent field (CASSCF) 41,42 wave function as reference. Three valence electrons are included in eleven active orbitals and three states (1 2 A′ , 2 2 A′ and 1 2 A″ ) of LiH 2 are set equal weight in state-averaged CASSCF calculations. In both the CASSCF and MRCI calculations, the WCVQZ and the AVQZ basis sets are employed for the Li and H atoms, respectively. For the two states (1 2 A′ , 2 2 A′ ) of LiH 2 , 30796 ab initio adiabatic energy points are obtained to construct the diabatic energies. The transition dipole moments mentioned above are also calculated.
All ab initio calculations in this work are performed by MOLPRO 43 package.
Fitting the diabatic PESs. The diabatic potential energies are obtained by the diabatization scheme mentioned above from the raw ab initio adiabatic energies. All of the diabatic potential energies are fitted by NN method. The back-propagation NN model is chosen in the fitting procedure. To ensure the quality of the fitting, there were two hidden layers in our NN method and 15 neurons were included in each hidden layer. There are 322 free parameters in the NN. The cross validation method was employed to detect and prevent over fitting of the training procedure. In this work, the ab initio energy points were separated into three sets, called training set, testing set and validation set. The training set was used for adapting the weights of the neural network and the validation set was used for early stopping of the training process. Moreover, the PESs have been checked by scanning the profiles for guaranteeing the non-physical behavior does not exist. A single neuron is a basic unit of the NN. Every neuron receives a set of input signals {x i } from the last layer and emits an output signal y i , which can be expressed as: where b j ( j = 1, 2, 3… .N) are the bias and w i are the connection weights between the two hidden layers. The transfer function f(x) is a critical factor for the quality of fitting the PES. The form of the transfer function in this work can be written as: x x x x In order to include the permutation symmetry induced by the two H atoms, the permutation invariant polynomials 44,45 are used in our fitting work. The final NN expansion can be presented as: Reaction Dynamics. On the new diabatic PESs, TDWP calculation was performed for the Li(2p) + H 2 (X 1 Σ + g ) → LiH(X 1 Σ + ) + H reaction. The TDWP method is a powerful tool to calculate initial state selected reactive collision and has been applied widely in many reactions 14,[46][47][48][49][50] . A brief introduction of TDWP method which used here is given in this work and the detailed discussions can be found in relevant literature 48,49 . The TDWP method used in this work is based on the reactant coordinate based (RCB) approach and the method was developed by Sun et al. 49 . The basic theory RCB is to propagate an initial wave packet in the reactant Jacobi coordinates as in an initial state-selected total reaction probability calculation to obtain a scattering wave function for the product channels directly before the wave packet is absorbed in the product region. This method could treat the diabatic effect efficiently for studying non-adiabatic reaction, such as Cl + H 2 reaction 51 . For a non-adiabatic reaction which take place over two coupled PESs, the Hamiltonian of the system is a 2 × 2 matrix. The off-diagonal elements of the matrix were employed to evaluate the transition probability between the diagonal elements. The Hamiltonian can be written as where R is the distance from the Li atom to the mass center of the H 2 molecule and r is the bond length of the H 2 molecule. μ R and μ r are the corresponding reduced masses. Ĵ and ˆj are the angular momentum operators of LiH 2 system and reactant diatom molecule. V denotes the 2 × 2 matrix representation of the potential energies. The RCB method is used to extract the state-to-state S-matrix