We report a theoretical study of the frequency shift (redshift) of the stretching fundamental transition of an H_{2} molecule confined inside the small dodecahedral cage of the structure II clathrate hydrate and its dependence on the condensed-phase environment. In order to determine how much the hydrate water molecules beyond the confining small cage contribute to the vibrational frequency shift, quantum five-dimensional (5D) calculations of the coupled translation-rotation eigenstates are performed for H_{2} in the $v=0$ and $v=1$ vibrational states inside spherical clathrate hydrate domains of increasing radius and a growing number of water molecules, ranging from 20 for the isolated small cage to over 1900. In these calculations, both H_{2} and the water domains are treated as rigid. The 5D intermolecular potential energy surface (PES) of H_{2} inside a hydrate domain is assumed to be pairwise additive. The H_{2}–H_{2}O pair interaction, represented by the 5D (rigid monomer) PES that depends on the vibrational state of H_{2}, $v=0$ or $v=1$, is derived from the high-quality *ab initio* full-dimensional (9D) PES of the H_{2}–H_{2}O complex [P. Valiron *et al.*, J. Chem. Phys. **129**, 134306 (2008)]. The H_{2} vibrational frequency shift calculated for the largest clathrate domain considered, which mimics the condensed-phase environment, is about 10% larger in magnitude than that obtained by taking into account only the small cage. The calculated splittings of the translational fundamental of H_{2} change very little with the domain size, unlike the H_{2} *j* = 1 rotational splittings that decrease significantly as the domain size increases. The changes in both the vibrational frequency shift and the *j* = 1 rotational splitting due to the condensed-phase effects arise predominantly from the H_{2}O molecules in the first three complete hydration shells around H_{2}.

## I. INTRODUCTION

Hydrogen clathrate hydrates are inclusion compounds, in which one or more hydrogen molecules are encapsulated inside closely packed polyhedral cavities within the crystalline three-dimensional (3D) framework of hydrogen-bonded water molecules.^{1–3} Simple hydrogen clathrate hydrates, having only hydrogen molecules as guests, were first identified by Dyadin *et al.*^{4} and subsequently studied in detail by Mao *et al.*^{5} They adopt the classical structure II (sII);^{1,2,5} its cubic unit cell has sixteen small cages, each consisting of 20 H_{2}O molecules and denoted as 5^{12} due to their 12 pentagonal faces, as well as eight large cages, each formed by 28 H_{2}O molecules and denoted as 5^{12}6^{4}, reflecting the presence of 12 pentagonal and 4 hexagonal faces. It has been established 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}

There has been a great deal of interest in hydrogen clathrate hydrates in recent years, motivated by their potential as efficient and environmentally friendly materials for hydrogen storage.^{1,2,7–10} In addition, hydrogen hydrates provide a rare opportunity for combined high-level experimental and theoretical investigations of the quantum dynamical effects and spectroscopic signatures associated with the entrapment of light molecules inside the large and small hydrate cages. This nanoscale confinement results in the quantization of the translational center-of-mass (c.m.) degrees of freedom of the guest molecule(s), which are coupled by the confining potential to their quantized rotational degrees of freedom. For a single hydrogen molecule in the cages of the sII clathrate hydrate, the salient features of its translation-rotation (TR) eigenstates have been revealed by the quantum 5D calculations.^{11–13} These include the splittings of both the translational fundamental and the *j* = 1 and *j* = 2 rotational levels of H_{2} caused to the anisotropies of the cage environment, the former with respect to the translational motion of the c.m. of H_{2}, and the latter with respect to its angular orientation within the cage. These features have been observed in the inelastic neutron scattering (INS)^{14,15} and the Raman spectra^{16–18} of the binary tetrahydrofuran (THF) + H_{2}/HD/D_{2} sII clathrate hydrate. The quantum TR dynamics of up to four H_{2} molecules in the large hydrate cage have been studied at *T* = 0 K by means of the diffusion Monte Carlo (DMC) calculations^{19} and at elevated temperatures (*T* = 25–200 K) using path-integral molecular dynamics (PIMD) simulations.^{20} In addition, fully quantal calculations of the TR eigenstates have been performed for two^{21,22} and four^{23} H_{2} molecules inside the large clathrate hydrate cage. Recently, diffusion was studied in the structure II clathrate hydrate for both H_{2} and D_{2}. The study concluded that quantum effects play an important role in the free-energy profiles for H_{2}.^{24} This result further indicates the significant role of quantum effects in this system.

