We report the results of calculations pertaining to the HH intramolecular stretching fundamentals of (p-H2)2 encapsulated in the large cage of structure II clathrate hydrate. The eight-dimensional (8D) quantum treatment assumes rotationless (j = 0) H2 moieties and a rigid clathrate structure but is otherwise fully coupled. The (H2)2-clathrate interaction is constructed in a pairwise-additive fashion, by combining the ab initio H2–H2O pair potential for flexible H2 and rigid H2O [D. Lauvergnat et al., J. Chem. Phys. 150, 154303 (2019)] and the six-dimensional (6D) H2–H2 potential energy surface [R. J. Hinde, J. Chem. Phys. 128, 154308 (2008)]. The calculations are performed by first solving for the eigenstates of a reduced-dimension 6D “intermolecular” Hamiltonian extracted from the full 8D Hamiltonian by taking the H2 moieties to be rigid. An 8D contracted product basis for the solution of the full problem is then constructed from a small number of the lowest-energy 6D intermolecular eigenstates and two discrete variable representations covering the H2-monomer internuclear distances. Converged results are obtained already by including just the two lowest intermolecular eigenstates in the final 8D basis of dimension 128. The two HH vibrational stretching fundamentals are computed for three hydrate domains having an increasing number of H2O molecules. For the largest domain, the two fundamentals are found to be site-split by ∼0.5 cm−1 and to be redshifted by about 24 cm−1 from the free-H2 monomer stretch frequency, in excellent agreement with the experimental value of 26 cm−1. A first-order perturbation theory treatment gives results that are nearly identical to those of the 8D quantum calculations.

Hydrogen clathrate hydrates are inclusion compounds in which one or more hydrogen molecules are confined inside closely packed polyhedral cavities within the three-dimensional (3D) crystalline framework created by hydrogen-bonded water molecules.1–3 Simple hydrogen clathrate hydrates, with only hydrogen molecules as guests, first identified by Dyadin et al.4 and subsequently studied in more detail by Mao et al.,5 adopt the classical structure II (sII).1,2,5 Its unit cell is cubic, comprised of two types of cages: (a) 16 small cages, each consisting of 20 H2O molecules and denoted 512 due to its 12 pentagonal faces, and (b) eight large cages, each formed by 28 H2O molecules arranged in 12 pentagonal and 4 hexagonal faces and therefore denoted as 51264. The small cage has been shown to accommodate only one H2 molecule, while up to four H2 molecules can be encapsulated in the large cage.6 Hydrogen clathrate hydrates have been the subject of a great deal of research in recent years, owing to their potential as hydrogen storage materials that could be both economical and environmentally friendly.1,2,7–10

The dynamics and spectroscopy of hydrogen molecules entrapped inside the cages of the clathrate hydrate are dominated by strong quantum effects to a degree seen only for light molecules inside fullerenes.11 These quantum effects have multiple sources. Nanoscale confinement in the hydrate cages gives rise to the quantized translational center-of-mass (c.m.) degrees of freedom (DOFs) of the guest molecule(s) (particle-in-a-box effect). They are coupled by the confining potential of the hydrate cage to the also quantized rotational DOFs of the hydrogen molecule(s). The translation-rotation (TR) energy level structure that results is sparse due to the low molecular mass of H2/HD/D2, their large rotational constants, and the small size of the hydrate cavities. For a single hydrogen molecule in the cages of the sII clathrate hydrate, the salient features of its TR eigenstates, 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 calculations12–15 and rigorous computations of the corresponding INS spectra.15–19 These features have also been observed experimentally in the INS20,21 and the Raman spectra22–24 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 has been studied at T = 0 K by means of the diffusion Monte Carlo (DMC) calculations25 and at elevated temperatures (T = 25–200 K) using path-integral molecular dynamics (PIMD) simulations.26 In addition, fully quantal calculations of the TR eigenstates have been performed for two27,28 and four29 H2 molecules inside the large clathrate hydrate cage. In all these calculations, both the H2 molecules and the hydrogen-bonded clathrate hydrate framework were treated as rigid. Recently,30 this constraint was relaxed partially, by performing quantum 5D calculations of the TR levels of (rigid) H2 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.

Another spectroscopic manifestation of the encapsulation of hydrogen molecules in the cages of clathrate hydrates, particularly relevant for this study, is the shift in the frequency of the H2 intramolecular stretching vibration away from that in the gas phase. It is readily observable in the Raman spectra of the binary tetrahydrofuran (THF) + H2 sII hydrate, where the large cages are completely occupied by the THF while the small cages are singly occupied by H2 and also in those of simple sII hydrates in which H2 molecules are the only guests.10,22,23 The vibrational frequencies of H2 molecules encapsulated in the sII hydrates are always lower than, i.e., redshifted relative to, the gas-phase H2. The largest redshift, −34 cm−1, is observed in the Raman spectra of the THF + H2 sII hydrate and can be assigned unambiguously to the singly H2 occupied small cage.10,22,23 The same redshift of −34 cm−1 appearing in the Raman spectra of the simple sII hydrate is therefore attributed to H2 in the small cage.

Also observed in the Raman spectra of the simple II hydrate are the bands redshifted by −26, −18, and −11 cm−1, respectively,10,22,23 that must represent contributions from the large cages whose H2 occupancy ranges between two and four. Thus, the frequency shifts are very sensitive to the number of H2 molecules confined in the cage. However, interpreting them in terms of a particular H2 occupancy of the large cages has turned out to be nontrivial, and the initial attempts22 proved to be incorrect. The subsequent elaborate experiments led to the assignment of these three redshifts to double, triple, and quadruple H2 occupancies of the large cages, respectively.23 

In the case of sII hydrogen hydrates, it was possible to assign the observed frequency shifts to different H2 occupancies of the small and large clathrate cages based on detailed experimental data alone. However, in general, e.g., molecular hydrogen in metal-organic frameworks (MOFs),31,32 extracting the information encoded in the measured vibrational frequency shifts regarding the H2 occupancies of the cavities of nanoporous materials, and other structural as well as dynamical aspects of the entrapped H2, demands theoretical methods capable of accurate calculation of the frequency shifts. This task is highly challenging for two reasons. The first is posed by the high dimensionality of the problem. Even if the host cavity is treated as rigid, the dimensionality of the calculations is 6nD, n being the number of encapsulated H2 molecules, when treated as flexible. 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 H2-clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H2 molecules. In addition, both interactions must include the dependence on the H2 intramolecular stretch coordinate and its coupling to the intermolecular DOFs. Second, dynamical quantum effects and anharmonicities in both intra- and intermolecular DOFs are significant, particularly at the low temperatures of the Raman spectroscopy measurements. Consequently, these key features have to be described correctly by any first-principles theoretical method whose goal is to determine accurate intramolecular vibrational frequency shifts of encapsulated hydrogen molecules.

A number of approaches, relying on different approximations, have been utilized to address this fundamental and difficult problem. In some of them, the H2 molecules encapsulated in the isolated small or large hydrate cages were treated as static, frozen in the geometry corresponding to the minimum energy of the system.33–35 This leaves out nuclear quantum effects, especially the averaging over the large-amplitude intermolecular vibrations of the guest H2 molecules. This problem was also approached by combining classical molecular dynamics (MD) and PIMD simulations with electronic structure calculations at the density functional theory (DFT) (B3LYP) and MP2 levels.36 The H2 vibrational frequencies calculated in 1D for the H2 intermolecular coordinates taken from many snapshots of the MD simulations covered a broad distribution of frequencies that extended to that of the free H2 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 H2 vibrational spectra were calculated by Fourier transforming the H–H bond length autocorrelation function.37 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 H2.

