Three states global fittings with improved long range: singlet and triplet states of H⁺₃

Abstract

Full dimensional analytical fits of the coupled potential energy surfaces for the three lower singlet and triplet adiabatic states of H⁺₃ are developed, providing analytic derivatives and non-adiabatic coupling matrix elements. The fits are highly accurate and include an improved description of the long range interactions, including new terms for the description of the long range in the diatomic fits and the atom–diatom dissociation channels. The fits are based on the DIM formalism including three body terms in Hamiltonian matrix elements, each of them obeying S₂ permutational symmetry, where the positive charge is placed in either of the three hydrogen atoms, but the full system obeys S₃ permutational symmetry, invariant under all permutations of the nuclei. The ab initio points used in the fitting are obtained from a complete basis set extrapolation, made for all electronic states. Total root mean square errors of the fits are 27 and 12 cm⁻¹, for the singlet and triplet states, respectively. The errors in the channels are lower than 2 cm⁻¹ and 6 cm⁻¹ for the H + H⁺₂ and H⁺ + H₂ channels respectively. The new fits have been used to calculate the rovibrational bound states of the lowest singlet and lowest triplet states showing very good agreement with previous calculations in the literature.

---

1 Introduction

H⁺₃ plays a central role in astrochemistry, being the most abundant ionic molecule in the Universe and the main protonator, leading to the formation of many hydrides, which constitutes the first step of chemistry in the interstellar medium.^1–3

In the past two decades much effort has been made experimentally trying to determine the rovibrational bands in the lab including higher and higher energies and vibrational hot bands⁴ of the ground singlet electronic state. The same amount of effort has been made by theoreticians in improving the description of the ground electronic state,^5–8 aiming to reproduce the experimental rovibrational energy levels and bands using local⁹ and global potential energy surfaces. In addition, many works have been devoted to the description of the ground and excited electronic states of H⁺₃ molecule in the attempt of obtaining a global surface that could be used for spectroscopic calculations as well as reactivity simulations. The description of ground and excited states of H⁺₃ is crucial for reactivity simulations. The ground state dissociates into H₂ + H⁺ while the first excited state dissociates into H⁺₂ + H. Both states cross when H₂ intermolecular distance increases, so for vibrational states of H₂ with v ≥ 3 the collision simulations would require the excited state dissociating into H⁺₂ + H. In addition, the description of the H⁺₂ + H channel is of particular interest since the formation of H₂ by collisions between H⁺₂ and H are very important in the Early Universe, where cosmic grains are not present. Another key point is the description of the long range interactions. Most of the works including the long range interactions have only considered the ground electronic states.^8,10

The description of the lowest^11–13 and first excited triplet^14,15 states have been the main goal of several studies with special attention to the degeneracy at the conical intersection between the two first triplet states.^16,17 An extension to the standard single-valued double many-body expansion was developed by the group of Prof. Varandas^12,14,16 for the representation of the double-valued potential energy surface ensuring the degeneracy at the intersection line. This same method^15,17 has been used recently to analyse the effect of the geometrical phase, at the conical intersection, on the vibrational state calculations of the lowest triplet states showing a change of sign after a full circle in the hyperspherical angle ϕ.

In this work we use a formalism very similar to that developed by Viegas and coworkers,^6,7,18 based on the diatomic-in-molecules (DIM) approach to represent the electronic matrix. The 3 × 3 DIM based matrix of singlet and triplet states can be used in the triatomic-in-molecules (TRIM) formalism used for the fitting of larger hydrogen-like molecules as H⁺₄¹⁹ and H⁺₅.²⁰ The better the fitting of the H⁺₃ electronic states are, the better the description of the channels of H⁺₃ + H and H⁺₃ + H₂ would be, crucial for the simulation of astrophysical processes.

In this paper we provide with a multisheeted potential fit of very high accuracy up to energies above the complete dissociation including the long range interactions, for the three lowest singlet and triplet states. Extra terms are added in the description of the long range of the diatomic molecules (H₂, H⁺₂) and for the long range of the two triatomic channels (H₂ + H⁺ and H⁺₂ + H). These terms not only add accuracy in the long and medium range regions of the potential surface, but they also help to restrict the area of the surface in which the many body fitting applies, thus allowing higher global accuracy.

In Section II, devoted to the description of the theoretical methodology, are described the ab initio calculations and the fitting procedure. Section III presents an analysis of the potential energy surfaces separated in individual subsections for singlet and triplet states. The definite test that proves the accuracy of the fit near the minimum is the calculation of vibrational states, those are presented and compared with previous calculations in Section IV.

2 Theoretical methodology

2.1 Ab initio calculations

In this work we calculate new energy points for the three first singlet states of C_s symmetry (¹A′) and for the three first triplet states (for a global idea of the energetics of the system see Fig. 1). The ab initio calculations use the Dunning augmented correlation-consistent polarised basis sets, aug-cc-pVXZ, where X runs from 2 to 6, to estimate the complete basis set (CBS) energies. The total energies are fitted using an extrapolation scheme based in the results of Dunning and coworkers^21–23 and probed and used already in the ground state of this system.²⁴ For the ground state of H⁺₃ a single two point extrapolation scheme, using aug-cc-pV5Z and aug-cc-pV6Z basis set, appeared to give accurate and reliable CBS energies according to the eqn (3) of Velilla and coworkers.²⁴ In this work the extrapolation procedure has been used in singlet excited states as well as the triplet states. A total of 16 000 energy ab initio points are calculated for the six states using Molpro²⁵ package of programs. The points are calculated using the complete active space self-consistent field/multireferece configuration interaction (CASSCF/MRCI) scheme using C_s symmetry. The active space consists of 2 electrons in 5 orbitals, 4 orbitals of A′ symmetry and 1 orbital of A′′ symmetry. This active space for this system produces results very close to the full CI. In Table 1 are listed some asymptotic and stationary points using the aug-cc-pV6Z basis set and the extrapolation to CBS procedure, both are compared with the most accurate values from literature. The fits presented below are made using the extrapolated values.

Fig. 1 Stationary points along the minimum energy path for the three singlet states and the two lowest triplet states, the third triplet does not present a global minimum. The energy is given in wavenumbers considering the zero of energy in the global dissociation. The structure of the molecule in each stationary point is given by balls and sticks representations.

