## Abstract

We report calculations of the infrared shifts for the water dimer, as obtained from the recent *ab initio* fully flexible HBB2 potential of Bowman and co-workers. The rovibrational calculations, which formally are 12-dimensional plus overall rotation, were performed within the [6+6]*d* adiabatic separation which decouples the ‘fast’ intramolecular modes from the ‘slow’ intermolecular ones. Apart from this decoupling, each set of modes is treated in a fully variational approach. The intramolecular motion was described in terms of Radau coordinates, using the *f*-embedding formulation of Wei & Carrington, and neglecting the rovibrational Coriolis coupling terms. Within this adiabatic approximation, the intermolecular motion is handled in a similar way as for rigid monomers, except for the rotational constants *B*'s, averaged over intramolecular modes that depend now on the intermolecular geometry. Comparison with experimental data shows an excellent overall agreement.

## 1. Introduction

Numerous efforts have been devoted to the role of hydrated complexes in the atmosphere [1–3] as they have been proposed to contribute to solar absorption, although their precise role is not yet fully understood. The most important hydrated complex is the water dimer (WD) (H_{2}O)_{2} as it has been invoked in several atmospheric processes, such as excess absorption of solar radiation [2,4–7], the water continuum absorption in the far-IR [8,9], homogeneous nucleation of water into droplets and ice [10,11] and catalysis of important chemical reactions [12,13], such as acid rain formation. As its *in situ* experimental study is inherently made difficult owing to the presence of water itself, many theoretical approaches have been undertaken in order to address its spectroscopic properties.

The cornerstone of such theoretical studies is the existence of a potential energy surface describing the WD. It is now well understood [14,15] that a meaningful comparison with experiment requires an extremely accurate description of this potential, at the *ab initio* level. For example, Cencek *et al.* [16] recently produced such a rigid monomer potential (6d) calculated using second-order perturbation theory with Møller–Plesset decomposition of the Hamiltonian (MP2) and the coupled cluster method with single, double and non-iterative triple excitations [CCSD(T)]. This rigid potential, based on vibrationally averaged monomer geometries, leads to a dimer vibration–rotation–tunnelling (VRT) spectrum with the best agreement so far with experiments.

In the present study, we focus on the infrared shifts of the WD that constitute a signature of hydrogen bonding [17]. They are defined as the shifts in the fundamental frequencies of the monomers, when embedded in the dimer. Obviously, such calculations require a flexible monomer potential energy surface, which is 12-dimensional. In a series of papers, Bowman and co-workers [18–20] produced successive generations of such a flexible potential at the [CCSD(T)] level and using the aug-cc-pVTZ basis set. These potentials were expressed as polynomials invariant under all permutations of like atoms, fitted to 30 000 distorted geometries. We consider here only the latest version HBB2 [20]. Table 1 and figure 1 present the resulting properties of the WD at its equilibrium geometry, as compared with the benchmark calculations of Tschumper *et al*. [21]. It can be seen that the main characteristics are very well reproduced by the HBB2 potential.

The outline of this paper is as follows. In §2, we introduce the adiabatic formulation which allows us to handle flexible monomers. Section 3 presents the resulting microwave and far-IR spectra, as well as the infrared shifts, and compares with experiments. Finally, some conclusions are drawn in §4.

## 2. Flexible formulation

In this section, we present the [6+6]*d* adiabatic approximation that we use in order to explicitly deal with flexible monomers. We then describe how this formulation can be recast into a pseudo-rigid calculation for each excited intramolecular state, but using intermolecular geometry-dependent rotational constants on each monomer.

### (a) Adiabatic decoupling

Using the cluster formulation of Gatti & Iung [22], the quantum exact Hamiltonian operator for the fully flexible WD can be formally written as
2.1
where ** Q** stands for the intermolecular coordinates,

**q**

^{X}is the intramolecular coordinates of monomer

*X*,

**j**

_{AB}=

**j**

_{A}+