Experimentally, the TR dynamics of hydrogen molecules in nanoporous materials can be probed most selectively utilizing the inelastic neutron scattering (INS) spectroscopy. However, reliable interpretation and assignment of such INS spectra is not possible without a theory that is capable of simulating them accurately. With this in mind, we have developed the quantum methodology enabling rigorous calculation of the INS spectra of a hydrogen molecule inside a nanocavity of arbitrary shape.^{25–27} With this novel approach, we have computed the INS spectra of H_{2} and HD in two binary D_{2}O clathrates: one with the cubic sII structure^{25,28,29} and the other having the hexagonal structure (sH).^{29,30} In both, a large promoter molecule resides in the large cages making them unavailable to H_{2}, while a single H_{2}/HD occupies the small cages, and in the case of the sH clathrate also the medium cages. Semiquantitative agreement was obtained between the simulated and experimental INS spectra, allowing the assignment of the latter in terms of various TR excitations of the guest molecule.

Another spectroscopic manifestation of the encapsulation of H_{2} molecules inside the cages of clathrate hydrates is the shift in the frequency of the H_{2} intramolecular stretching vibration away from that in the gas phase. This vibrational frequency shift is caused in part by the dependence of the H_{2}–hydrate interactions on the vibrational state, ground or excited, of H_{2}, that modifies the gap between the initial and final states of the spectroscopic transitions. In addition, vibrational frequency shifts reflect and are a sensitive probe of the TR dynamics of the confined chromophore, which in turn is affected by the geometry and chemical nature of the confining environment, presence of other guest molecules, temperature, and pressure. Vibrational frequencies of the confined H_{2} molecules have been measured primarily by Raman spectroscopy of the binary 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,16,17} The vibrations of the guest H_{2} molecules have always been found to have lower frequencies, i.e., redshifted, relative to the gas phase. The redshift is the largest, −34 cm^{−1}, for H_{2}, in the small cage.^{17} Raman spectroscopy measurements of simple hydrogen hydrates identified H_{2} vibrations with frequency shifts of −34 cm^{−1} attributed to the single H_{2} occupation of the small cage, while the frequency shifts of −26, −18, and −11 cm^{−1} have been associated with the large cavities occupied by two, three, and four H_{2} molecules, respectively.^{17}

For hydrogen hydrates, it has been possible to uniquely assign the observed frequency shifts to different H_{2} occupancies of the small and large clathrate cages, through very careful experiments that involved multiple cycles of heating and quenching of the sample and measurements of the released H_{2}.^{17} But even for this system, the redshifts from earlier Raman spectroscopic measurements^{16} were mis-assigned with respect to the cage occupancies by H_{2}.^{17} In general, extracting the information that vibrational frequency shifts contain regarding the cage occupancies and other structural, chemical, and dynamical aspects of a nanoconfined system requires theoretical methods capable of reliably calculating the frequency shifts. Accomplishing this task for H_{2} molecules is a challenging task for several reasons: (i) It is necessary to have an accurate potential energy surface (PES) describing the H_{2}-clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H_{2} molecules as well. (ii) The PES has to include the dependence on the H_{2} intermolecular stretch coordinate and its coupling to other intermolecular degrees of freedom. (iii) Due to the light mass of H_{2} and the low temperatures at which the Raman spectroscopy measurements are typically performed, dynamical quantum effects play a significant role and must be included in the treatment. (iv) For a meaningful comparison with the experimental H_{2} vibrational frequency shifts, it is not sufficient to consider only the nearest water molecules of the confining small or large cages, but also the framework water molecules further away.

Several calculations of the vibrational frequencies and frequency shifts of H_{2} molecules in sII clathrate hydrates have been reported in recent years. In some of them,^{31–33} one or more H_{2} molecules encapsulated in the isolated small or large hydrate cages are treated as static, frozen in the geometry corresponding to the minimum energy of the system. One of the shortcomings of this approach is that it completely leaves out nuclear quantum effects, in particular, the averaging over the large-amplitude intermolecular vibrations of the guest H_{2} molecule(s), which is important at low temperatures. In addition, since only isolated clathrate cages are considered, the effects of the condensed-matter environment are entirely neglected. Plattner and Meuwly^{34} have calculated the H_{2} vibrational frequency shifts in clathrate hydrates by combining classical molecular dynamics (MD) and path integral MD simulations with electronic structure calculations at the density functional theory (DFT) (B3LYP) and MP2 levels, for a system of 2 × 2 × 2 unit cells. H_{2} vibrational frequencies were calculated for the H_{2} coordinates taken from many snapshots of the MD simulations, resulting in a broad range of frequencies that extended to that of the free H_{2} at 4155 cm^{−1}. The maxima of these frequency distributions are in a reasonable agreement with the corresponding experimental values. However, it should be noted that since both the DFT and MP2 methods gave a too high vibrational frequency for free H_{2}, a scaling factor had to be introduced in all the calculations. Finally, Futera *et al.*^{35} have performed classical MD simulations within the DFT framework for an sII hydrate unit cell. The H_{2} vibrational spectra were calculated by Fourier transforming the H–H bond length autocorrelation function. This treatment does not account for the quantum effects. Moreover, the computed vibrational spectra are shifted by 100-150 cm^{−1} to higher frequencies relative to the experimental results, and above the stretch fundamental of free H_{2}.

