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) H_{2} 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* H_{2}–H_{2}O pair potential for flexible H_{2} and rigid H_{2}O [D. Lauvergnat *et al.*, J. Chem. Phys. **150**, 154303 (2019)] and the six-dimensional (6D) H_{2}–H_{2} 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 H_{2} 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 H_{2}-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 H_{2}O 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-H_{2} 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.

## I. INTRODUCTION

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 H_{2}O molecules and denoted 5^{12} due to its 12 pentagonal faces, and (b) eight large cages, each formed by 28 H_{2}O molecules arranged in 12 pentagonal and 4 hexagonal faces and therefore denoted as 5^{12}6^{4}. The small cage has been shown to accommodate only one H_{2} molecule, while up to four H_{2} molecules can be encapsulated in the large cage.^{6} Hydrogen clathrate hydrates have 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 H_{2}/HD/D_{2}, 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 calculations^{12–15} and rigorous computations of the corresponding INS spectra.^{15–19} These features have also been observed experimentally in the INS^{20,21} and the Raman spectra^{22–24} of the binary tetrahydrofuran (THF) + H_{2}/HD/D_{2} sII clathrate hydrate. The quantum TR dynamics of up to four H_{2} molecules in the large hydrate cage has been studied at *T* = 0 K by means of the diffusion Monte Carlo (DMC) calculations^{25} 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 two^{27,28} and four^{29} H_{2} molecules inside the large clathrate hydrate cage. In all these calculations, both the H_{2} 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) H_{2} in the small sII hydrate cage, while taking into account the quantum delocalization of the proton nuclei of the framework water molecules arising from their hindered rotations about the fixed positions of their O atoms.

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 H_{2} intramolecular stretching vibration away from that in the gas phase. It is readily observable in the Raman spectra of the binary tetrahydrofuran (THF) + H_{2} sII hydrate, where the large cages are completely occupied by the THF while the small cages are singly occupied by H_{2} and also in those of simple sII hydrates in which H_{2} molecules are the only guests.^{10,22,23} The vibrational frequencies of H_{2} molecules encapsulated in the sII hydrates are always lower than, i.e., redshifted relative to, the gas-phase H_{2}. The largest redshift, −34 cm^{−1}, is observed in the Raman spectra of the THF + H_{2} sII hydrate and can be assigned unambiguously to the singly H_{2} occupied small cage.^{10,22,23} The same redshift of −34 cm^{−1} appearing in the Raman spectra of the simple sII hydrate is therefore attributed to H_{2} 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 H_{2} occupancy ranges between two and four. Thus, the frequency shifts are very sensitive to the number of H_{2} molecules confined in the cage. However, interpreting them in terms of a particular H_{2} occupancy of the large cages has turned out to be nontrivial, and the initial attempts^{22} proved to be incorrect. The subsequent elaborate experiments led to the assignment of these three redshifts to double, triple, and quadruple H_{2} 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 H_{2} 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 H_{2} occupancies of the cavities of nanoporous materials, and other structural as well as dynamical aspects of the entrapped H_{2}, 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 6*n*D, *n* being the number of encapsulated H_{2} 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 H_{2}-clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H_{2} molecules. In addition, both interactions must include the dependence on the H_{2} 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 H_{2} 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 H_{2} 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 H_{2} vibrational frequencies calculated in 1D for the H_{2} intermolecular coordinates taken from many snapshots of the MD simulations covered a broad distribution of frequencies that extended to that of the free H_{2} at 4155 cm^{−1}. Their maxima agree reasonably well with the experiment after a scaling factor was introduced in the calculations. Finally, classical MD simulations within the DFT framework were performed for an sII hydrate unit cell, and the H_{2} vibrational spectra were calculated by Fourier transforming the H–H bond length autocorrelation function.^{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 H_{2}.

Recently, Powers *et al.*^{38} have calculated the frequency shift of H_{2} 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 Ar_{n}HF clusters.^{39–42} The H_{2} 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 H_{2} in a (rigid) hydrate domain that depend on the vibrational state of H_{2}, $v$ = 0, or $v$ = 1, respectively. These 5D PESs were constructed using the 5D (rigid-monomer) pair potential for the interaction of H_{2} in the ground and first excited vibrational states, respectively, with H_{2}O, obtained by averaging the full-dimensional (9D) *ab initio* PES of H_{2}–H_{2}O by Valiron *et al.*^{43} over the vibrational ground state wave function of H_{2}O and the vibrational wave functions of H_{2} for $v$ = 0 and $v$ = 1, respectively. This approach rests on the assumption of dynamical decoupling between the H_{2} intramolecular vibration and the TR modes, well-justified by their large energy separation. The H_{2} vibrational frequency shift of ∼−44 cm^{−1} calculated for the largest clathrate domain considered, with 1945 H_{2}O molecules, 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 Bowman^{44} have performed diffusion Monte Carlo (DMC) calculations of the vibrational frequency shift of H_{2} encapsulated in the (rigid) small cage of the sII hydrate, without and with surrounding water molecules, for the PES that included *ab initio* 3-body H_{2}–H_{2}O–H_{2}O interactions, in addition to the 2-body H_{2}–H_{2}O interactions. 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 Bowman^{44} is well-suited for ground-state calculations, but already the first excited state poses a challenge arising from the need to locate the node in the wave function, which is generally unknown (unless it can be determined from symmetry considerations^{45}). Therefore, the calculations for the first excited vibrational state of the caged H_{2} were done in the fixed-node approximation, applying the “adiabatic” method of McCoy and co-workers^{46} 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 H_{2}, HD, and D_{2} molecule entrapped in the (rigid) small cage of the sII hydrate that extended to the first excited ($v$ = 1) vibrational state of H_{2}. 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 H_{2} 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 H_{2} 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 H_{2} 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 H_{2} 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 H_{2} 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 H_{2} moieties. The $(H2)2$-clathrate interaction potential is generated in a pairwise-additive fashion, by combining the *ab initio* H_{2}–H_{2}O pair potential for flexible H_{2} and rigid H_{2}O^{47} and the 6D H_{2}–H_{2} potential energy surface.^{49} The quantum 8D calculations are performed for three clathrate hydrate domains of increasing size, the largest of which contains 98 H_{2}O molecules. They yield the vibrational stretching fundamentals of the caged $(p-H2)2$ and their redshifts from the free-H_{2} 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.

## II. GENERAL APPROACH

As stated in the Introduction, for *n* (flexible) H_{2} molecules inside a nanocavity, taken to be rigid, the dimensionality of the quantum bound-state treatment is 6*n*D. 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*-H_{2} 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 H_{2} molecules to 8D for two rotationless (*j* = 0) H_{2} monomers. This approximation is made rather commonly, e.g., in the well known Silvera-Goldman isotropic effective pair potential for H_{2}.^{50} Moreover, as demonstrated in Sec. IV A for a single *p*-H_{2} 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) H_{2}-monomer moieties, can be written as

