We report the first fully coupled quantum six-dimensional (6D) bound-state calculations of the vibration-translation-rotation eigenstates of a flexible H2, HD, and D2 molecule confined inside the small cage of the structure II clathrate hydrate embedded in larger hydrate domains with up to 76 H2O molecules, treated as rigid. Our calculations use a pairwise-additive 6D intermolecular potential energy surface for H2 in the hydrate domain, based on an ab initio 6D H2–H2O pair potential for flexible H2 and rigid H2O. They extend to the first excited (v = 1) vibrational state of H2, along with two isotopologues, providing a direct computation of vibrational frequency shifts. We show that obtaining a converged v = 1 vibrational state of the caged molecule does not require converging the very large number of intermolecular translation-rotation states belonging to the v = 0 manifold up to the energy of the intramolecular stretch fundamental (≈4100 cm−1 for H2). Only a relatively modest-size basis for the intermolecular degrees of freedom is needed to accurately describe the vibrational averaging over the delocalized wave function of the quantum ground state of the system. For the caged H2, our computed fundamental translational excitations, rotational j = 0 → 1 transitions, and frequency shifts of the stretch fundamental are in excellent agreement with recent quantum 5D (rigid H2) results [A. Powers et al., J. Chem. Phys. 148, 144304 (2018)]. Our computed frequency shift of −43 cm−1 for H2 is only 14% away from the experimental value at 20 K.

Hydrogen clathrate hydrates are crystalline inclusion compounds where one or more hydrogen molecules are encapsulated inside the cavities, or cages, created by the three-dimensional (3D) framework of hydrogen-bonded water molecules.1–3 Simple hydrogen clathrate hydrates, which have only hydrogen molecules as guests, were first identified by Dyadin et al.4 and later characterized by Mao et al. in more detail.5 They have the classical structure II (sII),1,2,5 whose cubic unit cell consists of two types of cages. One of them is the small dodecahedral cage, sixteen per unit cell, each comprising 20 H2O molecules forming 12 pentagonal faces, hence designated 512. The second type are the large cages, eight per unit cell, in which 28 H2O molecules are arranged in 12 pentagonal and 4 hexagonal faces and therefore denoted 51264. Experiments have shown that the small cage can accommodate only one H2 molecule, while up to four H2 molecules can be encapsulated in the large cage.6 Hydrogen clathrate hydrates have attracted a great deal of interest in recent years, owing to their potential as economical and environmentally friendly hydrogen storage materials.1,2,7–10

Moreover, hydrogen molecules entrapped in the clathrate hydrate cages constitute fascinating and unconventional chemical systems whose dynamics and spectroscopy are thoroughly dominated by strong quantum effects, to a degree matched only by light molecules inside fullerenes.11 The pronounced quantum effects have multiple sources; one of them is the quantization of the translational center-of-mass (c.m.) degrees of freedom (DOFs) of the guest molecule(s) due to the nanoscale confinement in the clathrate cage, small or large (particle-in-a-box effect). The confining potential of the hydrate cage couples the quantized translational DOFs to the also quantized rotational DOFs of the hydrogen molecule(s). The resulting translation-rotation (TR) energy level structure is sparse, owing the low molecular mass of H2/HD/D2, their large rotational constants, and the small size of the hydrate cavities. The salient features of the TR eigenstates of a single hydrogen molecule in the cages of the sII clathrate hydrate, notably the splittings of both the translational fundamental and rotational levels, as well as their manifestations in the inelastic neutron scattering (INS) spectra, have been characterized by Bačić and co-workers through quantum 5D bound-state calculations12–15 and rigorous computations of the corresponding INS spectra.15–19 Quantum TR dynamics of multiple hydrogen molecules in the large hydrate cage has been investigated by means of the diffusion Monte Carlo (DMC)20 and path-integral molecular dynamics (PIMD) simulations21 and also by fully coupled eigenstate-resolved calculations.22–24 In all these calculations, the hydrogen-bonded clathrate hydrate framework was treated as rigid. In a recent study,25 this constraint was relaxed partially, by performing quantum 5D calculations of the TR levels of 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.

Besides giving rise to the TR energy level structure, the encapsulation of hydrogen molecules in the cages of clathrate hydrates results in the shift in the frequency of the H2 intramolecular stretching vibration away from that in the gas phase. This frequency shift is readily observable in the Raman spectra of the binary tetrahydrofuran (THF) + H2 sII hydrate, where the large cages are completely occupied by the THF while the small cages are singly occupied by H2, and simple sII hydrates in which H2 molecules are the only guests.10,26,27 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,26,27 The same redshift of −34 cm−1 appearing in the Raman spectra of the simple sII hydrate is, therefore, also attributed to H2 in the small cage.

The Raman spectra of the simple II hydrate also display bands redshifted by − 26, −18, and −11 cm−1, respectively,10,26,27 that must represent contributions from the large cages whose H2 occupancy ranges between two and four. However, associating each of these redshifts with a particular H2 occupancy of the large cages proved to be nontrivial. Initially, the redshifts of −26, −18, and −11 cm−1 were interpreted in terms of triply, doubly, and singly occupied large cages, respectively.26 Subsequent very careful experiments that involved multiple cycles of heating and quenching of the sample and measurements of the amounts of H2 released in each led to the essentially opposite assignment of these three redshifts to double, triple, and quadruple occupancies of the large cages, respectively.27 The observed trend in the H2 redshift can be understood in terms of the interplay between two kinds of interactions.27 One of them, the attractive interaction between H2 and the cage, softens the intramolecular stretch potential of H2 and lowers its vibrational frequency relative to the gas-phase. As the large-cage occupancy of the large cage increases, the tighter packing and the largely repulsive H2–H2 interactions lead to the increasing vibrational frequency of H2 and the decreasing redshift. The fact that the H2 vibrational frequency is redshifted even for the highest, quadruple occupancy of the large cage suggests that the attractive guest-host interaction always remains dominant over the repulsive H2–H2 interactions.