In this paper, we present the results of the quantum 5D bound-state calculations of the vibrational frequency shift, and the TR eigenstates, of the H_{2} molecule in the small dodecahedral cage of the sII clathrate hydrate. The impact of the condensed phase is explored by performing these calculations for H_{2} inside several spherical clathrate domains of increasing radius and a growing number of H_{2}O molecules in them, from 20 in the isolated small cage to more than 1900 in the largest domain considered, and observing how (a) the frequency shift of H_{2} and (b) the splittings of translational fundamental and the rotational *j* = 1 triplet vary with the domain size. The same approach was used by us recently to elucidate the impact of the condensed-matter environment on the TR dynamics and INS spectra of H_{2} in the sII hydrate.^{36} In the quantum calculations, both H_{2} and the condensed-phase domains are treated as rigid. The 5D intermolecular PES for H_{2} in a clathrate domain is constructed in a pairwise additive fashion, utilizing an *ab initio* calculated 5D (rigid-monomer) H_{2}–H_{2}O pair potential. This pair potential depends on the vibrational state of H_{2}, $v=0$ or $v=1$, thus enabling the calculation of the H_{2} vibrational frequency shifts.

## II. METHODOLOGY

Accurate calculation of vibrational frequency shifts requires taking into account quantum effects, in particular, the large amplitude motion (LAM) vibrations present in the system considered. In systems where the chromophore is a light molecule, weakly bound to or encapsulated within a much heavier host, the most important quantum dynamical contributions to the vibrational frequency shift of the chromophore come from the intermolecular LAM vibrations of the guest molecule. This was clearly demonstrated by our earlier theoretical treatment of the vibrational frequency shifts of the HF stretch in various isomers of Ar_{n}HF clusters.^{37–39} In these calculations, the heavy Ar_{n} subunit was treated as rigid (i.e., its vibrational modes were frozen), while the quantum dynamics of the five remaining coupled LAM vibrations of (rigid) HF relative to Ar_{n} were treated with high accuracy. The pairwise additive 5D intermolecular PES employed included the dependence on the HF vibrational quantum numbers $v=0$ and $v=1$. For Ar_{n}HF, the HF vibrational frequency shift was obtained as the difference between the ground state energies $Ev=1$ and $Ev=0$ from two separate quantum 5D bound state calculations of the intermolecular vibrational levels, on the intermolecular PESs for HF $v=1$ and $v=0$, respectively. The frequency shifts calculated in this way agreed much better with the experimental results available at that time for *n* = 1–4 clusters^{40} than those that we obtained previously in a static approximation, with HF fixed at the equilibrium geometry.^{41} This vividly illustrated the vital importance of including the LAM vibrations of the chromophore in the frequency shift calculations. In fact, the high accuracy of our calculations of Ar_{n}HF vibrational frequency shifts^{39,42} enabled Nauta and Miller to assign their measured frequency shifts to various isomers of Ar_{n}HF clusters.^{43}

Adopting the approach outlined above, the vibrational frequency shift Δ*ν*_{n} for H_{2} inside the clathrate hydrate domain containing *n* water molecules can be computed as

where $En,iv$ is the ground translation-rotation state energy (*i* = 1) obtained from quantum 5D bound state calculations of the intermolecular vibrational energy levels of (rigid) H_{2} in the clathrate domain, on pairwise-additive PESs defined below, which depend on the vibrational states of H_{2}, $v\u2009=\u20091$ and $v\u2009=\u20090$, respectively. These calculations are five-dimensional since the hydrogen-bonded framework of water molecules encapsulating H_{2} is taken to be rigid. A large body of quantum calculations of the TR dynamics and INS spectra of H_{2} molecules in clathrate hydrates^{11–13,20,26,28,36} has shown that this assumption does not introduce significant error in the results.

