H2, HD, and D2 in the small cage of structure II clathrate hydrate: Vibrational frequency shifts from fully coupled quantum six-dimensional calculations of the vibration-translation-rotation eigenstates

We report the first fully coupled quantum six-dimensional (6D) bound-state calculations of the vibration-translation-rotation eigenstates of a flexible H2, HD, and D2 molecule confined inside the small cage of the structure II clathrate hydrate embedded in larger hydrate domains with up to 76 H2O molecules, treated as rigid. Our calculations use a pairwise-additive 6D intermolecular potential energy surface for H2 in the hydrate domain, based on an ab initio 6D H2–H2O pair potential for flexible H2 and rigid H2O. They extend to the first excited (v = 1) vibrational state of H2, along with two isotopologues, providing a direct computation of vibrational frequency shifts. We show that obtaining a converged v = 1 vibrational state of the caged molecule does not require converging the very large number of intermolecular translationrotation states belonging to the v = 0 manifold up to the energy of the intramolecular stretch fundamental (≈4100 cm−1 for H2). Only a relatively modest-size basis for the intermolecular degrees of freedom is needed to accurately describe the vibrational averaging over the delocalized wave function of the quantum ground state of the system. For the caged H2, our computed fundamental translational excitations, rotational j = 0→ 1 transitions, and frequency shifts of the stretch fundamental are in excellent agreement with recent quantum 5D (rigid H2) results [A. Powers et al., J. Chem. Phys. 148, 144304 (2018)]. Our computed frequency shift of −43 cm−1 for H2 is only 14% away from the experimental value at 20 K. Published under license by AIP Publishing. https://doi.org/10.1063/1.5090573