Recently, Powers et al.38 have calculated the frequency shift of H2 inside the small cage of the sII hydrate, isolated, and surrounded by spherical hydrate domains of increasing size, allowing the investigation of the effects of the condensed-phase environment. The approach employed was that developed earlier by Bačić and co-workers for the purpose of computing the HF stretch frequency shift in ArnHF clusters.39–42 The H2 frequency shift was obtained by means of the quantum 5D bound-state calculations of the coupled TR eigenstates on a pair of effective pairwise-additive intermolecular PESs for rigid H2 in a (rigid) hydrate domain that depend on the vibrational state of H2, v = 0, or v = 1, respectively. These 5D PESs were constructed using the 5D (rigid-monomer) pair potential for the interaction of H2 in the ground and first excited vibrational states, respectively, with H2O, obtained by averaging the full-dimensional (9D) ab initio PES of H2–H2O by Valiron et al.43 over the vibrational ground state wave function of H2O and the vibrational wave functions of H2 for v = 0 and v = 1, respectively. This approach rests on the assumption of dynamical decoupling between the H2 intramolecular vibration and the TR modes, well-justified by their large energy separation. The H2 vibrational frequency shift of ∼−44 cm−1 calculated for the largest clathrate domain considered, with 1945 H2O molecules, that 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,22 −37 cm−1, and at 76 K,23 −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.38 

Motivated in part by this suggestion, as well as by other considerations, Qu and Bowman44 have performed diffusion Monte Carlo (DMC) calculations of the vibrational frequency shift of H2 encapsulated in the (rigid) small cage of the sII hydrate, without and with surrounding water molecules, for the PES that included ab initio 3-body H2–H2O–H2O interactions, in addition to the 2-body H2–H2O interactions. For the largest hydrate domain considered, the inclusion of the 3-body interactions resulted in the shift of −40 ± 4 cm−1, in good agreement with the experiment. The DMC method employed by Qu and Bowman44 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 considerations45). Therefore, the calculations for the first excited vibrational state of the caged H2 were done in the fixed-node approximation, applying the “adiabatic” method of McCoy and co-workers46 to find the position of the node.

In our recent study,47 we presented the results of the first fully coupled quantum 6D calculations of the vibration-translation-rotation (VTR) eigenstates of a single flexible H2, HD, and D2 molecule entrapped in the (rigid) small cage of the sII hydrate that extended to the first excited (v = 1) vibrational state of H2. Prior to this work, it has been a widely held opinion that for molecular systems which have both high-frequency intramolecular vibrational mode(s) and low-frequency intermolecular vibrations, such as H2 in hydrate cages and hydrogen-bonded and van der Waals (vdW) molecular complexes, rigorous calculation of fundamental (and overtone) excitation(s) of their intramolecular vibrational mode(s) would be an extremely difficult and prohibitively costly task. The main source of the difficulty was seen to be the very large number of highly excited intermolecular vibrational eigenstates in the manifold of the intramolecular ground state lying below the energy of the intramolecular vibrational excitation(s), and the assumption that they all have to be converged in order to accurately compute the latter.

However, we demonstrated that, contrary to the above expectation, accurate computation of the intramolecular stretch fundamental of the entrapped H2 at ≈ 4100 cm−1 required having only a modest number of converged TR states in the v = 0 manifold up to at most 400–450 cm−1 above the ground state and none within several thousand wave numbers of the intramolecular fundamental.47 Our explanation for this most surprising finding was that although the number of highly excited intermolecular v = 0 TR states in the vicinity of the H2 stretch fundamental is large, their coupling to it is extremely weak. Consequently, not including them in the calculations has a negligible effect on the accuracy with which the intramolecular vibrational excitation is calculated. Of course, this resulted in a major reduction in the basis size required and transformed a computationally very demanding task to one that was readily tractable.47 

Exploiting the valuable insight gained in the above study, we recently developed a very efficient method for full-dimensional and fully coupled quantum calculations of intramolecular vibrational fundamentals and overtones of weakly bound molecular dimers.48 Its application to the 6D problem of (HF)2 produced results in excellent agreement with those in the literature, but with a fraction of the basis sets required by the other methods.

In this paper, we go one step further and report the results of the calculations of the intramolecular stretch fundamentals of (p-H2)2 confined inside the large cage of sII clathrate hydrate. These challenging calculations, like those for (HF)2,48 are made possible by what we learned in our investigation of a single H2 in the small clathrate hydrate cage.47 They also borrow from the computational methodology implemented in the (HF)2 calculations.48 The quantum treatment is in 8D (instead of 12D) since it assumes that the H2 moieties are rotationless (j = 0) and that the hydrate framework is rigid. Apart from that, the approach is fully coupled. The 8D product contracted basis includes a very small number of the lowest-energy eigenstates of the reduced-dimension 6D intermolecular Hamiltonian obtained from the full 8D Hamiltonian by fixing the bond lengths of the H2 moieties. The (H2)2-clathrate interaction potential is generated in a pairwise-additive fashion, by combining the ab initio H2–H2O pair potential for flexible H2 and rigid H2O47 and the 6D H2–H2 potential energy surface.49 The quantum 8D calculations are performed for three clathrate hydrate domains of increasing size, the largest of which contains 98 H2O molecules. They yield the vibrational stretching fundamentals of the caged (p-H2)2 and their redshifts from the free-H2 stretch frequency that for the largest domain considered are in excellent agreement with the measured values. We also formulate a first-order perturbation theory (PT) treatment, whose results are extremely close to those from the quantum 8D calculations.

This paper is organized as follows: The general approach is described in Sec. II, and the computational details are described in Sec. III. In Sec. IV, we present and discuss the results. Section V summarizes the work and outlines possible directions of further research.

As stated in the Introduction, for n (flexible) H2 molecules inside a nanocavity, taken to be rigid, the dimensionality of the quantum bound-state treatment is 6nD. Therefore, computing the coupled vibration-translation-rotation (VTR) eigenstates of (H2)2 in a cage represents a 12D problem. This poses a serious computational challenge, especially when one is interested in the excited intramolecular stretching vibrational eigenstates of (H2)2. In order to alleviate this problem, in the following, the two caged p-H2 molecules are treated as rotationless, i.e., in ground (j = 0) rotational state, and thus effectively as spherical particles. This reduces the dimensionality of the problem to be solved from 12D for two rotating H2 molecules to 8D for two rotationless (j = 0) H2 monomers. This approximation is made rather commonly, e.g., in the well known Silvera-Goldman isotropic effective pair potential for H2.50 Moreover, as demonstrated in Sec. IV A for a single p-H2 in the large clathrate hydrate cage, the error introduced by the j = 0 approximation is not large, and it decreases further as the size of the hydrate domain increases.

The 8D Hamiltonian for the coupled VRT motions of (H2)2 inside a large cage, isolated or within a larger sII hydrate domain, assumed to be rigid, and for rotationless (j = 0) H2-monomer moieties, can be written as

(1)

Here, Ri is the position vector of the center of mass (c.m.) of the ith H2 moiety measured with respect to a cage-fixed axis system with origin at cage center, i2 is the Laplacian associated with that vector, VH2-domain(4D) is the interaction energy between a rotationless H2 moiety and the hydrate domain, ri is the distance between the two H nuclei in the ith H2, VHH is the intramolecular PES of monomer H2 in its ground electronic state, R12 ≡ |R1R2|, VH2H2 is the isotropic interaction energy between the two H2 moieties, m = 2.015 65 amu is the mass of H2, and μ = 0.503 912 5 amu is the reduced mass of H2.

The approach we take in order to solve for the eigenstates and eigenvalues of the 8D Ĥ in Eq. (1) adopts some of the key elements of our recent full-dimensional treatment of the intramolecular vibrational eigenstates of weakly bound molecular dimers,48 specifically (HF)2. We first rewrite Ĥ by separating out its six-dimensional (6D) “intermolecular” portion from the remainder

(2)

