## Abstract

After a short historical account of the theory of the ion, two *ab initio* methods are reviewed that allow the computation of the ground-state potential energy surface (PES) of in the Born–Oppenheimer (BO) approximation, with microhartree or even sub-microhartree accuracy, namely the R12 method and the method of explicitly correlated Gaussians. The BO-PES is improved by the inclusion of relativistic effects and adiabatic corrections. It is discussed how non-adiabatic effects on rotation and vibration can be simulated by corrections to the moving nuclear masses. The importance of the appropriate analytic fit to the computed points of the PES for the subsequent computation of the rovibronic spectrum is addressed. Some recent extensions of the computed PES in the energy region above the barrier to linearity are reviewed. This involves a large set of input geometries and the correct treatment of the dissociation asymptotics, including the coupling with the first excited singlet state. Some comments on this state as well as on the lowest triplet state of are made. The paper ends with a few remarks on the ion .

## 1. Introduction

, consisting of three nuclei and two electrons, is the simplest *polyatomic* molecule. While very accurate quantum mechanical calculations on 2-electron *atoms* (the He isoelectronic series) and 2-electron *diatomic* molecules (H_{2} or HeH^{+}) were published already in the pioneer days of quantum mechanics (Hylleraas 1929; James & Coolidge 1933), has resisted an accurate quantum mechanical study much longer. In table 1, the progress of accuracy as a function of time is documented. One clearly distinguishes three classes of methods: (i) self-consistent field (SCF), (ii) configuration interaction (CI)—including that in terms of *pair natural orbitals* (PNO), and (iii) explicitly correlated methods such as Gaussian geminals (GG) or R12, to be explained later in this section.

Nowadays, the potential energy surface (PES) of is available with an error of about 0.1 *μE*_{h} or 0.02 cm^{−1} (*submicrohartree accuracy*) (Cencek *et al*. 1998), which is much better than what is usually called *spectroscopic accuracy*, by which one means errors of the order of 1 cm^{−1}.

In order to achieve this high accuracy, very sophisticated quantum chemical methods are required. One must take care of *electron correlation* at a high level, namely using an ansatz for the electronic wave function, that depends explicitly on the interelectronic distance *r*_{12}. There are essentially two approaches available, namely

the R12-method (Kutzelnigg 1985; Kutzelnigg & Klopper 1991; Klopper & Noga 2003), which is a configuration-interaction (CI) method, augmented by terms linear in

*r*_{12}, andthe method of exponentially correlated Gaussians (GG; Rychlewski

*et al*. 1999; Bukowski*et al*. 2003; Rychlewski & Komasa 2003), in which the basis consists of Gaussians depending on , and .

One must also take care of relativistic effects. For light atoms such as H these are very small, so that the lowest order of perturbation theory (Kutzelnigg 1989) is actually sufficient.

The experimental information available for is mainly in its infrared (IR) spectrum (Oka 1980, 2000), such that theory is challenged to reproduce and interpret the vibration–rotation term levels. On the lowest non-trivial theoretical level (i.e. beyond the harmonic approximation) this requires the knowledge of the PES in a sufficiently large neighbourhood of the equilibrium structure, from which the rovibrational term levels and the corresponding nuclear wave functions are accessible in the Born–Oppenheimer (BO) approximation. The first step to go beyond this approximation consists in the evaluation of the *adiabatic corrections*. To achieve full agreement between theory and experiment, and to make unambiguous assignments, *non-adiabatic corrections* are necessary as well. These are much more difficult, and represent at present the most serious challenge, if one strives for very high accuracy.

There is increasing interest in understanding what happens if one approaches energies close to the linear structure of (Gottfried *et al*. 2003) and the avoided crossing with the excited surface corresponding to the dissociation channel +H (Jaquet 2003).

## 2. Previous calculations on

The simplest model to describe is the Hückel theory (Hückel 1932; Kutzelnigg 1978). The Hückel matrices for in a cyclic (equilateral triangular) and an open (symmetric linear) structure are respectively(2.1)with the eigenvalues (*α*+2*β*, *α*−*β*, *α*−*β*) and . For the binding energy (with respect to H^{+}+2 H, for which the Hückel energy is 2*α*) one gets for the two structures (with the lowest molecular orbital (MO) doubly occupied)(2.2)with the cyclic structure more strongly bound. Comparing this with the Hückel binding energy 2*β* for H_{2}, one finds that an extra proton is bound to H_{2} with about the same energy as the two H-atoms are bound in H_{2}. Both this prediction and that of a cyclic structure for were confirmed by more sophisticated calculations and experiment. The Hückel theory even accounts for the bond lengths in the two structures of (Kutzelnigg 1978) via a bond-length/bond-order relation.