. The second type are the large cages, eight per unit cell, in which 28 H 2 O molecules are arranged in 12 pentagonal and 4 hexagonal faces and therefore denoted 5 12 6 4 . Experiments have shown that the small cage can accommodate only one H 2 molecule, while up to four H 2 molecules can be encapsulated in the large cage. 6 Hydrogen clathrate hydrates have attracted a great deal of interest in recent years, owing to their potential as economical and environmentally friendly hydrogen storage materials. 1,2,[7][8][9][10] Moreover, hydrogen molecules entrapped in the clathrate hydrate cages constitute fascinating and unconventional chemical systems whose dynamics and spectroscopy are thoroughly dominated by strong quantum effects, to a degree matched only by light molecules inside fullerenes. 11 The pronounced quantum effects have multiple sources; one of them is the quantization of the translational center-of-mass (c.m.) degrees of freedom (DOFs) of the guest molecule(s) due to the nanoscale confinement in the clathrate cage, small or large (particle-in-a-box effect). The confining potential of the hydrate cage couples the quantized translational DOFs to the also quantized rotational DOFs of the hydrogen molecule(s). The resulting translation-rotation (TR) energy level structure is sparse, owing the low molecular mass of H 2 /HD/D 2 , their large rotational constants, and the small size of the hydrate cavities. The salient features of the TR eigenstates of a single hydrogen molecule in the cages of the sII clathrate hydrate, notably the splittings of both the translational fundamental and rotational levels, as well as their manifestations in the inelastic neutron scattering (INS) spectra, have been characterized by Bačić and co-workers through quantum 5D bound-state calculations [12][13][14][15] and rigorous computations of the corresponding INS spectra. [15][16][17][18][19] Quantum TR dynamics of multiple hydrogen molecules in the large hydrate cage has been investigated by means of the diffusion Monte Carlo (DMC) 20 and path-integral molecular dynamics (PIMD) simulations 21 and also by fully coupled eigenstate-resolved calculations. [22][23][24] In all these calculations, the hydrogen-bonded clathrate hydrate framework was treated as rigid. In a recent study, 25 this constraint was relaxed partially, by performing quantum 5D calculations of the TR levels of H 2 in the small sII hydrate cage, while taking into account the quantum delocalization of the proton nuclei of the framework water molecules arising from their hindered rotations about the fixed positions of their O atoms.
Besides giving rise to the TR energy level structure, the encapsulation of hydrogen molecules in the cages of clathrate hydrates results in the shift in the frequency of the H 2 intramolecular stretching vibration away from that in the gas phase. This frequency shift is readily observable in the Raman spectra of the binary tetrahydrofuran (THF) + H 2 sII hydrate, where the large cages are completely occupied by the THF while the small cages are singly occupied by H 2 , and simple sII hydrates in which H 2 molecules are the only guests. 10,26,27 The vibrational frequencies of H 2 molecules encapsulated in the sII hydrates are always lower than, i.e., redshifted, relative to, the gas-phase H 2 . The largest redshift (−34 cm −1 ) is observed in the Raman spectra of the THF + H 2 sII hydrate and can be assigned unambiguously to the singly H 2 occupied small cage. 10,26,27 The same redshift of −34 cm −1 appearing in the Raman spectra of the simple sII hydrate is, therefore, also attributed to H 2 in the small cage.
The Raman spectra of the simple II hydrate also display bands redshifted by − 26, −18, and −11 cm −1 , respectively, 10,26,27 that must represent contributions from the large cages whose H 2 occupancy ranges between two and four. However, associating each of these redshifts with a particular H 2 occupancy of the large cages proved to be nontrivial. Initially, the redshifts of −26, −18, and −11 cm −1 were interpreted in terms of triply, doubly, and singly occupied large cages, respectively. 26 Subsequent very careful experiments that involved multiple cycles of heating and quenching of the sample and measurements of the amounts of H 2 released in each led to the essentially opposite assignment of these three redshifts to double, triple, and quadruple occupancies of the large cages, respectively. 27 The observed trend in the H 2 redshift can be understood in terms of the interplay between two kinds of interactions. 27 One of them, the attractive interaction between H 2 and the cage, softens the intramolecular stretch potential of H 2 and lowers its vibrational frequency relative to the gas-phase. As the large-cage occupancy of the large cage increases, the tighter packing and the largely repulsive H 2 -H 2 interactions lead to the increasing vibrational frequency of H 2 and the decreasing redshift. The fact that the H 2 vibrational frequency is redshifted even for the highest, quadruple occupancy of the large cage suggests that the attractive guesthost interaction always remains dominant over the repulsive H 2 -H 2 interactions.
In the case of sII hydrogen hydrates, it has been possible to assign with confidence the observed frequency shifts to different H 2 occupancies of the small and large clathrate cages largely guided by the experimental data. But, in general, e.g., molecular hydrogen in metal-organic frameworks (MOFs), 28,29 reliable decoding of the information contained in the vibrational frequency shifts regarding the H 2 occupancies of the cavities of nanoporous materials, and other structural as well as dynamical aspects of the entrapped H 2 , requires theoretical methods capable of reliably calculating the frequency shifts. This is a highly challenging task, for two reasons. First, the problem is inherently high-dimensional. Even if the hydrate framework is treated as rigid, the dimensionality of the calculations is 6nD, where n is the number of H 2 molecules considered; thus, for n = 1 − 4, one has to be able to deal with the problem whose dimensionality ranges from 6D to 24D. This requires having accurate highdimensional potential energy surfaces (PESs) that incorporate the H 2 -clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H 2 molecules. Both interactions must include the dependence on the H 2 intermolecular stretch coordinate and its coupling to the intermolecular degrees of freedom. Second, dynamical quantum effects and anharmonicities in both intra-and intermolecular DOFs play a significant role, particularly at the low temperatures of the Raman spectroscopy measurements. Consequently, these key features have to be fully accounted for in any first-principles theoretical method aiming to generate accurate frequency shifts of encapsulated hydrogen molecules.
Within the past decade, a number of approaches, involving a variety of approximations, have been taken to address this fundamental and difficult problem. In some of them, the H 2 molecules encapsulated in the isolated small or large hydrate cages were taken to be frozen in the geometry corresponding to the minimum energy of the system. [30][31][32] As a result, nuclear quantum effects are left out, in particular, the averaging over the large-amplitude intermolecular vibrations of the guest H 2 molecules. In addition, since only The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp isolated clathrate cages are considered, the effects of the condensedmatter environment are unaccounted for. This problem has also been treated through a combination of classical molecular dynamics (MD) and PIMD simulations with electronic structure calculations at the density functional theory [DFT (B3LYP)] and MP2 levels. 33 The H 2 vibrational frequencies calculated in 1D for the H 2 intermolecular coordinates taken from many snapshots of the MD simulations covered a broad distribution of frequencies that extended to that of the free H 2 at 4155 cm −1 . Their maxima agree reasonably well with the experiment, after a scaling factor was introduced in the calculations. Finally, classical MD simulations within the DFT framework were performed for an sII hydrate unit cell, and the H 2 vibrational spectra were calculated by Fourier transforming the H-H bond length autocorrelation function. 34 This classical treatment does not account for the quantum effects. Moreover, it gives the vibrational spectra that are shifted by 100-150 cm −1 to higher frequencies relative to the experimental results and above the stretch fundamental of free H 2 . Very recently, Powers et al. 35 have calculated the frequency shift of H 2 inside the small cage of the sII hydrate, isolated and embedded in spherical hydrate domains of increasing size, in order to investigate the effect of the condensed-phase environment. The approach employed was developed earlier by Bačić and co-workers for the purpose of computing the HF stretch frequency shift in ArnHF clusters. [36][37][38][39] The H 2 frequency shift was obtained by means of the quantum 5D bound-state calculations of the coupled TR eigenstates on a pairwise-additive intermolecular PES for rigid H 2 in a (rigid) hydrate domain, which depends on the vibrational state of H 2 , v = 0 or v = 1. This 5D PES was constructed using the 5D (rigidmonomer) pair potential for the interaction of H 2 in the ground and first excited vibrational states, respectively, with H 2 O, obtained by averaging the full-dimensional (9D) ab initio PES of H 2 -H 2 O by Valiron et al. 40 over the vibrational ground state wave function of H 2 O and the vibrational wave functions of H 2 for v = 0 and v = 1, respectively. Implicit in this approach is the assumption of dynamical decoupling between the H 2 intramolecular vibration and the TR modes, well-justified by their large energy separation. The H 2 vibrational frequency shift of ∼−44 cm −1 calculated for the largest clathrate domain considered, with 1945 H 2 O molecules, which mimics the condensed-phase environment, was about 10% larger in magnitude than that obtained for the isolated small cage. This 0 K value agrees well with the frequency shifts measured at 20 K 26 (−37 cm −1 ) and at 76 K 27 (−34 cm −1 ). It was suggested that improving further the agreement with the experiment may require including many-body interactions, three-body, in particular, missing from the pairwise-additive intermolecular PES employed. 35 Motivated in part by this suggestion and other considerations, Qu and Bowman 41  The DMC method employed by Qu and Bowman 41 is wellsuited for ground-state calculations, but already the first excited state poses a challenge arising from the need to locate the node in the wave function, which is generally unknown (unless it can be determined from symmetry considerations 42 ). The calculations for the first excited vibrational state of the caged H 2 were done in the fixed-node approximation, applying the "adiabatic" method of McCoy and co-workers 43 to find the position of the node. However, determining the correct nodal surface in a six-dimensional (6D) system is very difficult, virtually impractical. Therefore, Qu and Bowman made the approximation that the node is located entirely on the H-H intramolecular stretch coordinate and is independent of the TR coordinates of H 2 , thereby reducing the search for its position to 1D. This is equivalent to the assumption that the intra-and intermolecular coordinates of the caged H 2 are decoupled, justified by the large energy separation between the two types of modes. 41 Thus, the quantum methodologies employed in the two recent computations of the H 2 vibrational frequency shift in the small sII hydrate cage, 35,41 although entirely different, both rely on the approximation of no coupling between the high-energy intramolecular vibrational mode of H 2 and its low-energy TR modes. There is no reason to doubt its validity for this system (in both approaches) and the accuracy of the results it yields. Still, one can ask whether it is possible to perform quantum 6D calculation of the bound states of H 2 in the small cage up to, and including, the energy of the first excited (v = 1) intramolecular vibrational state (the stretch fundamental), around 4100 cm −1 , treating the intra-and intermolecular (TR) degrees of freedom as fully coupled, i.e., not invoking their separability. After all, fully coupled full-dimensional (6D) quantum calculations of the vibrational levels of nonrigid molecular systems, such as (HF) 2 , 44 (HCl) 2 , 45 and CO on Cu(100), 46 have been feasible for some 25 years. In some cases, e.g., (HCl) 2 47 and CO on Cu(100), 46 the quantum 6D calculations yielded the energies of the intramolecular stretch fundamental(s), and thus their shifts from the respective gas-phase values.
Nevertheless, it has been generally thought that for molecular systems which have both high-frequency intramolecular mode(s) and low-frequency intermolecular vibrations, such as H 2 in the small sII hydrate cage, and hydrogen-bonded and van der Waals (vdW) complexes, rigorous calculation of fundamental excitation(s) of their intramolecular mode(s), e.g., the v = 1 vibrational state of the encapsulated H 2 that is 6D for a rigid cage, would be an extremely difficult and prohibitively costly task. The main source of the difficulty was the assumption that the very large number of highly excited intermolecular vibrational eigenstates in the manifold of the intramolecular ground state below the energy of the intramolecular excitation(s) all have to be converged in order to compute accurate fundamental intramolecular excitation(s). In this paper, we demonstrate that certainly for the intramolecular stretch fundamental of H 2 (HD, D 2 ) inside the small cage of the sII hydrate, and its frequency shift, this widely held view is not correct, making this problem entirely tractable.
We present the results of the fully coupled quantum 6D calculations of the vibration-translation-rotation (VTR) eigenstates of a single flexible H 2 , HD, and D 2 molecule entrapped in the (rigid) small cage of the sII hydrate. We show that computing the converged energy of the first excited (v = 1) intramolecular vibrational The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp state of the caged H 2 (and the isotopologues) at ≈4100 cm −1 requires converging only the TR states in the v = 0 manifold up to at most 400-450 cm −1 above the ground state. Guided by our previous studies, 15,35 quantum 6D calculations of the coupled VTR eigenstates, which extend to the v = 1 state and its frequency shift away from the gas-phase value, are performed for H 2 encapsulated inside the spherical sII hydrate domains of increasing radius treated as rigid. The 6D intermolecular PES for a flexible hydrogen molecule inside the hydrate domain utilized in these calculations is constructed in a pairwise additive fashion, based on an ab initio 6D H 2 -H 2 O pair potential. The TR eigenstates and vibrational frequency shifts computed for H 2 , HD, and D 2 are compared with the available experimental data, as well as the results of the quantum 5D calculations in Ref. 35. This paper is organized as follows: Methodology is described in Sec. II. In Sec. III, we present and discuss the results. Section IV summarizes the work and outlines possible extensions.