Here, **R**_{i} is the position vector of the center of mass (c.m.) of the *i*th H_{2} moiety measured with respect to a cage-fixed axis system with origin at cage center, $\u2207i2$ is the Laplacian associated with that vector, $VH2-domain(4D)$ is the interaction energy between a rotationless H_{2} moiety and the hydrate domain, *r*_{i} is the distance between the two H nuclei in the *i*th H_{2}, *V*_{HH} is the intramolecular PES of monomer H_{2} in its ground electronic state, *R*_{12} ≡ |**R**_{1} − **R**_{2}|, $VH2H2$ is the isotropic interaction energy between the two H_{2} moieties, *m* = 2.015 65 amu is the mass of H_{2}, and *μ* = 0.503 912 5 amu is the reduced mass of H_{2}.

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

where the constant $r\xaf$ is an intramolecular HH distance close to the bottom of the *V*_{HH}(*r*_{i}) well, and

and

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

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

The *Ĥ*_{inter} matrix in this basis is diagonal with elements given by

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

and

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 H_{2} 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.

## III. COMPUTATIONAL DETAILS

### A. Potential energy surfaces

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

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 H_{2} #*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 H_{2}–H_{2}O *ab initio* pair potential of Valiron *et al.*^{43} by fixing the intramolecular coordinates of the H_{2}O moiety to their ground-state values. The angle brackets in Eq. (11)—$\u27e8\u2026\u27e9\omega i$—signify taking the average over all *ω*_{i}. This averaging, appropriate for *j* = 0 H_{2}, was performed by Lebedev quadrature on a grid of 26 points. Three different spherical domains of increasing size and growing number of H_{2}O 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 H_{2} in the small cage of the sII hydrate. The smallest domain ($NH2O=28$) consists of just the H_{2}O 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 H_{2}O molecules with O atoms within 9 Å of the cage center.