The Hückel theory also predicts a stable excited state of with the two lowest MOs each singly occupied, and with the energy minimum in the symmetric linear geometry corresponding to the electronic configuration 1*σ*_{g}1*σ*_{u}, a binding energy of with respect to +H and a bond length of approximately 2.5*a*_{0}, again using the bond-length/bond-order relation (Kutzelnigg 1978). This state is higher in energy than H_{2} +H^{+}. In the Hückel theory, only *valence-excited states* and no Rydberg states are possible. The Hückel theory does not distinguish between a singlet and a triplet state, but for the triplet only the stability with respect to the lowest *triplet dissociation limit* (+H) matters.

The first *ab initio* calculations on the ground state were performed by Coulson (1935) in terms of MO theory and Hirschfelder *et al*. (1936) and Hirschfelder (1938) at valence-bond level with Slater type orbital (STO) basis sets. Coulson assumed an equilateral triangular structure and got a bond length of 0.85 Å, in good agreement with the now accepted value of 0.873 Å, while Hirschfelder did not impose a *D*_{3h} structure and predicted a distorted triangular equilibrium structure. Calculations striving at quantitatively acceptable results started around 1964 using Gaussian basis sets.

That the progress with accurate calculations on was much slower than for the isoelectronic systems He or H_{2} is related to the fact that special coordinate systems that simplified the 1- or 2-centre problem (Hylleraas 1929; James & Coolidge 1933) have not been available for a 3-centre system.

In an early study of both the triangular and the linear structure of , including electron correlation (Kutzelnigg *et al*. 1967), the significant geometry dependence of the correlation energy was pointed out (−0.039*E*_{H} for the triangular and −0.047*E*_{H} for the linear structure), which made pure SCF calculations obsolete. is the prototype of a 3-centre-2-electron bond, and it is certainly noteworthy that this bond is about twice as strong as the corresponding 2-centre-2-electron bond.

H_{2} has a substantial *proton affinity* (418 kJ mol^{−1}) and must hence be classified as a *base*. However, compared to other *basic* molecules in the gas phase, e.g. H_{2}O (670 kJ mol^{−1}) or ions, H_{2} is an extremely weak base. Correspondingly, is a very strong acid. By far, much higher proton affinities than for isolated molecules are found for hydrogen-bonded liquids, such as water.

The study of the entire potential surface of in view of an understanding of the rovibrational spectrum started probably with Carney & Porter (1974).

An important step in advance was a paper by Meyer *et al*. (1986). It was based on traditional CI calculations (*not* including *r*_{12}-dependent terms).

A special feature of this paper has been a very elegant selection of the points of the PES to be computed that was later followed by other authors (Röhse *et al*. 1994; Jaquet *et al*. 1998).

Meyer *et al*. have chosen 69 points of the PES, labelled by three integers *n*_{a}, *n*_{x} and *n*_{y}, which are related to the Morse-type symmetry adapted deformation coordinates (MSADC) *S*_{a}, *S*_{x} and *S*_{y}(2.3)(2.4)(2.5)with(2.6)

These 69 points define a part of the whole PES from the minimum up to well above the barrier to linearity. The functional form is chosen such that it can be used up to the asymptotic region (H^{+}+H_{2}, H+H+H^{+}), but the most reliable part includes the energy region up to approximately 20 000 cm^{−1} above the minimum.

The choice *β*=1.3 in equation (2.6) is recommended for a region close to the minimum, while *β*=1.0 has to be used in the asymptotic region close to dissociation.

The essential progress in accuracy came with the use of explicitly correlated wave functions of the R12 and GG type to be described subsequently.

A study by Anderson (1992) deserves to be mentioned, since it has demonstrated that the Monte-Carlo (MC) method is able to achieve up to microhartree accuracy for systems like . Anderson (1992) has only published calculations for the equilibrium geometry of . MC studies never played a role for the PES of .

## 3. The R12-method

Although the R12-method is applicable to systems with an arbitrary number of electrons (Klopper & Noga 2003), here, it is sufficient to formulate it for a 2-electron system (Röhse *et al*. 1993). The quantum chemical standard method to include electron correlation is configuration interaction (CI). For a 2-electron system, one can factor out the spin, and for a singlet state, the spin-free CI function is of the simple form(3.1)where this expansion converges slowly. Actually, if *L* is the highest angular quantum number *l* of the *ϕ*_{p} included in the expansion (3.1), then the error goes as (*L*+1)^{−3}. This is owing to the fact that the exact 2-electron wave function has a *correlation cusp* at the point of coalescence of two electrons (Kato 1957; Kutzelnigg 2003). There, the wave function is linear in *r*_{12}, which cannot be described by a CI wave function. A considerable speed-up of the convergence is obtained if one replaces equation (3.1) by(3.2)where *Φ*(1, 2) is a simple product-type *reference function*, e.g. an eigenfunction of the bare-nuclear Hamiltonian, i.e. the Hamiltonian for two non-interacting electrons in the field of two protons. This is actually the most convenient choice and preferable to a Hartree–Fock reference state. Now the error goes ideally as (*L*+1)^{−7}. The implementation of this ansatz requires in addition to the traditional electron repulsion integrals(3.3)also the integrals (with *T* the operator of the kinetic energy)(3.4)which are not really difficult. The details are found elsewhere (Klopper & Noga 2003). The improved convergence of CI-R12 with respect to conventional CI is documented in table 2.