where the constant r¯ is an intramolecular HH distance close to the bottom of the VHH(ri) well, and

(3)
(4)

and

(5)

Second, as in Ref. 48, we diagonalize the 6D Ĥinter in Eq. (3) to obtain the intermolecular eigenfunctions ∣κ⟩ and associated eigenvalues Eκinter. Third, we construct a contracted 8D basis, consisting of states ∣κ, γ1, γ2⟩, as the direct product of the Ninter lowest-energy ∣κ⟩ and two one-dimensional (1D) Morse discrete variable representations (DVRs) covering the r1 and r2 coordinates, respectively, and consisting of NMorse functions apiece

(6)

Finally, we compute the matrix of Ĥ in this basis and diagonalize.

The Ĥinter matrix in this basis is diagonal with elements given by

(7)

The matrix elements of the other pieces of Ĥ are easily written as follows:

(8)
(9)

and

(10)

The efficiency of this approach in calculating accurate intramolecular fundamental frequencies for the (H2)2@cage species depends crucially on the assumption of weak coupling between the intramolecular vibrational fundamentals of the confined (H2)2 and its intermolecular vibrational excitations. If this assumption is fulfilled, then the low-energy intermolecular vibrational manifold is not particularly sensitive to the intramolecular vibrational excitations of the caged H2 moieties and can be converged rapidly for both the ground and excited intramolecular vibrations using a compact contracted 6D basis of rigid-monomer intermolecular vibrational eigenstates. Moreover, weak coupling between the intra- and intermolecular vibrational DOFs of the confined (H2)2 is the key prerequisite for the fulfillment of the expectation, stemming from our recent work on weakly bound species,47,48 that intramolecular vibrational levels can be obtained accurately from full-dimensional quantum calculations that converge only a fraction of intermolecular vibrational eigenstates with the energies far below those of intramolecular vibrational fundamentals (and overtones) of interests. Based on our study of (HF)2,48 in the weak-coupling regime, it should be possible to compute accurate intramolecular fundamentals of the (H2)2@cage species by including only a small number of 6D (rigid-monomer) intermolecular vibrational eigenfunctions of Ĥinter in Eq. (3) in the 8D basis that span a range of energies much lower those of the intramolecular vibrational excitations. Our results discussed below bear out this expectation to a very high degree.

The four-dimensional VH2-domain(4D)(Ri,ri) PES describing the interaction of j = 0 H2 moiety #i with the clathrate-water moieties is taken as

(11)

where the sum is over the water moieties in the sII hydrate domain chosen, whose number is NH2O, Ξk denotes the (fixed) coordinates of water #k, and ωi denotes the two angles that fix the orientation of H2 #i with respect to the cage-fixed axis system. The 6D pair potential Vh,wk(2b) has the same meaning as in our recent work.47 It is obtained from the full 9D H2–H2O ab initio pair potential of Valiron et al.43 by fixing the intramolecular coordinates of the H2O moiety to their ground-state values. The angle brackets in Eq. (11)ωi—signify taking the average over all ωi. This averaging, appropriate for j = 0 H2, was performed by Lebedev quadrature on a grid of 26 points. Three different spherical domains of increasing size and growing number of H2O molecules (NH2O) were considered, all extracted from the 3 × 3 × 3 supercell of the sII hydrate similar to that done by Powers et al.15 in studies of H2 in the small cage of the sII hydrate. The smallest domain (NH2O=28) consists of just the H2O moieties comprising the hexakaidecahedral large cage of the sII hydrate. The next larger domain (NH2O=44) encompasses all such moieties in the sII hydrate structure whose O atoms lie within 7.5 Å of the center of the large cage. The largest domain considered (NH2O=98) includes all those H2O molecules with O atoms within 9 Å of the cage center.

For VHH, we use the identical one-body term, Vh(1b), that we employed previously in the study of H2 inside the small sII hydrate cage.47 It comes from the paper by Bowman et al.51 and is based on the work by Schwenke.52 This potential gives rise to a vibrational fundamental frequency of 4161.2 cm−1 for free-monomer H2, given the value of μ that we use.

For VH2H2, we use the isotropic part of the 6D H2-H2 PES reported by Hinde.49 This surface extends only to R12 = 4.25 bohrs on the low side. As such, for VH2H2 values corresponding to R12 < 4.25 bohrs on the 6D grid that we use for the R1, R2 coordinates (see Subsection III B), we substitute the VH2H2 value corresponding to R12 = 4.25 bohrs. We anticipate minimal errors due to this approximation since the potential at these small R12 values sufficiently exceeds that at the potential minimum that low-energy intermolecular wave functions should have negligible amplitude at these distances.

For Ĥinter, we use Eq. (3) with r¯=0.74 Å, a value close to that corresponding to the minimum of VHH. We diagonalize the operator in a 6D product basis consisting of 1D sine DVRs53 covering the six Cartesian coordinates associated with R1 and R2. These 1D DVRs are all of identical form, with each derived from the 15 lowest-energy eigenfunctions of a particle in a box whose interior ranges from −2.7 to +2.7 Å. Diagonalization was affected by using the Chebyshev version54 of filter diagonalization.55 This procedure involves the repeated application of the Hamiltonian on a random initial state vector. We symmetrized the initial state vector such that it was either symmetric or antisymmetric with respect to the H2 interchange prior to implementing the filter diagonalization algorithm. Thus, the complete diagonalization of Ĥinter consisted of separate calculations yielding “even” (∣κ+⟩) and then the corresponding “odd” (∣κ⟩) intermolecular eigenfunctions. [While the odd intermolecular states are not physical on their own, they can combine with (H2)2 intramolecular states that are also antisymmetric with respect to H2 interchange to produce physically realizable states. Hence, they must be included in the construction of the full 8D basis.]

Operation with Ĥinter on the state vector was accomplished by matrix-on-vector multiplication. The matrix elements of the pieces of the kinetic-energy operator in Ĥinter were easily obtained by first evaluating them analytically in the basis of 1D particle-in-box eigenfunctions and then transforming the results to the DVR representation. Potential-energy matrix elements are diagonal in the 6D DVR basis and are given by the value of the potential at the relevant DVR point.

Table I lists properties of the 10 lowest-energy even/odd eigenfunction pairs obtained by the diagonalization of Ĥinter for the NH2O=98 domain. The results for the domains with NH2O=28 and 44 are very similar except for overall near-constant shifts in the absolute energies of the levels to higher values. This low-energy level structure is qualitatively similar to that previously computed28 for (p-H2)2 in the large cage of sII hydrate, for a different PES and a symmetrized cage geometry, and the H2 moieties taken to be rigid. That is, the excitations are in the range of tens of wave numbers and correspond essentially to hindered rotations of the (H2)2 pseudodiatom. Notably, the interchange-tunneling splittings of the states in Table I are all very small, and the properties of the two states in each tunneling pair are almost identical. Apparently, the barrier to H2 interchange is large compared to the excitation energies of these low-energy states.

TABLE I.

Properties of the low-energy eigenstates of Ĥinter for the sII clathrate hydrate domain with NH2O=98. (Energies in cm−1 and distances in bohr.)

