A new potential energy surface for the H2S system and dynamics study on the S(1D) + H2(X1Σg+) reaction

We constructed a new global potential energy surface (PES) for the electronic ground state (1A′) of H2S based on 21,300 accurate ab initio energy points over a large configuration space. The ab initio energies are obtained from multireference configuration interaction calculations with a Davidson correction using basis sets of quadruple zeta quality. The neural network method is applied to fit the PES, and the root mean square error of fitting is small (1.68 meV). Time-dependent wave packet studies for the S(1D) + H2(X1Σg+) → H(2S) + SH(X2Π) reaction on the new PES are conducted to study the reaction dynamics. The calculated integral cross sections decrease with increasing collision energy and remain fairly constant within the high collision energy range. Both forward and backward scatterings can be observed as expected for a barrierless reaction with a deep well on the PES. The calculated integral cross sections and differential cross sections are in good agreement with the experimental results.

quantum-mechanical (QM) scattering method at a collision energy of 2.24 kcal/mol. They found that the DCS result exhibits forward/backward symmetry, which is characteristic of an insertion reaction. Later, Banares et al. 16 conducted QM and QCT calculations at collision energies of 2.24 kcal/mol and 3.96 kcal/ mol, and the theoretical results were used to simulate the experimental data. Based on Ho PES, the isotope effect of the title reaction was studied 18,21 using QM and QCT methods. In 2009, Song et al. reported a realistic global PES with the double many-body expansion (DMBE) method, which is widely used to construct PESs [27][28][29] , based on 1984 energy points calculated at the MRCI level using the full valence complete active space (FVCAS) reference function 12 . The PES is referred to as the DMBE/CBS PES. The intramolecular isotope effect obtained by exploratory dynamics calculations on the DMBE PES reproduces the experimental result 12 . In the same year, Song et al. provided another PES for this reactive system using a DMBE-scaled external correlation (SEC) method based on 1972 energy points 13 , with the PES is referred to as the DMBE/SEC PES. Hankel et al. 20 presented exact quantum ICSs and DCS for the title reaction using DMBE/CBS PES and DMBE/SEC PES, and they found that the results obtained on the DMBE/CBS PES were smaller due to an unphysical barrier. The influence of the barrier on the DMBE/ CBS PES at low energies was investigated 9 in detail through a comparion with the results obtained on the Ho PES, and the two dynamical behaviors were very different. Jambrina et al. 22 calculated the cumulative reaction probabilities of the title reaction using three different theoretical approaches: time-independent QM, QCT and statistical QCT. The agreement among the three methods was good indicating that the title reaction can be considered statistical.
The reaction of S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) can occur within the very low collision energy range due to exothermicity and the absence of a barrier along the reactive path. Recently, reaction dynamics in the low temperature regime have attracted a great deal of attention, because the information on reactive processes within the very low kinetic energy range can be explored by newly developed experimental techniques. In 2010, Berteloite et al. reported the absolute rate coefficients down to 5.8 K and ICSs to the collision energy as low as 0.68 meV obtained by kinetics and crossed-beam experiments, and they found that the cumulative contribution of various partial waves can explain the behavior of the excitation function 6 . In 2012, the ICSs for the two channels of the S( 1 D) + HD reaction were obtained through crossed-beam experiments down to 0.46 meV. 7 On the reactive path of the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction, there is a deep and wide well in which the long-life complex can form at a low collision energy; thus, it can be expected that the dynamics results are extremely sensitive to the PES. Therefore, the accuracy of the PES is crucial for dynamics calculations of the title reaction because significant error in reaction dynamics may arise from a small inaccuracy on the PES. Moreover, for the barrierless reaction, the long-range interaction plays an important role, especially within the very low collision energy range. Hence, a correct description of the potential energy in the long-range region is of great importance to the dynamics calculation. In this work, our major goal is to construct an accurate global PES for the H 2 S ( 1 A′ ) system, which can meet the requirements above. To achieve this objective, we calculated a mass of ab initio energy points in a large region of configuration space, and used the neural network (NN) method to fit the new PES. Then, time-dependent wave packet (TDWP) calculations were conducted for the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction to study the reactive mechanism and verify the new PES.

