A globally accurate potential energy surface and quantum dynamics calculations on the Be(1S) + H2(v 0 = 0, j 0 = 0) → BeH + H reaction

The reactive collision between Be atom and H2 molecule has received great interest both experimentally and theoretically due to its significant role in hydrogen storage, astrophysics, quantum chemistry and other fields, but the corresponding dynamics calculations have not been reported. Herein, a globally accurate ground-state BeH2 PES is represented using the neural network strategy based on 12371 high-level ab initio points. On this newly constructed PES, the quantum time-dependent wave packet calculations on the Be(1S) + H2(v 0 = 0, j 0 = 0) → BeH + H reaction are performed to study the microscopic dynamics mechanisms. The calculated results indicate that this reaction follows the complex-forming mechanism near the reactive threshold, whereas a direct H-abstraction process gradually plays the dominant role when the collision energy is large enough. The newly constructed PES can be used for further dynamics calculations on the BeH2 reactive system, such as the rovibrational excitations and isotopic substitutions of the H2 molecule, and the presented dynamics data would be of importance in experimental research at a finer level.


Introduction
In recent decades, the interactions between beryllium atom and hydrogen molecules have been of great attention because of their significance in astrophysics, hydrogen storage, quantum chemistry and other fields. On the one hand, the collision product BeH 2 molecule presents the fundamental and technological interest in potential applications, such as the nuclear materials and rocket fuel technology [1,2], owing to its small mass and large hydrogen-to-metal mass ratios. Moreover, the molecular BeH 2 , with a simple electronic structure, has become an excellent candidate for testing new computational methods for quantum chemistry [3][4][5][6]. On the other hand, the further product BeH molecule in the collision process of Be + H 2 is a popular testing target for the electronic OPEN ACCESS EDITED BY structure calculations in open-shell systems [7,8]. In addition, BeH is also an important interstellar molecule, which has been identified in stars and comets [9,10].
Various experimental studies on the BeH 2 system have been implemented [11][12][13][14][15][16][17][18]. Tague and Andrews first detected the BeH 2 in molecular form by using infrared spectroscopy and the matrix isolation technique [11]. In their experiment, the pulsed laser evaporated Be atoms react with the hydrogen, and the primary product BeH and BeH 2 are largely favored compared with the other four more complex product molecules of Be 2 H, HBeHBeH, HBe(H) 2 BeH and HBeBeH.
[12] synthesized the gaseous BeH 2 molecule using an electrical discharge facility, which is verified by infrared emission spectroscopy. Their study concluded that the stable BeH 2 is a linearly symmetric molecule with the BeH bond length of 1.334 Å. The highprecision infrared emission spectra of the BeH 2 and BeD 2 molecules were measured by [14]. The antisymmetric stretching modes and some hot bands of the two molecules were studied and the spectroscopic data were accurately determined. In their later study [17], the new vibrationrotation hot bands of the BeH 2 molecule were analyzed, and an accurate value was obtained for the frequency of the bending vibrational mode.
In the theoretical aspect, numerous ab initio calculations on the BeH 2 molecule have also been performed [19][20][21][22][23][24][25][26]. Martin and Lee [19] accurately calculated the quartic force field of BeH 2 using the CCSD(T) method, and the obtained spectroscopic constants are consistent with the corresponding experimental measurements. Hrenar et al. reported the first potential energy surface (PES) of the BeH 2 molecule used for the vibrational configuration-interaction calculations by a multilevel scheme [23]. Their calculated results can reproduce the experimental values of the gas phase measurements and matrix isolation. The ground-state equilibrium structure and PES of BeH 2 were calculated utilizing the CCSD(T) method combined with the cc-pVTZ through cc-pV6Z basis sets by Koput and Peterson [24]. Furthermore, the rovibrational energy levels of BeH 2 and its isotopic variations of BeD 2 and BeHD were accurately calculated by a variational method. The newest PES of the BeH 2 system was constructed by Li and Roy [25] utilizing the three-dimensional spline interpolation over 6,864 energy points with the internally contracted multi-reference configuration interaction (icMRCI)/ aug-cc-pV5Z level. On this PES, the spectral constants of the BeH 2 and BeD 2 molecules were accurately calculated and the corresponding data of the BeHD molecule were predicted.
Although the BeH 2 system has received great attention both experimentally and theoretically, most of those studies focused on its structural and spectral properties, and the dynamics mechanisms of the Be + H 2 reaction process have not been reported up to now. In theory, the most reliable approach for obtaining the accurate dynamics information of a chemical reaction is to implement rigorous quantum scattering calculations on a globally high-precision PES [27,28]. The previous PESs of the BeH 2 system are extremely reliable and  accurate for describing the BeH 2 complex, whereas they are not  suitable for the reaction dynamics calculations since some key  regions where the reaction could reach are not included. Therefore, constructing a global and accurate BeH 2 PES is a crucial premise for studying the microscopic dynamics mechanisms of this reactive system.
Herein, a high-fidelity ground-state BeH 2 PES is represented based on a mass of high-precision ab initio energy points and the permutation invariant polynomialneural network (PIP-NN) scheme [29,30]. Moreover, the quantum dynamics calculations at the state-resolved level for the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction are carried out by the time-dependent wave packet (TDWP) method [31, 32] on this newly constructed PES. The computational details and the characteristics of the PES are given in Section 2. Section 3 displays the calculated dynamics results and the relevant discussion of the dynamics mechanisms for the title reaction and Section 4 concludes this work.
2 Ground-state BeH 2 potential energy surface

