We present a rigorous and comprehensive theoretical treatment of the vibrational dynamics of benzene–H2O 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 H2O 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 H2O 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 H2O and HDO complexes, respectively, allowing for direct diagonalization. These calculations characterize in detail the H2O/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.

Size-selected small benzene–(H2O)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 Raman7 spectroscopy. The microwave and Raman spectroscopies probed the two dimers in the ground intramolecular vibrational state, providing information about their vibrationally averaged geometries1,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–H2O 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 H2O 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–H2O 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 microwave1,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 C6 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 H2O in the π hydrogen bonding with benzene is limited to the 1:1 complex. In larger benzene–(H2O)n clusters, with n ≥ 2, benzene resides on the exterior, or the surface, of the (H2O)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–H2O, while no such combination bands are observed in the symmetric-stretch region.3,6 For both benzene–H2O 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 H2O/HDO intramolecular vibrations. The features arising from the LAMs are present also in the Raman spectra of benzene–H2O, 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–H2O 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–H2O (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 H2O on the C6 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–H2O intermolecular interaction, for rigid monomers, was also characterized by means of density functional theory (DFT) calculations.11 For benzene–(H2O)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 computed9,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–(H2O)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–(H2O)n clusters with n ≥ 2 where, as mentioned above, the (H2O)n cluster donates a single free OH group to the π hydrogen bond, and the cluster rigidity increases with n. However, benzene–H2O/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–H2O 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 H2O 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–H2O 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–H2O/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 H2O/HDO as rigid. A rigorous approach to the H2O intramolecular vibrations coupled to the intermolecular vibrational states of the dimer requires fully coupled 9D quantum bound-state calculations for flexible H2O and rigid benzene, extending to the energies of the H2O 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 H2O 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–H2O/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 H2O 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 H2O/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 H2, HD, and D2 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 H2 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 H2 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 H2 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., H2O2, CH4, 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 approach20 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 (H2)2 in the large hydrate cage21 and in the fully coupled 9D quantum treatment of the intra- and intermolecular vibrational levels of H2O in (rigid) C60,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 H2O 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 above20 and also builds on the earlier quantum 6D (rigid-monomer) treatment of the intermolecular vibrations of benzene–H2O.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., H2O), 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 G12 molecular symmetry group. Their dimensions, 1360 and 1680, for each G12 irrep of benzene–H2O 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–H2O/HDO in the OH-stretch fundamental region. Also calculated is the Raman spectrum of benzene–H2O 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.

The paper is organized as follows: Computational methodology is described in Sec. II. In Sec. III, we present and discuss the results. Section IV presents the conclusions.

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 work17 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^×Ŷ. We can then write the classical nuclear kinetic energy of this system (apart from overall translational energy) as

T=|pd|22μd+i(jiB)22Ii+TWF,
(1)

where μd is the reduced mass of the complex, mwatermbenzene/(mwater + mbenzene), pd 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 ith axis, Ii is the moment of inertia about the benzene principal axis parallel to i, and TWF 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

Ĥ=d22μd+i=X^,Ŷ,Z^L^i2(d)+(ĵiW(ω))2+2L^i(d)ĵiW(ω)2Ii+T^WF(ω,q)+Vtot(d,ω,q),
(2)

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, d2 is the Laplacian associated with d, L^i is the operator corresponding to the component, measured along the ith 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, ĵiW is the operator corresponding to the component of water’s angular momentum measured along the ith CF axis,

T^WF(ω,q)=T^vrWF(ω,q)+T^corWF(ω,q)+T^vibWF(q)
(3)

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

Vtot(d,ω,q)=VBW(d,ω,q)+VintraWF(q),
(4)

where the first term on the rhs of Eq. (4) is the atom–atom pairwise-additive benzene–water interaction potential from the work of Karlström et al.10 while the second term is the 3D PES for the isolated, gas-phase water molecule from the study of Mizus et al.28 

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κinter) and those of Ĥintra (which we denote as |γ⟩, with eigenenergies Eγintra), we build up a contracted product basis |κ, γ⟩ ≡|κ⟩|γ⟩, which we then use to diagonalize the full Hamiltonian. We define Ĥinter as

Ĥinterd22μd+i=X,Y,ZL^i2(d)+(ĵiW(ω))2+2L^i(d)ĵiW(ω)2Ii+T^vrWF(ω;q0)+Vinter(d,ω),
(5)

where

Vinter(d,ω)VBW(d,ω;q0)
(6)

and q0 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

ĤintraT^vibWF(q)+Vadap(q),
(7)

where the “dimer-adapted” 3D intramolecular potential, Vadap, is the hydrogen-exchange-symmetrized version of

Vadap(q)VintraWF(q)+VBW(d0,ω0,q).
(8)

Here, d0 and ω0 are fixed values of d and ω. We choose these values to correspond to those at, or near, the minimum of VB-W when the q are fixed to q0. The resulting Vadap thus includes some of the effects of the benzene–water interaction on the water vibrational PES.

This leaves

