We present a method for efficient calculation of intramolecular vibrational excitations of H2O inside C60, together with the low-energy intermolecular translation-rotation states within each intramolecular vibrational manifold. Apart from assuming rigid C60, 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 H2O@C60 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 H2O 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 H2O are blueshifted relative to those of the gas-phase H2O, the two stretch modes much more so than the bend. Excitation of the bend mode affects the energies of the low-lying H2O 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 H2O@C60 and ab initio calculation of a high-quality 9D potential energy surface for this endohedral complex, neither of which is presently available.

H2O@C60 is a remarkable supramolecular complex where a polar molecule, H2O, is encapsulated inside a highly symmetric, nonpolar all-carbon interior of C60. 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. H2O@C60 is a member of the family of light-molecule endofullerenes (LMEFs), which have molecules such as H2, HF, and CH4, characterized by small masses and large rotational constants, encapsulated inside the cages of C60, C70, and other fullerenes.1,2 These inclusion compounds have all been synthesized3–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 H2O@C60 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 H2O, and also H2, 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 1H 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 H2O (and H2), having total nuclear spins I = 0 and 1, respectively. The rotational states of H2O are conventionally labeled with the asymmetric top quantum numbers jkakc; for para-H2O, ka + kc has even parity, while for ortho-H2O, ka + kc has odd parity.11 The restrictions placed by the Pauli principle on the rotational quantum numbers of the spin isomers of H2O make its already sparse TR level structure inside the fullerenes even sparser and push the TR dynamics deeper in the quantum regime.

H2O@C60, where the fullerene cage encloses a dipolar H2O molecule (just as HF in C605,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 C60 cages and the possible emergence of dipole-ordered phases and ferroelectricity at low temperatures.14 Various aspects of the quantized rotational dynamics of H2O@C60 have been probed by NMR spectroscopy.15–20 The most comprehensive information about the TR eigenstates of H2O@C60 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-H2O 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 H2O@C60. We calculated the TR eigenstates of para- and ortho-H2O in an (isolated) C60 cage with Ih symmetry by means of fully coupled quantum 6D calculations, treating both H2O and C60 as rigid.22 The computed TR eigenstates allowed tentative assignments of a number of transitions observed in the experimental INS spectra of H2O@C6021 that were not assigned previously. What also emerged from these calculations22 was that the TR level structure of H2O@C60 exhibits all of the key qualitative features that have been identified earlier for H2@C60.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 C60, 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 H2O 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 H2O and the rotational angular momentum j of H2O couple to give the total angular momentum λ = l + j, with λ = l + j, l + j − 1, …, |lj|. 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 H2O@C60, for flexible H2O, but limited to its ground intramolecular vibrational state, and rigid, isolated C60 with Ih 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 H2O and C60 were taken to be rigid.22 

The INS spectra have also revealed that for ortho-H2O inside C60, in solid C60 and at low temperatures, its 101 ground state, which would be threefold degenerate if the cage had the Ih symmetry of isolated C60, 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 H2O must be lower than Ih. This lowering of symmetry, commonly referred to as symmetry breaking, was observed earlier at low temperatures in solid H2@C6027,28 and subsequently also for solid HF@C60.5 

Since our calculations above, and those of Prosmiti et al.,25 were performed for H2O inside an environment having Ih 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 101 ground state of ortho-H2O inside C60 predicted by the quadrupolar model for the dominant P orientation29 are in excellent agreement with experimental results obtained at low temperatures for this system in solid C60.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-H2 and the j = 1 rotational level of HF in C60.29,30

An important and readily observable consequence of the encapsulation of H2O in C60 for which, conspicuously, no experimental results have been reported so far is the shifts in the frequencies of the H2O 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 H231 and HF5 in C60 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 H2O in C60 and their frequency shifts have been reported. In one, harmonic stretching frequencies of the caged H2O were computed at the Hartree-Fock (HF) level, for the static, optimized structure of H2O.32 They gave blueshifts of 41 and 30 cm−1 for the symmetric and asymmetric stretch modes of H2O, respectively. In an attempt to go beyond the static treatment, the IR spectrum of H2O@C60 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 H2O inside the fullerene. Then, at each configuration, the 3D vibrational levels of H2O 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 H2O inside C60 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 H2O stretching modes predicted by the three theoretical treatments32–34 are somewhat surprising, in view of the substantial redshifts observed for both H2 and HF in C60.5,31 It is conceivable that the larger size of H2O, relative to H2 and HF, can result in increased repulsive interaction with the C60 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 H2O in C6032 also gave a slight blueshift for HF in C60, in disagreement with the substantial redshift (−170.5 cm−1) measured for this endohedral complex.5 Therefore, whether the vibrational modes of H2O in C60 are blue- or redshifted is a question that cannot be considered fully settled yet.

As emphasized earlier in this section, all DOFs of H2O inside C60, 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 H2O, and their frequency shifts, coupled to the intermolecular TR motions, requires fully coupled 9D quantum bound-state calculations for flexible H2O and rigid C60 extending to the energies of the H2O vibrational fundamentals and overtones. The earlier quantum treatments of the TR dynamics of H2O@C6022,25,29,30,35 did not, and could not, accomplish this either because H2O was treated as rigid22,29,30,35 or because the calculations were restricted to H2O in the ground intramolecular vibrational state.25 

The frequencies of the intramolecular vibrational fundamentals of H2O [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 phase36] are at least an order of magnitude higher than those calculated by us22 for the intermolecular modes, c.m. translation (162.1 cm−1), and lowest rotational levels of ortho- and para-H2O in C60, 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 (H2O@C60 can be viewed as a special case, where one of the monomers, C60, encloses the other, H2O). 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 H2O 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 H2, HD, and D2 inside the clathrate hydrate cage37 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 H2 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 H2 intramolecular fundamental. We explained this unexpected finding in terms of extremely weak coupling between the intramolecular stretch of the guest H2 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 H2 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 degrees37 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 H2O2, CH4, 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 H2O molecule in (rigid) C60, 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 H2O 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 H2O stretch and bend fundamentals. No ab initio-calculated 9D potential energy surface (PES) for H2O inside C60 is available (not even in 6D, for rigid H2O). Therefore, the 9D PES employed in these calculations is obtained by combining the high-accuracy 3D PES for the isolated H2O molecule36 with the intermolecular PES represented by pairwise-additive atom-atom interactions between H2O and C60 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 H2O in (rigid) C60.22 

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.

We consider a flexible H2O molecule encaged in an isolated and rigid C60 with Ih symmetry. It would be straightforward, at least formally, to include the effects of symmetry breaking29,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 H2O 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 H2O@C60 can be written as

Ĥ(R,ω,q)=R22M+T^vr(ω,q)+T^cor(ω,q)+T^v(q)+Vtot(R,ω,q),
(1)

where M is the mass of H2O. 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^,Ŷ,Z^) with origin at the center of the C60 cage. We choose these axes such that is along a C5 symmetry axis of the icosahedral cage and Ŷ is along one of the C2 symmetry axes, X^=Ŷ×Z^. For the components of R, we work with the spherical coordinates (R, β, α), with R ≡ |R|, cosβR^Z^, and (cosα,sinα)(R^X^,R^Ŷ)/sinβ. 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 (R1, R2, Θ).44–46 We also use Radau bisector-z embedding47 to define the water MF frame, with the MF ŷ axis taken to be parallel to the cross product of Radau vectors R1 ×R2 (i.e., ŷ is perpendicular to the molecular plane) and x^=ŷ×z^.

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 H2O inside C60 (Vtot). This potential term can be written as

