Due to their importance in planetary ionospheres and interstellar chemistry, ion– molecule reactions have been the subject of a vast amount of studies.^{[1, 2]} From the experimental point of view, since it is quite simple to accelerate ions, molecular beam experiments can be carried out relatively easier on such a kind of reaction over a wide range of relative translational energies. However, it is usually difficult to explore the reaction dynamics studies theoretically on such a kind of reaction system, since the energies of ground and excited potential energy surfaces (PESs) are frequently degenerate or close, and the electronically nonadiabatic processes easily occur.^{[3– 5]} Among a series of ion– molecule reactions, the reaction is an exceptional case, which mainly takes place adiabatically on the lowest energy quadruplet (1⁴A″ in C_s symmetry) PES over a quite wide range of collision energy.^[3] Thus, the reaction which is regarded as a prototype of moderately exothermic ion– molecule reactions involving an H atom transfer, has received plenty of studies both experimentally and theoretically.

Experimentally, the rate constants for reactions in which a series of positive ions (O⁺ , N⁺ , Ar⁺ , , CO⁺ , and N₂O⁺ ) abstract a hydrogen atom from H₂ were measured by Fehsenfeld et al.^[6] at 300  K in a pulsed, flowing-after glow reaction tube. By using the ion cyclotron resonance and selected ion flow tube techniques separately, the rate constants for the reaction O⁺ + H₂ → OH⁺ + H were obtained by Kim and Smith et al., ^[7] respectively. Both the rate constants and branching fractions were reported by Viggiano et al.^[8] for the reaction of O⁺ with HD as a function of average center-of-mass kinetic energy at three temperatures: 93, 300, and 509  K. For the state-selected reactions of O⁺ (⁴S, ²D, ²P) + H₂/D₂, the absolute total cross sections were measured over the collision energy range of 0.02  eV– 12  eV.^[9] The total cross sections of O⁺ (⁴S, ²D, ²P)+ HD reaction were measured by Sunderlin and Armentrout.^[10] The velocity vector distributions for the reactive and nonreactive scattering of O⁺ by H₂, HD, and D₂ were reported by Gillen, ^{[11, 12]} which showed that such a reaction and its isotopic variations proceed from ground-state reactants to predominantly electronic ground state products via a direct interaction. An interesting conclusion can be drawn from the above experimental studies, which is that the electronically nonadiabatic processes, such as charge transfer and dissociative charge transfer, are important only at high collision energies.

Amongst the vast quantity of theoretical work, employing the quasiclassical trajectory (QCT) approach, Gonzá lez et al.^{[13, 14]} investigated the ion– molecule reaction of O⁺ (⁴S) with and its isotopic variants (D₂ and HD), based on an empirical surface derived from experimental and ab initio information. Using the diatomics-in-molecules (DIM) model, Whitton et al.^[15] constructed the PES of H₂O⁺ (X⁴A″ ). A more accurate global PES of H₂O⁺ (X⁴A″ ) was established by Martí nez et al.^[4] using many-body expansion (MBE). The properties of the most relevant stationary points of this reaction were analyzed in detail. Subsequently, in order to investigate the dynamics of the O⁺ + H₂ reaction and its isotopic variants D₂ and HD, they performed the dynamics studies using methods of QCT, ^[16] time-dependent real wave packet quantum dynamics^[17] and time-independent close-coupling hyperspherical (CCH) quantum dynamics^[18] based on such PES. Moreover, a lot of other dynamics studies were also performed using this global PES. By using centrifugal sudden (CS) approximation, time-dependent wave packet calculations have been performed by Kł os et al.^[19] for the O⁺ + H₂, and exact total reflection probabilities at the total angular momentum J = 0 and approximate ones for J > 0 have been calculated. Zhang et al.^[20] studied the influence of coriolis coupling (CC) on the O⁺ + H₂, showing CC plays an important role in this reaction. The dynamics of the exchange reaction H + OH⁺ were also investigated by Zhang et al.^[21] using the time-dependent wave packet quantum method under CS approximation. The results of dynamics obtained by theoretical calculations can be compared well with the experimental ones.

