We report a theoretical study of the frequency shift (redshift) of the stretching fundamental transition of an H2 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 H2 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 H2 and the water domains are treated as rigid. The 5D intermolecular potential energy surface (PES) of H2 inside a hydrate domain is assumed to be pairwise additive. The H2–H2O pair interaction, represented by the 5D (rigid monomer) PES that depends on the vibrational state of H2, v=0 or v=1, is derived from the high-quality ab initio full-dimensional (9D) PES of the H2–H2O complex [P. Valiron et al., J. Chem. Phys. 129, 134306 (2008)]. The H2 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 H2 change very little with the domain size, unlike the H2j = 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 H2O molecules in the first three complete hydration shells around H2.

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 H2O molecules and denoted as 512 due to their 12 pentagonal faces, as well as eight large cages, each formed by 28 H2O molecules and denoted as 51264, reflecting the presence of 12 pentagonal and 4 hexagonal faces. It has been established that the small cage can accommodate only one H2 molecule, while up to four H2 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 H2 caused to the anisotropies of the cage environment, the former with respect to the translational motion of the c.m. of H2, 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 spectra16–18 of the binary tetrahydrofuran (THF) + H2/HD/D2 sII clathrate hydrate. The quantum TR dynamics of up to four H2 molecules in the large hydrate cage have been studied at T = 0 K by means of the diffusion Monte Carlo (DMC) calculations19 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 two21,22 and four23 H2 molecules inside the large clathrate hydrate cage. Recently, diffusion was studied in the structure II clathrate hydrate for both H2 and D2. The study concluded that quantum effects play an important role in the free-energy profiles for H2.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 H2 and HD in two binary D2O clathrates: one with the cubic sII structure25,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 H2, while a single H2/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 H2 molecules inside the cages of clathrate hydrates is the shift in the frequency of the H2 intramolecular stretching vibration away from that in the gas phase. This vibrational frequency shift is caused in part by the dependence of the H2–hydrate interactions on the vibrational state, ground or excited, of H2, 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 H2 molecules have been measured primarily by Raman spectroscopy of the binary THF + H2 sII hydrate, where the large cages are completely occupied by the THF, while the small cages are singly occupied by H2, and simple sII hydrates in which H2 molecules are the only guests.10,16,17 The vibrations of the guest H2 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 H2, in the small cage.17 Raman spectroscopy measurements of simple hydrogen hydrates identified H2 vibrations with frequency shifts of −34 cm−1 attributed to the single H2 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 H2 molecules, respectively.17 

For hydrogen hydrates, it has been possible to uniquely assign the observed frequency shifts to different H2 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 H2.17 But even for this system, the redshifts from earlier Raman spectroscopic measurements16 were mis-assigned with respect to the cage occupancies by H2.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 H2 molecules is a challenging task for several reasons: (i) It is necessary to have an accurate potential energy surface (PES) describing the H2-clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H2 molecules as well. (ii) The PES has to include the dependence on the H2 intermolecular stretch coordinate and its coupling to other intermolecular degrees of freedom. (iii) Due to the light mass of H2 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 H2 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 H2 molecules in sII clathrate hydrates have been reported in recent years. In some of them,31–33 one or more H2 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 H2 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 Meuwly34 have calculated the H2 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. H2 vibrational frequencies were calculated for the H2 coordinates taken from many snapshots of the MD simulations, resulting in a broad range of frequencies that extended to that of the free H2 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 H2, 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 H2 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 H2.

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 H2 molecule in the small dodecahedral cage of the sII clathrate hydrate. The impact of the condensed phase is explored by performing these calculations for H2 inside several spherical clathrate domains of increasing radius and a growing number of H2O 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 H2 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 H2 in the sII hydrate.36 In the quantum calculations, both H2 and the condensed-phase domains are treated as rigid. The 5D intermolecular PES for H2 in a clathrate domain is constructed in a pairwise additive fashion, utilizing an ab initio calculated 5D (rigid-monomer) H2–H2O pair potential. This pair potential depends on the vibrational state of H2, v=0 or v=1, thus enabling the calculation of the H2 vibrational frequency shifts.

