Encapsulation of a Water Molecule inside C60 Fullerene: The Impact of Confinement on Quantum Features

We introduce an efficient quantum fully coupled computational scheme within the multiconfiguration time-dependent Hartree (MCTDH) approach to handle the otherwise extremely costly computations of translational–rotational–vibrational states and energies of light-molecule endofullenes. Quantum calculations on energy levels are reported for a water molecule inside C60 fullerene by means of such a systematic approach that includes all nine degrees of freedom of H2O@C60 and does not consider restrictions above them. The potential energy operator is represented as a sum of natural potentials employing the n-mode expansion, along with the exact kinetic energy operator, by introducing a set of Radau internal coordinates for the H2O molecule. On the basis of the present rigorous computations, various aspects of the quantized intermolecular dynamics upon confinement of H2O@C60 are discussed, such as the rotational energy level splitting and the significant frequency shifts of the encapsulated water molecule vibrations. The impact of water encapsulation on quantum features is explored, and insights into the nature of the underlying forces are provided, highlighting the importance of a reliable first-principles description of the guest–host interactions.


INTRODUCTION
Compounds in which atoms or small molecules are trapped in the cavity of fullerenes, such as C 60 or C 70 , are known as endofullerenes or endohedral fullerenes and offer the opportunity to explore unusual patterns of the entrapped species. The guest confinement inside the fullerene cages influences both the host and the guest, leading to various applications ranging from fundamental research to applied technology in medicine, material science, and electronics, in particular single-molecule transistors for quantum computing. 1−10 In view of exciting properties and the prospects of applications, as new energy storage materials, among others in nanoscience, these species were revealed to exhibit some unique and unexpected properties, which triggered their extensive investigation from both theoretical and experimental sides. 2, 11−31 A step advance in endofullerene science has been made with the advent of a novel, multistep organic synthesis procedure known as molecular surgery. 9−12,29,32−34 Such a procedure consists of a series of chemical reactions for creating an open C 60 cage, in which the guest molecule can be inserted, followed by the closing of the cage with the guest trapped inside. Various molecules, such as H 2 , HD, HF, H 2 O, and CH 4 , have been encapsulated into C 60 in this manner to date, and their properties have been experimentally investigated by X-ray diffraction, inelastic neutron scattering (INS), far-infrared spectroscopy (FIR), and nuclear magnetic resonance (NMR) spectroscopy techniques. 14,16−18 Such confined light molecules exhibit dominant quantum effects as a result of strong couplings between translational and rotational motions. Further, H 2 @C 60 , H 2 O@C 60 , and CH 4 @C 60 endohedral fullerenes are of particular interest, as they also exhibit nuclear spin isomerism, introducing new quantum features, and only certain combinations of nuclear spin and rotational states are allowed by the Pauli exclusion principle.
H 2 O@C 60 endofullerene, considered as polar C 60 , enclosing a dipolar water molecule in its cage, is particularly intriguing for a variety of reasons. The co-crystallized structure of H 2 O@ C 60 with (NiOEP) 2 has been determined by single-crystal Xray diffraction analysis, 12 and it has been reported 12 that the oxygen atom of the encapsulated H 2 O molecule is located at the center of the C 60 cage, with the position of the O−H bonds to be directed toward the Ni ions through the C 60 cage. Further, on the basis of UV/vis, Fourier transform infrared (FTIR), and NMR spectra analysis of H 2 O@C 60 and empty C 60 systems, it has been concluded 12 that the inclusion of the H 2 O molecule does not affect the structure of the C 60 cage, indicating no detectable electronic interaction between water and the cage. However, theoretical molecular dynamics (MD) and electronic structure studies 13 38 In addition, a theoretical investigation on the covalent binding of H 2 O to C 60 upon compression of H 2 O@C 60 endofullerene has revealed 15 that the pressure-induced intercavity chemical reaction can take place depending on the direction of the compression, leading to the formation of endohedral covalent fullerene, if applied in the direction of the two opposite pentagons. In addition, earlier INS, FTIR spectroscopy, and cryogenic NMR experiments 14,18,27,30,43 and more recent INS spectra 17 in highly pure samples of solid H 2 O@C 60 have revealed energy splitting in the ground state of the encapsulated ortho-water, raising the 3fold degeneracy into single and doublet states (see refs 44 and 26 therein), associated with symmetry breaking of the water environment.
From a theoretical perspective, the investigation of such unexpected properties of water endofullerene requires consideration of quantum treatments in both electronic and nuclear degrees of freedom, and a variety of theoretical models and scenarios have been elaborated 21,22,25,28,45,46 to date to explore the experimental observations. In view of the complexity of the problem, most of the previous studies 22,44 have highlighted the importance of fully coupled quantum treatments of both nuclear and electronic degrees of freedom of such nine-dimensional (9D) guest−host systems. Thus, our current investigation 22,47−53 aims to contribute to the field by proposing a systematic protocol to crosscheck first-principles methodologies on noncovalent guest−host interactions and an efficient procedure within the multiconfiguration time-dependent Hartree (MCTDH) framework 54,55 to treat the nuclear quantum dynamics of any nanoconfined light molecule. Here, our efforts are focused on fully coupled quantum treatments employing available three-dimensional (3D) and six-dimensional (6D) potentials for the isolated and enclathrated H 2 O molecule, respectively. So, the problem is conveniently discussed in two stages: (a) rotational, translational, and vibrational quantization of a nanoconfined light−heavy−light molecule through a novel computational implementation of an efficient Radau coordinate Hamiltonian representation, allowing to converge the 9D fully coupled computations of highly excited vibrational states of H 2 O@C 60 and (b) the effect of water molecule encapsulation on its quantum spectral properties due to interfullerene interactions.
The paper is organized as follows. In Section 2, we describe our approach and outline the computational details employed. In Section 3, we present benchmark results on rotational and novel translational and vibrational energy values, as well as comparisons, in cases that are available with previously reported values. Final conclusions and future directions are drawn in Section 4.