The TR dynamics of H_{2} inside a clathrate domain are described in terms of the five coordinates **q** ≡ {*x*, *y*, *z*, *θ*, *φ*}; *x*, *y*, and *z* are the coordinates of the c.m. of H_{2} in the Cartesian frame fixed in the cage, while the two polar angles *θ* and *φ* specify its orientation relative to this frame. Since the overall system is a crystalline solid, the clathrate domains considered can be safely assumed to be infinitely heavy (in comparison with H_{2}) and non-rotating. Then the 5D TR Hamiltonian for the confined H_{2} molecule is^{11}

where $rv$ is the H_{2} bond length in the vibrational state $v$ (0 or 1), *m*_{H} is the atomic mass of hydrogen (1.008 amu), and the angular momentum operator $L^2$ of the diatomic is expanded in terms of the conjugate momenta of *θ* and *φ*. $VH2\u2212domain(v)(q)$ in Eq. (2) is the 5D intermolecular PES, defined below, for the interaction of H_{2} ($v\u2009=\u20090$ or $v\u2009=\u20091$) with the H_{2}O molecules within the water domain considered. The H_{2} bond lengths for $v\u2009=\u20090$, 1 are *r*_{0} = 1.419 bohrs and *r*_{1} = 1.457 bohrs, respectively, corresponding to the rotational constants *B*_{0} = 59.322 cm^{−1} and *B*_{1} = 56.260 cm^{−1} of the free H_{2} in $v=0$ and $v=1$ vibrational states.

The Smolyak sparse-grid technique^{44} implemented in ElVibRot^{45} is utilized to compute the 5D TR eigenstates of the Hamiltonian in Eq. (2). This method avoids the need for constructing a direct-product basis and grid and has been recently investigated by Avila and Carrington^{46–48} for the calculation of vibrational energy levels of semi-rigid molecules. Later, the same approach was also proposed by Lauvergnat and Nauts^{49–51} and has been used to compute the torsional levels of methanol in full dimensionality (12D).^{50} In this work, we present an application of this technique to the calculation of TR eigenstates. In this case, the basis functions are judiciously selected, similarly to the pruned basis set scheme of Avila and Carrington,^{48} and the single large direct-product grid is replaced by a sum of small direct-product grids. A major advantage of the Smolyak sparse grids is the reduction of the number of grid points with respect to a direct product approach. While useful in the present 5D calculations, this feature gains in importance for systems with six or more degrees of freedom. The primitive grids associated with the harmonic oscillator bases for the *x*, *y* and *z* coordinates, as well as the spherical harmonics for the coordinates *θ* and *φ*, are the Gauss-Hermite quadrature and the Lebedev grid points, respectively. For each translational coordinate, the largest numbers of basis functions and grid points are 21 and 24, respectively. For the spherical harmonics, the largest *j*_{max} is 10 and the largest number of Lebedev grid points is 170. These parameters give rise to the basis set with dimension 12 053 (for the TR levels of both *ortho*- and *para*-H_{2} in a single calculation) and 282 702 grid points; the TR energy levels from the quantum 5D calculations are converged to within 0.01 cm^{−1}. However, calculations with a smaller basis set (4065 basis functions) and grid (90 570 grid points) give TR energy levels with errors smaller than 0.1 cm^{−1} with respect to converged energy levels. For comparison, in a direct-product approach (basis set and grid) with the previous parameters, the total numbers of basis functions and grid points would be 1 120 581 and 23 500 080, respectively. Therefore, our approach gives a basis set of about 93 times smaller than the direct-product one and a Smolyak grid of 8 times smaller than the direct-product one.

The clathrate hydrate domains utilized in these calculations are constructed by the procedure described in Ref. 36. The starting point is the 3 × 3 × 3 supercell of the sII hydrate, having 27 primitive unit cells of the crystal. The O atoms of the water molecules in the supercell are placed at their crystallographic positions.^{52} Since the H atoms are configurationally disordered, different random distributions of the framework water protons are generated, each consistent with the Bernal-Fowler ice rules.^{53} One proton configuration with a negligible dipole moment is selected for use in all further calculations. Next, the energy of the empty supercell with the chosen proton distribution is minimized using the lattice dynamics program GULP,^{54} by allowing the (essentially rigid) H_{2}O molecules to change their orientation slightly, while keeping the O atoms fixed at their crystallographic positions.

