We present a method for efficient calculation of intramolecular vibrational excitations of H_{2}O inside C_{60}, together with the low-energy intermolecular translation-rotation states within each intramolecular vibrational manifold. Apart from assuming rigid C_{60}, this nine-dimensional (9D) quantum treatment is fully coupled. Following the recently introduced approach [P. M. Felker and Z. Bačić, J. Chem. Phys. **151**, 024305 (2019)], the full 9D vibrational Hamiltonian of H_{2}O@C_{60} is partitioned into two reduced-dimension Hamiltonians, a 6D one for the intermolecular vibrations and another in 3D for the intramolecular degrees of freedom, and a 9D remainder term. The two reduced-dimension Hamiltonians are diagonalized, and their eigenvectors are used to build up a product contracted basis in which the full vibrational Hamiltonian is diagonalized. The efficiency of this methodology derives from the insight of our earlier study referenced above that converged high-energy intramolecular vibrational excitations of weakly bound molecular complexes can be obtained from fully coupled quantum calculations where the full-dimensional product contracted basis includes only a small number of intermolecular vibrational eigenstates spanning the range of energies much below those of the intramolecular vibrational states of interest. In this study, the eigenstates included in the 6D intermolecular contacted basis extend to only 410 cm^{−1} above the ground state, which is much less than the H_{2}O stretch and bend fundamentals, at ≈3700 and ≈1600 cm^{−1}, respectively. The 9D calculations predict that the fundamentals of all three intramolecular modes, as well as the bend overtone, of the caged H_{2}O are blueshifted relative to those of the gas-phase H_{2}O, the two stretch modes much more so than the bend. Excitation of the bend mode affects the energies of the low-lying H_{2}O rotational states significantly more than exciting either of the stretching modes. The center-of-mass translational fundamental is virtually unaffected by the excitation of any of the intramolecular vibrational modes. Further progress hinges on the experimental measurement of the vibrational frequency shifts in H_{2}O@C_{60} and *ab initio* calculation of a high-quality 9D potential energy surface for this endohedral complex, neither of which is presently available.

## I. INTRODUCTION

H_{2}O@C_{60} is a remarkable supramolecular complex where a polar molecule, H_{2}O, is encapsulated inside a highly symmetric, nonpolar all-carbon interior of C_{60}. This molecular container provides a unique environment in which the water molecule is unable to form any hydrogen bonds, a highly unusual situation for a species with a strong propensity for hydrogen bonding. H_{2}O@C_{60} is a member of the family of light-molecule endofullerenes (LMEFs), which have molecules such as H_{2}, HF, and CH_{4}, characterized by small masses and large rotational constants, encapsulated inside the cages of C_{60}, C_{70}, and other fullerenes.^{1,2} These inclusion compounds have all been synthesized^{3–6} utilizing the approach known as molecular surgery,^{7–9} in which a hole is created in the fullerene cage through a sequence of chemical reactions, allowing the insertion of the guest molecule into the cavity, followed by reforming the pristine fullerene shell with the guest molecule permanently trapped in its interior.

A conspicuous feature of H_{2}O@C_{60} and other LMEFs, which has motivated a large number of experimental and theoretical studies, is the dominance of quantum effects in the dynamics and spectroscopy of the guest molecules,^{1,2,10} particularly for the low temperatures (typically ranging from 1.5 K to about 30 K) at which the spectroscopic measurements on LMEFs are usually carried out. These species exemplify many of the fundamental, textbook principles of quantum mechanics to the degree and with a clarity that are unmatched among molecular systems amenable to experimental investigations. The quantum effects arise from several sources. One of them is the quantization of the translational center-of-mass (c.m.) degrees of freedom (DOFs) of the encapsulated molecules due to their confinement inside the fullerene (particle-in-a-box effect). Since this confinement is tight and the guest molecules have low molecular masses, the energy differences between the translational eigenstates are large relative to *kT* (*k* being the Boltzmann constant). The same holds for the quantized rotational states of the light molecules due to their large rotational constants. Moreover, the confining potential of the fullerene interior couples the quantized translational and rotational DOFs of the guest molecule, giving rise to a sparse and intricate translation-rotation (TR) energy level structure.^{1,10}

Molecules such as H_{2}O, and also H_{2}, having two symmetrically equivalent H atoms, exhibit the phenomenon of nuclear spin isomerism, which introduces new and prominent quantum features in their TR dynamics inside the fullerene cages. The two identical ^{1}H nuclei are fermions since their nuclear spin is 1/2. In order to satisfy the Pauli principle, the total molecular wave function of the caged molecule, its spatial and spin components, must be antisymmetric with respect to the exchange of the two fermions. This requirement entangles the spin and spatial quantum states in a particular way, leading to nuclear spin isomers, *para* and *ortho*, of H_{2}O (and H_{2}), having total nuclear spins *I* = 0 and 1, respectively. The rotational states of H_{2}O are conventionally labeled with the asymmetric top quantum numbers $jkakc$; for *para*-H_{2}O, *k*_{a} + *k*_{c} has even parity, while for *ortho*-H_{2}O, *k*_{a} + *k*_{c} has odd parity.^{11} The restrictions placed by the Pauli principle on the rotational quantum numbers of the spin isomers of H_{2}O make its already sparse TR level structure inside the fullerenes even sparser and push the TR dynamics deeper in the quantum regime.

H_{2}O@C_{60}, where the fullerene cage encloses a dipolar H_{2}O molecule (just as HF in C_{60}^{5,12}), in its crystalline form constitutes in effect a cubic dipolar lattice.^{13} This provides a unique opportunity for investigating the many-body correlations among the highly quantum molecular dipoles confined in the C_{60} cages and the possible emergence of dipole-ordered phases and ferroelectricity at low temperatures.^{14} Various aspects of the quantized rotational dynamics of H_{2}O@C_{60} have been probed by NMR spectroscopy.^{15–20} The most comprehensive information about the TR eigenstates of H_{2}O@C_{60} has emerged from the low-temperature inelastic neutron scattering (INS) spectra,^{15,21} which exhibit numerous peaks corresponding to the transitions out of the ground states of *para*- and *ortho*-H_{2}O to a broad range of excited TR states of these two nuclear spin isomers.

These experimental findings have motivated a number of theoretical investigations of the quantum TR dynamics and spectroscopy of H_{2}O@C_{60}. We calculated the TR eigenstates of *para*- and *ortho*-H_{2}O in an (isolated) C_{60} cage with *I*_{h} symmetry by means of fully coupled quantum 6D calculations, treating both H_{2}O and C_{60} as rigid.^{22} The computed TR eigenstates allowed tentative assignments of a number of transitions observed in the experimental INS spectra of H_{2}O@C_{60}^{21} that were not assigned previously. What also emerged from these calculations^{22} was that the TR level structure of H_{2}O@C_{60} exhibits all of the key qualitative features that have been identified earlier for H_{2}@C_{60}.^{1,10,23,24} The purely translational eigenstates can be assigned with the quantum numbers of the 3D isotropic harmonic oscillator (HO) suggested by the high symmetry of C_{60}, the principal quantum number *n* = 0, 1, 2, …, and the orbital angular momentum quantum number *l* = *n*, *n* − 2, …, 1 or 0, for odd and even *n*, respectively. For the H_{2}O rotational excitations, the asymmetric top quantum numbers $jkakc$ are employed. In the case of the TR eigenstates that are excited both translationally and rotationally, the orbital angular momentum **l** of the c.m. of H_{2}O and the rotational angular momentum **j** of H_{2}O couple to give the total angular momentum ** λ** =

**l**+

**j**, with

*λ*=

*l*+

*j*,

*l*+

*j*− 1, …, |

*l*−

*j*|. The values of

*l*are those allowed for the quantum number

*n*of the 3D isotropic HO. The TR states having the same quantum numbers

*n*and $jkakc$ are split into as many distinct levels as there are different values of

*λ*, each with the degeneracy 2

*λ*+ 1.

^{22}Prosmiti

*et al.*

^{25}have performed fully coupled 9D quantum calculations of the TR energy levels of H

_{2}O@C

_{60}, for flexible H

_{2}O, but limited to its ground intramolecular vibrational state, and rigid, isolated C

_{60}with

*I*

_{h}symmetry, employing the multiconfiguration time-dependent Hartree (MCTDH) approach. The computed TR levels differed very little from those obtained in the quantum 6D calculations where both H

_{2}O and C

_{60}were taken to be rigid.

^{22}