COMPUTATIONAL DETAILS
In this study, the methodology already developed in ref 22 is exploited for any light−heavy−light encapsulated molecule to describe the quantum mechanical features of H 2 O@C 60 endohedral fullerene. We demonstrate the robustness of the method on systems that cannot be described through simple models by carrying out high-accuracy calculations of the molecular levels of the system.
2.1. Coordinate System and Hamiltonian Operator. In Figure 1, we display the 9D coordinate system used to represent the Hamiltonian operator of the fully coupled system, consisting of a flexible H 2 O molecule in rigid C 60 cages. The origin of the Cartesian xyz system is fixed at the center of the fullerene cage with the z-axis chosen to lie along the C 5 symmetry axis of C 60 , while the y-axis coincides with one of the C 2 symmetry axes. The position of the water molecule's center of mass inside the fullerene cage is described by means of the spherical coordinates (R, β, α) with respect to the xyz system, where R is the center of the mass vector of the water molecule, and β and α are its polar and azimuthal angles, respectively. The BF coordinate system x BF y BF z BF refers to the space-fixed xyz system through the set of Euler (ϕ, θ, χ) angles (cf. Figure 1), with its origin at the center of mass of the water molecule. An appropriate set of coordinates to describe a triatomic light−heavy−light molecule, 56−58 as H 2 O, is the internal Radau coordinates ( 1 , 2 ,γ), where 1 ⎯→ ⎯ and 2 ⎯ → ⎯⎯ are the position vectors of the hydrogen atoms with respect to the Radau canonical point, 56 and γ the angle between them. The z BF axis was chosen parallel to the 1 ⎯→ ⎯ vector, whereas 2 ⎯ → ⎯⎯ lies over the x BF −z BF plane.
The Hamiltonian operator, Ĥ= T(q) + V̂(q), was deduced in the set of q R , , , , , , , , In previous works of the H 2 O@C 60 system, 22,44,60 the potential term, V, was built up employing the sum-of-potentials approach, with the intermolecular interaction between the water molecule and C 60 , plus the intramolecular water p o t e n t i a l , , with radial coordinates in Å, and angular ones in rad. The C 60 cage has 10 C 2 -axes, and we should note that each of the symmetric minimum energy configurations is located along them, at n 5 α = π , with n = 1−10.