Ab initio calculations
The energy points of the BeH 2 system at the 1 1 A′ state are calculated using the icMRCI method [33,34] with the Davidson correction (+Q). The molecular orbitals are optimized by the complete active space self-consistent field (CASSCF) method [35,36] before the MRCI calculations are carried out. The CASSCF orbitals are determined by the state-averaged calculations with equal weight for the 1 1 A′, 2 1 A′, 1 1 A″ and 2 1 A″ states. The active space is composed of nine active orbitals (8a′ + 1a″). The aug-cc-pV5Z basis set [37] is used for both the two different atoms. The energies calculated for the symmetrical configuration of Be-H 2 is defined by 0.8 ≤ R HH /a 0 ≤ 8.0, 0.1 ≤ R Be-HH /a 0 ≤ 16.0, 0 ≤ θ ≤ π/2, and the configuration of H-BeH is constructed by 2.0 ≤ R BeH /a 0 ≤ 10.0, 0.1 ≤ R H-BeH /a 0 ≤ 16.0, 0 ≤ θ′ ≤ π, Here, the ab initio calculations are performed utilizing Molpro 2012 software [38].

Permutation invariant polynomialneural network fitting
The ground-state BeH 2 PES can be expressed by the summing of the two-body potentials and three-body potential: where R i (i = 1, 2, 3) are the bond length of Be-H a , H a -H b and Be-H b , respectively. A switch function f(R) is used to get a better representation in the asymptotic areas of the PES, and its form is written as: where R d and R w are the central position and the constant of switch strength, respectively. The two-body potentials are obtained by a feedforward NN structure, which consists of two hidden layers with five neurons. A total of 69 and 53 ab initio points are calculated to fit the potential energy curves (PECs) of HH and BeH, respectively, and the corresponding root mean square error (RMSE) are 0.036 and 0.385 meV. Figure 1 shows that the fitting PECs of H 2 (X 1 Σ g + ) and BeH(X 1 Σ + ) molecules can pass through the center of each ab initio point. To further demonstrate the accuracy of the two-body potentials, Table 1 displays that the spectroscopic constants of the two diatomic molecules determined on the analytical PECs are in good agreement with the corresponding experimental data [39][40][41], suggesting the presented PES are sufficiently accurate for representing the reactant and product channels when the dynamics calculations are carried out. The global ground-state BeH 2 PES is represented by the PIP-NN strategy [29,30], which can rigorously assure that the constructed PES satisfies the exchange symmetry of the two hydrogen atoms, and this scheme has been widely and successfully applied to lots of molecular systems [42-51]. First, the fundamental invariants can be expressed as: where α is a constant between 0 and 1, and here the value of α is set as 0.2. Second, the symmetrized polynomial vector G = {G i } is constructed as: Finally, G is normalized as the input of the NN model: where G i,max and G i,min are the maximum and minimum values of G i , respectively. The NN model used for constructing the global CaH 2 + PES consists of two hidden layers with 12 neurons. The hyperbolic tangent function and linear function are used as the transfer functions φ in the 1-2, 2-3 layers, and 3-4 layers, respectively.