κ±ΔEaX¯b(ΔX)cȲY)Z¯(ΔZ)R¯12(ΔR12)d
0+ 0.000 −0.026 (0.648) −0.335 (0.407) 0.045 (0.533) 6.155 (0.573) 
0 0.000 −0.026 (0.648) −0.335 (0.407) 0.045 (0.533) 6.155 (0.573) 
1+ 10.361 0.465 (0.731) −0.311 (0.408) −0.096 (0.501) 6.121 (0.562) 
1 10.361 0.465 (0.731) −0.311 (0.408) −0.096 (0.501) 6.121 (0.562) 
2+ 22.268 0.376 (0.561) −0.197 (0.372) 0.761 (0.610) 6.016 (0.549) 
2 22.267 0.376 (0.561) −0.197 (0.372) 0.761 (0.610) 6.016 (0.549) 
3+ 42.034 0.054 (0.677) −0.328 (0.407) 0.400 (0.650) 6.102 (0.571) 
3 42.036 0.054 (0.677) −0.328 (0.407) 0.399 (0.650) 6.102 (0.571) 
4+ 47.833 0.325 (0.639) −0.254 (0.424) 0.294 (0.614) 6.074 (0.554) 
4 47.847 0.325 (0.639) −0.254 (0.423) 0.295 (0.614) 6.074 (0.554) 
5+ 53.822 0.536 (0.594) −0.188 (0.423) 0.173 (0.632) 6.045 (0.550) 
5 53.833 0.535 (0.594) −0.187 (0.423) 0.172 (0.632) 6.045 (0.550) 
6+ 55.512 0.298 (0.714) −0.201 (0.415) −0.027 (0.582) 6.095 (0.561) 
6 55.501 0.320 (0.703) −0.199 (0.414) −0.017 (0.602) 6.090 (0.561) 
7+ 55.679 −0.008 (0.675) −0.160 (0.479) 0.117 (0.649) 6.095 (0.568) 
7 55.654 −0.027 (0.676) −0.161 (0.481) 0.107 (0.631) 6.099 (0.568) 
8+ 61.361 0.319 (0.706) −0.198 (0.465) 0.204 (0.680) 6.079 (0.562) 
8 61.350 0.320 (0.705) −0.197 (0.465) 0.204 (0.679) 6.078 (0.562) 
9+ 63.695 0.104 (0.734) −0.207 (0.447) 0.109 (0.786) 6.070 (0.561) 
9 63.685 0.100 (0.734) −0.208 (0.447) 0.107 (0.785) 6.070 (0.561) 
κ±ΔEaX¯b(ΔX)cȲY)Z¯(ΔZ)R¯12(ΔR12)d
0+ 0.000 −0.026 (0.648) −0.335 (0.407) 0.045 (0.533) 6.155 (0.573) 
0 0.000 −0.026 (0.648) −0.335 (0.407) 0.045 (0.533) 6.155 (0.573) 
1+ 10.361 0.465 (0.731) −0.311 (0.408) −0.096 (0.501) 6.121 (0.562) 
1 10.361 0.465 (0.731) −0.311 (0.408) −0.096 (0.501) 6.121 (0.562) 
2+ 22.268 0.376 (0.561) −0.197 (0.372) 0.761 (0.610) 6.016 (0.549) 
2 22.267 0.376 (0.561) −0.197 (0.372) 0.761 (0.610) 6.016 (0.549) 
3+ 42.034 0.054 (0.677) −0.328 (0.407) 0.400 (0.650) 6.102 (0.571) 
3 42.036 0.054 (0.677) −0.328 (0.407) 0.399 (0.650) 6.102 (0.571) 
4+ 47.833 0.325 (0.639) −0.254 (0.424) 0.294 (0.614) 6.074 (0.554) 
4 47.847 0.325 (0.639) −0.254 (0.423) 0.295 (0.614) 6.074 (0.554) 
5+ 53.822 0.536 (0.594) −0.188 (0.423) 0.173 (0.632) 6.045 (0.550) 
5 53.833 0.535 (0.594) −0.187 (0.423) 0.172 (0.632) 6.045 (0.550) 
6+ 55.512 0.298 (0.714) −0.201 (0.415) −0.027 (0.582) 6.095 (0.561) 
6 55.501 0.320 (0.703) −0.199 (0.414) −0.017 (0.602) 6.090 (0.561) 
7+ 55.679 −0.008 (0.675) −0.160 (0.479) 0.117 (0.649) 6.095 (0.568) 
7 55.654 −0.027 (0.676) −0.161 (0.481) 0.107 (0.631) 6.099 (0.568) 
8+ 61.361 0.319 (0.706) −0.198 (0.465) 0.204 (0.680) 6.079 (0.562) 
8 61.350 0.320 (0.705) −0.197 (0.465) 0.204 (0.679) 6.078 (0.562) 
9+ 63.695 0.104 (0.734) −0.207 (0.447) 0.109 (0.786) 6.070 (0.561) 
9 63.685 0.100 (0.734) −0.208 (0.447) 0.107 (0.785) 6.070 (0.561) 
a

Relative to E0 = −1028.736 cm−1.

b

X¯, Ȳ, and Z¯ are expectation values of the Cartesian components of (R1 + R2)/2.

c

ΔXX2¯(X¯)2 with analogous definitions for ΔY and ΔZ.

d

R¯12 is the expectation value of R12. ΔR12R122¯(R¯12)2.

Further insight into the low-energy intermolecular level structure can be gained by considering the single-H2 c.m. probability densities (PDs)

(12)

