We present a rigorous and comprehensive theoretical treatment of the vibrational dynamics of benzene–H_{2}O and benzene–HDO dimers, where the quantum bound-state calculations of the coupled intra- and intermolecular vibrational states of the dimers are complemented by the quantum simulations of their infrared (IR) and Raman spectra utilizing the computed eigenstates. Apart from taking benzene to be rigid, the methodology for the nine-dimensional (9D) vibrational quantum calculations introduced in this study is fully coupled. The approach yields the intramolecular vibrational fundamentals and the bend (*ν*_{2}) overtone of H_{2}O and HDO in the complex, together with the low-lying intermolecular vibrational states in each of the intramolecular vibrational manifolds considered. Following the recently introduced general procedure [P. M. Felker and Z. Bačić, J. Chem. Phys. **151**, 024305 (2019)], the full 9D vibrational Hamiltonian of the dimer is divided into a 6D intermolecular Hamiltonian, a 3D intramolecular Hamiltonian, and a 9D remainder term. A 9D contracted product basis is constructed from the low-energy eigenstates of the two reduced-dimension Hamiltonians, and the full vibrational dimer Hamiltonian is diagonalized in it. The symmetry present in the dimers is exploited to reduce the Hamiltonian matrix to a block diagonal form. Guided by the findings of our earlier study referenced above, the 6D intermolecular contracted bases for each symmetry block include only 40 eigenstates with energies up to about 225 cm^{−1}, far below the stretch and bend fundamentals of H_{2}O and HDO, which range between 1400 cm^{−1} and 3800 cm^{−1}. As a result, the matrices representing the symmetry blocks of the 9D Hamiltonian are small for the high-dimensional quantum problem, 1360 and 1680 for the H_{2}O and HDO complexes, respectively, allowing for direct diagonalization. These calculations characterize in detail the H_{2}O/HDO intramolecular vibrations, their frequency shifts, and couplings to the large-amplitude-motion intermolecular vibrational sates. The computed IR spectra of the two complexes in the OH-stretch region, as well as the intermolecular Raman spectra, are compared to the experimental spectra in the literature.

## I. INTRODUCTION

Size-selected small benzene–(H_{2}O)_{n} clusters with *n* = 1–7, as well as the benzene–HDO dimer, formed in a supersonic expansion, have been intensely investigated experimentally by means of microwave,^{1,2} infrared (IR),^{3–6} and Raman^{7} spectroscopy. The microwave and Raman spectroscopies probed the two dimers in the ground intramolecular vibrational state, providing information about their vibrationally averaged geometries^{1,2} and intermolecular vibrations,^{7} while the IR spectra showed the coupled intra- and intermolecular vibrational excitations in the region of the OH stretch fundamental. Prior to these studies in molecular beams, the IR spectra of the 1:1 complexes of benzene with water in solid argon were recorded and analyzed.^{8} Very recently, a far-IR (FIR) spectroscopic investigation of the intermolecular vibrations of the benzene–water dimer in neon and argon matrices was reported.^{9} There are multiple reasons for the attention given to these clusters. First, benzene is a prototypical aromatic molecule, and its microsolvation by a variable number of water molecules opens a venue for achieving a quantitative, molecular-level understanding of the interactions of aromatic hydrocarbons with water as the solvent. Second, benzene is incorporated into proteins through the side chains of aromatic amino acids (e.g., tyrosine, phenylalanine, and tryptophan), and their interactions with water play an important role in shaping the tertiary structure of the proteins through the hydrophobic effect. Third, the interaction between benzene and water is unusual, involving the so-called *π* hydrogen bond, in which water donates one or both of its hydrogens to the delocalized *π*-electron cloud above the benzene plane, rather than to a localized electron pair on the proton acceptor molecule.

Considerable experimental effort has been directed at the benzene–H_{2}O and benzene–HDO dimers.^{1–3,5–7,9} A strong motivation for it has been the expectation that this will enable the determination of an accurate potential energy surface (PES) of the dimer, preferably a nine-dimensional (9D) one, for flexible H_{2}O and rigid benzene, through interaction with high-level theory. Such a dimer PES would provide an essential two-body term in the many-body expansion for the noncovalent interactions in larger benzene–water clusters and ultimately the condensed phase. Another major impetus was provided by the realization that benzene–H_{2}O and benzene–HDO are highly nonrigid, which emerged first from the matrix isolation study.^{8} This fluxional character was confirmed and characterized more precisely by the subsequent microwave^{1,2} and IR spectra.^{3,6} The picture that emerged from these studies has the vibrationally averaged position of the center of mass (c.m.) of the water molecule on the *C*_{6} axis of benzene, with both hydrogens oriented preferentially toward the ring in the *π* hydrogen-bonding configurations. The water molecule undergoes multiple intermolecular large-amplitude motions (LAMs) already in the ground vibrational state, in particular the internal rotation about the sixfold axis and the facile in-plane torsional (or rocking) motion that exchanges its two H atoms, in such a way that, on the vibrationally averaged basis, the complex retains the sixfold symmetry of the benzene ring. In benzene–HDO, HDO preferentially forms the *π* hydrogen bond via the D atom.^{5,8} The participation of both H atoms of H_{2}O in the *π* hydrogen bonding with benzene is limited to the 1:1 complex. In larger benzene–(H_{2}O)_{n} clusters, with *n* ≥ 2, benzene resides on the exterior, or the surface, of the (H_{2}O)_{n} cluster, which is *π* hydrogen bonded to benzene through a single free OH group.^{4,5}

The intermolecular LAMs are responsible not only for the vibrationally averaged sixfold symmetry of the benzene–water dimer. The low-energy intermolecular levels associated with them give rise to a number of transitions in the OH stretch region of the IR spectrum, which in the absence of the LAMs would be negligibly weak.^{3,6} Thus, in the OH-stretch region of benzene–HDO, five transitions are observed originating in the zero-point-energy level, the OH-stretch fundamental, and four higher-energy bands corresponding to the combined excitations of the OH stretch and various low-frequency intermolecular vibrations.^{3,6} Interestingly, these combination bands have significantly greater intensities than the OH-stretch fundamental. At least six combination bands of similar origin are present in the asymmetric-stretch region of the IR spectrum of benzene–H_{2}O, while no such combination bands are observed in the symmetric-stretch region.^{3,6} For both benzene–H_{2}O and benzene–HDO, these combination bands, with excitation energies up to ≈80 cm^{−1} relative to the respective intramolecular fundamentals, provide further quantitative evidence regarding the LAMs and their coupling to the H_{2}O/HDO intramolecular vibrations. The features arising from the LAMs are present also in the Raman spectra of benzene–H_{2}O, as four prominent intermolecular bands in the low-energy region below 60 cm^{−1},^{7} and in the FIR spectra of the matrix-isolated benzene–water complex.^{9} However, these FIR spectra do not extend below ≈200 cm^{−1}.

Theoretical investigations of benzene–H_{2}O have a long history. The six-dimensional (6D) intermolecular potential energy surface (PES) for this dimer was first calculated by Karlström *et al.*,^{10} at the *ab initio* self-consistent field (SCF) CI level and treating both monomers as rigid. Suzuki *et al.*^{1} complemented their microwave spectroscopy investigations of benzene–H_{2}O (and isotopologues) with a scan of the intermolecular (rigid monomer) PES of the dimer, using the MP2 level of theory. The results obtained from their electronic structure calculations are in qualitative agreement with experiment: c.m. of H_{2}O on the *C*_{6} axis of benzene, with both H atoms oriented toward the plane of benzene, and the potential that is very flat over the benzene surface. The benzene–H_{2}O intermolecular interaction, for rigid monomers, was also characterized by means of density functional theory (DFT) calculations.^{11} For benzene–(H_{2}O)_{n} clusters with *n* ranging from 1 to 10, the optimized minimum-energy cluster structures, their binding energies, the harmonic frequencies, and IR intensities have been computed^{9,12,13} for the purpose of comparison between the calculated and measured IR spectra of these clusters in the OH-stretch fundamental region. The OH-stretch spectra of benzene–(H_{2}O)_{n} clusters, *n* = 2–7, were modeled by utilizing the local-mode approach as well.^{14} In both the harmonic and local-mode calculations, the cluster structures are taken to be static, thus excluding the low-frequency intermolecular vibrations. Nevertheless, they can be valuable for interpreting the IR spectra of benzene–(H_{2}O)_{n} clusters with *n* ≥ 2 where, as mentioned above, the (H_{2}O)_{n} cluster donates a single free OH group to the *π* hydrogen bond, and the cluster rigidity increases with *n*. However, benzene–H_{2}O/HDO is highly fluxional, and theoretical approaches assuming the minimum-energy structure of the dimer and small harmonic displacement from it are not appropriate. The LAMs must be included in any realistic theoretical treatment of the dynamical and spectroscopic features of these complexes.

Three 6D (rigid-monomer) quantum calculations of the LAMs of benzene–H_{2}O have been reported.^{15–17} One of them was a diffusion Monte Carlo simulation,^{16} providing information regarding the vibrationally averaged properties of the dimer in the ground vibrational state only. In the quantum 6D approach of Linse,^{15} the intermolecular vibrational levels of rigid H_{2}O inside the external potential of benzene molecules are calculated, but under the simplifying assumption that benzene is infinitely heavy, thereby leaving out a great deal of angular momentum coupling present in the complex. The only dynamically rigorous quantum 6D (rigid-monomer) calculation of the intermolecular vibrational states of benzene–H_{2}O is that of Kim *et al.*,^{17} which proved to be successful in identifying some significant features of the low-energy states of the complex.

A comprehensive, quantitative treatment of the vibrational level structure of benzene–H_{2}O/HDO, which can account for the combination bands prominent in the measured IR and Raman spectra as well as the OH-stretch vibrational frequency shifts, cannot be achieved while treating H_{2}O/HDO as rigid. A rigorous approach to the H_{2}O intramolecular vibrations coupled to the intermolecular vibrational states of the dimer requires fully coupled 9D quantum bound-state calculations for flexible H_{2}O and rigid benzene, extending to the energies of the H_{2}O vibrational fundamentals and overtones. The realization of this approach has to overcome two serious obstacles. One of them is the need to deal efficiently and accurately with nine coupled anharmonic vibrational degrees of freedom (DOFs), of which the six intermolecular DOFs have large-amplitude character. The second one is described in the following.

The measured frequencies of the intramolecular vibrational fundamentals of H_{2}O in the gas phase for the symmetric (*ν*_{1}) and asymmetric stretch (*ν*_{3}) modes and the bend (*ν*_{2}) are 3657.01 cm^{−1}, 3755.9 cm^{−1}, and 1594.7 cm^{−1}, respectively.^{18} Likewise, for the gas-phase HDO, the measured fundamental frequencies of *ν*_{OD}, *ν*_{2}, and *ν*_{OH} modes are 2723.7 cm^{−1}, 1403.5 cm^{−1}, and 3707.5 cm^{−1}, respectively.^{18} In contrast, the fundamentals of virtually all intermolecular modes of benzene–H_{2}O/HDO are below 100 cm^{−1}. As a result of such a large disparity between the frequencies of the intra- and intermolecular vibrational modes, typical for weakly bound molecular dimers, hundreds, if not thousands, of very highly excited intermolecular states in the intramolecular ground-state manifold lie below the energies of the H_{2}O and HDO intramolecular vibrational excitations considered. Until our recent work,^{19–22} the prevailing view was that for weakly bound molecular dimers, fully coupled quantum bound-state calculations of their intramolecular vibrational fundamentals and overtones require converging either all intermolecular vibrational states below the energies of these intramolecular excitations or at least a dense set of highly excited intermolecular states energetically close to them. If this were true, then the intramolecular vibrational fundamentals of H_{2}O/HDO in the complexes with benzene, coupled to the intermolecular states in these manifolds, would be virtually intractable for rigorous quantum treatment.