A. Clathrate hydrate domains and the ab initio 6D H 2 -H 2 O pair potential
The three spherical sII clathrate hydrate domains used in this work are identical to those generated previously by Bačić and coworkers 15 and employed in our recent quantum 5D H 2 frequency shift calculations. 35 They are carved out of the 3 × 3 × 3 supercell of the sII hydrate. 15 The three domains of increasing size and number of H 2 O molecules N have the cutoff radii set to (a) 5.0 Å, enclosing only the N = 20 water molecules of the small dodecahedra cage itself; (b) 7.5 Å, encompassing additional 20 H 2 O molecules, for the total of N = 40 water molecules; and (c) 9.0 Å, encompassing N = 76 water molecules that form the first three complete solvation shells around H 2 . 15 In our previous study, 35 the largest hydrate domain considered was much larger and included 1945 water molecules. However, that study also showed that the frequency shift computed for this domain differs by only 3% from that obtained for the N = 76 domain. Therefore, in the present work, we do not go beyond the sII hydrate domain containing N = 76 water molecules. In the bound-state calculations, the domains are taken to be rigid.
For the 6D intermolecular PES of flexible H 2 inside the domain with N water molecules, denoted V H 2 −domain , only one-and twobody terms of its many-body expansion are retained is the one-body term for the intramolecular stretching coordinate (r) of the hydrogen molecule. For it, we use the corresponding one-body term in the many-body representation of the PES for H 2 in the sII hydrate by Bowman and co-workers. 48 The second term in Eq. (1) represents the summation over the two-body interactions V Ξi) between the hydrogen molecule and each of the N water molecules in the domain. The coordinates of H 2 are q h ≡ {R, ω, r}, where R is the vector pointing from the center of the confining small cage, which is also the origin of the space-fixed (SF) Cartesian frame attached to the cage, to the center of mass (c.m.) of H 2 , while ω ≡ (θ, ) are the polar and azimuthal angles, respectively, which specify the orientation of H 2 relative to the SF frame. The position vector R can be expressed either in terms of the SF Cartesian coordinates {x, y, z} 12,35 or the spherical polar coordinates {R, }, where R ≡ ||R||, and ≡ (Θ, Φ) are the polar and azimuthal angles of R relative to the SF axes. 49 Ξi denotes the coordinates of the ith water molecule in the domain; these are fixed since the domains are assumed to be rigid.
The 6D pair potential 40 In this 9D PES, the flexibility of H 2 and H 2 O monomers is included as a correction to the rigid-monomer dimer 5D PES, computed as a Taylor expansion around the equilibrium geometries of monomers. The 9D PES V08 was obtained by combining standard coupled-cluster with single, double, and with the perturbative contributions of the triples excitations [CCSD(T)] calculations with the explicitly correlated CCSD(T)-R12 calculations, and it is expected to provide the currently most accurate description of the H 2 -H 2 O interaction, with an accuracy of a few cm −1 in the region of the van der Waals minimum. 40 In this work, the 9D H 2 -H 2 O PES was reduced to the 6D pair potential, for flexible H 2 and rigid H 2 O, by fixing the intramolecular coordinates of H 2 O to their values in the ground vibrational state (OH bond length = 1.843 a 0 and the HOH bending angle = 104.41 ○ ). The accuracy of this procedure is comparable to that of averaging the 9D PES over the vibrational ground-state wave function of H 2 O. Finally, the intermolecular coordinates employed in the V08 potential are transformed numerically to the coordinates used for the V