In the next step, we define several spherical clathrate hydrate domains of increasing radius, with a growing number of H_{2}O molecules, around the small dodecahedral cage at the center of the supercell in the following manner.^{36} Based on the radial distribution function of the O atoms, several radial distances from the center of the central small cage are chosen so as to enclose successive complete well-defined shells of water molecules: (a) The spherical domain for the cutoff radius set to 5.0 Å consists of the 20 water molecules of the small cage itself. (b) When the radius is increased to 7.5 Å, additional 20 water molecules are included that are directly hydrogen-bonded to the small cage, generating a spherical domain with the total of 40 water molecules. (c) For the cutoff radius of 9.0 Å, yet another shell of 36 water molecules is enclosed, giving rise to a sphere containing 76 water molecules. (d) The largest spherical domain possible for our supercell has a radius of 25.65 Å, with a total of 1945 water molecules in it. As in Ref. 36, this domain constitutes our condensed phase environment for the small cage enclosing the H_{2} molecule. The structures of the four domains (a)–(d) above are displayed in Fig. 1 of Ref. 36. The smallest domain (a), i.e., the small cage itself, is also shown in Fig. 1, while the structure of the domain (c) with 76 water molecule is displayed in Fig. 2.

As stated earlier, the 5D intermolecular PES for the interaction between the encapsulated H_{2} in the vibrational states $v=0$ and $v=1$ and the water molecules within the clathrate domain considered, denoted as $VH2\u2212domain(v)$, is assumed to be pairwise additive. For *n* H_{2}O molecules within the domain (with *n* ranging from 20 to 1945), $VH2\u2212domain(v)$ can be written as

where **q** ≡ {*x*, *y*, *z*, *θ*, *φ*} are the coordinates of H_{2}, $VH2\u2212H2O(v)$ is the 5D pair interaction between H_{2}$v=0$, 1 and a framework H_{2}O molecule, and the index $w$ runs over the water molecules of the domain, whose coordinates $\Xi w$ are fixed.

The 5D (rigid monomer) two-body term $VH2\u2212H2O(v)$ in Eq. (3) is obtained from the recent accurate full-dimensional *ab initio* PES of H_{2}–H_{2}O by Valiron *et al.*^{55} This PES depends explicitly on the nine internal coordinates of the dimer, five of which are the intermolecular coordinates and four are the intramolecular coordinates for the distortion of the two monomers (see also Ref. 56). The 9D PES was obtained by combining standard coupled-cluster with single, double and perturbative contributions of the triple excitations [CCSD(T)] calculations with explicitly correlated CCSD(T)-R12 calculations, and it is expected to provide the most accurate description of the H_{2}O–H_{2} interaction currently available, to within a few wave numbers in the region of van der Waals minimum.^{55} The 5D (rigid monomer) PES $VH2\u2212H2O(v)$ ($v=0$, 1) is obtained by averaging the 9D PES over the vibrational ground state wave function of H_{2}O by Tennyson,^{57} and the vibrational wave functions of H_{2} for $v=0$ and $v=1$, respectively. The latter were computed using the Fourier Grid Hamiltonian (FGH) method^{58,59} on the H_{2} potential of Schwartz and Le Roy.^{60} Once the wave function averaging has been performed, $VH2\u2212H2O(v)(v=0,1)$ is numerically converted to atomic Cartesian coordinates for use in the quantum 5D bound-state calculations.

In order to verify the correctness of the bound state code using the TR Hamiltonian in Eq. (2) and the Smolyak sparse-grid scheme, we have also performed a rigid-body DMC (RB-DMC) calculation of the TR ground state of H_{2} ($v\u2009=\u20090$) in the small dodecahedral water cage. The simulations were performed with the Xdmc code developed by Benoit^{61} (see also Ref. 62 for implementation details). In this study, we used 2000 replicas, a stabilization period of 13 × 100 cycles, with Δ*τ* = 40 a.u., and an averaging phase of 200 × 100 cycles with Δ*τ* = 15 a.u. The converged value obtained is −742.5 ± 0.7 cm^{−1} and should be compared to the value −741.89 cm^{−1} computed using our Smolyak sparse-grid scheme. The very good agreement between the two values within the error bars confirms the accuracy of the Smolyak sparse-grid scheme employed. The 3D H_{2} c.m. probability distribution in the small dodecahedral water cage from the RB-DMC calculations is shown in Fig. 1.

## III. RESULTS AND DISCUSSIONS