The INS spectra have also revealed that for *ortho*-H_{2}O inside C_{60}, in solid C_{60} and at low temperatures, its 1_{01} ground state, which would be threefold degenerate if the cage had the *I*_{h} symmetry of isolated C_{60}, is in fact split by about 4.19 cm^{−1} into a nondegenerate lower-energy state and a doubly degenerate higher-energy state.^{15,21,26} The implication of this was that the symmetry of the environment felt by the guest H_{2}O must be lower than *I*_{h}. This lowering of symmetry, commonly referred to as symmetry breaking, was observed earlier at low temperatures in solid H_{2}@C_{60}^{27,28} and subsequently also for solid HF@C_{60}.^{5}

Since our calculations above, and those of Prosmiti *et al.*,^{25} were performed for H_{2}O inside an environment having *I*_{h} symmetry, they could not address the issue of symmetry breaking. For this purpose, we recently developed a first-principles electrostatic quadrupolar model of symmetry breaking,^{29,30} which is free from any adjustable parameters. The magnitude of the splitting and the 1:2 splitting pattern of the 1_{01} ground state of *ortho*-H_{2}O inside C_{60} predicted by the quadrupolar model for the dominant P orientation^{29} are in excellent agreement with experimental results obtained at low temperatures for this system in solid C_{60}.^{15,21,26} The same level of agreement with the experiment regarding the splitting magnitude and pattern was achieved for the ground state of *ortho*-H_{2} and the *j* = 1 rotational level of HF in C_{60}.^{29,30}

An important and readily observable consequence of the encapsulation of H_{2}O in C_{60} for which, conspicuously, no experimental results have been reported so far is the shifts in the frequencies of the H_{2}O intramolecular vibrations away from those of the isolated molecule in the gas phase. These frequency shifts are expected to be appreciable, given that the intramolecular stretch fundamentals measured in the infrared (IR) spectra of H_{2}^{31} and HF^{5} in C_{60} are shifted to lower frequencies, i.e., redshifted, by −98.8 and −170.5 cm^{−1}, respectively. Only a few theoretical treatments of the intramolecular vibrations of H_{2}O in C_{60} and their frequency shifts have been reported. In one, harmonic stretching frequencies of the caged H_{2}O were computed at the Hartree-Fock (HF) level, for the static, optimized structure of H_{2}O.^{32} They gave blueshifts of 41 and 30 cm^{−1} for the symmetric and asymmetric stretch modes of H_{2}O, respectively. In an attempt to go beyond the static treatment, the IR spectrum of H_{2}O@C_{60} was computed by a mixed quantum-classical approach,^{33} where classical molecular-dynamics (MD) simulations were performed first in order to generate an ensemble of configurations averaged over the intermolecular, translational, and rotational motions of H_{2}O inside the fullerene. Then, at each configuration, the 3D vibrational levels of H_{2}O were calculated quantum mechanically, thus generating the lineshape of the IR spectrum. The symmetric stretch mode was found to be blueshifted by 13–19 cm^{−1}, while the asymmetric stretch mode exhibited a small blueshift or none at all.^{33} Finally, the frequency shifts of the vibrational modes of H_{2}O inside C_{60} were computed by means of classical MD simulations,^{34} yielding blueshifts of 14 and 15 cm^{−1} for the symmetric and asymmetric stretch modes, respectively.

The blueshifts of the H_{2}O stretching modes predicted by the three theoretical treatments^{32–34} are somewhat surprising, in view of the substantial redshifts observed for both H_{2} and HF in C_{60}.^{5,31} It is conceivable that the larger size of H_{2}O, relative to H_{2} and HF, can result in increased repulsive interaction with the C_{60} cage, leading to the blueshift of the stretching frequencies. However, it should be kept in mind that the level of *ab initio* theory that yielded blueshifts of the vibrational stretches for H_{2}O in C_{60}^{32} also gave a slight blueshift for HF in C_{60}, in disagreement with the substantial redshift (−170.5 cm^{−1}) measured for this endohedral complex.^{5} Therefore, whether the vibrational modes of H_{2}O in C_{60} are blue- or redshifted is a question that cannot be considered fully settled yet.

As emphasized earlier in this section, all DOFs of H_{2}O inside C_{60}, those pertaining to its intramolecular vibrations, (weakly) hindered rotation, and the c.m. translation, are quantized, and the effects of this, and of the coupling among the DOFs, leave clear fingerprints in the INS and NMR spectra. A rigorous approach to the intramolecular vibrational excitations of H_{2}O, and their frequency shifts, coupled to the intermolecular TR motions, requires fully coupled 9D quantum bound-state calculations for flexible H_{2}O and rigid C_{60} extending to the energies of the H_{2}O vibrational fundamentals and overtones. The earlier quantum treatments of the TR dynamics of H_{2}O@C_{60}^{22,25,29,30,35} did not, and could not, accomplish this either because H_{2}O was treated as rigid^{22,29,30,35} or because the calculations were restricted to H_{2}O in the ground intramolecular vibrational state.^{25}

The frequencies of the intramolecular vibrational fundamentals of H_{2}O [3657.05, 3755.93, and 1594.74 cm^{−1}, for the symmetric (*ν*_{1}) and asymmetric stretch (*ν*_{3}) modes and the bend (*ν*_{2}), respectively, in the gas phase^{36}] are at least an order of magnitude higher than those calculated by us^{22} for the intermolecular modes, c.m. translation (162.1 cm^{−1}), and lowest rotational levels of *ortho*- and *para*-H_{2}O in C_{60}, 23.80 and 37.16 cm^{−1}, respectively. Such a large disparity between the frequencies of the intramolecular and intermolecular vibrational modes is typical for weakly bound, hydrogen-bonded (HB) and van der Waals (vdW) molecular dimers (H_{2}O@C_{60} can be viewed as a special case, where one of the monomers, C_{60}, encloses the other, H_{2}O). As a result, hundreds, if not a 1000 or more, of highly excited TR states in the intramolecular ground-state manifold lie below the energies of the H_{2}O intramolecular vibrational excitations considered. Prior to our recent work,^{12,37,38} it has been generally assumed that accurate full-dimensional and fully coupled quantum bound-state calculations of weakly bound molecular dimers encompassing the intramolecular vibrational fundamentals and overtones of the monomers require converging either all intermolecular vibrational eigenstates lying below the energies of these intramolecular excitations or at least a dense set of highly excited intermolecular vibrational states close to them in energy. This makes such calculations extremely difficult, and in many instances prohibitively costly, especially if the monomers are polyatomic.

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^{37} revealed that this widely held view is not correct, making the task much easier than previously thought. In that study, we demonstrated that for accurate computation of the first excited ($v$ = 1) vibrational state of the caged H_{2} at ≈4100 cm^{−1}, it sufficed to have converged only a modest number of TR states in the $v$ = 0 zero manifold up to at most 400–450 cm^{−1} above the ground state and none within several thousand wave numbers of the H_{2} intramolecular fundamental. We explained this unexpected finding in terms of extremely weak coupling between the intramolecular stretch of the guest H_{2} and the large number of highly excited $v$ = 0 intermolecular TR states in its vicinity. Consequently, not including the latter in the calculations affects negligibly the accuracy with which the H_{2} intramolecular vibrational excitation is calculated. This greatly reduced the dimension of the basis set, transforming what would be a formidable computation into a readily tractable one.^{37}

This insight led us to the development of a very efficient method for full-dimensional and fully coupled quantum calculations of intramolecular vibrational excitations of weakly bound molecular dimers, together with the intermolecular vibrational states in each intramolecular vibrational manifold.^{38} In this approach, the vibrational Hamiltonian of the dimer is partitioned into two reduced-dimension Hamiltonians, a rigid-monomer one for the intermolecular vibrations and the other for all intramolecular vibrational degrees of freedom, and a remainder. The eigenstates of the two reduced-dimension Hamiltonians are used to build up a product contracted basis for the diagonalization of the full vibrational Hamiltonian. Such a basis is particularly convenient for taking advantage of what we learned previously about the weak coupling between the intramolecular and intermolecular degrees^{37} as it enables including only a certain number (small, as it turns out) of lowest-energy intramolecular and intermolecular vibrational eigenstates in the final full-dimensional direct-product basis.^{38} The application of this approach to the 6D vibration-rotation-tunneling eigenstates of (HF)_{2}, a prototypical HB dimer, gave results in excellent agreement with those in the literature while using just a fraction of the basis-set size required by the other methods.^{38} 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 cage,^{39} which would otherwise be virtually intractable.