Table 1 Values of the energies calculated at MRCI level using the aug-cc-pV6Z (aV6Z) basis set and the extrapolated to complete basis set (CBS) limit for hydrogen atom, diatomic molecules of H₂ and H⁺₂ and the triatomic system H⁺₃ at the value of the minimum. Values are compared with the most accurate values in literature. Distances are given in Bohr and energies in Hartrees

System	r e	aV6Z	CBS	Literature
H (^2S)	—	−0.499999	−0.500003	−0.500000
H2 (X^1Σ^+g)	1.4	−1.174360	−1.744473	−1.744476
H^+2 (X^2Σ^+g)	2.0	−0.602632	−0.602643	−0.602635
H^+3 (1^1A′)	1.65	−1.343718	−1.343837	−1.343836

---

2.2 Three states global fitting procedure

The DIM approximation has been used as an starting point for several hydrogen-like systems, as H⁺₃,^18,26 H⁺₄ ¹⁹ and H⁺₅.²⁰ In particular for the H⁺₃ system it has been used to fit very accurately the ground state²⁶ and the first three singlet states.^10,18,27 The three states fit used in this work is very similar to the one proposed by Viegas and coworkers¹⁸ in which extra terms are added inside the DIM matrix at the diabatic level, to improve the accuracy in certain regions of the surfaces, and the adiabatic states are obtained after diagonalisation. Although both works are similar, there are several differences between them. We list here the main differences:

1. The full system must obey the S₃ permutational symmetry. However, each element must obey S₂ permutational symmetry since each term of the matrix corresponds to an electron hole in one hydrogen atom. In Viegas and coworkers work they claim that the same three body term must be added into the diagonal elements of the matrix to fulfill this condition, but we show here that the S₃ permutation symmetry for the full system is obeyed even when different three body terms are added into the diagonal elements of the matrix.

2. We add several terms into the off-diagonal elements that fulfill the S₂ permutational symmetry. The three body terms added into the off-diagonal elements ensure a very accurate adiabatic description of the non adiabatic coupling matrix elements (NACME).

3. The two-body long range term is not enough to ensure an accurate description of the intermediate and long distances regions. A new three-body long range term has been added into the fitting for the improvement of the long range of the singlet states.

The full 3 × 3 electronic matrix describing the diabatic states of H⁺₃ is written as

where the H̄_DIM is the DIM matrix and V̄^(3B)_ABC is the three-body matrix. The DIM matrix elements for singlet and triple states are very similar, and are written only once in a compact form, denoting in parenthesis the difference corresponding to triplet states, as^11,26where all terms V⁽²⁾ are the fits of the potential curves of the respective diatomic molecules. The analytical representation of the diatomic terms follows the same scheme as the used in other hydrogen-like cations^10,19,26,28where ρ(R) = R exp(−βR) are Rydberg functions, and the term V_long is the long range interaction where R̃ is been modified to eliminate the divergence for R = 0, avoiding the use of switching functions, by

R = R + R₀ exp(−(R − R_e)) (4)

with R₀ being a distance parameter and R_e the equilibrium distance.

The long range interaction, for the ground and first excited states of H₂, is expanded in an asymptotically converging series²⁹

V_long(R) = −6.499027R⁻⁶ − 124.4R⁻⁸ − (1135 + 2150)R⁻¹⁰ − 3986R⁻¹¹ − (+)0.818 R^5/2 exp(−2R) + ⋯ (5)

The long range interaction between a proton an a H atom is represented by the spherical charge-induced dipole and charge-induced quadrupole dispersion,³⁰ already used^10,19,20,28 in these hydrogen cation type of systems. An extra term (J) is included with respect to our last fit of the ground state of H⁺₃ ¹⁰where J = −(2/e)Re^−R and is positive for the gerade state and negative fo the ungerade. This term takes into account the exchange splitting of the interaction energy at large interatomic distances, due to the approximation of this energy as a series of the inverse powers of R.^31–34

The second term of eqn (1) is the sum of the three body terms and the long range contributions. Each element of the V̂_ABC matrix is expressed as

V^ABC_ij = V^(3)k_ij + V_long(R,r,θ) × δ_ij. (7)

The long range term in eqn (7), only added to the diagonal terms, depends of the channel considered, H₂ + H⁺ or H⁺₂ + H. This long range term is a function of Jacobi coordinates since these coordinates are better adapted to the interaction between an atom and a diatomic molecule.¹⁰ Considering the channel H₂(¹⁽³⁾Σ⁺_g(u)) + H⁺, the long range interaction includes two contributions

wherewhere P₂ and P₄ are the Legendre polynomials, r′ = r exp(−r), R′ = R + RLIM exp(−α(R − R_e)) (RLIM = 20, α = 1 and R_e = 1.4 for the singlets and RLIM = 10, α = 0.5 and R_e = 4.0 for the Triplets). R is the distance between the centre-of-mass of the diatomic molecule (H₂ or H⁺₂ depending of the channel) and the nucleus (H⁺ or H depending of the channel), r is the interatomic distance, and θ is the angle between R and r. In the first long range term of eqn (9) the term in R⁻³ is the anisotropic charge-quadrupole dispersion energy and the term in R⁻⁴, which has spherical and anisotropic components, is the charge-induced dipole dispersion energy, α₀(r), α_‖(r) and α_⊥(r) are the average, parallel and perpendicular polarizabilities of H₂, respectively, and Q₂(r) is the quadrupole moment of H₂. The variation of H₂ polarizabilities and the quadrupole moment with the r internuclear distance is calculated numerically as the energy derivative using the extrapolated CBS energies. This magnitudes are as well calculated for the terms required in the DIM matrix elements of the triple states.

The Legendre polynomial of fourth order, P₄, is added to correct the angle dependency that is introduced by the DIM approximation at intermediate distances. The values of the coefficients are optimised fitting the points of only the ground state of H⁺₃. These optimised coefficients are maintained constant during the three states fitting.

The long range interaction for the channel H⁺₂ + H is the same for the fitting of singlet and triplet states. The term includes two contributions, one describing the interaction between a proton and an H-atom (the two terms of eqn (6)) and a multiple expansion as the one in eqn (9).

In the DIM matrix, each diagonal element i corresponds to a charge in atom i, that imposes a permutation symmetry on the V^ABC_ij elements added to the matrix. Thus the relevant three body elements are