**j**

_{B}is the coupled internal rotational angular momentum and

**=**

*J***j**

_{AB}+

**is the total angular momentum (**

*L***is the relative angular momentum between the monomers' centres of mass). For further manipulations, we express the vibration–rotation kinetic energy operator (KEO) of each monomer in the following matricial form: 2.2 where is the conjugate momenta of the internal coordinates (**

*L**q*

_{1},

*q*

_{2},

*q*

_{3}), and is the rotational angular momentum components.

To handle this 12-dimensional problem, one can resort to an *adiabatic* separation between the ‘fast’ intramolecular coordinates {**q**^{A},**q**^{B}} and the ‘slow’ intermolecular coordinates **Q**, as applied before by Klopper *et al*. [23,24] in a [4+2]*d* treatment of the HF dimer, and more recently by us for the WD [25]. More specifically, at fixed intermolecular geometry **Q**, we solve for the six-dimensional intramolecular subsystem that provides the adiabatic potentials governing the intermolecular motion. As will be discussed later on, we will further neglect the rovibrational Coriolis coupling terms ** σ** entering equation (2.2)
2.3
leading to the partitioning of the total Hamiltonian according to
2.4
2.5
and
2.6
where and (centrifugal + Coriolis) stand for the first and last terms of equation (2.1), respectively.

### (b) Intramolecular calculation

Within this adiabatic approximation, one has to solve for the intramolecular subsystem at fixed intermolecular geometry **Q**
2.7
where **N** is a composite index representing the six quantum numbers associated to the two monomers, and *Φ*_{N}(**q**^{A},**q**^{B};**Q**) the intramolecular eigenstate at fixed geometry. For convergence efficiency, we express the operators in terms of Radau coordinates [26] (*r*_{1},*r*_{2},*ϑ*)
2.8
which are ideally suited to the H_{2}O molecule as shown by Light and co-workers [27]. We use the *f*-embedding formulation of Wei & Carrington [28,29], choosing the *z*-axis of each monomer frame as the vector which bisects the two Radau vectors, which results in a ** Γ** matrix given by
2.9
where . The resulting Coriolis matrix

**entering equation (2.2) displays only zero elements except for**

*σ**σ*

_{ϑx}=

*σ*

_{xϑ}2.10 which vanishes at equilibrium (

*r*

_{1}=

*r*

_{2}). In the work presented here, this term has been neglected which allows for the decoupling of the intermolecular calculations from the intramolecular ones.

The *V*^{ad}_{N}(**Q**) potential results from a six-dimensional calculation (equation (2.7)), to be performed for every six-dimensional intermolecular **Q**-geometry (≃5×10^{5} **Q**-points). In order to make this efficient, we used the following two step procedure at a given geometry **Q**.

(i) Defining the instantaneous intramolecular optimized geometry from 2.11 and its associated energy , we first compute the vibrational states of monomer

*X*in the field of monomer*Y*≠*X*frozen at its optimized geometry as defined above 2.12 by means of a sequential truncation-reduction scheme [30,31].*E*^{X}_{nX}(**Q**) represents the energy of mode**n**_{X}with respect to the instantaneous energy minimum. The six-dimensional zero-order intramolecular states are then defined as 2.13 with associated energy 2.14 being the zero point energy of an isolated monomer. This corresponds to the [3+3]*d*adiabatic formulation in which one ignores the potential coupling between intramolecular modes, and, for example,*V*^{[3+3]d}_{0}(**Q**) reflects the change in zero point energies of both monomers when the intermolecular geometry**Q**varies. With this definition, the dissociation limit corresponds to zero provided that .(ii) The residual potential term Δ

*V*2.15 can be taken into account within a full variational treatment, the*V*^{6d}_{N}(**Q**) potentials being defined from the eigenvalues of the operator 2.16 the*Φ*_{N}'s being the associated eigenvectors, which corresponds to the 6*d*adiabatic formulation 2.17