Vtot(R,ω,q)=VH2O-cage(R,ω,q)+Vintra(q),
(2)

where the first term represents the interaction between the water moiety and the interior of C60, 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 H2O@C60 was constructed in the same manner.

The explicit forms of T^vr, T^cor, and T^v are given by47 

T^vr(ω,q)=12B1(R1)+B2(R2)ĵx21+cosΘ+ĵy22+ĵz21cosΘ+B1(R1)B2(R2)ĵxĵz+ĵzĵx2sinΘ,
(3)
T^cor(ω,q)B1(R1)B2(R2)p^ΘHĵy,
(4)

and

T^v(q)=i=1212μi2Ri2B1(R1)+B2(R2)×1sinΘΘsinΘΘ,
(5)

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

pΘH=iΘ+12cotΘ
(6)

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

As to the 9D PES Vtot 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 H2O with each atom of C60, with the parameterization (MD1) by Aluru et al.,34 and the C60 nuclear coordinates assumed in our previous work.22 For Vintra, we use the high-quality 3D PES of the isolated H2O (and the computer code to generate it) from the work of Mizus et al.36 

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κinter, of Ĥinter and those, |γ⟩ and Eγ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

Ĥinter(R,ω;q0)R22M+T^vr(ω;q0)+Vinter(R,ω),
(7)

where q0 ≡ (R1,0, R2,0, Θ0) are constants set to values of the water vibrational coordinates close to those of the equilibrium geometry of the monomer, and Vinter is defined as

Vinter(R,ω)=VH2O-cage(R,ω,q0).
(8)

We also define

Ĥintra(q)T^v(q)+Vadap(q),
(9)

where the “C60-adapted” 3D intramolecular potential, Vadap, is defined as the proton-exchange-symmetrized version of

Vadap(q)Vintra(q)+VH2O-cage(R0,ω0,q).
(10)

Here, R0 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 (R1,0, R2,0, Θ0), one obtains a 3D H2O vibrational potential that includes some of the effects (perhaps the bulk thereof) of the H2O–C60 interaction on the vibrational PES of H2O.

With Eq. (1) and the definitions of Eqs. (7) and (9), one has

ΔĤ(R,ω,q)=ΔT^vr(ω,q)+ΔV(R,ω,q)+T^cor(ω,q),
(11)

where

ΔT^vr(ω,q)=T^vr(ω,q)T^vr(ω;q0)
(12)

and

ΔV(R,ω,q)VH2O-cage(R,ω,q)VH2O-cage(R,ω;q0)+Vintra(q)Vadap(q).
(13)

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

κ,γ|Ĥ|κ,γ=(Eκinter+Eγintra)δκ,κδγγ+κ,γ|ΔĤ|κ,γ.
(14)

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

1. Primitive intermolecular basis

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

|n,l,ml,j,mj,k=|n,l,ml|j,mj,k,
(15)