Our idea to divide the coordinates of molecular dimers into two groups, intramolecular and intermolecular, respectively, generate contracted eigenvector bases for both, and use them to form the final product contracted basis has similarities with the work of Carrington and co-workers on the vibrational levels of polyatomic molecules, such as H_{2}O_{2}, CH_{4}, vinylidene, and $CH5+$.^{40–43} They defined two reduced-dimension Hamiltonians, for the stretch and bend groups of coordinates, respectively, by freezing the coordinates in the other group at their equilibrium values. Their eigenvectors served as the contracted stretch and bend basis functions, respectively. The strong stretch-bend coupling in these molecules precluded the possibility of including 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 this paper, we present a computational methodology, based on the approach outlined above,^{38} for performing, with modest effort, fully coupled 9D quantum calculations of the intramolecular vibrational frequencies of a flexible H_{2}O molecule in (rigid) C_{60}, together with the intermolecular TR states within each intramolecular vibrational manifold considered. In this treatment, the effects of the averaging over the quantized TR motions on the H_{2}O intramolecular vibrational fundamentals (and overtone, for the *ν*_{2} mode) and their frequency shifts are accounted for rigorously. Careful symmetry considerations result in reducing the matrix representation of the full 9D Hamiltonian to a block-diagonal form, of which ten symmetry blocks need to be diagonalized separately. These matrices are small, the largest having the dimension of 816, primarily because only a very small number of the lowest-energy (symmetry-adapted) eigenstates of the reduced-dimension 6D intermolecular Hamiltonian, at most 24, needed to be included in the 9D bases for the different symmetry blocks. The energies of these 6D eigenstates extend to only 410 cm^{−1} above the ground state, far below the energies of the H_{2}O stretch and bend fundamentals. No *ab initio*-calculated 9D potential energy surface (PES) for H_{2}O inside C_{60} is available (not even in 6D, for rigid H_{2}O). Therefore, the 9D PES employed in these calculations is obtained by combining the high-accuracy 3D PES for the isolated H_{2}O molecule^{36} with the intermolecular PES represented by pairwise-additive atom-atom interactions between H_{2}O and C_{60} parameterized by Aluru and co-workers.^{34} The latter was used by us previously in the quantum 6D calculations of the TR eigenstates of rigid H_{2}O in (rigid) C_{60}.^{22}

## II. COMPUTATIONAL METHODOLOGY

### A. Hamiltonian

We consider a flexible H_{2}O molecule encaged in an isolated and rigid C_{60} with *I*_{h} symmetry. It would be straightforward, at least formally, to include the effects of symmetry breaking^{29,30} on the intramolecular vibrational excitations. However, we suspect that the dominant contribution to this would come from the dependence of the quadrupole moment of H_{2}O on the intramolecular coordinates, which is not known. Therefore, we postpone the inclusion of symmetry breaking for now.

In this case, the 9D vibrational Hamiltonian of H_{2}O@C_{60} can be written as

where *M* is the mass of H_{2}O. The nine coordinates on which this operator depends are as follows. First, there are three coordinates associated with the position vector, **R**, of the water center of mass (c.m.) measured with respect to a cage-fixed (CF) Cartesian axis system ($X^,\u0176,Z^$) with origin at the center of the C_{60} cage. We choose these axes such that *Ẑ* is along a *C*_{5} symmetry axis of the icosahedral cage and *Ŷ* is along one of the *C*_{2} symmetry axes, $X^=\u0176\xd7Z^$. For the components of **R**, we work with the spherical coordinates (*R*, *β*, *α*), with *R* ≡ |**R**|, $cos\u2061\beta \u2261R^\u22c5Z^$, and $(cos\u2061\alpha ,sin\u2061\alpha )\u2261(R^\u22c5X^,R^\u22c5\u0176)/sin\u2061\beta $. Second, the three coordinates represented by *ω* are the Euler angles (*ϕ*, *θ*, *χ*) that fix the orientation of the water moiety’s molecule-fixed (MF) axis system relative to the CF frame. Third, the three coordinates represented by **q** are the vibrational coordinates of the water moiety. For these, we choose the Radau coordinates (*R*_{1}, *R*_{2}, Θ).^{44–46} We also use Radau bisector-z embedding^{47} to define the water MF frame, with the MF *ŷ* axis taken to be parallel to the cross product of Radau vectors **R**_{1} ×**R**_{2} (i.e., *ŷ* is perpendicular to the molecular plane) and $x^=\u0177\xd7z^$.

The terms in *Ĥ* correspond, respectively, from left to right in Eq. (1), to the kinetic energy of the c.m. translational motion of the water, the vibration-rotation kinetic energy of the water $(T^vr)$, the Coriolis contribution to the kinetic energy of the water $(T^cor)$, the vibrational kinetic energy of the water $(T^v)$, and the 9D PES of flexible H_{2}O inside C_{60} (*V*_{tot}). This potential term can be written as

where the first term represents the interaction between the water moiety and the interior of C_{60}, while the second term is the 3D PES for the isolated, gas-phase water molecule. The 9D PES used by Prosmiti *et al.*^{25} in their 9D quantum calculations of H_{2}O@C_{60} was constructed in the same manner.

The explicit forms of $T^vr$, $T^cor$, and $T^v$ are given by^{47}

and

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

is a Hermitian operator connected with the angular momentum of the bending vibration in water.^{47}

As to the 9D PES *V*_{tot} in Eq. (2), $VH2O-cage$ is represented by the sum over the pairwise interactions, modeled using the Lennard-Jones (LJ) 12-6 potentials, of each atom of H_{2}O with each atom of C_{60}, with the parameterization (MD1) by Aluru *et al.*,^{34} and the C_{60} nuclear coordinates assumed in our previous work.^{22} For *V*_{intra}, we use the high-quality 3D PES of the isolated H_{2}O (and the computer code to generate it) from the work of Mizus *et al.*^{36}

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

To solve for the eigenstates of *Ĥ*, we use the general procedure described by us elsewhere.^{38} Briefly, we divide the operator up into three parts: a 6D intermolecular Hamiltonian *Ĥ*_{inter}, a 3D intramolecular Hamiltonian *Ĥ*_{intra}, and a 9D remainder term Δ*Ĥ*. We then solve for the low-energy eigenfunctions, |*κ*⟩, and eigenvalues, $E\kappa inter$, of *Ĥ*_{inter} and those, |*γ*⟩ and $E\gamma intra$, of *Ĥ*_{intra}. From these reduced-dimension eigenfunctions, we construct a contracted product basis composed of functions of the form |*κ*, *γ*⟩ ≡ |*κ*⟩|*γ*⟩ and then use that basis to diagonalize the full *Ĥ*.

We define

where **q**_{0} ≡ (*R*_{1,0}, *R*_{2,0}, Θ_{0}) are constants set to values of the water vibrational coordinates close to those of the equilibrium geometry of the monomer, and *V*_{inter} is defined as

We also define

where the “C_{60}-adapted” 3D intramolecular potential, *V*_{adap}, is defined as the proton-exchange-symmetrized version of

Here, **R**_{0} and *ω*_{0} are fixed values of **R** and *ω*. By choosing these values to correspond to the minimum of $VH2O-cage$ when the Radau coordinates are fixed to (*R*_{1,0}, *R*_{2,0}, Θ_{0}), one obtains a 3D H_{2}O vibrational potential that includes some of the effects (perhaps the bulk thereof) of the H_{2}O–C_{60} interaction on the vibrational PES of H_{2}O.

where

and

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

The main task in diagonalizing *Ĥ* is computing the matrix elements of Δ*Ĥ*.

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

#### 1. Primitive intermolecular basis

For the problem of diagonalizing *Ĥ*_{inter}, we employ a primitive basis composed of functions of the form

where |*n*, *l*, *m*_{l}⟩ are normalized eigenfunctions of the isotropic 3D harmonic oscillator Hamiltonian^{48,49} and |*j*, *m*_{j}, *k*⟩ are symmetric-top rotational eigenfunctions.^{50} The former are of the form

(*l* = 0, 1, 2 … ; *n* = *l*, *l* + 2, …, *n*_{max}; |*m*_{l}| = 0, 1, …, *l*), where the radial function, $Fn,l(M\Omega )(R)$, depends on the choice of the product of the oscillator mass (*M*) and its angular frequency (Ω), and $Ylml$ is a spherical harmonic. The latter are normalized Wigner rotation matrix elements,

#### 2. Symmetry considerations