B. Quantum 6D diffusion Monte Carlo calculations
Although the focus of this study is on the excited VTR eigenstates, we also use the diffusion Monte Carlo (DMC) method to compute in 6D the VTR ground-state energy of flexible H 2 , HD, and D 2 inside the rigid water domains. This approach simulates a diffusion process in imaginary time on a given PES. The general DMC approach has been described in detail in Ref. 50 and here we use a standard (i.e., not rigid body) formulation of the algorithm for the caged molecule, while the cage itself remains fixed.
The simulations are performed using a revised parallelized version of the XDMC code developed by Benoit 51 (see also Ref. 50 for implementation details). For each simulation in this study, we use 1000 replicas, stabilization periods of 61 500 cycles (H 2 ), 80 900 cycles (HD), and 108 300 cycles (D 2 ) with ∆τ = 5 a.u. and an averaging phase of 1000 × 100 cycles with ∆τ = 1 a.u.

C. Quantum 6D calculations of the coupled vibration-translation-rotation eigenstates
The 6D Hamiltonian for the coupled VTR motions of a vibrating diatomic molecule AB, which in this study corresponds to H 2 and its isotopologues HD and D 2 , inside a rigid clathrate hydrate domain, can be written aŝ where mAB and µAB are the total mass and the reduced mass of AB, respectively, while R, ω, and r were defined in Sec. II A. ∇ 2 is the Laplacian associated with R, andĵ 2 is the operator associated with the square of the rotational angular momentum of AB. For the isotopic masses of hydrogen, the values m H = 1.008 g mol −1 and m D = 2.0141 g mol −1 were used.