where |n, l, ml⟩ are normalized eigenfunctions of the isotropic 3D harmonic oscillator Hamiltonian48,49 and |j, mj, k⟩ are symmetric-top rotational eigenfunctions.50 The former are of the form

|n,l,ml=Fn,l(MΩ)(R)Ylml(β,α)
(16)

(l = 0, 1, 2 … ; n = l, l + 2, …, nmax; |ml| = 0, 1, …, l), where the radial function, Fn,l(MΩ)(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,

|j,mj,k=2j+18π2[Dmjkj(ϕ,θ,χ)]*.
(17)

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 p2π/5 radians (p an integer) about transforms αα + p2π/5 and ϕϕ + p2π/5 while leaving Ĥinter (and Ĥ) unchanged. Thus, only those |n, l, ml, j, mj, 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 (ml + mj) mod 5. Hence, we factorize the intermolecular basis into five different sets characterized by values of (ml + mj) 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 p2π/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

Ê*|n,l,ml,j,mj,k=(1)ml+j+k|n,l,ml,j,mj,k.
(18)

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 (ml + mj) and by parity-filtering does not in a given calculation produce intermolecular eigenstates that all transform according to the same irreducible representation (irrep) of Ih. This is because these symmetry operations are elements of a subgroup (D5d) of Ih, and, in general, an irrep of Ih correlates with several different irreps of D5d. Thus, for example, an even-parity state with (ml + mj) mod 5 = 0 belongs to either the A1g or A2g irrep of D5d. However, an A1g state in D5d could be either an Ag or an Hg state in Ih. An A2g state in D5d could be either a T1g or a T2g state in Ih. The point is that the calculations employing even(odd)-parity, (ml + mj) mod 5 = 0 basis states yield eigenstates that belong to the Ag(Au), T1g(T1u), T2g(T2u), and Hg(Hu) irreps of Ih. Similarly, calculations employing even(odd)-parity, (ml + mj) mod 5 = ±1 basis states yield eigenstates that belong to the Ih irreps T1g(T1u), Gg(Gu), and Hg(Hu). Finally, even/odd-parity basis states that have (ml + mj) mod 5 = ±2 produce eigenstates belonging to the T2g(T2u), Gg(Gu), and Hg(Hu) irreps of Ih. One sees that it is not necessary to do calculations for all five (ml + mj) symmetry-factored basis-state sets to get the energies and Ih-symmetry assignments of all the levels in a given energy range: It is only necessary to do three such calculations, for example, for (ml + mj) 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 Ihand (b) are composed of basis states that have the same value of (ml + mj) 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 Ih irreps.

3. Computational details

To diagonalize Ĥinter, we make use of the Chebyshev version51 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,

|ψinter=I|II|ψinter,
(19)

where I is a meta-index denoting a particular set of the indices n, l, ml, j, mj, k. Such operations take the form

Ĥinter|ψinter=IĤinter|II|ψ=I|III|Ĥinter|II|ψinter.
(20)

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

I|T^inter|In,l,ml,j,mj,k|R22M+T^vrω;q0×|n,l,ml,j,mj,k=δllδml,mlδjjδmjmjn,l,ml|R22M|n,l,mlδkk+j,mj,k|T^vr(ω;q0)|j,mk,kδnn
(21)

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, Vinter(R, ω; q0), we transform the state function coordinate-by-coordinate to a 6D grid representation (Rp, (cos β)q, (cos θ)r, αs, ϕt, χu). We then multiply it at each grid point by the value of Vinter at that point. Finally, we transform the result back to the |n, l, ml, j, mj, k⟩ basis. The grid we use consists of the following. Rp=up/(MΩ), p = 1, …, Np, where up are Gauss-associated-Laguerre quadrature points. (cos β)q, q = 1, …, Nq, and (cos θ)r, r = 1, …, Nr, 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)/Nt, t = 1, …, Nt, and χu = 2π(u − 1)/Nu, u = 1, …, Nu, is included. However, owing to the C5 symmetry about the axis, one need only include one-fifth of the full range of Fourier points for the α angle: αs = 2π(s − 1)/Ns, s = 1, …, Ns/5, with Ns chosen to be a multiple of 5 (supplementary material). Thus, the full 6D grid comprises 1/5 × NpNqNrNsNtNu points. Full details pertaining to the procedure by which operation with Vinter 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 Vinter, the Ih 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, …, nmax and j = 0, 1, …, jmax were included for each value of (ml + mj) 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), mH = 1836.153 a.u. (the mass of each of the H nuclei), mO = 29 156.377 a.u. (the mass of the oxygen nucleus), and R1,0 = R2,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 formulation44 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.]

TABLE I.

Basis-set and grid parameters associated with the diagonalization of Ĥinter.