These two different definitions of the adiabatic potential will be compared in §2*c*.

In a previous paper [32], we discussed the ambiguity of defining excited intramolecular states in a simple way owing to the presence of multiple minima, and the resulting *G*_{16} symmetry of the system. Here, we briefly recall this discussion, and refer the reader to the paper of Fraser [33] for a more rigorous treatment of this case.

(i) One option is to identify the donor (D) and the acceptor (A) monomers at each intermolecular geometry

**Q**, and to define in such a way D- or A-excited adiabatic potentials. The acceptor can for example be identified as the monomer displaying the shortest*d*_{O…H}distance between its oxygen atom and a hydrogen atom located on the opposite monomer. Such a definition retains the full*G*_{16}symmetry but implies that upon the donor–acceptor interchange motion (which corresponds to bringing the complex from one minimum to another one) the excitation energy simultaneously migrates from one monomer to the other one. Actually, high-resolution spectra [34] show that the corresponding tunnelling splitting is reduced by one order of magnitude, which means that the donor–acceptor interchange is hindered in the excited state.(ii) A second possibility is to stipulate that, whatever is the intermolecular geometry

**Q**, the excitation energy stays on the same monomer, say A or B. In that case, upon interchange, there is no migration of the excitation energy, but only a slight modification owing to the different natures (D or A) of the two monomers. In such a scheme, the*G*_{16}symmetry is lost as the energy is not conserved upon interchange (permutation of the two monomers), and one must use the*G*_{8}subgroup instead.

We will compare these two approximations in §2*c* and show that they lead to close results for the WD.

### (c) Six-dimensional intermolecular calculations

Within the framework defined so far, the dimer full wave function *Ψ*_{N}(**q**^{A},**q**^{B},**Q**) associated to intramolecular state *Φ*_{N} is a solution of
2.18
Neglecting non-adiabatic coupling terms, this wave function can be written as
2.19
and projecting onto the *Φ*_{N}-adiabatic function, one obtains the *intermolecular* equation
2.20
where the notation means averaging over the **q**^{X}-intramolecular coordinates, that is
2.21
In order to simplify the formulation, we diagonalize the resulting *Γ*^{X(N)} matrix, which allows us to define *effective* rotational constants as its eigenvalues.

Equation (2.20) formally corresponds to the rigid formulation of Brocks *et al.* [35], but it uses an adiabatic potential reflecting the total intramolecular excitation energy of the monomers as a function of the intermolecular geometry **Q**. Depending on the method ([3+3]*d* or 6*d*) used to define this adiabatic potential as described in §2*b*, one terms these formulations [6+[3+3]]*d* or [6+6]*d* respectively. A second, important, difference stems from the **Q**-dependence of the rotational matrices *Γ*^{X(N)} of each monomer as shown by equation (2.21). Explicitly handling these **Q**-dependent terms makes the intermolecular energy levels calculation about one order of magnitude more expensive than in a rigid case. We then investigated the approximation of averaging the *Γ*^{X(N)} matrices over the Euler angles
2.22
where *α* is some constant (approx. 10^{2} arb. units), while explicitly retaining the -dependence. This approximation allows us to recast the flexible formulation into a rigid one, except for the rotational matrices which depend then on the separation . Within this averaging approximation, that is
2.23
we will consider the two possibilities of either ignoring the angular dependence of the *Γ*^{X(N)} matrices
2.24
or subsequently retrieving this angular dependence
2.25
by first-order perturbation theory
2.26
We will show in §3 that this latter correction essentially gives the exact transition energies obtained from equation (2.20).

A detailed description of the implementation of the intermolecular calculations has been given previously [36,37] in the case of rigid monomers. Here, for the sake of clarity, we briefly only recall it.