where LMAX denote the number of three body terms, LMAX = 2 in this study for both, singlets and triplets. The three body terms V_ABC(R_AB,R_AC,R_BC) are sums of products of Rydberg functions,³⁵ where the generic terms in these expressions are given by where n + m + p ≤ M, and d_nmp = d_mnp to ensure the S₂ permutational symmetry of the three-body terms, meaning that V^L_ABC(x,y,z) = V^L_ABC(y,x,z).

All the elements of eqn (11) fulfill the S₂ symmetry explained above, one of the main differences between the fitting in this work and the fitting of Viegas and coworkers.¹⁸

In this work we generate the analytic full 3 × 3 diabatic electronic matrix, in eqn (1), and the first derivatives as a function of the internuclear distances. The derivatives of the adiabatic potentials and the non-adiabatic couplings among all the states can then be obtained using the Hellmann–Feynman theorem, given in eqn (8) and (9) of our previous work.²⁸ All of them are analytical, depending on the analytical terms of the diabatic matrix and the eigenvectors obtained in the diagonalisation. All these derivatives with respect to the internuclear distances can be transformed to any other set of coordinates applying the chain rule.

3 Results

3.1 Analytical fit of singlet states

The three lowest singlet states of the H⁺₃ molecule present several features that makes the fitting very challenging. The system presents an avoided crossing in the entrance channel due to the swapping of the positive charge between the initial fragments. There are saddle points appearing at symmetric linear geometries and the first and second states are degenerate for D_3h geometries. The multisheeted fitting presented in this work can reproduce well the avoided crossings, degeneracies and symmetry of the system, intrinsically included into the DIM matrix.

The channels considering the approaching of an atom to a diatomic molecule (H + H⁺₂(X³Σ⁺_g) or H⁺ + H₂(X¹Σ⁺_g)) in the three singlet states lie below the complete dissociation (H + H + H⁺) of the system. In this fit are included energies up to −0.9E_H, energy above the global dissociation. Adding points at higher energies ensure that main features of the three states are well described at very short distances. Considering energies only below −0.9E_H, the number of points used for the fitting are 15 331, 11 564 and 2555 for the first, second and third singlet states respectively. The order of the expansion in products of modified Rydberg functions is M = 10. The points up to 14 000 cm⁻¹ above the global minimum of the ground state have been weighted with a value of 25 to make sure that the vibrational states of the global minimum present errors not higher than 1 cm⁻¹. In Fig. 2 are plotted the residuals as a function of the energy where the zero is placed at the global minimum of the ground state.

Fig. 2 Residuals (algebraic deviations) of the three singlet PESs vs. the energy (in cm⁻¹). The black dashed vertical line specifies the energy of the complete dissociation (H + H + H⁺). The horizontal black lines are the rms deviation (27 cm⁻¹).

The errors are shown in Table 2 in intervals of energy and Jacobi distance R for the present fit and the DIM matrix with the long range terms. For the upper states we only include the total error. The errors in intervals of the Jacobi R coordinate include points in which the Jacobi r coordinate varies within values 1.2 a₀ ≤ r ≤ 1.6 a₀, giving an idea of the error in the entrance channel H₂ + H⁺. The error in intervals of energy shows a small rms for energies up to the energy of the linear saddle points, about 14 000 cm⁻¹, the error increases along with the energy increasing. More detailed characterization of the three sheets would show that the higher errors do not necessary means a poor fit, since it reproduces all the stationary and crossing points correctly. Including in the fit simultaneously the three states add problems such as: larger convergence times and higher number of parameters. Higher degrees means larger number of parameters which would slow down the calculation of the potential energy in the dynamical calculations. The first excited state has an rms error around 24 cm⁻¹, similar to that of the ground state, while the error of the third state is slightly higher, of 48⁻¹. For extra details on the rms in the long range region (H₂ + H⁺) see section A in the ESI.†

Table 2 RMS obtained for several intervals of distances and energies using the DIM approach with the new long-range terms and the global three-state fitting. For the two excited states states is only given the total rms

	DIM + LR	3-State fit
	rms/cm^−1	rms/cm^−1
Distances/a0
3.0 ≤ R ≤ 7.0	1955.5	13.3
7.0 ≤ R ≤ 15.0	8.9	3.9
15.0 ≤ R ≤ ∞	0.07	0.07
Energies (1^1A′)
Minimum to linearity	3255.4	3.1
Minimum to H2 + H^+	2446.8	17.0
H2 + H^+ to H^+2 + H	837.2	19.0
H^+2 + H to H + H + H^+	1664.8	29.8
Energies (2^1A′)
All points	2856.2	23.6
Energies (3^1A′)
All points	646.1	47.9

---

The stationary points for the three singlet electronic states are given in Table 3 for ab initio energies, the multi-sheet fitted values and some previous results from the literature.^18,36,37 The results published by Röhse and coworkers used an optimised basis set, (16s10p8d)/[10s8p3d], in configuration interaction calculation with explicit linear r₁₂-term (CISD-R12), equivalent to a full CI for a two-electron system. The absolute and relative errors claimed by the authors is of a few micro-hartrees in the relevant regions. The other set of results published by Viegas et al.¹⁸ applied the same CASSCF/MRCI procedure followed in this work using MOLPRO package.²⁵ However, they used the correlation consistent cc-pV5Z basis set of Dunning. In their work they added 187 points of very high accuracy, calculated by other groups.^38,39 The results published by Wolniewicz³⁷ are one of the most precise calculations for the molecule of H₂ which include relativistic and adiabatic corrections.

Table 3 Energy comparison of the stationary points present in the singlet electronic states. Distances are given in Bohr and energies are given in Hartrees