For *V*_{HH}, we use the identical one-body term, $Vh(1b)$, that we employed previously in the study of H_{2} 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 H_{2}, given the value of *μ* that we use.

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

### B. Diagonalization of *Ĥ*_{inter}

For *Ĥ*_{inter}, we use Eq. (3) with $r\xaf=0.74$ Å, a value close to that corresponding to the minimum of *V*_{HH}. We diagonalize the operator in a 6D product basis consisting of 1D sine DVRs^{53} covering the six Cartesian coordinates associated with **R**_{1} and **R**_{2}. 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 version^{54} 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 H_{2} 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 H_{2} 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 computed^{28} for $(p-H2)2$ in the large cage of sII hydrate, for a different PES and a symmetrized cage geometry, and the H_{2} 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 H_{2} interchange is large compared to the excitation energies of these low-energy states.

κ^{±}
. | ΔE^{a}
. | $X\xaf$^{b} $(\Delta X)$^{c}
. | $\u0232$(ΔY)
. | $Z\xaf(\Delta Z)$ . | $R\xaf12\u2009(\Delta 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) |

κ^{±}
. | ΔE^{a}
. | $X\xaf$^{b} $(\Delta X)$^{c}
. | $\u0232$(ΔY)
. | $Z\xaf(\Delta Z)$ . | $R\xaf12\u2009(\Delta 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 *E*_{0} = −1028.736 cm^{−1}.

^{b}

$X\xaf$, $\u0232$, and $Z\xaf$ are expectation values of the Cartesian components of (**R**_{1} + **R**_{2})/2.

^{c}

$\Delta X\u2261X2\xaf\u2212(X\xaf)2$ with analogous definitions for Δ*Y* and Δ*Z*.

^{d}

$R\xaf12$ is the expectation value of *R*_{12}. $\Delta R12\u2261R122\xaf\u2212(R\xaf12)2$.

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

associated with the states. [Note that although *ρ*_{κ} in Eq. (12) is for H_{2} #1, the same density function pertains to H_{2} #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 H_{2} moieties localized in two regions [Fig. 1(a)]. This ground-state behavior is different from that reported in previous studies^{25,28} in which a different PES was assumed and in which delocalization of the H_{2} 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.

### C. Diagonalization of *Ĥ*

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 *N*_{tot} = *N*_{inter} × *N*_{Morse} × *N*_{Morse} = 20 × 8 × 8 = 1280 states. Each Morse DVR was constructed^{56} by first solving the vibrational Schrödinger equation for a particle of mass *μ* = 0.504 amu moving in the Morse potential

with *D* = 0.1744 hartree, *α* = 1.027 64 bohr^{−1}, and *r*_{e} = 1.402 01 bohr that is similar to that of the free H_{2} monomer. The DVR was then obtained by diagonalizing the matrix of $z\u2261e\u2212\alpha (ri\u2212re)$ 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 H_{2} vibrational states associated with the *V*_{HH} 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 *R*_{12} values for different (**R**_{1}, **R**_{2}) 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 15^{6} 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 *N*_{tot} = 20 × 8 × 8, were used to test convergence. All the basis sets differed only in respect to the number of intermolecular states (*N*_{inter}) 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.

### D. Intramolecular stretch fundamental of a single H_{2} in the large cage of sII hydrate

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 H_{2} in that cage are relevant as well. First, since it is straightforward to perform full-dimensional calculations on the single-H_{2} 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 H_{2} occupancy of the large cage obviously requires the single-H_{2} data point.

We have performed calculations on *p*-H_{2} inside the large cage similar to those described in our previous work^{47} on H_{2} in the small cage of the sII clathrate hydrate. We use filter diagonalization to diagonalize the full 6D vibrational Hamiltonian

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

This potential involves the same *V*_{HH} and the same H_{2}–H_{2}O interaction, $Vh,wk(2b)(R,\omega ,r)$, that we use for the two-H_{2} system [see Eq. (11)] and that we used in Ref. 47, except that there is no averaging over H_{2} 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(\omega )$ to cover the H_{2} 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.

## IV. RESULTS AND DISCUSSION

### A. Single H_{2} in the large sII hydrate cage

Table II lists results from the vibrational calculations on a single *p*-H_{2} in the large cage. The results pertain to all three clathrate hydrate domains and are given for both the rotationless (*j*_{max} = 0) basis and for the basis that allows for H_{2} rotational motion (*j*_{max} = 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 studies^{38,47} of one H_{2} in the *small* cage of the sII clathrate hydrate. It is readily rationalized in that the H_{2}–H_{2}O interactions relevant to the hydrate inclusion compounds tend to be predominantly attractive and such as to lower the H_{2} intramolecular fundamental frequency. Hence, the more such interactions, the lower the ground-state energy and the smaller the frequency of the H_{2} stretch fundamental(s). Of course, one expects these effects to saturate as $NH2O\u2192\u221e$, given that the long-range H_{2}–H_{2}O interactions fall off considerably faster than (distance)^{−3}.

$NH2O$ . | 28 . | 44 . | 98 . |
---|---|---|---|

E_{0} | −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 |

$NH2O$ . | 28 . | 44 . | 98 . |
---|---|---|---|

E_{0} | −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 *j*_{max} = 0 results are not fully converged with respect to the H_{2} rotational basis. However, the differences between the *j*_{max} = 0 and *j*_{max} = 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 (*j*_{max} = 0) to the large (*j*_{max} = 4) basis. This trend suggests that in the bulk clathrate hydrate, the difference between the *j*_{max} = 0 and *j*_{max} = 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 H_{2}-hydrate interaction as $NH2O$ increases. An analogous angular anisotropy decrease with $NH2O$ (as reflected in a decrease in the magnitude of the H_{2} *j* = 1 rotational-level splitting) was also found to hold for H_{2} in the small cage.^{38,47} The trend is easily understood once one recognizes that the angular anisotropy in the H_{2}-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-H_{2} 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 H_{2} in the small hydrate cage.^{47} This is in keeping with the relative ground-state energies in the two cases (the small-cage H_{2} ground state is stabilized by approximately 200 cm^{−1} more than the large-cage one relative to monomer H_{2}), 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.

### B. $(H2)2$ in the large hydrate cage

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 *N*_{inter} = 2 (the two lowest-energy states in Table I) and *N*_{tot} = 128 and the other to *N*_{inter} = 20 (all the states in Table I) and *N*_{tot} = 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 ($|\psi \nu a$ and $|\psi \nu 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 (*N*_{tot} = 128) and the corresponding ones for the large basis (*N*_{tot} = 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 H_{2} moieties.

$NH2O$ . | 28 . | 44 . | 98 . | Expt. . |
---|---|---|---|---|

E_{0} | −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 |

$NH2O$ . | 28 . | 44 . | 98 . | Expt. . |
---|---|---|---|---|

E_{0} | −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 *N*_{inter} = 20 6D intermolecular basis. To do so, we calculate for 8D eigenstate ∣*ψ*_{n}⟩ the quantities

For ∣*ψ*_{0}⟩, the only 6D intermolecular state that contributes appreciably is ∣0^{+}⟩, and $P0+(\psi 0)=0.9986$. For the intramolecular stretch fundamental states, $|\psi \nu a$ and $|\psi \nu b$, just *two* lowest-energy 6D intermolecular states in Table I—∣0^{+}⟩ and ∣0^{−}⟩—dominate in their contributions and $P0+(\psi \nu a/b)+P0\u2212(\psi \nu 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 $|\psi \nu a$ and $|\psi \nu 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 (*N*_{inter} = 2) and *N*_{tot} = 128, which is incapable of describing *any* excited intermolecular vibrational states of $(p-H2)2$. Even the largest intermolecular basis employed, with *N*_{inter} = 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 H_{2} in the small cage of sII hydrate^{47} 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-H_{2} results in Sec. IV A.

Finally from Table III, one notes that the calculated fundamental frequencies are close to that of the measured Raman band^{23} 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 H_{2}O molecules is only ≈8% smaller by magnitude than the measured value^{23} 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 H_{2} stretch fundamental in going from H_{2}@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 H_{2} 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 H_{2} 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 $|\psi \nu a$ and $|\psi \nu 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 H_{2} moiety in both the $|\psi \nu a$ and the $|\psi \nu b$ states. To do this, we projected out of each 8D eigenstate $|\psi \nu x$, *x* = *a*, *b*, the 6D function

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

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

Table IV gives the dominant coefficients $\u27e8\kappa ,1,0|\psi \nu x\u27e9$ and $\u27e8\kappa ,0,1|\psi \nu x\u27e9$ in Eq. (17) computed for *x* = *a*/*b* for the *N*_{inter} = 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 $\rho \nu 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, $\rho \nu 0(0,0)(R1)$, which is essentially $\rho 0+(R1)$ in Fig. 1(a), one sees that the “excitation locale” for $|\psi \nu a$ matches one of the high-density regions found for ∣*ψ*_{0}⟩, whereas that for $|\psi \nu 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 H_{2} sites within the clathrate cage.

x
. | $\u27e80+,1,0|\psi \nu x\u27e9$ . | $\u27e80+,0,1|\psi \nu x\u27e9$ . | $\u27e80\u2212,1,0|\psi \nu x\u27e9$ . | $\u27e80\u2212,0,1|\psi \nu x\u27e9$ . | Sum of squares . |
---|---|---|---|---|---|

a | 0.5844 | 0.5844 | 0.3952 | −0.3952 | 0.9955 |

b | −0.3941 | −0.3941 | 0.5840 | −0.5840 | 0.9927 |

x
. | $\u27e80+,1,0|\psi \nu x\u27e9$ . | $\u27e80+,0,1|\psi \nu x\u27e9$ . | $\u27e80\u2212,1,0|\psi \nu x\u27e9$ . | $\u27e80\u2212,0,1|\psi \nu x\u27e9$ . | Sum of squares . |
---|---|---|---|---|---|

a | 0.5844 | 0.5844 | 0.3952 | −0.3952 | 0.9955 |

b | −0.3941 | −0.3941 | 0.5840 | −0.5840 | 0.9927 |

### C. First-order perturbation theory treatment is highly accurate

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 $|\psi \nu a$ and $|\psi \nu 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

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

We then compute (a) the ground-state energy by means of nondegenerate first-order PT and (b) the energies of the $|\psi \nu 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 fundamentals^{23} of $(H2)3$ and $(H2)4$ clusters in the large cage of sII clathrate.

$NH2O$ . | 28 . | 44 . | 98 . | Expt. . |
---|---|---|---|---|

E_{0} | −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 |

$NH2O$ . | 28 . | 44 . | 98 . | Expt. . |
---|---|---|---|---|

E_{0} | −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 |

## V. CONCLUSIONS

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) H_{2} 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* H_{2}–H_{2}O pair potential^{47} for flexible H_{2} and rigid H_{2}O and the 6D H_{2}–H_{2} 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 dimers^{48} and the key insights gained from the quantum 6D calculations of the VRT levels of H_{2}, HD, and D_{2} 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 H_{2} 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 H_{2}-monomer internuclear distances.

The approximation of treating the H_{2} moieties as rotationless (*j* = 0) is validated in part by performing quantum 6D calculations of the intramolecular stretch fundamental for a single *p*-H_{2} in the large sII hydrate cage for a rotationless (*j*_{max} = 0) basis and also for the basis that allows for the H_{2} rotation (*j*_{max} = 4). We find that differences between the *j*_{max} = 0 and *j*_{max} = 4 results, small to begin with, decrease further with the increase in the number of H_{2}O molecules in the domains. For the largest hydrate domain, the H_{2} vibrational fundamental changes by just −0.7 cm^{−1} in going from the small (*j*_{max} = 0) to the large (*j*_{max} = 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 H_{2}O 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 H_{2} moieties may now be within reach. What also remains to be explored in the future is the role of nonadditive, three-body H_{2}–H_{2}O–H_{2}O interactions, missing from the pairwise-additive interaction potential employed in this study.

Looking ahead, our recent accurate quantum calculations for a single *p*-H_{2} in the small cage of sII clathrate hydrate^{47} 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 H_{2} 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 previous^{47,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 dimer^{48}) 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 papers^{47,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.

## ACKNOWLEDGMENTS

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.