In the case of sII hydrogen hydrates, it has been possible to assign with confidence the observed frequency shifts to different H2 occupancies of the small and large clathrate cages largely guided by the experimental data. But, in general, e.g., molecular hydrogen in metal-organic frameworks (MOFs),28,29 reliable decoding of the information contained in the vibrational frequency shifts regarding the H2 occupancies of the cavities of nanoporous materials, and other structural as well as dynamical aspects of the entrapped H2, requires theoretical methods capable of reliably calculating the frequency shifts. This is a highly challenging task, for two reasons. First, the problem is inherently high-dimensional. Even if the hydrate framework is treated as rigid, the dimensionality of the calculations is 6nD, where n is the number of H2 molecules considered; thus, for n = 1 − 4, one has to be able to deal with the problem whose dimensionality ranges from 6D to 24D. This requires having accurate high-dimensional potential energy surfaces (PESs) that incorporate the H2–clathrate interactions and, in the case of multiple occupancy, the interactions among the guest H2 molecules. Both interactions must include the dependence on the H2 intermolecular stretch coordinate and its coupling to the intermolecular degrees of freedom. Second, dynamical quantum effects and anharmonicities in both intra- and intermolecular DOFs play a significant role, particularly at the low temperatures of the Raman spectroscopy measurements. Consequently, these key features have to be fully accounted for in any first-principles theoretical method aiming to generate accurate frequency shifts of encapsulated hydrogen molecules.

Within the past decade, a number of approaches, involving a variety of approximations, have been taken to address this fundamental and difficult problem. In some of them, the H2 molecules encapsulated in the isolated small or large hydrate cages were taken to be frozen in the geometry corresponding to the minimum energy of the system.30–32 As a result, nuclear quantum effects are left out, in particular, the averaging over the large-amplitude intermolecular vibrations of the guest H2 molecules. In addition, since only isolated clathrate cages are considered, the effects of the condensed-matter environment are unaccounted for. This problem has also been treated through a combination of classical molecular dynamics (MD) and PIMD simulations with electronic structure calculations at the density functional theory [DFT (B3LYP)] and MP2 levels.33 The 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.34 This classical treatment does not account for the quantum effects. Moreover, it gives the vibrational spectra that are shifted by 100–150 cm−1 to higher frequencies relative to the experimental results and above the stretch fundamental of free H2.

Very recently, Powers et al.35 have calculated the frequency shift of H2 inside the small cage of the sII hydrate, isolated and embedded in spherical hydrate domains of increasing size, in order to investigate the effect of the condensed-phase environment. The approach employed was developed earlier by Bačić and co-workers for the purpose of computing the HF stretch frequency shift in ArnHF clusters.36–39 The H2 frequency shift was obtained by means of the quantum 5D bound-state calculations of the coupled TR eigenstates on a pairwise-additive intermolecular PES for rigid H2 in a (rigid) hydrate domain, which depends on the vibrational state of H2, v = 0 or v = 1. This 5D PES was 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.40 over the vibrational ground state wave function of H2O and the vibrational wave functions of H2 for v = 0 and v = 1, respectively. Implicit in this approach is 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, which mimics the condensed-phase environment, was about 10% larger in magnitude than that obtained for the isolated small cage. This 0 K value agrees well with the frequency shifts measured at 20 K26 (−37 cm−1) and at 76 K27 (−34 cm−1). It was suggested that improving further the agreement with the experiment may require including many-body interactions, three-body, in particular, missing from the pairwise-additive intermolecular PES employed.35 

Motivated in part by this suggestion and other considerations, Qu and Bowman41 have included ab initio 3-body H2–H2O–H2O interactions, in addition to the 2-body H2–H2O interactions, in their 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 largest hydrate domain considered having 172 H2O water molecules, calculations based on 2-body interactions only yielded the frequency shift of about −26 cm−1, while the inclusion of the 3-body interactions resulted in the shift of −40 ± 4 cm−1, in good agreement with the experiment, in particular, −37 cm−1 at 20 K.26 

The DMC method employed by Qu and Bowman41 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 considerations42). 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-workers43 to find the position of the node. However, determining the correct nodal surface in a six-dimensional (6D) system is very difficult, virtually impractical. Therefore, Qu and Bowman made the approximation that the node is located entirely on the H–H intramolecular stretch coordinate and is independent of the TR coordinates of H2, thereby reducing the search for its position to 1D. This is equivalent to the assumption that the intra- and intermolecular coordinates of the caged H2 are decoupled, justified by the large energy separation between the two types of modes.41 

Thus, the quantum methodologies employed in the two recent computations of the H2 vibrational frequency shift in the small sII hydrate cage,35,41 although entirely different, both rely on the approximation of no coupling between the high-energy intramolecular vibrational mode of H2 and its low-energy TR modes. There is no reason to doubt its validity for this system (in both approaches) and the accuracy of the results it yields. Still, one can ask whether it is possible to perform quantum 6D calculation of the bound states of H2 in the small cage up to, and including, the energy of the first excited (v = 1) intramolecular vibrational state (the stretch fundamental), around 4100 cm−1, treating the intra- and intermolecular (TR) degrees of freedom as fully coupled, i.e., not invoking their separability. After all, fully coupled full-dimensional (6D) quantum calculations of the vibrational levels of nonrigid molecular systems, such as (HF)2,44 (HCl)2,45 and CO on Cu(100),46 have been feasible for some 25 years. In some cases, e.g., (HCl)247 and CO on Cu(100),46 the quantum 6D calculations yielded the energies of the intramolecular stretch fundamental(s), and thus their shifts from the respective gas-phase values.

Nevertheless, it has been generally thought that for molecular systems which have both high-frequency intramolecular mode(s) and low-frequency intermolecular vibrations, such as H2 in the small sII hydrate cage, and hydrogen-bonded and van der Waals (vdW) complexes, rigorous calculation of fundamental excitation(s) of their intramolecular mode(s), e.g., the v = 1 vibrational state of the encapsulated H2 that is 6D for a rigid cage, would be an extremely difficult and prohibitively costly task. The main source of the difficulty was the assumption that the very large number of highly excited intermolecular vibrational eigenstates in the manifold of the intramolecular ground state below the energy of the intramolecular excitation(s) all have to be converged in order to compute accurate fundamental intramolecular excitation(s). In this paper, we demonstrate that certainly for the intramolecular stretch fundamental of H2 (HD, D2) inside the small cage of the sII hydrate, and its frequency shift, this widely held view is not correct, making this problem entirely tractable.