r 12	r 13	r 23	CBS	Fit	Literature
a Röhse et al.^36 b Viegas et al.^18 c Wolniewicz.^37
Ground state (X^1A′)
Global minimum (D3h)
1.65	1.65	1.65	−1.343837	−1.343839	−1.343834^a
Linear saddle points (D∞h)
1.54	1.54	3.08	−1.278687	−1.278756	−1.278682^a
H2 + H^+ channel
1.40	∞	∞	−1.174473	−1.174473	−1.174476^c
Avoided crossing between 1st and 2nd states
2.5	10.0	10.0	−1.095046	−1.095094	−1.093938^c
First excited state (2^1A′)
Global minimum
1.992	7.5	9.492	−1.103053	−1.103015	−1.102894^b
Saddle points (C2v)
2.0	8.16	8.16	−1.103041	−1.103042	−1.102882^b
Saddle points (D∞h)
5.423	5.423	10.846	−1.001782	−1.001743	−1.001319^b
H^+2 + H channel
2.0	∞	∞	−1.102647	−1.102647	−1.102625^b
Second excited state (3^1A′)
Global minimum
11.04	11.04	11.04	−1.000197	−1.000198	−1.000097^b
Global dissociation (H + H + H^+)
∞	∞	∞	−1.000008	−1.000008	−1.000003^c

---

The ground singlet state has a global minimum in a D_3h geometry, the extrapolated ab initio value used in this work is 0.66 cm⁻¹ lower that the one from Röhse et al.³⁶ This state presents three equivalent saddle points for linear configurations 14298.77 cm⁻¹ above the minimum, underestimating the best from literature³⁶ by 0.44 cm⁻¹. The dissociation energy into the H₂(X²Σ⁺_g) + H⁺, coinciding with the minimum of the H₂ molecules, is 37170.99 cm⁻¹ above the global minimum, overestimating the exact dissociation energy³⁷ in less than 1 cm⁻¹. As it has been already mentioned, the two lowest singlet states present an avoided crossing when the inter-atomic distance of H₂ (X²Σ⁺_g) and H⁺₂ (X ³Σ⁺_g) is equal to 2.5 a₀, explained by the swapping of charge between both diatomic systems (H₂ and H⁺₂). This crossing appears 17432.21 cm⁻¹ above the H₂(X²Σ⁺_g) + H⁺ dissociation channel.

The first excited state presents the global minimum in a linear configuration with two of the atoms in a shorter distance than the equilibrium of H⁺₂ and the third atom far away. This state dissociates into the H⁺₂ (X²Σ⁺_g) + H(²S) channel, 22615.76 cm⁻¹ below the complete dissociation. There is a shallow well, for a linear geometry, just 89.11 cm⁻¹ below the dissociation limit corresponding to H⁺₂ + H. The depth of this well is 32.04 cm⁻¹ deeper than that reported by Viegas et al.¹⁸ The three equivalent minima are separated by two types of saddle points, three of C_2v symmetry only few wavenumbers above the minimum and three of D_∞h 23 000 cm⁻¹ above the minimum.

The first and second excited states cross for all D_3h configurations of H⁺₃. The minimum of the third single state lies along the intersection line with a depth of 41.48 cm⁻¹ with respect to global dissociation (H + H + H⁺).

In Fig. 3 are plotted two-dimensional cuts of the potential energy surfaces using stereographic projections of hyperspherical coordinates.⁴⁰ The zero of energy is placed in the complete dissociation (H + H + H⁺). The left/middle/right panels correspond to ρ = 2.18, 6.6 and 10.2 a₀, respectively. The minimum of the well of the ground electronic state is at ρ = 2.18 a₀, and the avoided crossing between the first two singlet states, appearing for Jacobi R distances larger than 7–8 a₀ (see Fig. 4) is clearly apparent at ρ = 10.2 a₀. The ground state and the first excited state present an avoided crossing (Fig. 4) at configurations where two H atoms have an internuclear distance of 2.5 a₀ and the third H atom is far away. To better understand the stereographic projections of Fig. 3, the D_3h, C_2v and D_∞h geometries are included in the figure using balls and sticks representation. In the first column (ρ = 2.18 a₀) the global minimum of the ground electronic state, with a D_3h geometry appears in the centre of the plot. For this configuration the upper states are much higher in energy and present an avoided crossing. The avoided crossing between the two excited states appears for all the D_3h configurations, as it has been mentioned before (see Fig. 3). For longer values of hyperradius (second and third columns) the two excited surfaces approach to the ground state, and within each other. In Fig. 4 are represented several one-dimensional cuts of the ground and first excited singlet states showing the approaching of the surfaces using Jacobi coordinates. In this figure are plotted several values of R Jacobi coordinate showing that the avoided crossing appears for R > 7 a₀.

Fig. 3 Stereographic projections (β = sin(θ)cos(φ − 30) and γ = sin(θ) sin(φ − 30)) of hyperspherical coordinates⁴⁰ of the potential energy of the three singlet states for three constant values of ρ. The top/middle/bottom contourplots correspond to the 1¹A′/2¹A′/3¹A′ states, respectively.

Fig. 4 Avoided crossing between the ground and first excited electronic singlet states. The representation uses the two distances Jacobi coordinates with a fixed value of the angle in 90 degrees.

The full dimensional NACMEs, obtained analytically using the Hellmann–Feynman theorem,²⁸ are specially important at conical intersections (where they diverge) and at avoided crossings. The potential energy curves and NACMEs are plotted in Fig. 5 together with the ab initio points as a function of the hyperspherical angle φ at a fixed hyperradius ρ = 10 Bohr for three different values of hyperspherical angle θ, specified in every plot. The transformation between bond coordinates and hyperspherical coordinates assumed in this work is the one given by Billing and coworkers,⁴⁰ and the derivatives and NACMEs in any of the hyperspherical coordinates can be obtained using the chain rule.

Fig. 5 Potential energy curves (left column) and NACMEs (right column) of the three singlet states as a function of the hyperspherical⁴⁰ angle coordinate ϕ for a fixed value of ρ = 10 a.u. The value of θ angle is specified in each plot. Points correspond to ab initio points and the line is the analytical fit and/or NACME. The zero of energy corresponds to H + H + H⁺.

The zero of energies is displaced so that the ground state is zero for θ = 30^o and φ = 0^o. The agreement between the ab initio points atnd the fitting is very good, even when energies above −0.9E_H (0.1E_H in the figure) are not included in the fitting. The right column representing the NACMEs shows as well a good agreement. The position and magnitude agree well and the main disagreements, in the NACMEs between the two upper states, represented by blue empty circles and blue lines, occur at the cusps of the NACMEs where the analytical functions slightly underestimate the value. These disagreements could be due to several things, the NACMEs calculated with MOLPRO package are obtained by means of finite differences method, and it has problems at the divergence occurring at conical intersections. On the other side the fitting procedure, although is expected to be correct at these points thanks to the DIM approximations, is as well very sensitive because the surfaces present changes that can be challenging for the optimisation of the three body terms.