The intermolecular KEO of equation (2.20) leads to simple matrix elements in the overall basis set
2.27
where {|*n*〉} is an appropriate basis for the interfragment distance , a Wigner basis set describing the rotation of each monomer, and {|*J*,*K*,*M*〉} the Wigner basis set associated to the overall rotation of the complex. The basis can be projected onto the different irreducible representations (Irreps) *Γ* of the *G*_{16} (or *G*_{8}) molecular symmetry group governing the energy levels [38]
The most compact representation for the potential energy is the six-dimensional grid , where *φ*=*φ*^{A}−*φ*^{B}. This grid is *restricted* to points where the potential energy is lower than some threshold , and only non-symmetry-equivalent points are computed, which typically correspond to a dimension of *ca* 5×10^{5}. This expensive step of the whole calculation, as the six-dimensional intramolecular system has to be solved for each point of this grid, is fully parallelized by means of the Open MP protocol. Energy levels and eigenstates are then obtained from an iterative Lanczos [39] procedure, which is also parallelized. Transformations between the spectral and grid representations, as required to act the potential operator on the successive Lanczos vectors, are performed by means of a six-dimensional pseudo-spectral scheme involving a three-dimensional fast Fourier transform.

## 3. Results

In this section, we present the results obtained from the new potential for the microwave (MW) and far-IR spectra as well as the IR shifts, and compare with available experimental data. But first, we assess the accuracy of the intramolecular basis sets used, and test different approximations presented in §2.

### (a) Convergence tests

All the calculations were performed with Wigner basis sets, on each monomer, up to *j*=11 for the MW and far-IR spectra and to *j*=10 for the IR shifts (see equation (2.27)), and a primitive radial basis set of 20 sine functions spanning the range [4.2,10] (in bohrs), contracted to nine functions by means of the Harris–Engerholm–Gwinn procedure [40]. These specifications lead to a convergence better than 0.3 per cent on the positions of the energy levels with respect to the ground state, and 1 per cent on the splitting values.

Intramolecular calculations used on each monomer a local adiabatic biharmonic basis set {*H*_{m}(*r*_{1};*ϑ*_{s})×*H*_{n}(*r*_{2};*ϑ*_{s})} defined at each sampled value of the Radau *ϑ* angle grid. The monomer intramolecular basis consisted of seven harmonic functions along each radius *r*_{i} and 19 sampled angles *ϑ*_{s}. The three-dimensional resulting eigenvalues/eigenvectors are then obtained by means of a sequential truncation-reduction scheme [30,31]. Convergence with respect to a larger basis of nine harmonic functions and 25 sampled angles was tested.

We have shown in §2*c* that, within the adiabatic approximation, one can recast the 12-dimensional flexible calculation into a pseudo-rigid one using an adiabatic potential. However, the resulting **Q**-dependence of the rotational matrices ** Γ** renders the Lanczos iterative diagonalization scheme about one order of magnitude more costly because one has to switch from the spectral to the grid representation many times in order to evaluate the effect of the rotational operators . The perturbative approach defined by equations (2.22)–(2.26) allows one to retrieve the exact transition energies and splittings within 0.03 per cent, and the IR shifts within 0.01 cm

^{−1}. It should be noted that ignoring the final perturbation step, equation (2.26), increases these errors to 1 per cent and 0.3 cm

^{−1}, respectively.

As a final test, we compared the dissociation energy *D*_{0}=1096 cm^{−1} obtained using the specifications given above with the experimental value of 1105±10 cm^{−1} recently measured by Reisler and co-workers [41] by velocity map imaging. Another calculation obtained by Shank *et al.* [20] using a diffusion Monte Carlo (DMC) method gives the value 1103±4 cm^{−1}. The very small discrepancy (0.6%) between these two calculations, both within the experimental error, might be due to the adiabatic scheme and the different approximations invoked in it, as the DMC approach is in principle exact.

### (b) Microwave and far-infrared vibration–rotation–tunnelling levels and splittings