However, our recent fully coupled quantum 6D calculations of the vibration–translation–rotation (VRT) eigenstates of H_{2}, HD, and D_{2} inside the clathrate hydrate cage,^{19} taken to be rigid, led to a very surprising finding: the first excited ($v$ = 1) vibrational state of the caged H_{2} at ≈4100 cm^{−1} could be calculated accurately while having converged only a modest number of translation-rotation (TR) states in the $v$ = 0 zero manifold up to at most 400–450 cm^{−1} above the ground state and none several thousand wave numbers below the H_{2} intramolecular fundamental. Evidently, their absence affects negligibly the accuracy with which this high-energy intramolecular excitation is calculated. This was explained in terms of extremely weak coupling between the $v$ = 1 state of H_{2} and the highly excited $v$ = 0 intermolecular TR states close in energy. The resulting large reduction in the dimension of the basis set transformed a formidable computation into one that is readily tractable.^{19}

We immediately recognized the relevance of this insight for a broad class of weakly bound, van der Waals and hydrogen bonded (HB) molecular dimers in which the frequencies of the monomer intramolecular vibrations are typically at least an order of magnitude higher than those of the intermolecular modes. This led us to develop a very efficient method for high-dimensional, fully coupled quantum calculations of intramolecular vibrational excitations of such complexes, together with the intermolecular vibrational states in each intramolecular vibrational manifold.^{20} This approach involves the partitioning of the vibrational Hamiltonian of the dimer into two reduced-dimension Hamiltonians, a rigid-monomer one for the intermolecular vibrations and the other for all intramolecular vibrational DOFs, and a remainder. The eigenstates of these two reduced-dimension Hamiltonians are used to build up a product contracted basis in which the full vibrational Hamiltonian of the dimer is diagonalized. Such a basis is particularly convenient for taking advantage of what we learned previously about the weak coupling between the intra- and intermolecular degrees,^{19} as it enables including only a small number of lowest-energy intra- and intermolecular vibrational eigenstates in the final product contracted basis.^{20}

For polyatomic molecules, e.g., H_{2}O_{2}, CH_{4}, and $CH5+$, Carrington et al.^{23–26} implemented similar ideas of dividing the internal coordinates into two groups, stretch and bend, and used the eigenvectors of the corresponding reduced-dimension Hamiltonians in the final product contracted basis. However, due to strong stretch–bend coupling in these molecules, it was not possible to include only a small number of contracted functions for one group, bend or stretch, in the final product contracted basis, in contrast to the weakly bound molecular dimers.

In the application to the full-dimensional (6D) vibration–rotation–tunneling eigenstates of (HF)_{2}, a prototypical HB dimer, the results of our approach^{20} are in excellent agreement with those in the literature, while using just a fraction of the basis-set size required by the other methods. Subsequently, the same approach was utilized for rigorous, yet remarkably efficient, 8D quantum calculations of the intramolecular stretch fundamentals and frequency shifts of (H_{2})_{2} in the large hydrate cage^{21} and in the fully coupled 9D quantum treatment of the intra- and intermolecular vibrational levels of H_{2}O in (rigid) C_{60},^{22} which would otherwise be prohibitively costly.

In this paper, we present the computational methodology for performing, with modest effort, rigorous, fully coupled 9D quantum calculations of the intramolecular vibrational excitations of flexible H_{2}O and HDO weakly bound to rigid benzene, together with the low-energy intermolecular vibrational states within each intramolecular vibrational manifold. The method incorporates our approach outlined above^{20} and also builds on the earlier quantum 6D (rigid-monomer) treatment of the intermolecular vibrations of benzene–H_{2}O.^{17} The bound-state problem is formulated in a way that is particularly convenient for weakly bound complexes where one molecule (e.g., benzene) is significantly larger than the other (e.g., H_{2}O), and as a result, the intermolecular interaction shows pronounced anisotropy. Symmetry factorization of the 9D basis reduces the matrix of the full 9D Hamiltonian to a block-diagonal form, each of the blocks corresponding to one of the irreducible representations (irreps) of the *G*_{12} molecular symmetry group. Their dimensions, 1360 and 1680, for each *G*_{12} irrep of benzene–H_{2}O and benzene–HDO, respectively, are small enough to allow for direct diagonalization. The resulting 9D eigenstates are used to compute the IR spectra of benzene–H_{2}O/HDO in the OH-stretch fundamental region. Also calculated is the Raman spectrum of benzene–H_{2}O in the low-frequency region corresponding to transitions between intermolecular vibrational states within the intramolecular vibrational ground-state manifold of the dimer. Comparison is made between the measured IR and Raman spectra of these species in the literature.

## II. COMPUTATIONAL METHODOLOGY

### A. Vibrational Hamiltonian of the dimer

Among the numerous ways to formulate the problem of vibrational states in a weakly bound molecular complex, an approach wherein the complex-fixed (CF) axis frame is defined in terms of the nuclear positions of one of the component moieties is particularly convenient when that moiety is significantly larger than the other component, and the intermolecular interaction, as a result, is significantly anisotropic. This is the situation that obtains for the benzene–water (BW) dimer. As such, in an earlier work^{17} on the intermolecular vibrational states of this complex (in the rigid-monomer approximation), we started with a CF frame affixed to the benzene moiety, derived an intermolecular vibrational Hamiltonian in terms of coordinates referred to that frame, and computed intermolecular states by numerical diagonalization of that Hamiltonian. We extend that approach here to the BW complex in which the benzene moiety is still taken to be rigid, but the water is treated as fully flexible.

We label the C nuclei of the benzene #1 to 6 in order about the benzene ring. We define the CF frame such that its origin is at the c.m. of the complex, its $X^$ axis is parallel to the vector going from the c.m. of the benzene moiety to C(1), its *Ŷ* axis is parallel to the vector going from the c.m. of benzene to the middle of the C(2)–C(3) bond, and its *Ẑ* axis is $X^\xd7\u0176$. We can then write the classical nuclear kinetic energy of this system (apart from overall translational energy) as

where *μ*_{d} is the reduced mass of the complex, *m*_{water}*m*_{benzene}/(*m*_{water} + *m*_{benzene}), **p**_{d} is the momentum conjugate to **d**, the vector going from the c.m. of benzene to that of the water moiety, the index *i* runs over the CF-axis directions, $jiB$ is the component of the rotational angular momentum of the benzene measured along the *i*th axis, *I*_{i} is the moment of inertia about the benzene principal axis parallel to *i*, and *T*^{WF} is the kinetic energy of the water moiety due to the motion of its nuclei referred to a water-fixed (WF) frame.

Starting with Eq. (1) and proceeding with a development analogous to that mentioned in Sec. III A of Ref. 17, we obtain the quantum vibrational Hamiltonian of the BW complex as

where *ω* ≡ (*ϕ*, *θ*, *χ*) denotes the Euler angles that fix the orientation of the WF frame relative to the CF frame, **q** denotes the vibrational coordinates of the water moiety, $\u2207d2$ is the Laplacian associated with **d**, $L^i$ is the operator corresponding to the component, measured along the *i*th axis, of the orbital angular momentum associated with the rotation of the c.m.s of the benzene and water moieties about the c.m. of the dimer, $\u0135iW$ is the operator corresponding to the component of water’s angular momentum measured along the *i*th CF axis,

is the operator associated with the rotational–vibrational kinetic energy of the water monomer^{27} (explicit expressions for $T^vrWF$, $T^corWF$, and $T^vibWF$ are given in Sec. II C), and *V*_{tot}(**d**, *ω*, **q**) is the 9D PES of the benzene–water dimer. We take *V*_{tot} to have the form

### B. General procedure for calculating the eigenstates of *Ĥ*

To solve for the eigenstates of *Ĥ* in Eq. (2), we follow the general procedure of Ref. 20 and divide it into a reduced dimension (6D) “intermolecular” Hamiltonian *Ĥ*_{inter}, a reduced dimension (3D) “intramolecular” Hamiltonian *Ĥ*_{intra}, and a remainder term Δ*Ĥ*. From the low-energy eigenstates of *Ĥ*_{inter} (which we denote as |*κ*⟩, with eigenenergies $E\kappa inter$) and those of *Ĥ*_{intra} (which we denote as |*γ*⟩, with eigenenergies $E\gamma intra$), we build up a contracted product basis |*κ*, *γ*⟩ ≡|*κ*⟩|*γ*⟩, which we then use to diagonalize the full Hamiltonian. We define *Ĥ*_{inter} as

where

and **q**_{0} denotes vibrational coordinates of the water moiety fixed to constant values near to those characterizing the minimum-energy geometry of the water monomer. For *Ĥ*_{intra}, we have

where the “dimer-adapted” 3D intramolecular potential, *V*_{adap}, is the hydrogen-exchange-symmetrized version of

Here, **d**_{0} and *ω*_{0} are fixed values of **d** and *ω*. We choose these values to correspond to those at, or near, the minimum of *V*_{B-W} when the **q** are fixed to **q**_{0}. The resulting *V*_{adap} thus includes some of the effects of the benzene–water interaction on the water vibrational PES.

This leaves

where

and

The matrix elements of *Ĥ* in the 9D |*κ*, *γ*⟩ basis are thus given by

### C. Coordinates and kinetic-energy operators

