The interaction of OH(X²Π) with H₂: Ab initio potential energy surfaces and bound states

There has been considerable interest in the dynamics of the OH–H₂ system. One important reason is that the OH + H₂ → H₂O + H reaction is one of the simplest four-atom reactions.^1,2 The OH(X²Π)–H₂ system is also of astrophysical importance.^3,4

Loomis and Lester have stabilized the weakly bound OH(X²Π)⋅⋅⋅H₂/D₂ complexes,⁵ reporting binding energies of, for H₂, 54 cm⁻¹ and, for D₂, >66 cm⁻¹, which they assigned to the OH–ortho-H₂ and OH–para-D₂ complexes, respectively.⁶ Subsequently, Lester and co-workers have observed and analyzed the rotational structure of these complexes.^7–9 

Andresen et al.^10,11 first determined relative inelastic cross sections of OH j = 3/2 F₁ (with equal populations in both Λ-doublet levels) in collisions with H₂ and D₂. Later, Schreel and ter Meulen¹² measured Λ-doublet resolved inelastic cross sections for OH–H₂, using an electrostatic hexapole as a state selector for OH. Recently, Kirste et al.¹³ used a Stark decelerator¹⁴ to measure the energy dependence of inelastic cross sections for transitions out of the upper j = 3/2 F₁f Λ-doublet level of OH in collisions with D₂.

Here we report two accurate sets of OH(X²Π)–H₂ potential energy surfaces (PES’s) for the study of both bound states and inelastic collisions. Because of its orbital degeneracy, two PES’s are needed to describe the interaction of OH with a closed-shell collision partner.¹⁵ 