MCTDH Computational Setup.
Within the MCTDH scheme, 55 the time-dependent Schrodinger equation is solved by expanding the wave function in a sum of products of the time-dependent low-dimensional basis of so-called single-particle functions (SPFs), ϕ(Q, t), so Ψ(q 1 , ..., where f and p are the number of degrees of freedom and the number of particles or combined modes, respectively, while A J is the MCTDH expansion coefficient, with J being the composite index (j 1 , ..., j p ), and Φ J is a Hartree product of SPFs, ϕ(q, t). The SPFs are in turn represented by linear combinations of time-independent primitive basis functions or discrete variable representation (DVR) grids. The solution of the MCTDH equations of motion requires the computation of the mean-field matrices at every time step. Thus, if the Hamiltonian operator is expressed by a sum of products of monoparticle operators, the so-called product structure, then multidimensional integrals can be solved efficiently, allowing the treatment of high-dimensional problems. The kinetic energy operators are in the sum-ofproducts form, while potential energy operators must be expanded in the product form representation. Within the MCTDH code, 55 the POTFIT approach 62 is the default procedure to transform low-dimensional PESs to the product form, while for treating systems of larger dimensionality, alternatives, such as the n-mode representation, 63 the multigrid POTFIT, 64 Monte Carlo POTFIT (MCPOTFIT), 65 or those reported more recently, such as the Monte Carlo canonical polyadic decomposition (MCCPD) 66 and the rectangular collocation MCTDH (RC-MCTDH) 67 methods are used.
We used an n-mode representation of the potential, 63 as the grid of at least 10 9 points required for the 9D potential is too large to be treated with the POTFIT algorithm. 62 By considering the strong coupling between the 9D coordinates for the H 2 O@C 60 system, we adopted a 7-mode combination scheme as, , with V (0) being the reference configuration potential value, V (1) the intramode terms, while the remaining V (n) are the 2-up to 7-mode correlation terms.
As previously discussed, for affordable 9D calculations, we proceeded to truncate the full 7-mode expansion by considering the following n-mode selective representations of the V 9D potential. Thus, V 9D (Q 1 , ..., takes into account the internal coordinates of the water molecule, as well as the position of its center of mass, with the V 6D (Q 1 ,Q 2 ,Q 3 ,Q 4 ) = V 5D (Q 2 ,Q 3 ,Q 4 ) + Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article ΔV 6D R (Q 1 ,Q 2 ,Q 3 ,Q 4 ) term being the exact potential representation of the 6D system, and ΔV 6D Qi (Q i ,Q 2 ,Q 3 ,Q 4 ) = V 6D (Q i ,Q 2 ,Q 3 ,Q 4 ) − V 5D (Q 2 ,Q 3 ,Q 4 ) with i = 1, 5, 6, 7 being the corresponding terms in the expansion, keeping fixed at their equilibrium values the independent degrees of freedom at each potential term. According to the n-mode representation, the resulting nine-dimensional potential was approximated as a sum of potentials of five and six dimensions, and the POTFIT algorithm, 62 as implemented in the MCTDH code, 55 was applied to each of these potential terms to represent them as a sum of operators of one dimension (natural potentials). Some details on the convergence of such approaches are given in Figure S1 and Table S1 (see the SI). In particular. in Table S1, we list the POTFIT parameters (number of natural potentials and the relevant regions of the potential) used for each of the n-mode terms together with the corresponding root-meansquare (RMS) errors of the potential fits, while in Figure S1 (see the SI), we show the convergence of the V 9D expansion as a function of n-mode expansion terms over a total number of 10 9 configurations in the range of each of the nine internal coordinates, as defined in Table S2 (see the SI). One can see that 90% of the configurations have a mean absolute error (MAE) of less than 2 cm −1 , while for the remaining of them, it does not exceed 4 cm −1 . Once the Hamiltonian operator was set up, we should also define the parameters of the grid for its representation. Taking into account the singularities of the KEO, we have carefully chosen the type and range of the DVR basis sets employed in each coordinate. Thus, the underlying primitive basis sets used in the MCTDH calculations are the harmonic oscillator (HO) and radial form solutions (rHO) for radial coordinates, and Legendre (Leg) and exponential (exp) DVR functions for the angular ones depending on the coordinate range, except for the γ angle where the restricted Legendre-type (Leg/R) DVR functions were employed, that allow choosing specific angular range intervals. In Table S2, we list the grid parameters (range, number, and type of primitive basis set functions) for each of the nine coordinates in the MCTDH calculations for obtaining the rotational/translational/vibrational states of H 2 O@C 60 endofullerene, while Table S3 summarizes the number of SPFs employed in each combination mode in the present computations.

RESULTS AND DISCUSSION
3.1. Benchmarking Rotational State Calculations. For calculating rotationally excited states, we employed the block improved relaxation procedure (BIR) 68 implemented in the Heidelberg MCTDH package. 55 We initiated with a block of 42 initial vectors corresponding to the number of the desired states that are simultaneously converged collectively, within 10 −3 cm −1 , using the same set of SPFs for all of them. This set of 42 eigenstates includes the ground and another 10 (multiple degenerated) lower rotationally excited states. Converging to higher excited states has become more difficult as the excitation energy increases and couplings with translational motions are present. Thus, here, we assigned pure rotational states, and in Table  1, we show their energy values together with their degeneracy, g, and corresponding 2D probability density distributions of one of the degenerated states. Each H 2 O rotational excited state is labeled as J K a K c , according to quantum numbers of total rotational angular momentum J, while K a and K c are related to the angular momentum for prolate and oblate symmetric tops, respectively. J K a K c states have a 2J + 1 degeneracy and can be further classified into para-H 2 O (total nuclear spin I = 0, and even parity for K a + K c ) and ortho-H 2 O (I = 1 and K a + K c of odd parity) states. The ground state of para-H 2 O is singlet 0 00 , while the ground ortho-H 2 O state is triplet 1 01 . We found the ground state energy for H 2 O@C 60 at E H 2 O@C 60 0 = −1967.37 cm −1 , which agrees very well with the value of −1966.64 cm −1 reported recently from quantum calculations 44 and in our previous work 22 using a set of internal valence coordinates instead of the Radau grid. Apart from the ground state, small differences are also observed for the rotationally excited levels. Such subtle deviations can be attributed to the differences between the spatial positioning of the grids used in each set of coordinates in the present and earlier 9D MCTDH calculations 22 or to different molecular parameters and potentials used for the water molecule. 44 Also, Table 1 shows the energy levels obtained for the isolated H 2 O molecule in the present 6D MCTDH calculations in comparison with those reported in the literature 60 for a freely rotating rigid water molecule. One can see that the energies of the lower rotational levels of the encapsulated para-and ortho-H 2 O are close (within about 0.5 cm −1 ) to those of the isolated freely rotating H 2 O, indicating a rather weakly hindered rotation of the trapped molecule. However, as the excitation energy increases, the observed energy differences also increase, with values up to 2.3 cm −1 between those of the isolated and encapsulated water molecules, for the higher rotational states in the present work (see Figure 3). In addition, for J ≥ 3, rotational levels of the isolated water molecule appear with the 2J + 1 degeneracy, while the presence of the icosahedral, I h symmetry, C 60 cage, breaks such degeneracy. For example, in the case of J = 3, the 7-fold-degenerated 3 03 levels of the isolated water molecule at an energy of 136.50 cm −1 are split into a pair of 3-folddegenerated, namely 3 03 (a), and 4-fold-degenerated, namely 3 03 (b), levels at energies of 135.29 and 137.11 cm −1 , respectively, for the encapsulated H 2 O molecule. The energy splitting of the 3 03 levels induced by the nonsphericity of C 60 is 1.8 cm −1 , indicating a rather small deviation. In the same vein, by comparing the corresponding probability density distributions between the isolated 3 03 and trapped 3 03 (a) water molecule rotational states, we identified some tiny differences, mainly in the ϕ/θ coordinates, which could indicate some preferent orientation of the encapsulated water molecule into the cage.
3.2. Novel Translational and Vibrational State Calculations. As earlier discussed, the rotational and translational degrees of freedom of the water molecule inside C 60 are strongly coupled, with a rigorous approach requiring also their full coupling to the vibrational water modes. However, we should emphasize that intermolecular rotational/translational frequencies are much lower than those of intramolecular vibrational ones; thus, a large number (several hundred) of highly excited rotational−translational states are lying between vibrational excited states, which makes computations rather expensive, and in cases, when the excitation energy increases, even prohibitive. Thus, so far, such issues on fully 9D coupled rotational−translational− vibrational excitations have not been widely addressed by theory, and only a few attempts have been reported in the literature 22,44,47,69,70 for nanoconfined systems. Indeed, the computationally most challenging objective was to compute specific vibrational states of the water molecule inside the C 60 cage that allows identifying frequency shifts in H 2 O vibrational modes, indicating further possible observable consequences of its encapsulation.
To perform such calculations of high-lying translational and vibrational states of H 2 O@C 60 , we employed the improved relaxation (IR) method, 71 as implemented in the Heidelberg MCTDH code. 55 The energy relaxation method uses imaginary time propagation of the wave function to generate the ground state of a system. The relaxation can be accelerated and, also, excited states can be computed if the MCTDH coefficient vector is not determined by relaxation but by diagonalization. Thus, we choose an initial state that has a reasonable overlap with the eigenstate of interest. Then, the Hamiltonian is diagonalized by a Davidson routine within the SPFs, followed by the selection of the specific eigenstate, the calculations of the mean fields, the relaxation of the SPFs, and so on, until convergence is achieved. For the calculation of the ground state, the lowest energy eigenvector is selected, while for any individual translationally or vibrationally excited state, the eigenvector that has the largest overlap with the initial state is chosen. Converging to excited states is numerically more demanding and depends mainly on how well separated are the states of interest. Then, we can, in principle, converge to any of them by selecting the appropriate initial wave function.
Both translational and vibrational energy levels are computed from the 9D Hamiltonian. A set of lock calculations (whose algorithm is implemented in the MCTDH package 55 ) were performed along with the IR method. With the lock keyword, a calculation converges to the eigenvector with the largest overlap with the initial wave function; so, setting a proper set of 60/20 translational-/vibrational-like states, initially obtained from BIR MCTDH calculations with a small and selective basis set of SPFs, it was possible to achieve convergence toward the actual translational or vibrational eigenenergies and eigenstates of each system. Types and numbers of DVR basis functions, together with the mode assignation, can be consulted in Tables S2 and S3, respectively. Such calculations for the first five pure translationally excited states of H 2 O@C 60 endofullerene were performed, and their energies with respect to the ground state energy, together with representative probability density distributions in the translational Cartesian coordinates, are plotted in Figure 4. Translationally excited states were assigned according to the principal quantum number n = n x + n y + n z of the 3D isotropic harmonic oscillator, and each of them has a g n = (n + 1)(n + 2)/2 degeneracy. The first translationally excited state of para-H 2 O, labeled as n = 1, is a triplet (g = 3 with one quantum in each n x /n y /n z coordinate) state at an energy of 163.61 cm −1 above the ground E 0 state, corresponding to the fundamental translational excitation. This value agrees well with those previously reported at 162.08 and 163.04 cm −1 from 6D and 9D quantum calculations, 44,60 respectively. Higher excited levels with two to five quanta of translational excitation are also computed from the present 9D MCTDH calculations and assigned to n = 2 (g = 5/1), 3 (g = 7/3), 4 (g = 9/5/1), and 5 (g = 11/7/3) states at energies of 325.97/330.43, 491.37/ 498.88, 660.03, and 851.74 cm −1 , respectively. The energy splitting found between the degenerate translationally excited levels, e.g., 4.5 and 7.5 cm −1 in the n = 2 and 3 levels, respectively, reflects the anharmonicity in this mode. In a previous study, 60 values have been reported only for the second n = 2 excited states from a 6D rigid rotor model, corresponding to l = 2 and 0 at energies of 325.97 (g = 5) and 328.21 (g = 1) cm −1 , respectively. One can see that our results are in accordance with these values, while the small differences can mainly be attributed to the fully coupled 9D treatment developed in the present study.
As before, a set of lock IR MCTDH calculations were performed to calculate the vibrational energy levels of H 2 O@ C 60 using the 9D Hamiltonian. Types and numbers of DVR basis functions, together with the mode assignation, can be consulted in Tables S2 and S3 (see the SI). The computed vibrational energies for the encapsulated water molecule in the C 60 cavity are listed in Table 2, together with the corresponding probability density distributions in ( 1 ,γ) and ( 1 , 2 ) planes. One can see that convergence was achieved for specific vibrationally excited states of H 2 O@C 60 endofullerene; in total, 13 levels were computed up to energies of about 7500 cm −1 above the ground state energy. In turn, on the basis of their probability distributions, each state was assigned as ν 1 ν 2 ν 3 , using the vibrational quantum numbers of the isolated water molecule, with ν 1 and ν 3 corresponding to the symmetric and asymmetric OH stretching modes, respectively, while ν 2 to the HOH bending one. Thus, apart from the fundamentals, we were also able to assign vibrationally excited levels with up to 4 quanta in the bending mode, as well as several combinations with stretching modes. In addition, in Table 2, we also present vibrational energies, obtained here from 3D MCTDH calculations for an isolated water molecule. We found that such 3D results are in accordance with the values reported previously using the same H 2 O PES. 61,72 Just recently, quantum full 9D calculations have been reported 44 for the lower four vibrationally excited states of H 2 O@C 60 . We found that the present values agree very well with those energies, within less than 1 cm −1 .
By comparing now our energy values from the 9D and 3D calculations (see Figure 5) for the encapsulated and isolated water molecules, respectively, we found that both ν 1 and ν 3 OH vibrational fundamental frequencies and overtones of the H 2 O molecule undergo appreciable upshifts upon its encapsulation, while the ν 2 bending mode exhibits a much smaller one. In particular, OH stretching fundamentals of the C 60 -encapsulated H 2 O molecule are blue-shifted from those of the isolated gas-phase water molecule by 24 and 25 cm −1 , while the 010 HOH bending fundamental and its 020 overtone are blue-shifted by only 2.5 and 3.8 cm −1 , respectively, much less than the stretching ones. The latter values are in accordance with previously reported values of 1−6 cm −1 from PBE0/6-311++G(d,p) harmonic frequency optimization calculations 38 and 2.19 and 3.38 cm −1 from the quantum 9D calculations by Felker and Bacǐc. 44 So far, frequency blue shifts of the stretching H 2 O modes have been reported by different theoretical treatments. 35,36,38,39,44 In particular, normal mode analysis from Hartree−Fock optimizations 35   Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article which are smaller for the symmetric mode than the asymmetric mode. Further, instantaneous vibrational analysis applied to molecular dynamics simulation configurations at the B3LYP/ cc-pVDZ(STO-3G) level of theory has reported blue shifts of 13−19 cm −1 only for the symmetric stretching mode, whereas no shift has been found for the asymmetric one. 36 Finally, more recently, OH bond vibrational frequency blue shifts from the gas-phase water, of similar values of 14 and 15 cm −1 for the symmetric and asymmetric OH stretching modes, respectively, have been reported 39 from classical MD simulations using the same DFT-SAPT force-field parameters 39 as those in the present work. Our estimated blue shift values are compared quite well with those given by Farimani et al. 39 and are very close (within 1 cm −1 ) to the values of 24.2 and 23.8 cm −1 reported by Felker and Bacǐc. 44 However, we should note that in the absence of any experimental evidence on such effects upon water molecule's encapsulation into the C 60 cage, together with the contradictory results on the values of the vibrational frequency shifts, one should consider to carefully evaluate the potential models employed in the present and most recent studies, on the basis of first-principles guest−host interactions taking into account many-body contributions.

SUMMARY AND CONCLUSIONS
The effect of nanoconfinement on the quantum features of the encapsulated molecules is still an open topic, with relevance in various molecular processes in several scientific areas. In an effort to contribute to the field, a light endohedral fullerene, such as H 2 O@C 60 , was modeled by means of an exact full 9D Hamiltonian, deduced for a system of internal Radau coordinates to describe the encapsulated water molecule, which correlates all of the degrees of freedom and considers anisotropies in the rigid C 60 cage and the extended H 2 O molecule interaction. One specific focus of the present study was to develop a systematic and computationally efficient The corresponding values obtained for an isolated water molecule from 3D MCTDH calculations are also presented for comparison purposes, as well as the frequency shifts, δ ν , between the encapsulated and isolated water vibration modes. Plots of probability density distributions, in Radau coordinates, are shown, with / 1 2 bond lengths in the range of 0.7−1.3 Å and angle γ in 0.9−3.0 rad. procedure to treat quantum mechanically such a highly coupled system within the MCTDH framework. Thus, we have shown that such a quantum computational scheme could be utilized for rigorous and tractable computations of any light−heavy−light encapsulated molecule. The calculations yield well-converged intermolecular low-lying rotational and translational states, as well as intramolecular vibrational fundamentals and overtones up to higher energies of 7500 cm −1 of the H 2 O@C 60 system, and the obtained results for the encapsulated water molecule were compared with those for the isolated molecule.
The results reported can be established as the first reference values in the study of the quantum features of such an endofullerene. In some cases, we have verified, while in others, we have managed to expand the previously reported values in higher energies. We have found small energy splitting in higher rotational levels of the encapsulated water molecule that reflects the nonsphericity of the C 60 cage, while such splitting in translationally excited levels was attributed to anharmonicities in this mode. Finally, we have predicted significant blue shifts in the vibrational frequencies of the encapsulated water with respect to those of the isolated gas-phase molecule, as an important consequence of its encapsulation in the C 60 cage.
However, we do realize that conclusions in the present study, as well as in previous ones, are constrained by the limitation of the sum-of-potentials PES used. Although the results are numerically exact for the PES employed, there is evidence from electronic structure calculations of significant noncovalent interactions between the water molecule and the C 60 cage. Importantly, the lack of such many-body contributions could introduce considerable uncertainty for the PES's quality. To date, there is no ab initio/first-principlesbased interaction potential for such endofullerenes. Thus, our efforts should be directed on the accurate description of such guest−host interactions, through well-converged reference training data sets from wave function and/or reliable DFT approaches, for developing predictive data-driven potential/ force-field models in a general automated fashion by sampling the most representative and diverse set of configurations. Following current challenges in the field, in the near future, we aspire to address by proper machine-learned modeling these noncovalent interactions and employ such improved (highquality) models to carry on further molecular computer simulations, providing decisive results on fully coupled intramolecular and intermolecular excitations and conclusive insights into the impact of the guest−host forces on the quantized dynamics of the encapsulated water molecule.
Kinetic energy operator in the 9D coordinates; details on the convergence of the n-mode potential expansion and the multidimensional POTFIT approach ( Figure S1 and Table S1); and specific computational details for carrying out the MCTDH calculations (Tables S2 and  S3