Results
Characteristics of the potential energy surface. Table 1 shows the equilibrium structural parameters of the new H 2 S PES. The experimental results 30 and the previous theoretical results from two analytical PESs 12,13 are also listed in the table. For both the geometry and the energy, the agreement between the calculated and experimental results is good, and our results are closer to the experimental data. The contour maps for the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction in internal coordinates at four S-H-H angles (45°, 90°, 135°, and 180°) are shown in Figure 1. In each map, there are two valleys: the left valley corresponds H( 2 S) + SH(X 2 Π ), and the valley at the bottom corresponds to S( 1 D) + H 2 (X 1 Σ g + ).
There are wells at three angles (45°, 90° and 135°), and the wells are deeper for larger angles. For the angle of 180°, there is no well but there is a barrier, indicating that the title reaction with a sufficiently low collision energy cannot occur through the S-H-H collinear path. For a better understanding of the characteristics of the PES, the minimum energy paths (MEPs) of the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH (X 2 Π ) reaction at the four S-H-H angles are shown in Fig. 2 Figure 3(a) shows an energy plot for a sulfur atom moving around a H 2 molecule with a fixed bond length at the equilibrium distance. The zero energy is set as the energy in the configuration when the sulfur atom is far from the H 2 molecule. There is a well 1.58 eV deep at x = 0.0 angstrom, y = 1.37 angstrom. For the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction, as the sulfur atom moves slowly to the H 2 molecular, it is attracted to the well. Therefore, it is expected that the insertion reaction plays an important role in the title reaction. Figure 3(b) shows an energy plot for a hydrogen atom moving around a SH molecule of which the bond length is fixed at the equilibrium distance and the zero energy is set as the energy in the configuration when the hydrogen atom is far from the SH molecule. There is a deep well (4.15 eV) close to the sulfur atom at of x = 0.74 angstrom, y = 1.33 angstrom. Figure 4 presents the reaction probabilities (J = 0, 10, 20 and 30) of

Molecular reaction dynamics.
reaction calculated by the TDWP method as a function of the collision energy. As shown in Fig. 4, a large number of sharp and dense oscillations are found in the    , which also has a deep well and a small exoergicity, has similar highly oscillatory feature. However, for the insertion reactions 32  In the TDWP calculation, the maximum total angular momentum quantum number is 45, and the threshold energy for J = 45 is above 0.2 eV. Therefore, the collision energy of 0.2 eV is set as the upper limit of the collision energy in the calculations of the ICSs. Figure 5 shows the total ICS of the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction calculated by the TDWP method. Within the very low collision energy range, the calculated ICSs have high values and decrease rapidly with increasing collision energy, which is a typical feature of exothermic reactions without a barrier. With increasing collision energy, the rate of decrease becomes lower and the values of the ICSs remain fairly constant within the high collision energy range. Oscillations are present on the curve of the calculated ICSs, which is a sign of quantum resonance. For the H 2 S system, there is a sufficiently deep potential well on the PES to support a mass of quantum states, which leads to the resonance effect. For comparison, three previous results are also shown in Fig Fig. 6. The unit in the measurement is an arbitrary unit (a.u.); thus, it is necessary to scale the experimental results to each of the theoretical results. Each scale coefficient is determined to minimize the root mean square error (RMSE) between the experimental data and the corresponding theoretical results. As shown in the figure, the agreement between our calculated results and the experimental results is strikingly good for the range of the collision energy, and the ICSs obtained from the Ho PES also match the experimental data. For the DMBE/CBS PES, the results are too flat to reproduce the experimental results. The DMBE/SEC PES results agree with the measurement results at high collision energy, but at low collision energy (the first two data points), the calculated values are significantly smaller than the experimental results.
To further verify the new PES, we calculated the total DCSs of the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction according to the TDWP method at the two collision energies (2.24 kcal/mol and 3.96 kcal/mol) displayed in Fig. 7. The experimental results 5 and several theoretical results 16,20 are also shown in the figure. Because the theoretical results obtained from DMBE/CBS and DMBE/SEC PESs are relative in the paper of Hankel and co-workers 20 , the previous theoretical results have been scaled separately to make the value at 90° equal to our result for comparison. For the experimental results and all of the theoretical results, both forward and backward scattering can be observed in the reaction as expected for a reaction with a deep well and no barrier. The DCSs are not strictly backward-forward symmetric. At the collision energy of 2.24 kcal/mol, a slight forward preference is found in the experimental results, and the DCSs obtained from the new PES reproduce this feature well. However, other previous theoretical results, particularly the DCSs 16 obtained from Ho PES, show an opposite preference. For the collision energy of 3.96 kcal/mol, the forward preference in the experiment results is more significant, and our results match the preference well. The DCSs obtained from DMBE/CBS and DMBE/SEC PESs show a significant backward preference, contrary to the experimental results.