With a (30s, 20p, 12d, 9f) basis a ground-state energy *E*=−1.343835 *E*_{h} at the (equilateral triangular) equilibrium geometry has been obtained, with microhartree accuracy (Röhse *et al*. 1993). Omitting the *f*-functions in the basis, one gets errors of the PSE of a few microhartrees. The size of the s, p, d is not so critical. The published PES (Röhse *et al*. 1994) was performed with a (16s, 10p, 8d) basis contracted to (10s, 8p, 6d), but for some points additional calculations with a contracted (10s, 8p, 6d, 4f) basis were performed. More recent calculations on an extended region of the PES (Polyansky *et al*. 2000; Jaquet 2002, 2003) were done with the latter basis.

A generalization of the ansatz (3.2) that contains more than one reference function multiplied by *r*_{12} (Gdanitz & Röhse 1995) has been used in order to account for excited singlet states together with the ground state.

## 4. The method of exponentially correlated Gaussians

In the GG method (Rychlewski *et al*. 1999) for a 2-electron system in a singlet state, one makes the following ansatz for the spin-free electronic wave function (Cencek *et al*. 1998):(4.1)

(4.2)

All nonlinear parameters *α*_{1i}, *α*_{2i} and *β*_{i}, as well as the reference positions *A*_{i} and *B*_{i} are optimized consecutively in various cycles. The optimization is stopped, when the gain in one cycle was less than 10^{−8} Hartree. The expansion length was *K*=1300. Again we must refer to the original reference (Cencek *et al*. 1998) for details. The energy obtained at the equilibrium geometry *E*=−1.343835518 *E*_{h} was compared with the estimate for the exact value *E*=−1.343835624 *E*_{h} (Cencek *et al*. 1995), i.e. the error is of the order 0.1 *μE*_{h} or 0.02 cm^{−1}.

Unlike with the R12 method, one does, with GGs, not describe the correlation cusp, i.e. the linear dependence of the wave function on *r*_{12} for *r*_{12}→0 correctly. It is therefore not immediately understandable why one gets an excellent convergence behaviour. The key to this is obviously related to the satisfactory convergence of the expansion of one-electron-functions in a Gaussian basis (Kutzelnigg 1994).

Precursors of a refined application of the GG method were studies by Salmon & Poshusta (1973) and Alexander *et al*. (1990). Owing to a less sophisticated optimization of the nonlinear parameters, the results were less accurate than those of Cencek *et al*. (1995, 1998).

## 5. Relativistic corrections

Relativistic corrections have been evaluated by means of direct perturbation theory (DPT) (Kutzelnigg 1989; Cencek 2003), to be more precise, by stationary direct perturbation theory (S-DPT) (Kutzelnigg 1996), although to the leading order there is hardly a difference between DPT and S-DPT. The result of DPT differs from that of the well-known Breit–Pauli approximation (Bethe & Salpeter 1957) in a correction term *Δ*_{DPT} that takes care of a faster convergence of the overall result with extension of the basis, but which is not very important in the present case, since the nonlinear parameters are well optimized. The relativistic corrections involve expectation values of *c*^{−2} times the following operators:(5.1)(5.2)(5.3)(5.4)(5.5)

There are speculations whether even the Lamb shift (Pyykkö *et al*. 2001) ought to be taken care of for a high-precision study of the IR spectrum of (tables 3 and 4).

## 6. Adiabatic corrections

In the spirit of the BO approximation, one first constructs a PES for clamped nuclei and one considers then the motion of the nuclei in this potential. This is not rigorous, because the nuclei have finite masses, which require corrections to the BO separation of nuclear and electronic motions. Two types of corrections have to be considered:

the

*adiabatic*or*diagonal*corrections, andthe

*non-adiabatic*or*non-diagonal*corrections.

In the traditional evaluation of adiabatic corrections, one first separates the centre-of-mass motions (translation and rotation) from the relative motion and then evaluates the adiabatic corrections in terms of relative nuclear coordinates. This leads to rather complicated expressions. Handy *et al*. (1986) and Handy & Lee (1986) have proposed a much simpler approach in which the centre-of-mass motion is not separated off. One of the present authors has given an *a posteriori* justification (Kutzelnigg 1997) of this approach (see also Herman & Asgharian 1966; Watson 1973) that has been applied to (Cencek *et al*. 1998) and related systems (Cencek & Kutzelnigg 1997). The adiabatic correction *E*_{ad} is essentially the expectation value of the kinetic energy of the nuclei, evaluated with the electronic wave function (exploiting the parametric dependence of the latter on the nuclear coordinates)(6.1)

