We report the first fully coupled quantum six-dimensional (6D) bound-state calculations of the vibration-translation-rotation eigenstates of a flexible H_{2}, HD, and D_{2} molecule confined inside the small cage of the structure II clathrate hydrate embedded in larger hydrate domains with up to 76 H_{2}O molecules, treated as rigid. Our calculations use a pairwise-additive 6D intermolecular potential energy surface for H_{2} in the hydrate domain, based on an *ab initio* 6D H_{2}–H_{2}O pair potential for flexible H_{2} and rigid H_{2}O. They extend to the first excited ($v$ = 1) vibrational state of H_{2}, 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 translation-rotation states belonging to the $v$ = 0 manifold up to the energy of the intramolecular stretch fundamental (≈4100 cm^{−1} for H_{2}). 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 H_{2}, 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 H_{2}) results [A. Powers *et al.*, J. Chem. Phys. **148**, 144304 (2018)]. Our computed frequency shift of −43 cm^{−1} for H_{2} is only 14% away from the experimental value at 20 K.

## I. INTRODUCTION

Hydrogen clathrate hydrates are crystalline inclusion compounds where one or more hydrogen molecules are encapsulated inside the cavities, or cages, created by the three-dimensional (3D) framework of hydrogen-bonded water molecules.^{1–3} Simple hydrogen clathrate hydrates, which have only hydrogen molecules as guests, were first identified by Dyadin *et al.*^{4} and later characterized by Mao *et al.* in more detail.^{5} They have the classical structure II (sII),^{1,2,5} whose cubic unit cell consists of two types of cages. One of them is the small dodecahedral cage, sixteen per unit cell, each comprising 20 H_{2}O molecules forming 12 pentagonal faces, hence designated 5^{12}. 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–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–15} and rigorous computations of the corresponding INS spectra.^{15–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–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 guest-host 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 6*n*D, 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 high-dimensional 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–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 isolated clathrate cages are considered, the effects of the condensed-matter 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 Ar_{n}HF clusters.^{36–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 (rigid-monomer) 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} have included *ab initio* 3-body H_{2}–H_{2}O–H_{2}O interactions, in addition to the 2-body H_{2}–H_{2}O interactions, in their diffusion Monte Carlo (DMC) calculations of the vibrational frequency shift of H_{2} encapsulated in the (rigid) small cage of the sII hydrate, without and with surrounding water molecules. For the largest hydrate domain considered having 172 H_{2}O water molecules, calculations based on 2-body interactions only yielded the frequency shift of about −26 cm^{−1}, while the inclusion of the 3-body interactions resulted in the shift of −40 ± 4 cm^{−1}, in good agreement with the experiment, in particular, −37 cm^{−1} at 20 K.^{26}

The DMC method employed by Qu and Bowman^{41} is well-suited 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 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.

## II. COMPUTATIONAL METHODOLOGY

### 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 co-workers^{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 $VH2\u2212domain$, only one- and two-body terms of its many-body expansion are retained

Here, $Vh(1b)(r)$ 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 $Vh,wi(2b)\u2261VH2\u2212H2O(qh,\Xi 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

*i*th water molecule in the domain; these are fixed since the domains are assumed to be rigid.

The 6D pair potential $VH2\u2212H2O(qh,\Xi i)$ is derived from the accurate full-dimensional (9D) *ab initio* H_{2}–H_{2}O pair potential V08 of Valiron *et al.*^{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 $Vh,wi(2b)(qh,\Xi i)$ term in Eq. (1).

### 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 as

where *m*_{AB} 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.

#### 1. 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, $\u2202\u2202\theta $ and $\u2202\u2202\varphi $. For this choice of {

*x*,

*y*,

*z*,

*θ*,

*ϕ*,

*r*} coordinates, we have used the Smolyak scheme approach

^{52}introduced by Avila and Carrington

^{53–55}and also proposed by Lauvergnat and Nauts.

^{56–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 $SLSrep$,

where $S\u2113ii$ represents the *i*th primitive basis or grid. The parameter *ℓ*_{i} defines the size of this primitive basis, *nb*_{i}(*ℓ*_{i}), or grid, *nq*_{i}(*ℓ*_{i}), as shown in Table I and |*L*| = *∑ℓ*_{i}. The size of the nondirect product grid or basis $SLSrep$ in Eq. (3) is determined through the parameter *L*_{S}. In this study, five (*n* = 5) types of primitive basis sets are required: 4 harmonic-oscillator (HO) basis sets for the description of the vibrational and translational DOFs of AB (associated with the *r*, *x*, *y*, *z* coordinates) and spherical harmonics $Yjm(\theta ,\varphi )$ for the rotational DOFs of AB (coordinates *θ* and *ϕ*). The corresponding primitive grids are, respectively, the Gauss-Hermite quadrature for the HO basis sets and the Lebedev grid points for the spherical harmonics. In order to minimize the number of basis functions, the HO basis sets are scaled such that the arguments of the basis are $ui=si\u22c5(Qi\u2212Qi0)$, where *s*_{i} and $Qi0$ are the scaling parameter and the center, respectively, of the *i*th basis set. In this work, $Qi0$ has the values of 1.41 bohrs for the intramolecular stretch coordinate and zero for the translational DOFs, for all three isotopologues. For the translational DOFs (*i* = 1 − 3), the HO scaling parameter *s*_{i} is chosen to be 1.2 for H_{2} and 1.3 for HD and D_{2}. For the intramolecular stretching mode (*i* = 4), the scaling parameter *s*_{4} has the values 4.4, 4.7, and 5.2 for H_{2}, HD, and D_{2}, respectively.

ℓ_{i}
. | 0 . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . |
---|---|---|---|---|---|---|---|---|---|

Translation (HO) | |||||||||

nb_{i} | 1 | 3 | 5 | 7 | 9 | 11 | 13 | 15 | 17 |

Vibration (HO) | |||||||||

nb_{i} | 1 | 4 | 7 | 10 | 13 | 16 | 19 | 22 | 25 |

Rotation ($Yjm$) | |||||||||

j_{max} | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |

nb_{5} | 1 | 4 | 9 | 16 | 25 | 36 | 49 | 64 | 81 |

nq_{5} (Lebedev) | 6 | 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 |

ℓ_{i}
. | 0 . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . |
---|---|---|---|---|---|---|---|---|---|

Translation (HO) | |||||||||

nb_{i} | 1 | 3 | 5 | 7 | 9 | 11 | 13 | 15 | 17 |

Vibration (HO) | |||||||||

nb_{i} | 1 | 4 | 7 | 10 | 13 | 16 | 19 | 22 | 25 |

Rotation ($Yjm$) | |||||||||

j_{max} | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |

nb_{5} | 1 | 4 | 9 | 16 | 25 | 36 | 49 | 64 | 81 |

nq_{5} (Lebedev) | 6 | 6 | 14 | 26 | 38 | 50 | 74 | 86 | 110 |

#### 2. 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

where *n* = 0, 1, …, *n*_{max}, *l* = *n*, *n* − 2, …, ≥ 0, |*m*_{l}| = 0, 1, …, *l*, *j* = 0, 1, …, *j*_{max}, |*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 *β*/*m*_{AB}, $\u27e8\theta ,\varphi |j,m\u27e9\u2261Yjm(\theta ,\varphi )$, 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 *r*_{e} are constants chosen so that $VMorse(r)\u2243Vh(1b)(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, *r*_{e} = 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 $\u22022\u2202r2$ are diagonal in all the basis-set indices except *γ*, and the $\u27e8r\gamma \u2032|\u22022\u2202r2|r\gamma \u27e9$ can be straightforwardly obtained by numerical transformation of the matrix elements from the Morse-eigenvector 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 $VH2-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.

#### 3. 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 high-energy $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 *L*_{S} values for the basis set (*L*_{B} = 6) and for the grid (*L*_{G} = 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 *L*_{B} to 7 and *L*_{G} 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 *L*_{B} = 6 and *L*_{B} = 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 (*n*_{max}, *l*_{max}, *j*_{max}, *γ*_{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 filter-diagonalization calculations is that, owing primarily to the low *j*_{max} = 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 Ar_{n}HF clusters.^{36–39}

## III. RESULTS AND DISCUSSION

Figure 1 displays the 3D isosurface plot of the spatial distribution of the H_{2} molecule encapsulated in the small cage of the sII hydrate, from the 6D DMC calculation. It is clear from it that the wave function of the caged H_{2} is delocalized already in the ground state. The 3D spatial distribution is nearly spherical, reflecting the weakly hindered rotation of H_{2} inside the cage. The DMC-calculated ground-state energies of H_{2}, HD, and D_{2} in the small cage are shown in Table II. Their positive values result from the relative magnitudes of two contributions: (i) the interaction between H_{2} and the cage: this contribution is negative since the reference energy corresponds to H_{2} at a large distance from the cage; and (ii) the zero-point energy (ZPE) of the H_{2} intramolecular vibration, about 2179 cm^{−1} for the free H_{2}. For H_{2}, the ground-state energy from the DMC calculation (1441 ± 10 cm^{−1}) compares favorably with the *E*_{0} value in Table III (1438.3 cm^{−1}) computed for H_{2} in the small cage using the Smolyak scheme approach. Good agreement between these two results provides additional confirmation of the accuracy of the Smolyak scheme approach.

Method . | H_{2}
. | HD . | D_{2}
. |
---|---|---|---|

DMC | 1441 ± 10 | 1127 ± 9 | 770 ± 7 |

Smolyak | 1438.3 ± 0.1 | 1129.9 ± 0.1 | 771.2 ± 0.1 |

Method . | H_{2}
. | HD . | D_{2}
. |
---|---|---|---|

DMC | 1441 ± 10 | 1127 ± 9 | 770 ± 7 |

Smolyak | 1438.3 ± 0.1 | 1129.9 ± 0.1 | 771.2 ± 0.1 |

. | . | N
. | ||
---|---|---|---|---|

. | Expt. . | 20 . | 40 . | 76 . |

E_{0} | … | 1438.3 | 1401.7 | 1382.6 |

Translations | ||||

I | 71.0 | 66.8 (66.78) | 66.5 (66.36) | 66.2 (66.13) |

II | 80.2 | 76.1 (76.02) | 75.5 (75.35) | 75.3 (75.12) |

III | 101.1 | 93.3 (93.13) | 92.5 (92.34) | 92.2 (92.06) |

Rotations | ||||

j = 1 | ||||

I | 110.0 | 85.4 (86.61) | 93.5 (94.68) | 97.7 (98.95) |

II | 116.5 | 121.2 (122.18) | 121.1 (122.17) | 118.4 (119.40) |

III | 122.1 | 147.6 (148.65) | 140.2 (141.39) | 137.7 (138.85) |

ν | 4120.9 | 4119.4 | 4119.2 | |

Δν | −34 (76 K)/−37 (20 K) | −40.2 (−39.81) | −41.7 (−41.29) | −41.9 (−42.30) |

. | . | N
. | ||
---|---|---|---|---|

. | Expt. . | 20 . | 40 . | 76 . |

E_{0} | … | 1438.3 | 1401.7 | 1382.6 |

Translations | ||||

I | 71.0 | 66.8 (66.78) | 66.5 (66.36) | 66.2 (66.13) |

II | 80.2 | 76.1 (76.02) | 75.5 (75.35) | 75.3 (75.12) |

III | 101.1 | 93.3 (93.13) | 92.5 (92.34) | 92.2 (92.06) |

Rotations | ||||

j = 1 | ||||

I | 110.0 | 85.4 (86.61) | 93.5 (94.68) | 97.7 (98.95) |

II | 116.5 | 121.2 (122.18) | 121.1 (122.17) | 118.4 (119.40) |

III | 122.1 | 147.6 (148.65) | 140.2 (141.39) | 137.7 (138.85) |

ν | 4120.9 | 4119.4 | 4119.2 | |

Δν | −34 (76 K)/−37 (20 K) | −40.2 (−39.81) | −41.7 (−41.29) | −41.9 (−42.30) |

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 $Vh(1b)(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–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.

. | H_{2}
. | HD . | D_{2}
. | |||
---|---|---|---|---|---|---|

. | Expt. . | Theory . | Expt. . | Theory . | . | Theory . |

Translations | ||||||

I | 71.0 | 66.2 | 53.2 | 48.9 | 39.8 | |

II | 80.2 | 75.3 | 58.7 | 54.7 | 46.4 | |

III | 101.1 | 92.2 | 70.6 | 79.5 | 60.3 | |

Rotations | ||||||

j = 1 | ||||||

I | 110.0 | 97.7 | 87.6 | 62.6 | 40.6 | |

II | 116.5 | 118.4 | 93.4 | 91.4 | 61.9 | |

III | 122.1 | 137.7 | 98.5 | 106.9 | 77.7 | |

ν | 4119.2 | 3595.5 | 2962.9 | |||

ν_{free} | 4161.1 | 3631.8 | 2993.2 | |||

Δν | −34 (76 K)/−37 (20 K) | −41.9 | −36.3 | −25 (76 K) | −30.3 |

. | H_{2}
. | HD . | D_{2}
. | |||
---|---|---|---|---|---|---|

. | Expt. . | Theory . | Expt. . | Theory . | . | Theory . |

Translations | ||||||

I | 71.0 | 66.2 | 53.2 | 48.9 | 39.8 | |

II | 80.2 | 75.3 | 58.7 | 54.7 | 46.4 | |

III | 101.1 | 92.2 | 70.6 | 79.5 | 60.3 | |

Rotations | ||||||

j = 1 | ||||||

I | 110.0 | 97.7 | 87.6 | 62.6 | 40.6 | |

II | 116.5 | 118.4 | 93.4 | 91.4 | 61.9 | |

III | 122.1 | 137.7 | 98.5 | 106.9 | 77.7 | |

ν | 4119.2 | 3595.5 | 2962.9 | |||

ν_{free} | 4161.1 | 3631.8 | 2993.2 | |||

Δν | −34 (76 K)/−37 (20 K) | −41.9 | −36.3 | −25 (76 K) | −30.3 |

Table IV displays the results for H_{2}, HD, and D_{2} encapsulated inside the small cage embedded within the sII hydrate domain having 76 water molecules, from the quantum 6D calculations employing the Smolyak scheme approach. The computational results shown for H_{2} and the isotopologues include the fundamental translational excitations and the rotational *j* = 0 → 1 transitions in the $v$ = 0 state, and the fundamental stretch frequency *ν*, and its shift Δ*ν*. The pertinent experimental data are included for comparison. As explained in Sec. II A, the sII hydrate domain with 76 H_{2}O molecules, which complete the first three hydration shells surrounding the hydrogen molecule,^{15} is chosen primarily because the magnitude of the H_{2} frequency shift computed previously (in 5D) for this domain^{35} is only 3% smaller than the shift calculated for the largest domain comprising 1945 H_{2} molecules, intended to mimic the bulk sII hydrate. Thus, the domain of this size captures virtually all of the condensed-phase effect on the vibrational frequency shift and on the TR excitations.^{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} that are missing from this pairwise-additive PES. Therefore, it would be interesting to repeat these quantum 6D calculations for the H_{2}-hydrate PES that, in addition to the 2-body H_{2}–H_{2}O interactions, would incorporate the 3-body H_{2}–H_{2}O–H_{2}O interactions, such as those computed by Bowman and co-workers,^{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 $Vh(1b)(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} the frequency shift measured at 20 K is about −37 cm^{−1},^{26} compared to −34 cm^{−1} at 76 K.^{27} Since the quantum 6D results are for 0 K, it is more appropriate to compare our computed value for H_{2} (−43 cm^{−1}) to the experimental shift at 20 K (−37 cm^{−1}). In that case, the measured shift is only 14% smaller in magnitude than the theoretical result.

Since the computed fully coupled quantum 6D vibrational frequency shifts are essentially numerically exact for the PES employed, the residual discrepancies between theory and experiment can be attributed primarily to the deficiencies in the description of the H_{2}-hydrate interaction. These can stem from the inaccuracy of the *ab initio* 6D H_{2}–H_{2}O pair potential and, more likely, the lack of three-body terms in the H_{2}-hydrate PES. The recent DMC calculations of Qu and Bowman^{41} of the vibrational frequency shift of H_{2} encapsulated in the (rigid) small cage of the sII hydrate surrounded by additional water molecules did include both the 2-body H_{2}–H_{2}O and 3-body H_{2}–H_{2}O–H_{2}O interactions. When only the 2-body interactions were considered, the frequency shift of about 26 cm^{−1} was obtained. Including the 3-body interactions yielded the shift of −40 ± 4 cm^{−1}, in good agreement with the measured shift of 37 cm^{−1} at 20 K.^{26}

## 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 gas-phase 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 gas-phase 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–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 ground-state 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–66} as is the wealth of spectroscopic data about its intra- and intermolecular excitations.

## ACKNOWLEDGMENTS

Z.B. is grateful to the National Science Foundation for its partial support of this research through Grant No. CHE-1566085. D.M.B. acknowledges the Viper High Performance Computing facility of the University of Hull and its support team.

## REFERENCES

*ElVibRot*: Quantum Dynamics Code, http://pagesperso.lcp.u-psud.fr/lauvergnat/ElVibRot/ElVibRot.html,