However, the PES constructed by Martí nez et al.^[4] was obtained by fitting around 600 ab initio points calculated at CCSD(T) level with cc-pVQZ basis set. We are also aware that the dissociation energy of OH⁺ (X³∑ ⁻ ) was calculated from such PES to be 5.1268  eV, while the experimental value is 5.301  eV, showing a relatively larger difference of 0.1742  eV (∼ 4.0172  kcal· mol^{− 1}). Compared with the most recent result (5.2816  eV) obtained by Paniagua, ^[5] the difference is calculated to be 0.1544  eV (∼ 3.5606  kcal· mol^{− 1}). Thus, we are motivated to constructing more accurate global analytical PES of H₂O⁺ (X⁴A″ ) by fitting more accurate ab initio energies to be calculated with multireference configuration interaction (MRCI)^[22] level and larger basis set aug-cc-pVQZ^{[23, 24]} (denoted as AVQZ). In this work, a realistic global analytical PES for H₂O⁺ (1⁴A″ ) is reported, which is obtained by fitting 2690 ab initio points that were calculated at MRCI^[22] level of theory, using the full valence complete active space (FVCAS)^[25] reference wave function and the Dunning' s basis set AVQZ. The resulting ab initio energies have been subsequently corrected, using the double many-body expansion-scaled external correlation (DMBE-SEC)^[26] method to extrapolate to the one-electron complete basis set (CBS) limit.^[27] The PES so obtained shows the correct behavior at all dissociation channels while providing a realistic representation of the surface features at all interatomic separations.

This paper is organized as follows. Section  2 presents the ab initio calculations employed in the present work. Section  3 describes the analytical representation of the H₂O⁺ (X⁴A″ ) PES. Main topographical features of the current PES are discussed in Section  4. Finally, Section  5 summarizes the concluding remarks.

In order to map the H₂O⁺ (X⁴A″ ) PES, a grid of raw 2690 ab initio energies has been calculated, which is defined by 1.2 ≤ R_H2/a₀ ≤ 3.8, 1.0 ≤ r_{O⁺ − H2}/a₀ ≤ 15.0, and 0.0 ≤ γ /deg ≤ 90 for O⁺ − H₂ channel; while 1.5 ≤ R_OH⁺ /a₀ ≤ 4.0, 2.0 ≤ r_{H− OH⁺} /a₀ ≤ 15.0, and 0.0 ≤ γ /deg ≤ 180 for H− OH⁺ channel. As described in Ref.  [28]–   [30], R, r, and γ are the atom– diatom Jacobi coordinates. By using the FVCAS^[25] as reference, all ab initio calculations have been carried out at the MRCI^[22] level. The AVQZ basis set of Dunning has been employed, with the calculations being performed using the Molpro^[31] package. C_s point group symmetry is employed in the ab initio calculation, which holds two irreducible representations, namely, A′ andA″ . For H₂O⁺ (X⁴A″ ), five A′ and oneA″ symmetry molecular orbitals (MOs) are determined as the active space, amounting to a total of 172 (110A′ + 62A″ ) configuration state functions. The molecular orbitals, which correlate with 2s, 2p atomic orbitals of the O⁺ ion and the 1s orbital of each H atom, are utilized to determine the valence space. Such a procedure provides a good description of bond formation/breaking which takes place in the reactive process, in which case the CCSD(T) usually fails.

Aiming at accounting for the excitations beyond singles and doubles and, most importantly, for the incompleteness of the basis set, the resulting ab initio energies have subsequently been modeled by scaling the external correlation employing the DMBE-SEC method developed by Varandas.^[26] According to the DMBE-SEC scheme, the total energy can be written as the following form

where

where the first terms of Eqs.  (2) and (3) sum over all the diatomic components with i = AB, BC, AC, and R = {R_AB, R_AC, R_BC} in and is a collective variable of all internuclear distances. and in Eq.  (3) are then expressed as

As is done in the previous work, ^{[32– 34]} the scaling factors of diatoms in Eq.  (4) are chosen to reproduce the corresponding bond dissociation energies of the i-th diatom, while the scaling factor F⁽³⁾ in Eq.  (5) is fixed at the average of three two-body scaling factors. For the AVQZ basis set, such a procedure leads to the two-body and three-body scaling factors as: and F⁽³⁾= 0.8602.