The Smolyak scheme approach with ELVIBROT
In most, although not all, of the calculations performed in this study, the AB c.m. position vector R in the Hamiltonian in Eq. (2) is expressed in terms of the Cartesian coordinates {x, y, z}, and ω ≡ (θ, ). Furthermore, the operatorĵ 2 of Eq. (2) is expanded in terms of partial derivative operators, ∂ ∂θ and ∂ ∂φ . For this choice of {x, y, z, θ, , r} coordinates, we have used the Smolyak scheme approach 52 introduced by Avila and Carrington [53][54][55] and also proposed by Lauvergnat and Nauts. [56][57][58] More recently, it has been used to calculate the energy levels of H 2 in a clathrate hydrate. 25,35 In the Smolyak scheme, the single large direct-product basis or grid is replaced by a sum of small direct-products, denoted as S rep L S , where S i i represents the ith primitive basis or grid. The parameter i defines the size of this primitive basis, nbi( i), or grid, nqi( i), as shown in Table I

Filter diagonalization in a direct-product basis
An alternative approach employed in this study is to compute the VTR eigenstates of the 6D Hamiltonian in Eq. (2) in selected regions of the energy spectrum, utilizing the Chebyshev variant 60 of filter diagonalization, 61 together with the direct-product basis described below. The AB c.m. position vector R in Eq. (2) is now expressed in terms of the spherical polar coordinates {R, }, with R ≡ ||R||, and ≡ (Θ, Φ), so that the complete set of coordinates is {R, Θ, Φ, θ, , r}. The matrix representation of the Hamiltonian is formed in a basis consisting of the product functions n, l, m l , j, m, rγ⟩ ≡ n, l, m l ⟩ j, m⟩ rγ⟩, where n = 0, 1, . . ., nmax, l = n, n − 2, . . ., ≥ 0, |m l | = 0, 1, . . ., l, j = 0, 1, . . ., jmax, |m| = 0, 1, . . ., j, and γ = 1, . . ., γmax. Here, are eigenfunctions of the 3D isotropic HO (e.g., see the supplementary material in Ref. 49) having the angular frequency β/mAB, ⟨θ, φ j, m⟩ ≡ Y m j (θ, φ), and the |rγ⟩ constitute a discrete variable representation (DVR) 62 derived from the eigenfunctions of a 1D oscillator of mass µAB moving in a Morse potential of the form where D, α, and re are constants chosen so that V Morse (r) ≃ V (1b) h (r) in Eq. (1). The specific parameters that we have used in conjunction with this basis (all in atomic units) are as follows: D = 0.1744, α = 1.02764, re = 1.40201, and β = 2.9889. We note that because of the Pauli principle, j can be either even (para-H 2 ) or odd (ortho-H 2 ).
The Chebyshev variant of filter diagonalization requires the repeated application ofĤ 6D in Eq. (2) on an initial, random state vector. This is readily accomplished by matrix-vector multiplication for the kinetic-energy portion of Eq. (2). The ∇ 2 andĵ 2 /r 2 parts ofĤ 6D have analytic matrix elements in the basis of Eq. (4). The matrix elements of ∂ 2 ∂r 2 are diagonal in all the basis-set indices except γ, and the ⟨r γ ′ ∂ 2 ∂r 2 rγ⟩ can be straightforwardly obtained by numerical transformation of the matrix elements from the Morseeigenvector representation to the Morse-DVR one. To operate with the potential-energy portion of Eq. (2), we first transform the state vector to a grid representation |Rρ, (Θ,Φ) ξ , (θ, )η, rγ⟩, where the Rρ (ρ = 1, . . ., Nρ) are associated-Laguerre quadrature points and the (Θ, Φ) ξ (ξ = 1, . . ., N ξ ) and the (θ, )η (η = 1, . . ., Nη) are Lebedev quadrature points. We then multiply the state vector at each grid point with the value of V H 2 -domain at that grid point. Finally, we transform the result back to the |n, l, m l , j, m, rγ⟩ representation. In this study, the grid sizes Nρ, N ξ , and Nη are 10, 110, and 38, respectively.