FIGURE 1
Comparison of the ab initio data and the NN fitting results of the PECs of H 2 (X 1 Σ g + ) and BeH(X 1 Σ + ).

FIGURE. 2
Distribution of the PIP-NN fitting errors of the ground-state BeH 2 PES.
Frontiers in Physics frontiersin.org The finally analytical expansion of the final PES can be presented as: where y represents the normalized potential energy. The connecting weight w and bias b between the adjacent two layers are iteratively optimized by the Levenberg-Marquardt algorithm [52]. Here, a total of 12371 molecular configurations are picked out to take part in the PIP-NN fitting, which are randomly classified into 90% training data, 5% testing data, and 5% validation data to avoid over-fitting. The parameters of w and b of the analysis PES are determined by the training data; the testing data are used to evaluate the generalization performance of the trained PES and the training should stop immediately when the testing error starts to rise; the validation data can be used for the initial assessment and adjustment of the NN model. The distribution of the fitting errors of the ground-state BeH 2 PES is plotted in Figure 2. This figure shows that the constructed PES can keep small fitting errors in the whole energy area. The overall RMSE of the PIP-NN PES is only 1.972 meV, and the energy points with an absolutely fitting error less than 0.005 eV can reach 97.2% of all the selected configurations, implying the fitting PES is globally accurate and suitable for performing the reaction dynamics studies on the BeH 2 system.

Topographic characteristics of potential energy surface
Figures 3A,B display the contour plots of the PIP-NN PES at the D ∞h and C 2v symmetries, respectively. Excellent exchange symmetry of the PIP-PES is displayed in Figure 3A. There is a deep well with an energy minimum of −6.382 eV below the asymptotic H-Be-H at R 1 = R 3 = 2.515 a 0 , and it is also the global minimum (GM) of the ground-state BeH 2 PES, which has been demonstrated in the previous theoretical and experimental FIGURE 3 Contour plots of the ground-state BeH 2 PES at the (A) D ∞h and (B) C 2v symmetries.   25]. At the relatively low collision energy, the Be atom collides with the H 2 molecule with the remarkable elongation of the HH bond, and the BeH product is formed by the dissociation of the collinear BeH 2 molecule. For panel (B), a saddle point structure with the energy value of −2.156 eV is presented at R 2 = 2.269 a 0 , R Be-HH = 2.561 a 0 , which is corresponding to the transition state (TS) of the BeH 2 system and dominates the collision process of the H-exchange path of H a + BeH b → H b + BeH a . The valley at R 2 = 1.401 a 0 corresponds to the Be( 1 S) + H 2 channel, and the GM is also shown at R 2 = 5.030 a 0 , R Be-HH = 0 a 0 since the D ∞h configuration is a limitation of the C 2v symmetry. Table 2 lists the structures, energy values and vibrational frequencies of the GM and TS for the ground-state BeH 2 calculated at the PIP-NN PES, and the available experimental and ab initio values are also presented. The energy values are relative to the Be( 1 S) + H 2 asymptotic channel. The newly constructed PES can accurately reproduce the geometries and the corresponding energies of the two stationary points, and the vibrational frequency v 2 is consistent with the experimental [14, 17] and extremely high-precision ab initio data [25] well. There   Figure 4. It is clear that the constructed PES is smooth in the entire configuration space, and there is no nonphysical structure for each angle, suggesting the over-fitting behavior does not exist during the fitting PES. For the PES at every angle, the bottom valley is the Be( 1 S) + H 2 channel, and the left valley is corresponding to the BeH + H channel. The energy of the bottom channel is lower than the left channel, indicating that the Be( 1 S) + H 2 → BeH + H reaction is endothermic. For the angles of 45°and 90°, the reactant and product channels are connected by a barrier structure, which is generated by the avoid crossing behavior of the 2 1 A′ state. The energy value of the barrier is higher than the energy of the product channel, implying the larger collision energy is needed to initiate this reaction by the collision approach with a relatively small Be-H-H approaching angle. In addition, a potential well with the depth of 2.496 eV relative to the BeH + H asymptotic region is shown when the approaching angle is at 45°, and many bound states or quasibound states can be supported by this well. For the angles of 135°a nd 180°, no well or barrier exists in the PES, thus the title reaction proceeds via a direct H-abstraction process when the collision angle becomes larger enough. Figure 5 shows the minimum energy paths (MEPs) of the Be( 1 S) + H 2 → BeH + H reaction at four Be-H-H approaching angles (45°, 90°, 135°, 180°), calculated by scanning the groundstate BeH 2 PES with the fixed angle shown in Figure 4 at different coordinates to obtain the energy minimum. In addition, the global MEP generated by scanning the whole PES is also given in this plot, which plays the dominant role in  Figure 6A displays the contour plot of the ground-state BeH 2 PES in the case of the Be atom moving around the H 2 molecule fixed at its equilibrium distance. It is clear that the Be atom is always repelled by the H 2 molecule, so initiating the title reaction is difficult when the HH bond is stabilized at its equilibrium structure. As shown in Figure 3A, when the HH bond is elongated 5.030 a 0 , there exist the attractive interactions between the Be atom and the H 2 molecule, and a stable BeH 2 complex is formed. A similar map to Figure 6A but for a H atom moving around the BeH fixed at its equilibrium distance is displayed in Figure 6B. Different from the case of Figure 6A, it appears the attractive interactions between the H atom and BeH molecule, and the well around the Be atom is deeper than the well around another H atom, suggesting that this H atom prefers to get out from the side of H atom of BeH in the product region.