It is a mass-dependent correction to the PES. The most noticeable adiabatic effect in is found for the isotopomers H_{2}D^{+} and H, in which the *E*-type fundamental *ν*_{2}, degenerate in , is split owing to the asymmetry of the masses (Jaquet *et al*. 1998). If one ignores the adiabatic corrections, the splitting *ν*_{3}−*ν*_{2} is calculated as 128.7 and 111.6 cm^{−1}, respectively. With adiabatic corrections, one obtains 129.5 and 110.2 cm^{−1} in good agreement with the experimental values of 129.6 and 110.3 cm^{−1} (table 5).

It has often been recommended to interpret the atomic energy unit (Hartree) as a *reduced* unit (Shull & Hall 1959), the equivalent of which in standard units (e.g. kJ mol^{−1} or cm^{−1}) depends on the nucleus (by choosing the reduced mass rather than the genuine electronic mass). This choice makes sense only for atoms, and it should *not be used* if one takes care of adiabatic corrections in the way indicated here (Kutzelnigg 1997). This would lead to a double-counting.

## 7. Non-adiabatic effects

For a correct evaluation of non-adiabatic effects, one would have to consider the coupling of several electronic states. Non-adiabatic effects cannot simply be taken care of by a correction to the PES. For the H_{2} molecule, both the adiabatic and the non-adiabatic effects on the vibration frequencies are known rather accurately. The two effects are of comparable magnitude (Jaquet *et al*. 1998). For the fundamental vibrational frequency of H_{2}, the adiabatic correction is −1.4 cm^{−1}, and the non-adiabatic one −0.9 cm^{−1}, the corresponding corrections for HD are −1.0 and −0.6 cm^{−1}, for D_{2} −0.8 and −0.3 cm^{−1}. So the non-adiabatic correction is roughly half the adiabatic one.

There is some indication that the origin of the non-adiabatic correction to vibrational frequencies is the partial participation of the electrons in the vibrational motion, for the analogous effect on the rotation frequencies the participation of the electrons in rotation (Bunker & Moss 1977; Watson 1980; Moss 1996). One can simulate this effect to first-order by a change of the moving mass. For a full participation of the electrons, one must—in neutral systems—replace the nuclear mass that enters the BO approximation by the atomic mass. For this would be the mass of an electron plus two-thirds of that of an electron. There are indications that the effective vibrational and rotational masses are different, e.g. that the participation of the electrons in the rotation is more quenched than that in the vibration.

We do not know much about the magnitude of non-adiabatic corrections for the rovibrational transitions in . We can only conclude from a comparison of the best computed values, including adiabatic corrections, with the experimental values, how large the non-adiabatic corrections might be. They appear to be not much larger than a few tenths of a cm^{−1}, and are hence significantly smaller than for H_{2}, even if one relates them to the smaller vibrational frequencies in .

It is noteworthy that excellent agreement between theoretical and experimental term values (errors of a few hundredth to a few tenths of a cm^{−1}) is achieved if one uses the atomic mass for vibration and the nuclear mass for rotation (Jaquet 1999).

In this section, we have only cared about what one may call *second-order* non-adiabatic effects. They are present in regions where the considered PES is *well separated* energetically from the PES of other states. In the region of avoided crossings, the much larger *first-order* non-adiabatic effects must be considered (see §8).

## 8. Characterization of the ground state

The equilibrium geometry of the singlet state is an equilateral triangle (symmetry *D*_{3h}) with the equilibrium distance *r*_{e}=1.6500 *a*_{0}=0.8731 Å. The binding energy *D*_{e} with respect to 2 H+H^{+} is 0.3438 *E*_{h}=902.6 kJ mol^{−1} or 0.1694 *E*_{h}=444.8 kJ mol^{−1} with respect to H_{2}+H^{+}. That with respect to +H is higher by 0.071840 *E*_{h}=188.6 kJ mol^{−1}.

There are two vibrational frequencies *ν*_{1} and *ν*_{2} of symmetry *A*_{1} and *E*, respectively. Only the degenerate *ν*_{2} is IR-active. The (computed) harmonic frequencies are 3437.5 (*A*_{1}) and 2775.0 cm^{−1} (*E*), while the respective directly observable lowest transition frequencies are 3178.3 and 2521.3 cm^{−1}. These differences indicate a strong anharmonicity, and the necessity to treat the nuclear motion quantum mechanically.

Overtones and combination terms are classified by the occupation of *ν*_{1} and *ν*_{2}, and a vibrational angular quantum number *l*_{2}. For a detailed discussion of quantum numbers for rovibronic levels, see Lindsay & McCall (2001) and Watson (1984, 2000).