3.2 Analytical fit of triplet states

The fit of the triplet states includes ab initio points up to values above the dissociation energy, −0.98E_H. Higher energies are not included, as in the case of the singlet states, because all stationary points lie below the energy chosen, in addition the upper sheet is repulsive what increases the error of the lower sheets when higher energies are considered. The order of the expansion of the three body terms (eqn (11)) is M = 6. In Fig. 6 are plotted the residuals as a function of the energy, where the zero of energies is kept at the minimum of the ground (singlet) electronic state. The total rms deviation (12 cm⁻¹) is showed using the horizontal black lines and the global dissociation with the vertical line. The errors obtained for the lowest triplet state remain lower than 50 cm⁻¹ and the errors for the first triplet state increase up to maximum errors of 200 cm⁻¹ for energies above the complete dissociation. For a detailed description of the long range interaction in the lowest Triplet state see section B in the ESI.†

Fig. 6 Residuals (algebraic deviations) of the three triplet PESs vs. the energy (in cm⁻¹). The black dashed vertical line specifies the energy of the complete dissociation (H + H + H⁺).

In Table 4 we present the rms as a function of the Jacobi coordinate R for the lowest sheet, to show the error in the entrance channel (H⁺₂ + H). The rms for the lowest triplet state remains lower than 8 cm⁻¹ for all intervals considered and lower than 0.5 cm⁻¹ for the long-range region. Only the two lowest sheets are included in the table since all the points of the third state lie out of the maximum energy considered in the fit. Since the two upper singlet state dissociate directly into the complete dissociation, no rms for intervals is included for the second triplet state, the rms of the all points included in the fit is 25 cm⁻¹.

Table 4 RMS obtained for several intervals of energies with the DIM approximation plus the long-range terms and the global three-state fitting. The rms for the total number of points in the fit includes points up to −0.98E_H

	DIM + LR	3-State fit
	rms/cm^−1	rms/cm^−1
Distances/a0
3.0 ≤ R ≤ 7.0	2893.2	3.5
8.0 ≤ R ≤ 15.0	4.4	6.7
15.0 ≤ R ≤ ∞	0.02	0.4
Energies (1^3A′)
Minimum to H^+2 + H	3563.0	4.7
H^+2 + H to H + H + H^+	3201.9	7.4
All points	3445.5	7.1
Energies (2^3A′)
All points	459.5	25.1

---

A list of the stationary points of the two lowest triplet states is given in Table 5. The third triplet states presents no minima or stationary points for energies below the global dissociation. In the table are included the CBS extrapolated values, the multi-sheet fit values of this work and the ones obtained from literature.¹⁶ The results published by Varandas and coworkers for the two lowest triplet states are based on the single-valued double many-body expansion (DMBE) theory, applied to a grid of points calculated with MOLPRO package using the cc-pV5Z basis set. The lowest triplet state presents three equivalent minima at D_∞h configurations 52528.18 cm⁻¹ above the minimum of the ground singlet state (2957.2 cm⁻¹ below the dissociation channel H⁺₂ + H), 10 cm⁻¹ lower than the one from literature. The saddle points between these minima, of C_2v symmetry, appear 2550.3 cm⁻¹ above the linear minimum, 64 cm⁻¹ lower than the one provided by Varandas et al. The dissociation into the H⁺₂ + H is placed 2957 cm⁻¹ above the minimum, only 2 cm⁻¹ lower than the one from literature.

Table 5 Energy comparison of the stationary points present in the two triplet electronic states. Distances are given in Bohr and energies in Hartrees

r 12	r 13	r 23	CBS	Fit	Literature^a
a Varandas et al.^16
First triplet state (1^3A′)
Global minimum (D∞h)
2.456	2.456	4.912	−1.116121	−1.116103	−1.116061
Saddle points (C2v)
1.99	5.286	5.286	−1.104501	−1.104570	−1.104234
H^+2 + H channel
2.000	∞	∞	−1.102647	−1.102647	−1.102598
Avoided crossing between 1st and 2nd states
3.654	3.654	3.654	−1.034843	−1.034842	−1.034590
First triplet excited State (2^3A′)
Global minimum (D3h)
3.654	3.654	3.654	−1.034843	−1.034842	−1.034590
Global dissociation (H + H + H^+)
∞	∞	∞	−1.000008	−1.000008	—

---

The minimum of the upper sheet, of D_3h symmetry, is produced by the conical intersection between the first and second triplet states. Both states approach to each other creating a conical intersection for all equilateral triangle configurations. The global minimum of the second triplet state appears in this type of configuration, 14881 cm⁻¹ above the H⁺₂ + H channel. This state dissociates directly into the global dissociation channel (H + H + H⁺), the channel is 7645.40 cm⁻¹ above the minimum. A global idea of the energetic of the full system, including singlet and triplet states, can be seen in Fig. 1.

The stationary points of Table 5 are shown graphically in Fig. 9 using stereographic projections of hyperspherical coordinates. The zero of energies is placed in the global dissociation (H + H + H⁺). The top row presents a three dimensional representation of the three states, which are plotted below using the two-dimensional stereographic representations. The three columns represent three cuts of the hyperradius ρ of the potential energy surface. The first column, ρ = 4.53 a.u. corresponds to the minimum of the lowest triplet state, appearing in this plots in the vertices of the imaginary triangle (see Section 3.1). The saddle points joining the minima appear in the sides of the triangle. The value of ρ in the minimum of the second triplet state is 4.80 a.u., no apparent differences show with respect to the value of ρ for the minimum of the lowest triplet state, this is the reason why not a separate plot is added. The minimum of the second sheet can be easily seen in the three dimensional plots, where the two opposite cones form the conical intersection producing the global minimum of the second triplet state. In the two dimensional representations this minimum appears in the middle of the plot in the second row. The other values of ρ included in the figure show global features of the three states as the approaching of the surfaces with the increasing of the ρ value. However, there are no points in the third triplet state for any value of ρ that lie below the global dissociation. The first and second triplets form the double cone creating the conical intersection that can be seen for all values of ρ.