Quantum dynamics calculations
For most of the triatomic and some tetratomic reactive systems, the quantum TDWP method [31, 32, 53-55] is a high-efficiency and accurate tool for calculating the dynamics data. The full-dimensional quantum dynamics calculations of the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction are carried out on this newly constructed PIP-NN PES by the TDWP method for understanding the state-resolved dynamics mechanisms. The Coriolis coupling effect is included in the quantum TDWP calculations. Here, only the main equations in the TDWP calculations are displayed below. The Hamiltonian of the title reaction can be expressed as: where µ r and µ R are the reduced masses associated with r and R in the Jacobi coordinate, respectively. J and j express the total Be (  Total propagation time 20,000 a.u.
Frontiers in Physics frontiersin.org angular momentum quantum number of BeH 2 and rotational angular momentum quantum number of H 2 , respectively. The initial wave packet consists of a Gaussian type wave function, a rovibrational eigenfunction of H 2 , and an eigenfunction of the total angular momentum, written as: To avoid the reflection of wave packet at the grid edge, the absorption potential used in the TDWP calculations is defined as: where C i and x i (i = a, b) represent the strength and positions of the absorption potential, respectively. Here, the time evolution of the wave packet is realized by the split operator scheme [56] and using the reactant coordinate-based method [57, 58] to extract the state-resolved S-matrix. The rovibrationally stateresolved reaction probability obtained by the S-matrix is expressed as: The state-resolved integral cross sections (ICSs) are calculated by summing the probabilities of all the calculated partial wave J: where k v0j0 is the momenta in the entrance channel. The stateresolved differential cross sections (DCSs) can be obtained by the following equation: whered J KK0 (ϑ) represents the reduced Wigner matrix, andϑexpress the scattering angle.
In this work, the reactant H 2 molecule is set at its groundrovibrational state of v 0 = 0, j 0 = 0, and the number of partial waves is calculated up to 65, which can obtain the convergent ICS and DCS up to the collision energy of 5.0 eV. In Table 3, the main parameters determined by many times tests of convergence in the TDWP calculations are listed.
The collision energy dependence of total reaction probabilities for the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction with four partial waves (J = 0, 20, 40 and 50) are presented in Figure 7A. For J = 0, the curve exhibits relatively dense oscillation structures, which are attributed to the potential well on the reactive path. The title reaction is dominated by the global MEP and there is a well with the depth of 1.632 eV, resulting in obvious quantum resonances because numerous bound and quasi-bound states can be formed in the well. As the increase of J values, the reactive threshold becomes larger and the oscillations are gradually weakened. This is because the increasing centrifugal barrier reduces and even smooths the effective potential well, and the other collision channels shown in Figure 5 are opened, causing the amplitudes of oscillations on the reaction probability curves become less pronounced. Figure 7B shows the collision energy dependence of total ICS for the title reaction. The total ICS value increases monotonically with the increase of collision energy, which is consistent with the characteristic of an endothermic reaction. Compared to the reaction probabilities, there is no oscillation structures on the ICS curve due to the superposition of all the calculated partial waves.
To understand the dynamics mechanisms of the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction at the state-to-state level, the rovibrationally state-resolved ICSs of the product BeH molecule at four collision energies (3.0, 4.0, 4.5, and 5.0 eV) are shown in Figure 8. For the collision energy of 3.0 eV, the BeH molecule only can be excited to the lowest three vibrational states, but the maximum of the rotational quantum number can reach j′ = 21 at v′ = 0, and the peak value of the rovibrationally state-resolved ICS is located at v′ = 0, j′ = 16. The presented vibrationally cold and rotationally hot distribution conforms to the complexforming mechanism. More rovibrational states become available with the increase of collision energy, and there is a population inversion of the vibrational quantum number. For the collision energy of 5.0 eV, the product BeH molecule can populate at very high rovibrational states (v′ = 10, j′ = 28), suggesting more collision energy is effectively transformed into the internal energy of the product molecule. The contributions of high-order partial waves are larger and more reaction paths are gradually opened as the collision energy increases, thus the lifetime of the forming BeH 2 complex becomes shorter and the title reaction prefers a direct H-abstraction process when the collision energy is large enough.
To study the dynamics process of the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction more intuitively by giving the angular distribution of the product molecule, Figure 9 presents the total DCSs varying with the scattering angle and collision energy. It is clear that the peak values of the angular distribution are located at the two extreme angles (0°and 180°) and the forward-backward symmetric DCSs are displayed when the collision energy is slightly larger than the reactive threshold, which is due to the forming of a BeH 2 complex supported by the potential well on the global MEP. With the increase of collision energy, the product BeH molecule increasingly prefers the forward scattering, showing an obviously non-statistical behavior. It also can be explained by the increasing contributions of the centrifugal barriers and more open reactive paths without a well at large collision energy. The calculated results of the total DCS further imply that the title reaction follows the complex-forming mechanism near the reactive threshold, whereas s direct H-abstraction process gradually plays a dominant role at high collision energy.

FIGURE 9
Total DCS of the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction as a function of scattering angle and collision energy calculated by the TDWP method on the ground-state BeH 2 PES.
Frontiers in Physics frontiersin.org

Conclusion
In this paper, a globally accurate ground-state BeH 2 PES is structured using the PIP-NN scheme based on 12371 ab initio points calculated at the icMRCI + Q/aug-cc-pCV5Z level. The PES can accurately reproduce the original ab initio data in each region, and the global fitting RMSE is only 1.972 meV. The molecular constants of H 2 (X 1 Σ g + ) and BeH(X 1 Σ + ) calculated on the PES are consistent with the corresponding experimental data, and the PES can reproduce the characteristics of stationary points well. The GM and TS of the ground-state BeH 2 correspond to the D ∞h and C 2v symmetries, respectively. The topographic features of the PES are described in detail. On this newly constructed PES, the dynamics calculations are performed on the Be( 1 S) + H 2 (v 0 = 0, j 0 = 0) → BeH + H reaction at the state-to-state level by the quantum TDWP method for understanding the microscopic mechanisms. The endothermicity of the title reaction determined by the PES is 2.716 eV. There exist obvious oscillation structures on the curves of reaction probabilities since the well on the global MEP can support numerous bound and quasi-bound states, and the total ICS increases monotonically with the increase of collision energy. The rovibrationally state-resolved ICSs present vibrationally cold and rotationally hot distribution at relatively low collision energy, and the product BeH molecule can populate at very high rovibrational states. The total DCSs are forward-backward symmetric when the collision energy is slightly larger than the reactive threshold, but only the forward scatting is presented at high collision energy. The dynamics results indicate that the title reaction follows the complex-forming mechanism near the reactive threshold, whereas a direct H-abstraction process gradually plays the dominant role at high collision energy. Further dynamics studies for this reaction system can be carried out on the presented PES, such as the effects of rovibrational excitations and isotopic substitutions of the H 2 molecule, and the dynamics data calculated in this paper would be of importance in the experimental studies on the title reaction.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.