Obtaining converged first excited vibrational state of H 2 and the isotopologues
As outlined in the Introduction, it has been generally assumed that a converged fully coupled quantum 6D calculation of the highenergy v = 1 vibrational state of the caged H 2 would necessarily involve converging a very large number of the VTR states in the v = 0 manifold lying below it. This, of course, would require diagonalization of a prohibitively large matrix of the 6D VTR Hamiltonian, making the task virtually intractable.
However, convergence tests performed for the quantum 6D calculations utilizing the Smolyak scheme approach of the VTR levels of H 2 inside the small dodecahedral cage with 20 water molecules reveal an entirely different picture. The calculations using different LS values for the basis set (LB = 6) and for the grid (LG = 7), which generate 8246 basis functions and 460 000 grid points, yielded 4120.9 cm −1 as the energy of the first excited (v = 1) vibrational state of H 2 (for the TR DOFs in the ground state). The results shown in Tables III and IV are obtained utilizing this basis. Increasing LB to 7 and LG to 8 gives 17 900 basis functions and 1 167 282 grid points. Despite more than doubling the basis set size, the energies of the v = 1 vibrational state calculated for LB = 6 and LB = 7 differ by less than 0.1 cm −1 , indicating its "high degree" of convergence. In contrast, comparison of the results of the two calculations shows that the highly excited TR states in the v = 0 manifold, close in energy to the v = 1 state, are not converged at all. In fact, only the v = 0 TR states with excitation energies up to 400-450 cm −1 , far below the v = 1 state, are converged to within 1 cm −1 .
Thus, what emerges from these calculations is the unexpected result that obtaining a well-converged energy of the v = 1 vibrational state of H 2 does not require having converged high-lying v = 0 TR states in its vicinity and below. This suggests that the latter are very weakly coupled to the H 2 stretch fundamental, and therefore, it confirms the validity of our previous 5D results. 35 Quantum 6D calculations on the same system using filter diagonalization and a direct-product basis described above confirm this finding and the conclusion, and go one step further. The basis set parameters (nmax, lmax, jmax, γmax) ranging from (7,7,4,8) to (10,6,4,11) give rise to basis sets ranging in size from 14 400 to 30 030. Although differing in size by more than a factor of two, these basis sets, when used in the quantum 6D calculations, all give the energies of the v = 1 vibrational state that are to within 0.1 cm −1 of each other and converge on 4121.1 cm −1 . This result is very close to 4120.9 cm −1 , the value obtained for the v = 1 state with Smolyak scheme approach.
A very interesting feature of the results of the filterdiagonalization calculations is that, owing primarily to the low jmax = 4, the computed TR states of the v = 0 manifold extend only up to excitation energies of 800 to 1500 cm −1 , depending on the size of the basis, some 2500 cm −1 below the energy of the v = 1 state. Thus, there are no v = 0 TR states close in energy to the v = 1 state. Moreover, only the TR states up to about 200 cm −1 above the ground state are well-converged. Despite that, the calculated energy of the H 2 stretch fundamental is virtually identical, to within 0.2 cm −1 , to that obtained by the Smolyak scheme approach, which does generate TR states of the v = 0 manifold, albeit not converged, in the neighborhood of the v = 1 state.
The surprising conclusion that emerges from the Smolyak scheme approach and filter-diagonalization calculations above is that the converged first excited state of H 2 stretch can be obtained without (a) converging all TR states in the v = 0 manifold up to its energy, or (b) having any highly excited v = 0 TR states at all within a couple of thousands of wave numbers. This finding points out to extremely weak coupling between the v = 1 vibrational state of H 2 and the high-lying v = 0 TR states. We do not have a formal theoretical explanation for this weak coupling at the present time. However, the disparity between the nodal patterns of the states involved, completely irregular for the highly excited TR states vs a smooth one, with a single node for the v = 1 state, are likely to figure prominently in any theoretical model.
Both the Smolyak scheme approach and filter-diagonalization calculations also demonstrate that in order to compute a highly converged H 2 stretch fundamental, one (only) needs to use a basis for the intermolecular DOFs that can provide an accurate description of the vibrational averaging over the large-amplitude TR motions in the delocalized quantum ground state of the system. It should be stressed that this basis is much smaller than the one that would be required to get converged highly excited TR eigenstates in the vicinity of the v = 1 state. That accurate calculation of the vibrational shift in systems dominated by quantum effects must take into account averaging over the large-amplitude motions which was demonstrated first for the ArnHF clusters. [36][37][38][39] III. RESULTS AND DISCUSSION  In Table III, we report the energies of the fundamental translational excitations and the rotational j = 0 → 1 transitions of H 2 in the ground vibrational state (v = 0), as well as the frequency ν of the H 2 stretch fundamental (v = 1) and the corresponding frequency shift ∆ν, for three sII hydrate domains with the number of water molecules N equal to 20 (isolated small cage), 40, and 76. These results are from the quantum 6D calculations employing the Smolyak scheme approach. Shown for comparison are the corresponding results from the quantum 5D calculation for the same hydrate domains in Ref. 35, in which the H 2 molecule is treated as rigid, and the available experimental data pertaining to these quantities. 26,27,63 Even a cursory inspection of Table III reveals a striking agreement between the results of the fully coupled quantum 6D calculations with those from the quantum 5D, rigid-H 2 treatment, for all three domains. For the translational excitations, the agreement is better than 0.5 cm −1 , while the 6D and 5D rotational excitations agree to about 1 cm −1 .
Excellent agreement between the 6D and 5D calculations in Table III extends to the H 2 vibrational frequency shift as well. In our quantum 6D calculations, the frequency shift ∆ν for a given hydrate domain is obtained as the difference between the frequency ν of the H 2 stretch fundamental computed in 6D for this domain and the free-H 2 stretch frequency ν free evaluated for the one-body potential V (1b) h (r) in Eq. (1) (the ν free values for H 2 , HD, and D 2 are given in Table IV). For the three hydrate domains considered, the difference between the frequency shifts computed in 6D and 5D is very small, less than 0.5 cm −1 . This confirms the remarkably high accuracy of the quantum 5D method for computing vibrational frequency shifts used in Ref. 35 and earlier. [36][37][38][39] Moreover, the fact that all results of the 6D and 5D calculations, translational and rotational excitations and frequency shifts, agree so exceedingly well points to a high degree of decoupling between the high-frequency H 2 intramolecular stretch vibration and the low-frequency intermolecular TR modes. Table IV displays the results for H 15 The splitting of both the translational fundamental and j = 0 → 1 transition for all three isotopologues into three components, evident in Table IV (and also in Table III), has been discussed previously, 12,13,15,35 and is caused by the anisotropies, radial and angular, respectively, of the cage environment. For the isotopologues for which the experimental data are available, H 2 and HD, the calculated splittings are, in general, substantially larger than the measured values, especially for the rotational j = 0 → 1 transition. The implication is that the pairwise-additive 6D PES employed overestimates the anisotropies of the H 2 -hydrate interaction, angular anisotropy, in particular. This has been attributed to nonadditive many-body interactions 13 41,48 and see what degree this would improve the agreement between the calculated and measured TR excitation energies.
The theoretical frequency shifts ∆ν reported in Table IV for H 2 , HD, and D 2 are obtained as the difference between the frequency ν of the stretch fundamental calculated in 6D for each of the isotopologues inside the domain with 76 H 2 O molecules and its gas-phase value ν free (in Table IV) evaluated for the one-body potential V (1b) h (r) in Eq. (1). The vibrational frequency shifts computed in this way are −41.9 cm −1 for H 2 , −36.3 cm −1 for HD, and −30.3 cm −1 for D 2 . As expected, the shifts decrease in magnitude with the increasing mass of the isotopologue. However, the ratio |∆ν/ν| is virtually constant for H 2 , HD, and D 2 and equal to 0.010. In other words, |∆ν| ≈ 1% of ν for the three isotopologues.
As mentioned earlier, the frequency shift from the quantum 5D calculations for H 2 in the small cage within the domain with 76 water molecules is ∼3% smaller in magnitude than the shift computed for a very large domain encompassing 1945 water molecules aimed to mimic the bulk sII hydrate. 35 If we assume that the same relationship holds for the frequency shifts obtained with the quantum 6D calculations in this study, and for all three isotopologues, then the theoretical frequency shifts for H 2 , HD, and D 2 in the small cage of sII hydrate are −43, −37, and −31 cm −1 , respectively.
The vibrational frequency shift measured at 76 K for H 2 in the small cage of sII hydrate 27 (−34 cm −1 ) is about 21% smaller by magnitude than the (extrapolated) theoretical value of −43 cm −1 , while the shift measured for D 2 , also at 76 K 27 (−25 cm −1 ) is about 24% smaller in magnitude than the theoretical result of −31 cm −1 . Experimental data regarding the frequency shift of HD in the sII hydrate are not available. The agreement between theory and experiment is satisfactory, given that the shifts are computed rigorously and from first-principles, with no adjustable parameters. The agreement improves if the temperature dependence of the frequency shift measured for H 2 (but not HD and D 2 so far) in the small sII hydrate cage 26 is taken into account. As pointed out by Qu and Bowman,41