Employing the many-body expansion method proposed by Murrell and Carter, ^[35] the analytical functional form of H₂O⁺ (X⁴A″ ) PES can be expressed as

where is the energy for isolated atom, V⁽²⁾(R_i) represents the two-body term, and V⁽³⁾ is the three-body term which becomes zero at all dissociation limits. In the present case, the H₂O⁺ (X⁴A″ ) PES follows the dissociation scheme as

Thus, the one-body term equals zero since all possible dissociation channels of H₂O⁺ (X⁴A″ ) correlate with the atoms in their ground states (i.e., H(²S) and O⁺ (⁴S)). The two-body terms involve the and OH⁺ (X³∑ ⁻ ), with their analytical expression being represented by the formalism developed by Aguado and Paniagua, ^{[36, 37]} which was recently applied to construct the X²∑ ⁺ and A²Π states of CN, ^[38] can be written as the sum of the short- and long-range potentials,

where

which warrants the diatomic potentials tend to infinite value when R_AB→ 0. The long-range potential which tends to zero as R_AB→ ∞ takes the following expression .

The linear parameters a_i (i = 1, 2, … , n) and the nonlinear parameters β _i (i = 1, 2) in Eqs.  (9) and (10) are obtained by fitting the ab initio potential energies of the diatoms. All the coefficients used to model the diatomic potentials for OH⁺ (X³∑ ⁻ ) and are gathered in Table  1. While, the properties of both the OH⁺ and H₂ potential energy curves (PECs), including equilibrium geometries, dissociation energies, vibrational frequencies, and spectroscopic constants are gathered in Table  2, which agree excellently with other theoretical and experimental values. Shown in Fig.  1 are the OH⁺ (X³∑ ⁻ ) and PECs, which mimic the ab initio potential energies nicely, showing accurate and smooth behavior both in short and long regions.

Table 1.

Table 1.

Table 1. Fitted parameters of two-body energy term in Eq.  (10). OH+ (X3∑ − ) a0 9.9132747650 × 100 0.7905374169 × 100 a1 − 5.2367539406 × 100 − 4.4604940414 × 100 a2 1.2129990387 × 102 5.3346576691 × 101 a3 − 2.1271813965 × 103 − 7.3640325928 × 102 a4 2.2541205078 × 104 6.9652080078 × 103 a5 − 1.4416690625 × 105 − 4.2423710938 × 104 a6 5.4555181250 × 105 1.5623592187 × 105 a7 − 1.1230468750 × 106 − 3.1305646875 × 105 a8 9.7031368750 × 105 2.5589892187 × 105 β 1 2.649242 1.693474 β 2 1.333239 1.047213

Table 1. Fitted parameters of two-body energy term in Eq.  (10).

Table 2.

Table 2.

Table 2. Spectroscopic constants of OH+ and H2 diatoms, with the unit of Re in a.u. (atomic unit), De in unit kcal/mol, and ω e, ω exe, α e, and β e in unit cm− 1. Re De ω e ω exe α e β e OH+ (X3∑ − ) This work 1.9430 5.2972 3058.71 54.745 0.446 16.819 Theor.[4] 1.9428 5.1268 3126.60 – – – Theor.[5] 1.9445 5.2818 3127 – – – Theor.[41] 1.9460 5.2020 3113.37 78.515 0.750 16.794 Exp.[42] 1.9443 5.301 3113.37a – – – This work 1.4025 4.7479 4395.22 126.119 2.221 60.735 Theor.[4] 1.4010 4.7321 4351.24 – – – Theor.[43] 1.4010 4.7477 4403.60 126.602 2.232 60.864 Theor.[44] 1.4013 4.7105 4389.66 121.560 3.162 60.826 Exp.[42] 1.4014 4.476 4395.20 117.99 2.993 60.809 Exp.[45] 1.4011 4.4781b 4401.21 121.33 3.062 60.853 aRef.  [46], bRef.  [47].

Table 2. Spectroscopic constants of OH+ and H2 diatoms, with the unit of Re in a.u. (atomic unit), De in unit kcal/mol, and ω e, ω exe, α e, and β e in unit cm− 1.

By subtracting the two-body energies from the total ab initio energies, for a given triatomic geometry, one can obtain the three-body energies, which are then fitted to the three-body term V⁽³⁾ written as the M-th order polynomial^{[36, 37, 39, 40]}