We present the results of the fully coupled quantum 6D calculations of the vibration-translation-rotation (VTR) eigenstates of a single flexible H2, HD, and D2 molecule entrapped in the (rigid) small cage of the sII hydrate. We show that computing the converged energy of the first excited (v = 1) intramolecular vibrational state of the caged H2 (and the isotopologues) at ≈4100 cm−1 requires converging only the TR states in the v = 0 manifold up to at most 400–450 cm−1 above the ground state. Guided by our previous studies,15,35 quantum 6D calculations of the coupled VTR eigenstates, which extend to the v = 1 state and its frequency shift away from the gas-phase value, are performed for H2 encapsulated inside the spherical sII hydrate domains of increasing radius treated as rigid. The 6D intermolecular PES for a flexible hydrogen molecule inside the hydrate domain utilized in these calculations is constructed in a pairwise additive fashion, based on an ab initio 6D H2–H2O pair potential. The TR eigenstates and vibrational frequency shifts computed for H2, HD, and D2 are compared with the available experimental data, as well as the results of the quantum 5D calculations in Ref. 35.

This paper is organized as follows: Methodology is described in Sec. II. In Sec. III, we present and discuss the results. Section IV summarizes the work and outlines possible extensions.

The three spherical sII clathrate hydrate domains used in this work are identical to those generated previously by Bačić and co-workers15 and employed in our recent quantum 5D H2 frequency shift calculations.35 They are carved out of the 3 × 3 × 3 supercell of the sII hydrate.15 The three domains of increasing size and number of H2O molecules N have the cutoff radii set to (a) 5.0 Å, enclosing only the N = 20 water molecules of the small dodecahedra cage itself; (b) 7.5 Å, encompassing additional 20 H2O molecules, for the total of N = 40 water molecules; and (c) 9.0 Å, encompassing N = 76 water molecules that form the first three complete solvation shells around H2.15 In our previous study,35 the largest hydrate domain considered was much larger and included 1945 water molecules. However, that study also showed that the frequency shift computed for this domain differs by only 3% from that obtained for the N = 76 domain. Therefore, in the present work, we do not go beyond the sII hydrate domain containing N = 76 water molecules. In the bound-state calculations, the domains are taken to be rigid.

For the 6D intermolecular PES of flexible H2 inside the domain with N water molecules, denoted VH2domain, only one- and two-body terms of its many-body expansion are retained

VH2domainqh=Vh(1b)(r)+iNVh,wi(2b)qh,Ξi.
(1)

Here, Vh(1b)(r) is the one-body term for the intramolecular stretching coordinate (r) of the hydrogen molecule. For it, we use the corresponding one-body term in the many-body representation of the PES for H2 in the sII hydrate by Bowman and co-workers.48 

The second term in Eq. (1) represents the summation over the two-body interactions Vh,wi(2b)VH2H2O(qh,Ξi) between the hydrogen molecule and each of the N water molecules in the domain. The coordinates of H2 are qh ≡ {R, ω, r}, where R is the vector pointing from the center of the confining small cage, which is also the origin of the space-fixed (SF) Cartesian frame attached to the cage, to the center of mass (c.m.) of H2, while ω ≡ (θ, ϕ) are the polar and azimuthal angles, respectively, which specify the orientation of H2 relative to the SF frame. The position vector R can be expressed either in terms of the SF Cartesian coordinates {x, y, z}12,35 or the spherical polar coordinates {R, Ω}, where R ≡ ||R||, and Ω ≡ (Θ, Φ) are the polar and azimuthal angles of R relative to the SF axes.49Ξi denotes the coordinates of the ith water molecule in the domain; these are fixed since the domains are assumed to be rigid.

The 6D pair potential VH2H2O(qh,Ξi) is derived from the accurate full-dimensional (9D) ab initio H2–H2O pair potential V08 of Valiron et al.40 In this 9D PES, the flexibility of H2 and H2O monomers is included as a correction to the rigid-monomer dimer 5D PES, computed as a Taylor expansion around the equilibrium geometries of monomers. The 9D PES V08 was obtained by combining standard coupled-cluster with single, double, and with the perturbative contributions of the triples excitations [CCSD(T)] calculations with the explicitly correlated CCSD(T)-R12 calculations, and it is expected to provide the currently most accurate description of the H2–H2O interaction, with an accuracy of a few cm−1 in the region of the van der Waals minimum.40 In this work, the 9D H2–H2O PES was reduced to the 6D pair potential, for flexible H2 and rigid H2O, by fixing the intramolecular coordinates of H2O to their values in the ground vibrational state (OH bond length = 1.843 a0 and the HOH bending angle = 104.41°). The accuracy of this procedure is comparable to that of averaging the 9D PES over the vibrational ground-state wave function of H2O. Finally, the intermolecular coordinates employed in the V08 potential are transformed numerically to the coordinates used for the Vh,wi(2b)(qh,Ξi) term in Eq. (1).

Although the focus of this study is on the excited VTR eigenstates, we also use the diffusion Monte Carlo (DMC) method to compute in 6D the VTR ground-state energy of flexible H2, HD, and D2 inside the rigid water domains. This approach simulates a diffusion process in imaginary time on a given PES. The general DMC approach has been described in detail in Ref. 50 and here we use a standard (i.e., not rigid body) formulation of the algorithm for the caged molecule, while the cage itself remains fixed.

The simulations are performed using a revised parallelized version of the XDMC code developed by Benoit51 (see also Ref. 50 for implementation details). For each simulation in this study, we use 1000 replicas, stabilization periods of 61 500 cycles (H2), 80 900 cycles (HD), and 108 300 cycles (D2) with Δτ = 5 a.u. and an averaging phase of 1000 × 100 cycles with Δτ = 1 a.u.

The 6D Hamiltonian for the coupled VTR motions of a vibrating diatomic molecule AB, which in this study corresponds to H2 and its isotopologues HD and D2, inside a rigid clathrate hydrate domain, can be written as

Ĥ6D=22mAB222μAB2r2+ĵ22μABr2+VH2domain(R,ω,r),
(2)

where mAB and μAB are the total mass and the reduced mass of AB, respectively, while R, ω, and r were defined in Sec. II A. ∇2 is the Laplacian associated with R, and ĵ2 is the operator associated with the square of the rotational angular momentum of AB. For the isotopic masses of hydrogen, the values mH = 1.008 g mol−1 and mD = 2.0141 g mol−1 were used.

1. The Smolyak scheme approach with ELVIBROT