In order to assess the global accuracy of the HBB2 potential, we first compare the experimental energy levels [42,43], and those obtained in this study. As the WD is a near symmetric top, the energy *o*_{1} or *o*_{2} of each subfork is customarily defined according to the convention
3.1
where *E*_{i}(*J*,*K*) means the average energy of levels and pertaining to this subfork. The values given in parentheses correspond to the donor–acceptor splitting of the subfork , the acceptor splitting *a*_{i}(*K*) being defined as |*o*_{1}(*K*)−*o*_{2}(*K*)|. We compare the experimental and calculated MW and far-IR spectra in figure 2, using the specifications given above. One difference with respect to our previous studies [16] concerns the reassignment of some energy levels to the OO stretch. From the associated eigenstates *ψ*^{0}_{n}(**Q**), we computed their mean relative kinetic energy . While for most of these states this quantity has a value in the range 30–45 cm^{−1}, OO stretch states are characterized by a value increased to *ca* 80–95 cm^{−1}. In figure 2, the (1) and (2) indices refer, respectively, to the and subforks originating from each (*J*,*K*) level and owing to tunnelling motions (see [16] for a detailed description). One can note that all origins *o*_{i}(*K*) are within 0.5–4 cm^{−1} of their experimental counterparts, and all splittings are reproduced within 30 per cent. Concerning the [6+[3+3]]*d* approximation, i.e. completely neglecting the six-dimensional Δ*V* (**q**^{A},**q**^{B};**Q**) correction term of equation (2.15), the energy levels *o*_{j}(*K*) are barely modified (approx. 0.2% at most), while this effect is slightly larger (approx. 1%) on the donor–acceptor splittings *i*_{j}(*K*).

### (c) Infrared shifts

The upper part of table 2 presents the fundamental frequencies of the isolated water monomer as given by the HBB2 potential. One can note that the stretching modes are in error by roughly 20 cm^{−1} with respect to the experimental values. The reason is that, at the time these calculations were performed, the new HBB/PS potential, which incorporates the empirically adjusted Partridge– Schwenke [51] monomer potential, was not available. However, the IR shifts that we define with respect to the *calculated* isolated monomer frequencies are expected not to depend too much on these discrepancies.

The lower part concerns the IR shifts as defined above. For all modes but *as*[*A*], we report the shifts associated to the lower transition. In the latter case, owing to the vibrational asymmetry of the upper state, it corresponds to the observed transition. First, the right two columns consider the two limiting models, *G*_{16} or *G*_{8}, where in the first case vibrational excitation hops from one monomer to the other one upon donor–acceptor interchange, while it stays on a given monomer in the second case. For the *G*_{8} model, the acceptor can be identified *a posteriori* as the monomer displaying the shortest hydrogen bond O…H. These results were obtained in the [6+[3+3]]*d* approximation, i.e. completely neglecting the six-dimensional Δ*V* residual correction term of equation (2.15). It can be seen that the two models give extremely close results, the differences being at most of the order of 1.5 cm^{−1}. Comparison of columns 2 and 3 allows us to estimate the importance of this residual term. It leads to minor effects (approx. 0.5 cm^{−1}) except for the O–H_{b}[*D*] mode where the change is 1.5 cm^{−1}. Concerning the *ss*[*A*] mode, it should be noted that its value was taken from helium droplets experiments [50] as no molecular beams data are available. For most of the modes, it should be emphasized that it is not possible to give an accurate comparison between experimental and calculated IR shifts as no resolved band was identified, except for the acceptor asymmetric stretch *as*[*A*] [34]. For these former modes, the experimental data essentially provide estimates of the IR shifts.