The vibrational zero-point energy is 19.87 *mE*_{h} or 4361.4 cm^{−1}. There is also a rotational zero-point energy (e.g. Porter 1982). For *para*-, i.e. for the three nuclear spins coupled to a total spin of *J*=1/2, this is 64.12 cm^{−1} from experiment (Lindsay & McCall 2001) and 64.12 cm^{−1} from our calculation. The respective values for *ortho*- are 86.96 and 86.97 cm^{−1}. Taking care of the zero-point energy and that of H_{2} (2180.2 cm^{−1}), we get the dissociation energy *D*_{0}=34923.6 cm^{−1}=417.98 kJ mol^{−1} of *para*- into H_{2}+H^{+}.

On the ground-state PES there is a saddle point with a symmetric linear structure (Röhse *et al*. 1994) and the HH-bond length *r*_{e}=1.5391*a*_{0} which lies 0.278682 *E*_{h}=731.693 kJ mol^{−1} below 2 H+H^{+} or 0.0651513 *E*_{h}=171.057 kJ mol^{−1} above the *D*_{3h} minimum.

The minimum energy path on the ground-state PES of corresponds to the dissociation of into H_{2}+H^{+} in its ground state. However, one cannot exclude the possibility of a dissociation into +H.

Let *R* be the Jacobi coordinate, which measures the distance between H and H_{2}, and *r* be the HH distance in the or H_{2} subunit. In figure 1, the energy of the lowest singlet states of for *R*=∞ is plotted as a function of *r*. The two lowest curves correspond to the ground state of H_{2} and that of shifted by 0.5 *E*_{h} (Tully & Preston 1971; Bauschlicher *et al*. 1973). For *r* smaller than or slightly bigger than *r*_{e} for H_{2}, H_{2}+H^{+} is lower, but for large *r*, the *charge transfer* state +H is lower. The two curves cross at *r*=2.5017*a*_{0}, but the adiabatic ground state of is, of course, always the state of the lowest energy. For *R* finite, but large, there is a splitting between the two adiabatic states, typical for an *avoided crossing*. This means that in a significant region of the PES, the ground state of changes, near *r*=2.5*a*_{0}, from the *diabatic* state characterized as H_{2} +H^{+} to the charge transfer type *diabatic* state characterized as +H. In the neighbourhood of this avoided crossing the adiabatic corrections become extremely large (Jaquet 1999, 2002) and a dynamic treatment requires to take care of both the states. There is an upper adiabatic state, which switches between the complementary diabatic states.

## 9. Excited states of

A state of a triatomic molecule is regarded as bound, if its PES has (at least) a minimum lower in energy than its fragments, and deep enough to accommodate a zero-point vibration.

Schaad & Higgs (1974) have concluded from CI calculations that has only two bound states, namely the singlet ground state (with *D*_{3h}-symmetry) and an excited linear triplet state, see also Ahlrichs *et al*. (1977). There have been many more recent theoretical studies on the bound triplet state of . For details, see Varandas *et al*. (2005) and the paper by Alijah (2006).

As already mentioned, even the Hückel theory predicts a bound excited state of with linear equilibrium geometry.

This state is much less demanding from the point of view of quantum chemistry. Triplet wave functions vanish at the point of coalescence of two electrons. There is only a cusp of higher order, and the error of a conventional partial wave expansion goes as (*L*+1)^{−5} rather than as (*L*+1)^{−3}, as it does for singlet states. Therefore, standard CI without explicit inclusion of *r*_{12}-dependent terms is sufficient.

In *D*_{3h} symmetry, the triplet state is characterized by the configuration *a*_{1}*e*. As an *E* state it is degenerate and subject to a Jahn–Teller splitting, there are actually two PES with a conical intersection at *D*_{3h} symmetry.

The lower sheet has its minimum for symmetric linear geometry as predicted by the Hückel theory, with *r*=2.4537*a*_{0} (Cernei *et al*. 2003). This state dissociates into H and in its ground state, which is the lowest triplet asymptote. The binding energy *D*_{e} is 0.116 106 *E*_{h} with respect to 2H+H^{+}, or 14 *mE*_{h}*=*3000 cm^{−1} with respect to H and . In *C*_{2v} geometry, there is a saddle point 0.104 075 *E*_{h} below 2H+H^{+} or 2 *mE*_{h} below the asymptote, and 0.012 031 *E*_{h}=2600 cm^{−1} above the minimum.

There is also a singlet counterpart of this triplet state (Tully & Preston 1971; Bauschlicher *et al*. 1973). In regions where the subunits H and are far from each other, the exchange integral is small and the singlet state should only be slightly above the triplet state. Like the triplet state, the excited singlet state should also have a Jahn–Teller splitting. R. Jaquet (2005, unpublished results) indicates that the excited singlet state has a shallow minimum (linear, with *R*=8.5*a*_{0}, *r*=2.0*a*_{0}), but that it is unlikely that it can accommodate a zero point vibrational level.

As already mentioned in the previous section, the excited singlet state has an avoided crossing with the ground state, near the dissociation limit.

As long as the very small spin–orbit interaction can be neglected, a real crossing of the lowest triplet with the lowest singlet state is possible. This requires further investigation.