IV. CONCLUSIONS
We performed fully coupled quantum 6D calculations of the vibration-translation-rotation (VTR) eigenstates of a flexible H 2 , HD, and D 2 molecule inside the small cage of the sII clathrate hydrate, taken to be rigid. These calculations utilized two different approaches, the Smolyak scheme approach 52-58 and the Chebyshev variant 60 of filter diagonalization, 61 together with the direct-product basis described in Ref. 49. It was demonstrated that with both approaches, it is entirely feasible to obtain a highly converged energy of the first excited (v = 1) intramolecular vibrational state of the caged diatomic molecule, and its frequency shift relative to the gasphase value, without excessive computational effort. What made this possible was the realization that to obtain the converged intramolecular stretch fundamental of the entrapped H 2 at ≈4100 cm −1 it sufficed to have converged only the TR states in the v = 0 manifold up to at most 400-450 cm −1 above the ground state, necessary for a proper description of the delocalized ground state of the system and the vibrational averaging over the large-amplitude TR motions. This led to the conclusion that the v = 1 intramolecular vibrational state is extremely weakly coupled to the highly excited v = 0 TR states.
Quantum 6D calculations of the coupled VTR eigenstates, including the v = 1 state and its frequency shift relative to the gasphase value, were performed for H 2 , HD, and D 2 encapsulated inside three spherical sII hydrate domains of increasing radius, treated as rigid. A pairwise-additive 6D intermolecular PES for H 2 inside the hydrate domain was employed in these calculations, constructed using the ab initio-based 40 6D H 2 -H 2 O pair potential, for flexible H 2 and rigid H 2 O. In addition, the VTR ground state of H 2 in the (rigid) small cage was determined by means of the 6D DMC simulations, to partly verify the correctness of the eigenstate-resolved calculations.
All results of the quantum 6D calculations for H 2 in the three hydrate domains considered agree extremely well with those from the quantum 5D, rigid-H 2 treatment, 35 demonstrating the high accuracy of the quantum 5D method for computing vibrational frequency shifts employed in Ref. 35 and earlier applications. [36][37][38][39] Comparison of the quantum 6D frequency shifts for H 2 and D 2 with the corresponding experimental results at 76 K 27 shows that the latter are 21% and 24% smaller in magnitude, respectively. The difference in the magnitudes of the calculated H 2 frequency shift and that measured for H 2 at 20 K 26 is only 14%. The agreement between theory and experiment is satisfactory, but clearly there is room for improvement.
The quantum 6D calculation of the vibrational frequency shift is rigorous and yields results that are virtually exact numerically for the PES employed. The only remaining dynamical approximation is treating the hydrogen-bonded water framework as rigid. However, the H 2 -hydrate interaction is weak, and moreover, the disparity between the masses of H 2 (and isotopologues) and the confining hydrate is large. As a result, the coupling of VTR motions of H 2 to the vibrations (phonons) of the host water framework is weak as well, and its neglect (by treating the hydrate as rigid) is not expected to introduce significant errors in the calculated frequency shifts. Consequently, the main source of the residual differences between the computed and experimental values have to be certain shortcomings in the pairwise-additive intermolecular PES for the H 2 -hydrate interaction, having to do with either the ab initio 6D H 2 -H 2 O pair potential or the absence of the nonadditive three-body interactions, or a combination of both. These possibilities will be investigated in the future.
It is likely that very weak coupling between the high-frequency intramolecular modes and the low-frequency intermolecular vibrations is the feature of other molecular systems, in particular, the weakly bound ones, e.g., hydrogen-bonded and van der Waals complexes mentioned in the Introduction. In that case, the fundamental excitations of their intramolecular modes, and frequency shifts, could be calculated accurately from full-dimensional, fully coupled quantum bound-state calculations, without converging the very large number of highly excited intermolecular vibrational eigenstates in the manifold of the intramolecular ground state. This could be achieved with a relatively small basis for the intermolecular DOFs capable of accurately describing the vibrational averaging over the large-amplitude intermolecular motions in the delocalized groundstate wave function, but not the high-lying intermolecular eigenstates. One attractive and challenging target for such a treatment is the HF dimer, for which high-quality full-dimensional (6D) PESs are available, 44,[64][65][66] as is the wealth of spectroscopic data about its intraand intermolecular excitations.