The focus of this study is on the H_{2} vibrational frequency shifts and their evolution with the increasing size of the sII clathrate hydrate domain. Nevertheless, for completeness, in Table I, we give the energies of the fundamental translational excitations and the rotational *j* = 0 → 1 transitions of H_{2}($v=0$) in the four sII clathrate hydrate domains considered, where *N* ranges from 20 (small cage) to 1945, calculated with our Smolyak sparse-grid scheme. They are compared with the data from the INS spectroscopic measurements of the binary H_{2}–THF clathrate hydrate.^{14}

. | . | n
. | |||
---|---|---|---|---|---|

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

$En,10$ . | … . | −741.89 . | −778.56 . | −798.32 . | −819.05 . |

Translation | |||||

I | 71.0 | 66.78 | 66.36 | 66.13 | 66.03 |

II | 80.2 | 76.02 | 75.35 | 75.12 | 74.94 |

III | 101.1 | 93.13 | 92.34 | 92.06 | 92.07 |

Rotation j = 1 | |||||

I | 110.0 | 86.61 | 94.68 | 98.95 | 100.63 |

II | 116.5 | 122.18 | 122.17 | 119.40 | 119.35 |

III | 122.1 | 148.65 | 141.39 | 138.85 | 136.75 |

. | . | n
. | |||
---|---|---|---|---|---|

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

$En,10$ . | … . | −741.89 . | −778.56 . | −798.32 . | −819.05 . |

Translation | |||||

I | 71.0 | 66.78 | 66.36 | 66.13 | 66.03 |

II | 80.2 | 76.02 | 75.35 | 75.12 | 74.94 |

III | 101.1 | 93.13 | 92.34 | 92.06 | 92.07 |

Rotation j = 1 | |||||

I | 110.0 | 86.61 | 94.68 | 98.95 | 100.63 |

II | 116.5 | 122.18 | 122.17 | 119.40 | 119.35 |

III | 122.1 | 148.65 | 141.39 | 138.85 | 136.75 |

As mentioned in the Introduction, both the translational fundamental and *j* = 0 → 1 transition of H_{2} in the small sII cage are split into three components each (denoted as I, II, and III), by the anisotropies, radial and angular, respectively, of the cage environment.^{11,12} For *n* = 20, the calculated overall splitting of the translational fundamental (the energy difference between the components I and III), 26.4 cm^{−1}, agrees well with the measured splitting of 30.1 cm^{−1}. This shows that the pairwise-additive PES used in this work describes accurately the radial anisotropy of the environment of H_{2} in the small cage. However, the computed energies of the individual components I-III are 3-9 cm^{−1} lower than the corresponding experimental values. We note that the calculated energies of the translational fundamentals I-III hardly change with the size of the clathrate domain. This is consistent with the findings of the study in Ref. 36, where the pairwise-additive H_{2}-clathrate PES was constructed using a semi-empirical H_{2}–H_{2}O pair potential,^{63} that the water molecules beyond the confining small cage contribute very little to the radial anisotropy responsible for the splitting of the translational fundamental.

The situation is different in the case of the rotational *j* = 0 → 1 transition. As Table I shows, the splitting of the *j* = 1 triplet (the energy difference between the *j* = 1 components I and III) calculated for *N* = 20, 62.1 cm^{−1}, is much larger than the experimental result of 12.1 cm^{−1}. The calculated *j* = 1 splitting decreases significantly with increasing *N*, unlike the splitting of the translational fundamental. Most of this decrease is caused by the 76 water molecules comprising the first three complete solvation shells around H_{2}. For *n* = 76, the calculated *j* = 1 splitting of 39.9 cm^{−1} is only 3.8 cm^{−1}(9.5%) greater than the value of 36.12 cm^{−1} computed for *n* = 1945. The same conclusion that the first three hydration shells around H_{2} capture the main effects of the condensed phase on the *j* = 1 splittings in the small sII hydrate cage was reached earlier in Ref. 36.