Considering two equivalent H–H distances, we have performed an optimisation of these distances for angles between (α) two vectors running from 0 to 180 degrees at CI level using the aug-cc-pVTZ basis set of Dunning.⁴¹ With the optimised values we calculate the extrapolated values and called the global fit, the results are plotted in Fig. 7. The global minimum of the lowest triplet state appears for a value of α = 180° (D_∞h), the conical intersection and minimum of the second triplet state at α = 60° (D_3h) and the saddle point of the lowest singlet state at α = 21.5° (C_2v). All features are well described by the global fit and the discrepancies obtained for the second and triplet states are due to the maximum energy value considered for the points included in the fit, −0.98E_H (0.02E_H in the figure, since the zero of energy is placed in the global dissociation). Nevertheless, the fit provides very good agreement for much higher energies, including the third triplet state which no ab initio points have been included in the fit.

Fig. 7 Cuts of the potential energy surface of the three triplet states, as a function of the angle between two equivalent relaxed H–H distances. Solid lines correspond to fitted values and dots to ab initio points. The zero of energíes is placed in the global dissociation H + H + H⁺.

The potential curves and NACMEs as a function of the hyperspherical angle ϕ, for a fixed value of the hyperradius ρ = 10 Bohr and three values of the hyperspherical angle θ are shown in Fig. 8. The value of ρ was chosen so that closeness of the surfaces would show larger values of NACMEs. For values of ρ near the conical intersection only the NACMEs in the intersection was appreciable, due to degeneracy lifting stated by Jahn–Teller theorem.⁴² The agreement of the fit in the potential energy curves and the NACMEs is very good and only some disagreements are obtained for the NACMEs between the second and third states for θ = 70° for the values ϕ = 60,180 and 240°. The disagreements are a consequence of the error in the fit due to the non included points for higher energies. The three peaks appearing in the ab initio NACME points, not reproduced by the analytical fit, are between values of both states higher than 22 000 cm⁻¹ above the global dissociation.

Fig. 8 Potential energy curves (left column) and NACMEs (right column) of the three triplet states as a function of the hyperspherical⁴⁰ angle coordinate ϕ for a fixed value of ρ = 10 a.u. The value of θ angle is specified in each plot. Points correspond to ab initio points and the line is the analytical fit and/or NACME.

Fig. 9 Stereographic projections (β = sin(θ) cos (φ − 30) and γ = sin(θ) sin(φ − 30)) of hyperspherical coordinates⁴⁰ of the potential energy of the three triple states for three constant values of ρ. In the top row there are the three-dimensional curves of the three states. The top/middle/bottom contourplots correspond to the 1³A′/2³A′/3 ³A′ states, respectively.

3.3 Rovibrational analysis

Rovibrational bound states on the adiabatic PES are calculated using hyperspherical coordinates with the same methods previously described and used for the singlet^24,26 and triplet states.¹¹ For the details with respect to the hyperspherical coordinates, symmetry and Hamiltonian/wavefunction representations we refer the reader to these publications.

3.3.1 Ground singlet state. We have calculated the rovibrational states of the ground singlet state. The minimum of the second singlet state does not support vibrational states below the H⁺₂ + H channel, and the third singlet state can be considered as repulsive. The rovibrational eigenvalues for J = 0 are listed in Table 5 and compared with the experimental⁴³ values available in the literature and theoretical calculations.^7,18 The classification of the eigenstates is made using the approximated quantum numbers v₁, v₂ and l, which correspond to the symmetric stretch quantum number, the degenerate vibration (a superposition of bending and antisymmetric stretch) and the vibrational angular momentum, respectively. The quantum numbers given in the table coincide for previous calculations²⁶ for the ground Singlet state of H⁺₃ and the experimental⁴³ and calculated values,⁷ all given in Table 6. The rovibrational energy levels listed in Table 5 include states up to energies above the barrier to linearity. The comparison with the other multi-sheet¹⁸ fit shows a very good agreement with differences lower than 3 cm⁻¹ for each symmetry. The agreement with the values published by Pavanello et al. remains lower than 4 cm⁻¹ for most of the states for all symmetries of the wavefunction with the exception of a few states with discrepancies higher than 5 cm⁻¹.

Table 6 Rovibrational eigenvalues for J = 0 of the ground singlet electronic state for the present multi-sheet fit (with zero-point energy, ZPE = 4362.1197 cm⁻¹), some experimental values,⁴³ for the multi-sheet fit of Viegas et al.¹⁸ (with ZPE = 4361.59 cm⁻¹) and the very accurate potential energy surface of the ground state of Pavanello et al.⁷

(v1,v^l2)	Γ	This work	Exp.^43	MSGPES^18	GPES^7
(1,0^0)	A 1′	3180.96		3179.32	3178.29
(0,2^0)		4777.10		4777.87	4778.15
(2,0^0)		6265.03		6265.09	6262.13
(0,3^3)		7283.32		7286.25	7285.56
(1,2^0)		7770.10		7770.79	7769.23
(0,4^0)		8999.52		9002.03	9001.57
(3,0^0)		9252.86		9256.88	9251.91
(1,3^3)		9967.41		9970.31	9968.94
(2,2^0)		10593.04		10596.63	10593.19
(0,5^3)		10920.72		10922.82	10923.36
(1,4^0)		11813.80		11816.95	11814.52
(4,0^0)		12144.28			12146.46
(0,6^0)		12377.37			12382.15
(2,3^3)		12587.53			12590.53
(3,2^0)		13286.84			13288.91
(1,5^3)		13401.70			13405.25
(0,7^3)		13719.95			13725.43
(2,4^0)		14195.50			14198.51
(1,6^0)		14663.36			14666.18
(0,8^0)		14904.23			14909.87
(4,0^0)		14934.84			14940.15
(0,6^6)		15075.84			15080.18
(3,3^3)		15163.65			15168.30
(0,3^3)	A 2′	7491.26	7492.91	7493.86	7492.78
(1,3^3)		10209.79		10212.20	10210.33
(0,5^3)		11527.20		11530.21	11529.24
(2,3^3)		12831.24			12832.17
(1,5^3)		13755.51			13756.61
(0,7^3)		14564.54			14566.28
(0,6^6)		15188.77			15190.71
(0,1^1)	E′	2520.88	2521.41	2520.84	2521.30
(0,2^2)		4996.77	4998.05	4997.43	4997.89
(1,1^1)		5555.78	5554.06	5555.45	5554.20
(0,3^1)		7004.33		7006.18	7005.97
(1,2^2)		7870.39		7871.61	7870.23
(2,1^1)		8489.09		8491.08	8488.01
(0,4^2)		9110.73		9113.54	9113.04
(1,3^1)		9651.87		9654.96	9653.70
(0,4^4)		9996.90		9999.10	9997.18
(2,2^2)		10644.76	10645.38	10648.17	10645.31
(0,5^1)		10859.73	10862.90	10862.46	10862.75
(3,1^1)		11322.17	11323.10	11327.69	11323.12
(0,5^5)		11656.08	11658.40	11659.89	11658.31
(1,4^2)		12077.64		12081.86	12079.40
(2,3^1)		12300.93	12303.36		12303.33
(0,6^2)		12473.80	12477.38		12477.39
(1,4^4)		12696.24			12697.40
(3,2^2)		13315.70			13318.19
(1,5^1)		13392.26			13395.2
(0,6^4)		13589.82			13592.25
(0,7^1)		13696.81	13702.37		13702.58
(4,1^1)		14051.62			14055.01
(1,5^5)		14215.51			14218.02
(2,4^2)		14475.26			14478.22
(1,6^2)		14888.78			14890.47
(3,3^1)		14896.84			14901.82
(0,8^2)		15117.62	15122.80		15122.64