We now specify our choices for the particular coordinates on which *Ĥ* depends. For the components of **d**, we use the cylindrical coordinates (*d*_{Z}, *ρ*, Φ), where *d*_{Z} is the cartesian component of **d** along the *Ẑ* axis, $\rho \u2261dX2+dY2$, and (cos Φ, sin Φ) ≡ (*d*_{X}/*ρ*, *d*_{Y}/*ρ*), with *d*_{X} and *d*_{Y} being the Cartesian components of **d** along the $X^$ and *Ŷ* axes, respectively. For the intra-water **q**, we use the Radau coordinates (*R*_{1}, *R*_{2}, Θ).^{29} Furthermore, we use Radau bisector-*z* embedding to define the WF frame.^{27} Specifically, we take the WF *Ẑ* axis to be parallel to $\u2212(R^1+R^2)$, where $R^1$ and $R^2$ are the unit vectors associated with the Radau vectors **R**_{1} and **R**_{2}. The WF *ŷ* axis is taken to be parallel to $R^1\xd7R^2$, and the WF $x^$ axis is given by *ŷ* × *ẑ*. With the CF and WF frames defined, the Euler angles *ω* = (*ϕ*, *θ*, *χ*) can be specified. We use the convention employed in Ref. 30 with CF being the “space-fixed” frame and WF being the “body-fixed” frame.

where

and

where *A*_{⊥} ≡ 1(2*I*_{X}) = 1/(2*I*_{Y}) and *A*_{∥} ≡ 1/(2*I*_{Z}) = *A*_{⊥}/2 are the rotational constants of rigid benzene.

and

where $Bi(Ri)\u22611/(2\mu iRi2)$, *μ*_{i} is the Radau mass factor associated with *R*_{i},^{29} the $\u0135kWF$ (*k* = *x*, *y*, *z*) are the operators associated with the components of the water moiety’s rotational angular momentum measured along the WF axes, and

Note that the $(\u2009\u0135XW,\u0135YW,\u0135ZW)$ appearing in Eqs. (15) and (16) correspond to components of the water angular momentum measured along the CF axes, while the $(\u2009\u0135xWF,\u0135yWF,\u0135zWF)$ in Eqs. (17) and (18) refer to such components measured along the WF axes. The former conform to the normal angular-momentum commutation relations, while the latter conform to the anomalous ones.

### D. Diagonalization of *Ĥ*_{inter}

#### 1. Primitive intermolecular basis set

We choose a “cylindrical-Wigner” basis^{17} to represent the matrix of *Ĥ*_{inter}. These basis states are of the form

where the |*d*_{Z,α}⟩, *α* = 1, … *N*_{α}, comprise a 1D Gauss–Hermite DVR covering the *d*_{Z} coordinate, the |$v$ , *l*⟩ (|*l*| = 0, 1, …, $vmax$ and $v$ = |*l*|, |*l*| + 2, …, $vmax$) are eigenfunctions of a degenerate 2D harmonic oscillator (HO) dependent on *ρ* and Φ, and the | *j*, *m*_{j}, *k*⟩ (*j* = 0, 1, …, *j*_{max}, *m*_{j}, *k* = −*j*, −*j* + 1, …, *j*) are symmetric-top rotational eigenfunctions dependent on the *ω*. This type of basis was employed in our earlier work on benzene–water^{17} and was described in detail there. In the current work, we take *N*_{α} = 20, $vmax$ = 10, and *j*_{max} = 11. The |*d*_{Z,α}⟩ were generated from 1D HO wavefunctions centered at *d*_{Z} = 7.6 bohr (*z*_{0} in Ref. 17) with the mass × angular-frequency parameter ($\gamma z2$ in Ref. 17) set to (1.428 778)^{2} a.u. The resulting 20 *d*_{Z} quadrature points ranged from 3.829 to 11.371 bohr. The 2D HO eigenfunctions had their mass × angular-frequency parameter ($\gamma \rho 2$ in Ref. 17) set to (1.058 354 4)^{2} a.u. We have assessed the adequacy of this basis with respect to the convergence of the eigenstates of both *Ĥ*_{inter} and, more importantly, *Ĥ*. As to the former, we find that the states are generally converged to within several tenths of cm^{−1} for intermolecular excitations less than 250 cm^{−1}. A similar level of convergence with respect to the primitive intermolecular basis pertains to the 9D eigenstates that we report below.

#### 2. Symmetry and the primitive basis

Symmetry factorization of the cylindrical-Wigner basis ultimately allows for a significant increase in computational efficiency with respect to the diagonalization of *Ĥ*_{inter} and, more critically, with respect to the 9D diagonalization of *Ĥ*. We choose to effect such factorization into symmetry blocks associated with the irreducible representations (irreps) of the *G*_{12} molecular symmetry group (isomorphic to the *D*_{6} point group with the character table given by the upper-left quadrant of Table I in Ref. 3). With this choice, we are able to employ the same code to handle both the H_{2}O and HDO versions of the single-sided BW complex.

The first step in effecting such symmetry factorization is to include only those basis states with the same value of (*l* + *m*_{j}) mod 6 in a given symmetry block. This divides the states into six sets corresponding to the *A* irreps [(*l* + *m*_{j}) mod 6 = 0], the *B* irreps [(*l* + *m*_{j}) mod 6 = ±3], the two components of the *E*_{1} irrep [(*l* + *m*_{j}) mod 6 = +1 and −1], and the two components of the *E*_{2} irrep [(*l* + *m*_{j}) mod 6 = +2 and −2]. With the values of *N*_{α}, $vmax$, and *j*_{max} that we use, these symmetry blocks are composed of 505 760, 506 240, 506 240, and 505 760 states for the *A*, *B*, *E*_{1} (either component), and *E*_{2} (either component) irreps, respectively.

The second step, which factors the *A* block into *A*_{1} and *A*_{2} blocks and the *B* block into *B*_{1} and *B*_{2} blocks, entails sorting states by the way they transform with respect to the operation (35)(26)^{*}. Rather than building such symmetry directly into the primitive basis, we impose it upon the random initial state vector that serves as the starting point in the iterative matrix-on-vector scheme, which we employ to diagonalize *Ĥ*_{inter}. Specifically, we project out of that state vector an eigenfunction of (35)(26)^{*} with the desired eigenvalue *ϵ* = ±1. Working with the latter state function ultimately produces eigenfunctions of *Ĥ*_{inter}, which are simultaneous eigenfunctions of (35)(26)^{*} with eigenvalue *ϵ*.

It is important to note that the same block diagonalization due to symmetry that the above accomplishes for *Ĥ*_{inter} carries over to the *Ĥ* matrix when expressed in the |*κ*, *γ*⟩ representation. This is because *Ĥ* is invariant to the operations of *G*_{12}, the |*γ*⟩ are unaffected by those operations, and the |*κ*⟩ transform like *G*_{12} irreps.

#### 3. Iterative diagonalization procedure

We use the Chebyshev version^{31} of filter diagonalization^{32} to diagonalize *Ĥ*_{inter}. This method requires successive application of the operator on an initial, randomly generated state function. Such operation is straightforward for the kinetic-energy portion of *Ĥ*_{inter} [Eqs. (5) and (13)–(17)]. Matrix representations of many of the pieces of that operator in the primitive intermolecular basis can be evaluated analytically. The only exceptions are the pieces that involve derivatives with respect to *d*_{Z}. For those, we first evaluate analytically the operator matrix in the 1D HO-eigenfunction representation that is isomorphic to the *d*_{Z} DVR and then transform that matrix numerically to the DVR representation. The entire procedure is analogous to that described in Ref. 17, except that, here, there are some additional kinetic-energy-operator terms that arise from $T^vrWF(\omega ;q0)$.

To operate with the potential term in Eq. (5) [i.e., *V*_{inter}(**d**, *ω*)], we use a procedure very similar to that described in Kim *et al.*^{17} Briefly, we transform the state function from the |*α*, $v$ , *l*, *j*, *m*_{j}, *k*⟩ representation to the quadrature-grid representation |*α*, *ρ*_{β}, Φ_{γ}, *θ*_{q}, *ϕ*_{r}, *χ*_{s}⟩. We multiply each term in the resulting grid expansion of the state function by the value of *V*_{inter}(**d**, *ω*) at the associated grid point. Finally, we back-transform the result to the |*α*, $v$ , *l*, *j*, *m*_{j}, *k*⟩ basis. Details pertaining to the types of quadrature points and transformation matrices employed in this procedure are given in Ref. 17.

We make here two modifications to the earlier potential-energy-operation algorithm to increase its efficiency. First, to accomplish the (*m*_{j}, *k*) ↔ (*ϕ*, *χ*) transformations, we use 2D fast Fourier transforms. Second, we have reduced by a factor of six, without loss of accuracy, the number of (Φ_{γ}, *ϕ*_{r}) quadrature points in the grid representation. Specifically, while the full range of Fourier points for *ϕ* is included [*ϕ*_{r} = (*r* − 1)2*π*/*N*_{r}, where *r* = 1, …, *N*_{r}], only the first one-sixth of those for Φ are [Φ_{γ} = (*γ* − 1)2*π*/*N*_{γ}, *γ* = 1, … *N*_{γ}/6]. This is possible because of the sixfold symmetry of *Ĥ*_{inter} (and *Ĥ*) about *Ẑ* and because of the (*l* + *m*_{j})-symmetry factorization of the basis.^{33} The size of the 6D grid that we use is *N*_{α}*N*_{β}*N*_{γ}*N*_{q}*N*_{r}*N*_{s}/6 = 20(12)(24)(12)(24)(24)/6 = 6, 635, 520.

One last point pertaining to computational efficiency should be mentioned. The matrix elements of *Ĥ*_{inter} in the primitive basis are all real, even though the basis functions themselves are complex. This is straightforward to show for the kinetic-energy matrix elements from their analytic expressions. As to the potential-energy matrix elements, one can readily show that the symmetry transformations

and

require them to be real as well. In consequence, if the initial state function in filter diagonalization is chosen to have real expansion coefficients in the primitive basis, the analogous expansion coefficients of the state functions produced by successive operations of *Ĥ*_{inter} in the filter diagonalization procedure will remain real. Working with real state vectors allows for a decrease in computational effort by about a factor of two relative to working with complex ones, and we take advantage of this in the present work.

### E. Diagonalization of *Ĥ*_{intra}

To completely specify *Ĥ*_{intra}, the values of **d**_{0} and *ω*_{0} that define *V*_{adap} [see Eq. (8)] need to be fixed. For each isotopologue, we treat two cases. The first is the “bare-monomer” case corresponding to |**d**_{0}| = *∞* and $Vadap=VintraWF$. The second is the “dimer-adapted” case, for which we take **d**_{0} ≡ (*d*_{Z,0}, *ρ*_{0}, Φ_{0}) = (6.13°, 0.35°, 20°) (distances are in bohr) and *ω*_{0} ≡ (*ϕ*_{0}, *θ*_{0}, *χ*_{0}) = (20°, 30.1°, 180°). These values correspond to a rigid-body B–W interaction energy *V*_{B-W}(**d**_{0}, *ω*_{0}, **q**_{0}) = −1076 cm^{−1} relative to the energy of the infinitely separated, rigid monomers. This energy is very close to the minimum of *V*_{inter}.

Our procedure for computing the eigenfunctions of *Ĥ*_{intra} is very similar to that described in Ref. 22 in connection with a study of H_{2}O@C_{60}. Briefly, for a given species, we use a primitive basis consisting of the product of two tridiagonal-Morse DVRs^{34} covering the *R*_{1} and *R*_{2} coordinates (|*η*_{1}⟩ and |*η*_{2}⟩, each consisting of nine functions) and a potential-optimized DVR^{35} (|*η*_{3}⟩, consisting of 15 functions) covering the Θ coordinate. The basis used here for H_{2}O is identical to the intramolecular basis of Ref. 22. The basis used for HDO was constructed in the same way as that for H_{2}O (Sec. V A of Ref. 22) but with Radau mass factors appropriate to the HDO species.

Matrix elements of *V*_{adap} are diagonal in the above bases, with each nonzero element equal to that of the potential function evaluated at the relevant 3D quadrature point. Kinetic-energy (i.e., $T^vibWF$) matrix elements were evaluated as described in Sec. V B of Ref. 22. Filter diagonalization of *Ĥ*_{intra} was used to compute the intramolecular eigenvectors.

### F. Diagonalization of *Ĥ*

As mentioned in Sec. II D 2, the matrix of *Ĥ* is block diagonal in the |*κ*, *γ*⟩ basis, with the different blocks corresponding to the different *G*_{12} irreps to which the various |*κ*⟩ belong. In the case of benzene–H_{2}O, for all irreps, we set *N*_{inter} = 40 and *N*_{intra} = 34. For benzene–HDO, *N*_{inter} = 40 and *N*_{intra} = 42 were employed for all irreps. Characteristics of all of these *N*_{inter} = 40 computed intermolecular states are available (see the supplementary material). This resulted in the sizes of the 9D symmetry-factored bases that we have used, *N*_{inter} × *N*_{intra}, equal to 1360 for each *G*_{12} irrep of the H_{2}O complex and 1680 for each one of the HDO complex. Note that all of these bases are small enough to allow for direct diagonalization. The convergence of the calculated 9D levels with respect to *N*_{inter} was tested for the *A*_{1}, *E*_{1}, and *E*_{2} irreps of the HDO complex by doubling it from 40 to 80. In all intramolecular vibrational manifolds considered, for the intermolecular levels within 85 cm^{−1} relative to the respective “bare” intramolecular state, this changed their absolute energies by at most 0.2–0.3 cm^{−1} and energy differences by at most several tenths of cm^{−1}. The energies of the intramolecular vibrational excitations change by only a few hundredths of cm^{−1}. Concerning *N*_{intra}, its value of 34 for the H_{2}O complex was taken from our study of H_{2}O@C_{60},^{22} where it was tested for convergence. This set includes all eigenstates of *Ĥ*_{intra} with energies less than, or equal to, the energy of the state corresponding to three quanta in the asymmetric stretching (*ν*_{3}) mode. For the HDO complex, *N*_{intra} = 42 was chosen because it includes all eigenstates of *Ĥ*_{intra} with energies less than, or equal to, the state with three quanta in the OH stretching mode (*ν*_{OH}). The information on all of the computed intramolecular states that we have used to construct the 9D bases for the H_{2}O and HDO complexes, respectively, is provided (see the supplementary material). All such states correspond to the dimer-adapted potential.

The bulk of the computational effort in diagonalizing *Ĥ* is that of computing the matrix elements of $\Delta T^vrW$, $T^corW$, and Δ*V* [Eq. (12)]. Since $\Delta T^vrW$ and $T^corW$ depend on only six coordinates each, their matrix elements are considerably less challenging to compute than those of the 9D Δ*V*. We evaluate the former by using procedures completely analogous to those described in Sec. II E [Eqs. (27)–(30)] of Ref. 22. However, to evaluate the Δ*V* matrix elements, we introduce here a modification to that of our previous work that renders the computation a factor of approximately three more efficient.

As before [Eq. (27) of Ref. 22], we ultimately use

where |*η*⟩ ≡ |*η*_{1}*η*_{2}*η*_{3}⟩ is the 3D intramolecular (water) discrete variable representation (DVR). However, here, in order to compute the matrix elements on the rhs of Eq. (23), we work with the |*κ*⟩ in the 6D grid representation rather than in the primitive-basis representation. That is, at the start of the calculation, we compute and store all the relevant |*κ*⟩ in the |**Q**⟩ ≡ |*α*, *ρ*_{β}, Φ_{γ}, *θ*_{q}, *ϕ*_{r}, *χ*_{s}⟩ representation by transformation from the |*α*, $v$ , *l*, *j*, *m*_{j}, *k*⟩ primitive-basis representation (see Sec. II D 3). This is feasible because of the rather small number (i.e., 40) of |*κ*⟩ that needs to be included in any given symmetry block of the 9D basis. With the ⟨**Q**|*κ*⟩ in hand for all *κ*, it is straightforward to compute by quadrature the Δ*V* matrix elements on the rhs of Eq. (23) for a given *η* and all *κ* and *κ*′,

### G. Calculation of infrared vibrational-transition intensities

To compute IR transition intensities between the eigenstates of *Ĥ*, we make use of the general approach outlined in Refs. 36–38. That is, we compute the vibrational transition-dipole matrix elements between initial state |*ψ*_{i}⟩ and final state |*ψ*_{f}⟩ with the dipole spherical-tensor components, $\mu KECF$ (*K* = −1, 0, +1), referred to an Eckart complex-fixed (ECF) frame,

We then use these to estimate transition intensities as

To compute the $\u27e8\mu KECF\u27e9if$, we need to transform transition-dipole components from the CF frame that we use to compute the 9D eigenfunctions (see Sec. II A) to the ECF frame. This is accomplished by the following equation:

where the $\mu MCF$ are spherical-tensor dipole-component operators referred to the CF frame, the $DM,K(1)(\Omega )$ are Wigner rotation matrix elements, and Ω denotes the Euler angles that describe the rotation of the CF frame into the ECF frame.

Now, in principle, the Ω and the $\mu MCF$ depend on all nine of the coordinates (**d**, *ω*, **q**). Hence, the matrix elements on the rhs of Eq. (27) are 9D integrals. However, two approximations reduce considerably the effort required to compute them. The first is the assumption that the Ω are independent of the water intramolecular coordinates **q**. The rationale for this approximation derives from the semi-rigid nature of the water moiety. The small excursions of its nuclei away from its equilibrium geometry should not appreciably affect the orientation of the ECF frame relative to the CF frame. With this approximation,

where we have expressed |*ψ*_{f}⟩ and |*ψ*_{i}⟩ in terms of the |*κ*, *γ*⟩ basis. The second approximation is to assume that the transition-dipole magnitudes are those of the |*γ*⟩ → |*γ*′⟩ *monomer* transitions and that, consequently, their CF-frame components vary only by virtue of the variation in the orientation of the WF frame with respect to the CF frame. This leads to

where the $\mu K\u2032W$ are dipole-operator spherical-tensor components referred to the WF frame. Equation (28) can then be written as

One sees that the main task in evaluating these matrix elements is the computation of 6D integrals over the intermolecular coordinates **Q** on which the |*κ*⟩ depend. We accomplish this by transforming the |*κ*⟩ to the grid representation described above in Secs. II D 3 and II F. We also compute and store the $DM,K(1)(\Omega )$ and $DM,K\u2032(1)(\omega )$ on that grid. To compute the former, we make use of Ref. 39 to determine how the Ω depend on **Q**. In this task, the reference geometry of the complex is taken to be that of benzene plus a mass equal to that of the water moiety located at *ρ* = 0 and *Z* = 6.2 bohr. The $DM,K\u2032(1)(\omega )$ are easily obtained since the *ω* comprise part of the 6D grid. The matrix elements on the rhs of Eq. (30) are then computed by quadrature,

We are interested here in modeling IR spectra corresponding to the OH-stretching fundamentals of the complexes. Aside from having the relevant 9D eigenvectors in hand, the other information that is required to accomplish this is monomer transition-dipole matrix elements. For the H_{2}O complex, the transition dipole for the *ν*_{1} fundamental is along the WF *ẑ* axis and that for the *ν*_{3} fundamental is along the WF $x^$ axis. In isolated H_{2}O, the relative magnitudes of these dipole vectors are in the ratio $|\mu \nu 1|/|\mu \nu 3|\u22431/20$.^{40} In the benzene–H_{2}O complex, there is some evidence that this ratio is considerably larger ($2/3$).^{6} We have computed spectra for both of these relative values, but the spectra that we present below correspond to the latter one. In the case of the HDO complex, only a single monomer transition-dipole vector is relevant to the OH stretching region. We take that vector to be aligned along the OH bond.

### H. Calculation of Raman spectra involving transitions between intermolecular vibrational states

We also model herein the Raman spectrum of the H_{2}O complex in the low frequency (0 cm^{−1}–100 cm^{−1}) region that corresponds to transitions between intermolecular vibrational states within the intramolecular vibrational ground-state manifold of the species. In this endeavor, as in the calculation of the IR spectra, we make use of the general approach described in Refs. 36–38, although the relevant operator now is the polarizability tensor rather than the electric-dipole vector. In particular, we compute the Raman intensity of the transition from vibrational state |*ψ*_{i}⟩ to vibrational state |*ψ*_{f}⟩ as

where

and the $[\alpha K(2)]ECF$ are spherical-tensor operators corresponding to the anisotropic part of the complex’s polarizability tensor referred to the Eckart complex-fixed frame. We further make the assumption that this polarizability tensor is contributed to only by the permanent polarizabilities of the monomer moieties in their ground vibronic states. With this assumption, the polarizability tensor of the complex referred to the CF frame is given by

where $[\alpha M(2)]B$ is the anisotropic polarizability spherical tensor of benzene in its ground state with components referred to the CF (i.e., benzene-fixed) frame and $[\alpha q(2)]W$ is that of water in its ground state with components referred to the WF frame. Transforming Eq. (34) to the Eckart frame and evaluating Eq. (33), one has

The matrix elements on the rhs of Eq. (35) can be evaluated in a manner similar to the evaluation of Eqs. (30) and (31).

For the computed Raman spectra reported below, we have used the polarizability tensors of benzene and water reported in Refs. 41 and 42, respectively. Notably, the benzene polarizability dominates in determining the Raman intensities; the anisotropy of that tensor is roughly two orders of magnitude larger than that of water moiety.

## III. RESULTS AND DISCUSSION

### A. Intermolecular vibrational eigenstates

Selected low-energy intermolecular eigenstates of benzene–H_{2}O and benzene–HDO, calculated in both 6D (rigid-monomer) and 9D, for water in the ground intramolecular vibrational state, with the molecular parameters listed in Table I, are presented in Tables II and III, respectively. Inspection of the two tables shows that for both complexes, the 9D level energies are slightly larger than their 6D counterparts by less than 2 cm^{−1} for the H_{2}O complex and by less than 3 cm^{−1} for the HDO complex (with the exception of the state *N* = 4 in Table III, where the difference is 4.9 cm^{−1}). In addition, not surprisingly, the results for the H_{2}O complex are very similar to those of the earlier rigid-monomer study;^{17} any differences between the two sets can be attributed to slightly different monomer geometries, inertial parameters, and basis sets. Among their salient features is a vibrationally averaged dimer geometry in which the c.m. of the water moiety is about 6.2 bohr from the center of the benzene ring and along the benzene C_{6} axis, with the H atoms of the water moiety directed toward the benzene (*π* hydrogen bonding).

. | A_{⊥} (cm^{−1})
. | μ_{d}
. | μ_{1}
. | μ_{2}
. | R_{1,0}
. | R_{2,0}
. | cos(Θ_{0})
. |
---|---|---|---|---|---|---|---|

H_{2}O | 0.189 08 | 26 680.3725 | 1731.9320 | 1731.9320 | 1.828 367 4 | 1.828 367 4 | −0.305 438 61 |

HDO | 0.189 08 | 27 878.6684 | 1736.0869 | 3276.2121 | 1.842 675 6 | 1.828 160 9 | −0.329 414 95 |

. | A_{⊥} (cm^{−1})
. | μ_{d}
. | μ_{1}
. | μ_{2}
. | R_{1,0}
. | R_{2,0}
. | cos(Θ_{0})
. |
---|---|---|---|---|---|---|---|

H_{2}O | 0.189 08 | 26 680.3725 | 1731.9320 | 1731.9320 | 1.828 367 4 | 1.828 367 4 | −0.305 438 61 |

HDO | 0.189 08 | 27 878.6684 | 1736.0869 | 3276.2121 | 1.842 675 6 | 1.828 160 9 | −0.329 414 95 |

. | Para
. | Ortho
. | . | . | Basis-state . |
---|---|---|---|---|---|

N
. | ΔE: 6D/9D
. | ΔE: 6D/9D
. | Irrep . | Assignment . | norm . |

1 | 0.00/0.00 | $A1+$ | … | 0.999 | |

2 | 16.32/16.53 | $E1\u2212$ | m(1) | 0.999 | |

3 | 35.36/36.07 | $E1+$ | ν_{XY} | 0.999 | |

4 | 41.14/43.09 | $A1\u2212$ | ν_{rock} | 0.998 | |

6 | 53.22/54.21 | $E2\u2212$ | ν_{XY} + m(1) | 0.999 | |

7 | 61.24/62.13 | $A1+$ | 2ν_{XY} | 0.998 | |

8 | 64.86/65.72 | $E2+$ | m(2) | 0.996 | |

10 | 69.42/71.27 | $E2+$ | 2ν_{XY} | 0.996 | |

13 | 77.87/79.17 | $A1+$ | ν_{Z} | 0.997 |

. | Para
. | Ortho
. | . | . | Basis-state . |
---|---|---|---|---|---|

N
. | ΔE: 6D/9D
. | ΔE: 6D/9D
. | Irrep . | Assignment . | norm . |

1 | 0.00/0.00 | $A1+$ | … | 0.999 | |

2 | 16.32/16.53 | $E1\u2212$ | m(1) | 0.999 | |

3 | 35.36/36.07 | $E1+$ | ν_{XY} | 0.999 | |

4 | 41.14/43.09 | $A1\u2212$ | ν_{rock} | 0.998 | |

6 | 53.22/54.21 | $E2\u2212$ | ν_{XY} + m(1) | 0.999 | |

7 | 61.24/62.13 | $A1+$ | 2ν_{XY} | 0.998 | |

8 | 64.86/65.72 | $E2+$ | m(2) | 0.996 | |

10 | 69.42/71.27 | $E2+$ | 2ν_{XY} | 0.996 | |

13 | 77.87/79.17 | $A1+$ | ν_{Z} | 0.997 |

N
. | ΔE: 6D/9D
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00/0.00 | A_{1} | … | 0.991 |

2 | 12.59/13.27 | E_{1} | m(1) | 0.991 |

3 | 35.01/36.42 | E_{1} | ν_{XY} | 0.978 |

4 | 36.38/41.28 | A_{1} | ν_{rock} | 0.971 |

5 | 47.00/48.25 | E_{2} | ν_{XY} + m(1) | 0.964 |

7 | 50.41/52.93 | E_{2} | m(2) | 0.953 |

8 | 56.46/58.16 | A_{1} | 2ν_{XY} | 0.971 |

11 | 68.98/71.78 | E_{2} | 2ν_{XY} | 0.971 |

16 | 77.69/79.68 | A_{1} | ν_{Z} | 0.964 |

N
. | ΔE: 6D/9D
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00/0.00 | A_{1} | … | 0.991 |

2 | 12.59/13.27 | E_{1} | m(1) | 0.991 |

3 | 35.01/36.42 | E_{1} | ν_{XY} | 0.978 |

4 | 36.38/41.28 | A_{1} | ν_{rock} | 0.971 |

5 | 47.00/48.25 | E_{2} | ν_{XY} + m(1) | 0.964 |

7 | 50.41/52.93 | E_{2} | m(2) | 0.953 |

8 | 56.46/58.16 | A_{1} | 2ν_{XY} | 0.971 |

11 | 68.98/71.78 | E_{2} | 2ν_{XY} | 0.971 |

16 | 77.69/79.68 | A_{1} | ν_{Z} | 0.964 |

The energy levels in Tables II and III are associated with the fundamentals, overtones, and combinations of four intermolecular modes, which are as follows.

(1) the internal-rotation (about the benzene *C*_{6} axis) mode, with internal-rotation constants *B* ≃ 16/13 cm^{−1} (H_{2}O/HDO) and with the internal-rotation energy levels given by *E*_{m} = *Bm*^{2} (*m* = 0, ±1, ±2, …) [we designate this mode as *m*(0) corresponding to *m* = 0, *m*(1) corresponding to *m* = ±1, *m*(2) corresponding to *m* = ±2, etc.], (2) a doubly degenerate in-plane translational (or “van der Waals bending”) mode, *ν*_{XY}, at about 35 cm^{−1}, (3) a nondegenerate rocking (or “in-plane torsional”) mode, *ν*_{rock}, at about 41 cm^{−1}, which involves the hindered rotation of the water about the out-of-plane (i.e., *c*) principal axis of that moiety, and (4) a nondegenerate out-of-plane translational (or “van der Waals stretching”) mode, *ν*_{Z}, at about 80 cm^{−1}.

Comparison of Table III with Table II indicates that the intermolecular results for the HDO complex are similar to those for the H_{2}O species. However, there are two significant qualitative differences in the two sets of results. These relate to vibrationally averaged geometries and the nature of the “rocking” vibration.

As to geometry, like the H_{2}O species, the HDO ground state has ⟨*d*_{Z}⟩ ≃ 6.2 bohr and ⟨*d*_{X}⟩ = ⟨*d*_{Y}⟩ = 0 bohr. However, it is only the D nucleus that interacts closely with the benzene moiety in the HDO species, whereas (by symmetry) both H nuclei interact equivalently with the benzene in the H_{2}O species. This point is illustrated graphically by the plots shown in Figs. 1(a) and 1(b). These are contour plots of the two-point probability-density functions,

where the angle brackets on the rhs of Eq. (36) denote an expectation value and $Z1\u2032$ and $Z2\u2032$ represent the *Z* (CF frame) coordinates of the two hydrogen nuclei of the water moiety. Figure 1(a) is a plot of *P*(*Z*_{1}, *Z*_{2}) for the *Ĥ*_{inter} ground state of the H_{2}O species. As one expects, it is symmetric with respect to nuclear interchange and is such that when one H nucleus is close to the benzene ring plane, the other is further away. The analogous plot for the HDO complex [Fig. 1(b)], however, reveals a marked bias for the D nucleus to be localized close to the ring plane and for the H nucleus to be simultaneously further away. Moreover, this bias characterizes the wavefunctions for most of the elementary intermolecular excitations in the HDO species. In particular, the lowest-energy internal rotation (*m*), *ν*_{XY}, and *ν*_{Z} excited states in the HDO complex all have the D and H nuclei positioned in a manner similar to that reflected in Fig. 1(b). The *ν*_{rock} excitation, however, reverses this bias in the HDO species. Figures 1(c) and 1(d) show *P*(*Z*_{1}, *Z*_{2}) plots for the intermolecular states that we assign to *ν*_{rock} for the H_{2}O and HDO complexes, respectively. These two plots both exhibit the expected presence (for *ν*_{rock}) of a nodal surface that separates geometries in which one hydrogen nucleus is close to the benzene and the other is further away from those in which those positions are reversed. Unlike the (necessarily) symmetrical plot for the H_{2}O complex, the plot for the HDO isotopologue is asymmetrical and is such that the H nucleus has a greater probability of being close to the benzene than does the D nucleus.

### B. Intramolecular vibrational eigenstates of H_{2}O and HDO in the dimers

Table IV displays results, computational and experimental (where available), that pertain to the intramolecular vibrational excitations of H_{2}O and HDO in the complexes with benzene and their frequency shifts caused by the complexation. Columns 2 and 3 show the four lowest-energy intramolecular vibrational states of the H_{2}O and HDO moieties from 3D quantum calculations for $VintraWF$, the 3D PES of the isolated H_{2}O,^{28} and *V*_{adap}, the dimer-adapted 3D H_{2}O PES, respectively. In the case of isolated H_{2}O, the shown computed level energies agree to within a tenth of a cm^{−1} with the measured ones.^{18} However, for HDO, the calculated intramolecular levels differ from the experimental values^{18} by +2.7 cm^{−1} for *ν*_{2}, +1.5 cm^{−1} for *ν*_{OD}, and −3.3 cm^{−1} for *ν*_{OH}. This is in all likelihood due to the fact that the 3D PES employed was refined, or optimized, starting from an *ab initio* PES, to reproduce with exceptional accuracy the measured vibrational levels of H_{2}O over a broad range of energies.^{28} This optimization probably made the PES slightly less accurate for HDO.

PES: . | $VintraWF$ (3D) . | V_{adap} (3D)
. | Δν^{3D}
. | V_{tot} (9D)
. | Δν^{9D}
. | Expt. . | Δν^{expt}
. |
---|---|---|---|---|---|---|---|

H_{2}O | |||||||

ν_{2} | 1594.74 | 1610.89 | +16.15 | 1615.16 | +20.42 | ||

2ν_{2} | 3151.62 | 3186.55 | +34.93 | 3196.62 | +45.00 | ||

ν_{1} | 3656.97 | 3645.95 | −11.02 | 3645.16 | −11.81 | 3634 | −23 |

ν_{3} | 3755.84 | 3733.77 | −22.07 | 3739.39 | −16.45 | 3731 | −25 |

HDO | |||||||

ν_{2} | 1406.25 | 1420.63 | +14.38 | 1422.84 | +16.59 | ||

ν_{OD} | 2725.19 | 2714.36 | −10.83 | 2707.09 | −18.10 | ||

2ν_{2} | 2788.57 | 2817.19 | +28.62 | 2822.32 | +33.75 | ||

ν_{OH} | 3704.21 | 3686.64 | −17.57 | 3694.73 | −9.48 | 3683 | −24 |

PES: . | $VintraWF$ (3D) . | V_{adap} (3D)
. | Δν^{3D}
. | V_{tot} (9D)
. | Δν^{9D}
. | Expt. . | Δν^{expt}
. |
---|---|---|---|---|---|---|---|

H_{2}O | |||||||

ν_{2} | 1594.74 | 1610.89 | +16.15 | 1615.16 | +20.42 | ||

2ν_{2} | 3151.62 | 3186.55 | +34.93 | 3196.62 | +45.00 | ||

ν_{1} | 3656.97 | 3645.95 | −11.02 | 3645.16 | −11.81 | 3634 | −23 |

ν_{3} | 3755.84 | 3733.77 | −22.07 | 3739.39 | −16.45 | 3731 | −25 |

HDO | |||||||

ν_{2} | 1406.25 | 1420.63 | +14.38 | 1422.84 | +16.59 | ||

ν_{OD} | 2725.19 | 2714.36 | −10.83 | 2707.09 | −18.10 | ||

2ν_{2} | 2788.57 | 2817.19 | +28.62 | 2822.32 | +33.75 | ||

ν_{OH} | 3704.21 | 3686.64 | −17.57 | 3694.73 | −9.48 | 3683 | −24 |

Column 4 gives the frequency shifts (Δ*ν*^{3D}) of the intramolecular states of H_{2}O and HDO in going from the bare-monomer to the dimer-adapted potential, computed in 3D. One notes the appreciable shifts obtained when the influence of the intermolecular interaction on the intramolecular states is included, if only approximately, by using the dimer-adapted potential. The stretching modes of H_{2}O (*ν*_{1} and *ν*_{3}) and HDO (*ν*_{OD} and *ν*_{OH}) are redshifted from the ones calculated for the respective isolated monomers; the redshifts are similar in magnitude for H_{2}O and HDO. In contrast, the H_{2}O/HDO bend fundamentals (*ν*_{2}) and overtones (2*ν*_{2}) are significantly blueshifted; again, these blueshifts are similar for H_{2}O and HDO.

Columns 5 and 6 of Table IV show the results of the quantum 9D calculations pertaining to the four lowest-energy intramolecular vibrational states of H_{2}O and HDO in the dimers and their frequency shifts (Δ*ν*^{9D}), respectively. The 9D frequency shifts are, in general, appreciably larger in magnitude than their 3D counterparts. The same trend was observed for H_{2}O@C_{60}.^{22} The exception to this is the *ν*_{3} mode of H_{2}O and the *ν*_{OH} mode of HDO; their 9D values are 25% and 45% *smaller* by magnitude, respectively, than the 3D results. Combined with the observation that for *ν*_{OD}, the magnitude of the 9D redshift is 68% *larger* than the 3D one, it is evident that increasing the dimensionality of the treatment from 3D to 9D has a particularly large effect on the *ν*_{OH} and *ν*_{OD} stretch modes of HDO, but in opposite directions regarding their frequency shifts.

The main difference between the 3D quantum calculations on the benzene-adapted H_{2}O/HDO PES and the 9D quantum calculations is that the former are performed for the fixed position and orientation of H_{2}O/HDO relative to benzene and consider explicitly only the three intramolecular DOFs of the water moiety, while the latter treat both the intra- and intermolecular DOFs of H_{2}O/HDO in full dimensionality and as fully coupled. Therefore, comparison of the 3D and 9D results provides insight into the impact of the LAM intermolecular vibrations, and averaging over them, on the magnitude of the intramolecular frequency shifts of H_{2}O/HDO. Clearly, any quantitative treatment of the latter has to take LAMs into account rigorously.

Finally, the last two columns of Table IV present the experimental intramolecular vibrational levels for H_{2}O and HDO complexes and their frequency shifts (Δ*ν*^{expt}), respectively. Since for both dimers, the IR spectra were measured in the OH-stretch fundamental region,^{6} experimental data are available only for the *ν*_{1} and *ν*_{3} modes of H_{2}O and the *ν*_{OH} mode of HDO. For both H_{2}O and HDO, the comparison of the frequency shifts measured for H_{2}O and HDO (Δ*ν*^{expt}) with those from the 9D quantum calculations (Δ*ν*^{9D}) reveals a semiquantitative agreement. The theory reproduces the observed direction of the frequency shifts, i.e., redshifts, and the fact that the *ν*_{3} redshift is somewhat larger by absolute value than that for the *ν*_{1} mode. However, the magnitudes of the 9D-calculated redshifts for H_{2}O and HDO are significantly smaller, by roughly a factor of two, than those of the corresponding measured values. This points to a deficiency in the 9D PES employed that it underestimates the strength of the coupling between the intramolecular vibrations of H_{2}O and HDO and the intermolecular DOFs of the two dimers.

### C. Effects of H_{2}O/HDO intramolecular excitations on the intermolecular vibrational levels

Tables V–VIII show the set of 9D-calculated low-energy intermolecular levels of benzene–H_{2}O, which appears in Table II, but now in the *ν*_{2}, 2*ν*_{2}, *ν*_{1}, and *ν*_{3} intramolecular manifolds, respectively. Likewise, Tables IX–XII display the set of 9D-calculated low-energy intermolecular levels of benzene–HDO appearing in Table III, but in the *ν*_{2}, *ν*_{OD}, 2*ν*_{2}, and *ν*_{OH} intramolecular manifolds, respectively. By comparing these results with those computed in 9D for the ground intramolecular vibrational states of H_{2}O and HDO complexes, respectively, in Tables II and Table III, we can discern how, and to what degree, excitation of different intramolecular modes affects the intermolecular vibrational states of the dimers.

N
. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.998 | |

2 | 16.87 | $E1\u2212$ | m(1) | 0.998 | |

3 | 35.69 | $E1+$ | ν_{XY} | 0.998 | |

4 | 42.11 | $A1\u2212$ | ν_{rock} | 0.997 | |

6 | 54.18 | $E2\u2212$ | ν_{XY} + m(1) | 0.999 | |

7 | 61.43 | $A1+$ | 2ν_{XY} | 0.997 | |

8 | 66.89 | $E2+$ | m(2) | 0.980 | |

10 | 70.52 | $E2+$ | 2ν_{XY} | 0.979 | |

13 | 78.35 | $A1+$ | ν_{Z} | 0.997 |

N
. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.998 | |

2 | 16.87 | $E1\u2212$ | m(1) | 0.998 | |

3 | 35.69 | $E1+$ | ν_{XY} | 0.998 | |

4 | 42.11 | $A1\u2212$ | ν_{rock} | 0.997 | |

6 | 54.18 | $E2\u2212$ | ν_{XY} + m(1) | 0.999 | |

7 | 61.43 | $A1+$ | 2ν_{XY} | 0.997 | |

8 | 66.89 | $E2+$ | m(2) | 0.980 | |

10 | 70.52 | $E2+$ | 2ν_{XY} | 0.979 | |

13 | 78.35 | $A1+$ | ν_{Z} | 0.997 |

N
. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.996 | |

2 | 17.29 | $E1\u2212$ | m(1) | 0.997 | |

3 | 35.21 | $E1+$ | ν_{XY} | 0.996 | |

4 | 40.68 | $A1\u2212$ | ν_{rock} | 0.992 | |

6 | 54.18 | $E2\u2212$ | ν_{XY} + m(1) | 0.996 | |

7 | 60.58 | $A1+$ | 2ν_{XY} | 0.993 | |

9 | 68.05 | $E2+$ | m(2) | 0.804 | |

11 | 69.85 | $E2+$ | 2ν_{XY} | 0.803 | |

13 | 77.40 | $A1+$ | ν_{Z} | 0.993 |

N
. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.996 | |

2 | 17.29 | $E1\u2212$ | m(1) | 0.997 | |

3 | 35.21 | $E1+$ | ν_{XY} | 0.996 | |

4 | 40.68 | $A1\u2212$ | ν_{rock} | 0.992 | |

6 | 54.18 | $E2\u2212$ | ν_{XY} + m(1) | 0.996 | |

7 | 60.58 | $A1+$ | 2ν_{XY} | 0.993 | |

9 | 68.05 | $E2+$ | m(2) | 0.804 | |

11 | 69.85 | $E2+$ | 2ν_{XY} | 0.803 | |

13 | 77.40 | $A1+$ | ν_{Z} | 0.993 |

. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.990 | |

2 | 16.39 | $E1\u2212$ | m(1) | 0.991 | |

3 | 36.02 | $E1+$ | ν_{XY} | 0.988 | |

4 | 40.21 | $A1\u2212$ | ν_{rock} | 0.942 | |

6 | 54.19 | $E2\u2212$ | ν_{XY} + m(1) | 0.990 | |

7 | 61.89 | $A1+$ | 2ν_{XY} | 0.986 | |

9 | 65.13 | $E2+$ | m(2) | 0.990 | |

11 | 71.18 | $E2+$ | 2ν_{XY} | 0.984 | |

13 | 79.59 | $A1+$ | ν_{Z} | 0.988 |

. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1+$ | … | 0.990 | |

2 | 16.39 | $E1\u2212$ | m(1) | 0.991 | |

3 | 36.02 | $E1+$ | ν_{XY} | 0.988 | |

4 | 40.21 | $A1\u2212$ | ν_{rock} | 0.942 | |

6 | 54.19 | $E2\u2212$ | ν_{XY} + m(1) | 0.990 | |

7 | 61.89 | $A1+$ | 2ν_{XY} | 0.986 | |

9 | 65.13 | $E2+$ | m(2) | 0.990 | |

11 | 71.18 | $E2+$ | 2ν_{XY} | 0.984 | |

13 | 79.59 | $A1+$ | ν_{Z} | 0.988 |

. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1\u2212$ | … | 0.859 | |

2 | 17.35 | $E1+$ | m(1) | 0.670 | |

3 | 36.68 | $A1+$ | ν_{rock} | 0.571 | |

4 | 36.75 | $E1\u2212$ | ν_{XY} | 0.782 | |

5 | 50.73 | $E2+$ | ν_{XY} + m(1) | 0.418 | |

7 | 63.91 | $E2\u2212$ | m(2) | 0.686 | |

8 | 63.40 | $A1\u2212$ | 2ν_{XY} | 0.603 | |

10 | 72.96 | $E2\u2212$ | 2ν_{XY} | 0.761 | |

12 | 79.58 | $A1\u2212$ | ν_{Z} | 0.578 |

. | para ΔE
. | ortho ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|---|

1 | 0.00 | $A1\u2212$ | … | 0.859 | |

2 | 17.35 | $E1+$ | m(1) | 0.670 | |

3 | 36.68 | $A1+$ | ν_{rock} | 0.571 | |

4 | 36.75 | $E1\u2212$ | ν_{XY} | 0.782 | |

5 | 50.73 | $E2+$ | ν_{XY} + m(1) | 0.418 | |

7 | 63.91 | $E2\u2212$ | m(2) | 0.686 | |

8 | 63.40 | $A1\u2212$ | 2ν_{XY} | 0.603 | |

10 | 72.96 | $E2\u2212$ | 2ν_{XY} | 0.761 | |

12 | 79.58 | $A1\u2212$ | ν_{Z} | 0.578 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.986 |

2 | 13.68 | E_{1} | m(1) | 0.986 |

3 | 36.33 | E_{1} | ν_{XY} | 0.969 |

4 | 41.57 | A_{1} | ν_{rock} | 0.959 |

5 | 48.64 | E_{2} | ν_{XY} + m(1) | 0.919 |

7 | 54.00 | E_{2} | m(2) | 0.900 |

8 | 57.82 | A_{1} | 2ν_{XY} | 0.963 |

11 | 71.64 | E_{2} | 2ν_{XY} | 0.960 |

16 | 79.12 | A_{1} | ν_{Z} | 0.953 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.986 |

2 | 13.68 | E_{1} | m(1) | 0.986 |

3 | 36.33 | E_{1} | ν_{XY} | 0.969 |

4 | 41.57 | A_{1} | ν_{rock} | 0.959 |

5 | 48.64 | E_{2} | ν_{XY} + m(1) | 0.919 |

7 | 54.00 | E_{2} | m(2) | 0.900 |

8 | 57.82 | A_{1} | 2ν_{XY} | 0.963 |

11 | 71.64 | E_{2} | 2ν_{XY} | 0.960 |

16 | 79.12 | A_{1} | ν_{Z} | 0.953 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.950 |

2 | 14.13 | E_{1} | m(1) | 0.944 |

3 | 37.40 | E_{1} | ν_{XY} | 0.894 |

4 | 46.89 | A_{1} | ν_{rock} | 0.727 |

5 | 49.37 | E_{2} | ν_{XY} + m(1) | 0.849 |

7 | 56.24 | E_{2} | m(2) | 0.793 |

8 | 61.97 | A_{1} | 2ν_{XY} | 0.705 |

11 | 73.87 | E_{2} | 2ν_{XY} | 0.835 |

16 | 80.79 | A_{1} | ν_{Z} | 0.861 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.950 |

2 | 14.13 | E_{1} | m(1) | 0.944 |

3 | 37.40 | E_{1} | ν_{XY} | 0.894 |

4 | 46.89 | A_{1} | ν_{rock} | 0.727 |

5 | 49.37 | E_{2} | ν_{XY} + m(1) | 0.849 |

7 | 56.24 | E_{2} | m(2) | 0.793 |

8 | 61.97 | A_{1} | 2ν_{XY} | 0.705 |

11 | 73.87 | E_{2} | 2ν_{XY} | 0.835 |

16 | 80.79 | A_{1} | ν_{Z} | 0.861 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.976 |

2 | 14.09 | E_{1} | m(1) | 0.896 |

3 | 36.22 | E_{1} | ν_{XY} | 0.950 |

4 | 41.45 | A_{1} | ν_{rock} | 0.925 |

5 | 48.96 | E_{2} | ν_{XY} + m(1) | 0.808 |

7 | 55.22 | E_{2} | m(2) | 0.685 |

8 | 57.25 | A_{1} | 2ν_{XY} | 0.957 |

12 | 71.46 | E_{2} | 2ν_{XY} | 0.932 |

16 | 78.55 | A_{1} | ν_{Z} | 0.927 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.976 |

2 | 14.09 | E_{1} | m(1) | 0.896 |

3 | 36.22 | E_{1} | ν_{XY} | 0.950 |

4 | 41.45 | A_{1} | ν_{rock} | 0.925 |

5 | 48.96 | E_{2} | ν_{XY} + m(1) | 0.808 |

7 | 55.22 | E_{2} | m(2) | 0.685 |

8 | 57.25 | A_{1} | 2ν_{XY} | 0.957 |

12 | 71.46 | E_{2} | 2ν_{XY} | 0.932 |

16 | 78.55 | A_{1} | ν_{Z} | 0.927 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.955 |

2 | 11.66 | E_{1} | m(1) | 0.953 |

3 | 33.34 | A_{1} | ν_{rock} | 0.937 |

4 | 35.02 | E_{1} | ν_{XY} | 0.935 |

5 | 44.87 | E_{2} | ν_{XY} + m(1) | 0.668 |

7 | 49.07 | E_{2} | m(2) | 0.655 |

9 | 58.04 | A_{1} | 2ν_{XY} | 0.909 |

12 | 69.20 | E_{2} | 2ν_{XY} | 0.919 |

17 | 79.36 | A_{1} | ν_{Z} | 0.946 |

N
. | ΔE
. | Irrep . | Assignment . | Basis-state norm . |
---|---|---|---|---|

1 | 0.00 | A_{1} | … | 0.955 |

2 | 11.66 | E_{1} | m(1) | 0.953 |

3 | 33.34 | A_{1} | ν_{rock} | 0.937 |

4 | 35.02 | E_{1} | ν_{XY} | 0.935 |

5 | 44.87 | E_{2} | ν_{XY} + m(1) | 0.668 |

7 | 49.07 | E_{2} | m(2) | 0.655 |

9 | 58.04 | A_{1} | 2ν_{XY} | 0.909 |

12 | 69.20 | E_{2} | 2ν_{XY} | 0.919 |

17 | 79.36 | A_{1} | ν_{Z} | 0.946 |

The four intermolecular vibrational modes that fall in the energy range considered up to 85 cm^{−1}, *m*, *ν*_{XY}, *ν*_{rock}, and *ν*_{Z}, can be, roughly speaking, divided into two groups: (*i*) those that primarily involve the motion of the c.m. of H_{2}O relative to the benzene plane—*ν*_{XY} and *ν*_{Z}, and (*ii*) those that primarily involve the LAM motions of the H/D atoms of H_{2}O and HDO—the internal rotation mode (*m*) and the rocking (or “in-plane torsional”) mode (*ν*_{rock}). Inspection of the results in Tables V–XII readily shows that, for both H_{2}O and HDO, the two intermolecular modes in the first group, *ν*_{XY} and *ν*_{Z}, change only slightly, by 1–2 cm^{−1}, from their respective values in the ground intramolecular vibrational state upon the excitation of any of the intramolecular modes.

In contrast, a mode in the second group, *ν*_{rock}, exhibits pronounced sensitivity to the excitation of the stretching modes of H_{2}O and HDO. For H_{2}O, in the *ν*_{3} manifold, its energy, 36.7 cm^{−1}, is 6.4 cm^{−1} (or 14.9%) lower than that in the ground intramolecular state, 43.1 cm^{−1}. Excitation of the *ν*_{1} mode has a smaller effect on the energy of *ν*_{rock}, lowering it by 2.9 cm^{−1} (or 7.3%) to 40.2 cm^{−1}. Changes are even larger and more intricate in the case of HDO. In the *ν*_{OH} manifold, the energy of *ν*_{rock}, 33.3 cm^{−1}, is 8.0 cm^{−1} (or 19%) *lower* than that in the ground intramolecular state, 41.3 cm^{−1}. Interestingly, excitation of *ν*_{OD} has the *opposite* effect on the energy of *ν*_{rock}, *increasing* it by 5.6 cm^{−1} (or 11.4%) relative to the value in the ground intramolecular state.

The other mode in the second group, *m*, is less affected by the excitation of H_{2}O and HDO stretching modes. For H_{2}O, both *ν*_{1} and *ν*_{3} excitations change the energy of its fundamental, *m*(1), by less than one cm^{−1} from the value it has in the intramolecular vibrational ground state, 16.5 cm^{−1}. However, for HDO, in the *ν*_{OH} manifold, the energy of *m*(1), 11.7 cm^{−1}, is 1.6 cm^{−1} (or 12%) lower than that in the ground intramolecular state, 13.3 cm^{−1}. The excitation of *ν*_{OD} has the opposite effect, increasing the *m*(1) energy by 0.8 cm^{−1} (or 6%) relative to the value in the ground intramolecular state.

Thus, in the HDO complex, excitations of the *ν*_{OH} and *ν*_{OD} modes have appreciable but *opposite* effects on the fundamentals of the *m* and *ν*_{rock} modes. In the *ν*_{OH} intramolecular manifold, the energies of both *ν*_{rock} and *m*(1) fundamentals are *lower* than their corresponding values in the ground intramolecular state. However, exciting the *ν*_{OD} mode increases the energies of the *ν*_{rock} and *m*(1) fundamentals relative to those in the ground intramolecular state.

Interestingly, the excitation of the intramolecular bend (*ν*_{2}) mode has a very small effect, at most a cm^{−1}, on any of the four intermolecular modes in both the H_{2}O and HDO complexes.

Finally, we turn our attention to the last column in Tables II, III, and V–XII, labeled “Basis-state norm” (BSN). The entries in this column measure the contribution that the dominant product inter/intra-basis state |*κ*, *γ*⟩ ≡ |*κ*⟩|*γ*⟩ makes to a given 9D eigenstate. For benzene–H_{2}O, all the low-lying intermolecular eigenstates in the *ν*_{2} and *ν*_{1} intramolecular vibrational manifolds (Tables V and VII) have BSNs greater than 0.9, meaning that these eigenstates are highly “pure,” i.e., totally dominated by a single |*κ*⟩|*γ*⟩ basis state. However, in Table VI, two states in the 2*ν*_{2} manifold, *N* = 9 and *N* = 11, assigned as *m*(2) and 2*ν*_{XY}, respectively, have BSNs of about 0.80. This is due to basis-state mixing arising from the near-degeneracy of these two intermolecular excitations within the 2*ν*_{2} manifold. Much more extensive basis-state mixing can be seen in Table VIII, for the *ν*_{3} manifold, where all of the intermolecular states shown have BSNs ranging from 0.42 to 0.86, indicating that they have strongly mixed character. An examination of the other 9D basis states that contribute significantly to these mixed eigenstates reveals that they are predominantly excited intermolecular states in the *ν*_{1} manifold. This significant inter-manifold coupling is possible because the origin of the *ν*_{1} manifold is only about 95 cm^{−1} lower in energy than that of *ν*_{3}. As such, the states in the *ν*_{1} manifold that come into near-degeneracy with those in the *ν*_{3} manifold are not markedly different with respect to their level of intermolecular excitation. Thus, coupling matrix elements can be appreciable. A similar situation also obtains for the HDO complex in regard to the *ν*_{OD} and the slightly higher-lying 2*ν*_{2} manifold.

### D. Near-infrared spectra

Figure 2 shows computed IR spectra corresponding to the OH stretching region of the H_{2}O–benzene complex. The upper spectrum is that arising from the *ortho* ground state of the complex as the initial state [i.e., the lowest-energy *E*_{1}, *m*(1) state], and the lower one is that arising from the *para* ground state [i.e., the overall, *m*(0), ground state]. In each spectrum, the red bands correspond to parallel transitions and the blue ones correspond to perpendicular ones. In regard to the former, there is one dominant band in each spectrum (“a” and “A” in upper and lower, respectively). These two, which are located within about 0.1 cm^{−1} of one another, correspond to the *ν*_{1} fundamental for the *ortho* and *para* species, and each is characterized by Δ*m* = 0. In regard to the perpendicular bands, there are two dominant ones (“b” and “c”) in the *ortho* spectrum and one (“B”) in the *para* spectrum. Each one of these bands arises from the water-moiety *ν*_{3} fundamental, but involves, in addition, a change of |Δ*m*| = 1 in the internal rotation state of the complex. Band “b” involves internal-rotation state changes *m*(1) → *m*(0). Band “c” involves the state changes *m*(1) → *m*(2). Finally, band “B” involves *m*(0) → *m*(1). Comparison of the aforementioned dominant bands in the spectra of Fig. 2 with the experimental results and analysis of Pribble *et al.*^{6} (Figs. 8 and 10 of that reference, respectively) reveals excellent agreement with respect to relative band positions and band assignments. From this agreement, in particular that associated with bands “b,” “B,” and “c,” one concludes that the low-energy internal rotation states in the intramolecular ground-state and *ν*_{3} manifolds of the complex appear to be modeled very well with the PES that we have used. To be specific, our bands “b” and “c” are separated by the computed energy difference between the *m*(0) and *m*(2) states in the *ν*_{3} manifold: 63.9 cm^{−1}. The corresponding experimental value [see Fig. 10(b) of Ref. 6] is about 61.5 cm^{−1}. The separation between band “b” and “B” in Fig. 2 corresponds to the computed sum of the *m*(1) − *m*(0) energy differences in the ground state and *ν*_{3} manifolds, which equals 33.9 cm^{−1}. The corresponding experimental value is about 35.5 cm^{−1}.

In addition to the large-intensity bands in Fig. 2, there are a number of smaller-intensity ones, the origins of which deserve comment. The largest such band is “C” in the bottom spectrum. It nominally corresponds to the transition from the *para* ground state to the *m*(2) internal rotation state built on *ν*_{1} + *ν*_{XY}. However, it gains its intensity by virtue of the mixing of this latter state with the *m*(1) internal rotation level of *ν*_{3}. As such, it serves as a good example of the possibility that higher-lying states in the *ν*_{1} manifold might acquire transition strength in the IR due to coupling with lower-lying “bright” states in the *ν*_{3} manifold. The “D” and “E” bands in the lower spectrum are examples of the fact that intermolecular-state changes other than those associated with internal rotation can accompany the OH stretching fundamentals and lead to appreciably intense IR bands. In the case of “D,” the final state is *ν*_{1} + *ν*_{XY}. In the case of “E” it is *ν*_{3} + *ν*_{rock}. Thus, both *ν*_{XY} and *ν*_{rock} can, in principle, contribute structure to the OH stretching region beyond that due to the dominant bands that involve |Δ*m*| = 1 internal-rotation state changes. Indeed, Pribble *et al.*^{6} have previously noted this possibility in conjunction with the rocking vibration. One expects such a structure to be especially sensitive to the details of the PES. As such, given the rather rudimentary nature of the PES we employ, it is not surprising that we cannot positively assign our computed smaller-intensity bands to corresponding features in the measured spectra.^{6} In fact, we present evidence below that more accurate treatments of both *ν*_{rock} and *ν*_{XY} than those that are available from the present PES are necessary to model observed spectra quantitatively.

Figure 3 shows the computed spectrum in the OH-stretching fundamental region for the HDO complex. The initial state corresponding to this spectrum is the *m*(0), *A*_{1} ground state of the complex, the only one expected to have significant population under the very-low-temperature conditions of a seeded supersonic expansion. As in Fig. 2, the red bands are parallel transitions and the blue bands are perpendicular ones. There are several points to note about the spectrum. First, the structure begins with the appreciably intense “A” band, which is the bare *ν*_{OH} fundamental. Second, however, the next and most intense band “B” involves *ν*_{OH} plus a |Δ*m*| = 1 internal-rotation level change. Finally, there are higher-frequency bands with significant intensity that correspond to final states in which the intermolecular excitations *ν*_{rock} and *ν*_{XY} are involved. For example, band “E,” which is the second-most intense band in the spectrum, involves the final state *ν*_{OH} + *ν*_{rock} + *m*(1). Similarly, band “C” involves the final state *ν*_{OH} + *ν*_{rock}. The involvement of *ν*_{XY} is reflected in band “D,” whose final state is *ν*_{OH} + *ν*_{XY}. Comparison of Fig. 3 with the measured spectrum of Fig. 5(a) (Ref. 6) reveals that, while there is not a quantitative match between the two, there are notable points of agreement. First, the two lowest-frequency bands in each spectrum are spaced by close to the same frequency (11.7 cm^{−1} in Fig. 3 and ∼10 cm^{−1} in the experiment). Second, both computed and measured spectra consist of five significantly intense bands within a range of about 50 cm^{−1}. Given these points of agreement, together with what we know about the nature of the bands in the computed spectrum, we offer the following provisional interpretation of the experimental spectrum.

First, since (a) the spacing of the first two bands in the computed spectrum corresponds to the *m*(1) − *m*(0) energy difference in the *ν*_{OH} manifold, (b) the results on the H_{2}O complex indicate that our calculations do well at modeling the low-energy internal rotation states of that complex, and since (c) the measured band spacing of the two lowest-frequency HDO bands is close to the computed spacing, we attribute these two measured bands to the same transitions as those that are associated with the corresponding computed ones. Second, as to the higher-frequency bands in the experimental spectrum, we assume that they arise from the same transitions that produce bands “C,” “D,” and “E” in the computed spectrum. Of course, this can only be correct if our calculations do not accurately model *ν*_{rock} and *ν*_{XY}. However, if *ν*_{rock} and *ν*_{XY} in the *ν*_{OH} manifold of the HDO complex are actually ∼21 cm^{−1} and ∼49 cm^{−1}, respectively (rather than the 33.3 cm^{−1} and 35.0 cm^{−1} values that arise from the PES we use), then the relative positions of the third, fourth, and fifth measured bands are consistent with transitions to the final states *ν*_{OH} + *ν*_{rock}, *ν*_{OH} + *ν*_{rock} + *m*(1), and *ν*_{OH} + *ν*_{XY}, respectively. While this conjuring up of new frequencies for the rocking and in-plane translational modes may seem to be *ad hoc*, there is evidence from the intermolecular Raman spectra that further supports this scenario, as we describe below.

### E. Intermolecular Raman spectra

Figure 4 shows computed intermolecular Raman spectra for the H_{2}O–benzene complex. The upper and lower spectra correspond, respectively, to the *ortho* and *para* ground states of the species as the initial state. The one dominant feature in the *para* spectrum, band “A” at 36.1 cm^{−1}, is due to the *ν*_{XY} fundamental. A second minor band in the same spectrum, “B” at 71.3 cm^{−1}, corresponds to a final state that is best described as *ν*_{rock} + *m*(1). In the *ortho* spectrum, band “a” at 37.7 cm^{−1} corresponds to the transition from *m*(1) to *ν*_{XY} + *m*(1) (i.e., the *ortho ν*_{XY} fundamental). The only other bands with appreciable intensity in the *ortho* spectrum appear at 26.6 cm^{−1} and 52.3 cm^{−1}. These two, bands “b” and “c,” arise from transitions in which the final states are *ν*_{rock} and *ν*_{XY} + *ν*_{rock}, respectively.

The measured intermolecular Raman spectrum for the H_{2}O complex, which has contributions from both *ortho* and *para* species, exhibits four significantly intense bands at about 8 cm^{−1}, 40 cm^{−1}, 44 cm^{−1}, and 52 cm^{−1}. Of these, the most intense one at 44 cm^{−1} is almost certainly due to the *ν*_{XY} fundamentals for the *para* and *ortho* species.^{7} The assignments of the other three bands have yet to be made. However, if one works under the assumption that the measured Raman bands arise from the same types of transitions that have significant Raman intensity in the computed spectra and if one assumes that our calculations model the low-energy internal rotation states in the ground-state manifold reasonably accurately (see above), then an internally consistent set of Raman assignments can be had if one takes *ν*_{rock} ≃ 24.5 cm^{−1} and *ν*_{XY} ≃ 44 cm^{−1} in the ground-state manifold of the H_{2}O complex. Given these frequencies, the highest-intensity bands “A” and “a” would appear unresolved near 44 cm^{−1}. Bands “b” and “c” would appear at ∼8 cm^{−1} and 52 cm^{−1}, respectively. Finally, band “B” would occur in the vicinity of 41 cm^{−1}. In short, the observed Raman spectrum is very-well modeled under this set of assumptions. Moreover, these values for *ν*_{rock} and *ν*_{XY} in the ground-state manifold of the H_{2}O complex are close to the analogous frequencies (i.e., 21 cm^{−1} and 49 cm^{−1}, respectively) in the *ν*_{OH} manifold of the HDO complex that successfully modeled the IR spectral features of the latter species (see Sec. III D).

## IV. CONCLUSIONS

We have presented an efficient computational methodology for rigorous quantum bound-state calculations of the coupled intra- and intermolecular vibrational states of benzene–H_{2}O and benzene–HDO dimers in which H_{2}O and HDO are treated as flexible and benzene as rigid. Apart from the assumption of rigid benzene, the quantum 9D treatment is fully coupled. The calculations yield H_{2}O/HDO intramolecular vibrational fundamentals and the bend (*ν*_{2}) overtone as well as the low-lying intermolecular vibrational states in each intramolecular vibrational manifold. The general procedure introduced in Ref. 20 is followed in which the full 9D vibrational Hamiltonian of the dimer is partitioned into a 6D intermolecular Hamiltonian, a 3D intramolecular Hamiltonian, and a 9D remainder term. The low-energy eigenstates of the two reduced-dimension Hamiltonians are used for constructing the 9D product basis in which the full Hamiltonian is diagonalized. Symmetry of the system is exploited to bring the Hamiltonian matrix to a block-diagonal form, where each block is associated with one of the irreducible representations of the molecular symmetry group employed (*G*_{12}). In addition, we have presented the methods for computing accurately the IR transition intensities between the 9D eigenstates of the complexes, as well their low-frequency Raman spectra.

The bound-state method developed in the course of the present investigation owes much of its efficiency to the finding made in our recent studies^{19–22} that converged intramolecular vibrational excitations in weakly bound molecular complexes can be obtained by means of fully coupled quantum calculations in which the contracted intermolecular vibrational eigenstates included in the full-dimensional basis span the energies far below those of the intramolecular vibrational excitations of interest. In the present work, for both the H_{2}O and HDO complexes, the 6D intermolecular contracted bases for each symmetry block include only 40 eigenstates whose energies extend up to about 225 cm^{−1} that is much less than the stretch and bend fundamental frequencies of H_{2}O and HDO, which range between 1400 cm^{−1} and 3800 cm^{−1}. This results in matrices representing the symmetry blocks of the 9D Hamiltonian that are small for the 9D quantum problem, 1360 and 1680 for the H_{2}O and HDO complexes, respectively, which allows for straightforward direct diagonalization.

What has emerged from the 9D quantum bound state calculations and the simulations of the IR and Raman spectra is a comprehensive description of the H_{2}O/HDO intramolecular excitations coupled to the low-energy LAM eigenstates of the two dimers within the intramolecular manifolds and their manifestations in the spectra. Numerous interesting facets of the highly quantum behavior of the H_{2}O and HDO complexes are revealed. Thus, plots of the probability density functions show that in the ground intermolecular state of the H_{2}O complex, both H atoms of water interact equivalently with the benzene, as deduced experimentally.^{1,8} For the HDO complex, the same plots show that only the D nucleus is localized close to the benzene plane while the H atom is further away, in accord with the experimental observations.^{5,8} However, upon excitation of the low-frequency *ν*_{rock} mode, this bias is reversed, and now, it is the H atom that has a greater probability of being close to the benzene than the H atom.

The 9D quantum calculations of the four lowest-energy intramolecular vibrational states of H_{2}O and HDO in the complexes show that the stretching modes of H_{2}O (*ν*_{1} and *ν*_{3}) and HDO (*ν*_{OD} and *ν*_{OH}) are redshifted, while the H_{2}O/HDO bend fundamentals (*ν*_{2}) and overtones (2*ν*_{2}) are blueshifted, relative to the ones calculated for the respective isolated monomers. For the *ν*_{1} and *ν*_{3} modes of H_{2}O and the *ν*_{OH} mode of HDO, the prediction of their redshift agrees with the experimental results, but the redshift magnitudes are roughly a factor of two smaller the corresponding experimental values. No experimental data are available regarding the frequency shifts of the *ν*_{2} mode of gas-phase benzene–H_{2}O and benzene–HDO.

The quantum 9D calculations provide rich data about the effects that the excitation of various intramolecular modes of H_{2}O and HDO has on the LAM low-frequency intermolecular states of the two dimers. We find that the intermolecular modes that primarily involve the motions of the c.m. of H_{2}O/HDO relative to the benzene plane, *ν*_{XY} and *ν*_{Z}, upon the excitation of any of the intramolecular modes change only slightly, by 1–2 cm^{−1}, from their respective values in the ground intramolecular vibrational state. In contrast, the energies of the intermolecular modes that execute LAM motions of the H/D atoms, the internal-rotation mode (*m*) and *ν*_{rock}, are sensitive to the excitation of the stretching modes of H_{2}O and HDO, and this manifests in intricate ways, particularly in the case of the HDO complex.

The most complete comparison between theory and experiment is achieved by contrasting the simulated and measured IR^{6} and Raman^{7} spectra of the dimers. The dominant bands in the computed IR spectra of the H_{2}O complex in the OH-stretch region are in excellent agreement with the experimental results and analysis of Pribble *et al.*,^{6} regarding the relative band positions and band assignments, implying that the low-energy internal-rotation states in the intramolecular ground-state and *ν*_{3} manifolds of the complex are described well by the PES employed. However, this does not appear to be the case for the *ν*_{XY} and *ν*_{rock} modes. As a result, a number of smaller-intensity bands cannot be assigned with confidence to particular features in the spectra, reflecting the limitations in the accuracy of the PES. In the case of the HDO complex, the IR spectrum computed in the region of the OH stretch, while not providing a quantitative match with the measured spectrum,^{6} agrees with it in a number of important features, allowing us to interpret it tentatively after making some reasonable assumptions about the energies of *ν*_{XY} and *ν*_{rock}. The computed intermolecular Raman spectrum of benzene–H_{2}O can be brought in correspondence with the measured Raman spectrum,^{7} and the latter interpreted, under similar assumptions about *ν*_{XY} and *ν*_{rock}.

The rigorous three-pronged theoretical treatment of benzene–H_{2}O and benzene–HDO, combining fully coupled 9D quantum calculations of their intra- and intermolecular vibrations, as well as the IR and Raman spectra, is essential for full assessment of the current PES and any other that may be developed in the future. Such development is highly desirable in view of the deficiencies of the PES employed identified in this investigation. Moreover, the methodology formulated in this study is appropriate for a broad class of weakly bound molecular dimers beyond benzene–water in which one of the two moieties is significantly larger than the other, resulting in a strongly anisotropic intermolecular potential, such as a flexible triatomic molecule bound to a polyacene, planar (e.g., naphthalene and anthracene) or curved (e.g., corannulene^{43}), treated as rigid. After some simplifications of the Hamiltonian, the same approach is directly applicable also to the 9D quantum vibrational dynamics of a triatomic molecule physisorbed on a solid surface. Investigations along these lines are under way.

## SUPPLEMENTARY MATERIAL

See Sec. S1 of the supplementary material for an extended compilation of results from the diagonalization of *Ĥ*_{inter} and *Ĥ*_{intra}.

## ACKNOWLEDGMENTS

Z.B. is grateful to the National Science Foundation for its partial support of this research through the Grant No. CHE-1566085. P.F. is grateful to Professor Daniel Neuhauser for his support.

## REFERENCES

The analogous situation obtains for the H_{2}O@C_{60} system in regard to a five-fold axis of symmetry. See Supplementary Material, Sec. IV of Ref. 22.