ΔĤΔT^vrWF(ω,q)+T^corWF(ω,q)+ΔV(d,ω,q),
(9)

where

ΔT^vrWF(ω,q)T^vrWF(ω,q)T^vrWF(ω;q0)
(10)

and

ΔV(d,ω,q)VBW(d,ω,q)VBW(d,ω;q0)VBW(d0,ω0,q).
(11)

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

κ,γ|Ĥ|κ,γ=(Eκinter+Eγintra)δκ,κδγ,γ+κ,γ|(ΔT^vrWF+T^corWF+ΔV)|κ,γ.
(12)

We now specify our choices for the particular coordinates on which Ĥ depends. For the components of d, we use the cylindrical coordinates (dZ, ρ, Φ), where dZ is the cartesian component of d along the axis, ρdX2+dY2, and (cos Φ, sin Φ) ≡ (dX/ρ, dY/ρ), with dX and dY being the Cartesian components of d along the X^ and Ŷ axes, respectively. For the intra-water q, we use the Radau coordinates (R1, R2, Θ).29 Furthermore, we use Radau bisector-z embedding to define the WF frame.27 Specifically, we take the WF axis to be parallel to (R^1+R^2), where R^1 and R^2 are the unit vectors associated with the Radau vectors R1 and R2. The WF ŷ axis is taken to be parallel to R^1×R^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.

Apart from T^WF, the kinetic-energy portion of Ĥ [see Eq. (2)] can now be written as17 

d22μd+iL^i2(d)2Ii+i(ĵiW(ω))22Ii+i2L^i(d)ĵiW(ω)2Ii=T^1+T^2+T^3,
(13)

where

T^1=12μd+AdZ22ρ2+1ρρ+1ρ22Φ2+AL^Z212μd+Aρ22dZ2+2AdZdZ+Aρρ1+2dZdZ,
(14)
T^2=A(ĵW)2+(AA)(ĵZW)2,
(15)

and

T^3=ĵXWsinΦĵYWcosΦ2AiρdZ+idZρ2AdZρ[ĵXWcosΦ+ĵYWsinΦ]2AĵZWL^Z,
(16)

where A ≡ 1(2IX) = 1/(2IY) and A ≡ 1/(2IZ) = A/2 are the rotational constants of rigid benzene.

For the components of T^WF [Eq. (3)], one has27 

T^vrWF(ω,q)=B1(R1)+B2(R2)2(ĵxWF)21+cosΘ+(ĵyWF)22+(ĵzWF)21cosΘ+B1(R1)B2(R2)2ĵxWFĵzWF+ĵzWFĵxWFsinΘ,
(17)
T^corWF(ω,q)B1(R1)B2(R2)p^ΘHĵyWF,
(18)

and

T^vibWF(q)=i=1212μi2Ri2B1(R1)+B2(R2)1sinΘΘsinΘΘ,
(19)