The basis of Eq. (15) and the particular choice that we have made for the CF axes allow for a straightforward classification of intermolecular basis states (and intermolecular eigenstates) with respect to symmetry. In particular, rotation of the water moiety by *p*2*π*/5 radians (*p* an integer) about *Ẑ* transforms *α* → *α* + *p*2*π*/5 and *ϕ* → *ϕ* + *p*2*π*/5 while leaving *Ĥ*_{inter} (and *Ĥ*) unchanged. Thus, only those |*n*, *l*, *m*_{l}, *j*, *m*_{j}, *k*⟩ states that have the same symmetry with respect to such rotations are capable of being coupled by *Ĥ*_{inter} (or *Ĥ*). It is easy to see that such states must have the same value of (*m*_{l} + *m*_{j}) mod 5. Hence, we factorize the intermolecular basis into five different sets characterized by values of (*m*_{l} + *m*_{j}) mod 5 = 0, ±1, and ±2. We solve the 6D intermolecular problem by diagonalizing separately the blocks of the *Ĥ*_{inter} matrix corresponding to these different sets. We also make use of the same type of symmetry factorization in Sec. II E to block diagonalize *Ĥ* in the |*κ*⟩|*γ*⟩ basis. This is possible because |*κ*⟩ have well-defined symmetry with respect to the *p*2*π*/5 *Ẑ*-axis rotations and *Ĥ* is invariant with respect to them.

The other symmetry that we make use of to block diagonalize the *Ĥ*_{inter} matrix (and that of *Ĥ*) is parity, the symmetric or antisymmetric transformation of states due to inversion of the coordinates of the water particles though the CF origin. The states of the intermolecular basis transform upon inversion, *Ê*^{*}, as

Rather than building parity directly into the primitive basis, we build it into the random initial intermolecular state vector that is used in our iterative matrix-on-vector diagonalization routine described in the computational details below. We do this by applying a parity projection operator to a random state vector expressed in the intermolecular basis. We then proceed with the diagonalization routine by using the parity-filtered result as the initial state vector in that process. In this way, we can choose to compute states of even or odd parity.

It is important to point out that the symmetry factorization effected by sorting with respect to (*m*_{l} + *m*_{j}) and by parity-filtering *does not* in a given calculation produce intermolecular eigenstates that all transform according to the same irreducible representation (irrep) of *I*_{h}. This is because these symmetry operations are elements of a subgroup (*D*_{5d}) of *I*_{h}, and, in general, an irrep of *I*_{h} correlates with several different irreps of *D*_{5d}. Thus, for example, an even-parity state with (*m*_{l} + *m*_{j}) mod 5 = 0 belongs to either the *A*_{1g} or *A*_{2g} irrep of *D*_{5d}. However, an *A*_{1g} state in *D*_{5d} could be either an *A*_{g} or an *H*_{g} state in *I*_{h}. An *A*_{2g} state in *D*_{5d} could be either a *T*_{1g} or a *T*_{2g} state in *I*_{h}. The point is that the calculations employing even(odd)-parity, (*m*_{l} + *m*_{j}) mod 5 = 0 basis states yield eigenstates that belong to the *A*_{g}(*A*_{u}), *T*_{1g}(*T*_{1u}), *T*_{2g}(*T*_{2u}), and *H*_{g}(*H*_{u}) irreps of *I*_{h}. Similarly, calculations employing even(odd)-parity, (*m*_{l} + *m*_{j}) mod 5 = ±1 basis states yield eigenstates that belong to the *I*_{h} irreps *T*_{1g}(*T*_{1u}), *G*_{g}(*G*_{u}), and *H*_{g}(*H*_{u}). Finally, even/odd-parity basis states that have (*m*_{l} + *m*_{j}) mod 5 = ±2 produce eigenstates belonging to the *T*_{2g}(*T*_{2u}), *G*_{g}(*G*_{u}), and *H*_{g}(*H*_{u}) irreps of *I*_{h}. One sees that it is not necessary to do calculations for all five (*m*_{l} + *m*_{j}) symmetry-factored basis-state sets to get the energies and *I*_{h}-symmetry assignments of all the levels in a given energy range: It is only necessary to do three such calculations, for example, for (*m*_{l} + *m*_{j}) mod 5 = 0, +1, and +2.

The above also has implications in respect to the full 9D calculations. In constructing the 9D bases, one can effect a symmetry factorization by grouping together those intermolecular eigenstates that (a) belong to the same irrep of *I*_{h} *and* (b) are composed of basis states that have the same value of (*m*_{l} + *m*_{j}) mod 5. The resulting |*κ*, *γ*⟩ 9D basis states that correspond to each such set can only be coupled by *Ĥ* to states in the same set. Thus, the matrix of *Ĥ* can be block-diagonalized into 32 blocks with the 9D basis that we use. Of these 32, only ten need to be set up and diagonalized to obtain the energies of the 9D levels corresponding to all of the *I*_{h} irreps.

#### 3. Computational details

To diagonalize *Ĥ*_{inter}, we make use of the Chebyshev version^{51} of filter diagonalization.^{52} The main task in this matrix-on-vector iterative method is the repeated application of the Hamiltonian onto a random initial state vector,

where *I* is a meta-index denoting a particular set of the indices *n*, *l*, *m*_{l}, *j*, *m*_{j}, *k*. Such operations take the form

We break these operations up into two parts by treating the kinetic- and potential-energy parts of *Ĥ*_{inter} separately. The kinetic-energy matrix elements

can all be readily evaluated analytically (supplementary material). Hence, the kinetic-energy portion of Eq. (20) is straightforward to handle by matrix-on-vector multiplication.