One may wonder whether has Rydberg states. This is related to the question whether there is a stable configuration of that can serve as the core for an outer electron. Our attempts to construct a PES for led to the conclusion that the absolute minimum occurs for an in its equilibrium and a proton at infinite distance. For any approach of the proton to , the energy rises owing to the Coulomb repulsion. This makes the existence of Rydberg states rather unlikely.

## 10. Analytic fit and accuracy of the PES

In order to evaluate vibrational and rotational term values, an analytic fit to the computed points of the potential surface is necessary. It is crucial that the fit is at least as accurate as the computed points and it interpolates smoothly between the latter, and does not introduce spurious oscillations. It should further extrapolate reasonably to asymptotic regions. Starting with Meyer *et al*. (1986), the MSADC (2.3)–(2.5) became standard, namely *S*_{a}, *S*_{e} and *ϕ*. The energy is then expanded as a polynomial(10.1)

One has the choice between a single fit for the total energy and individual fits for the various contributions to the energy, which are then added up. The latter choice is preferable. In the case of the adiabatic corrections, the separation of the fits (splitting these corrections into two parts) has the advantage that the mass dependence of these corrections is easily taken care of. For more details see Jaquet (2003).

The most important corrections to the BO surface are the adiabatic corrections. In the neighbourhood of the minimum, these vary along the PES between 90 and 140 cm^{−1}, with a contribution of electron correlation to these corrections between 10 and 14 cm^{−1}, but they become very large for *R*>4*a*_{0}. In spite of the smallness of these corrections and their small variation with the geometry (near the minimum), their neglect leads to errors of the order of 1 cm^{−1} in the rovibrational frequencies. Therefore, they are absolutely necessary (figure 2). Relativistic corrections to the PES are negative and of the order of only a few cm^{−1}.

Interesting and promising alternatives to the standard polynomial fit are the *reproducing kernel Hilbert space method*, introduced by Ho & Rabitz (1996) and applied by these authors to the points computed by Röhse *et al*. (1994), as well as the Shepard interpolation of Collins (2002).

Some time ago, when an *ab initio* PES of sufficient accuracy was not available, there have been attempts to *improve* the existing surfaces by introducing adjustable parameters to be fitted to experimental data (Dinelli *et al*. 1994). It is usually possible to adjust even a poor surface, (provided it contains sufficient parameters) such that it reproduces all experimental data (figure 3). This may help to make predictions, but it seldom leads to more physical insight. In order that such an empirical adjustment is meaningful, two conditions must be satisfied (figure 4).

The error in the rovibrational spectrum must be a direct consequence of the inaccuracy of the surface, thus it must not be owing to non-adiabatic effects that cannot be described by a single surface.

All assignments in the observed spectrum must be beyond doubt. Otherwise, one will try to change the parameters in the potential to compensate for a wrong assignment.

There is no guarantee that either of these conditions is satisfied. Moreover, in the present case, the available PES is so accurate that the still existing differences between computed and measured IR frequencies are hardly owing to errors in the PES. An empirical adjustment of a PES may still make sense for larger systems, where the computed PES are necessarily less accurate, but not for .

To account for the low-lying IR frequencies (up to approximately 8000 cm^{−1}), the 69 points of Meyer *et al*. (1986), together with the polynomial fit (10.1) and the program TRIATOM (Tennyson & Sutcliffe 1984; Tennyson & Miller 1989; Sutcliffe & Tennyson 1991) for the nuclear dynamics, turned out to be sufficient. The agreement between theoretical and experimental rovibronic frequencies is quite good (table 6). This even holds for intensities that have been evaluated recently (Jaquet 2002, 2003) for a selected frequency range (figure 3).

For the spectrum in the higher frequency region many more points are necessary. There appears to be the need to further test different fits and an alternative treatment of the nuclear dynamics, e.g. in terms of hyperspherical coordinates (Alijah *et al*. 1995) by means of DVR-methods (Bramley *et al*. 1994; Tennyson *et al*. 1995) or finite elements (Jaquet *et al*. 1997), especially in the region close to and above the barrier to linearity.

If one approaches the dissociation limit, it becomes necessary to consider the nuclear motion on the two singlet surfaces that have an avoided crossing for large *R*.

## 11. A few words on

Under the conditions in which is present, addition complexes (H_{2})_{n} are also stable. Among these, the ion plays a special role. While is the prototype of a 3-centre-2-electron bond, represents a 5-centre-4-electron bond. While the bond between a proton and a H_{2} molecule is particularly strong (*D*_{e}=169 *mE*_{h}=444 kJ mol^{−1}), the gain on attaching to another H_{2} is relatively small (*D*_{e}=13.7 *mE*_{h}=35.9 kJ mol^{−1}=3007 cm^{−1}).

Although has only four electrons, it is much harder than the 2-electron system from the point of view of quantum chemistry. Both the R12-method and the GG-method can be applied to , but the accuracy that can be achieved is much more modest than for . Therefore, calculations with explicitly correlated wave functions were performed only for the most important stationary points of (Müller & Kutzelnigg 2000).