But even for *n* = 1945, the *j* = 1 splitting of 36.12 cm^{−1} differs by about a factor of three from the experiment. Clearly, the pairwise-additive PES employed, based on the *ab initio* 5D H_{2}–H_{2}O pair potential, overestimates the angular anisotropy of the H_{2}-clathrate interaction, and this is not compensated by enlarging the clathrate domain included in the calculations. It is unlikely that the fault for this deficiency lies with this particular *ab initio* two-body H_{2}–H_{2}O interaction. The TR levels of H_{2} in the small dodecahedral cage calculated by us earlier,^{11} also on the pairwise additive 5D intermolecular PES constructed using another high-quality *ab initio* 5D (rigid-monomer) PES for H_{2}–H_{2}O,^{65} exhibited the same problems, in particular, the *j* = 1 splitting (57 cm^{−1}) much too large in comparison with the experiment. We have attributed these discrepancies to the nonadditive many-body interactions,^{12} which are missing from both pairwise-additive PESs. This view is supported by the results of our earlier calculations of the TR eigenstates of H_{2} in the small cage of the sII hydrate,^{12} using the H_{2}–H_{2}O pair interaction by Alavi *et al.*^{63} to construct the pairwise additive H_{2}-cage PES. This H_{2}–H_{2}O interaction is based on the empirical SPC/E H_{2}O–H_{2}O effective pair potential,^{64} which incorporates the many-body induced polarization implicitly, in an average mean-field sense. As a result, this H_{2}–H_{2}O pair interaction includes at least approximately the effects of the nonadditive water-water interactions on the H_{2}-water interaction potential. The *j* = 1 splitting calculated for this H_{2}-cage PES,^{12} 19.1 cm^{−1}, agrees much better with the experimental result of 12.1 cm^{−1} than the splittings computed for the two H_{2}-cage PESs constructed from the *ab initio* two-body H_{2}–H_{2}O potentials, thus lacking any many-body contributions.

The ground-state energies $En,11$ and $En,10$ for H_{2} in the vibrational states $v\u2009=\u20091$ and $v\u2009=\u20090$, respectively, from the quantum 5D calculations for four domains of increasing size are shown in Fig. 3. Both exhibit virtually the same decrease with the growing number of H_{2}O molecules. According to Eq. (1), the difference between $En,11$ and $En,10$ gives the computed vibrational frequency shift Δ*ν*_{n} for the domain with *n* water molecules.

The H_{2} vibrational frequency shifts for the domains of increasing size obtained in this way are given in Table II. The magnitude of the frequency shift increases appreciably with the domain size up to *n* = 76 and then continues to increase very slowly until *n* = 1945. This implies that the dominant contribution to the frequency shift comes from the 76 water molecules that form the first three complete hydration shells around H_{2}. In going from *n* = 76 to *n* = 1945, i.e., adding nearly 1900 H_{2}O molecules, the frequency shift increases in magnitude by just ∼3%. In addition, Table I shows that the magnitude of the frequency shift calculated for *n* = 1945 is only 9.6% larger than that obtained for *n* = 20, the small cage encapsulating H_{2}. The implication is that at least in this case, the vibrational frequency shift calculated by considering just H_{2} in the isolated small cage provides a good estimate, to within ∼10%, of the frequency shift that would be computed for the condensed-phase environment.

. | . | n
. | |||
---|---|---|---|---|---|

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

$En,11$ | … | −781.71 | −819.85 | −840.61 | −862.67 |

Δν_{n} | −34 | −39.81 | −41.29 | −42.30 | −43.62 |

. | . | n
. | |||
---|---|---|---|---|---|

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

$En,11$ | … | −781.71 | −819.85 | −840.61 | −862.67 |

Δν_{n} | −34 | −39.81 | −41.29 | −42.30 | −43.62 |

Extrapolating from the result for *n* = 1945 in Table II, our computed vibrational frequency shift for H_{2} in the small cage of sII clathrate hydrate is −44 cm^{−1}. The corresponding experimental value^{17} of −34 cm^{−1} is about 23% smaller. This level of agreement between theory and experiment is satisfactory in our opinion, keeping in mind that our frequency shifts have been calculated from first principles, without any adjustable parameters. In a recent theoretical study, Plattner and Meuwly,^{34} using MD simulations at the B3LYP and MP2 levels of electronic structure theory, have obtained vibrational frequency shifts for H_{2} in the small sII hydrate cage of −28 cm^{−1} at the B3LYP level and −32 cm^{−1} at the MP2 level. Good agreement of these results with experiment should be viewed with some caution, as the authors used a scaling factor, actually two, to bring the B3LYP and MP2 vibrational frequencies for free H_{2} in agreement with the experimental value.^{34}

## IV. CONCLUSIONS