where . C_jkl and are the linear and nonlinear parameters to be determined in the fitting procedure. The two expressions of Eq.  (11) are used to satisfy the permutation symmetry of the two H atoms. The constraints i + j + k ≠ i ≠ j ≠ k and i + j + k ≤ M are used to warrant that the three-body term becomes zero at all dissociation limits and when at least one of the internuclear distances is zero. M equals 10 in the present work, which results in a complete set of parameters to be determined amounting to a

Fig.  1.

	Figure Option ViewDownloadNew Window
	Fig.  1. Potential energy curves of OH+ (X3∑ − ) and The circles indicate the AVQZ/DMBE-SEC energies.

Fig.  2.

	Figure Option ViewDownloadNew Window
	Fig.  2. Deviations between ab initio energies and the fitted PES of H2O+ (X4A″ ).

6.99% of the total ab initio energies, which implies that the fitted PES is in good quality.

Table  3 compares the attributes (geometries, energies, and vibrational frequencies) of the relevant stationary points of the current PES with other theoretical results for H₂O⁺ (X⁴A″ ): C_{∞ v}[OHH]⁺ , C_2v[OH₂]⁺ , and C_{∞ v}[HOH]⁺ minima, and D_{∞ v}[HOH]⁺ transition state. Also showing in this table are the results obtained by optimizing at MRCI/AVQZ level of theory in the present work, and at CCSD(T)/cc-pVQZ level of theory by Martí nez et al.^[4] As is shown in Table  3, the C_{∞ v}[OHH]⁺ minimum located at R₁ = 2.2950a₀ (R_H2), R₂ = 2.1209a₀, and R₃= 4.4159a₀ (R_OH⁺ ). Comparing with the result calculated from analytical PES of Martí nez et al., ^[4] the derivations are only 0.0573a₀, 0.0336a₀, and 0.0237a₀, respectively. The well depth calculated from the present PES is 0.2080 E_h relative to the O⁺ (⁴S) + H(²S) + H(²S) asymptote, which is 0.006 E_h lower than the results of Martí nez et al.,   ^[4] since the present PES is obtained by fitting ab initio energies calculated with higher accuracy, namely MRCI/AVQZ. Comparing with the well depth of the other minima and transition state, the C_{∞ v}[OHH]⁺ minimum has the lowest energy, which is attributed to the global minimum of H₂O⁺ (X⁴A″ ). Both the well depth and geometry are in good accordance with the results obtained in the present work through optimization at MRCI/AVQZ level of theory, with the well depth differing by 0.0008 E_h and the maximal difference 0.0789a₀ of geometries occurs for OH bond length. The harmonic vibrational frequencies for the C_{∞ v}[OHH]⁺ minimum obtained by the present PES are 1746.27  cm^{− 1}, 863.03  cm^{− 1}, and 562.55  cm^{− 1}, which are also in good agreement with the results through optimization searches at the MRCI(Q)/AVQZ level (1747.92  cm^{− 1}, 796.11  cm^{− 1}, and 719.09  cm^{− 1}), and the results of analytical PES obtained by Martí nez et al.^[4] (1644.09  cm^{− 1}, 785.04  cm^{− 1}, and 765.04  cm^{− 1}).

Table 3.

Table 3.