---

3.3.2 Triplet states. In the previous sections we mentioned the conical intersection line that is formed for D_3h symmetries between the two lowest triplet states, of special importance for the second triplet states which global minimum is produced by this intersection. The more recent studies, including rovibrational eigenstates in the triplet states of H⁺₃, try to account for this degeneracy and the geometrical phase produced by the intersection. The rovibrational levels presented in this section are calculated in the two lowest diagonalised adiabatic states, since only those states possess minima that support rovibrational levels.

In Table 7 are presented the rovibrational states of the first triplet state (a ³A′) for the present fit and the latest found in literature.¹⁷ The potential representation used by Alijah et al. is made using the diabatization procedure proposed by Longuet–Higgins⁴⁴ to transform two adiabatic potential energy surfaces to a 2 × 2 diabatic matrix. In their paper they present the eigenstates considering either 1- or 2-states for the calculations. Because in our calculations we use the diagonalised electronic states we compare our results with the one-state calculation given in their work.¹⁷ The quantum numbers used for the assignment of the vibrational levels use the same procedure that the one proposed by Watson.⁴⁵ In hyperspherical coordinates, v₁ and v₃ corresponds to symmetric and antisymmetric stretching vibrations, and v^l₂ corresponds to the bending vibration, where l is the associated vibrational angular momentum. Considering all the above the agreement between both calculations is very good with errors lower than 4 cm⁻¹ for the states below the dissociation channel (H⁺₂ + H, 2957.20 cm⁻¹ above the minimum), discrepancy increases for higher vibrational states. These disagreements are attributed to the description of the long-range region of the potential energy surface. In Table 4 where the rms of the fit was presented for for several ranges of distances (i.e. long-range) and energies, the rms remains lower than 8 cm⁻¹ for all points in the lowest sheet. The effect of the 2-state calculation, studied by Alijah and coworkers, never lead to differences with the 1-state calculation larger than 2–3 cm⁻¹. The comparison between some previous calculations of the same authors¹⁶ using the same potential energy surface for the lowest sheet, shows that the geometrical phase does not produce changes in the eigenvalues larger than 3–4 cm⁻¹, except in some of the higher vibrational states, very near the conical intersection. For all this, we conclude that the difference between the present work and the results of Alijah et al. might arise from the analytical fit and particularly from the long-range description, of very high accuracy in the present work.

Table 7 Rovibrational eigenvalues for J = 0 of the lowest triplet electronic state for the present multi-sheet fit (with zero-point energy, ZPE = 1720.14 cm⁻¹) and one state calculation of Alijah et al.¹⁷ (with ZPE = 1722.043 cm⁻¹). Note that the bending mode (v₂) has vibrational symmetry A₁, vibronic symmetry A₂ (³A′ electronic state is of symmetry A₂), what means that v₂ = 1 implies |l| = 1, therefore v₂ = 1 does not exist for J = 0⁴⁵

(v1,v^l2,v3)	i	A 1′	i	A 2′	i	E′	A 1′ ^a	i	A 2′ ^a	i	E′ ^a
a Bound states computed with a two-surface potential calculated by Alijah et al.^17
(0,0^0,0)	0	0.0			0	0.0	0.0			0	0.0
(0,0^0,1)			0	738.93	1	736.93		0	738.49	1	738.49
(1,0^0,0)	1	974.21			2	974.21	975.05			2	975.06
(0,2^0,0)	2	1269.41			3	1269.48	1273.73			3	1273.79
(1,0^0,1)			1	1467.06	4	1466.96		1	1474.57	4	1474.4
(0,0^0,2)	3	1566.29			5	1566.42	1573.70			5	1573.8
(0,2^0,1)			2	1720.43	6	1718.21		2	1730.39	6	1728.4
(2,0^0,0)	4	1912.83			7	1913.33	1922.55			7	1923.1
(1,2^0,0)	5	1922.22			8	1934.50	1940.34	3		8	1951.0
…			3	1956.74	9	1954.26		4	1972.87	9	1970.7
…	6	2145.24			10	2111.48	2158.72			10	2136.9
…	7	2165.66			11	2152.99	2188.49			11	2166.4
…			4	2183.87	12	2238.51		5	2205.01	12	2251.0
…			5	2259.96	13	2246.28		6	2271.13	13	2259.3
…	8	2313.36			14	2337.17	2308.61			14	2334.9
…	9	2330.54			15	2344.14	2340.09			15	2357.0

---

4 Conclusions

In the present work we have obtained an analytical representation of the potential energy surfaces, derivatives and non-adiabatic coupling matrix elements of the three lowest singlet and triplet states of the system of H⁺₃. New terms have been included to improve the long range interaction in the diatomic systems, and both dissociation channels, H⁺ + H₂ and H + H⁺₂. Those terms are particularly important to properly describe the collision dynamics at low and ultracold temperatures. In addition, we have shown that the present fits also describe very accurately the bound rovibrational states of H⁺₃, in singlet and triplet states. The NACME obtained in the present work are also very accurate and reliable. For all these reasons we consider these fits as the most accurate nowadays available in the literature to describe the non-adiabatic and low temperature collision dynamics for this system, of particular relevance in astrochemistry.