associated with the states. [Note that although ρκ in Eq. (12) is for H2 #1, the same density function pertains to H2 #2 owing to interchange symmetry.] Isosurface plots of ρκ are presented in Fig. 1 for several of the low-energy states listed in Table I. From the figure, one notes that the intermolecular ground state, ∣0+⟩, has its H2 moieties localized in two regions [Fig. 1(a)]. This ground-state behavior is different from that reported in previous studies25,28 in which a different PES was assumed and in which delocalization of the H2 c.m. positions over four sites was found. As the other plots in Fig. 1 show, however, the states at higher energies show increasing degrees of angular delocalization, consistent with excitations involving hindered rotational motion.

FIG. 1.

Single-H2 center-of-mass probability densities for the intermolecular eigenfunctions: (a) ∣0+⟩, (b) ∣2+⟩, (c) ∣4+⟩, and (d) ∣9+⟩.

FIG. 1.

Single-H2 center-of-mass probability densities for the intermolecular eigenfunctions: (a) ∣0+⟩, (b) ∣2+⟩, (c) ∣4+⟩, and (d) ∣9+⟩.

Close modal

As per Eq. (6), the 8D basis for the diagonalization of Ĥ was constructed from states of Table I together with those constituting two 1D Morse DVRs of eight functions apiece. Thus, the largest 8D basis that we used consisted of Ntot = Ninter × NMorse × NMorse = 20 × 8 × 8 = 1280 states. Each Morse DVR was constructed56 by first solving the vibrational Schrödinger equation for a particle of mass μ = 0.504 amu moving in the Morse potential

(13)

with D = 0.1744 hartree, α = 1.027 64 bohr−1, and re = 1.402 01 bohr that is similar to that of the free H2 monomer. The DVR was then obtained by diagonalizing the matrix of zeα(rire) in the finite basis representation (FBR) consisting of the eight lowest-energy Morse eigenfunctions. We tested the Morse DVR by using it as the basis for the calculation of the monomer H2 vibrational states associated with the VHH potential employed herein. Results converged to within ∼ 0.1 cm−1 were obtained for the v = 0 and v = 1 states.

More can now be said about the calculation of the matrix elements that appear in Eqs. (8)–(10) and that contribute to the matrix of Ĥ in the 8D basis. The 1D matrix elements on the rhs of Eq. (8) can be readily evaluated by computing the matrix of the relevant operator in the Morse FBR basis and then transforming it to the DVR representation. The 3D and 6D matrix elements on the rhs of Eqs. (9) and (10), respectively, are easily obtained by quadrature by expressing the intermolecular eigenfunctions (i.e., ∣κ′⟩, ∣κ⟩) in terms of the primitive 6D DVR intermolecular basis. In regard to the evaluation of Eq. (10), we also note that many of the R12 values for different (R1, R2) pairs on the 6D DVR grid are identical (or almost so, i.e., to within 0.001 bohr). As such, rather than having to calculate 156 values of VH2H2 to compute Eq. (10), only 1175 such values were required.

We used direct diagonalization to compute the eigenvectors and eigenvalues of Ĥ. Several basis sets, apart from the one with Ntot = 20 × 8 × 8, were used to test convergence. All the basis sets differed only in respect to the number of intermolecular states (Ninter) included. In all cases, if one intermolecular state of an interchange-tunneling pair was included in the basis, then the other one was included as well. Erroneous results obtain if this is not done because (as we show below in Sec. IV B and Table IV) vibrational eigenstates can have significant contributions from both ∣κ+⟩ and ∣κ⟩ intermolecular states.

While our main concern here is with the intramolecular vibrational fundamentals of (H2)2 moieties in the large cage, there are two reasons why similar calculations pertaining to a single H2 in that cage are relevant as well. First, since it is straightforward to perform full-dimensional calculations on the single-H2 system,47 one can investigate the magnitude of the errors that arise from making the j = 0 approximation. Second, determining the trend in frequency shift with the H2 occupancy of the large cage obviously requires the single-H2 data point.

We have performed calculations on p-H2 inside the large cage similar to those described in our previous work47 on H2 in the small cage of the sII clathrate hydrate. We use filter diagonalization to diagonalize the full 6D vibrational Hamiltonian

(14)

where R is the c.m. position vector of the H2 moiety, ∇2 is the Laplacian associated with R, r is the HH internuclear distance, ω denotes the two angles that fix the H2 axis with respect to the cage axis system, ĵ2 is the operator associated with the square of the rotational angular momentum of the H2 moiety, and

(15)

This potential involves the same VHH and the same H2–H2O interaction, Vh,wk(2b)(R,ω,r), that we use for the two-H2 system [see Eq. (11)] and that we used in Ref. 47, except that there is no averaging over H2 orientation in Eq. (15). Clathrate hydrate domains identical to those defined in Sec. III A and employed in the (H2)2 calculations were also used in evaluating Eq. (15).

For the 6D bases, we used (a) the product of three 1D sine DVRs identical to those described in Sec. III B to cover the R degrees of freedom times, (b) a 1D Morse DVR identical to that described in Sec. III C to cover the r degree of freedom times, and (c) a spherical harmonic Yjm(ω) to cover the H2 rotational degrees of freedom. In one 6D basis, we only included the j = 0 spherical harmonic. In a second one, we included all spherical harmonics with j = 0, 2, and 4.

Table II lists results from the vibrational calculations on a single p-H2 in the large cage. The results pertain to all three clathrate hydrate domains and are given for both the rotationless (jmax = 0) basis and for the basis that allows for H2 rotational motion (jmax = 4). There are three significant trends to note from these results. First, both the ground-state energy and the frequency of the intramolecular stretch fundamental decrease with increasing NH2O although the decrease for the latter is very small. This behavior matches that which has previously been observed in computational studies38,47 of one H2 in the small cage of the sII clathrate hydrate. It is readily rationalized in that the H2–H2O interactions relevant to the hydrate inclusion compounds tend to be predominantly attractive and such as to lower the H2 intramolecular fundamental frequency. Hence, the more such interactions, the lower the ground-state energy and the smaller the frequency of the H2 stretch fundamental(s). Of course, one expects these effects to saturate as NH2O, given that the long-range H2–H2O interactions fall off considerably faster than (distance)−3.

TABLE II.

Quantum 6D frequencies ν (in cm−1) of the intramolecular stretch fundamental (v = 1) of a single p-H2 inside the large cage within three sII hydrate domains having different numbers of H2O molecules (NH2O). Also shown are their respective frequency shifts Δν relative to the free-H2 stretch fundamental. E0 (in cm−1) denotes the ground-state energy of H2 in the given domain from the quantum 6D calculations, relative to that of the free H2 monomer, computed for each of the three domains. Two values are displayed for every ν and Δν, as well as E0, the first corresponding to the jmax = 0 basis and the second to jmax = 4. For additional explanations, see the text.

NH2O284498
E0 −517.6, −523.8 −535.0, −539.1 −562.2, −563.7 
ν 4134.7, 4132.7 4133.9, 4132.5 4133.1, 4132.4 
Δν −26.5, −28.5 −27.3, −28.7 −28.1, −28.8 
NH2O284498
E0 −517.6, −523.8 −535.0, −539.1 −562.2, −563.7 
ν 4134.7, 4132.7 4133.9, 4132.5 4133.1, 4132.4 
Δν −26.5, −28.5 −27.3, −28.7 −28.1, −28.8 

Second, the jmax = 0 results are not fully converged with respect to the H2 rotational basis. However, the differences between the jmax = 0 and jmax = 4 values narrow considerably as NH2O increases. Indeed, for the largest hydrate domain considered (NH2O=98), the vibrational fundamental changes by only −0.7 cm−1 in going from the small (jmax = 0) to the large (jmax = 4) basis. This trend suggests that in the bulk clathrate hydrate, the difference between the jmax = 0 and jmax = 4 results would be negligible. The better j = 0 convergence with a larger domain can be attributed to a decrease in the angular anisotropy in the H2-hydrate interaction as NH2O increases. An analogous angular anisotropy decrease with NH2O (as reflected in a decrease in the magnitude of the H2j = 1 rotational-level splitting) was also found to hold for H2 in the small cage.38,47 The trend is easily understood once one recognizes that the angular anisotropy in the H2-cage interaction is due to the asymmetry in the spatial distribution of H-atoms in the hydrogen-bond framework of the hydrate. The effects of such asymmetry tend to get averaged away as more layers of water moieties are included in the larger hydrate domains. In any case, the single-H2 results indicate that the j = 0 approximation that we have made in the (H2)2 calculations causes convergence errors of small and acceptable magnitude.

Finally, we note that the magnitudes of the HH-frequency shifts in Table II are significantly less than those that we have computed (∼−42 cm−1) for H2 in the small hydrate cage.47 This is in keeping with the relative ground-state energies in the two cases (the small-cage H2 ground state is stabilized by approximately 200 cm−1 more than the large-cage one relative to monomer H2), and the fact that the redshift of the stretch fundamental tracks with ground-state stabilization. Indeed, the ratio of redshifts for the two cases closely matches the ratio of the ground-state stabilization energies.

Table III lists results from the diagonalization of Ĥ for the j = 0, (p-H2)2 system for two 8D bases whose dimensions differ by a factor of ten—one corresponding to Ninter = 2 (the two lowest-energy states in Table I) and Ntot = 128 and the other to Ninter = 20 (all the states in Table I) and Ntot = 1280—and for the three different hydrate domains. The results pertain to the ground-state (∣ψ0⟩) of the encaged (H2)2 and to its two interchange-symmetry-allowed v1 + v2 = 1 intramolecular excited vibrational states (|ψνa and |ψνb, respectively). From the table, one notes first that there is remarkably little difference between the results, intramolecular stretch fundamentals, and frequency shifts, obtained for the small basis (Ntot = 128) and the corresponding ones for the large basis (Ntot = 1280). Thus, the computed values appear to be very well converged. Moreover, this is strong evidence that the low-energy intermolecular excitations of the confined (H2)2 cluster are not appreciably affected by the intramolecular excitations of the H2 moieties.

TABLE III.

Energies (in cm−1) from the quantum 8D calculations (jmax = 0 basis), for the ground state E0 (relative to that of two free H2 monomers) and the two v1 + v2 = 1 intramolecular stretching excited states νa and νb of (p-H2)2 inside the large cage within three sII hydrate domains having different numbers of H2O molecules (NH2O). νa and νb denote the low- and high-energy intramolecular stretch fundamentals, respectively, while Δν is the shift of the average of νa and νb from the stretch fundamental of the free H2 monomer. For all computed quantities, the values shown in parentheses correspond to Ninter = 2 and all others to Ninter = 20. The measured value of intramolecular stretch fundamental for (H2)2 in the large sII hydrate cage is from Ref. 23; νa/b are not resolved experimentally. For additional explanations, see the text.

NH2O284498Expt.
E0 −961.74 (−961.72) −999.45 (−999.42) −1056.17 (−1056.15)  
νa 4139.24 (4139.28) 4138.26 (4138.31) 4137.34 (4137.39) 4135.5 
νb 4139.77 (4139.87) 4138.82 (4138.92) 4137.79 (4137.89) … 
Δν −21.70 (−21.63) −22.66 (−22.59) −23.64 (−23.56) −25.7 
NH2O284498Expt.
E0 −961.74 (−961.72) −999.45 (−999.42) −1056.17 (−1056.15)  
νa 4139.24 (4139.28) 4138.26 (4138.31) 4137.34 (4137.39) 4135.5 
νb 4139.77 (4139.87) 4138.82 (4138.92) 4137.79 (4137.89) … 
Δν −21.70 (−21.63) −22.66 (−22.59) −23.64 (−23.56) −25.7 

In this regard, it is instructive to examine the intermolecular basis-state composition of the pertinent 8D eigenfunctions computed by using the Ninter = 20 6D intermolecular basis. To do so, we calculate for 8D eigenstate ∣ψn⟩ the quantities

(16)

For ∣ψ0⟩, the only 6D intermolecular state that contributes appreciably is ∣0+⟩, and P0+(ψ0)=0.9986. For the intramolecular stretch fundamental states, |ψνa and |ψνb, just two lowest-energy 6D intermolecular states in Table I—∣0+⟩ and ∣0⟩—dominate in their contributions and P0+(ψνa/b)+P0(ψνa/b)=0.996. The implications of this are twofold: (1) The remaining nine higher-energy even/odd intermolecular eigenfunction pairs in Table I mix negligibly with |ψνa and |ψνb. (2) The 8D basis that we have used very efficiently overlaps the 8D eigenstates that have small degrees of intermolecular excitation.

These observations help explain a most striking result of this study that converged values of the intramolecular stretch fundamentals of the caged (p-H2)2 are obtained already with the smallest 6D intermolecular basis (Ninter = 2) and Ntot = 128, which is incapable of describing any excited intermolecular vibrational states of (p-H2)2. Even the largest intermolecular basis employed, with Ninter = 20, extends to only ∼64 cm−1 above the ground state, which is far below the energies of the (p-H2)2 intramolecular stretch fundamentals, around 4100 cm−1. At this energy, the density of the intermolecular vibrational states in the v1 + v2 = 0 (ground-state) intramolecular vibrational manifold of (p-H2)2 must be very high, but their coupling to the intramolecular vibrational excitations is negligible so that leaving them out of the treatment has no effect whatsoever on the calculated intramolecular stretch fundamentals. This is fully in line with the findings of our theoretical studies of H2 in the small cage of sII hydrate47 and (HF)2.48 For both systems, accurate intramolecular vibrational fundamentals (and overtones for HF dimer) were obtained although the highest-energy intermolecular vibrational levels computed lay thousands of wave numbers below the intramolecular vibrational excitations of interest.47,48

One also sees from Table III that there are clear trends in respect to the computed energies as a function of the clathrate domain. Both the ground-state energy and the vibrational fundamental frequencies decrease as NH2O increases. These trends match those referred to above in conjunction with the single-H2 results in Sec. IV A.

Finally from Table III, one notes that the calculated fundamental frequencies are close to that of the measured Raman band23 assigned to the intramolecular fundamental of (p-H2)2 in the large cage of the sII hydrate. The redshift of −23.64 cm−1 computed for the hydrate domain with 98 H2O molecules is only ≈8% smaller by magnitude than the measured value23 of −25.7 cm−1. While such quantitative agreement may be somewhat fortuitous, it is also notable that two important qualitative features of the experimental results are reproduced by the calculated ones. First, the computed νa/νb splitting is small (∼0.5 cm−1) and is thus consistent with the fact that only a single (p-H2)2 Raman band has been observed in the experimental study. Second, the observed trend of decreasing (by magnitude) redshift of the H2 stretch fundamental in going from H2@small cage to (H2)2@large cage is accurately reproduced by the present calculations together with those from our previous work.47 This general trend of decreasing redshift with increasing H2 occupancy of the hydrate cages has been discussed by us previously.47 It is worth mentioning here that from the quantum 6D bound state calculations of gas-phase (H2)2,49 the dimerization-induced H2 redshift is only 0.405 cm−1 and therefore contributes very little to the overall redshift of the caged (H2)2.

It is of some interest to get a more complete picture of the nature of the |ψνa and |ψνb intramolecular excitations than is available from their energies alone. In particular, are the excitations local within the dimer or delocalized? Equivalently, what portion of the νa/νb splitting is due to site splitting as opposed to excitation exchange? We have addressed these questions by computing the center-of-mass probability density corresponding to the vibrationally excited H2 moiety in both the |ψνa and the |ψνb states. To do this, we projected out of each 8D eigenstate |ψνx, x = a, b, the 6D function

(17)

where ∣v1, v2⟩ corresponds to an intramolecular excitation in which H2 moiety #1 has vibrational quantum number v1 and H2 moiety #2 has vibrational quantum number v2. ⟨v1, v2γ1, γ2⟩ were obtained by solving for the free H2-monomer vibrational states ∣vi⟩ in the ∣γi⟩ basis. We then computed the single-particle c.m. PD

(18)

This PD gives for each state the distribution in space of v = 1-excited H2 #1 for all possible positions of v = 0-excited H2 #2. [Given H2-interchange symmetry, the same functional form applies to the PD of v = 1-excited H2 #2 for all possible positions of v = 0-excited H2 #1—i.e., ρνx(0,1)(R2).]

Table IV gives the dominant coefficients κ,1,0|ψνx and κ,0,1|ψνx in Eq. (17) computed for x = a/b for the Ninter = 20 basis. From the sum of squares of these values (the right-most column of Table IV), it is clear that the two eigenstates correspond overwhelmingly to v1 + v2 = 1 intramolecular excitations as assigned. Calculation of ρνx(1,0)(R1) for the two states from such coefficients yields the plots shown in Fig. 2. The striking feature of these plots is the marked degree of spatial localization associated with each intramolecular excitation. Furthermore, if one compares these PDs with the ground state PD, ρν0(0,0)(R1), which is essentially ρ0+(R1) in Fig. 1(a), one sees that the “excitation locale” for |ψνa matches one of the high-density regions found for ∣ψ0⟩, whereas that for |ψνb matches the other such ground-state region. One concludes that, for the PES we have employed, the splitting of the intramolecular vibrational fundamentals in the encaged (H2)2 is dominated by effects due to inequivalent H2 sites within the clathrate cage.

TABLE IV.

Dominant values of κ,v1,v2|ψνx for x = a and b. For additional information, see Eq. (17) in the text.

x0+,1,0|ψνx 0+,0,1|ψνx0,1,0|ψνx0,0,1|ψνxSum of squares
a 0.5844 0.5844 0.3952 −0.3952 0.9955 
b −0.3941 −0.3941 0.5840 −0.5840 0.9927 
x0+,1,0|ψνx 0+,0,1|ψνx0,1,0|ψνx0,0,1|ψνxSum of squares
a 0.5844 0.5844 0.3952 −0.3952 0.9955 
b −0.3941 −0.3941 0.5840 −0.5840 0.9927 
FIG. 2.

Single-H2 center-of-mass probability densities corresponding to the vibrationally excited H2 moiety in the eigenstates (a) |ψνa and (b) |ψνb.

FIG. 2.

Single-H2 center-of-mass probability densities corresponding to the vibrationally excited H2 moiety in the eigenstates (a) |ψνa and (b) |ψνb.

Close modal

Given the results in Table IV, it is pertinent to make one last point in closing this section. Namely, if we were to choose an 8D basis consisting of just the four “zeroth-order” states, ∣κ, v1, v2⟩ = ∣ 0±, 1, 0⟩ and ∣0±, 0, 1⟩, and diagonalize Ĥ in that basis, we would capture the bulk of the eigenvector content of both intramolecular vibrational fundamentals |ψνa and |ψνb. Similarly, the character of the 8D ground state is almost entirely captured by a single state in the ∣κ, v1, v2⟩ basis: |⟨0+, 0, 0∣ψ0⟩|2 = 0.9984. The upshot is that computing the energies of these states, and hence the intramolecular stretch fundamental frequencies, would seem to fall well within the reach of first-order perturbation theory (PT). We have tested this as follows. We take the unperturbed Hamiltonian to be

(19)

This has ∣0+, 0, 0⟩, ∣0±, 1, 0⟩, and ∣0±, 0, 1⟩ as eigenstates. The perturbation then is

(20)

We then compute (a) the ground-state energy by means of nondegenerate first-order PT and (b) the energies of the |ψνx (x = a, b) vibrationally excited states by using degenerate first-order PT. The results for all three hydrate domains are given in Table V. A comparison with the quantum 8D results in Table III shows that while the absolute energies of the ground state from PT are ∼1 cm−1 higher, the frequencies of the intramolecular stretch fundamentals from the PT and variational approaches are nearly identical. These results on (H2)2 suggest that analogous implementations of perturbation theory might well be a viable way forward in respect to the even more challenging problem of computing the intramolecular vibrational fundamentals23 of (H2)3 and (H2)4 clusters in the large cage of sII clathrate.

TABLE V.

First-order perturbation theory energies (in cm−1), for the ground state E0 (relative to that of two free H2 monomers) and the two v1 + v2 = 1 intramolecular stretching excited states νa and νb of (p-H2)2 inside the large cage within three sII hydrate domains having different numbers of H2O molecules (NH2O). νa and νb denote the low- and high-energy intramolecular stretch fundamentals, respectively, while Δν is the shift of the average of νa and νb from the stretch fundamental of the free H2 monomer. The values shown in parentheses are from Table III, resulting from the quantum 8D calculations with Ninter = 20. The experimental value of Δν for (H2)2 in the large sII hydrate cage is from Ref. 23. For additional explanations, see the text.

NH2O284498Expt.
E0 −960.88 (−961.74) −998.52 (−999.45) −1054.78 (−1056.17)  
νa 4139.33 (4139.24) 4138.36 (4138.26) 4137.44 (4137.34)  
νb 4139.90 (4139.77) 4138.95 (4138.82) 4137.92 (4137.79)  
Δν −21.59 (−21.70) −22.55 (−22.66) −23.52 (−23.64) −25.7 
NH2O284498Expt.
E0 −960.88 (−961.74) −998.52 (−999.45) −1054.78 (−1056.17)  
νa 4139.33 (4139.24) 4138.36 (4138.26) 4137.44 (4137.34)  
νb 4139.90 (4139.77) 4138.95 (4138.82) 4137.92 (4137.79)  
Δν −21.59 (−21.70) −22.55 (−22.66) −23.52 (−23.64) −25.7 

We have performed fully coupled quantum 8D calculations of the intramolecular stretching fundamentals of (p-H2)2 confined inside the large cage of sII clathrate hydrate. The calculations assume rotationless (j = 0) H2 moieties and a rigid clathrate hydrate framework. The (H2)2-hydrate interaction potential employed is constructed in a pairwise-additive fashion by combining the ab initio H2–H2O pair potential47 for flexible H2 and rigid H2O and the 6D H2–H2 potential energy surface.49 The quantum 8D calculations are performed for three clathrate hydrate domains of increasing size with the number of water molecules in them ranging from 28 to 98.

The computational methodology incorporates important elements of the full-dimensional quantum treatment of the vibrational levels of weakly bound molecular dimers48 and the key insights gained from the quantum 6D calculations of the VRT levels of H2, HD, and D2 in the small cage of the sII hydrate.47 In this approach, we first solve for the eigenstates of a reduced-dimension 6D intermolecular Hamiltonian, obtained from the full 8D Hamiltonian by fixing the bond lengths of the two H2 moietes. An 8D contracted product basis for the solution of the full problem is then constructed from a small number of the lowest-energy 6D intermolecular eigenstates and two discrete variable representations (DVRs) covering the H2-monomer internuclear distances.

The approximation of treating the H2 moieties as rotationless (j = 0) is validated in part by performing quantum 6D calculations of the intramolecular stretch fundamental for a single p-H2 in the large sII hydrate cage for a rotationless (jmax = 0) basis and also for the basis that allows for the H2 rotation (jmax = 4). We find that differences between the jmax = 0 and jmax = 4 results, small to begin with, decrease further with the increase in the number of H2O molecules in the domains. For the largest hydrate domain, the H2 vibrational fundamental changes by just −0.7 cm−1 in going from the small (jmax = 0) to the large (jmax = 4) basis.

For the (H2)2 system, the two intramolecular fundamentals are computed to be site-split by about 0.5 cm−1 and to be red-shifted from the free-monomer stretch frequency by −23.64 cm−1 (for the hydrate domain with 98 H2O molecules), only ≈8% less than the value of −25.7 cm−1 observed experimentally.23 Such excellent agreement may be somewhat fortuitous in view of the complexity of the interactions in the system and the fact that no adjustable parameters are involved. Remarkably, an 8D basis as small as 128 functions, that includes just the two lowest-energy 6D intermolecular eigenstates, is found to be sufficient to achieve convergence of the fundamental frequencies to within about 0.1 cm−1. In addition, we find that first-order perturbation theory, based on 8D zeroth-order states consisting of products of intermolecular eigenstates and free-monomer vibrational eigenstates, produces results for the fundamental frequencies that match the variational results to within a few tenths of cm−1. The work thus provides reasons to believe that accurate calculations of intramolecular vibrational frequencies in inclusion compounds (clathrate hydrates, MOFs, and others) containing more than two H2 moieties may now be within reach. What also remains to be explored in the future is the role of nonadditive, three-body H2–H2O–H2O interactions, missing from the pairwise-additive interaction potential employed in this study.

Looking ahead, our recent accurate quantum calculations for a single p-H2 in the small cage of sII clathrate hydrate47 and (HF)2,48 as well as the ones for (p-H2)2 presented in this paper, allow us to draw conclusions that, we believe, have important implications for a much broader range of weakly bound molecular clusters. What these three systems, quite different in their nature, have in common is that they possess both high-energy intramolecular stretching vibrational mode(s), those of H2 or HF, whose fundamentals happen to lie around 4000 cm−1, and intermolecular vibrations with frequencies that are at least an order of magnitude lower. As a result, the intramolecular vibrational fundamentals are embedded in the dense set of highly excited intermolecular vibrational states belonging to the intramolecular ground-state manifold. The widely held view has been that all these high-lying intermolecular eigenstates have to be converged in order to obtain accurate intramolecular excitations, which would pose a formidable, nearly intractable computational problem.

However, our previous47,48 and present calculations have demonstrated that this is not the case for the three systems considered. Converged intramolecular stretch fundamentals (and overtones in the case of the HF dimer48) have been obtained by means of fully coupled quantum calculations that converged only a small number of low-lying intermolecular vibrational eigenstates in the intramolecular ground-state manifold with the energies several thousand wave numbers below those of the intramolecular vibrational fundamentals and overtones. The density of highly excited intermolecular vibrational states in the vicinity of the intramolecular excitations must be high. The fact that not including any of them in the calculations has a negligible effect on the accuracy with which the intramolecular fundamentals and overtones are computed, can be explained only by the extremely weak coupling between the intra- and intermolecular vibrational DOFs of the systems. This greatly reduces, by order(s) of magnitude, the total size of the basis set required and the computational effort involved. The large disparity between the energies of the intra- and intermolecular vibrational modes, and hence their weak coupling, present in hydrogen clathrate hydrates and (HF)2, is a general feature of weakly bound, hydrogen-bonded and van der Waals molecular clusters. This suggests that the computational approach developed in this study and in our earlier papers47,48 has a wide scope of applicability and will allow computing accurate intramolecular vibrational eigenstates, and the low-energy intermolecular vibrational states associate with them, of many weakly bound molecular dimers, and perhaps trimers as well, that would otherwise be beyond reach of rigorous quantum treatments.

Z.B. is grateful to the National Science Foundation for its partial support of this research through Grant No. CHE-1566085. P.F. is grateful to Professor Daniel Neuhauser for the generous sharing of his computer resources.

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.
Z.
Bačić
,
J. Chem. Phys.
149
,
100901
(
2018
).
12.
M.
Xu
,
Y.
Elmatad
,
F.
Sebastianelli
,
J. W.
Moskowitz
, and
Z.
Bačić
,
J. Phys. Chem. B
110
,
24806
(
2006
).
13.
M.
Xu
,
F.
Sebastianelli
, and
Z.
Bačić
,
J. Chem. Phys.
128
,
244715
(
2008
).
14.
M.
Xu
,
F.
Sebastianelli
, and
Z.
Bačić
,
J. Phys. Chem. A
113
,
7601
(
2009
).
15.
A.
Powers
,
O.
Marsalek
,
M.
Xu
,
L.
Ulivi
,
D.
Colognesi
,
M. E.
Tuckerman
, and
Z.
Bačić
,
J. Phys. Chem. Lett.
7
,
308
(
2016
).
16.
M.
Xu
,
L.
Ulivi
,
M.
Celli
,
D.
Colognesi
, and
Z.
Bačić
,
Phys. Rev. B
83
,
241403(R)
(
2011
).
17.
M.
Xu
,
L.
Ulivi
,
M.
Celli
,
D.
Colognesi
, and
Z.
Bačić
,
Chem. Phys. Lett.
563
,
1
(
2013
).
18.
D.
Colognesi
,
M.
Celli
,
L.
Ulivi
,
M.
Xu
, and
Z.
Bačić
,
J. Phys. Chem. A
117
,
7314
(
2013
).
19.
D.
Colognesi
,
A.
Powers
,
M.
Celli
,
M.
Xu
,
Z.
Bačić
, and
L.
Ulivi
,
J. Chem. Phys.
141
,
134501
(
2014
).
20.
L.
Ulivi
,
M.
Celli
,
A.
Giannasi
,
A. J.
Ramirez-Cuesta
,
D. J.
Bull
, and
M.
Zoppi
,
Phys. Rev. B
76
,
161401(R)
(
2007
).
21.
L.
Ulivi
,
M.
Celli
,
A.
Giannasi
,
A. J.
Ramirez-Cuesta
, and
M.
Zoppi
,
J. Phys.: Condens. Matter
20
,
104242
(
2008
).
22.
A.
Giannasi
,
M.
Celli
,
L.
Ulivi
, and
M.
Zoppi
,
J. Chem. Phys.
129
,
084705
(
2008
).
23.
T. A.
Strobel
,
E. D.
Sloan
, and
C. A.
Koh
,
J. Chem. Phys.
130
,
014506
(
2009
).
24.
A.
Giannasi
,
M.
Celli
,
M.
Zoppi
,
M.
Moraldi
, and
L.
Ulivi
,
J. Chem. Phys.
135
,
054506
(
2011
).
25.
F.
Sebastianelli
,
M.
Xu
, and
Z.
Bačić
,
J. Chem. Phys.
129
,
244706
(
2008
).
26.
A.
Witt
,
F.
Sebastianelli
,
M. E.
Tuckerman
, and
Z.
Bačić
,
J. Phys. Chem. C
114
,
20775
(
2010
).
27.
A.
Valdeś
and
G. J.
Kroes
,
J. Phys. Chem. C
116
,
21664
(
2012
).
28.
P. M.
Felker
,
J. Chem. Phys.
141
,
184305
(
2014
).
29.
P. M.
Felker
,
J. Chem. Phys.
138
,
174306
(
2013
).
30.
D.
Benoit
,
D.
Lauvergnat
, and
Y.
Scribano
,
Faraday Discuss.
212
,
533
(
2018
).
31.
S. A.
FitzGerald
,
K.
Allen
,
P.
Landerman
,
J.
Hopkins
,
J.
Matters
,
R.
Myers
, and
J. L. C.
Rowsell
,
Phys. Rev. B
77
,
224301
(
2008
).
32.
S. A.
FitzGerald
,
J.
Hopkins
,
B.
Burkholder
,
M.
Friedman
, and
J. L. C.
Rowsell
,
Phys. Rev. B
81
,
104305
(
2010
).
33.
J.
Wang
,
H.
Lu
, and
J. A.
Ripmeester
,
J. Am. Chem. Soc.
131
,
14132
(
2010
).
34.
J.
Wang
,
H.
Lu
,
J. A.
Ripmeester
, and
U.
Becker
,
J. Phys. Chem. C
114
,
21042
(
2010
).
35.
K. R.
Ramya
and
A.
Venkatnathan
,
J. Chem. Phys.
138
,
124305
(
2013
).
36.
N.
Plattner
and
M.
Meuwly
,
J. Chem. Phys.
140
,
024311
(
2014
).
37.
Z.
Futera
,
M.
Celli
,
L.
del Rosso
,
C. J.
Burnham
,
L.
Ulivi
, and
N. J.
English
,
J. Phys. Chem. C
121
,
3690
(
2017
).
38.
A.
Powers
,
Y.
Scribano
,
D.
Lauvergnat
,
E.
Mebe
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
148
,
144304
(
2018
).
39.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
6359
(
1994
).
40.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
10181
(
1994
).
41.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
103
,
1829
(
1995
).
42.
Z.
Bačić
,
J. Chem. Soc., Faraday Trans.
93
,
1459
(
1997
).
43.
P.
Valiron
,
M.
Wernli
,
A.
Faure
,
L.
Wiesenfeld
,
C.
Rist
,
S.
Kedzuch
, and
J.
Noga
,
J. Chem. Phys.
129
,
134306
(
2008
).
44.
C.
Qu
and
J. M.
Bowman
,
J. Phys. Chem. A
123
,
329
(
2019
).
45.
A.
Sarsa
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
Phys. Rev. Lett.
88
,
123401
(
2002
).
46.
H.-S.
Lee
,
J. M.
Herbert
, and
A. B.
McCoy
,
J. Chem. Phys.
110
,
5481
(
1999
).
47.
D.
Lauvergnat
,
P.
Felker
,
Y.
Scribano
,
D. M.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
150
,
154303
(
2019
).
48.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
151
,
024305
(
2019
).
49.
R. J.
Hinde
,
J. Chem. Phys.
128
,
154308
(
2008
).
50.
I. F.
Silvera
and
V. V.
Goldman
,
J. Chem. Phys.
69
,
4209
(
1978
).
51.
Z.
Homayoon
,
R.
Conte
,
C.
Qu
, and
J. M.
Bowman
,
J. Chem. Phys.
143
,
084302
(
2015
).
52.
D. W.
Schwenke
,
J. Chem. Phys.
89
,
2076
(
1988
).
53.
D. T.
Colbert
and
W. H.
Miller
,
J. Chem. Phys.
96
,
1982
(
1992
).
54.
V. A.
Mandelshtam
and
H. S.
Taylor
,
J. Chem. Phys.
106
,
5085
(
1997
).
55.
M. R.
Wall
and
D.
Neuhauser
,
J. Chem. Phys.
102
,
8011
(
1995
).
56.
H.
Wei
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
97
,
3029
(
1992
).