Table 3. Properties of stationary points on the fitted H2O+ (X4A″ ) PES a. Feature R1/a0 R2/a0 R3/a0 E/Eha ω sym ω asym ω bending C∞ v[OHH]+ minimum This workc 2.2950 2.1209 4.4159 – 0.2080 1746.27 863.03 562.55 MRCI/AVQZd 2.3011 2.1326 4.4337 – 0.2072 1747.92 796.11 719.09 Ab initioe 2.2396 2.1545 4.4391 – 0.2029 1643.42 767.59 763.00 Analytical fitf 2.2377 2.1545 4.3922 – 0.2029 1644.09 785.04 765.04 C2v[OH2]+ minimum This workc 1.4945 3.4241 3.4241 – 0.1954 823.19 326.29 3652.18 MRCI/AVQZd 1.4929 3.5563 3.5563 – 0.1970 743.66 556.33 3665.22 Ab initioe 1.5101 3.2585 3.2585 – 0.1929 744.11 108.33 3582.54 Analytical fitf 1.5101 3.5890 3.5890 – 0.1929 736.02 192.66 3401.79 C∞ v[HOH]+ minimum This workc 8.3805 6.4362 1.9443 – 0.1944 90.69 3050.63 99.00 MRCI/AVQZd 7.6157 5.6733 1.9424 – 0.1933 199.61 3127.53 65.39 Ab initioe 7.6108 5.6679 1.9429 – 0.2029 171.71 3131.67 44.00 Analytical fitf 8.4102 6.4579 1.9523 – 0.1913 286.26 3074.26 97.70 D∞ v[HOH]+ transition state This workc 5.4107 2.7053 2.7053 – 0.1230 1519.61 2507.02i 456.82 MRCI/AVQZd 5.3144 2.6572 2.6572 – 0.1265 1438.66 3618.64i 841.35 Ab initioe 5.2994 2.6497 2.6497 – 0.1215 1416.07 2759.75i 860.59 Analytical fitf 5.4922 2.7461 2.7461 – 0.1224 1204.39 1759.44i 675.72 aHarmonic frequencies in unit cm− 1.bEnergy relative to the three-atom limit O+ (4S) + H(2S) + H(2S).c This work, fitted by the current analytical PES.d This work, obtained through optimization at the MRCI/AVQZ level of theory.eAb initio results calculated at CCSD(T)/cc-pVQZ level of theory.[4]f From the analytical PES of Martí nez et al.[4]

Table 3. Properties of stationary points on the fitted H2O+ (X4A″ ) PES a.

The major topographical features of the H₂O⁺ (X⁴A″ ) PES reported in the present work are displayed in Figs.  3– 7. The salient features from these contour plots are some of the most relevant stationary points for the title system. Clearly, a smooth and correct behavior can be observed over the whole configuration space. Figure  3 shows a contour plot for the C_2v insertion of O⁺ into H₂ diatom from large atom– diatom separations along T-shaped geometries. In this process, when the O⁺ inserts into H₂ diatom, the ∠ HOH opens progressively, while the separation between O⁺ and the center of H– H shortens, and the bond length of H– H increases. It can be seen from this figure that there exists a local minimum (C_2v[OH₂]⁺ ) locating at x ≈ 1.5a₀ and y ≈ 3.3a₀ (R₁ = 1.4945a₀, R₂ = R₃ = 3.4241a₀), with the well depth lying at − 0.1954 E_h which is 8.5341  kcal· mol^{− 1} above the C_{∞ v}[OHH]⁺ minimum, and the harmonic vibrational frequencies being 1577, 562, and 1886  cm^{− 1}, respectively. The geometry, well depth, and harmonic vibrational frequencies compare well with the results obtained by Martí nez et al.^[4] Another notable feature shown in this figure is the D_{∞ v}[HOH]⁺ transition state lies at R₁ = 5.1406a₀ and R₂ = R₃ = 2.7053a₀, with the barrier height of − 0.1230 E_h which is 53.338  kcal· mol^{− 1} above C_{∞ v}[OHH]⁺ minimum and 44.804  kcal· mol^{− 1} above the C₂ v[OH₂]⁺ minimum. Comparing with the results of Martí nez et al.:^[4] the bond length differences are calculated to be 0.0816a₀ and 0.0408a₀, and the well depth predicted by the current PES is 3.7651  kcal· mol^{− 1} deeper. The harmonic vibrational frequencies are of D_{∞ v}[HOH]⁺ transition state are calculated to be 1519.61  cm^{− 1}, 2507.02i  cm^{− 1}, and 456.82  cm^{− 1}, with the MRCI/AVQZ results being 1438.66  cm^{− 1}, 3618.64i  cm^{− 1}, and 841.35  cm^{− 1}. In turn, the results obtained from the analytical PES by Martí nez et al.^[4] are 1204.39  cm^{− 1}, 1759.44i  cm^{− 1}, and 675.72  cm^{− 1}, with the CCSD(T)/cc-pvQZ results being 1416.07  cm^{− 1}, 2759.75i  cm^{− 1}, and 860.59  cm^{− 1}. As seen from Table  3, the results from the current PES show overall good agreement with the results of Martí nez et al.^[4]