In most, although not all, of the calculations performed in this study, the AB c.m. position vector R in the Hamiltonian in Eq. (2) is expressed in terms of the Cartesian coordinates {x, y, z}, and ω ≡ (θ, ϕ). Furthermore, the operator ĵ2 of Eq. (2) is expanded in terms of partial derivative operators, θ and ϕ. For this choice of {x, y, z, θ, ϕ, r} coordinates, we have used the Smolyak scheme approach52 introduced by Avila and Carrington53–55 and also proposed by Lauvergnat and Nauts.56–58 More recently, it has been used to calculate the energy levels of H2 in a clathrate hydrate.25,35 In the Smolyak scheme, the single large direct-product basis or grid is replaced by a sum of small direct-products, denoted as SLSrep,

SLSrep=L=[1,,n]LSn+1|L|LS(1)LS|L|Cn1LS|L|S11Snn,
(3)

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

TABLE I.

Numbers of primitive basis functions nbi and grid points nqi (i = 1, …, 4) as a function of the parameter i. For the translational (i = 1, …, 3) and vibrational (i = 4) degrees of freedom (DOFs), the number of grid points is equal to the number of primitive functions nbi. For the rotational DOF (i = 5), the number of Lebedev points is a function of i and is explicitly given and the quantum rotational number j is such that 0 ⩽ jjmax.