Conflicts of interest

There are no conflicts to declare.

References

1. T. Oka, Philos. Trans. R. Soc., A, 2012, 113, 8738–8761.
2. T. Oka, Chem. Rev., 2013, 370, 4991–5000.
3. J. Pelley, ACS Cent. Sci., 2019, 5, 741–744.
4. B. J. McCall, Philos. Trans. R. Soc. London, Ser. A, 2000, 358, 2385–2401.
5. M. Pavanello, W.-C. Tung, F. Leonarski and L. Adamowicz, J. Chem. Phys., 2009, 130, 074105.
6. O. L. Polyansky, A. Alijah, N. F. Zobov, I. I. Mizus, R. I. Ovsyannikov, J. Tennyson, L. Lodi, T. Szidarovszky and A. G. Császár, Philos. Trans. R. Soc., A, 2012, 370, 5014–5027.
7. M. Pavanello, L. Adamowicz, A. Alijah, N. F. Zobov, I. I. Mizus, O. L. Polyansky, J. Tennyson, T. Szidarovszky and A. G. Császár, J. Chem. Phys., 2012, 136, 184303.
8. I. I. Mizus, O. L. Polyansky, L. K. McKemmish, J. Tennyson, A. Alijah and N. F. Zobov, Mol. Phys., 2019, 117, 1663–1672.
9. S. M. Bijit Mukherjee and S. Adhikari, J. Phys.: Conf. Ser., 2016, 7, 012050.
10. L. Velilla, B. Lepetit, A. Aguado, J. A. Beswick and M. Paniagua, J. Chem. Phys., 2008, 129, 084307.
11. C. Sanz-Sanz, O. Roncero, C. Tablero, A. Aguado and M. Paniagua, J. Chem. Phys., 2001, 114, 2182.
12. M. Cernei, A. Alijah and A. J. C. Varandas, J. Chem. Phys., 2003, 118, 2637.
13. T. M. Ferreira, A. Alijah and A. J. C. Varandas, J. Chem. Phys., 2008, 128, 054301.
14. L. P. Viegas, M. Cernei, A. Alijah and A. J. C. Varandas, J. Chem. Phys., 2004, 120, 253–259.
15. L. P. Viegas, A. Alijah and A. J. C. Varandas, J. Phys. Chem. A, 2005, 109, 3307–3310.
16. A. J. Varandas, A. Alijah and M. Cernei, Chem. Phys., 2005, 308, 285–295.
17. A. Alijah and V. Kokoouline, Chem. Phys., 2015, 460, 43.
18. L. P. Viegas, A. Alijah and A. J. C. Varandas, J. Chem. Phys., 2007, 126, 074309.
19. C. Sanz-Sanz, O. Roncero, M. Paniagua and A. Aguado, J. Chem. Phys., 2013, 139, 184302.
20. A. Aguado, P. Barragán, R. Prosmiti, G. Delgado-Barrio, P. Villarreal and O. Roncero, J. Chem. Phys., 2010, 133, 024306.
21. D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys., 1993, 99, 1914.
22. K. A. Peterson, R. A. Kendall and T. H. Dunning, J. Chem. Phys., 1993, 99, 1930–1944.
23. K. A. Peterson, R. A. Kendall and T. H. Dunning, Jr., J. Chem. Phys., 1993, 99, 9790.
24. L. Velilla, M. Paniagua and A. Aguado, Int. J. Quantum Chem., 2010, 111, 387–399.
25. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schützet al., MOLPRO, version 2015.1, a package of ab initio programs, 2015, see http://www.molpro.net.
26. A. Aguado, O. Roncero, C. Tablero, C. Sanz-Sanz and M. Paniagua, J. Chem. Phys., 2000, 112, 1240.
27. S. Ghosh, S. Mukherjee, B. Mukherjee, S. Mandal, R. Sharma, P. Chaudhury and S. Adhikari, J. Chem. Phys., 2017, 147, 074105.
28. C. Sanz-Sanz, A. Aguado, O. Roncero and F. Naumkin, J. Chem. Phys., 2015, 139, 184302.
29. J. O. Hirschfelder and W. J. Meath, The Nature of Intermolecular Forces, Wiley, Chichester, 1967, vol. 12.
30. C. A. Coulson, Proc. - R. Soc. Edinburgh, Sect. A: Math. Phys. Sci., 1941, 161, 20.
31. P. Gniewek and B. Jeziorski, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 90, 022506.
32. P. Gniewek and B. Jeziorski, Mol. Phys., 2016, 14, 1176.
33. T. Holstein, J. Phys. Chem., 1952, 56, 832.
34. C. Herring, Rev. Mod. Phys., 1962, 34, 631.
35. A. Aguado, C. Tablero and M. Paniagua, Comput. Phys. Commun., 1998, 108, 259.
36. R. Röhse, W. Kutzelnigg, R. Jaquet and W. Klopper, J. Chem. Phys., 1994, 101, 2231–2243.
37. L. Wolniewicz, J. Chem. Phys., 1993, 99, 1851–1868.
38. W. Cencek, J. Rychlewski, R. Jaquet and W. Kutzelnigg, J. Chem. Phys., 1998, 108, 2831.
39. O. L. Polyansky, R. Prosmiti, W. Klopper and J. Tennyson, Mol. Phys., 2000, 98, 261–273.
40. G. D. Billing and J. T. Muckerman, J. Chem. Phys., 1989, 91, 6830–6838.
41. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007.
42. H. A. Jahn, E. Teller and F. G. Donnan, Proc. R. Soc. London, Ser. A, 1937, 161, 220–235.
43. T. Furtenbacher, T. Szidarovsky, E. Matyus, C. Fabri and A. G. Csaszar, J. Chem. Theory Comput., 2013, 9, 5471–5478.
44. H. C. Longuet-Higgins, Adv. Spectrosc., 1961, 2, 429–472.
45. J. K. Watson, Mol. Phys., 1970, 19, 465–487.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0cp04100a

---