Figure  4 shows a contour plot for both OH⁺ stretching in H– O– H⁺ linear configuration. The notable features shown in Fig.  4 are the two equivalent asymmetric hydrogen-bonded minima, which are connected by a barrier lying at R₂ = R₃ = 2.7053a₀. Figure  5 shows the contour plot for bond

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. Contour plot for C2v insertion of O+ into H2 diatom. Contours equally spaced by 0.0065Eh, starting at − 0.192Eh.

Fig.  4.

	Figure Option ViewDownloadNew Window
	Fig.  4. Contour plot for bond stretching in [H– O– H]+ linear configuration. Contours equally spaced by 0.004Eh, starting at − 0.1158Eh.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. Contour plot for bond stretching in [O– H– H]+ collinear configuration. Contours equally spaced by 0.0085Eh, starting at − 0.205Eh.

stretching of OH⁺ and HH in O– H– H⁺ collinear configuration. A minimum can be seen to locate at R₁ = 2.2950a₀, and R₂ = 2.1209a₀, which is attributed to be the C_{∞ v}[OHH]⁺ minimum. Also visible in this figure is O⁺ can approach HH diatom via a barrierless process from large atom– diatom separations along the O– H– H⁺ collinear configurations, thus agreeing with findings of Martí nez et al.^[4]

Figure  6 shows a contour plot for an H atom moving around an OH⁺ diatom with the bond length of OH⁺ fixed at equilibrium geometry R_OH+ = 1.943a₀, which lies along the X axis with the center of the bond fixed at the origin. One of the salient features is the C_{∞ v}[OHH]⁺ minimum exists when H approaches OH⁺ collinearly, which can also be seen in Fig.  5. In turn, figure  7 shows a contour plot for O⁺ atom moving around a fixed HH diatom with the bond length fixed at equilibrium geometry R_HH = 1.401a₀, which lies along the X axis

Fig.  6.

Figure Option ViewDownloadNew Window

Fig.  6. Contour plot for H atom with the coordinate (x, y) moving around a fixed OH+ diatom with the bond length fixed at ROH+ = 1.943a0, which lies along the X axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.0008Eh, starting at − 0.21Eh. The dashed lines are contours equally spaced by − 0.00006Eh, starting at − 0.1946Eh. Also shown are the calculated ab initio points for the H atom moving around OH+ diatom.

Fig.  7.

Figure Option ViewDownloadNew Window

Fig.  7. Contour plot for O+ with the coordinate (x, y) moving around a fixed H2 diatom in equilibrium geometry RH2 = 1.401a0, which lies along the X axis with the center of the bond fixed at the origin. Contours are equally spaced by 0.002Eh, starting at − 0.196Eh. The dashed area is contours equally spaced by − 0.0002Eh, starting at − 0.1745Eh. Also shown are the calculated ab initio points for the O+ atom moving around H2 diatom.

with the center of the bond fixed at the origin. As can be seen in this figure, there exist a C_2v and two equivalent C_{∞ v} minima in the process of O⁺ moving around the fixed HH diatom. Also shown in Figs.  6 and 7 are the ab initio energies utilized in the fitting procedure, which cover both the short and long regions, warranting the reliable behavior of the current PES. These two contour plots clearly reveal a smooth behavior both at short and at long-range regions.

A wealth of accurate MRCI/AVQZ energies have been calculated in the present work to map the PES of H₂O⁺ (X⁴A″ ). The resulting raw energies have been corrected subsequently using the DMBE-SEC method and then utilized in the fitting procedure to obtain the globally accurate PES, which is expected to be realistic over the entire configuration space. The various topographical features of the current PES have been examined in detail and compared with the previous theoretical results. The properties of several minima and transition states, including the geometries, energies, and vibrational frequencies, have been characterized on the current H₂O⁺ (X⁴A″ ) PES, which agree nicely with the topographical predictions of other theoretical work. Based on such features, it is concluded that the PES here reported is globally valid while accurately fitting the ab initio points utilized for its calibration. Thus, the present H₂O⁺ (X⁴A″ ) PES should be valuable for studying the kinetics and dynamics of reactions involving ground-state H₂O⁺ , such as which is currently in progress.