i012345678
Translation (HO)          
nbi 11 13 15 17 
Vibration (HO)          
nbi 10 13 16 19 22 25 
Rotation (Yjm         
jmax 
nb5 16 25 36 49 64 81 
nq5 (Lebedev) 14 26 38 50 74 86 110 
i012345678
Translation (HO)          
nbi 11 13 15 17 
Vibration (HO)          
nbi 10 13 16 19 22 25 
Rotation (Yjm         
jmax 
nb5 16 25 36 49 64 81 
nq5 (Lebedev) 14 26 38 50 74 86 110 

The desired 6D VTR eigenstates are obtained by direct diagonalization of the Hamiltonian in Eq. (2) in this basis and computed with the help of ELVIBROT.59 

2. Filter diagonalization in a direct-product basis

An alternative approach employed in this study is to compute the VTR eigenstates of the 6D Hamiltonian in Eq. (2) in selected regions of the energy spectrum, utilizing the Chebyshev variant60 of filter diagonalization,61 together with the direct-product basis described below. The AB c.m. position vector R in Eq. (2) is now expressed in terms of the spherical polar coordinates {R, Ω}, with R ≡ ||R||, and Ω ≡ (Θ, Φ), so that the complete set of coordinates is {R, Θ, Φ, θ, ϕ, r}. The matrix representation of the Hamiltonian is formed in a basis consisting of the product functions

|n,l,ml,j,m,rγ|n,l,ml|j,m|rγ,
(4)

where n = 0, 1, …, nmax, l = n, n − 2, …, ≥ 0, |ml| = 0, 1, …, l, j = 0, 1, …, jmax, |m| = 0, 1, …, j, and γ = 1, …, γmax. Here,

R,Θ,Φ|n,l,mlNnlRleβR2/2Lnl2(l+1/2)(βR)Ylml(Θ,Φ)
(5)

are eigenfunctions of the 3D isotropic HO (e.g., see the supplementary material in Ref. 49) having the angular frequency β/mAB, θ,ϕ|j,mYjm(θ,ϕ), and the |rγ⟩ constitute a discrete variable representation (DVR)62 derived from the eigenfunctions of a 1D oscillator of mass μAB moving in a Morse potential of the form

VMorse(r)=D1eα(rre)2,
(6)

where D, α, and re are constants chosen so that VMorse(r)Vh(1b)(r) in Eq. (1). The specific parameters that we have used in conjunction with this basis (all in atomic units) are as follows: D = 0.1744, α = 1.02764, re = 1.40201, and β = 2.9889. We note that because of the Pauli principle, j can be either even (para-H2) or odd (ortho-H2).

The Chebyshev variant of filter diagonalization requires the repeated application of Ĥ6D in Eq. (2) on an initial, random state vector. This is readily accomplished by matrix-vector multiplication for the kinetic-energy portion of Eq. (2). The ∇2 and ĵ2/r2 parts of Ĥ6D have analytic matrix elements in the basis of Eq. (4). The matrix elements of 2r2 are diagonal in all the basis-set indices except γ, and the rγ|2r2|rγ can be straightforwardly obtained by numerical transformation of the matrix elements from the Morse-eigenvector representation to the Morse-DVR one. To operate with the potential-energy portion of Eq. (2), we first transform the state vector to a grid representation |Rρ, (Θ,Φ)ξ, (θ,ϕ)η, rγ⟩, where the Rρ (ρ = 1, …, Nρ) are associated-Laguerre quadrature points and the (Θ, Φ)ξ (ξ = 1, …, Nξ) and the (θ, ϕ)η (η = 1, …, Nη) are Lebedev quadrature points. We then multiply the state vector at each grid point with the value of VH2-domain at that grid point. Finally, we transform the result back to the |n, l, ml, j, m, rγ⟩ representation. In this study, the grid sizes Nρ, Nξ, and Nη are 10, 110, and 38, respectively.

3. Obtaining converged first excited vibrational state of H2 and the isotopologues

As outlined in the Introduction, it has been generally assumed that a converged fully coupled quantum 6D calculation of the high-energy v = 1 vibrational state of the caged H2 would necessarily involve converging a very large number of the VTR states in the v = 0 manifold lying below it. This, of course, would require diagonalization of a prohibitively large matrix of the 6D VTR Hamiltonian, making the task virtually intractable.

However, convergence tests performed for the quantum 6D calculations utilizing the Smolyak scheme approach of the VTR levels of H2 inside the small dodecahedral cage with 20 water molecules reveal an entirely different picture. The calculations using different LS values for the basis set (LB = 6) and for the grid (LG = 7), which generate 8246 basis functions and 460 000 grid points, yielded 4120.9 cm−1 as the energy of the first excited (v = 1) vibrational state of H2 (for the TR DOFs in the ground state). The results shown in Tables III and IV are obtained utilizing this basis. Increasing LB to 7 and LG to 8 gives 17 900 basis functions and 1 167 282 grid points. Despite more than doubling the basis set size, the energies of the v = 1 vibrational state calculated for LB = 6 and LB = 7 differ by less than 0.1 cm−1, indicating its “high degree” of convergence. In contrast, comparison of the results of the two calculations shows that the highly excited TR states in the v = 0 manifold, close in energy to the v = 1 state, are not converged at all. In fact, only the v = 0 TR states with excitation energies up to 400–450 cm−1, far below the v = 1 state, are converged to within 1 cm−1.

Thus, what emerges from these calculations is the unexpected result that obtaining a well-converged energy of the v = 1 vibrational state of H2 does not require having converged high-lying v = 0 TR states in its vicinity and below. This suggests that the latter are very weakly coupled to the H2 stretch fundamental, and therefore, it confirms the validity of our previous 5D results.35 

Quantum 6D calculations on the same system using filter diagonalization and a direct-product basis described above confirm this finding and the conclusion, and go one step further. The basis set parameters (nmax, lmax, jmax, γmax) ranging from (7, 7, 4, 8) to (10, 6, 4, 11) give rise to basis sets ranging in size from 14 400 to 30 030. Although differing in size by more than a factor of two, these basis sets, when used in the quantum 6D calculations, all give the energies of the v = 1 vibrational state that are to within 0.1 cm−1 of each other and converge on 4121.1 cm−1. This result is very close to 4120.9 cm−1, the value obtained for the v = 1 state with Smolyak scheme approach.

A very interesting feature of the results of the filter-diagonalization calculations is that, owing primarily to the low jmax = 4, the computed TR states of the v = 0 manifold extend only up to excitation energies of 800 to 1500 cm−1, depending on the size of the basis, some 2500 cm−1below the energy of the v = 1 state. Thus, there are nov = 0 TR states close in energy to the v = 1 state. Moreover, only the TR states up to about 200 cm−1 above the ground state are well-converged. Despite that, the calculated energy of the H2 stretch fundamental is virtually identical, to within 0.2 cm−1, to that obtained by the Smolyak scheme approach, which does generate TR states of the v = 0 manifold, albeit not converged, in the neighborhood of the v = 1 state.

The surprising conclusion that emerges from the Smolyak scheme approach and filter-diagonalization calculations above is that the converged first excited state of H2 stretch can be obtained without (a) converging all TR states in the v = 0 manifold up to its energy, or (b) having any highly excited v = 0 TR states at all within a couple of thousands of wave numbers. This finding points out to extremely weak coupling between the v = 1 vibrational state of H2 and the high-lying v = 0 TR states. We do not have a formal theoretical explanation for this weak coupling at the present time. However, the disparity between the nodal patterns of the states involved, completely irregular for the highly excited TR states vs a smooth one, with a single node for the v = 1 state, are likely to figure prominently in any theoretical model.

Both the Smolyak scheme approach and filter-diagonalization calculations also demonstrate that in order to compute a highly converged H2 stretch fundamental, one (only) needs to use a basis for the intermolecular DOFs that can provide an accurate description of the vibrational averaging over the large-amplitude TR motions in the delocalized quantum ground state of the system. It should be stressed that this basis is much smaller than the one that would be required to get converged highly excited TR eigenstates in the vicinity of the v = 1 state. That accurate calculation of the vibrational shift in systems dominated by quantum effects must take into account averaging over the large-amplitude motions which was demonstrated first for the ArnHF clusters.36–39 

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

FIG. 1.

3D isosurface of the H2 c.m. probability distribution inside the small cage of the sII hydrate, from the 6D DMC simulation.

FIG. 1.

3D isosurface of the H2 c.m. probability distribution inside the small cage of the sII hydrate, from the 6D DMC simulation.

Close modal
TABLE II.

Ground-state energies E0 (in cm−1) of the H2 isotopologues inside the small cage of the sII hydrate comprising 20 water molecules, from the 6D DMC and Smolyak calculations.

MethodH2HDD2
DMC 1441 ± 10 1127 ± 9 770 ± 7 
Smolyak 1438.3 ± 0.1 1129.9 ± 0.1 771.2 ± 0.1 
MethodH2HDD2
DMC 1441 ± 10 1127 ± 9 770 ± 7 
Smolyak 1438.3 ± 0.1 1129.9 ± 0.1 771.2 ± 0.1 
TABLE III.

Comparison of the energies (in cm−1) of the fundamental translational excitations and the rotational j = 0 → 1 transitions of H2 in the ground (v = 0) vibrational state from quantum 6D calculations, in which the Smolyak scheme approach is employed, for three sII clathrate hydrate domains with the experimental results from Ref. 63 (in boldface). For a domain with N water molecules, where N = 20 corresponds to the isolated small cage, the calculated excitation energies are relative to the ground-state energy E0 of this domain from the quantum 6D calculations. Also shown are the quantum 6D frequencies ν (in cm−1) of the stretching fundamental (v = 1) of H2 in the three sII hydrate domains and their respective frequency shifts Δν. The experimental frequency shifts (in boldface) at 76 K and 20 K are from Refs. 27 and 26, respectively. The numbers in the brackets are the corresponding quantum 5D results reported in Ref. 35.

N
Expt.204076
E0 … 1438.3 1401.7 1382.6 
Translations     
71.0 66.8 (66.78) 66.5 (66.36) 66.2 (66.13) 
II 80.2 76.1 (76.02) 75.5 (75.35) 75.3 (75.12) 
III 101.1 93.3 (93.13) 92.5 (92.34) 92.2 (92.06) 
Rotations     
j = 1     
110.0 85.4 (86.61) 93.5 (94.68) 97.7 (98.95) 
II 116.5 121.2 (122.18) 121.1 (122.17) 118.4 (119.40) 
III 122.1 147.6 (148.65) 140.2 (141.39) 137.7 (138.85) 
ν  4120.9 4119.4 4119.2 
Δν −34 (76 K)/−37 (20 K) −40.2 (−39.81) −41.7 (−41.29) −41.9 (−42.30) 
N
Expt.204076
E0 … 1438.3 1401.7 1382.6 
Translations     
71.0 66.8 (66.78) 66.5 (66.36) 66.2 (66.13) 
II 80.2 76.1 (76.02) 75.5 (75.35) 75.3 (75.12) 
III 101.1 93.3 (93.13) 92.5 (92.34) 92.2 (92.06) 
Rotations     
j = 1     
110.0 85.4 (86.61) 93.5 (94.68) 97.7 (98.95) 
II 116.5 121.2 (122.18) 121.1 (122.17) 118.4 (119.40) 
III 122.1 147.6 (148.65) 140.2 (141.39) 137.7 (138.85) 
ν  4120.9 4119.4 4119.2 
Δν −34 (76 K)/−37 (20 K) −40.2 (−39.81) −41.7 (−41.29) −41.9 (−42.30) 

In Table III, we report the energies of the fundamental translational excitations and the rotational j = 0 → 1 transitions of H2 in the ground vibrational state (v = 0), as well as the frequency ν of the H2 stretch fundamental (v = 1) and the corresponding frequency shift Δν, for three sII hydrate domains with the number of water molecules N equal to 20 (isolated small cage), 40, and 76. These results are from the quantum 6D calculations employing the Smolyak scheme approach. Shown for comparison are the corresponding results from the quantum 5D calculation for the same hydrate domains in Ref. 35, in which the H2 molecule is treated as rigid, and the available experimental data pertaining to these quantities.26,27,63

Even a cursory inspection of Table III reveals a striking agreement between the results of the fully coupled quantum 6D calculations with those from the quantum 5D, rigid-H2 treatment, for all three domains. For the translational excitations, the agreement is better than 0.5 cm−1, while the 6D and 5D rotational excitations agree to about 1 cm−1.

Excellent agreement between the 6D and 5D calculations in Table III extends to the H2 vibrational frequency shift as well. In our quantum 6D calculations, the frequency shift Δν for a given hydrate domain is obtained as the difference between the frequency ν of the H2 stretch fundamental computed in 6D for this domain and the free-H2 stretch frequency νfree evaluated for the one-body potential Vh(1b)(r) in Eq. (1) (the νfree values for H2, HD, and D2 are given in Table IV). For the three hydrate domains considered, the difference between the frequency shifts computed in 6D and 5D is very small, less than 0.5 cm−1. This confirms the remarkably high accuracy of the quantum 5D method for computing vibrational frequency shifts used in Ref. 35 and earlier.36–39 Moreover, the fact that all results of the 6D and 5D calculations, translational and rotational excitations and frequency shifts, agree so exceedingly well points to a high degree of decoupling between the high-frequency H2 intramolecular stretch vibration and the low-frequency intermolecular TR modes.

TABLE IV.

Comparison of the energies (in cm−1) of the fundamental translational excitations and the rotational j = 0 → 1 transitions of H2, HD, and D2 in the ground (v = 0) vibrational state from quantum 6D calculations, in which the Smolyak scheme approach is employed, for the sII hydrate domain having 76 water molecules with the experimental results from Ref. 63 (in boldface). Also shown are the quantum 6D frequencies ν (in cm−1) of the stretching fundamentals (v = 1) of H2, HD, and D2 in the same sII hydrate domain and their respective frequency shifts Δν. The experimental frequency shifts (in boldface) for H2 and D2 at 76 K are from Ref. 27, while the shift for H2 at 20 K is from Ref. 26.

H2HDD2
Expt.TheoryExpt.TheoryTheory
Translations       
71.0 66.2 53.2 48.9  39.8 
II 80.2 75.3 58.7 54.7  46.4 
III 101.1 92.2 70.6 79.5  60.3 
Rotations       
j = 1       
110.0 97.7 87.6 62.6  40.6 
II 116.5 118.4 93.4 91.4  61.9 
III 122.1 137.7 98.5 106.9  77.7 
ν  4119.2  3595.5  2962.9 
νfree  4161.1  3631.8  2993.2 
Δν −34 (76 K)/−37 (20 K) −41.9  −36.3 −25 (76 K) −30.3 
H2HDD2
Expt.TheoryExpt.TheoryTheory
Translations       
71.0 66.2 53.2 48.9  39.8 
II 80.2 75.3 58.7 54.7  46.4 
III 101.1 92.2 70.6 79.5  60.3 
Rotations       
j = 1       
110.0 97.7 87.6 62.6  40.6 
II 116.5 118.4 93.4 91.4  61.9 
III 122.1 137.7 98.5 106.9  77.7 
ν  4119.2  3595.5  2962.9 
νfree  4161.1  3631.8  2993.2 
Δν −34 (76 K)/−37 (20 K) −41.9  −36.3 −25 (76 K) −30.3 

Table IV displays the results for H2, HD, and D2 encapsulated inside the small cage embedded within the sII hydrate domain having 76 water molecules, from the quantum 6D calculations employing the Smolyak scheme approach. The computational results shown for H2 and the isotopologues include the fundamental translational excitations and the rotational j = 0 → 1 transitions in the v = 0 state, and the fundamental stretch frequency ν, and its shift Δν. The pertinent experimental data are included for comparison. As explained in Sec. II A, the sII hydrate domain with 76 H2O molecules, which complete the first three hydration shells surrounding the hydrogen molecule,15 is chosen primarily because the magnitude of the H2 frequency shift computed previously (in 5D) for this domain35 is only 3% smaller than the shift calculated for the largest domain comprising 1945 H2 molecules, intended to mimic the bulk sII hydrate. Thus, the domain of this size captures virtually all of the condensed-phase effect on the vibrational frequency shift and on the TR excitations.15 

The splitting of both the translational fundamental and j = 0 → 1 transition for all three isotopologues into three components, evident in Table IV (and also in Table III), has been discussed previously,12,13,15,35 and is caused by the anisotropies, radial and angular, respectively, of the cage environment. For the isotopologues for which the experimental data are available, H2 and HD, the calculated splittings are, in general, substantially larger than the measured values, especially for the rotational j = 0 → 1 transition. The implication is that the pairwise-additive 6D PES employed overestimates the anisotropies of the H2-hydrate interaction, angular anisotropy, in particular. This has been attributed to nonadditive many-body interactions13 that are missing from this pairwise-additive PES. Therefore, it would be interesting to repeat these quantum 6D calculations for the H2-hydrate PES that, in addition to the 2-body H2–H2O interactions, would incorporate the 3-body H2–H2O–H2O interactions, such as those computed by Bowman and co-workers,41,48 and see what degree this would improve the agreement between the calculated and measured TR excitation energies.

The theoretical frequency shifts Δν reported in Table IV for H2, HD, and D2 are obtained as the difference between the frequency ν of the stretch fundamental calculated in 6D for each of the isotopologues inside the domain with 76 H2O molecules and its gas-phase value νfree (in Table IV) evaluated for the one-body potential Vh(1b)(r) in Eq. (1). The vibrational frequency shifts computed in this way are −41.9 cm−1 for H2, −36.3 cm−1 for HD, and −30.3 cm−1 for D2. As expected, the shifts decrease in magnitude with the increasing mass of the isotopologue. However, the ratio |Δν/ν| is virtually constant for H2, HD, and D2 and equal to 0.010. In other words, |Δν| ≈ 1% of ν for the three isotopologues.

As mentioned earlier, the frequency shift from the quantum 5D calculations for H2 in the small cage within the domain with 76 water molecules is ∼3% smaller in magnitude than the shift computed for a very large domain encompassing 1945 water molecules aimed to mimic the bulk sII hydrate.35 If we assume that the same relationship holds for the frequency shifts obtained with the quantum 6D calculations in this study, and for all three isotopologues, then the theoretical frequency shifts for H2, HD, and D2 in the small cage of sII hydrate are −43, −37, and −31 cm−1, respectively.

The vibrational frequency shift measured at 76 K for H2 in the small cage of sII hydrate27 (−34 cm−1) is about 21% smaller by magnitude than the (extrapolated) theoretical value of −43 cm−1, while the shift measured for D2, also at 76 K27 (−25 cm−1) is about 24% smaller in magnitude than the theoretical result of −31 cm−1. Experimental data regarding the frequency shift of HD in the sII hydrate are not available. The agreement between theory and experiment is satisfactory, given that the shifts are computed rigorously and from first-principles, with no adjustable parameters. The agreement improves if the temperature dependence of the frequency shift measured for H2 (but not HD and D2 so far) in the small sII hydrate cage26 is taken into account. As pointed out by Qu and Bowman,41 the frequency shift measured at 20 K is about −37 cm−1,26 compared to −34 cm−1 at 76 K.27 Since the quantum 6D results are for 0 K, it is more appropriate to compare our computed value for H2 (−43 cm−1) to the experimental shift at 20 K (−37 cm−1). In that case, the measured shift is only 14% smaller in magnitude than the theoretical result.

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

We performed fully coupled quantum 6D calculations of the vibration-translation-rotation (VTR) eigenstates of a flexible H2, HD, and D2 molecule inside the small cage of the sII clathrate hydrate, taken to be rigid. These calculations utilized two different approaches, the Smolyak scheme approach52–58 and the Chebyshev variant60 of filter diagonalization,61 together with the direct-product basis described in Ref. 49. It was demonstrated that with both approaches, it is entirely feasible to obtain a highly converged energy of the first excited (v = 1) intramolecular vibrational state of the caged diatomic molecule, and its frequency shift relative to the gas-phase value, without excessive computational effort. What made this possible was the realization that to obtain the converged intramolecular stretch fundamental of the entrapped H2 at ≈4100 cm−1 it sufficed to have converged only the TR states in the v = 0 manifold up to at most 400–450 cm−1 above the ground state, necessary for a proper description of the delocalized ground state of the system and the vibrational averaging over the large-amplitude TR motions. This led to the conclusion that the v = 1 intramolecular vibrational state is extremely weakly coupled to the highly excited v = 0 TR states.

Quantum 6D calculations of the coupled VTR eigenstates, including the v = 1 state and its frequency shift relative to the gas-phase value, were performed for H2, HD, and D2 encapsulated inside three spherical sII hydrate domains of increasing radius, treated as rigid. A pairwise-additive 6D intermolecular PES for H2 inside the hydrate domain was employed in these calculations, constructed using the ab initio-based40 6D H2–H2O pair potential, for flexible H2 and rigid H2O. In addition, the VTR ground state of H2 in the (rigid) small cage was determined by means of the 6D DMC simulations, to partly verify the correctness of the eigenstate-resolved calculations.

All results of the quantum 6D calculations for H2 in the three hydrate domains considered agree extremely well with those from the quantum 5D, rigid-H2 treatment,35 demonstrating the high accuracy of the quantum 5D method for computing vibrational frequency shifts employed in Ref. 35 and earlier applications.36–39 

Comparison of the quantum 6D frequency shifts for H2 and D2 with the corresponding experimental results at 76 K27 shows that the latter are 21% and 24% smaller in magnitude, respectively. The difference in the magnitudes of the calculated H2 frequency shift and that measured for H2 at 20 K26 is only 14%. The agreement between theory and experiment is satisfactory, but clearly there is room for improvement.

The quantum 6D calculation of the vibrational frequency shift is rigorous and yields results that are virtually exact numerically for the PES employed. The only remaining dynamical approximation is treating the hydrogen-bonded water framework as rigid. However, the H2-hydrate interaction is weak, and moreover, the disparity between the masses of H2 (and isotopologues) and the confining hydrate is large. As a result, the coupling of VTR motions of H2 to the vibrations (phonons) of the host water framework is weak as well, and its neglect (by treating the hydrate as rigid) is not expected to introduce significant errors in the calculated frequency shifts. Consequently, the main source of the residual differences between the computed and experimental values have to be certain shortcomings in the pairwise-additive intermolecular PES for the H2-hydrate interaction, having to do with either the ab initio 6D H2–H2O pair potential or the absence of the nonadditive three-body interactions, or a combination of both. These possibilities will be investigated in the future.

It is likely that very weak coupling between the high-frequency intramolecular modes and the low-frequency intermolecular vibrations is the feature of other molecular systems, in particular, the weakly bound ones, e.g., hydrogen-bonded and van der Waals complexes mentioned in the Introduction. In that case, the fundamental excitations of their intramolecular modes, and frequency shifts, could be calculated accurately from full-dimensional, fully coupled quantum bound-state calculations, without converging the very large number of highly excited intermolecular vibrational eigenstates in the manifold of the intramolecular ground state. This could be achieved with a relatively small basis for the intermolecular DOFs capable of accurately describing the vibrational averaging over the large-amplitude intermolecular motions in the delocalized ground-state wave function, but not the high-lying intermolecular eigenstates. One attractive and challenging target for such a treatment is the HF dimer, for which high-quality full-dimensional (6D) PESs are available,44,64–66 as is the wealth of spectroscopic data about its intra- and intermolecular excitations.

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

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.
F.
Sebastianelli
,
M.
Xu
, and
Z.
Bačić
,
J. Chem. Phys.
129
,
244706
(
2008
).
21.
A.
Witt
,
F.
Sebastianelli
,
M. E.
Tuckerman
, and
Z.
Bačić
,
J. Phys. Chem. C
114
,
20775
(
2010
).
22.
A.
Valdeś
and
G. J.
Kroes
,
J. Phys. Chem. C
116
,
21664
(
2012
).
23.
P. M.
Felker
,
J. Chem. Phys.
138
,
174306
(
2013
).
24.
P. M.
Felker
,
J. Chem. Phys.
141
,
184305
(
2014
).
25.
D. M.
Benoit
,
D.
Lauvergnat
, and
Y.
Scribano
,
Faraday Discuss.
212
,
533
(
2018
).
26.
A.
Giannasi
,
M.
Celli
,
L.
Ulivi
, and
M.
Zoppi
,
J. Chem. Phys.
129
,
084705
(
2008
).
27.
T. A.
Strobel
,
E. D.
Sloan
, and
C. A.
Koh
,
J. Chem. Phys.
130
,
014506
(
2009
).
28.
S. A.
FitzGerald
,
K.
Allen
,
P.
Landerman
,
J.
Hopkins
,
J.
Matters
,
R.
Myers
, and
J. L. C.
Rowsell
,
Phys. Rev. B
77
,
224301
(
2008
).
29.
S. A.
FitzGerald
,
J.
Hopkins
,
B.
Burkholder
,
M.
Friedman
, and
J. L. C.
Rowsell
,
Phys. Rev. B
81
,
104305
(
2010
).
30.
J.
Wang
,
H.
Lu
, and
J. A.
Ripmeester
,
J. Am. Chem. Soc.
131
,
14132
(
2010
).
31.
J.
Wang
,
H.
Lu
,
J. A.
Ripmeester
, and
U.
Becker
,
J. Phys. Chem. C
114
,
21042
(
2010
).
32.
K. R.
Ramya
and
A.
Venkatnathan
,
J. Chem. Phys.
138
,
124305
(
2013
).
33.
N.
Plattner
and
M.
Meuwly
,
J. Chem. Phys.
140
,
024311
(
2014
).
34.
Z.
Futera
,
M.
Celli
,
L.
del Rosso
,
C. J.
Burnham
,
L.
Ulivi
, and
N. J.
English
,
J. Phys. Chem. C
121
,
3690
(
2017
).
35.
A.
Powers
,
Y.
Scribano
,
D.
Lauvergnat
,
E.
Mebe
,
D.
Benoit
, and
Z.
Bačić
,
J. Chem. Phys.
148
,
144304
(
2018
).
36.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
6359
(
1994
).
37.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
101
,
10181
(
1994
).
38.
S.
Liu
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
J. Chem. Phys.
103
,
1829
(
1995
).
39.
Z.
Bačić
,
J. Chem. Soc., Faraday Trans.
93
,
1459
(
1997
).
40.
P.
Valiron
,
M.
Wernli
,
A.
Faure
,
L.
Wiesenfeld
,
C.
Rist
,
S.
Kedzuch
, and
J.
Noga
,
J. Chem. Phys.
129
,
134306
(
2008
).
41.
C.
Qu
and
J. M.
Bowman
,
J. Phys. Chem. A
123
,
329
(
2019
).
42.
A.
Sarsa
,
Z.
Bačić
,
J. W.
Moskowitz
, and
K. E.
Schmidt
,
Phys. Rev. Lett.
88
,
123401
(
2002
).
43.
H. S.
Lee
,
J. M.
Herbert
, and
A. B.
McCoy
,
J. Chem. Phys.
110
,
5481
(
1999
).
44.
D. H.
Zhang
,
Q.
Wu
,
J. Z. H.
Zhang
,
M.
von Dirke
, and
Z.
Bačić
,
J. Chem. Phys.
102
,
2315
(
1995
).
45.
Y.
Qiu
and
Z.
Bačić
,
J. Chem. Phys.
106
,
2158
(
1997
).
46.
A.
Bahel
and
Z.
Bačić
,
J. Chem. Phys.
111
,
11164
(
1999
).
47.
Y.
Qiu
,
J. Z. H.
Zhang
, and
Z.
Bačić
,
J. Chem. Phys.
108
,
4804
(
1998
).
48.
Z.
Homayoon
,
R.
Conte
,
C.
Qu
, and
J. M.
Bowman
,
J. Chem. Phys.
143
,
084302
(
2015
).
49.
P. M.
Felker
and
Z.
Bačić
,
J. Chem. Phys.
145
,
084310
(
2016
).
50.
D. M.
Benoit
and
D. C.
Clary
,
J. Chem. Phys.
113
,
5193
(
2000
).
51.
D. M.
Benoit
, XDMC, Quantum Diffusion Monte Carlo Code,
1999
.
52.
S. A.
Smolyak
,
Sov. Math. Dokl.
4
,
240
(
1963
).
53.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
131
,
174103
(
2009
).
54.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
134
,
054126
(
2011
).
55.
G.
Avila
and
T.
Carrington
,
J. Chem. Phys.
135
,
064101
(
2011
).
56.
D.
Lauvergnat
and
A.
Nauts
,
Phys. Chem. Chem. Phys.
12
,
8405
(
2010
).
57.
D.
Lauvergnat
and
A.
Nauts
,
Spectrochim. Acta. Part A
119
,
18
(
2014
).
58.
J.
Sarka
,
D.
Lauvergnat
,
V.
Brites
,
A. G.
Császár
, and
C.
Léonard
,
Phys. Chem. Chem. Phys.
18
,
17678
(
2016
).
59.
D.
Lauvergnat
, ElVibRot: Quantum Dynamics Code, http://pagesperso.lcp.u-psud.fr/lauvergnat/ElVibRot/ElVibRot.html,
2006
.
60.
V. A.
Mandelshtam
and
H. S.
Taylor
,
J. Chem. Phys.
106
,
5085
(
1997
).
61.
M. R.
Wall
and
D.
Neuhauser
,
J. Chem. Phys.
102
,
8011
(
1995
).
62.
H.
Wei
and
T.
Carrington
, Jr.
,
J. Chem. Phys.
97
,
3029
(
1992
).
63.
L.
Ulivi
,
M.
Celli
,
A.
Giannasi
,
A. J.
Ramirez-Cuesta
,
D. J.
Bull
, and
M.
Zoppi
,
Phys. Rev. B
76
,
161401(R)
(
2007
).
64.
M.
Quack
and
M. A.
Suhm
,
J. Chem. Phys.
95
,
28
(
1991
).
65.
Z.
Bačić
and
Y.
Qiu
, in
Advances in Molecular Vibrations and Collision Dynamics
, edited by
J. M.
Bowman
and
Z.
Bačić
(
JAI Press, Inc.
,
Stamford
,
1998
), Vol. 3, p.
183
.
66.
G. W. M.
Vissers
,
G. C.
Groenenboom
, and
A.
van der Avoird
,
J. Chem. Phys.
119
,
277
(
2003
).