where Bi(Ri)1/(2μiRi2), μi is the Radau mass factor associated with Ri,29 the ĵkWF (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

pΘHiΘ+12cotΘ.
(20)

Note that the (ĵXW,ĵYW,ĵZW) appearing in Eqs. (15) and (16) correspond to components of the water angular momentum measured along the CF axes, while the (ĵxWF,ĵyWF,ĵzWF) 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.

1. Primitive intermolecular basis set

We choose a “cylindrical-Wigner” basis17 to represent the matrix of Ĥinter. These basis states are of the form

|α,v,l,j,mj,k=|dZ,α|v,l|j,mj,k,
(21)

where the |dZ,α⟩, α = 1, … Nα, comprise a 1D Gauss–Hermite DVR covering the dZ 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, mj, k⟩ (j = 0, 1, …, jmax, mj, 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–water17 and was described in detail there. In the current work, we take Nα = 20, vmax = 10, and jmax = 11. The |dZ,α⟩ were generated from 1D HO wavefunctions centered at dZ = 7.6 bohr (z0 in Ref. 17) with the mass × angular-frequency parameter (γz2 in Ref. 17) set to (1.428 778)2 a.u. The resulting 20 dZ quadrature points ranged from 3.829 to 11.371 bohr. The 2D HO eigenfunctions had their mass × angular-frequency parameter (γρ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 G12 molecular symmetry group (isomorphic to the D6 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 H2O 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 + mj) mod 6 in a given symmetry block. This divides the states into six sets corresponding to the A irreps [(l + mj) mod 6 = 0], the B irreps [(l + mj) mod 6 = ±3], the two components of the E1 irrep [(l + mj) mod 6 = +1 and −1], and the two components of the E2 irrep [(l + mj) mod 6 = +2 and −2]. With the values of Nα, vmax, and jmax that we use, these symmetry blocks are composed of 505 760, 506 240, 506 240, and 505 760 states for the A, B, E1 (either component), and E2 (either component) irreps, respectively.

The second step, which factors the A block into A1 and A2 blocks and the B block into B1 and B2 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 G12, the |γ⟩ are unaffected by those operations, and the |κ⟩ transform like G12 irreps.

3. Iterative diagonalization procedure

We use the Chebyshev version31 of filter diagonalization32 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 dZ. For those, we first evaluate analytically the operator matrix in the 1D HO-eigenfunction representation that is isomorphic to the dZ 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(ω;q0).

To operate with the potential term in Eq. (5) [i.e., Vinter(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, mj, 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 Vinter(d, ω) at the associated grid point. Finally, we back-transform the result to the |α, v , l, j, mj, 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 (mj, 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π/Nr, where r = 1, …, Nr], 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 + mj)-symmetry factorization of the basis.33 The size of the 6D grid that we use is NαNβNγNqNrNs/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

(35)(26)*|α,v,l,j,mj,k=[|α,v,l,j,mj,k]*

and

(35)(26)*Vinter[(35)(26)*]1=Vinter
(22)

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.

To completely specify Ĥintra, the values of d0 and ω0 that define Vadap [see Eq. (8)] need to be fixed. For each isotopologue, we treat two cases. The first is the “bare-monomer” case corresponding to |d0| = and Vadap=VintraWF. The second is the “dimer-adapted” case, for which we take d0 ≡ (dZ,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 VB-W(d0, ω0, q0) = −1076 cm−1 relative to the energy of the infinitely separated, rigid monomers. This energy is very close to the minimum of Vinter.

Our procedure for computing the eigenfunctions of Ĥintra is very similar to that described in Ref. 22 in connection with a study of H2O@C60. Briefly, for a given species, we use a primitive basis consisting of the product of two tridiagonal-Morse DVRs34 covering the R1 and R2 coordinates (|η1⟩ and |η2⟩, each consisting of nine functions) and a potential-optimized DVR35 (|η3⟩, consisting of 15 functions) covering the Θ coordinate. The basis used here for H2O is identical to the intramolecular basis of Ref. 22. The basis used for HDO was constructed in the same way as that for H2O (Sec. V A of Ref. 22) but with Radau mass factors appropriate to the HDO species.

Matrix elements of Vadap 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.

As mentioned in Sec. II D 2, the matrix of Ĥ is block diagonal in the |κ, γ⟩ basis, with the different blocks corresponding to the different G12 irreps to which the various |κ⟩ belong. In the case of benzene–H2O, for all irreps, we set Ninter = 40 and Nintra = 34. For benzene–HDO, Ninter = 40 and Nintra = 42 were employed for all irreps. Characteristics of all of these Ninter = 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, Ninter × Nintra, equal to 1360 for each G12 irrep of the H2O 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 Ninter was tested for the A1, E1, and E2 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 Nintra, its value of 34 for the H2O complex was taken from our study of H2O@C60,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, Nintra = 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 H2O 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 ΔT^vrW, T^corW, and ΔV [Eq. (12)]. Since Δ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

κ,γ|ΔV|κ,γ=ηγ|ηη|γκ,η|ΔV|κ,η,
(23)

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, mj, 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 κ′,

κ,η|ΔV|κ,η=Qκ|QΔV(Q,η)Q|κ.
(24)

This step only requires computing and storing a 6D piece of ΔV. One then proceeds to cycle through all the η, computing for each a 6D piece of ΔV and evaluating Eq. (24) for all κ and κ′ until all the matrix elements required for the evaluation of Eq. (23) are computed.

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, μKECF (K = −1, 0, +1), referred to an Eckart complex-fixed (ECF) frame,

μKECFifψf|μKECF|ψi.
(25)

We then use these to estimate transition intensities as

IifK|μKECFif|2.
(26)

To compute the μKECFif, 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:

μKECFif=Mψf|DM,K(1)(Ω)μMCF|ψi,
(27)

where the μMCF are spherical-tensor dipole-component operators referred to the CF frame, the DM,K(1)(Ω) 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 μ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,

ψf|DM,K(1)(Ω)μMCF|ψi=κ,γκ,γψf|κ,γκ,γ|ψi×κ|DM,K(1)(Ω)γ|μMCF|γ|κ,
(28)

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

γ|μMCF|γ=K[DM,K(1)(ω)]*γ|μKW|γ,
(29)

where the μKW are dipole-operator spherical-tensor components referred to the WF frame. Equation (28) can then be written as

ψf|DM,K(1)(Ω)μMCF|ψi=Kγ,γγ|μKW|γκ,κψf|κ,γ×κ,γ|ψiκ|DM,K(1)(Ω)[DM,K(1)(ω)]*|κ.
(30)

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)(Ω) and DM,K(1)(ω) 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(1)(ω) 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,

κ|DM,K(1)(Ω)[DM,K(1)(ω)]*|κ  =Qκ|QQ|κDM,K(1)(Ω(Q))DM,K(1)ωQ*.
(31)

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 H2O 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 H2O, the relative magnitudes of these dipole vectors are in the ratio |μν1|/|μν3|1/20.40 In the benzene–H2O 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.

We also model herein the Raman spectrum of the H2O 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

SifK|[αK(2)]ECFif|2,
(32)

where

[αK(2)]ECFifψf|[αK(2)]ECF|ψi
(33)

and the [α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

[αM(2)]CF=[αM(2)]B+q[DM,q(2)(ω)]*[αq(2)]W,
(34)

where [α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 [α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

[αK(2)]ECFif=M[αM(2)]Bψi|DM,K(2)(Ω)|ψf+q[αq(2)]Wψi|DM,K(2)(Ω)[DM,q(2)(ω)]*|ψf.
(35)

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.

Selected low-energy intermolecular eigenstates of benzene–H2O 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 H2O 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 H2O 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 C6 axis, with the H atoms of the water moiety directed toward the benzene (π hydrogen bonding).

TABLE I.

Molecular parameters (a.u., except where otherwise noted) used in Ĥinter in Eq. (5) for benzene–H2O and benzene–HDO. Radau distances and mass factors conform to the convention of Ref. 29. For both isotopologues, the fixed Radau coordinates correspond to a rigid water geometry with OH/OD bond distances of 0.96 Å and a bond angle of 104.52°.

A (cm−1)μdμ1μ2R1,0R2,0cos(Θ0)
H20.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μ2R1,0R2,0cos(Θ0)
H20.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 
TABLE II.

Selected low-energy intermolecular states of benzene–H2O from 6D (rigid H2O) and 9D (flexible H2O) calculations. The latter are for H2O in the ground intramolecular vibrational state. “Basis-state norm” refers to the 9D results.

ParaOrthoBasis-state
NΔE: 6D/9DΔE: 6D/9DIrrepAssignmentnorm
0.00/0.00  A1+ … 0.999 
 16.32/16.53 E1 m(1) 0.999 
35.36/36.07  E1+ νXY 0.999 
 41.14/43.09 A1 νrock 0.998 
 53.22/54.21 E2 νXY + m(1) 0.999 
61.24/62.13  A1+ 2νXY 0.998 
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 
ParaOrthoBasis-state
NΔE: 6D/9DΔE: 6D/9DIrrepAssignmentnorm
0.00/0.00  A1+ … 0.999 
 16.32/16.53 E1 m(1) 0.999 
35.36/36.07  E1+ νXY 0.999 
 41.14/43.09 A1 νrock 0.998 
 53.22/54.21 E2 νXY + m(1) 0.999 
61.24/62.13  A1+ 2νXY 0.998 
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 
TABLE III.

Selected low-energy intermolecular states of benzene–HDO from 6D (rigid HDO) and 9D (flexible HDO) calculations. The latter are for HDO in the ground intramolecular vibrational state. “Basis-state norm” refers to the 9D results.

NΔE: 6D/9DIrrepAssignmentBasis-state norm
0.00/0.00 A1 … 0.991 
12.59/13.27 E1 m(1) 0.991 
35.01/36.42 E1 νXY 0.978 
36.38/41.28 A1 νrock 0.971 
47.00/48.25 E2 νXY + m(1) 0.964 
50.41/52.93 E2 m(2) 0.953 
56.46/58.16 A1 2νXY 0.971 
11 68.98/71.78 E2 2νXY 0.971 
16 77.69/79.68 A1 νZ 0.964 
NΔE: 6D/9DIrrepAssignmentBasis-state norm
0.00/0.00 A1 … 0.991 
12.59/13.27 E1 m(1) 0.991 
35.01/36.42 E1 νXY 0.978 
36.38/41.28 A1 νrock 0.971 
47.00/48.25 E2 νXY + m(1) 0.964 
50.41/52.93 E2 m(2) 0.953 
56.46/58.16 A1 2νXY 0.971 
11 68.98/71.78 E2 2νXY 0.971 
16 77.69/79.68 A1 ν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 C6 axis) mode, with internal-rotation constants B ≃ 16/13 cm−1 (H2O/HDO) and with the internal-rotation energy levels given by Em = Bm2 (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 H2O 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 H2O species, the HDO ground state has ⟨dZ⟩ ≃ 6.2 bohr and ⟨dX⟩ = ⟨dY⟩ = 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 H2O 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,

P(Z1,Z2)δ(Z1Z1)δ(Z2Z2),
(36)

where the angle brackets on the rhs of Eq. (36) denote an expectation value and Z1 and Z2 represent the Z (CF frame) coordinates of the two hydrogen nuclei of the water moiety. Figure 1(a) is a plot of P(Z1, Z2) for the Ĥinter ground state of the H2O 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(Z1, Z2) plots for the intermolecular states that we assign to νrock for the H2O 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 H2O 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.

FIG. 1.

Contour plots of the two-point probability density P(Z1, Z2) for four benzene–water intermolecular eigenfunctions. (a) Ground state of benzene–H2O. (b) Ground state of benzene–HDO. (For the HDO complex, Z1 corresponds to the H nucleus and Z2 corresponds to the D nucleus.) (c) νrock state of benzene–H2O. (d) νrock state of benzene–HDO. In (c) and (d), the contour corresponding to the largest local probability-density magnitude is labeled with that magnitude.

FIG. 1.

Contour plots of the two-point probability density P(Z1, Z2) for four benzene–water intermolecular eigenfunctions. (a) Ground state of benzene–H2O. (b) Ground state of benzene–HDO. (For the HDO complex, Z1 corresponds to the H nucleus and Z2 corresponds to the D nucleus.) (c) νrock state of benzene–H2O. (d) νrock state of benzene–HDO. In (c) and (d), the contour corresponding to the largest local probability-density magnitude is labeled with that magnitude.

Close modal

Table IV displays results, computational and experimental (where available), that pertain to the intramolecular vibrational excitations of H2O 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 H2O and HDO moieties from 3D quantum calculations for VintraWF, the 3D PES of the isolated H2O,28 and Vadap, the dimer-adapted 3D H2O PES, respectively. In the case of isolated H2O, 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 values18 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 H2O over a broad range of energies.28 This optimization probably made the PES slightly less accurate for HDO.

TABLE IV.

Energies of the four lowest-energy H2O and HDO excited intramolecular vibrational states in benzene–H2O/HDO computed for three PESs. VintraWF is the 3D PES of the isolated water molecule.28Vadap (3D) is the dimer-adapted water potential defined by Eq. (8) and the intermolecular geometrical parameters specified in Sec. II E. Vtot (9D) is defined in Eq. (4). Next to them are the calculated frequency shifts, Δν3D obtained for Vadap, and Δν9D obtained for Vtot. They are relative to the energies computed for VintraWF in the second column. Also shown, under Expt., are the available experimental intramolecular vibrational fundamentals for the H2O and HDO complexes (Ref. 6) and their frequency shifts Δνexpt relative to the gas-phase intramolecular fundamentals (Ref. 18). All quantities shown are in cm−1.

PES:VintraWF (3D)Vadap (3D)Δν3DVtot (9D)Δν9DExpt.Δνexpt
H2       
ν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)Vadap (3D)Δν3DVtot (9D)Δν9DExpt.Δνexpt
H2       
ν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 H2O 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 H2O (ν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 H2O and HDO. In contrast, the H2O/HDO bend fundamentals (ν2) and overtones (2ν2) are significantly blueshifted; again, these blueshifts are similar for H2O 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 H2O 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 H2O@C60.22 The exception to this is the ν3 mode of H2O 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 H2O/HDO PES and the 9D quantum calculations is that the former are performed for the fixed position and orientation of H2O/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 H2O/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 H2O/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 H2O 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 H2O and the νOH mode of HDO. For both H2O and HDO, the comparison of the frequency shifts measured for H2O 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 H2O 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 H2O and HDO and the intermolecular DOFs of the two dimers.

Tables V–VIII show the set of 9D-calculated low-energy intermolecular levels of benzene–H2O, 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 H2O 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.

TABLE V.

Selected low-energy states in the ν2 manifold of benzene–H2O from 9D calculations. ΔE values are relative to the ν2 fundamental frequency at 1615.16 cm−1.

Npara ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.998 
 16.87 E1 m(1) 0.998 
35.69  E1+ νXY 0.998 
 42.11 A1 νrock 0.997 
 54.18 E2 νXY + m(1) 0.999 
61.43  A1+ 2νXY 0.997 
66.89  E2+ m(2) 0.980 
10 70.52  E2+ 2νXY 0.979 
13 78.35  A1+ νZ 0.997 
Npara ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.998 
 16.87 E1 m(1) 0.998 
35.69  E1+ νXY 0.998 
 42.11 A1 νrock 0.997 
 54.18 E2 νXY + m(1) 0.999 
61.43  A1+ 2νXY 0.997 
66.89  E2+ m(2) 0.980 
10 70.52  E2+ 2νXY 0.979 
13 78.35  A1+ νZ 0.997 
TABLE VI.

Selected low-energy states in the 2ν2 manifold of benzene–H2O from 9D calculations. ΔE values are relative to the ν2 first-overtone frequency at 3196.62 cm−1.

Npara ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.996 
 17.29 E1 m(1) 0.997 
35.21  E1+ νXY 0.996 
 40.68 A1 νrock 0.992 
 54.18 E2 νXY + m(1) 0.996 
60.58  A1+ 2νXY 0.993 
68.05  E2+ m(2) 0.804 
11 69.85  E2+ 2νXY 0.803 
13 77.40  A1+ νZ 0.993 
Npara ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.996 
 17.29 E1 m(1) 0.997 
35.21  E1+ νXY 0.996 
 40.68 A1 νrock 0.992 
 54.18 E2 νXY + m(1) 0.996 
60.58  A1+ 2νXY 0.993 
68.05  E2+ m(2) 0.804 
11 69.85  E2+ 2νXY 0.803 
13 77.40  A1+ νZ 0.993 
TABLE VII.

Selected low-energy states in the ν1 manifold of benzene–H2O from 9D calculations. ΔE values are relative to the ν1 fundamental frequency at 3645.16 cm−1.

para ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.990 
 16.39 E1 m(1) 0.991 
36.02  E1+ νXY 0.988 
 40.21 A1 νrock 0.942 
 54.19 E2 νXY + m(1) 0.990 
61.89  A1+ 2νXY 0.986 
65.13  E2+ m(2) 0.990 
11 71.18  E2+ 2νXY 0.984 
13 79.59  A1+ νZ 0.988 
para ΔEortho ΔEIrrepAssignmentBasis-state norm
0.00  A1+ … 0.990 
 16.39 E1 m(1) 0.991 
36.02  E1+ νXY 0.988 
 40.21 A1 νrock 0.942 
 54.19 E2 νXY + m(1) 0.990 
61.89  A1+ 2νXY 0.986 
65.13  E2+ m(2) 0.990 
11 71.18  E2+ 2νXY 0.984 
13 79.59  A1+ νZ 0.988 
TABLE VIII.

Selected low-energy states in the ν3 manifold of benzene–H2O from 9D calculations. ΔE values are relative to the ν3 fundamental frequency at 3739.39 cm−1.

para ΔEortho ΔEIrrepAssignmentBasis-state norm
 0.00 A1 … 0.859 
17.35  E1+ m(1) 0.670 
36.68  A1+ νrock 0.571 
 36.75 E1 νXY 0.782 
50.73  E2+ νXY + m(1) 0.418 
 63.91 E2 m(2) 0.686 
 63.40 A1 2νXY 0.603 
10  72.96 E2 2νXY 0.761 
12  79.58 A1 νZ 0.578 
para ΔEortho ΔEIrrepAssignmentBasis-state norm
 0.00 A1 … 0.859 
17.35  E1+ m(1) 0.670 
36.68  A1+ νrock 0.571 
 36.75 E1 νXY 0.782 
50.73  E2+ νXY + m(1) 0.418 
 63.91 E2 m(2) 0.686 
 63.40 A1 2νXY 0.603 
10  72.96 E2 2νXY 0.761 
12  79.58 A1 νZ 0.578 
TABLE IX.

Selected low-energy states in the ν2 manifold of benzene–HDO from 9D calculations. ΔE values are relative to the ν2 fundamental frequency at 1422.84 cm−1.

NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.986 
13.68 E1 m(1) 0.986 
36.33 E1 νXY 0.969 
41.57 A1 νrock 0.959 
48.64 E2 νXY + m(1) 0.919 
54.00 E2 m(2) 0.900 
57.82 A1 2νXY 0.963 
11 71.64 E2 2νXY 0.960 
16 79.12 A1 νZ 0.953 
NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.986 
13.68 E1 m(1) 0.986 
36.33 E1 νXY 0.969 
41.57 A1 νrock 0.959 
48.64 E2 νXY + m(1) 0.919 
54.00 E2 m(2) 0.900 
57.82 A1 2νXY 0.963 
11 71.64 E2 2νXY 0.960 
16 79.12 A1 νZ 0.953 
TABLE X.

Selected low-energy states in the νOD manifold of benzene–HDO from 9D calculations. ΔE values are relative to the νOD fundamental frequency at 2707.09 cm−1.

NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.950 
14.13 E1 m(1) 0.944 
37.40 E1 νXY 0.894 
46.89 A1 νrock 0.727 
49.37 E2 νXY + m(1) 0.849 
56.24 E2 m(2) 0.793 
61.97 A1 2νXY 0.705 
11 73.87 E2 2νXY 0.835 
16 80.79 A1 νZ 0.861 
NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.950 
14.13 E1 m(1) 0.944 
37.40 E1 νXY 0.894 
46.89 A1 νrock 0.727 
49.37 E2 νXY + m(1) 0.849 
56.24 E2 m(2) 0.793 
61.97 A1 2νXY 0.705 
11 73.87 E2 2νXY 0.835 
16 80.79 A1 νZ 0.861 
TABLE XI.

Selected low-energy states in the 2ν2 manifold of benzene–HDO from 9D calculations. ΔE values are relative to the ν2 first-overtone frequency at 2822.32 cm−1.

NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.976 
14.09 E1 m(1) 0.896 
36.22 E1 νXY 0.950 
41.45 A1 νrock 0.925 
48.96 E2 νXY + m(1) 0.808 
55.22 E2 m(2) 0.685 
57.25 A1 2νXY 0.957 
12 71.46 E2 2νXY 0.932 
16 78.55 A1 νZ 0.927 
NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.976 
14.09 E1 m(1) 0.896 
36.22 E1 νXY 0.950 
41.45 A1 νrock 0.925 
48.96 E2 νXY + m(1) 0.808 
55.22 E2 m(2) 0.685 
57.25 A1 2νXY 0.957 
12 71.46 E2 2νXY 0.932 
16 78.55 A1 νZ 0.927 
TABLE XII.

Selected low-energy states in the νOH manifold of benzene–HDO from 9D calculations. ΔE values are relative to the νOH fundamental frequency at 3694.73 cm−1.

NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.955 
11.66 E1 m(1) 0.953 
33.34 A1 νrock 0.937 
35.02 E1 νXY 0.935 
44.87 E2 νXY + m(1) 0.668 
49.07 E2 m(2) 0.655 
58.04 A1 2νXY 0.909 
12 69.20 E2 2νXY 0.919 
17 79.36 A1 νZ 0.946 
NΔEIrrepAssignmentBasis-state norm
0.00 A1 … 0.955 
11.66 E1 m(1) 0.953 
33.34 A1 νrock 0.937 
35.02 E1 νXY 0.935 
44.87 E2 νXY + m(1) 0.668 
49.07 E2 m(2) 0.655 
58.04 A1 2νXY 0.909 
12 69.20 E2 2νXY 0.919 
17 79.36 A1 ν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 H2O relative to the benzene plane—νXY and νZ, and (ii) those that primarily involve the LAM motions of the H/D atoms of H2O 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 H2O 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 H2O and HDO. For H2O, 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 H2O and HDO stretching modes. For H2O, 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 H2O 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–H2O, 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.

Figure 2 shows computed IR spectra corresponding to the OH stretching region of the H2O–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 E1, 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.

FIG. 2.

Calculated infrared spectra in the region of the OH stretching fundamentals of benzene–H2O. Top: initial state is the ortho ground state. Bottom: initial state is the para ground state. In both spectra, red bars signify parallel bands and blue bars signify perpendicular bands. Labeled bands are referred to in the text.

FIG. 2.

Calculated infrared spectra in the region of the OH stretching fundamentals of benzene–H2O. Top: initial state is the ortho ground state. Bottom: initial state is the para ground state. In both spectra, red bars signify parallel bands and blue bars signify perpendicular bands. Labeled bands are referred to in the text.

Close modal

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), A1 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.

FIG. 3.

Calculated infrared spectrum in the region of the OH stretching fundamental of benzene–HDO. Initial state is the ground state. Red bars signify parallel bands and blue bars signify perpendicular bands. Labeled bands are referred to in the text.

FIG. 3.

Calculated infrared spectrum in the region of the OH stretching fundamental of benzene–HDO. Initial state is the ground state. Red bars signify parallel bands and blue bars signify perpendicular bands. Labeled bands are referred to in the text.

Close modal

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 H2O 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.

Figure 4 shows computed intermolecular Raman spectra for the H2O–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.

FIG. 4.

Calculated intermolecular Raman spectra of the benzene–H2O complex. The initial state for the upper spectrum is the ortho ground state and that for the lower one is the para ground state. Labeled bands are referred to in the text.

FIG. 4.

Calculated intermolecular Raman spectra of the benzene–H2O complex. The initial state for the upper spectrum is the ortho ground state and that for the lower one is the para ground state. Labeled bands are referred to in the text.

Close modal

The measured intermolecular Raman spectrum for the H2O 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 H2O 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 H2O 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).

We have presented an efficient computational methodology for rigorous quantum bound-state calculations of the coupled intra- and intermolecular vibrational states of benzene–H2O and benzene–HDO dimers in which H2O 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 H2O/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 (G12). 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 studies19–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 H2O 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 H2O 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 H2O 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 H2O/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 H2O and HDO complexes are revealed. Thus, plots of the probability density functions show that in the ground intermolecular state of the H2O 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 H2O and HDO in the complexes show that the stretching modes of H2O (ν1 and ν3) and HDO (νOD and νOH) are redshifted, while the H2O/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 H2O 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–H2O and benzene–HDO.

The quantum 9D calculations provide rich data about the effects that the excitation of various intramolecular modes of H2O 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 H2O/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 H2O 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 IR6 and Raman7 spectra of the dimers. The dominant bands in the computed IR spectra of the H2O 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–H2O 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–H2O 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., corannulene43), 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.

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

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.

1.
S.
Suzuki
,
P. G.
Green
,
R. E.
Bumgarner
,
S.
Dasgupta
,
W. A.
Goddard
 III
, and
G. A.
Blake
,
Science
257
,
942
(
1992
).
2.
H. S.
Gutowsky
,
T.
Emilsson
, and
E.
Arunan
,
J. Chem. Phys.
99
,
4883
(
1993
).
3.
A. J.
Gotch
and
T. S.
Zwier
,
J. Chem. Phys.
96
,
3388
(
1992
).
4.
R. N.
Pribble
and
T. S.
Zwier
,
Science
265
,
75
(
1994
).
5.
R. N.
Pribble
and
T. S.
Zwier
,
Faraday Discuss.
97
,
229
(
1994
).
6.
R. N.
Pribble
,
A. W.
Garrett
,
K.
Haber
, and
T. S.
Zwier
,
J. Chem. Phys.
103
,
531
(
1995
).
7.
P. M.
Maxton
,
M. W.
Schaeffer
, and
P. M.
Felker
,
Chem. Phys. Lett.
241
,
603
(
1995
).
8.
A.
Engdahl
and
B.
Nelander
,
J. Phys. Chem.
89
,
2860
(
1985
).
9.
J.
Andersen
,
R. W.
Larsen
,
J.
Ceponkus
,
P.
Uvdal
, and
B.
Nelander
,
J. Phys. Chem. A
124
,
513
(
2020
).
10.
G.
Karlström
,
P.
Linse
,
A.
Wallqvist
, and
B.
Jönsson
,
J. Am. Chem. Soc.
105
,
3777
(
1983
).
11.
S.
Li
,
V. R.
Cooper
,
T.
Thonhauser
,
A.
Puzder
, and
D. C.
Langreth
,
J. Phys. Chem.
112
,
9031
(
2008
).
12.
S. Y.
Fredericks
,
K. D.
Jordan
, and
T. S.
Zwier
,
J. Phys. Chem.
100
,
7810
(
1996
).
13.
M.
Prakash
,
K. G.
Samy
, and
V.
Subramanian
,
J. Phys. Chem. A
113
,
13845
(
2009
).
14.
D. P.
Tabor
,
R.
Kusaka
,
P. S.
Walsh
,
T. S.
Zwier
, and
E. L.
Sibert
 III
,
J. Phys. Chem. A
119
,
9917
(
2015
).
15.
P.
Linse
,
J. Comput. Chem.
9
,
505
(
1988
).
16.
J. K.
Gregory
and
D. C.
Clary
,
Mol. Phys.
88
,
33
(
1996
).
17.
W.
Kim
,
D.
Neuhauser
,
M. R.
Wall
, and
P. M.
Felker
,
J. Chem. Phys.
110
,
8461
(
1999
).
18.
M.
Chaplin
,
Water absorption spectrum
, available at http://www1.lsbu.ac.uk/water/water_vibrational_spectrum.html; 14 January 2019.
19.
D.
Lauvergnat
,
P. M.
Felker
,
Y.
Scribano
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
150
,
154303
(
2019
).
20.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
151
,
024305
(
2019
).
21.
P. M.
Felker
,
D.
Lauvergnat
,
Y.
Scribano
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
151
,
124311
(
2019
).
22.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
152
,
014108
(
2020
).
23.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
117
,
6923
(
2002
).
24.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
119
,
101
(
2003
).
25.
J. C.
Tremblay
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
125
,
094311
(
2006
).
26.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
129
,
234102
(
2009
).
27.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
146
,
104105
(
2017
).
28.
I. I.
Mizus
,
A. A.
Kyuberis
,
N. F.
Zobov
,
V. Y.
Makhnev
,
O. L.
Polyansky
, and
J.
Tennyson
,
Philos. Trans. R. Soc., A
376
,
20170149
(
2018
).
29.
B. T.
Sutcliffe
and
J.
Tennyson
,
Int. J. Quantum Chem.
39
,
183
(
1991
).
30.
R. N.
Zare
,
Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics
(
Wiley-Intersience
,
New York
,
1988
).
31.
V. A.
Mandelshtam
and
H. S.
Taylor
,
J. Chem. Phys.
106
,
5085
(
1997
).
32.
M. R.
Wall
and
D.
Neuhauser
,
J. Chem. Phys.
102
,
8011
(
1995
).
33.

The analogous situation obtains for the H2O@C60 system in regard to a five-fold axis of symmetry. See Supplementary Material, Sec. IV of Ref. 22.

34.
H.
Wei
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
97
,
3029
(
1992
).
35.
J.
Echave
and
D. C.
Clary
,
Chem. Phys. Lett.
190
,
225
(
1992
).
36.
P. M.
Maxton
,
M. W.
Schaeffer
,
S. M.
Ohline
,
W.
Kim
,
V. A.
Venturo
, and
P. M.
Felker
,
J. Chem. Phys.
101
,
8391
(
1994
).
37.
P. M.
Felker
,
D.
Neuhauser
, and
W.
Kim
,
J. Chem. Phys.
114
,
1233
(
2001
).
38.
P. M.
Felker
,
J. Chem. Phys.
114
,
7901
(
2001
).
39.
J. D.
Louck
and
H. W.
Galbraith
,
Rev. Mod. Phys.
48
,
69
(
1976
).
40.
F.
Culot
and
J.
Liévin
,
Phys. Scr.
46
,
502
(
1992
).
41.
M. R.
Battaglia
,
A. D.
Buckingham
, and
J. H.
Williams
,
Chem. Phys. Lett.
78
,
421
(
1981
).
42.
W. F.
Murphy
,
J. Chem. Phys.
67
,
5877
(
1977
).
43.
C.
Pérez
,
A. L.
Steber
,
A. M.
Rijs
,
B.
Temelso
,
G. C.
Shields
,
J. C.
Lopez
,
Z.
Kisiel
, and
M.
Schnell
,
Phys. Chem. Chem. Phys.
19
,
14214
(
2017
).

Supplementary Material