One of the reasons that limited the accuracy is that 4-electron correlation effects, although small, are not entirely negligible. We could only afford the levels CCSD(T)-R12 and CCSD[T]-R12 (which hardly differ among each other) (Müller & Kutzelnigg 2000). CCSD-R12 is not sufficient. The errors of the absolute energies of our calculations are of the order 100 microhartree, but the relative errors are much smaller, of the order 0.1 kJ mol^{−1}, i.e. significantly better than so-called chemical accuracy (Müller & Kutzelnigg 2000).

Early calculations at Hartree–Fock level have predicted that is a rather loose complex between and H_{2}, symmetric structures being much higher in energy. Later it was found (Ahlrichs 1975) that electron correlation strongly stabilizes the symmetric structures to the extent that in the first sophisticated study including electron correlation, the configuration of lowest energy had *D*_{2d} symmetry (Ahlrichs 1975). This prediction could not be confirmed at a still higher level of sophistication.

Ten stationary points of the PES have been characterized (Yamaguchi *et al*. 1987; Müller & Kutzelnigg 2000; Prosmiti *et al*. 2001; Xie *et al*. 2005), four of these are very close in energy (differences of less than 1 *mE*_{h} or 200 cm^{−1}). The other six points are higher by 7–13 *mE*_{h}.

There is no doubt that the structure of lowest energy *C*_{2v}(1)—and the only minimum—consists of an unit and an H_{2} perpendicular to the plane. The binding energy *D*_{e} with respect to and H_{2} (not taking care of zero-point energies) is 13.7*mE*_{h}=35.9 kJ mol^{−1}=3007 cm^{−1} (compared with 169*mE*_{h}=444 kJ mol^{−1}=37091 cm^{−1} for the binding of H^{+} to H_{2}). The next higher structure has *D*_{2d} symmetry, with a proton binding symmetrically to two H_{2} units. This is only 0.28*mE*_{h}=0.74 kJ mol^{−1}=61.5 cm^{−1} above *C*_{2v}(1), followed by *C*_{2v}(2), a planar structure 0.44*mE*_{h}=1.2 kJ mol^{−1}=96.6 cm^{−1} above *C*_{2v}(1) and *D*_{2h}, a symmetric planar structure 0.83*mE*_{h}=182.2 cm^{−1} above *C*_{2v}(1). The structures *C*_{2v}(1) and *C*_{2v}(2) have practically the same geometrical parameters for the subunits and H_{2} and also the same distance between them, they differ only in the *twist angle*. Similarly, the structures *D*_{2d} and *D*_{2h} differ significantly only in the twist angle.

Taking care of the change in the zero-point energy, which has to be done beyond the harmonic approximation (Kraemer *et al*. 1994) using a double-well potential for , we get a binding energy *D*_{0} for with respect to +H_{2} of 25.5 kJ mol^{−1}=213 cm^{−1}, and for of 28.4 kJ mol^{−1}. This can be compared with the recent results obtained by a Monte-Carlo treatment of the zero-point energy (Xie *et al*. 2005), namely 6.33 kcal mol^{−1}=26.5 kJ mol^{−1} for (figure 5, table 7).

The dynamics of has been studied by Kraemer *et al*. (1994) and Spirko *et al*. (1997), and more recently by Moyano & Collins (2003) and Xie *et al*. (2005).

Obviously, is a very floppy molecule, in which all protons are easily scrambled. The Longuet-Higgins group (Longuet-Higgins 1963) has the order 240.

The harmonic vibrational frequencies have been evaluated for the *C*_{2v}(1) structure (Prosmiti *et al*. 2001). The two lowest frequencies are 202 and 506 cm^{−1}. Their zero-point levels lie high above the energetically closest stationary points *D*_{2d}, *C*_{2v}(2) and *D*_{2h}. This is an indication that the four lowest stationary points have no physical reality and that the vibrational ground state is highly symmetric with practically free internal rotation.

One must further consider that along the path through a double well, the zero-point energy of the vibrations perpendicular to this path changes, which manifests itself as a change of the effective isomerization barrier.

## Discussion

B. T. Sutcliffe (*ULB, B 1050 Bruxelles Belgium*). Can we be sure that the triplet state does not have a correlation cusp?

W. Kutzelnigg. The answer consists of two parts.

Kato's cusp relation connects the angular average of the first derivative of the wave function with respect to the internuclear distance at the point of coalescence of two electrons to the value of the wave function at this point

Since for a triplet state, both the wave function and the first derivative vanish at the point of coalescence, the cusp relation is fulfilled trivially in the form 0=0.

If one expands the wave function in the neighbourhood of the coalescence of two electrons in powers of

*r*_{12}, one gets for the spherical average

So even triplet states have a generalized correlation cusp. For details see Kutzelnigg & Morgan III (1992).