Discussion
In this work, we report a new global PES for the electronic ground state of the H 2 S system. A mass of ab initio energies over a large configuration space are calculated at the MRCI level with the aug-cc-pVQZ basis sets plus the Davidson correction. The diatomic potential and the three-body term are fitted by the NN method. For the diatomic potential, the precision of the fitting is extremely high, and the RMSEs are 5.1 × 10 −5 eV for ( ) V HH 2 and 1.5 × 10 −5 eV for ( ) V SH 2 . The overall RMSE of the new PES is only 1.68 meV. For the equilibrium structural parameters of the H 2 S system, the agreement between the calculated and the experimental results is good, and our results are closer to the experimental data than those obtained from two previous analytical PESs 12,13 . Based on the new PES, the TDWP calculation is performed for the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction. The reaction probabilities, ICSs and DCSs of the title reaction are obtained. A large number of sharp oscillations are found in the reaction probabilities, especially at low collision energies. The deep well on the H 2 S PES and small exoergicity of the title reaction are the major causes of the highly oscillatory feature. To further verify the new PES, we have compared our dynamics results (ICSs and DCSs) with experimental 4,5 and theoretical results 16,20,21 . The ICSs obtained from our new PES agree with the experimental results 4 for the entire range of collision energy. The results 20 from two DMBE PESs do not match the experimental data exactly. The ICSs 21 from the Ho PES reproduce the experimental data well, but the DCSs 16 obtained from Ho PES at the collision energy of 2.24 kcal/mol show a significant backward preference, contrary to the experimental results. In the experimental results, a slight forward preference is found at the collision energy of 2.24 kcal/mol, and the forward preference becomes more significant at the collision energy of 3.96 kcal/mol. The DCSs on our new PES reproduce this forward preference. The close agreement of our results with the experimental data shows that the new PES is sufficiently accurate to describe the dynamics of the title reaction. In addition to the high accuracy, another advantage of the new PES is the large configuration space employed in the calculations for the ab initio energy points. The points for our new PES extend to larger distances (18.5 Å) than other previous PESs (less than 6.5 Å). Thus, the new PES provides a good description of not only short-range regions but also long-range regions. The title reaction is barrierless and exothermic; thus, the interaction within the long-range regions is important, especially in a low-temperature regime. Based on the qualities of the new PES described above, we believe that the new PES provides a good description of the reaction dynamics in the low-temperature regime, and additional research will be conducted in the future.