nmaxMΩ (a.u.)jmaxNp (#R)Nq ()Nr ()Ns (5 × )Nt ()Nu ()
24.38 12 10 10 20 20 20 
nmaxMΩ (a.u.)jmaxNp (#R)Nq ()Nr ()Ns (5 × )Nt ()Nu ()
24.38 12 10 10 20 20 20 
TABLE II.

Number of symmetry-adapted basis functions for the diagonalization of Ĥinter.

(ml + mj) mod 5 ±1 ±2 
No. of |n, l, ml, j, mj, k⟩ 43 794 43 874 43 954 
(ml + mj) mod 5 ±1 ±2 
No. of |n, l, ml, j, mj, 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 work22 (column three of Table III) on the rigid-monomer intermolecular states of H2O@C60, 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 Ih. 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).

TABLE III.

Selected low-energy intermolecular eigenstates of H2O@C60 from 6D (rigid H2O) and 9D (flexible H2O) calculations. The latter are for H2O in the ground intramolecular vibrational state.

NΔE6D (this work)ΔE6D (Ref. 22)ΔE9D (this work)IrrepAssignment
0.0a 0.0b 0.0c Ag(s) Ground state 
23.94 23.74 23.82 T1u(a) 101 
36.13 36.50 36.73 T1u(s) 111 
41.30 41.85 42.03 T1g(a) 110 
70.47 69.85 70.08 Hg(s) 202 
78.75 78.55 78.99 Hg(a) 212 
94.27 94.61 94.91 Hu(s) 211 
130.04 132.0 132.67 Hu(a) 221 
131.29 133.27 133.94 Hg(s) 220 
10 136.90 135.57 135.99 T2u(a) 303(a
11 137.73 136.38 136.94 Gu(a) 303(b
14 162.11 162.08 163.04 T1u(s) n = 1 
NΔE6D (this work)ΔE6D (Ref. 22)ΔE9D (this work)IrrepAssignment
0.0a 0.0b 0.0c Ag(s) Ground state 
23.94 23.74 23.82 T1u(a) 101 
36.13 36.50 36.73 T1u(s) 111 
41.30 41.85 42.03 T1g(a) 110 
70.47 69.85 70.08 Hg(s) 202 
78.75 78.55 78.99 Hg(a) 212 
94.27 94.61 94.91 Hu(s) 211 
130.04 132.0 132.67 Hu(a) 221 
131.29 133.27 133.94 Hg(s) 220 
10 136.90 135.57 135.99 T2u(a) 303(a
11 137.73 136.38 136.94 Gu(a) 303(b
14 162.11 162.08 163.04 T1u(s) n = 1 
a

E0 = −1986.31 cm−1 relative to separated rigid C60 and H2O.

b

E0 = −1987.57 cm−1 relative to separated rigid C60 and H2O.

c

E0 = −1966.64 cm−1 relative to separated rigid C60 and flexible H2O.

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

|η1,η2,η3|R1,η1|R2,η2|(cosΘ)η3.
(22)

The radial DVRs, |Ri,ηi, are tridiagonal Morse DVRs54 generated by using a Morse potential

V(Ri)=D[ea(RiRe)1]2,
(23)

with the mass factor μi = 1731.932 a.u., D = 0.25 hartree, a = 0.90 bohr−1, and Re = 1.8 bohrs. The bend DVR, |(cosΘ)η3, is a potential-optimized DVR (PODVR)55 generated by first diagonalizing T^v+Vintra for fixed values of R1 = R2 = 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 Ri 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 Vadap = Vintra [i.e., we chose R0 = 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 H2O. Second, we chose R0 = 0 and ω0 = (ϕ0, θ0, χ0) = (0, 0, 0) in computing Vadap from Eq. (10). These values for the R and ω coordinates correspond to those at, or very close to, the minimum of Vinter for the R1,0, R2,0, and Θ0 values that we use. The intramolecular eigenvectors corresponding to this C60-adapted intramolecular potential are those that we ultimately use in the 9D basis.

For both versions of Vadap, 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 Vadap 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,

12μiη1η2η3|2Ri2|η1η2η3=12μiδηj,ηjδη3,η3ηi|2Ri2|ηi,i,j=1,2,ji
(24)

and

η1η2η3|[B1(R1)+B2(R2)]1sinΘΘsinΘΘ|η1η2η3  =δη1η1δη2η2[B1(R1,η1)+B2(R2,η2)]   ×η3|1sinΘΘsinΘΘ|η3.
(25)

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 derived54 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).

With the intermolecular and intramolecular eigenstates computed, we construct the symmetry-adapted 9D bases |κ, γ⟩ by including the 34 lowest-energy intramolecular states, |γ⟩, corresponding to ΔEintra ≤ 11 300 cm−1, and all |κ⟩ with the appropriate symmetry having ΔEinter ≤ 410 cm−1. By “appropriate symmetry,” we refer to the way in which the |κ⟩ transform with respect to (a) the operations of Ih and (b) the specific C5 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 Nintra = 34 intramolecular states |γ⟩, the number of intermolecular states |κ⟩ (Ninter) in the basis is obtained by dividing the basis-set size for each symmetry block in Table IV by 34.

TABLE IV.

Number of symmetry-adapted basis functions for the diagonalization of Ĥ.

(ml + mj), mod 5:(Ih irreps): 0: (Ag, T1g, T2g, Hg+1: Gg 0: (Au, T1u, T2u, Hu+1:(Gu
# of |κ, γ⟩ (170, 340, 408, 816) 646 (34, 510, 578, 714) 680 
(ml + mj), mod 5:(Ih irreps): 0: (Ag, T1g, T2g, Hg+1: Gg 0: (Au, T1u, T2u, Hu+1:(Gu
# 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,

κ,γ|ΔĤ|κ,γ=κ,γ|ΔTvr(ω,q)+ΔV(R,ω,q)+T^cor(ω,q)|κ,γ.
(26)

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 Carrington56) for all κ′, κ, and η and then transform them to the |κ, γ⟩ basis,

κγ|Ô|κ,γ=ηγ|ηη|γκ,η|Ô|κ,η,
(27)

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

Ô|κ,η=I|I,ηII,η|Ô|I,ηI|κ
(28)

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 κ,γ|T^cor(ω,q)|κ,γ 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 κ,η1η2η3|T^cor|κ,η1η2η3 for all κ′, κ, η3, η3, η1, and η2. We then obtain the desired matrix elements from

κ,γ|T^cor(ω,q)|κ,γ=η1,η2η3,η3γ|η1η2η3η1η2η3|γ×κ,η1η2η3|T^cor|κ,η1η2η3.
(29)

To evaluate κ,η1η2η3|T^cor|κ,η1η2η3, we first compute

T^cor|κ,η=I,η3|I,η1η2η3II,η1η2η3|T^cor|I,η1η2η3I|κ
(30)

for a given κ and ηη1, η2, η3. We then contract T^cor|κ,η with the bras ⟨κ′, η1η2η3| to obtain κ,η1,η2,η3|T^cor|κ,η1,η2,η3 for all values of κ′ and η3. We repeat the calculation of Eq. (30) and the subsequent contractions for all the other values of κ to obtain κ,η1,η2,η3|T^cor|κ,η1,η2,η3 for all values of κ′, κ, and η3 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 Nintra, it was tested for the Ag 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 Ninter, we looked at all symmetry blocks except for the Au block. For some of these, we roughly doubled Ninter from the values that pertain to the bases listed in Table IV. Specifically, we changed Ninter to 13 from 5 for Ag block, to 24 from 10 for T1g, to 22 from 12 for T2g, to 26 from 15 for T1u, and to 30 from 17 for T2u. For the G and H blocks, we roughly halved the Ninter values pertaining to Table IV: to 13 from 24 for Hg, to 11 from 19 for Gg, to 12 from 21 for Hu, and to 11 from 20 for Gu. 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.

Tables V–IX present the main results of this study. Those given in Table V pertain to the central topics of the investigation, H2O intramolecular vibrational excitations, and their frequency shifts due to the confinement inside C60. 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 H2O and for the 3D PES of the isolated H2O36 (Vintra) and the C60-adapted 3D H2O PES (Vadap) defined in Eq. (10), respectively. The energies computed in this work for the case of isolated H2O 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 H2O intramolecular vibrational fundamentals and the 2ν2 overtone calculated on the C60-adapted 3D PES are all blueshifted (Δν3D) from the ones for the isolated H2O, the two OH stretching fundamentals (ν1 and ν3) significantly so, by ≈21 cm−1, much more than the H2O bend fundamental (ν2) and overtone (2ν2). The latter are blueshifted by only 1.4 and 3.1 cm−1, respectively.

TABLE V.

Energies (cm−1) of the four lowest-energy H2O intramolecular vibrational states computed for three PESs. Vintra is the 3D PES of the isolated H2O.36Vadap (3D) is defined in Eq. (10). Vtot (9D) is defined in Eq. (2). Also shown are the frequency shifts (cm−1), Δν3D obtained for Vadap, and Δν9D obtained for Vtot. They are relative to the energies computed for Vintra, in the second column. For additional explanation, see the text.

PESVintra (3D)Vadap (3D)Δν3DVtot (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 
PESVintra (3D)Vadap (3D)Δν3DVtot (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 
TABLE VI.

Select low-energy intermolecular TR states of H2O@C60 in the ν1 manifold, from 9D quantum calculations. Their energies ΔE (cm−1) are relative to the ν1 fundamental at 3681.14 cm−1.

paraortho
NΔEΔEIrrepAssignment
 23.45 T1u(a) 101 
36.00  T1u(s) 111 
 41.23 T1g(a) 110 
68.97  Hg(s) 202 
 77.61 Hg(a) 212 
93.32  Hu(s) 211 
 130.34 Hu(a) 221 
131.62  Hg(s) 220 
14 164.37  T1u(s) n = 1 
paraortho
NΔEΔEIrrepAssignment
 23.45 T1u(a) 101 
36.00  T1u(s) 111 
 41.23 T1g(a) 110 
68.97  Hg(s) 202 
 77.61 Hg(a) 212 
93.32  Hu(s) 211 
 130.34 Hu(a) 221 
131.62  Hg(s) 220 
14 164.37  T1u(s) n = 1 
TABLE VII.

Select low-energy intermolecular TR states of H2O@C60 in the ν3 manifold, from 9D quantum calculations. Their energies ΔE (cm−1) are relative to the ν3 fundamental at 3779.65 cm−1.

paraortho
NΔEΔEIrrepAssignment
23.48  T1u(s) 101 
 35.47 T1u(a) 111 
40.81  T1g(s) 110 
 68.82 Hg(a) 202 
76.99  Hg(s) 212 
 92.74 Hu(a) 211 
128.26  Hu(s) 221 
 129.86 Hg(a) 220 
14  164.19 T1u(a) n = 1 
paraortho
NΔEΔEIrrepAssignment
23.48  T1u(s) 101 
 35.47 T1u(a) 111 
40.81  T1g(s) 110 
 68.82 Hg(a) 202 
76.99  Hg(s) 212 
 92.74 Hu(a) 211 
128.26  Hu(s) 221 
 129.86 Hg(a) 220 
14  164.19 T1u(a) n = 1 
TABLE VIII.

Select low-energy intermolecular TR states of H2O@C60 in the ν2 manifold, from 9D quantum calculations. Their energies ΔE (cm−1) are relative to the ν2 fundamental at 1596.93 cm−1.

paraortho
NΔEΔEIrrepAssignment
 23.93 T1u(a) 101 
39.86  T1u(s) 111 
 45.38 T1g(a) 110 
70.54  Hg(s) 202 
 82.13 Hg(a) 212 
98.68  Hu(s) 211 
 144.91 Hu(a) 221 
146.00  Hg(s) 220 
14 163.22  T1u(s) n = 1 
paraortho
NΔEΔEIrrepAssignment
 23.93 T1u(a) 101 
39.86  T1u(s) 111 
 45.38 T1g(a) 110 
70.54  Hg(s) 202 
 82.13 Hg(a) 212 
98.68  Hu(s) 211 
 144.91 Hu(a) 221 
146.00  Hg(s) 220 
14 163.22  T1u(s) n = 1 
TABLE IX.

Select low-energy intermolecular TR states of H2O@C60 in the 2ν2 manifold, from 9D quantum calculations. Their energies ΔE (cm−1) are relative to the 2ν2 overtone at 3155.00 cm−1.

paraortho
NΔEΔEIrrepAssignment
 25.02 T1u(a) 101 
45.13  T1u(s) 111 
 50.84 T1g(a) 110 
71.96  Hg(s) 202 
 87.40 Hg(a) 212 
104.52  Hu(s) 211 
 162.05 Hu(a) 221 
162.86  Hg(s) 220 
14 164.40  T1u(s) n = 1 
paraortho
NΔEΔEIrrepAssignment
 25.02 T1u(a) 101 
45.13  T1u(s) 111 
 50.84 T1g(a) 110 
71.96  Hg(s) 202 
 87.40 Hg(a) 212 
104.52  Hu(s) 211 
 162.05 Hu(a) 221 
162.86  Hg(s) 220 
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 H2O 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 C60-adapted H2O PES and the 9D quantum calculations is that the former are performed for the fixed position and orientation of H2O inside C60 and considering explicitly only its three intramolecular DOFs, while the latter treat both the intramolecular and intermolecular DOFs of H2O in full dimensionality and as fully coupled. From the comparison of the 3D and 9D frequency shifts, one can conclude that (a) C60-adapted 3D H2O PES includes the bulk of the effects of the interaction between H2O and the C60 interior and (b) averaging rigorously over the TR motions and their couplings to the intramolecular DOFs of H2O, 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 H2O held fixed. In the hindsight, result (b) is not surprising. C60 provides an environment with weak angular anisotropy, and tight confinement, for the guest H2O 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 H2O, 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 H2O, listed in the fourth column (ΔE9D) of Table III, in order to discern how excitation of different intramolecular modes affects the TR levels of H2O@C60. 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 H2O 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 H2O 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 H2O 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 H2O moiety inside C60, 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 C60 as rigid, the 9D quantum methodology for computing fully coupled intramolecular and intermolecular excitations of the entrapped H2O 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 H2O are blueshifted due to the confinement in C60. 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 C60 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.

We have presented an efficient method for computing accurate intramolecular vibrational fundamentals and overtones of H2O inside the fullerene C60. This 9D quantum treatment assumes rigid C60 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 Ih irreps. Following the recently introduced approach,38 the full 9D vibrational Hamiltonian of H2O@C60 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 studies37–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 H2O 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 H2O in C60, 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 C60-adapted 3D PES of H2O, in which the position and orientation of the molecule in C60 are fixed, and (b) 9D quantum calculations, where all intramolecular and intermolecular DOFs of H2O 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 H2O. 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 C60 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-H2O dimer.

Finally, the 9D quantum calculations show that the low-energy rotational states of H2O 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 H2O 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 treatments32–34 that the confinement in C60 results in blueshifts of the intramolecular vibrational modes of H2O, 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 H2O in (rigid) C60 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 H2O@C60, introduced in this work, will stimulate advances in both directions.

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 Vinter on a state function. Section S4 demonstrates that the matrix elements of Vinter in the |n, l, ml, j, mj, k⟩ basis are all real. Section S5 presents the 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 Grant No. CHE-1566085. P.F. is grateful to Professor Daniel Neuhauser for the generous sharing of his computer resources.

1.
Z.
Bačić
,
J. Chem. Phys.
149
,
100901
(
2018
).
2.
M. H.
Levitt
,
Philos. Trans. R. Soc., A
371
,
20120429
(
2013
).
3.
K.
Komatsu
,
M.
Murata
, and
Y.
Murata
,
Science
307
,
238
(
2005
).
4.
K.
Kurotobi
and
Y.
Murata
,
Science
333
,
613
(
2011
).
5.
A.
Krachmalnicoff
,
R.
Bounds
,
S.
Mamone
,
S.
Alom
,
M.
Concistrè
,
B.
Meier
,
K.
Kouřil
,
M. E.
Light
,
M. R.
Johnson
,
S.
Rols
,
A. J.
Horsewill
,
A.
Shugai
,
U.
Nagel
,
T.
Rõõm
,
M.
Carravetta
,
M.
Levitt
, and
R. J.
Whitby
,
Nat. Chem.
8
,
953
(
2016
).
6.
S.
Bloodworth
,
G.
Sitinova
,
S.
Alom
,
S.
Vidal
,
G. R.
Bacanu
,
S. J.
Elliott
,
M. E.
Light
,
J. M.
Herniman
,
G. J.
Langley
,
M. H.
Levitt
, and
R. J.
Whitby
,
Angew. Chem., Int. Ed.
58
,
5038
(
2019
).
7.
Y.
Rubin
,
Chem. Eur. J.
3
,
1009
(
1997
).
8.
Y.
Rubin
,
Top. Curr. Chem.
199
,
67
(
1999
).
9.
Y.
Rubin
,
T.
Jarrosson
,
G. W.
Wang
,
M. D.
Bartberger
,
K. N.
Houk
,
G.
Schick
,
M.
Saunders
, and
R. J.
Cross
,
Angew. Chem., Int. Ed.
40
,
1543
(
2001
).
10.
Z.
Bačić
,
M.
Xu
, and
P. M.
Felker
,
Adv. Chem. Phys.
163
,
195
(
2018
).
11.
P. R.
Bunker
and
P.
Jensen
,
Molecular Symmetry and Spectroscopy
, E-book edition (
NRC Research Press
,
Ottawa, Ontario, Canada
,
2006
).
12.
M.
Xu
,
P. M.
Felker
,
S.
Mamone
,
A. J.
Horsewill
,
S.
Rols
,
R. J.
Whitby
, and
Z.
Bačić
,
J. Phys. Chem. Lett.
10
,
5365
(
2019
).
13.
S.
Aoyagi
,
N.
Hoshino
,
T.
Akutagawa
,
Y.
Sado
,
R.
Kitaura
,
H.
Shinohara
,
K.
Sugimoto
,
R.
Zhang
, and
Y.
Murata
,
Chem. Commun.
50
,
524
(
2014
).
14.
J.
Cioslowski
and
A.
Nanayakkara
,
Phys. Rev. Lett.
69
,
2871
(
1992
).
15.
C.
Beduz
,
M.
Carravetta
,
J. Y.-C.
Chen
,
M.
Concistre
,
M.
Denning
,
M.
Frunzi
,
A. J.
Horsewill
,
O. G.
Johannessenn
,
R.
Lawler
,
X.
Lei
,
M. H.
Levitt
,
Y.
Li
,
S.
Mamone
,
Y.
Murata
,
U.
Nagel
,
T.
Nishida
,
J.
Ollivier
,
S.
Rols
,
T.
Rõõm
,
R.
Sarkar
,
N. J.
Turro
, and
Y.
Yang
,
Proc. Natl. Acad. Sci. U. S. A.
109
,
12894
(
2012
).
16.
S.
Mamone
,
M.
Concistre
,
E.
Carignani
,
B.
Meier
,
A.
Krachmalnicoff
,
O.
Johannessen
,
X.
Lei
,
Y.
Li
,
M.
Denning
,
M.
Carravetta
,
K.
Goh
,
A. J.
Horsewill
,
R. J.
Whitby
, and
M. H.
Levitt
,
J. Chem. Phys.
140
,
194306
(
2014
).
17.
B.
Meier
,
S.
Mamone
,
M.
Concistre
,
J.
Alonso-Valdesueiro
,
A.
Krachmalnicoff
,
R. J.
Whitby
, and
M. H.
Levitt
,
Nat. Commun.
6
,
8112
(
2015
).
18.
S. J.
Elliot
,
C.
Bengs
,
K.
Kouřil
,
B.
Meier
,
S.
Alom
,
R. J.
Whitby
, and
M. H.
Levitt
,
ChemPhysChem
19
,
251
(
2018
).
19.
B.
Meier
,
K.
Kouřil
,
C.
Bengs
,
H.
Kouřilová
,
T. C.
Barker
,
S. J.
Elliot
,
S.
Alom
,
R. J.
Whitby
, and
M. H.
Levitt
,
Phys. Rev. Lett.
120
,
266001
(
2018
).
20.
H.
Suzuki
,
M.
Nakano
,
Y.
Hashikawa
, and
Y.
Murata
,
J. Phys. Chem. Lett.
10
,
1306
(
2019
).
21.
K. S. K.
Goh
,
M.
Jiménez-Ruiz
,
M. R.
Johnson
,
S.
Rols
,
J.
Ollivier
,
M. S.
Denning
,
S.
Mamone
,
M. H.
Levitt
,
X.
Lei
,
Y.
Li
,
N. J.
Turro
,
Y.
Murata
, and
A. J.
Horsewill
,
Phys. Chem. Chem. Phys.
16
,
21330
(
2014
).
22.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
144
,
201101
(
2016
).
23.
M.
Xu
,
F.
Sebastianelli
,
Z.
Bačić
,
R.
Lawler
, and
N. J.
Turro
,
J. Chem. Phys.
128
,
011101
(
2008
).
24.
M.
Xu
,
F.
Sebastianelli
,
Z.
Bačić
,
R.
Lawler
, and
N. J.
Turro
,
J. Chem. Phys.
129
,
064313
(
2008
).
25.
Á.
Valdés
,
O.
Carrillo-Bohórquez
, and
A. R.
Prosmiti
,
J. Chem. Theory Comput.
14
,
6521
(
2018
).
26.

A. J. Horsewill has informed us that the observed 101 splitting in H2O@C60 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 

27.
Y.
Kohama
,
T.
Rachi
,
J.
Jing
,
Z.
Li
,
J.
Tang
,
R.
Kumashiro
,
S.
Izumisawa
,
H.
Kawaji
,
T.
Atake
,
H.
Sawa
,
Y.
Murata
,
K.
Komatsu
, and
K.
Tanigaki
,
Phys. Rev. Lett.
103
,
073001
(
2009
).
28.
S.
Mamone
,
M. R.
Johnson
,
J.
Ollivier
,
S.
Rols
,
M. H.
Levitt
, and
A. J.
Horsewill
,
Phys. Chem. Chem. Phys.
18
,
1998
(
2016
).
29.
P. M.
Felker
,
V.
Vlček
,
I.
Hietanen
,
S.
FitzGerald
,
D.
Neuhauser
, and
Z.
Bačić
,
Phys. Chem. Chem. Phys.
19
,
31274
(
2017
).
30.
Z.
Bačić
,
V.
Vlček
,
D.
Neuhauser
, and
P. M.
Felker
,
Faraday Discuss.
212
,
547
(
2018
).
31.
M.
Ge
,
U.
Nagel
,
D.
Hüvonen
,
T.
Rõõm
,
S.
Mamone
,
M. H.
Levitt
,
M.
Carravetta
,
Y.
Murata
,
K.
Komatsu
,
J. Y.-C.
Chen
, and
N. J.
Turro
,
J. Chem. Phys.
134
,
054507
(
2011
).
32.
O.
Shameema
,
C. N.
Ramachandran
, and
N.
Sathyamurthy
,
J. Phys. Chem. A
110
,
2
(
2006
).
33.
K.
Yagi
and
D.
Watanabe
,
Int. J. Quantum Chem.
109
,
2080
(
2009
).
34.
A. B.
Farimani
,
Y.
Wu
, and
N. R.
Aluru
,
Phys. Chem. Chem. Phys.
15
,
17993
(
2013
).
35.
E.
Rashed
and
J. L.
Dunn
,
Phys. Chem. Chem. Phys.
21
,
3347
(
2019
).
36.
I. I.
Mizus
,
A. A.
Kyuberis
,
N. F.
Zubov
,
V. Y.
Makhnev
,
O. L.
Polyansky
, and
J.
Tennyson
,
Philos. Trans. R. Soc., A
376
,
20170149
(
2018
).
37.
D.
Lauvergnat
,
P. M.
Felker
,
Y.
Scribano
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
150
,
154303
(
2019
).
38.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
151
,
024305
(
2019
).
39.
P. M.
Felker
,
D.
Lauvergnat
,
Y.
Scribano
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
151
,
124311
(
2019
).
40.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
117
,
6923
(
2002
).
41.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
119
,
101
(
2003
).
42.
J. C.
Tremblay
and
T.
Carrington
,
J. Chem. Phys.
125
,
094311
(
2006
).
43.
X.-G.
Wang
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
129
,
234102
(
2009
).
44.
B. R.
Johnson
and
W. P.
Reinhardt
,
J. Chem. Phys.
85
,
4538
(
1986
).
45.
Z.
Bačić
,
D.
Watt
, and
J. C.
Light
,
J. Chem. Phys.
89
,
947
(
1988
).
46.
B. T.
Sutcliffe
and
J.
Tennyson
,
Int. J. Quantum Chem.
39
,
183
(
1991
).
47.
X.-G.
Wang
and
T.
Carrington
,
J. Chem. Phys.
146
,
104105
(
2017
).
49.
See, for example, supplementary material, Sec. II, from
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
145
,
084310
(
2016
).
50.
R. N.
Zare
,
Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics
(
Wiley-Intersience
,
New York
,
1988
).
51.
V. A.
Mandelshtam
and
H. S.
Taylor
,
J. Chem. Phys.
106
,
5085
(
1997
).
52.
M. R.
Wall
and
D.
Neuhauser
,
J. Chem. Phys.
102
,
8011
(
1995
).
53.
B.
Poirier
,
J. Chem. Phys.
143
,
101104
(
2015
).
54.
H.
Wei
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
97
,
3029
(
1992
).
55.
J.
Echave
and
D. C.
Clary
,
Chem. Phys. Lett.
190
,
225
(
1992
).
56.
X.-G.
Wang
and
T.
Carrington
,
J. Chem. Phys.
148
,
074108
(
2018
).
57.
G.
Herzberg
,
Molecular Spectra and Molecular Structure. III. Electronic Spectra and Electronic Structure of Polyatomic Molecules
(
Van Nostrand
,
Princeton
,
1966
).

Supplementary Material