To treat the dynamics of the OH–H₂ system, it is most convenient to describe the interaction of the collision partner in terms of the two components (with body-frame projections Λ = ±1 of the electronic orbital angular momentum) of OH(⁠

In this Cartesian basis, the wave functions for the supermolecular system of OH and a structureless atom are eigenfunctions of the same reflection operator, and correspond to two states of differing reflection symmetry (A^′ and A^′′). The electronic energies of these two states correspond to interaction potentials which are conventionally labelled

In general the tetratomic OH–H₂ system is nonplanar. There will be still two electronic states, which we can express as linear combinations of the two Cartesian basis functions:

Two strategies have been previously employed to treat this: The first is to sample only OH–H₂ geometries that have a plane of symmetry, in which case the wave functions for the two electronic states have differing reflection symmetries. This allows the two PES’s, which correspond to the lowest states of A^′ and A^′′ symmetry, to be calculated separately using conventional quantum chemical approaches. This strategy was adopted by Kochanski and Flower (KF),¹⁷ who used self-consistent field (SCF) calculations with a perturbation theory approximation to include electron correlations, and later by Miller et al. (MCKW)¹⁸ in coupled electron pair approximation (CEPA) calculations.

Alternatively, one can obtain the mixing angle directly in the ab initio calculations. This strategy was implemented by Offer and van Hemert (OvH),¹⁵ using the multi-configuration self-consistent field (MCSCF) method with corrections for the dispersion interaction.

The PES’s mentioned above have been tested by comparison with experimental studies of rotationally inelastic scattering and the determination of OH–H₂ bound state energies. Both the KF and MCKW PES’s gave state-to-state inelastic cross sections in good agreement with the experiments by Andresen et al., although the agreement was better for cross sections in which OH remains in its initial (F₁) spin-orbit manifold.^18–20 Cross sections computed from the OvH PES’s agree well with the later scattering experiments of Schreel and ter Meulen.^12,21 The KF and MCKW PES’s predict dissociation energies (⁠

It is probable that these disagreements reflect inaccuracies in the PES’s, either because the ab initio calculations sampled a small number of geometries, or because of a simplistic (by today's standards) treatment of electron correlation.

The two strategies mentioned above have their own shortcomings. In the first approach, one samples a very limited number of orientations. The second tack requires ab initio methods that can compute excited-state energies accurately, which is less cost-effective. In this paper, we present two high-quality sets of OH–H₂ PES’s, one based on a multi-reference configuration interaction [MRCISD+Q(Davidson)] method and the other based on an explicitly correlated coupled cluster method with single-, double-, and (pertubatively) triple-excitations [RCCSD(T)-F12a]. In a similar vein, Dobbyn et al.²³ combined RCCSD(T) and MRCI methods to construct PES’s to describe the Cl(²P) + HCl → ClH + Cl(²P) reaction.

This paper is organized as follows: In Sec. II, we present the theory for the analytical expansion of the OH–H₂ PES’s in a definite-Λ basis. In Sec. III, we describe our treatment in developing the PES’s. We show that despite its incomplete angular sampling, the CCSD(T) PES’s are more accurate. Features of the PES’s and comparison with previous PES’s are also illustrated in Sec. III. We present bound-state calculations for the OH–H₂ complex with our PES’s in Sec. IV. A discussion in Sec. V concludes this paper. Further technical details are available in the supplementary material.²⁴ 

The theoretical framework for the PES’s describing the interaction of a diatomic molecule in a ²Π electronic state, like OH(X²Π), and another diatomic molecule in a ¹Σ⁺ electronic state has been described previously.^15,19,25 In this section, we summarize and clarify the previous treatments.

We use the coordinate system of Fig. 1 to describe the OH–H₂ complex. The OH moiety lies in the XZ plane. We keep the OH and HH bond distances frozen. The interaction then depends on the two polar angles θ_O and θ_H, the azimuthal angle ϕ_H of the HH moiety, and the distance R between the centers-of-mass of the OH and HH moieties.

As discussed in greater detail in the supplementary materials,²⁴ we follow our earlier work on the interaction between a ²Π molecule and a closed shell atom.¹⁶ For the more general case of a ²Π molecule and a closed shell diatom the elements of the 2 × 2 matrix of the interaction potential in the definite-Λ basis can be expressed in the body frame (BF) expansion

where (l₁m₁, l₂m₂|lm) is a Clebsch-Gordan coefficient, and

Since the Hamiltonian is invariant with respect to the chirality of a rotation around

Wormer and co-workers have used the invariance of the Hamiltonian with respect to time-reversal to show that only terms with even values of (l₁ + l₂ + l) occur in the expansion of the diagonal matrix elements with Λ^′ = Λ.²⁵ Hence the electronic Hamiltonian will be invariant with respect to a reflection in the XZ plane in Fig. 1 (ϕ_H → −ϕ_H). Under this constraint, the sum of exponentials in Eq. (2) is proportional to cos (mϕ_H), which allows us to write V_d as

Since the diagonal elements of a Hermitian matrix must be real, and since all the terms in this equation except the

The reader can easily deduce from the work of Wormer and co-workers²⁵ [see their Eqs. (30) and (31)] that the expansion coefficients in Eq. (1) obey the symmetry relation

Because the matrix of the Hamiltonian must be Hermitian, Eq. (4) implies that the

We designate the off-diagonal

The off-diagonal potential for the interaction of a linear molecule in a Π electronic state with a closed-shell diatomic exhibits the same lowering of an underlying symmetry, reflected by the appearance of expansion terms with odd (l₁ + l₂ + l), as the interaction of a nonlinear polyatomic with a closed-shell diatomic.^27,28 Consequently, in fitting our ab initio interaction energies, we did include the flexibility of both even and odd values of (l₁ + l₂ + l) for the off-diagonal matrix elements.

The individual {l₁, l₂, l} basis functions in the expansions in Eqs. (2) and (5) are not normalized with respect to integration over θ_O, θ_H, and ϕ_H. Comparable expansions in orthonormal bases are

and

Examination of the magnitude of these

Ab initio calculations give electronically adiabatic states. For a four-atom ²Π–¹Σ⁺ system in geometries in which there is a reflection plane that lies in the XZ or YZ plane, the electronic wave functions will belong to different irreducible representations A^′ and A^′′ in C_s symmetry. We will designate these states as |Π_x⟩ and |Π_y⟩, in which the Cartesian index designates the singly-filled π lone-pair orbitals on the OH molecule, i.e.,

Consider now an arbitrary OH–H₂ geometry, in which there are no symmetry elements save the identity operator (C₁ symmetry). In an ab initio calculation, it is possible to rotate the orthogonal π molecular orbitals so that one lies in the XZ plane and the other, perpendicular to this plane. In our specific application, the molecular orbitals from a complete active space, self-consistent field (CASSCF) calculation are rotated to achieve maximum overlap with the CASSCF orbitals for the OH–H₂ system at the same geometry except for ϕ_H = 0.

We use these rotated orbitals to define two state functions

where χ is the so-called “mixing angle.” This latter depends on the OH–H₂ geometry, just as the potential energy surfaces.

The 2 × 2 basis in terms of the rotated π orbitals defines what are called “quasi-diabatic” (or, more simply, “diabatic”) states.^30,31 In this diabatic basis, the matrix of the interaction Hamiltonian

where

As with the one-electron molecular orbitals, the definite-Λ states are linear combinations of the Cartesian states, namely,

Wormer and co-workers have shown [see Eq. (4) of Ref. 25] that in the definite-Λ basis the diagonal potential is

and the off-diagonal potential is

For any OH–H₂ geometry in which the XZ plane is a plane of reflection symmetry V_xy and V_yx vanish. For these geometries, the diagonal potentials are identical to the adiabatic potentials V_l and V_u, and it follows that

and

Equations (13) and (14) are equivalent to the “sum” and “difference” PES’s usually used in describing the interaction between a ²Π atom and a structureless atom.^16,29,32

The symmetry of the potential with respect to interchange of the two hydrogen atoms and to reflection in the XZ plane (Fig. 1) allows us to restrict the ab initio calculations to a subspace of the complete coordinate space grid. As a consequence, in the ab initio calculation of the potential energy surfaces, we can restrict ourselves to the subspace 0 ⩽ θ_H ⩽ π/2 and 0 ⩽ ϕ_H ⩽ π.

When the OH molecule lies along the Z axis, θ_O = 0 or π. The two adiabatic potential energies, V_l and V_u, will obviously be independent of ϕ_H.

When the H₂ molecule is in the XZ plane, i.e., ϕ_H = 0, or θ_H = 0 or π,

When the OH molecule lies along the Z axis, the OH–HH system is always planar. However, the CASSCF orbitals are rotated to achieve maximum overlap with those from a calculation in which the XZ plane (Fig. 1) defines the plane of reflection symmetry. Because the system is planar, there is no additional coupling between the two states. In the case where θ_O = 0 the mixing angle χ is equal to ϕ_H, while for θ_O = π we have χ = −ϕ_H.^15,25

Finally, in linear geometries (both polar angles 0 or π), the degeneracy of the two Π states of OH is not lifted. The two adiabatic states of the complex are degenerate.

We used the MOLPRO program suite³³ to determine the two lowest adiabatic energies of the OH–H₂ complex, E_l and E_u, as well as the mixing angle χ. For each geometry, we first performed a CASSCF calculation, with 1 core orbital and 7 active orbitals, and the augmented correlation-consistent polarized quadruple zeta (aug-cc-pVQZ) basis.³⁴ As discussed in Sec. II E, the resulting orbitals were rotated to maximize the overlap with CASSCF orbitals from a calculation with ϕ_H = 0 and the other coordinates unchanged.

With configuration state functions built from these diabatic orbitals, we then performed an internally contracted, multi-reference configuration interaction calculation including all single and double excitations (ic-MRCISD) to determine the two adiabatic energies. The CI calculation included 7 active orbitals and 210 external orbitals. We used the cluster corrected energies of Langhoff and Davidson.³⁵ These MRCI calculations were carried out on a four-dimensional (R, θ_O, θ_H, ϕ_H) grid (see Fig. 1).

We assumed that the OH and HH bond distances remain fixed at the average value of r in the v = 0 vibrational levels [r₀ = 1.8509 a₀ for OH (Ref. 36) and r₀ = 1.448736 a₀ for H₂ (Ref. 37)]. The MRCI method is not size-extensive, i.e., the calculated energy of the OH-H₂ complex does not precisely equal the sum of the energies of OH and H₂ when the distance R goes to infinity. Therefore, we obtained the interaction energy at a given geometry from the expression

where Ω denotes the three angles. Finally, the quasi-diabatization procedure in MOLPRO was used to compute χ(R, Ω) from analysis of the CI coefficients.^38,39 We did not use the alternative, but computationally more cumbersome, approach of diabatization by integration of calculated nonadiabatic coupling matrix elements.

The angular grid included 845 orientations, with 11 values of θ_O defined by cos θ_O = −1 to 1 in steps of 0.2, 6 values of θ_H defined by cos θ_H = 0 to 1 in steps of 0.2, and 13 values of ϕ_H from 0 to π in steps of π/12. This choice of grid provides a uniform sampling of the differential solid angle sin θ₁dθ₁sin θ₂dθ₂ dϕ₂.⁴⁰ The OH–H₂ distance was spanned by 32 values of R ranging from 3.5 a₀ to 16 a₀.

We found that the angular dependence of the computed interaction energy is not smooth beyond R = 16 a₀. We discarded 112 of the 27456 points, corresponding to small values of R, for which there was unacceptable scatter in the calculated values of V_o.

We also carried out restricted explicitly correlated coupled-cluster calculations with single-, double-, and (perturbatively) triple-excitations [RCCSD(T)-F12a].^41,42 The advantage of a coupled-cluster method is that one can recover more of the correlation energy. The disadvantage is that it is difficult to determine the energy of the second adiabatic state, except in high-symmetry geometries where the wave functions for the two states belong to different irreducible representations.

We performed RCCSD(T)-F12a calculations with the augmented correlation-consistent triple-zeta basis (aug-cc-pVTZ)^34,43,44 supplemented with mid-bond functions (3s3p2d2f1g1h) at the following high-symmetry OH–H₂ geometries: (1) ϕ_H = 0 or 180°, and a grid including 11 values of θ_O (0°, 20°, 40°, 60°, 80°, 90°, 100°, 120°, 140°, 160°, and 180°) and 7 values of θ_H (0°, 22.5°, 30°, 45°, 60°, 67.5°, and 90°); (2) ϕ_H = 90°, θ_O = 0° or 180°, and θ_H with each of the 7 values listed above; and (3) ϕ_H = 90°, θ_H = 0° or 90°, and θ_O with each of the 11 values listed above excluding 0° and 180°. We note that many of the geometries listed above are equivalent, e.g., all geometries in set (2) and geometries with θ_H = 0 in set (3) are equivalent to some geometries in set (1). However, we used the interaction energies of all these geometries in fitting the PES’s. This choice of angles gives rise to 186 orientations, of which 131 are unique.

For all these geometries the existence of a plane of symmetry allows us to determine RCCSD(T) energies for the lowest states of both A^′ and A^′′ symmetry. These energies correspond to V_xx and V_yy depending on the chosen plane of symmetry. The potential in the definite-Λ basis was computed with Eqs. (13) and (14). The grid in the OH–HH distance R consisted of 39 values ranging from 3.5 a₀ to 30 a₀.

We used the scaled triples correction⁴² as implemented in MOLPRO.³³ In these RCCSD(T) calculations, we used the counterpoise method to estimate the basis set superposition error (BSSE),⁴⁵ in which we subtract the energies of the fragments computed in the (supermolecular) atomic orbital basis. The interaction energy is

To check the accuracy of our MRCI and RCCSD(T) calculations, we carried out additional computations, with various basis sets and methodologies, for T-shaped geometry at several values of the intermolecular separation R, ranging from the repulsive region to that of the van der Waals well. The two sections of Table I present the results of these additional calculations.

We can draw a number of conclusions from Table I. In the MRCI calculations, there is little difference in the interaction energies when only two electronic states [the two components of OH(X²Π)], rather than the first three states [OH(X²Π + A²Σ⁺)], are optimized. Second, CBS extrapolation with an active space including the O 3s and 3p orbitals yields only a slightly smaller value for the repulsive potential. Finally, use of an explicitly correlated F12 basis gives results similar to the CBS extrapolation.

Calculations with full triples (CCSDT)⁴⁶ give results very similar to those with perturbative inclusion of triples. Unfortunately, it is not possible to perform CCSDT calculations with larger basis sets, nor does a F12 version of the CCSDT code currently exist.

Expanding the calculations to allow core-valence correlation has little effect upon the computed interaction energies. Also, explicit inclusion of scalar relativistic effects (using the Douglas-Kroll Hamiltonian)⁴⁷ also does not significantly affect the results; this is expected since oxygen is a light atom.

In our calculations, we kept the OH and H₂ bond distances fixed at r₀. A better approximation would be to calculate the PES at several values of the bond distances and then average over the ground vibrational wave functions for both species. We see in Table I that vibrationally averaging over r(H₂) has little effect. Vibrational averaging over r(OH) leads to only a slight increase in the repulsive interaction.

Present computational resources do not permit full-triples calculations (CCSDT), or those involving vibrational averaging, for the more than ∼6000 points used to determine the RCCSD(T)-F12 potential. Notwithstanding, from Table I we would not anticipate a significant difference in the calculated interaction energies.

The phase of the CI wave functions is not defined uniquely in the CI calculation. Consequently the signs of the mixing angle χ, as well as of the real and/or imaginary part of the off-diagonal potential, may vary inconsistently from point to point. To compensate for this uncertainty in sign, we first fitted our CCSD(T) ab initio potential to a small set of angular basis functions. This fitted potential was then used to resolve ambiguities in the signs of the MRCI off-diagonal potentials. Manual inspection confirmed smoothness in the off-diagonal PES. Unfortunately, for a few geometries, V_o is very close to zero (|V_o| < 0.05 cm⁻¹), preventing the determination of a consistent sign for these points.

The fitting procedure was as follows: We first set upper limits of 9 and 6 for, respectively, l₁ and l₂. The subsequent expansion contains 126 terms for V_d and 198 terms for V_o. We then carried out a least-squares fit of the sign-corrected MRCI interaction energies at R = 5 a₀ (a point in the repulsive region of the PES where we do not anticipate that the important expansion coefficients will be changing sign) to Eqs. (6) and (7). All terms with expansion coefficients of absolute value greater than 1.5 cm⁻¹ were retained. This corresponds to 39 terms for V_d and 15 terms for V_o. In the latter case, 3 of the terms have odd values of (l₁ + l₂ + l).

We used this (39/15)-term fitting scheme and Eqs. (6) and (7) to fit both the MRCI and the CCSD(T) PES’s. Figure 2 displays the root-mean-square deviation (RMS) of the fits to V_d and V_o as a function of R. Overall, we see the fit is good, and the RMS is well below 1 cm⁻¹ for R ⩾ 6 a₀. Because, in general the off-diagonal potential is smaller in magnitude, the RMS for V_o is smaller. We will discuss this further in Sec. III F.

At short range, V_d of the CCSD(T) potential has a larger rms. This is reasonable since the terms we used in the fit were optimized for the MRCI PES’s. For V_o and V_d at long range, the CCSD(T) PES’s have smaller RMS. The RMS for the MRCI PES’s generally remains a constant (0.05–0.1 cm⁻¹) for R ⩾ 8 a₀. This is a measure of the precision of our MRCI calculations.

To estimate the systematic error of the CCSD(T) PES’s due to the restriction of the angular sampling to high symmetry points, we also created test V_d and V_o PES’s using only the MRCI points whose molecular geometry has a plane of symmetry [referred to in the following as the MRCI-HS (high-symmetry) PES’s]. Out of the 858 orientations for each R, there are 153 that have a plane of symmetry. Fitting just these points, by means of the same procedure as described above, we obtain the MRCI-HS PES’s.

The three sets of OH–H₂ PES’s, in the form of expansion coefficients as a function of R, can be obtained from the supplementary material.²⁴ 

The calculated PES’s were limited to R ⩾ 3.5 a₀, which is sufficiently small for scattering and bound-state calculations. Inside this value of R, the expansion coefficients are held constant to their values at 3.5 a₀. We note that a few expansion coefficients, of relatively small magnitude, do not vary smoothly for R ⩽ 4.25 a₀.

The CCSD(T) PES’s were extrapolated by a least-squares fit of a few ab initio points at large R to the form

Here

The MRCI PES’s were extrapolated by a slightly different method since the ab initio points are less precise at long-range. We employed a multipole expansion method similar to that described in Ref. 25 to model the long-range electrostatic interaction in terms of the multipole moments of the separate OH and H₂ molecules. The relevant theory as well as the calculated multipole moments are included in the additional technical details which are available in the supplementary material.²⁴ 

We included in this expansion only the dipole, quadrupole, and octupole moments of OH, and the quadrupole moment of H₂. For the diagonal PES, we can thus estimate the OH-dipole/H₂-quadrupole (B₁₂₃, R⁻⁴), the quadrupole-quadrupole (B₂₂₄, R⁻⁵), and the OH-octupole/H₂-quadrupole (B₃₂₅, R⁻⁶) interactions. Since there is no dipole moment between the Λ = ±1 and Λ = ∓1 states of OH, for the off-diagonal PES we included just the quadrupole-quadrupole (F₂₂₄, R⁻⁵) and the OH-octupole/H₂-quadrupole (F₃₂₅, R⁻⁶) terms.

For these multipole-multipole interactions, the

We obtain the final MRCI expansion coefficients by smoothly switching, as R approaches 16 a₀, between the fit to the ab initio points at short range and the extrapolated expansion coefficients at long range, for the {l₁l₂l₃} = 123, 224, and 325 terms. For the other {l₁l₂l} terms with known inverse R dependences we extrapolated by smoothly switching at R = 16 a₀ to Eq. (17), or, if not, to an exponential decay.

In the case of the CCSD(T) expansion coefficients, which are smoothly varying out to R = 30 a₀, we did not switch to the known multipole-multipole terms. For R > 30 a₀ we just extrapolated by Eq. (17) for the terms with a known inverse R dependence, or by an exponential decay.

The MRCI-HS PES’s were extrapolated with the same procedure and parameters as for the full MRCI PES’s.

Figure 3 shows several radial cuts of V_d for the MRCI and CCSD(T) PES’s at θ_H = 90° and ϕ_H = 0°. These correspond to planar geometries with H₂ perpendicular to R. Three values of θ_O (0°, 90°, and 180°) were chosen. For comparison, we also show comparable radial cuts of the CEPA MCKW PES’s (adapted from Ref. 18).

The MRCI and the CCSD(T) PES’s agree very well. Of the three orientations plotted, only the θ_O = 0° curve shows a significant well, which corresponds to a T-shaped geometry with the H-end of the OH pointing to H₂. The minima on these curves are also the global minima of V_d, whose positions and depths are tabulated in Table II. The geometry of these global minima are a consequence of both the dipole-quadrupole interaction and the quadrupole-quadrupole interaction. The signs of both quadrupole moments are positive, which implies that there is an excess of positive charge at both ends of each molecule and of negative charge in the middle. Therefore, the quadrupole-quadrupole interaction favors a T-shaped geometry.

Both minima are significantly more attractive than predicted by the MCKW PES’s (see Table II). This is not unexpected, since the MCKW ab initio calculations were based on the CEPA method which gives a less complete recovery of electron correlation and, hence, a shallower van der Waals well. Second, the monomer bond lengths were slightly different. Here, for both OH and H₂ we used the average value of r in the ground vibrational level, namely, ⟨r⟩_{v = 0} for both OH and H₂, as recommended by Faure et al.⁵⁰ Since ⟨r⟩_{v = 0} is larger than r_e, the polarizability and the electrostatic moments of the monomer fragments will be larger, so that the long-range interaction will be stronger.

We display the larger fitted expansion coefficients as a function of R in Figs. 4 and 5 for the MRCI, MRCI-HS, and CCSD(T) PES’s. The difference between the expansion coefficients of the MRCI PES’s and MRCI-HS PES’s are insignificant. Consequently, restricting the angular sampling in the CCSD(T) ab initio calculations to geometries with a plane of symmetry will not lead to significant error in the fitted CCSD(T) PES.

As a justification of the use of only the points with high symmetry, we performed some full close-coupling scattering calculations for several rotationally inelastic transitions of OH in collisions of H₂ j = 0, j = 1, and j = 2, at collision energies ranging from 70 to 150 cm⁻¹. We found that the state-to-state cross sections computed with the MRCI and the MRCI-HS PES’s differ by no more than ∼2%, comparable to the precision of the scattering calculations (∼1%). For calibration, the maximum difference between cross sections computed on the MRCI and CCSD(T) PES’s is 8%. Also, we see from Table II that OH⋅⋅⋅H₂ dissociation energies computed with MRCI and MRCI-HS PES’s differ by only ∼0.3 cm⁻¹. We will further show the OH⋅⋅⋅H₂ bound state energies and spectroscopic constants computed from these two sets of PES’s in Sec. IV.

We thus assert that our CCSD(T) set of PES’s, determined at points with a plane of symmetry, is the most accurate currently available for OH–H₂. Inclusion of triple excitations allows a more complete recovery of electron correlation than in the MRCISD+Q calculations. Use of the CCSD(T)-F12a method with an aug-cc-pVTZ basis is expected to be more accurate than a standard (without explicit correlation) CCSD(T) calculation with a larger aug-cc-pVQZ basis.⁴² In addition, the CCSD(T) calculations are faster. Also, whenever there is a plane of symmetry, the diabatic coupling V_xy vanishes, so that the determination of the diagonal V_d and off-diagonal V_o PES’s from the V_xx and V_yy Cartesian PES’s is straightforward by means of Eqs. (13) and (14).

The top panels of Figs. 4 and 5 show coefficients with l = l₁ + l₂, l₁ ⩾ 1, and l₂ ⩾ 2, to which the multipole-multipole electrostatic interactions make significant contributions at long range. For the diagonal potential, the B₁₂₃ (OH-dipole/H₂-quadrupole, R⁻⁴ dependence) and the B₂₂₄ (quadrupole-quadrupole, R⁻⁵ dependence) terms dominate the PES at medium to long range. In this range the only comparable non-multipole contribution is the isotropic B₀₀₀ term, which is a reflection of the R⁻⁶ dispersion interaction. As seen in the lower panel of Fig. 5, the largest non-multipole contribution is the F₂₀₂ term, which is a difference between the dispersion interaction between H₂ and OH in its

We discussed in Sec. II C, the appearance of terms in the expansion of the off-diagonal potential V_o corresponding to an odd sum of the indices {l₁, l₂, l}. For the electrostatic interactions between the two fragments, the indices obey the relation l₁ + l₂ = l, so they do not give rise to terms with odd (l₁ + l₂ + l). The induction, dispersion, and overlap interactions are not subject to this restriction. Thus, for these interactions there arise terms with l₁ + l₂ ≠ l, both with odd and even (l₁ + l₂ + l). For this reason, and guided by Green's investigation of similar terms in the expansion of the interaction between a non-linear polyatomic and a diatomic,²⁷ we included terms with odd (l₁ + l₂ + l) in our fit. The largest of these was the F₂₂₃ term (see Fig. 5). As can be seen, this term is quite small except at small values of R.

At small R, the anisotropy of both the diagonal and the off-diagonal PES’s is more complex than predicted by a simple multipole-multipole model. A number of expansion terms make important contributions. When H₂ is in the lowest rotational level (j = 0), its wave function is spherically symmetric. The OH–H₂(j = 0) PES is a function of only R and θ_O, and can be obtained by an equal-weight averaging over all orientations (θ_H and ϕ_H) of the H₂ moiety. Figure 6 compares the OH–H₂(j = 0) V_d PES with the V_sum PES’s for OH–He⁵¹ and OH–Ne.⁵² 

The anisotropies of these PES’s are very similar. All three have global minima at θ_O = 0° and R ∼ 6.5 a₀. Although OH–H₂ is isoelectronic with OH–He, the OH–H₂(j = 0) minimum is more than three times deeper than the OH–He minimum, and the interaction is more repulsive at short range. Both differences are reasonable. For H₂ in j = 0, the charge distribution is spherical so that there will be no multipole moments. Hence, the dominant attractive contribution is the dispersion interaction. Because the charge distribution of H₂ is more polarizable than that of He, the OH–H₂ attractive interaction will be larger. Similarly the size of an H₂ molecule, even averaged over orientation, will be larger than that of the isoelectronic He atom. Thus, we expect a somewhat steeper repulsive wall. In reality, the OH–H₂(j = 0) PES’s appears more similar to that of OH–Ne.

Figure 7 displays contour plots of the off-diagonal OH–H₂(j = 0) PES and the V_dif PES’s for OH–He and OH–Ne. The V_o (or

Figure 8 presents contour plots of V_d and V_o in the OH molecule frame similar to plots given by OvH (Fig. 3(a) of Ref. 15). In these plots, the center of mass of the OH molecule defines the origin. The OH molecule lies on the z axis with the H atom on the positive direction. The OH–H₂ separation is (x² + y²)^1/2, and the orientation of the OH molecule in the space-fixed frame defined in Fig. 1 is θ_O = tan ⁻¹(−x/z). The orientation of the H₂ molecule is fixed at {θ_H, ϕ_H} = {27.464°, 45°}. Overall, our PES’s and those of OvH¹⁵ show a very similar dependence on the OH–HH distances and on the orientation of the OH moiety. For both sets of PES’s, V_d has an attractive well on the oxygen side of OH, and |V_o| is slightly asymmetric about the x = 0 axis.

Our CCSD(T) calculations sampled only high-symmetry geometries while the OvH PES’s are fits to ab initio calculations which sampled non-planar geometries. Thus the great qualitative similarity in the OvH and CCSD(T) PES’s confirms that sampling only high-symmetry geometries will provide an excellent description of the OH–H₂ system. The OvH PES’s were computed with a SCF + dispersion method, which yields a deeper well at this H₂ orientation compared with our more accurate method.

Contour plots showing the dependence of V_d and V_o on θ_O and θ_H for ϕ_H = 0° and 90° at R = 6 a₀ are presented in Figs. 9 and 10, respectively. At medium to long range, the interaction between the OH dipole (negative at the O end) and the HH quadrupole (positive at the ends, negative in the middle) will dominate the interaction. The interaction will be most attractive when the OH approaches the H₂ in T-shape geometry with the H-end of the OH closer to the H₂ (θ_O = 0°, θ_H = 90°, T-shape OH–HH). In addition to this global minimum, there are local minima in linear HO–HH geometry, in which the negative end of the OH dipole approaches the positive ends of the HH quadrupole (θ_O = 180°, θ_H = 90°, linear HO–HH).

Figure 10 presents contour plots of the off-diagonal potential V_o. The minima and the maxima of ℜ(V_o) and ℑ(V_o) are all located at θ_O = 90°. At this value of θ_O, one lobe of the π_x orbital of OH points to H₂, and thus will maximize the difference between the energies of the

When ϕ_H = 90°, exp (imϕ_H) = i^m and cos (mϕ_H) vanishes unless m is even. Since the index l₂ is also even (because of the permutation symmetry of the HH moiety), one can show that the diagonal potential [Eq. (3)] and the real part of the off-diagonal potential [Eq. (5)] is unchanged when θ_H is replaced by π − θ_H. In contrast, the imaginary part of the off-diagonal potential changes sign when θ_H → π − θ_H. Thus V_d and ℜ(V_o) should be symmetric, while ℑ(V_o) should be antisymmetric, about θ_H = 90°. This prediction is borne out by the right panel of Fig. 9, and the middle and right panels of Fig. 10.

Finally, we compare in Fig. 11 at R = 15 a₀ contour plots of the contribution to V_d of just the terms which correlate at long-range to the lower-order multipole-multipole interactions (B₁₂₃, B₂₂₄, and B₃₂₅) and the full fitted CCSD(T) PES’s. The electrostatic interactions were computed from the multipole moments listed in Table I of the supplementary material.²⁴ Overall, consideration of just the electrostatic interactions represents well the angular dependence and magnitude of the PES. Notwithstanding, the CCSD(T) PES’s is less repulsive at all orientations than that predicted by with the electrostatic interactions alone. We infer that there are still non-trivial contributions from induction and dispersion at this large value of R. We thus believe that a full calculation of the long-range potential will be more accurate than relying on electrostatics alone.

We note that Wormer et al.²⁵ used electrostatic interactions only in extrapolating their OH–HCl PES. Notably, though, the OH-HCl PES is dominated by the strong R⁻³ dipole-dipole interaction which is absent for OH–H₂.

Although we did not use the multipole moments in extrapolating the CCSD(T) PES, fits of the CCSD(T) ab initio points at large R actually give quite similar results. For example, from the multipole moments listed in the supplementary materials,²⁴ we predict B₁₂₃ = 2.86 × 10⁵(R/a₀)⁻⁴ cm⁻¹, compared to a coefficient of 2.66 obtained by direct fit of the ab initio points. This good agreement justifies our direct fit of the CCSD(T) points to Eq. (17).

The Hamiltonian for the OH–H₂ system with fixed bond lengths can be written as

where μ is the OH–H₂ reduced mass,

We construct a set of basis functions from the product of stretching functions and angular functions and use a set of distributed Gaussians to span the R-space.⁵³ The construction of the rotational basis functions, which can also be used for scattering calculations, is discussed in detail in the supplementary material.²⁴ 

We employed the following parameters for

The expression for the matrix elements of the interaction potential in this angular basis is given in the supplementary material.²⁴ We used our Hibridon⁵⁶ suite of programs, recently extended to handle the OH–H₂ system, to construct and diagonalize the Hamiltonian matrix to determine the energies of the OH–H₂ bend-stretch states. The bound-state energies were also computed independently using another program described in Ref. 57. The results obtained from the two calculations are practically identical.

We list in Tables III and IV the energies of the bend-stretch levels of the OH–H₂ complex calculated with the CCSD(T) PES’s for total angular momenta J = 1/2 and 3/2. In computing these energies, we included all angular basis functions with j₁ ⩽ 13/2, and j₂ ⩽ 5 and ⩽ 4 for ortho and para H₂, respectively. The stretch basis (the expansion in the OH–H₂ distance) was 41 equally spaced Gaussian functions spanning the range 4 a₀ ⩽ R ⩽ 16 a₀. The computed bound state energies are converged to ∼0.01 cm⁻¹.

We have inspected the wave functions for each bound level and assigned OH–H₂ stretching quantum numbers v_s and the body-fixed projection P of the angular momentum of the complex J. P was assigned using the vector model of angular momenta similar to the method described in Ref. 22.

Overall, we find that the computed bound level energies are in qualitative agreement with the values reported by Miller et al.:^18,22 The number of bound levels is similar, the zero point energies are large, and for both J = 1/2 and 3/2 the binding energies for OH–ortho-H₂ are greater than for the para–H₂ complex.

The stronger binding with ortho-H₂ is a general phenomenon for all H₂ complexes.^{18,22,58–63} The para–H₂ molecule is spherically symmetric and cannot orient to sample the most attractive orientations (at least without mixing in the higher j = 2 rotational level). Alternatively, because of this spherical symmetry, the j = 0 state of H₂ lacks a quadrupole moment. Thus, the long-range quadrupole-quadrupole attraction does not contribute in the interaction of OH with H₂(j = 0).

In Sec. III F 1, we discussed whether restricting the angular sampling to geometries with a plane of symmetry would lead to significant error in the fitted PES. The differences between the expansion coefficients for the MRCI and MRCI-HS PES’s were found to be insignificant. Moreover, there were only minimal differences between state-to-state cross sections computed with these two PES’s. As a further check on the validity of considering interaction energies only for geometries with a plane of symmetry, we compare the energies of bend-stretch levels computed with these two PES’s. Table V presents comparison of the energies of the lower J = 1/2 and 3/2 bound levels computed with the MRCI and MRCI-HS PES’s. We see that the energies computed with the MRCI PES’s are slightly lower than the MRCI-HS values, but only by a fraction of a wavenumber (ΔE < 0.3 cm⁻¹). This provides additional justification for using only geometries that have a plane of symmetry in our fit. This restriction permits the use of RCCSD(T) calculations of the interaction energies.

We see in Tables III–V that there are significant parity splittings in the energies of the bend-stretch levels. Green and Lester⁶⁴ and Dubernet et al.⁶⁵ have employed a perturbation theory to determine the J-dependence of the bend-stretch energies and parity splittings, which can be expressed for P = 1/2 bend levels as

where B is the rotational constant of the complex and p is the parity splitting parameter of that particular bend level. As expected from the analysis in the above cited studies,^64,65 the parity splittings are much larger for P = 1/2 levels than for P = 3/2 levels.

We have fit the J and parity dependence of the computed energy of the lowest bend-stretch level (P = 1/2) of OH–ortho-H₂ to Eq. (19). Table VI compares the rotational constant and parity splitting parameter for this bend-stretch level, determined both experimentally and theoretically. As expected from the small differences between the MRCI and MRCI-HS bend-stretch energies, the derived spectroscopic constants computed with these two PES’s are in good agreement. We also see that there are only slight differences in the parameters computed with the CCSD(T) and MRCI PES’s. By contrast, the agreement of our computed spectroscopic constants with those reported experimentally,⁷ which were obtained from fitting J = 1/2 and J = 3/2 level energies derived from the spectrum, is less satisfactory. Our values for B and p are, respectively, smaller and larger than the experimental values. The differences lie outside the experimental error bars. The spectroscopic constants computed with the less sophisticated MCKW PES’s are in better agreement with experiment. We note that the Coriolis coupling in the J = 3/2, P = 1/2 levels we used in deriving spectroscopic constants are weak.

Calculation of the bend-stretch energies allows us to predict the zero-point corrected dissociation energy D₀. These values are presented in Table II. Our predicted well depths (D_e) are ca. 31 cm⁻¹ larger than predicted for the MCKW PES’s. However, this difference is compensated somewhat by the larger curvature of our PES’s in the region of the well (see the blue lines in Fig. 3). Thus the zero-point-corrected dissociation energies (D₀) are only 8 cm⁻¹ (MRCI) and 11 cm⁻¹ [CCSD(T)] larger than the values computed for the MCKW PES’s.

The D₀ values obtained from our CCSD(T) PES’s, 36.1 cm⁻¹ and 53.7 cm⁻¹ for OH–para-H₂ and OH–ortho-H₂, respectively, are both in good agreement with the MRCI values (Table II). We note that our computed D₀ for OH–ortho-H₂ is in excellent agreement with the experimental value of 54 cm⁻¹.⁶ We have also computed the bound states for OH–ortho-D₂ and OH–para-D₂ with the CCSD(T) PES’s, and obtained D₀ values of 51.4 cm ⁻¹ and 67.5 cm⁻¹. These values are ∼7 cm⁻¹ larger than those obtained with the MCKW PES.⁹ The lower limit of the OH–para-D₂ binding energy determined via electronic spectroscopy is 66 cm⁻¹,⁶ which is again in excellent agreement with our prediction.

We have presented here two sets of OH(X²Π)–H₂ PES’s, calculated using the MRCISD+Q (Davidson) and the CCSD(T)-F12a methods. Comparison of the full MRCISD+Q PES’s with a fit to a subset of points determined at geometries with at least one plane of reflection symmetry shows that sampling based on this subset of points provides a very accurate, but much less computationally demanding, representation of the OH–H₂ PES’s. For these points with symmetry we can then use the CCSD(T)-F12a method, which is more accurate and, in addition, computationally faster. We expect this strategy to be useful in developing PES’s for similar systems involving a ²Π molecule and a ¹Σ⁺ molecule, such as NO(X²Π)–H₂.

We have performed calculations of the bend-stretch levels of the binary OH–H₂ complex for both the para and ortho nuclear spin modifications. The dissociation energies (D₀) predicted by fits to the MRCI and the CCSD(T) points differ by less than 2 cm⁻¹, and agree extremely well with the experimental estimate.⁶ This agreement is a measure of the accuracy of our PES’s, especially for the depth and shape of the attractive well.

Comparison of inelastic OH–H₂ scattering cross sections with experiment would constitute a further test of our PES’s. This comparison would be most sensitive to the PES’s near the onset of the repulsive well. Collisions of OH with H₂ have been well studied.^8,12 Most recently, use of a Stark decelerator^13,66 has permitted the determination of the relative magnitude of state-resolved OH(j → j^′) transitions in collisions with D₂ and H₂ over a wide range of collision energies. Comparable scattering calculations are in progress and will be reported soon.

M.H.A. and P.J.D. are grateful for support from the U.S. National Science Foundation (Grant No. CHE-1313332). J.K. acknowledges XSEDE.org for computational time awarded on Lonestar and Gordon high-performance computing facilities through Grant No. CHE-130120. Parts of the calculations were also performed on the Deepthought cluster at the University of Maryland.