Methods
Potential Energy Surface. Ab initio calculations. The ab initio calculations were conducted at the MRCI level with the CASSCF reference wave functions with the MOLPRO package 33 . The aug-cc-pVQZ basis set of Dunning was employed for H and S atoms. In the arrangement channel S( 1 D) + H 2 , the supermolecules are correlated with three states (1A′ , 2A′ , 1A″ ). Therefore, all of the three states were assigned equal weight factors in the state-averaged complete active space self-consistent field (SA-CASSCF) calculation. The wave function of SA-CASSCF was used in the internally contracted MRCI method as a reference. In addition, the higher excitations were included using a multi-reference Davidson correction (+ Q). Eight valence electrons were in 10 active orbitals (9a′ + 1a″ ) composed of two 1s orbitals on H and 2s, 2p, 3s and 3p orbitals on S. Furthermore, the remaining orbitals were doubly occupied in all ab initio calculations. To achieve a high-accuracy PES, a large number of configurations (~21,300) are chosen to calculate the ab initio energy. Due to the importance of the long-range potential at the low collision energy range, a large region of configuration space is applied in the calculation: the S-H 2 region is defined by 0.37 ≤ R HH /Å ≤ 18.5, 0 ≤ R S-HH /Å ≤ 18.5, 0 ≤ θ/degree ≤ 90, and the H-SH region is defined by 0.79 ≤ R SH /Å ≤ 18.5, 0 ≤ R H-SH /Å ≤ 18.5, 0 ≤ θ/degree ≤ 180. For the diatomic potential of H 2 and SH, 141 and 146 points are calculated using the method described above, respectively.
Fitting the potential energy surface. The global analytical surface of the triatomic system can be presented as where R is a collective variable of all internuclear distances, ( ) V n 2 (n = HH, SH a , SH b ) is the diatomic potential, R n is the distance between two atoms, and ( ) V SHH 3 is the three-body term. In this work, the diatomic potential and the three-body term are fitted by the NN 34 method, which is inspired by the central nervous system of animals. The basic unit of NN is the neuron, and as a synapse, the neuron receives input signals and emits an output signal. The output signal y can be written as follows: 1, …, N) is the input signals, w i is the connection weight, b is a bias, and φ is a transfer function. The structure of the NN is important to the computational efficiency and fitting precision; thus, we used a series of tests to determine the structures. For the two-body potential of HH, two hidden layers are included in the NN, and there are six neurons in each neuron. Two hidden layers are used to fit the different parts of the SH two-body potential, and the structures of the hidden layers are (5,5) and (4,4), respectively. The precision of the fitting is extremely high and the RMSEs are 5.1 × 10 −5 eV for ( ) V HH 2 and 1.5 × 10 −5 eV for ( ) V SH 2 . For the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction, the product molecule SH(X 2 Π ) dissociates to the ground-state atoms S( 3 P) and H( 2 S). Due to the transition of the state of the sulfur atom, the energy changes dramatically at certain configurations, forming cusp structures. Most of these cusp structures are far from the strong interaction region, and the energies near cusp structures are relatively high. Therefore, the region near the cusp structures has a slight impact on the calculation of the reaction dynamics, especially at the low collision energy range. However, cusp structures cause troubles in the fitting process because the fitting errors near cusp structures are relatively large, which may reduce the accuracy of the PES. To solve this problem, Zyubin et al. 10 used an anti-cusp to compensate, whereas Song et al. introduced a switching function 12,13 . In the present work, the configuration space is divided into two parts. Most of the configurations with relatively low energy exist in one part named Part 1, and most of the cusp structures exist in another part named Part 2. The NN fitting is performed for the two parts separately. With the divided method, the fitting precision improves greatly for Part 1 and Part 2. The high fitting precision for Part 1 is important for the calculation of the reaction dynamics, especially at low collision energies because most of the interaction regions are included in Part 1. The low-order permutation invariant polynomials (PIPs) 35,36 are applied in the fitting process to solve the problem of the adaptation of the permutation symmetry. For each part, more than twenty NN PESs are obtained, and several NN PESs (seven NN PESs for Part 1 and four NN PESs for Part 2) with the least fitting errors are selected to structure an average PES. The overall RMSE of the average PES is small (only 1.68 meV).
Dynamical Calculations. Based on the new PES, the dynamics of the S( 1 D) + H 2 (X 1 Σ g + ) → H( 2 S) + SH(X 2 Π ) reaction are studied by the TDWP method. The TDWP method is a powerful tool with which to study reaction dynamics and has been applied in our previous studies [37][38][39] . In this paper, the TDWP method is introduced in brief, and more details can be found in the relevant literature 40,41 . In the body fixed representation, the reactant Jacobi coordinates are used. The Hamiltonian can be written as where R and r are defined as the distances of S-H 2 and H-H, respectively. μ R and μ r are the reduced masses associated with R and r coordinates, respectively. Ĵ is the total angular momentum operator, and ˆj is the angular momentum operator of reactant diatom molecule H 2 . V is the potential energy of the H 2 S system. The reactant coordinate based method 42 is used to extract the state-to-state S-matrix. The state-to-state reaction probability can be obtained by in which θ is the scattering angle. In this work, the ground rovibrational state (v 0 = 0, j 0 = 0) of the reactant is considered. Numerous tests were conducted to determine the optimal numerical parameters shown in Table 2.