A. Alijah (*University of Coimbra, PT 3004-534 Coimbra, Portugal*). For the electronic triplet state, the calculation of the electronic energies still requires considerable effort with traditional methods. We are far away from matching sub-microhartree accuracy as you have achieved for the electronic singlet state. *r*_{12}-type calculations would be useful even for the triplet state.

W. Kutzelnigg. The error of a partial wave expansion (PWE) or of a conventional CI goes as ∼(*L*+1)^{−3} for an ordinary singlet state and ∼(*L*+1)^{−5} for a triplet state, where *L* means the highest angular momentum included in the AO basis. The corresponding error for the R12 method goes ideally as ∼(*L*+1)^{−7} in either case. So, one gains less for a triplet than for a singlet, but one still gains something. The remark that for triplet states traditional CI calculations may be sufficient, was made to explain the audience why the R12 method has so far not been applied to triplet states, but nevertheless acceptable results were obtained.

N. C. Handy (*Cambridge University, CB2 IEW Cambridge, UK*). Professor Kutzelnigg observed that energetic results using correlated Gaussians appeared to be superior to those using linear *r*_{12}. Does not the following rationalize this observation?

The factor shows that the value of the wave function at (and therefore near to) electron coalescence is unimportant.

As

*r*_{12}increases, from coalescence, higher powers of*r*_{12}, , …, become important to describe the shape of the cusp.As

*r*_{12}becomes infinite, linear*r*_{12}diverges and is therefore incorrect, whereas a correlated Gaussian tends to zero.

These three reasons suggest that the use of correlated Gaussians should be superior to the use of linear *r*_{12}. Finally, I am sure that these arguments have been presented previously.

W. Kutzelnigg. It is generally not true that correlated Gaussians are superior to wave functions with linear *r*_{12} terms. It is, on the contrary, surprising that correlated Gaussians perform so well although they do not describe the correlation cusp correctly. They do perform well, only if one uses them in a brute-force way, with a huge number of non-linear variational parameters with a highly sophisticated technique for their optimization. For a small number of parameters, wave functions with term linear (and of higher order) in *r*_{12} are always superior.

In order to understand the situation, one must distinguish between two classes of expansion methods (for details see Kutzelnigg (1994)): (i) those of the type of Fourier series (to which the partial wave expansion, and hence CI belong), and (ii) those which can be rationalized as a discretized integral transformation. (A prototype of this is the expansion of an STO type wave function such as that for the ground state of the H atom in a Gaussian basis.) Methods of type (i) are extremely sensitive to singularities (such as the cusp) of the wave function to be expanded, and display an inverse power-type convergence. Methods of type (ii)—to which also the correlated Gaussian ansatz belongs—are much less sensitive to singularities and obey a convergence law of the typefor *N* the expansion length.

The remarks (i)–(iii) made above are rather irrelevant. The smallness of the domain in which the singularity really matters does not affect the rate of convergence of a Fourier-type expansion. The divergence of the factor *r*_{12} at infinity is more than overcompensated by the exponential decay of the wave function for either *r*_{1}→∞ or *r*_{2}→∞. Recently, Tew & Klopper (2005) have used *r*_{12} multiplied with a damping factor for large *r*_{12} that lead to some, but not a spectacular improvement.

J. Hinze (*University of Bielefeld, D 33615 Bielefeld, Germany*).

Do you not have information of the PES of the excited singlet state of out of your CI calculations of the ground state? Unpublished?

Can we not think of as hydrogen-bonded to H

_{2}? An asymmetric hydrogen bond with a very low barrier for the ‘proton’ motion between the two H_{2}molecules (65 cm^{−1}).

W. Kutzelnigg.

Usually one does get information on excited states of the same symmetry if one performs a CI calculation for the ground state. If one uses the R12 ansatz with a reference function multiplied by

*r*_{12}, one is biased towards the ground state. However, as mentioned in §9, one can use a multi-reference generalization of the R12 ansatz to account for the excited state. Preliminary results for a local minimum of the first excited singlet state are mentioned in §9 as well.One refers to a hydrogen bond, or hydrogen bridge, if a proton binds two bases, e.g. OH

^{−}and H_{2}O to the H_{2}O dimer or two F^{−}to FHF^{−}. The former case represents an asymmetric and the latter a symmetric H-bond. In this sense, can be viewed as consisting of two H_{2}bound by an H-bond. This bond is asymmetric, but not very far from symmetric. While conventional H-bonds are 3-centre bonds, the one in is 5-centre bond.

## Acknowledgments

This work was supported by ‘Fonds der Chemie’. We thank all colleagues involved in this project, especially Wojciech Cencek, Wim Klopper, Jacek Komasa, Hendrik Müller and Robert Röhse.

## Footnotes

One contribution of 26 to a Discussion Meeting Issue ‘Physics, chemistry and astronomy of ’.

- © 2006 The Royal Society