The 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 directions of future research.

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 ArnHF clusters.37–39 In these calculations, the heavy Arn 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 Arn 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 ArnHF, 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 clusters40 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 ArnHF vibrational frequency shifts39,42 enabled Nauta and Miller to assign their measured frequency shifts to various isomers of ArnHF clusters.43 

Adopting the approach outlined above, the vibrational frequency shift Δνn for H2 inside the clathrate hydrate domain containing n water molecules can be computed as

Δνn=En,11En,10,
(1)

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) H2 in the clathrate domain, on pairwise-additive PESs defined below, which depend on the vibrational states of H2, v=1 and v=0, respectively. These calculations are five-dimensional since the hydrogen-bonded framework of water molecules encapsulating H2 is taken to be rigid. A large body of quantum calculations of the TR dynamics and INS spectra of H2 molecules in clathrate hydrates11–13,20,26,28,36 has shown that this assumption does not introduce significant error in the results.

The TR dynamics of H2 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 H2 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 H2) and non-rotating. Then the 5D TR Hamiltonian for the confined H2 molecule is11 

Ĥ=2212mH2x2+2y2+2z2+2mHrv22φ2+1tanθθ+1sin2θ2θ2+VH2domain(v)(q),
(2)

where rv is the H2 bond length in the vibrational state v (0 or 1), mH 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 φ. VH2domain(v)(q) in Eq. (2) is the 5D intermolecular PES, defined below, for the interaction of H2 (v=0 or v=1) with the H2O molecules within the water domain considered. The H2 bond lengths for v=0, 1 are r0 = 1.419 bohrs and r1 = 1.457 bohrs, respectively, corresponding to the rotational constants B0 = 59.322 cm−1 and B1 = 56.260 cm−1 of the free H2 in v=0 and v=1 vibrational states.

The Smolyak sparse-grid technique44 implemented in ElVibRot45 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 Carrington46–48 for the calculation of vibrational energy levels of semi-rigid molecules. Later, the same approach was also proposed by Lauvergnat and Nauts49–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 jmax 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-H2 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) H2O 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 H2O 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 H2 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.

FIG. 1.

3D isosurface of the H2 c.m. probability distribution from the RB-DMC calculations of H2 (v = 0) in the ground translation-rotation state inside the small dodecahedral clathrate hydrate cage.

FIG. 1.

3D isosurface of the H2 c.m. probability distribution from the RB-DMC calculations of H2 (v = 0) in the ground translation-rotation state inside the small dodecahedral clathrate hydrate cage.

Close modal
FIG. 2.

The structure of the sII clathrate hydrate domain including the small dodecahedral cage (20 H2O molecules, red) and additional 56 H2O molecules, 20 (blue) + 36 (green), around it, with the total of 76 H2O molecules.

FIG. 2.

The structure of the sII clathrate hydrate domain including the small dodecahedral cage (20 H2O molecules, red) and additional 56 H2O molecules, 20 (blue) + 36 (green), around it, with the total of 76 H2O molecules.

Close modal

As stated earlier, the 5D intermolecular PES for the interaction between the encapsulated H2 in the vibrational states v=0 and v=1 and the water molecules within the clathrate domain considered, denoted as VH2domain(v), is assumed to be pairwise additive. For n H2O molecules within the domain (with n ranging from 20 to 1945), VH2domain(v) can be written as

VH2domain(v)(q)=w=1nVH2H2O(v)(q,Ξw),
(3)

where q ≡ {x, y, z, θ, φ} are the coordinates of H2, VH2H2O(v) is the 5D pair interaction between H2v=0, 1 and a framework H2O molecule, and the index w runs over the water molecules of the domain, whose coordinates Ξw are fixed.