We have performed quantum 5D bound-state calculations of the TR eigenstates for H_{2} molecules in the vibrational states $v\u2009=\u20090$ and $v\u2009=\u20091$ inside the small dodecahedral cage of the sII clathrate hydrate, embedded in spherical hydrate domains of increasing size, with the number of H_{2}O molecules *n* in them ranging from 20, to 40, 76, and 1945. The primary goal of this study has been to elucidate from first principles the degree to which the hydrate water molecules beyond the confining small cage contribute to the H_{2} vibrational frequency shift (redshift), and how this influence varies with the number of water molecules included and their distance from the central small cage. For this purpose, a pairwise additive 5D intermolecular PES for rigid H_{2} inside a hydrate domain, also treated as rigid, was employed, which depends on the vibrational state of H_{2}, $v\u2009=\u20090$ or $v\u2009=\u20091$. It is constructed from the 5D (rigid monomer) pair potential for the interaction of H_{2} ($v\u2009=\u20090$, 1) with H_{2}O, derived from an accurate full-dimensional (9D) *ab initio* PES of H_{2}–H_{2}O by Valiron *et al.*^{55} The TR eigenstates were calculated with an efficient Smolyak sparse-grid technique that avoids the direct-product structure required in the approaches using the multidimensional discrete variable representation.^{66,67}

Based on these calculations, the magnitude of the H_{2} vibrational frequency shift increases appreciably with the growing domain size up to *n* = 76 and then much more slowly until *n* = 1945. 76 water molecules form the first three complete hydration shells around the H_{2} molecule. The contribution they make to the frequency shift is dominant since further increase in the number of water molecules from *n* = 76 to *n* = 1945 causes the frequency shift to increase by just ∼3%. The overall impact of the hydrate condensed phase environment on the computed H_{2} vibrational frequency shift is rather modest since its magnitude calculated for the largest domain with *n* = 1945 is only 9.6% larger than that obtained for *n* = 20, the water molecules forming the small cage encapsulating H_{2}.

Our calculations have allowed us also to investigate how the clathrate-induced splittings of the translational fundamental and the *j* = 1 rotational level of the guest H_{2} vary with the domain size. The splitting of the translational fundamental calculated for *n* = 20 is in good agreement with the experimental result, although the energies of the individual components are slightly lower than the measured values. The splitting remains virtually constant over the entire range of the domain sizes considered, consistent with the observation made in an earlier study,^{36} that the H_{2}O molecules beyond the first solvation shell (the small cage) make a very small contribution to the radial anisotropy that gives rise to the splitting of the translational fundamental.

By contrast, increasing the domain size leads to significant decrease in the *j* = 1 rotational splitting. Most of the decrease takes place between *n* = 20 and *n* = 76, i.e., the completion of the third hydration shell, in agreement with the observations made in Ref. 36. The *j* = 1 splitting calculated for *n* = 76 is only 9.5% larger than that for *n* = 1945. The latter is still about a factor of three greater than the experimental value, evidencing that the pairwise-additive PES employed in this study overestimates the actual angular anisotropy of the H_{2}-clathrate interaction.

Based on the domain-size dependence of the H_{2} vibrational frequency shift and the result obtained for the largest domain, our computed frequency shift for H_{2} in the small cage of the sII hydrate is ∼−44 cm^{−1}. The corresponding experimental result,^{17} −34 cm^{−1}, is ∼23% smaller. Given that our first-principles calculations involve no adjustable or scaling parameters, the level of agreement with experiment is good. Clearly, there is room for improvement. The PES employed is pairwise-additive and as such does not include nonadditive many-body interactions. A more sophisticated *ab initio* PES, such as the one by Homayoon *et al.*,^{68} may be considered for future investigations. Another aspect of the PES that needs to be examined, and possibly improved, is the coupling between the H_{2} intramolecular stretch vibration and the intermolecular degrees of freedom, which is directly responsible for the frequency shift.

A natural extension of this work is the development of the fully quantal methodology for the calculations the vibrational frequency shifts when multiple H_{2} molecules are confined inside nanocavities, such as the large clathrate hydrate cages. The obvious starting points are some of the existing methods capable of treating the quantum dynamics of more than one nanoconfined H_{2} molecule.^{20–23} Currently, these methods deal with diatomics in their ground vibrational states, but it should be possible to generalize them to include excited intermolecular stretching vibrations as well. Work on this is in progress in our laboratories.

## ACKNOWLEDGMENTS

Y.S. thanks Alexandre Faure for fruitful discussions. Z.B. is grateful to the National Science Foundation for its partial support of this research through the Grant No. CHE-1566085.

## REFERENCES

*Xdmc*, Quantum Diffusion Monte Carlo code,