To operate with the potential energy, *V*_{inter}(**R**, *ω*; **q**_{0}), we transform the state function coordinate-by-coordinate to a 6D grid representation (*R*_{p}, (cos *β*)_{q}, (cos *θ*)_{r}, *α*_{s}, *ϕ*_{t}, *χ*_{u}). We then multiply it at each grid point by the value of *V*_{inter} at that point. Finally, we transform the result back to the |*n*, *l*, *m*_{l}, *j*, *m*_{j}, *k*⟩ basis. The grid we use consists of the following. $Rp=up/(M\Omega )$, *p* = 1, …, *N*_{p}, where *u*_{p} are Gauss-associated-Laguerre quadrature points. (cos *β*)_{q}, *q* = 1, …, *N*_{q}, and (cos *θ*)_{r}, *r* = 1, …, *N*_{r}, are Gauss-Legendre quadrature points. The 3D portion of the grid associated with the azimuthal angles is a Fourier grid. For that grid the full range of points for *ϕ*_{t} = 2*π*(*t* − 1)/*N*_{t}, *t* = 1, …, *N*_{t}, and *χ*_{u} = 2*π*(*u* − 1)/*N*_{u}, *u* = 1, …, *N*_{u}, is included. However, owing to the *C*_{5} symmetry about the *Ẑ* axis, one need only include one-fifth of the full range of Fourier points for the *α* angle: *α*_{s} = 2*π*(*s* − 1)/*N*_{s}, *s* = 1, …, *N*_{s}/5, with *N*_{s} chosen to be a multiple of 5 (supplementary material). Thus, the full 6D grid comprises 1/5 × *N*_{p}*N*_{q}*N*_{r}*N*_{s}*N*_{t}*N*_{u} points. Full details pertaining to the procedure by which operation with *V*_{inter} is effected are provided (supplementary material).

One important point concerning Eq. (20) is that the matrix elements of *Ĥ*_{inter} in the primitive basis that we have chosen are all real numbers, despite the basis functions themselves being complex. The proof of this is straightforward. The analytic expressions for the kinetic-energy matrix elements show them to be real. As to those of *V*_{inter}, the *I*_{h} symmetry of the cage and our choice of CF axes require that they be real, as well (supplementary material). The consequence of the fact that ⟨*I*′|*Ĥ*_{inter}|*I*⟩ are real is that if one chooses the initial state-vector amplitudes ⟨*I*|*ψ*_{inter}⟩ to be real, the analogous amplitudes of the state vectors obtained by successive operations of *Ĥ*_{inter} will remain real. This reduces by about a factor 2 the computational effort required to affect the diagonalization of *Ĥ*_{inter}. It also produces intermolecular eigenvectors, |*κ*⟩, composed of real rather than complex elements. This property of |*κ*⟩, in turn, reduces significantly the effort required in the diagonalization of *Ĥ* relative to what it would be if complex amplitudes were involved.

The specific basis-set and grid parameters that we used in the diagonalization of *Ĥ*_{inter} are given in Table I. All possible basis functions corresponding to *n* = 0, 1, …, *n*_{max} and *j* = 0, 1, …, *j*_{max} were included for each value of (*m*_{l} + *m*_{j}) mod 5. The total size of the basis for each such value is given in Table II. The molecular parameters relevant to *Ĥ*_{inter} were taken as *M* = 32 838.683 a.u. (the mass of the water molecule), *m*_{H} = 1836.153 a.u. (the mass of each of the H nuclei), *m*_{O} = 29 156.377 a.u. (the mass of the oxygen nucleus), and *R*_{1,0} = *R*_{2,0} = 1.828 367 355 bohrs, cos Θ_{0} = −0.305 438 6 (the fixed Radau coordinates). These Radau coordinates correspond to a rigid water geometry in which the OH bond lengths are each 0.960 Å, and the HOH bond angle is 104.520°. [Note that the Radau-coordinate convention that we use is that of Sutcliffe and Tennyson (S&T).^{46} An entirely equivalent formulation^{44} defines the radial coordinates and their associated *μ*_{i} differently than do S&T. As such, for the same rigid water geometry that we use, the radial coordinates and mass factors in that alternative formulation would have different numerical values than the ones that we quote here.]

n_{max}
. | MΩ (a.u.)
. | j_{max}
. | N_{p} (#R)
. | N_{q} (#β)
. | N_{r} (#θ)
. | N_{s} (5 × #α)
. | N_{t} (#ϕ)
. | N_{u} (#χ)
. |
---|---|---|---|---|---|---|---|---|

8 | 24.38 | 9 | 12 | 10 | 10 | 20 | 20 | 20 |

n_{max}
. | MΩ (a.u.)
. | j_{max}
. | N_{p} (#R)
. | N_{q} (#β)
. | N_{r} (#θ)
. | N_{s} (5 × #α)
. | N_{t} (#ϕ)
. | N_{u} (#χ)
. |
---|---|---|---|---|---|---|---|---|

8 | 24.38 | 9 | 12 | 10 | 10 | 20 | 20 | 20 |

(m_{l} + m_{j}) mod 5 | 0 | ±1 | ±2 |

No. of |n, l, m_{l}, j, m_{j}, k⟩ | 43 794 | 43 874 | 43 954 |

(m_{l} + m_{j}) mod 5 | 0 | ±1 | ±2 |

No. of |n, l, m_{l}, j, m_{j}, k⟩ | 43 794 | 43 874 | 43 954 |

The second column of Table III presents limited results for some of the lowest-energy eigenstates of *Ĥ*_{inter}. These can be compared with those from our previous work^{22} (column three of Table III) on the rigid-monomer intermolecular states of H_{2}O@C_{60}, from which we take the assignments listed in the rightmost column. Any small differences between this work and the earlier one can be attributed to slightly different rigid-water geometries and rotational constants. The irreps listed in column five of Table III correspond to the $Ih(12)$ group,^{53} which adds proton-exchange symmetry (*s* = symmetric and *a* = antisymmetric with respect to proton exchange) to the operations of *I*_{h}. We emphasize that we have not built exchange symmetry into the primitive intermolecular basis in this work but instead have analyzed the computed eigenvectors of *Ĥ*_{inter} to ascertain their transformation properties due to exchange. A complete compilation of results pertaining to states obtained from the diagonalization of *Ĥ*_{inter} and used subsequently to construct the 9D bases is available (supplementary material).

N
. | ΔE_{6D} (this work)
. | ΔE_{6D} (Ref. 22)
. | ΔE_{9D} (this work)
. | Irrep . | Assignment . |
---|---|---|---|---|---|

1 | 0.0^{a} | 0.0^{b} | 0.0^{c} | $Ag(s)$ | Ground state |

2 | 23.94 | 23.74 | 23.82 | $T1u(a)$ | 1_{01} |

3 | 36.13 | 36.50 | 36.73 | $T1u(s)$ | 1_{11} |

4 | 41.30 | 41.85 | 42.03 | $T1g(a)$ | 1_{10} |

5 | 70.47 | 69.85 | 70.08 | $Hg(s)$ | 2_{02} |

6 | 78.75 | 78.55 | 78.99 | $Hg(a)$ | 2_{12} |

7 | 94.27 | 94.61 | 94.91 | $Hu(s)$ | 2_{11} |

8 | 130.04 | 132.0 | 132.67 | $Hu(a)$ | 2_{21} |

9 | 131.29 | 133.27 | 133.94 | $Hg(s)$ | 2_{20} |

10 | 136.90 | 135.57 | 135.99 | $T2u(a)$ | 3_{03}(a) |

11 | 137.73 | 136.38 | 136.94 | $Gu(a)$ | 3_{03}(b) |

14 | 162.11 | 162.08 | 163.04 | $T1u(s)$ | n = 1 |

N
. | ΔE_{6D} (this work)
. | ΔE_{6D} (Ref. 22)
. | ΔE_{9D} (this work)
. | Irrep . | Assignment . |
---|---|---|---|---|---|

1 | 0.0^{a} | 0.0^{b} | 0.0^{c} | $Ag(s)$ | Ground state |

2 | 23.94 | 23.74 | 23.82 | $T1u(a)$ | 1_{01} |

3 | 36.13 | 36.50 | 36.73 | $T1u(s)$ | 1_{11} |

4 | 41.30 | 41.85 | 42.03 | $T1g(a)$ | 1_{10} |

5 | 70.47 | 69.85 | 70.08 | $Hg(s)$ | 2_{02} |

6 | 78.75 | 78.55 | 78.99 | $Hg(a)$ | 2_{12} |

7 | 94.27 | 94.61 | 94.91 | $Hu(s)$ | 2_{11} |

8 | 130.04 | 132.0 | 132.67 | $Hu(a)$ | 2_{21} |

9 | 131.29 | 133.27 | 133.94 | $Hg(s)$ | 2_{20} |

10 | 136.90 | 135.57 | 135.99 | $T2u(a)$ | 3_{03}(a) |

11 | 137.73 | 136.38 | 136.94 | $Gu(a)$ | 3_{03}(b) |

14 | 162.11 | 162.08 | 163.04 | $T1u(s)$ | n = 1 |

^{a}

*E*_{0} = −1986.31 cm^{−1} relative to separated rigid C_{60} and H_{2}O.

^{b}

*E*_{0} = −1987.57 cm^{−1} relative to separated rigid C_{60} and H_{2}O.

^{c}

*E*_{0} = −1966.64 cm^{−1} relative to separated rigid C_{60} and flexible H_{2}O.

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

#### 1. Intramolecular basis

For the intramolecular problem, we use a product basis of three 1D discrete variable representation (DVR) functions covering each of the vibrational coordinates. These are of the form

The radial DVRs, $|Ri,\eta i\u2009$, are tridiagonal Morse DVRs^{54} generated by using a Morse potential

with the mass factor *μ*_{i} = 1731.932 a.u., *D* = 0.25 hartree, *a* = 0.90 bohr^{−1}, and *R*_{e} = 1.8 bohrs. The bend DVR, $|(cos\Theta )\eta 3\u2009$, is a potential-optimized DVR (PODVR)^{55} generated by first diagonalizing $T^v+Vintra$ for fixed values of *R*_{1} = *R*_{2} = 1.828 367 355 5 bohrs in a basis composed of 40 Gauss-Legendre DVRs covering the Θ coordinate. In the basis of the resulting 1D eigenvectors, the matrix of cos Θ was then constructed and diagonalized to obtain the PODVR and relevant operator matrices in the PODVR basis. The ultimate intramolecular basis that we employed consists of two identical DVRs of nine functions apiece for the *R*_{i} coordinates and 15 PODVR functions for the Θ coordinate, for a total of 1215 functions.

#### 2. Computational details

We solved for the eigenstates of *Ĥ*_{intra} [Eq. (9)] for two different potentials. First, we took *V*_{adap} = *V*_{intra} [i.e., we chose **R**_{0} = *∞* in Eq. (10)]. We did this so as to allow us to check our procedure by comparing results with those of Mizus *et al.*^{36} for isolated H_{2}O. Second, we chose **R**_{0} = 0 and *ω*_{0} = (*ϕ*_{0}, *θ*_{0}, *χ*_{0}) = (0, 0, 0) in computing *V*_{adap} from Eq. (10). These values for the **R** and *ω* coordinates correspond to those at, or very close to, the minimum of *V*_{inter} for the *R*_{1,0}, *R*_{2,0}, and Θ_{0} values that we use. The intramolecular eigenvectors corresponding to this C_{60}-adapted intramolecular potential are those that we ultimately use in the 9D basis.

For both versions of *V*_{adap}, we diagonalized *Ĥ*_{intra} by using the Chebyshev version of filter diagonalization. Operation with the potential on state functions expressed in the basis of Eq. (22) is trivial given that the *V*_{adap} matrix is diagonal in that basis with matrix elements equal to the value of the potential function evaluated at the corresponding 3D DVR quadrature point. To operate with the kinetic-energy portion of *Ĥ*_{intra}, one needs to evaluate the matrix elements,

and

The 1D matrix elements appearing in Eq. (24) are readily obtained by evaluating analytically the matrix elements of the radial kinetic energy operator in the tridiagonal Morse basis from which the |*η*_{i}⟩ radial DVRs are derived^{54} and then transforming that matrix to the DVR representation. Similarly, matrix elements of the bend operator in Eq. (25) are easily obtained analytically in a basis of ordinary Legendre functions.^{47} Transformation of that operator matrix to the PODVR representation is straightforward with the Legendre-to-PODVR transformation matrix in hand.

The results pertaining to all of the intramolecular states used in the construction of the 9D basis are tabulated in their entirety (supplementary material).

### E. Diagonalization of *Ĥ*

With the intermolecular and intramolecular eigenstates computed, we construct the symmetry-adapted 9D bases |*κ*, *γ*⟩ by including the 34 lowest-energy intramolecular states, |*γ*⟩, corresponding to Δ*E*^{intra} ≤ 11 300 cm^{−1}, and all |*κ*⟩ with the appropriate symmetry having Δ*E*^{inter} ≤ 410 cm^{−1}. By “appropriate symmetry,” we refer to the way in which the |*κ*⟩ transform with respect to (a) the operations of *I*_{h} and (b) the specific *C*_{5} operations corresponding to rotations about the *Ẑ* axis. The basis-set sizes for select symmetry blocks are given in Table IV. Since the basis sets for all symmetry blocks include *N*_{intra} = 34 intramolecular states |*γ*⟩, the number of intermolecular states |*κ*⟩ (*N*_{inter}) in the basis is obtained by dividing the basis-set size for each symmetry block in Table IV by 34.

(m_{l} + m_{j}), mod 5:(I_{h} irreps): | 0: (A_{g}, T_{1g}, T_{2g}, H_{g}) | +1: G_{g} | 0: (A_{u}, T_{1u}, T_{2u}, H_{u}) | +1:(G_{u}) |

# of |κ, γ⟩ | (170, 340, 408, 816) | 646 | (34, 510, 578, 714) | 680 |

(m_{l} + m_{j}), mod 5:(I_{h} irreps): | 0: (A_{g}, T_{1g}, T_{2g}, H_{g}) | +1: G_{g} | 0: (A_{u}, T_{1u}, T_{2u}, H_{u}) | +1:(G_{u}) |

# of |κ, γ⟩ | (170, 340, 408, 816) | 646 | (34, 510, 578, 714) | 680 |

For any given symmetry block in the 9D basis, the major task in diagonalizing *Ĥ* is to compute the matrix elements,

To evaluate the matrix elements of the first two terms on the rhs of Eq. (26) our general strategy derives from the fact that the relevant operators are diagonal in the |*η*⟩ ≡ |*η*_{1}*η*_{2}*η*_{3}⟩ intramolecular primitive basis. As such, defining *Ô* ≡ Δ$Tvr$ + Δ*V*, we evaluate first the quantities ⟨*κ*′, *η*|*Ô*|*κ*, *η*⟩ (similar to the *F* matrix of Carrington^{56}) for all *κ*′, *κ*, and *η* and then transform them to the |*κ*, *γ*⟩ basis,

where ⟨*η*|*γ*⟩ are available from the diagonalization of *Ĥ*_{intra}. ⟨*κ*′, *η*|*Ô*|*κ*, *η*⟩ are evaluated as follows. We first compute

for a given *η* and *κ*. This step is essentially equivalent to the iterative step in the diagonalization of *Ĥ*_{inter}, i.e., Eq. (20), and is carried out in a completely analogous way. We then contract *Ô*|*κ*, *η*⟩ with the bras ⟨*κ*′, *η*| to obtain ⟨*κ*′, *η*|*Ô*|*κ*, *η*⟩ for all *κ*′ and for the initial *κ* and *η*. We repeat the calculation of Eq. (28) and the subsequent contractions for all *κ* values until ⟨*κ*′, *η*|*Ô*|*κ*, *η*⟩ are computed for all *κ*′, *κ*, and the initial *η* value. Finally, we repeat this entire process for every value of *η*. Note that for each value of *η* in this procedure, the function Δ*V*(**R**, *ω*, *η*) has to be computed on a 6D (**R**, *ω*) grid. We use the same grid for this that we employ in the diagonalization of *Ĥ*_{inter} (Sec. II C 3). This function need only be stored in memory while Eq. (28) is being computed for the pertinent *η* value. In this way, storage of Δ*V* on a full 9D grid is avoided.

The calculation of $\u27e8\kappa \u2032,\gamma \u2032|T^cor(\omega ,q)|\kappa ,\gamma \u27e9$ matrix elements is slightly more complicated than the above because the operator is diagonal in *η*_{1} and *η*_{2} but not in *η*_{3}. We use the following procedure. We first evaluate the quantities $\u27e8\kappa \u2032,\eta 1\eta 2\eta 3\u2032|T^cor|\kappa ,\eta 1\eta 2\eta 3\u27e9$ for all *κ*′, *κ*, $\eta 3\u2032$, *η*_{3}, *η*_{1}, and *η*_{2}. We then obtain the desired matrix elements from

To evaluate $\u27e8\kappa \u2032,\eta 1\eta 2\eta 3\u2032|T^cor|\kappa ,\eta 1\eta 2\eta 3\u27e9$, we first compute

for a given *κ* and *η* ≡ *η*_{1}, *η*_{2}, *η*_{3}. We then contract $T^cor|\kappa ,\eta \u2009$ with the bras ⟨*κ*′, *η*_{1}*η*_{2}$\eta 3\u2032$| to obtain $\u27e8\kappa \u2032,\eta 1,\eta 2,\eta 3\u2032|T^cor|\kappa ,\eta 1,\eta 2,\eta 3\u27e9$ for all values of *κ*′ and $\eta 3\u2032$. We repeat the calculation of Eq. (30) and the subsequent contractions for all the other values of *κ* to obtain $\u27e8\kappa \u2032,\eta 1,\eta 2,\eta 3\u2032|T^cor|\kappa ,\eta 1,\eta 2,\eta 3\u27e9$ for all values of *κ*′, *κ*, and $\eta 3\u2032$ and for the initial value of *η*. Finally, we repeat this whole process for every value of *η*.

With all the matrix elements of Δ*Ĥ* computed, we construct the *Ĥ* matrix per Eq. (14) for a given symmetry block. The matrices are small enough that they can be easily diagonalized by direct diagonalization. The convergence of the calculated levels was tested as follows. For *N*_{intra}, it was tested for the *A*_{g} block only, by reducing it from 34 to 21. This seemed to be sufficient since all the symmetry blocks involve the same set of intramolecular states |*γ*⟩. To test the effect of changing *N*_{inter}, we looked at all symmetry blocks except for the *A*_{u} block. For some of these, we roughly doubled *N*_{inter} from the values that pertain to the bases listed in Table IV. Specifically, we changed *N*_{inter} to 13 from 5 for *A*_{g} block, to 24 from 10 for *T*_{1g}, to 22 from 12 for *T*_{2g}, to 26 from 15 for *T*_{1u}, and to 30 from 17 for *T*_{2u}. For the G and H blocks, we roughly halved the *N*_{inter} values pertaining to Table IV: to 13 from 24 for *H*_{g}, to 11 from 19 for *G*_{g}, to 12 from 21 for *H*_{u}, and to 11 from 20 for *G*_{u}. In all instances, this changed the results by only a few hundredths of cm^{−1} for the low-energy states in each of the intramolecular vibrational manifolds.

## III. RESULTS AND DISCUSSION

Tables V–IX present the main results of this study. Those given in Table V pertain to the central topics of the investigation, H_{2}O intramolecular vibrational excitations, and their frequency shifts due to the confinement inside C_{60}. The second and third columns of Table V list the results of the 3D quantum calculations, using *Ĥ*_{intra} in Eq. (9), for the four lowest-energy intramolecular vibrational states of H_{2}O and for the 3D PES of the isolated H_{2}O^{36} (*V*_{intra}) and the C_{60}-adapted 3D H_{2}O PES (*V*_{adap}) defined in Eq. (10), respectively. The energies computed in this work for the case of isolated H_{2}O agree to within a few hundredths to tenths of cm^{−1} with the values reported in Table 1 in the work of Mizus *et al.*^{36} Furthermore, one notes that the three H_{2}O intramolecular vibrational fundamentals and the 2*ν*_{2} overtone calculated on the C_{60}-adapted 3D PES are all blueshifted (Δ*ν*^{3D}) from the ones for the isolated H_{2}O, the two OH stretching fundamentals (*ν*_{1} and *ν*_{3}) significantly so, by ≈21 cm^{−1}, much more than the H_{2}O bend fundamental (*ν*_{2}) and overtone (2*ν*_{2}). The latter are blueshifted by only 1.4 and 3.1 cm^{−1}, respectively.

PES . | V_{intra} (3D)
. | V_{adap} (3D)
. | Δν^{3D}
. | V_{tot} (9D)
. | Δν^{9D}
. |
---|---|---|---|---|---|

ν_{2} | 1594.74 | 1596.14 | 1.40 | 1596.93 | 2.19 |

2ν_{2} | 3151.62 | 3154.72 | 3.10 | 3155.00 | 3.38 |

ν_{1} | 3656.97 | 3678.09 | 21.12 | 3681.14 | 24.17 |

ν_{3} | 3755.84 | 3776.91 | 21.07 | 3779.65 | 23.81 |

PES . | V_{intra} (3D)
. | V_{adap} (3D)
. | Δν^{3D}
. | V_{tot} (9D)
. | Δν^{9D}
. |
---|---|---|---|---|---|

ν_{2} | 1594.74 | 1596.14 | 1.40 | 1596.93 | 2.19 |

2ν_{2} | 3151.62 | 3154.72 | 3.10 | 3155.00 | 3.38 |

ν_{1} | 3656.97 | 3678.09 | 21.12 | 3681.14 | 24.17 |

ν_{3} | 3755.84 | 3776.91 | 21.07 | 3779.65 | 23.81 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.45 | $T1u(a)$ | 1_{01} | |

2 | 36.00 | $T1u(s)$ | 1_{11} | |

3 | 41.23 | $T1g(a)$ | 1_{10} | |

4 | 68.97 | $Hg(s)$ | 2_{02} | |

5 | 77.61 | $Hg(a)$ | 2_{12} | |

6 | 93.32 | $Hu(s)$ | 2_{11} | |

7 | 130.34 | $Hu(a)$ | 2_{21} | |

8 | 131.62 | $Hg(s)$ | 2_{20} | |

14 | 164.37 | $T1u(s)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.45 | $T1u(a)$ | 1_{01} | |

2 | 36.00 | $T1u(s)$ | 1_{11} | |

3 | 41.23 | $T1g(a)$ | 1_{10} | |

4 | 68.97 | $Hg(s)$ | 2_{02} | |

5 | 77.61 | $Hg(a)$ | 2_{12} | |

6 | 93.32 | $Hu(s)$ | 2_{11} | |

7 | 130.34 | $Hu(a)$ | 2_{21} | |

8 | 131.62 | $Hg(s)$ | 2_{20} | |

14 | 164.37 | $T1u(s)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.48 | $T1u(s)$ | 1_{01} | |

2 | 35.47 | $T1u(a)$ | 1_{11} | |

3 | 40.81 | $T1g(s)$ | 1_{10} | |

4 | 68.82 | $Hg(a)$ | 2_{02} | |

5 | 76.99 | $Hg(s)$ | 2_{12} | |

6 | 92.74 | $Hu(a)$ | 2_{11} | |

7 | 128.26 | $Hu(s)$ | 2_{21} | |

8 | 129.86 | $Hg(a)$ | 2_{20} | |

14 | 164.19 | $T1u(a)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.48 | $T1u(s)$ | 1_{01} | |

2 | 35.47 | $T1u(a)$ | 1_{11} | |

3 | 40.81 | $T1g(s)$ | 1_{10} | |

4 | 68.82 | $Hg(a)$ | 2_{02} | |

5 | 76.99 | $Hg(s)$ | 2_{12} | |

6 | 92.74 | $Hu(a)$ | 2_{11} | |

7 | 128.26 | $Hu(s)$ | 2_{21} | |

8 | 129.86 | $Hg(a)$ | 2_{20} | |

14 | 164.19 | $T1u(a)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.93 | $T1u(a)$ | 1_{01} | |

2 | 39.86 | $T1u(s)$ | 1_{11} | |

3 | 45.38 | $T1g(a)$ | 1_{10} | |

4 | 70.54 | $Hg(s)$ | 2_{02} | |

5 | 82.13 | $Hg(a)$ | 2_{12} | |

6 | 98.68 | $Hu(s)$ | 2_{11} | |

7 | 144.91 | $Hu(a)$ | 2_{21} | |

8 | 146.00 | $Hg(s)$ | 2_{20} | |

14 | 163.22 | $T1u(s)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 23.93 | $T1u(a)$ | 1_{01} | |

2 | 39.86 | $T1u(s)$ | 1_{11} | |

3 | 45.38 | $T1g(a)$ | 1_{10} | |

4 | 70.54 | $Hg(s)$ | 2_{02} | |

5 | 82.13 | $Hg(a)$ | 2_{12} | |

6 | 98.68 | $Hu(s)$ | 2_{11} | |

7 | 144.91 | $Hu(a)$ | 2_{21} | |

8 | 146.00 | $Hg(s)$ | 2_{20} | |

14 | 163.22 | $T1u(s)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 25.02 | $T1u(a)$ | 1_{01} | |

2 | 45.13 | $T1u(s)$ | 1_{11} | |

3 | 50.84 | $T1g(a)$ | 1_{10} | |

4 | 71.96 | $Hg(s)$ | 2_{02} | |

5 | 87.40 | $Hg(a)$ | 2_{12} | |

6 | 104.52 | $Hu(s)$ | 2_{11} | |

7 | 162.05 | $Hu(a)$ | 2_{21} | |

8 | 162.86 | $Hg(s)$ | 2_{20} | |

14 | 164.40 | $T1u(s)$ | n = 1 |

. | para
. | ortho
. | . | . |
---|---|---|---|---|

N
. | ΔE
. | ΔE
. | Irrep . | Assignment . |

1 | 25.02 | $T1u(a)$ | 1_{01} | |

2 | 45.13 | $T1u(s)$ | 1_{11} | |

3 | 50.84 | $T1g(a)$ | 1_{10} | |

4 | 71.96 | $Hg(s)$ | 2_{02} | |

5 | 87.40 | $Hg(a)$ | 2_{12} | |

6 | 104.52 | $Hu(s)$ | 2_{11} | |

7 | 162.05 | $Hu(a)$ | 2_{21} | |

8 | 162.86 | $Hg(s)$ | 2_{20} | |

14 | 164.40 | $T1u(s)$ | n = 1 |

The last two columns of Table V show the results of the 9D quantum calculations and the energies of the four lowest-energy intramolecular vibrational states of H_{2}O and their frequency shifts (Δ*ν*^{9D}), respectively. The 9D frequency shifts are somewhat larger than their 3D counterparts, by 14.4% for *ν*_{1}, 13.0% for *ν*_{3}, and 56% for *ν*_{2}. Despite this large increase for the bend mode, its blueshift remains much smaller than those for the two stretching modes. The main difference between the 3D quantum calculations on the C_{60}-adapted H_{2}O PES and the 9D quantum calculations is that the former are performed for the fixed position and orientation of H_{2}O inside C_{60} and considering explicitly only its three intramolecular DOFs, while the latter treat both the intramolecular and intermolecular DOFs of H_{2}O in full dimensionality and as fully coupled. From the comparison of the 3D and 9D frequency shifts, one can conclude that (a) C_{60}-adapted 3D H_{2}O PES includes the bulk of the effects of the interaction between H_{2}O and the C_{60} interior and (b) averaging rigorously over the TR motions and their couplings to the intramolecular DOFs of H_{2}O, performed in the 9D quantum calculations, has a modest effect on the magnitude of the frequency shifts, at least for the PES employed, increasing them by about 14% for the stretch modes relative to those computed for H_{2}O held fixed. In the hindsight, result (b) is not surprising. C_{60} provides an environment with weak angular anisotropy, and tight confinement, for the guest H_{2}O so that the effects of averaging over the quantized TR motions should not be large.

Tables VI–VIII show select low-energy TR states in the *ν*_{1}, *ν*_{3}, and *ν*_{2} intramolecular manifolds of H_{2}O, respectively, from 9D quantum calculations. Table IX shows the same for the 2*ν*_{2} manifold. These results can be compared with those computed in 9D for the ground intramolecular vibrational state of H_{2}O, listed in the fourth column (Δ*E*_{9D}) of Table III, in order to discern how excitation of different intramolecular modes affects the TR levels of H_{2}O@C_{60}. The picture that emerges from the inspection of these tables is interesting. First, the c.m. translational fundamental *n* = 1 at ≈164 cm^{−1} changes very little upon the excitation of any of the intramolecular modes. On the other hand, while the H_{2}O rotational states in the *ν*_{1} and *ν*_{3} intramolecular stretching manifolds have energies that are only slightly lower (at most 2–4 cm^{−1}) than their counterparts in the intramolecular ground-state manifold, the energies of the rotational states of H_{2}O in the intramolecular bend fundamental (*ν*_{2}) and overtone (2*ν*_{2}) manifolds are appreciably higher than those in ground-state manifold. These rotational results can be understood by considering the way in which the rotational constants of the H_{2}O moiety change with intramolecular vibrational state.^{57} Those of the monomer *ν*_{1} and *ν*_{3} states are slightly smaller (∼3% or less) than those of the vibrational ground state. In contrast, the *A* rotational constant of the monomer *ν*_{2} state is about 3.3 cm^{−1} (∼12%) larger, and that of the 2*ν*_{2} state is about 6.6 cm^{−1} (∼24%) larger, than that of the ground state. The vibrational-state changes in rotational constants, together with the assumption of essentially free rotation of the H_{2}O moiety inside C_{60}, can account nearly quantitatively for the intramolecular-state variations in rotational intervals reported in Tables VI–IX.

In the absence of relevant spectroscopic data, it is difficult to assess the accuracy of the theoretical results above, especially the computed frequency shifts. Apart from treating C_{60} as rigid, the 9D quantum methodology for computing fully coupled intramolecular and intermolecular excitations of the entrapped H_{2}O introduced in this paper is rigorous, much more so than the earlier treatments of this problem.^{32–34} The calculated energy levels are highly converged and thus numerically exact for the PES employed. Therefore, the only significant source of uncertainty regarding the results is the quality of this PES. All the calculations to date, including ours in this paper, predict that the intramolecular vibrational modes of H_{2}O are blueshifted due to the confinement in C_{60}. However, whether the intramolecular excitations are blue- or redshifted depends on a subtle balance between the attractive and repulsive interactions of the guest molecule with the confining C_{60} cage. It is not clear whether the PESs utilized so far capture correctly all the competing effects, including the dipole-induced dipole attractive guest-host interaction, since none of them was calculated at a high level of *ab initio* electronic structure theory.

## IV. CONCLUSIONS

We have presented an efficient method for computing accurate intramolecular vibrational fundamentals and overtones of H_{2}O inside the fullerene C_{60}. This 9D quantum treatment assumes rigid C_{60} but is otherwise fully coupled. The calculations also yield low-lying intermolecular TR states within each intramolecular manifold. Symmetries present in this system are exploited, reducing the Hamiltonian matrix to a block-diagonal form, of which ten blocks need to be diagonalized separately in order to obtain the 9D energy levels corresponding to all of the *I*_{h} irreps. Following the recently introduced approach,^{38} the full 9D vibrational Hamiltonian of H_{2}O@C_{60} is divided up into a 6D intermolecular Hamiltonian, a 3D intramolecular Hamiltonian, and a 9D remainder term. The eigenvectors of the two reduced-dimension Hamiltonians are used to construct 9D product contracted (symmetry-adapted) bases in which the symmetry blocks of the full vibrational Hamiltonian are diagonalized. The key to the efficiency of this method is the finding of our recent studies^{37–39} that converged intramolecular vibrational excitations in weakly bound complexes can be obtained by means of fully coupled quantum calculations in which the (6D) intermolecular vibrational eigenstates included in the full-dimensional (9D) basis span the range of energies far below those of the intramolecular vibrational excitations of interest. In the present work, the 6D intermolecular (symmetry adapted) contracted bases include only a small number of eigenstates, at most 24, with energies up to 410 cm^{−1}, which is much less than the H_{2}O stretch and bend fundamental frequencies, at ≈3700 and ≈1600 cm^{−1}, respectively. This results in very small matrices representing the symmetry blocks of the 9D Hamiltonian that need to be diagonalized, the largest with the dimension of 816.

The four lowest-energy intramolecular vibrational levels of H_{2}O in C_{60}, corresponding to the fundamental (*ν*_{2}) and the overtone (2*ν*_{2}) of the water bend, as well as the fundamentals of the *ν*_{1} and *ν*_{3} stretching modes are calculated in two ways: (a) 3D quantum calculation for a C_{60}-adapted 3D PES of H_{2}O, in which the position and orientation of the molecule in C_{60} are fixed, and (b) 9D quantum calculations, where all intramolecular and intermolecular DOFs of H_{2}O are treated as fully coupled. According to both the 3D and the 9D quantum calculations, all four intramolecular vibrational levels considered are blueshifted relative to the ones for the gas-phase H_{2}O. The fundamentals of the two stretching modes, *ν*_{1} and *ν*_{3}, are blueshifted significantly, by ≈24 cm^{−1} in the 9D calculations. The 9D blueshifts of the bend fundamental (*ν*_{2}) and overtone (2*ν*_{2}) are much smaller, 2.19 and 3.38 cm^{−1}, respectively. The frequency shifts from the 9D quantum treatment are somewhat larger than those from the 3D quantum calculations by 14.4% for *ν*_{1}, 13.0% for *ν*_{3}, and 56% for *ν*_{2}. From this, one can conclude that rigorous averaging over the quantum TR motions in the 9D calculations contributes moderately to the frequency shifts for this complex. This is due to the exceptionally high symmetry of C_{60} and hence weak angular anisotropy of the interaction with the guest molecule, as well as its tight confinement. Vibrational averaging is expected to play a more important role in weakly bound molecular complexes with lower symmetry, such as the benzene-H_{2}O dimer.

Finally, the 9D quantum calculations show that the low-energy rotational states of H_{2}O in the *ν*_{1} and *ν*_{3} stretching manifolds differ relatively little, by no more than 4 cm^{−1}, from the corresponding ones in the intramolecular ground-state manifold. However, excitation of the *ν*_{2} bend mode causes larger changes in the energies of the rotational states considered, raising them by up to 12 cm^{−1} for *ν*_{2} manifold and up to 29 cm^{−1} in the 2*ν*_{2} manifold, relative to those in the ground-state manifold. This behavior is consistent with the known vibrational-state dependences of the rotational constants of the H_{2}O moiety. In contrast, the c.m. translational fundamental *n* = 1 is virtually insensitive to excitation of the intramolecular modes, stretch and bend.

Given the qualitative agreement between the 3D and 9D quantum calculations in this study and the earlier and less rigorous treatments^{32–34} that the confinement in C_{60} results in blueshifts of the intramolecular vibrational modes of H_{2}O, one might be tempted to view this question as resolved (apart from the exact values of the shifts). However, this is not necessarily the case. As discussed earlier in the paper, vibrational frequency shifts of confined molecules (inside fullerenes, clathrate hydrates,^{37,39} and other nanocavities) are rather subtle quantities, and it is not obvious that the PESs used thus far are accurate enough to produce reliable values. What is missing are a 9D PES for H_{2}O in (rigid) C_{60} from high-level *ab initio* electronic structure calculations and also spectroscopic measurements of the frequency shifts for comparison with theory. We hope that the existence of the methodology for fully coupled quantum 9D calculations of intramolecular vibrational excitations in H_{2}O@C_{60}, introduced in this work, will stimulate advances in both directions.

## SUPPLEMENTARY MATERIAL

See sec. S1 of the supplementary material for the evaluation of the kinetic-energy matrix elements in the primitive intermolecular basis. Section S2 shows that fivefold truncation of the *α*, *ϕ* 2D grid is possible due to the fivefold rotational symmetry about *Ẑ*. Section S3 describes in detail the operation with *V*_{inter} on a state function. Section S4 demonstrates that the matrix elements of *V*_{inter} in the |*n*, *l*, *m*_{l}, *j*, *m*_{j}, *k*⟩ basis are all real. Section S5 presents the 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 Grant No. CHE-1566085. P.F. is grateful to Professor Daniel Neuhauser for the generous sharing of his computer resources.

## REFERENCES

A. J. Horsewill has informed us that the observed 1_{01} splitting in H_{2}O@C_{60} is such that the lower energy level is nondegenerate and the upper one is doubly degenerate. Due to a misprint, the opposite ordering was erroneously reported in the work of Goh *et al*.^{21}