The 5D (rigid monomer) two-body term VH2H2O(v) in Eq. (3) is obtained from the recent accurate full-dimensional ab initio PES of H2–H2O 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 H2O–H2 interaction currently available, to within a few wave numbers in the region of van der Waals minimum.55 The 5D (rigid monomer) PES VH2H2O(v) (v=0, 1) is obtained by averaging the 9D PES over the vibrational ground state wave function of H2O by Tennyson,57 and the vibrational wave functions of H2 for v=0 and v=1, respectively. The latter were computed using the Fourier Grid Hamiltonian (FGH) method58,59 on the H2 potential of Schwartz and Le Roy.60 Once the wave function averaging has been performed, VH2H2O(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 H2 (v=0) in the small dodecahedral water cage. The simulations were performed with the Xdmc code developed by Benoit61 (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 H2 c.m. probability distribution in the small dodecahedral water cage from the RB-DMC calculations is shown in Fig. 1.

The focus of this study is on the H2 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 H2(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 H2–THF clathrate hydrate.14 

TABLE I.

Comparison of the calculated energies (in cm−1) of the fundamental translational excitations and the rotational j = 0 → 1 transitions of H2 (v=0) in four sII clathrate hydrate domains with the experimental results from Ref. 14 (in boldface). For a domain with N water molecules, the calculated excitation energies are relative to the ground translation-rotation state energy En,1v=0 of this domain from the quantum 5D calculations.

n
Expt.2040761945
En,10−741.89−778.56−798.32−819.05
Translation 
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 
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.2040761945
En,10−741.89−778.56−798.32−819.05
Translation 
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 
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 H2 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 H2 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 H2-clathrate PES was constructed using a semi-empirical H2–H2O 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 H2. 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 H2 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 H2–H2O pair potential, overestimates the angular anisotropy of the H2-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 H2–H2O interaction. The TR levels of H2 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 H2–H2O,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 H2 in the small cage of the sII hydrate,12 using the H2–H2O pair interaction by Alavi et al.63 to construct the pairwise additive H2-cage PES. This H2–H2O interaction is based on the empirical SPC/E H2O–H2O effective pair potential,64 which incorporates the many-body induced polarization implicitly, in an average mean-field sense. As a result, this H2–H2O pair interaction includes at least approximately the effects of the nonadditive water-water interactions on the H2-water interaction potential. The j = 1 splitting calculated for this H2-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 H2-cage PESs constructed from the ab initio two-body H2–H2O potentials, thus lacking any many-body contributions.

The ground-state energies En,11 and En,10 for H2 in the vibrational states v=1 and v=0, 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 H2O 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.

FIG. 3.

The ground translation-rotation state energy En,1v for H2 in the vibrational states v=1 and v=0, respectively, from the quantum 5D calculations for the domains as a function of n−1, with n = 20, 40, 76, and 1945. The dashed line is a linear regression of the data points.

FIG. 3.

The ground translation-rotation state energy En,1v for H2 in the vibrational states v=1 and v=0, respectively, from the quantum 5D calculations for the domains as a function of n−1, with n = 20, 40, 76, and 1945. The dashed line is a linear regression of the data points.

Close modal

The H2 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 H2. In going from n = 76 to n = 1945, i.e., adding nearly 1900 H2O 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 H2. The implication is that at least in this case, the vibrational frequency shift calculated by considering just H2 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.

TABLE II.

H2 vibrational frequency shifts Δνn (in cm−1) from the quantum 5D calculations in this work, for clathrate domains having n water molecules. The experimental result for H2 in the small cage of the sII hydrate is from Ref. 17. En,11 (in cm−1) is the quantum 5D calculated ground-state energy for H2 (v=1) inside the domain of size n.

n
Expt.2040761945
En,11 … −781.71 −819.85 −840.61 −862.67 
Δνn −34 −39.81 −41.29 −42.30 −43.62 
n
Expt.2040761945
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 H2 in the small cage of sII clathrate hydrate is −44 cm−1. The corresponding experimental value17 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 H2 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 H2 in agreement with the experimental value.34 

We have performed quantum 5D bound-state calculations of the TR eigenstates for H2 molecules in the vibrational states v=0 and v=1 inside the small dodecahedral cage of the sII clathrate hydrate, embedded in spherical hydrate domains of increasing size, with the number of H2O 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 H2 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 H2 inside a hydrate domain, also treated as rigid, was employed, which depends on the vibrational state of H2, v=0 or v=1. It is constructed from the 5D (rigid monomer) pair potential for the interaction of H2 (v=0, 1) with H2O, derived from an accurate full-dimensional (9D) ab initio PES of H2–H2O 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 H2 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 H2 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 H2 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 H2.

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 H2 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 H2O 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 H2-clathrate interaction.

Based on the domain-size dependence of the H2 vibrational frequency shift and the result obtained for the largest domain, our computed frequency shift for H2 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 H2 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 H2 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 H2 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.

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.

1.
W. L.
Mao
,
C. A.
Koh
, and
E. D.
Sloan
,
Phys. Today
60
(
10
),
42
(
2007
).
2.
V. V.
Struzhkin
,
B.
Militzer
,
W. L.
Mao
,
H. K.
Mao
, and
R. J.
Hemley
,
Chem. Rev.
107
,
4133
(
2007
).
3.
E. D.
Sloan
,
Clathrate Hydrates of Natural Gases
(
Marcel Dekker
,
New York
,
1998
).
4.
Y. A.
Dyadin
,
E. G.
Larionov
,
A. Y.
Manakov
,
F. V.
Zhurko
,
E. Y.
Aladko
,
T. V.
Mikina
, and
V. Y.
Komarov
,
Mendeleev Commun.
9
,
209
(
1999
).
5.
W. L.
Mao
,
H. K.
Mao
,
A. F.
Goncharov
,
V. V.
Struzhkin
,
Q.
Guo
,
J.
Hu
,
J.
Shu
,
R. J.
Hemley
,
M.
Somayazulu
, and
Y.
Zhao
,
Science
297
,
2247
(
2002
).
6.
K. A.
Lokshin
,
Y.
Zhao
,
D.
He
,
W. L.
Mao
,
H. K.
Mao
,
R. J.
Hemley
,
M. V.
Lobanov
, and
M.
Greenblatt
,
Phys. Rev. Lett.
93
,
125503
(
2004
).
7.
W. L.
Mao
and
H. K.
Mao
,
Proc. Natl. Acad. Sci. U. S. A.
101
,
708
(
2004
).
9.
Y. H.
Hu
and
E.
Ruckenstein
,
Angew. Chem., Int. Ed.
45
,
2011
(
2006
).
10.
T. A.
Strobel
,
K. C.
Hester
,
C. A.
Koh
,
A. K.
Sum
, and
E. D.
Sloan
, Jr.
,
Chem. Phys. Lett.
478
,
97
(
2009
).
11.
M.
Xu
,
Y.
Elmatad
,
F.
Sebastianelli
,
J. W.
Moskowitz
, and
Z.
Bačić
,
J. Phys. Chem. B
110
,
24806
(
2006
).
12.
M.
Xu
,
F.
Sebastianelli
, and
Z.
Bačić
,
J. Chem. Phys.
128
,
244715
(
2008
).
13.
M.
Xu
,
F.
Sebastianelli
, and
Z.
Bačić
,
J. Phys. Chem. A
113
,
7601
(
2009
).
14.
L.
Ulivi
,
M.
Celli
,
A.
Giannasi
,
A. J.
Ramirez-Cuesta
,
D. J.
Bull
, and
M.
Zoppi
,
Phys. Rev. B
76
,
161401(R)
(
2007
).
15.
L.
Ulivi
,
M.
Celli
,
A.
Giannasi
,
A. J.
Ramirez-Cuesta
, and
M.
Zoppi
,
J. Phys.: Condens. Matter
20
,
104242
(
2008
).
16.
A.
Giannasi
,
M.
Celli
,
L.
Ulivi
, and
M.
Zoppi
,
J. Chem. Phys.
129
,
084705
(
2008
).
17.
T. A.
Strobel
,
E. D.
Sloan
, and
C. A.
Koh
,
J. Chem. Phys.
130
,
014506
(
2009
).
18.
A.
Giannasi
,
M.
Celli
,
M.
Zoppi
,
M.
Moraldi
, and
L.
Ulivi
,
J. Chem. Phys.
135
,
054506
(
2011
).
19.
F.
Sebastianelli
,
M.
Xu
, and
Z.
Bačić
,
J. Chem. Phys.
129
,
244706
(
2008
).
20.
A.
Witt
,
F.
Sebastianelli
,
M. E.
Tuckerman
, and
Z.
Bačić
,
J. Phys. Chem. C
114
,
20775
(
2010
).
21.
A.
Valdeś
and
G. J.
Kroes
,
J. Phys. Chem. C
116
,
21664
(
2012
).
22.
P. M.
Felker
,
J. Chem. Phys.
141
,
184305
(
2014
).
23.
P. M.
Felker
,
J. Chem. Phys.
138
,
174306
(
2013
).
24.
J. R.
Cendagorta
,
A.
Powers
,
T. J. H.
Hele
,
O.
Marsalek
,
Z.
Bačić
, and
M. E.
Tuckerman
,
Phys. Chem. Chem. Phys.
18
,
32169
(
2016
).
25.
M.
Xu
,
L.
Ulivi
,
M.
Celli
,
D.
Colognesi
, and
Z.
Bačić
,
Phys. Rev. B
83
,
241403(R)
(
2011
).
26.
M.
Xu
and
Z.
Bačić
,
Phys. Rev. B
84
,
195445
(
2011
).
27.
M.
Xu
,
L.
Ulivi
,
M.
Celli
,
D.
Colognesi
, and
Z.
Bačić
,
Chem. Phys. Lett.
563
,
1
(
2013
).
28.
D.
Colognesi
,
M.
Celli
,
L.
Ulivi
,
M.
Xu
, and
Z.
Bačić
,
J. Phys. Chem. A
117
,
7314
(
2013
).
29.
D.
Colognesi
,
A.
Powers
,
M.
Celli
,
M.
Xu
,
Z.
Bačić
, and
L.
Ulivi
,
J. Chem. Phys.
141
,
134501
(
2014
).
30.
M.
Celli
,
A.
Powers
,
D.
Colognesi
,
M.
Xu
,
Z.
Bačić
, and
L.
Ulivi
,
J. Chem. Phys.
139
,
164507
(
2013
).
31.
J.
Wang
,
H.
Lu
, and
J. A.
Ripmeester
,
J. Am. Chem. Soc.
131
,
14132
(
2010
).
32.
J.
Wang
,
H.
Lu
,
J. A.
Ripmeester
, and
U.
Becker
,
J. Phys. Chem. C
114
,
21042
(
2010
).
33.
K. R.
Ramya
and
A.
Venkatnathan
,
J. Chem. Phys.
138
,
124305
(
2013
).
34.
N.
Plattner
and
M.
Meuwly
,
J. Chem. Phys.
140
,
024311
(
2014
).
35.
Z.
Futera
,
M.
Celli
,
L.
del Rosso
,
C. J.
Burnham
,
L.
Ulivi
, and
N. J.
English
,
J. Phys. Chem. C
121
,
3690
(
2017
).
36.
A.
Powers
,
O.
Marsalek
,
M.
Xu
,
L.
Ulivi
,
D.
Colognesi
,
M. E.
Tuckerman
, and
Z.
Bačić
,
J. Phys. Chem. Lett.
7
,
308
(
2016
).
37.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
6359
(
1994
).
38.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
10181
(
1994
).
39.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
103
,
1829
(
1995
).
40.
A.
McIlroy
,
R.
Lascola
,
C. M.
Lovejoy
, and
D. J.
Nesbitt
,
J. Phys. Chem.
95
,
2636
(
1991
).
41.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
100
,
7166
(
1994
).
42.
J. M.
Hutson
,
S.
Liu
,
J. W.
Moskowitz
, and
Z.
Bačić
,
J. Chem. Phys.
111
,
8378
(
1999
).
43.
K.
Nauta
and
R. E.
Miller
,
J. Chem. Phys.
115
,
10138
(
2001
).
44.
S. A.
Smolyak
,
Sov. Math. Dokl.
4
,
240
(
1963
).
45.
D.
Lauvergnat
, ElVibRot: Quantum dynamics code http://pagesperso.lcp.u-psud.fr/lauvergnat/ElVibRot/ElVibRot.html, 2006.
46.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
131
,
174103
(
2009
).
47.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
134
,
054126
(
2011
).
48.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
135
,
064101
(
2011
).
49.
D.
Lauvergnat
and
A.
Nauts
,
Phys. Chem. Chem. Phys.
12
,
8405
(
2010
).
50.
D.
Lauvergnat
and
A.
Nauts
,
Spectrochim. Acta, Part A
119
,
18
(
2014
).
51.
J.
Sarka
,
D.
Lauvergnat
,
V.
Brites
,
A. G.
Császár
, and
C.
Léonard
,
Phys. Chem. Chem. Phys.
18
,
17678
(
2016
).
52.
C. Y.
Jones
,
S. L.
Marshall
,
B. C.
Chakoumakos
,
C. J.
Rawn
, and
Y.
Ishii
,
J. Phys. Chem. B
107
,
6026
(
2003
).
53.
J. D.
Bernal
and
R. H.
Fowler
,
J. Chem. Phys.
1
,
515
(
1933
).
54.
J. D.
Gale
,
J. Chem. Soc., Faraday Trans.
93
,
629
(
1997
).
55.
P.
Valiron
,
M.
Wernli
,
A.
Faure
,
L.
Wiesenfeld
,
C.
Rist
,
S.
Kedzuch
, and
J.
Noga
,
J. Chem. Phys.
129
,
134306
(
2008
).
56.
A.
Faure
,
L.
Wiesenfeld
,
M.
Wernli
, and
P.
Valiron
,
J. Chem. Phys.
123
,
104309
(
2005
).
57.
O. L.
Polyansky
,
A. G.
Császár
,
S. V.
Shirin
,
N. F.
Zobov
,
J.
Tennyson
,
P.
Barletta
,
D. W.
Schwenke
, and
P. J.
Knowles
,
Science
299
,
539
(
2003
).
58.
C. C.
Marston
and
G. G.
Balint-Kurti
,
J. Chem. Phys.
91
,
3571
(
1989
).
59.
G. G.
Balint-Kurti
and
C. L.
Ward
,
Comput. Phys. Commun.
67
,
285
(
1991
).
60.
C.
Schwartz
and
R. J.
LeRoy
,
J. Mol. Spectrosc.
121
,
420
(
1987
).
61.
D. M.
Benoit
, Xdmc, Quantum Diffusion Monte Carlo code,
1999
.
62.
D. M.
Benoit
and
D. C.
Clary
,
J. Chem. Phys.
113
,
5193
(
2000
).
63.
S.
Alavi
,
J. A.
Ripmeester
, and
D. D.
Klug
,
J. Chem. Phys.
123
,
024507
(
2005
).
64.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
65.
M. P.
Hodges
,
R. J.
Wheatley
,
G. K.
Schenter
, and
A. H.
Harvey
,
J. Chem. Phys.
120
,
710
(
2004
).
66.
Z.
Bačić
and
J. C.
Light
,
Annu. Rev. Phys. Chem.
40
,
469
(
1989
).
67.
J. C.
Light
and
T.
Carrington
, Jr.
Adv. Chem. Phys.
114
,
263
(
2000
).
68.
Z.
Homayoon
,
R.
Conte
,
C.
Qu
, and
J. M.
Bowman
,
J. Chem. Phys.
143
,
084302
(
2015
).