Except for the donor O–H_{b}[*D*] and acceptor asymmetric *as*[*A*] modes, the calculated IR shifts agree remarkably well with the experimental estimated values. To understand these discrepancies, it should be first kept in mind that these shifts actually correspond to changes in large frequencies, of the order of 3600–3700 cm^{−1}, which makes their accurate determination challenging. The second point concerns the functional form used to fit the *ab initio* data. Contrary to other flexible potentials, which separate out the contribution of one-body terms
the HBB2 potential was directly fitted in terms of the 15 interatomic distances {*r*_{ij}}
Such a functional form, which does not rely on a quasi-exact monomer potential reference, certainly makes more difficult the accurate description of excited monomers, as shown in table 2 for the fundamental frequencies of isolated monomers. In particular, this effect should be more pronounced for O–H_{b}[*D*] and *as*[*A*] modes which correspond to highly distorted monomer geometries. Recent MULTIMODE calculations by Shank *et al.* [20] comparing the HBB2 and HBB2/PS potentials, the latter one incorporating the empirically adjusted Partridge–Schwenke monomer potential [48], did show improved IR shifts for these two modes.

### (d) (D_{2}O)_{2} results

These calculations were conducted with a larger intermolecular basis (up to *j* = 13) in order to converge the very small shifts of the ground state levels caused by bifurcation tunnelling, and within the approximate [6+[3+3]]*d* treatment. The reason is that, owing to the more localized nature of the acceptor and donor intramolecular states resulting from the heavier deuterated hydrogens, coupling between these states is much less effective. As an example, the change in zero point energy for the dimer at its equilibrium geometry is only 0.06 cm^{−1} when considering the full variational [6+6]*d* approach as compared with the approximate [6+[3+3]]*d* one.

Figure 3 presents a comparison of the calculated VRT levels and splittings with the experimental data [42–44,52–54]. The ground state levels and the very small splittings are well reproduced except for the *o*_{2}(0) level which displays a slow convergence with respect to the basis set size. One can note that the acceptor wag and acceptor torsion levels are particularly well described, within 0.2–0.8 cm^{−1} of the experimental values, with all associated splittings but one being off by 20 per cent at most. Discrepancies are somewhat larger concerning the other modes.

## 4. Discussion

We have tested the spectroscopic properties of the new, fully flexible, *ab initio* HBB2 potential of Bowman and co-workers [20]. Such a calculation constitutes a challenging task as it corresponds formally to a 12-dimensional problem. We have used the general cluster formulation of Gatti & Iung [22] and described each monomer in terms of *f*-embedded Radau coordinates, as worked out by Wei & Carrington [28,29], which results in a very simple expression of the KEO. The large separation in frequencies between the ‘*fast*’ intramolecular modes and the ‘*slow*’ intermolecular ones allowed us to invoke an adiabatic [6+6]*d* approximation where one has to solve the intramolecular subsystem at each six-dimensional intermolecular geometry. The whole calculation is thus recast in a way similar to the rigid monomers case, but using rotational constants *B*'s explicitly depending on the geometry. We showed that this dependence, which increases considerably the computational effort, could actually be retrieved by a perturbation correction at a much lower cost comparable to the one required for rigid monomers. Calculations were enhanced by using a straightforward, Open MP, parallel implementation which greatly reduced the large amount of computing time required. The [6+[3+3]]*d* approximation, which neglects the potential coupling between intramolecular modes, was shown to have a minor effect on spectroscopic properties, much smaller than the errors with respect to experimental values.

This new HBB2 potential constitutes a definite improvement for the description of IR shifts over the previous HBB version [18,19,32]. As compared with experiments, it leads to an accuracy somewhat better than the results obtained from the Sapt-5s*f* potential [55]. The newer HBB2/PS version, which corresponds to asymptotically correct one-body terms, should be tested in order to check if it still improves these shifts. Finally, it might be necessary to explicitly consider a non-adiabatic formulation, the impact of which has never been estimated. It would in particular allow for a rigorous description of the intramolecular levels, instead of resorting to the two limiting models (*G*_{16} or *G*_{8}) invoked here.

## Acknowledgements

Dr F. Gatti is acknowledged for helpful discussions.

## Footnotes

One contribution of 17 to a Theo Murphy Meeting Issue ‘Water in the gas phase’.

- This journal is © 2012 The Royal Society