Path Integral Molecular Dynamics Study of Small H2 Clusters in the

Oct 12, 2010 - Chemistry and Courant Institute of Mathematical Sciences, New York ... Institute of Theoretical and Computational Science, East China N...
0 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. C 2010, 114, 20775–20782

20775

Path Integral Molecular Dynamics Study of Small H2 Clusters in the Large Cage of Structure II Clathrate Hydrate: Temperature Dependence of Quantum Spatial Distributions† Alexander Witt,‡ Francesco Sebastianelli,§ Mark E. Tuckerman,*,§,| and Zlatko Bacˇic´*,§,⊥ Lehrstuhl fu¨r Theoretische Chemie, Ruhr-UniVersita¨t Bochum, D-44780 Bochum, Germany, Department of Chemistry and Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York 10003, United States, and State Key Laboratory of Precision Spectroscopy and Department of Physics, Institute of Theoretical and Computational Science, East China Normal UniVersity, Shanghai 200062, China ReceiVed: July 27, 2010; ReVised Manuscript ReceiVed: September 27, 2010

We report a path integral molecular dynamics (PIMD) study of the temperature dependence of the spatial distribution of two and four H2 molecules inside the large cage of the structure II clathrate hydrate. The PIMD calculations were performed at five temperatures ranging from 25 to 200 K. Their results were combined with those from an earlier diffusion Monte Carlo (DMC) study of this system at T ) 0 K [Sebastianelli, F.; Xu, M.; Bacˇic´, Z. J. Chem. Phys. 2008, 129, 244706]. The spatial distribution of the confined H2 molecules at each of the temperatures considered was characterized with the help of several one-dimensional (1D) and three-dimensional (3D) distribution functions of suitably chosen intermolecular coordinates, generated by the PIMD and DMC calculations. The 1D distribution that proved to be the most strongly temperature dependent, and also the most revealing about the structural properties of the system as a function of temperature, involves the H2-cage center-H2 angle. In the case of four caged H2 molecules, this angular distribution provides clear evidence that between 50 and ∼100 K the system undergoes a qualitative change. At 50 K and below, the system is fully localized in the global minimum of the intermolecular potential, corresponding to a tetrahedral configuration of H2 molecules with a unique orientation relative to the cage frame. At temperatures of 75-100 K and higher, nearly degenerate local minima ∼200 cm-1 above the global minimum become accessible and are increasingly sampled by the system. The 3D spatial distributions also show this growing delocalization above 75-100 K. Our findings are in accord with the localization-delocalization transition observed experimentally to occur at 50 K for four D2 molecules in the large cage [Lokshin, K. A.; Zhao, Y.; He, D.; Mao, W. L.; Mao, H. K.; Hemley, R. J.; Lobanov, M. V.; Greenblatt, M. Phys. ReV. Lett. 2004, 93, 125503]. I. Introduction Clathrate hydrates are crystalline solids where guest molecules of different sizes are encapsulated inside close-packed polyhedral cavities of the three-dimensional (3D) host lattice formed from hydrogen-bonded water molecules.1-3 There are three common clathrate hydrate structures, and which of them is formed depends strongly on the size of the guest molecule(s) essential for the stability of the clathrates.1-3 For a long time, the prevailing opinion was that hydrogen molecules were too small, and interacted too weakly with the water molecules, to stabilize the clathrate hydrate framework. However, this widely held view had to be revised in the past decade, when it was demonstrated that molecular hydrogen does form clathrate hydrates.4,5 Simple hydrogen hydrates, having the hydrogen molecules as the sole guests, adopt the classical structure II (sII),1,2 whose unit cell is cubic, with 136 water molecules organized in the hydrogen-bonded framework comprised of two types of cages: 16 dodecahedral (512), or “small” cages, with 12 pentagonal faces, and 8 hexakaidecahedral (51264), or “large” †

Part of the “Mark A. Ratner Festschrift”. * To whom correspondence should be addressed. E-mail: mark.tuckerman@ nyu.edu (M.E.T.), [email protected] (Z.B.). ‡ Ruhr-Universita¨t Bochum. § Department of Chemistry, New York University. | Courant Institute of Mathematical Sciences, New York University. ⊥ East China Normal University.

cages, having 12 pentagonal and 4 hexagonal faces. Initially, based on the Raman spectroscopic measurements, it was estimated that each small cage is occupied by two H2 molecules, whereas four H2 molecules occupy the large cage,5 equivalent to 5.3 mass % of molecular hydrogen in the material, which is close to the US DOE 2010 target value of 5.5 mass %. This generated a great deal of interest in hydrogen hydrates as potential hydrogen storage materials.1,2,6-9 A medium for hydrogen storage made of pure water would be environmentally very friendly, safe, and economical. Following the work of Mao et al.,5 Lokshin et al.10 performed neutron diffraction studies of the pure D2 clathrate hydrate, confirming its sII crystal structure. However, unlike Mao et al.,5 they found only one D2 molecule in the small cages over the entire range of temperatures and pressures tested. In the large cage, Lokshin et al.10 observed four D2 molecules at 200 MPa and up to ∼180 K, confirming the estimate of Mao et al.5 But at ambient pressure, as the sample was heated above ∼80 K, the number of D2 molecules in the large cavity gradually decreased from four to two, until the decomposition temperature of ∼160 K was reached.10 With only single hydrogen molecule occupancy of the small cage, and the maximum quadruple occupancy of the large cage, the hydrogen storage capacity is reduced to ∼3.8 mass %. The main obstacle to the practical applications of simple H2 hydrates are the rather extreme conditions, high pressure and low temperature, required for their formation, as well as the low temperatures necessary

10.1021/jp107021t  2010 American Chemical Society Published on Web 10/12/2010

20776

J. Phys. Chem. C, Vol. 114, No. 48, 2010

to store them. It has been shown that the addition of a second, larger promoter molecule such as tetrahydrofuran (THF) reduces drastically the formation pressure from 200 to 5 MPa at 280 K.11,12 However, this is achieved at the expense of the hydrogen storage capacity, since in the binary THF + H2 hydrate all of the large cages are filled with the THF molecules, leaving only the small clathrate cages, which only hold one H2 per cage, available for the hydrogen molecules. Bringing the conditions under which hydrogen hydrates are formed closer to ambient conditions without sacrificing too much of their capacity for hydrogen storage remains a major challenge in this field. In addition to their potential as hydrogen storage materials, hydrogen hydrates provide an exceptional opportunity for investigating the dynamical and spectroscopic manifestations of the trapping of one or more light vib-rotors inside nanocavities of different sizes and shapes. The confinement leads to the discretization of the translational degrees of freedom of the guest molecules, in addition to the already quantized rotational states. Moreover, inside the cavity, the translational “rattling” motion of the center of mass (cm) of the caged molecule and its overall rotation are strongly coupled. Both the discrete translational and rotational eigenstates are well separated in energy, due to the small mass and large rotational constant of the hydrogen molecule, as well as the nanoscale size of the confining spaces. Consequently, the translation-rotation (T-R) energy level structure is sparse. It is even sparser in the case of the homonuclear isotopologues H2 and D2, where the symmetry constraints on their total wave functions result in the existence of two quite distinct species: para-H2 (p-H2), with the total nuclear spin I ) 0, and ortho-D2 (o-D2) with I ) 0 or 2, having even-j rotational states only (j ) 0, 2, 4, ...), and ortho-H2 (oH2), and para-D2 (p-D2), both with I ) 1, with exclusively odd-j rotational states (j ) 1, 3, 5, ...). Finally, the zero-point energy (ZPE) of the coupled T-R motions is substantial relative to the well depth of the interaction potential. As a result of these combined quantum effects, hydrogen molecules encapsulated in the clathrate cages constitute a highly quantum mechanical system, particularly at the low temperatures at which the hydrogen hydrates are synthesized and at which the experiments on them are typically carried out. Therefore, quantitative theoretical predictions regarding the energetics, spectroscopy, and dynamics of hydrogen hydrates at low temperatures, the number of trapped hydrogen molecules, and their spatial distribution inside the cavities, demand a fully quantum mechanical treatment, solving the multidimensional Schro¨dinger equation for the coupled T-R motions of the guest hydrogen molecule(s). In a series of recent papers, we have reported the first rigorous investigations of the quantum T-R dynamics of hydrogen molecules encapsulated in the small and large cages of the sII clathrate hydrate.13-17 In our theoretical approach, the three translational and the two rotational degrees of freedom of a confined hydrogen molecule are treated explicitly, as fully coupled, without any dynamical approximation, while the clathrate cages are assumed to be rigid, and their geometries determined in the X-ray diffraction experiments18 are employed. For a single confined H2 (HD, D2) molecule, our studies have revealed the salient features of its T-R dynamics, in particular the subtle connections between the symmetries of the small and large cages on one hand, and the patterns of the translationally excited energy levels and the quantum numbers required for their assignment on the other. When only the heavy oxygen atoms of the framework water molecules are taken into account, the three principal moments of inertia of the large cage are

Witt et al. identical, while for the small cage two principal moments of inertia are the same and the third one is larger. Therefore, the large cage is essentially a spherical top, and the small cage has the lower symmetry of a (slightly) oblate symmetric top. This difference in the symmetries of the two cages leaves unmistakable fingerprints in the translational energy level structure. In the large cage, the triple degeneracy of the fundamental translational excitation of the caged H2 (HD, D2) molecule is left virtually intact, evidence of the high symmetry of the confining environment.17 In contrast, inside the small cage, the H2 (HD, D2) translational fundamental is split into the (almost) doubly degenerate in-plane (xy) excitation and a nondegenerate perpendicular (z) excitation,15 as expected for the symmetrictop environment whose two axes are symmetrically equivalent and one is distinct. Moreover, inspection of the patterns of the translational excitations and the plots of the translational components of the wave functions showed that in the large cage they can be assigned using the quantum numbers n and l of the 3D isotropic harmonic oscillator (HO),17 whereas in the small cage the appropriate quantum numbers are V and l of the 2D isotropic HO for the in-plane (xy) modes and the Cartesian quantum number Vz for the perpendicular z-mode excitations.15 In both the small and the large cavities, the 2j + 1 degeneracy of the j ) 1 and j ) 2 levels is lifted due to the angular anisotropy of the H2-cage potential energy surface (PES).15,17 The splittings of the j ) 1 and j ) 2 multiplets in the two cages are very similar with respect to their patterns, magnitudes, and the energies of the individual components. The j ) 1 triplet is split into three distinct levels, and the j ) 2 quintuplet has one level in the center of the pattern and the other four appear as two pairs of closely spaced levels separated by ∼1 cm-1, with one pair lying energetically above and the other below the central component.15,17 For a hydrogen molecule in the small cage, our results15 regarding the patterns and magnitudes of the splittings of the translational fundamental and the j ) 1 triplet of H2 and HD, as well as the energy of the translational overtone excitation, are in near-quantitative agreement with the inelastic neutron scattering (INS) studies of THF + H2 and THF + HD hydrates.19,20 The distinctive splitting pattern of the j ) 2 quintuplet of H2, HD, and D2 which we predicted15 has been observed in the rotational Raman spectra of THF + H2,21,22 THF + HD,21 and THF + D2 hydrates.22 Our finding17 that the j ) 1 and j ) 2 levels of the hydrogen molecules exhibit almost identical splittings in the small and large cages is consistent with the observation that the S0(0) (j ) 0 f 2) bands in the rotational Raman spectra measured for simple H2 hydrate and the binary THF + H2 hydrate show striking similarity in their frequencies, widths, shapes, and internal structure, for low H2 occupancy of the large cage of simple H2 hydrate.22 The dimensionality of the quantum bound-state calculations of the T-R eigenstates is 5n, where n is the number of H2 molecules confined inside the (rigid) nanocavities, when the H2 bond length is held fixed; otherwise it is 6n. The 5D calculations for a single encapsulated molecule are readily performed. However, obtaining the excited T-R energy levels of two confined H2 molecules already requires quantum 10D calculations, which are exceedingly difficult and demanding in terms of both computation time and memory requirements. The methodology for solving this most challenging problem is under development in our group. In the meantime, we have performed diffusion Monte Carlo (DMC) calculations of the ground-state properties, energetics and the vibrationally averaged geometries, of one and two p-H2 molecules in the small cage,15 and of up

Temperature Dependence of Quantum Spatial Distributions to five p-H2 and o-D2 molecules inside the large cage.16 Although limited to the ground state of the systems considered, the DMC calculations produced valuable results which could be directly compared with the experimental data. Thus, in the case of the small cage, for one caged H2 molecule, the computed groundstate energy is negative, implying that it is truly bound and stable relative to the hydrogen molecule at a large distance outside the cage (the zero of energy). When the second molecule is placed inside the small cage, the ground-state energy takes a high positive value, essentially ruling out the double occupancy of H2 as energetically highly unfavorable.15 This agrees with the recent experiments showing only one H2 molecule in the small cage.10,23 For the large cage, the DMC ground-state energies of (pH2)n and (o-D2)n clusters are negative and show a weak dependence on n in the range n ) 1-4. However, upon addition of the fifth p-H2 or o-D2 molecule the ground-state energy of the cluster sharply rises to high positive values. This means that up to four H2 or D2 molecules are stabilized by the encapsulation in the large cage. The energetic penalty for adding the fifth hydrogen molecule is very high, making the quintuple H2/D2 occupancy of the large cage extremely improbable.16 The predicted maximum large cage occupancy of four H2 or D2 molecules is in agreement with the results of the neutron diffraction experiments10 for temperatures below 70 K under ambient pressure. For n ) 4, the DMC-calculated averages of the D2-D2 separation and the distance between the D2 molecules and the cage center agree with the experimental results10 to better than 10%. In addition, the vibrationally averaged tetrahedral geometry of the (o-D2)4 cluster from the DMC calculations16 reproduced the experimental finding10 that below 50 K each vertex of the (D2)4 tetrahedron is oriented toward the center of one of the four hexagonal faces of the large cage. The DMC calculations, although very informative, provide a quantum treatment of the properties of nanoconfined molecular hydrogen at 0 K only. In this paper, path integral molecular dynamics is employed to calculate several different quantum spatial distributions, both 1D and 3D, of two and four hydrogen molecules in the large cage, for temperatures ranging from 25 to 200 K. These calculations yield a quantitative description of the temperature dependence of the spatial distribution of small clusters of H2 molecules inside the large cage. The main motivation for this work comes from the already mentioned neutron diffraction study by Lokshin et al.,10 in which they also reported a localization-delocalization transition for the four D2 molecules residing in the large cage. As described above, below 50 K the four tetrahedrally arranged D2 molecules are localized in a unique orientation relative to the cage framework, in which they point at the centers of the hexagonal faces of the large cage. As the temperature is increased above 50 K, D2 molecules become increasingly delocalized, approaching a nearly uniform distribution on a spherical shell at 200 K.10 We have qualitatively explained this localization-delocalization transition in terms of certain features of the potential energy landscape of the (H2)4 cluster in the large cage.16 The path integral molecular dynamics calculations reported in this paper provide the much needed quantitative treatment of this interesting phenomenon. II. Computational Details The aim of a path integral molecular dynamics (PIMD) calculation is to generate the equilibrium properties of a system at temperature T described by a canonical density matrix Fˆ ∝ ˆ is the Hamiltonian. The ˆ ), where β ) 1/kBT, and H exp(-βH partition function Q(N, β) of the system is just the trace of the

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20777 ˆ )]. (Since the present density matrix: Q(N, β) ) Tr[exp(-βH study concerns confined systems, Q has no explicit volume dependence.) The Feynman path integral formulation of equilibrium quantum statistical mechanics is derived by carrying out the trace in the basis of coordinate eigenstates. Given a system of N quantum particles with coordinate operators ˆ N obeying ordinary Boltzmann statistics, the partition ˆ 1, ..., R R function can be expressed as ˆ

Q(N, β) ) Tr[e-βH] )

∫dNR〈R1 · · · RN|e-βH|R1 · · · RN〉 ˆ

(1) ˆ ) Kˆ + Assuming the Hamiltonian has the standard form H ˆ , where Kˆ and U ˆ are the usual kinetic and potential energy U operators, respectively, the Trotter theorem is used to factorize ˆ +U ˆ ] ) limPf∞[exp(-βK ˆ/ the Boltzmann operator as exp[-β(K ˆ /P)]P. This is then substituted into eq 1, and P) exp(-βU complete sets of coordinate and momentum eigenstates are inserted between each of the factors of exp(-βKˆ/P) and ˆ /P). The resulting momentum integrations can be exp(-βU performed analytically, and the path integral expression ultimately becomes

[∏( ) N

Q(N, β) ) lim

Pf∞

{ ∑[∑ P

N

s)1

I)1

exp -β

I)1

MIP

2πβp

3P/2

2

∫dR

(1) (P) I · · · dRI

]

×

]}

1 1 M ω2 (R(s+1) - RI(s))2 + U(R1(s), ..., R(s) N) 2 I P I P

(2)

where MI, I ) 1, ..., N are the particle masses, ωP ) P1/2/βp, and RI(s), s ) 1, ..., P indexes the imaginary time path for the Ith particle. The trace requires that the paths satisfy RI(P+1) ) RI(1). As was shown by Wolynes and Chandler,24 eq 2 is isomorphic to the classical partition function of a set of N cyclic polymer chains, each containing P pseudoparticles, usually called “beads”, with a harmonic nearest-neighboring coupling in each chain and interacting via the N-particle potential (1/ (s) P)U(R(s) 1 , ..., RN ). The potential can be seen to act between beads on each chain having the same imaginary time index s. In principle, the multidimensional integral in eq 2 could be sampled, for a finite value of P, in a straightforward Monte Carlo scheme or via molecular dynamics with a fictitious Hamiltonian of the form P

HP )

{ [ N

∑ ∑ s)1

I)1

2 (P(s) I ) + 2MI

N

2 - R(s) ∑ 21 MIωP2 (R(s+1) I I ) I)1

] }

1 (s) U(R(s) 1 , ..., RN ) P

+

(3)

however, such straightforward schemes are known to suffer from severe ergodicity problems in the limit of large P. In particular, Hall and Berne25 showed that molecular dynamics trajectories generated with eq 3 remain close to the invariant tori of the stiff harmonic coupling term. Tuckerman et al.26 showed that this problem could be solved by introducing a change of variables into eq 2 of the form

20778

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Witt et al.

(1) u(1) I ) RI

(s) u(s) I ) RI -

(s - 1)R(s+1) - R(1) I I , s ) 2, ..., P s

(4)

referred to as a staging transformation.26,27 This transformation has the effect of diagonalizing the harmonic coupling term: P

P

s)1

s)1

1 )2 ) ∑ M(s) ω2 (u(s))2 ∑ 21 MIωP2 (R(s)I - R(s+1) I 2 I P I

(5)

where MI(1) ) 0 and MI(s) ) sMI/(s - 1), s ) 2, ..., P. Thus, in order to ensure that each staging variable moves on the same time scale, a molecular dynamics scheme is constructed using the fictitious Hamiltonian P

HP )

{ [ N

∑ ∑ s)1

I)1

2 (P(s) I )

˜ (s) 2M I

N

+

]

∑ 21 M(s)I ωP2 (u(s)I )2 I)1

+

}

1 (s) U(R(s) 1 (u1), ..., RN (uN)) P

(6)

˜ I(1) ) MI and M ˜ I(s) ) MI(s), ˜ I(s) are M where the ficitious masses M (s) and RI (uI) denotes the inverse of the staging transformation. This choice of fictitious masses ensures that the fundamental frequency in the now uncoupled harmonic term is the same for each staging mode. When eq 6 is used, molecular dynamics trajectories are generated in terms of the staging variables, and in order to ensure optimal sampling, a Nose´-Hoover chain thermostat28 is coupled to each Cartesian component of each staging variable. More details of path integral molecular dynamics calculations are given in ref 29. It was also argued in ref 26 and elsewhere30,31 that PIMD simulations can be performed in terms of normal mode variables as well, although they are not employed here. In the present study, PIMD calculations were performed in staging variables for two and four molecules in the large cage at temperatures of 25, 50, 75, 100, and 200 K. Each H2 molecule is composed of two quantum hydrogen atoms discretized into cyclic paths containing 256, 128, 128, 128, and 32 beads, at 25, 50, 75, 100, and 200 K, respectively. In these simulations, water molecules are kept frozen and, therefore, act as an external potential on the H2 molecules. The positions of the O and H atoms of the framework water molecules were taken from the earlier DMC study of hydrogen molecules in the large cage.16 Our past investigations of one or more hydrogen molecules in the small and large cages of the sII clathrate hydrate13-17 have demonstrated that treating the cages as rigid captures accurately the key features of the quantum T-R dynamics of the guest molecules, and yields results that are in near-quantitative agreement with a variety of recent spectroscopic measurements. The H atoms of the water molecules forming the cage are configurationally disordered. The number of distinct hydrogenbonding (H-B) arrangements exceeds 30 000 already for the small cage32 and is undoubtedly much greater for the large cage. Quantum computations in our previous studies have shown that the main results depend very weakly on the particular H-B topology.13-17 As before,16 the interactions between the guest hydrogen molecules and the water nanocage, as well as those between the confined hydrogen molecules, were assumed to be pairwise

additive. The anisotropic H2-H2O and H2-H2 pair potentials employed were identical to those used in the DMC calculations of molecular hydrogen clusters in the large cage.16 The H2-H2O pair interaction is due to Alavi et al.,33 and the H2-H2 pair potential is described by the ab initio 4D PES by Diep and Johnson.34 Detailed description of the total intermolecular PES for n H2 molecules in the large cage can be found in ref.16 In this PES the H2 molecules are treated as rigid, which complicates the path-integral treatment. Therefore, in order to circumvent this technicality, we introduced a harmonic bond between the hydrogen atoms of H2 with a spring coupling of 4.177 × 105 K/Å2 and equilibrium bond length of 0.74 Å. We tested the value of the spring constant against runs in which the spring constant was doubled and found the results to be insensitive to an increase in the bond stiffness. Each path integral trajectory was equilibrated for 10 ps followed by production runs of 2, 1.6, 1.6, 1.5, and 1.0 ns at 200, 100, 75, 50, and 25 K, respectively, for four H2 molecules in the cage and 2.0, 1.6, 1.8, 1.9, and 1.2 ns at 200, 100, 75, 50, and 25 K, respectively, for two H2 molecules in the cage. All simulations are carried using a time step of 0.125 fs using the PINY_MD code.35 III. Results and Discussion In presenting the results, we will focus on three 1D distribution functions in particular. These are P(R), where R is the distance between the cage center (cc) and the center of mass (cm) of an H2 molecule, P(r), where r is the distance between the centers of mass of pairs of H2 molecules, and P(Θ), where Θ is the cm1-cc-cm2 angle, cm1 and cm2 being the centers of mass of the individual H2 molecules in a pair. For each of these distributions, the raw histograms generated from the PIMD trajectories and the DMC were divided by appropriate Jacobian factors [R2 for P(R), r2 for P(r), and sin(Θ) for P(Θ)] to produce the final probability distributions. The distributions are then normalized such that ∫0∞ R2P(R) dR ) 1, ∫0∞ r2P(r) dr ) 1, and ∫0π sin(Θ)P(Θ) dΘ ) 1. For additional visual insight, we also present full 3D spatial distributions. A. Two H2 Molecules in the Large Cage. Figure 1 shows the distributions P(R) (top), P(r) (middle), and P(Θ) (bottom) for two H2 molecules in the large cage generated from PIMD calculations at the five chosen temperatures in the range 25-200 K, and from DMC calculations at T ) 0 K (ref 16). It is evident that at all temperatures except the highest (200 K), P(R) is virtually insensitive to temperature and peaks around 3 au. Mutual repulsion keeps the two H2 molecules outside the central region of the cage with the radius R of about 1.5 au, and they reside within a shell whose width is ∼2 au. Only the 200 K distribution deviates perceptibly from the T ) 0 K DMC result. At the highest temperature, P(R) peaks at a slighly larger R value, and the distribution is somewhat broader. This provides evidence that thermal effects cause the hydrogen molecules to penetrate further in the directions of both the cage center and the water molecules of the cage wall. The peak of the distribution P(r) shown in Figure 1 (middle) gradually moves to smaller r values with increasing temperature, indicating that thermal excitation allows the H2 molecules to come closer to each other. As the temperature is decreased, the P(r) distributions approach that for T ) 0 K obtained using the DMC method.16 The most interesting of the three distributions is the angular distribution P(Θ) displayed in Figure 1 (bottom), which exhibits a significant temperature dependence. At all temperatures its maximum is at 180°, corresponding to a collinear geometry of the two H2 molecules and the cage center. For higher temper-

Temperature Dependence of Quantum Spatial Distributions

Figure 1. Distribution functions of the H2-cage center distance (top), H2-H2 distance (middle), and H2-cc-H2 angle (bottom) for two H2 molecules in the large cage generated from path integral molecular dynamics at the five chosen temperatures, and from diffusion Monte Carlo (DMC) at T ) 0 K.16

Figure 2. The H2-cage center-H2 angle distribution sin (Θ)P(Θ) for two H2 molecules in the large cage from path integral molecular dynamics at the five chosen temperatures, and from diffusion Monte Carlo (DMC) at T ) 0.16 The inclusion of the Jacobian factor highlights the peak of the T ) 0 K distribution at 157°, close to the global minimum geometry for which the angle is Θ ) 156.8°.

atures, the distribution becomes wider and extends to smaller angles Θ, implying that strongly bent H2-cc-H2 configurations become increasingly accessible; in particular, the 200 K distribution, although small, is nonzero for angles as small as 90°. Further insights into the vibrationally averaged geometry of two H2 molecules are gained by plotting the distribution sin(Θ)P(Θ) at different temperatures, shown in Figure 2. As a result of weighting by sin (Θ), this distribution goes to zero at 180° (and also at 0°) and has a maximum for Θ < 180°. At T ) 0 K, the maximum of the weighted distribution from the DMC calculations is around 157°. Interestingly, this is very close to the equilibrium value of Θ, 156.8°, for the global minimum on the PES of two H2 molecules in the large cage. Thus, this

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20779

Figure 3. Isosurfaces of the 3D spatial distribution for two H2 molecules in the large cage generated using path integral molecular dynamics at the five chosen temperatures. The golden isosurface corresponds to 25% of the maximum isocontour value, and the white surface corresponds to 5% of the maximum value.

structural feature of the PES is clearly reflected in the quantum mechanical distribution despite strong quantum effects. However, even at T ) 0 K the angular distribution spans a wide range of angles, demonstrating that the system is highly fluxional already in the ground state. Increasing the temperature causes the peak to gradually shift away from 157° to smaller angles and broadens the distribution further. At 200 K, the distribution has the peak around 130°, and its angular range exceeds 100°. Considering the three distributions together, the conclusion is that the radial distributions P(R) and P(r) display much weaker temperature variations than P(Θ) [and sin (Θ)P(Θ)]. This implies that the PES is significantly stiffer along the two radial coordinates R and r than in the angular coordinate Θ. The H2 molecules are essentially confined to a spherical shell, which is ∼2 au wide and centered at R ∼ 3 au. There is only a very slight deviation from this at 200 K. As the temperature is increased, the most probable angle decreases while the distribution itself becomes increasingly broad, revealing a growing degree of angular delocalization of the molecules within the spherical shell. This is clearly confirmed by the 3D spatial distributions of Figure 3. B. Four H2 Molecules in the Large Cage. For four H2 molecules Figure 4 displays the distributions P(R) (top), P(r) (middle), and P(Θ) (bottom), from the PIMD calculations at the five chosen temperatures ranging from 25 to 200 K, and also from DMC calculations at T ) 0 K.16 Comparison of Figure 4 (top) with Figure 1 (top) shows that the distribution P(R) for four H2 molecules in the cage is displaced to substantially larger R values relative to P(R) for two H2 molecules at the same temperature. Thus, at T ) 0 K the maximum of P(R) is at 3.0 au for two H2 molecules, and it shifts to 3.7 au for four H2 molecules.16 Therefore, as more H2 molecules are added to the cage, somewhat paradoxically, the volume of empty space at

20780

J. Phys. Chem. C, Vol. 114, No. 48, 2010

Figure 4. Distribution functions of the H2-cage center distance (top), H2-H2 distance (middle), and H2-cc-H2 angle (bottom) for four H2 molecules in the large cage generated from path integral molecular dynamics at the five chosen temperatures, and from diffusion Monte Carlo (DMC) at T ) 0 K.16

the center of the cage, from which the guest molecules are excluded, increases. At higher temperatures, P(R) for four H2 molecules shifts as a whole to somewhat smaller values of R, indicating their ability to get closer to the center of the cage. The temperature dependence of the distribution P(r) in Figure 4 (middle) is qualitatively the same as that of the corresponding distribution for two H2 molecules shown in Figure 1 (middle); the peak of P(r) appears at decreasing r values as the temperature increases, and its position at a given temperature is nearly the same for four and two H2 molecules. As in the case of two H2 molecules, the most interesting distribution is the angular distribution P(Θ) shown in Figure 4(bottom). For all temperatures except the highest, 200 K, it peaks around 109.5°. This value of the H2-cc-H2 angle is an unambiguous sign that at lower temperatures the four H2 molecules are preferentially located at the vertices of a tetrahedron. The angular distributions at 25 and 50 K are indistinguishable from each other, within error bars (the difference in the peak heights is 0.0005, which is comparable to the estimated error in the difference), and from that at 0 K, in contrast to the distributions at 75 and 100 K which are distinct and broadened. Moreover, the differences in the distributions at 50, 75, and 100 K, although not large, lie well outside estimated error bars. This strongly suggests that between 50 and ∼100 K the system of four confined H2 molecules undergoes a qualitative structural change. We know from our earlier DMC study16 that at 0 K the (H2)4 tetrahedron has a specific, well-defined orientation relative to the water framework of the large cage, with each H2 molecule pointing toward the center of one of the four hexagonal faces of the large cage. This configuration also corresponds to the global minimum on

Witt et al.

Figure 5. Isosurfaces of the 3D spatial distribution for four H2 in the large cage generated using path integral molecular dynamics at the five chosen temperatures. The golden isosurface corresponds to 25% of the maximum isocontour value, and the white surface corresponds to 5% of the maximum value.

the PES of this system. Some 200 cm-1 above the global minimum there is a group of lowest-lying nearly degenerate local minima; they are also tetrahedral, but each is differently oriented inside the cage.16 It is therefore natural to infer that at 50 K and below the four H2 molecules are entirely localized in the neighborhood of the global minimum and maintain their unique orientation, whereas around 75-100 K the higher-energy local minima are accessible as well and are sampled by the H2 molecules. This also accounts for the change in the appearance of the angular distributions at these higher temperatures, apparent in Figure 4 (bottom). The confirmation for this interpretation of the temperatureinduced changes in the angular distribution comes from the 3D spatial distributions shown in Figure 5. The isosurfaces in the figure at 25 and 50 K depict a localized tetrahedral structure, whose orientation with respect to the cage is identical to that of (H2)4 tetrahedron at 0 K. On the other hand, from the 3D spatial distributions at 75 and 100 K it is evident that other basins on the PES, energetically above the global minimum, are becoming populated at elevated temperatures. At 200 K, the system is clearly able to explore considerably more of the PES, and is substantially delocalized. This is consistent with the much greater width of the 200 K angular distribution in Figure 4 (bottom), and the shift of its peak to a smaller angle. The overall conclusion from Figures 5 and Figure 4 (bottom) is that as the temperature is raised from 50 to ∼100 K, the system does undergoe a distinct change in the behavior. Below 50 K, it is localized within the global minimum basin, and the distributions are therefore very close the T ) 0 results of the DMC calculations. At temperatures higher than 50 K, we start to see the system probe the higher-energy local minima, whose population gradually increases with the raising temperature. Our results are in accord with the experimental observation by Lokshin et al.10 of a localization-delocalization transition, which four hydrogen (D2) molecules in the large cage exhibit

Temperature Dependence of Quantum Spatial Distributions at 50 K. Below 50 K their neutron diffraction data were best described in terms of localized guest molecules at the corners of a tetrahedron, directed toward the centers of the four hexagonal faces of the cage. Above 50 K, fitting the experimental data required using a linear combination of the (tetrahedrally) localized D2 model and a delocalized model of D2 rotating freely on the surface of a sphere, with the fraction of delocalized D2 molecules growing linearly with increasing temperature. IV. Conclusions We have reported path integral molecular dynamics (PIMD) calculations of the temperature dependence of the quantum spatial distributions of two and four H2 molecules in the large cage of sII clathrate hydrate. The cage was treated as rigid. Its geometry and the H2-cage and H2-H2 interaction potentials were taken from a recent DMC study of this system.16 The vibrationally averaged spatial configuration of the encapsulated H2 molecules was characterized by means of several 1D and 3D distribution functions generated from the PIMD calculations performed at 25, 50, 75, 100, and 200 K. They were combined with the T ) 0 K results of the DMC calculations.16 For both two and four H2 molecules, the two radial distributions, one for the distance between an H2 molecule and the cage center (cc) and another for the distance between pairs of H2 molecules, depend rather weakly on temperature and do not depart much from the corresponding T ) 0 K distributions from the DMC calculations. The exception are the distributions at the highest temperature considered, 200 K, for which the difference from the 0 K results is significant. On the basis of these distributions, the H2 molecules are completely excluded from a large portion of the central region of the cage and confined to a spherical shell whose width is about 2 au. The 1D distribution that we found to be the most sensitive to temperature variations, and also the most indicative of the structural changes in the system as a function of temperature, is that involving the H2-cc-H2 angle Θ. For two H2 molecules, when weighted by the Jacobian sin (Θ), the angular distribution at T ) 0 K peaks around Θ ) 157°, which is very close to the value of this angle at the minimumenergy configuration of two H2 molecules in the large cage. Thus, the most probable angle from the DMC calculations virtually coincides with its equilibrium value in the global minimum of the PES. With increasing temperatures, the maximum of the angular distribution shifts to smaller angles; at 200 K it is around 130°, reflecting a growing preference for strongly bent H2-cc-H2 configurations. At all temperatures, the angular distributions are broad, indicating large amplitude motions in the angular coordinate. In the case of four caged H2 molecules, the maximum of the angular distribution is at 109.5° at all temperatures except 200 K, where it appears at a slightly smaller angle. This implies a predominantly tetrahedral configuration of the guest molecules. The angular distributions at 0, 25, and 50 K overlap completely and are visibly different from those at 75 and 100 K (and 200 K). This suggests that between 50 and ∼100 K, a qualitative change takes place in the system. At 50 K and below, the four H2 molecules are localized within the global minimum and are preferentially located at the corners of a tetrahedron, which are oriented toward the centers of the four hexagonal faces of the large cage. At temperatures of 75-100 K and above, nearly degenerate local minima lying ∼200 cm-1 above the global minimum16 start being sampled as well, with their population

J. Phys. Chem. C, Vol. 114, No. 48, 2010 20781 increasing with the raising temperature. These local minima also correspond to (H2)4 tetrahedra, but each has a different orientation relative to the cage framework. The 3D spatial distributions also show very clearly that above 50 K the system is able to gradually explore increasing portions of the PES, and thus becomes more fluxional. Our results are consistent with the experimental finding of Lokshin et al.,10 that at 50 K the four hydrogen (D2) molecules in the large cage undergo a transition from being localized in a tetrahedral geometry with a unique orientation with respect to the cage to a state whose degree of delocalization within a spherical shell grows with increasing temperature. Acknowledgment. Z.B. is grateful to the National Science Foundation for partial support of this research. The computational resources used in this work were funded in part by the NSF MRI grant CHE-0420870. Acknowledgment is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research. M.E.T. acknowledges support from NSF CHE-0704036. A.W. acknowledges financial support from the German Academic Exchange Service (Deutscher Akademischer Austausch Dienst, DAAD). The PIMD simulations were carried out using computational recources of Bovilab@RUB and Rechnerverbund-NRW (Dortmund). References and Notes (1) Mao, W. L.; Koh, C. A.; Sloan, E. D. Phys. Today 2007, 60 (10), 42. (2) Struzhkin, V. V.; Militzer, B.; Mao, W. L.; Mao, H. K.; Hemley, R. J. Chem. ReV. 2007, 107, 4133. (3) Sloan, E. D. Clathrate Hydrates of Natural Gases; Marcel Dekker: New York, 1998. (4) Dyadin, Y. A.; Larionov, E. G.; Manakov, A. Y.; Zhurko, F. V.; Aladko, E. Y.; Mikina, T. V.; Komarov, V. Y. MendeleeV Commun. 1999, 9, 209. (5) Mao, W. L.; Mao, H. K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Science 2002, 297, 2247. (6) Mao, W. L.; Mao, H. K. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 708. (7) Schu¨th, F. Nature 2005, 434, 712. (8) Hu, Y. H.; Ruckenstein, E. Angew. Chem., Int. Ed. 2006, 45, 2011. (9) Strobel, T. A.; Hester, K. C.; Koh, C. A.; Sum, A. K.; Sloan, E. D., Jr. Chem. Phys. Lett. 2009, 478, 97. (10) Lokshin, K. A.; Zhao, Y.; He, D.; Mao, W. L.; Mao, H. K.; Hemley, R. J.; Lobanov, M. V.; Greenblatt, M. Phys. ReV. Lett. 2004, 93, 125503. (11) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Science 2004, 306, 469. (12) Lee, H.; Lee, J.-W.; Kim, D. Y.; Park, J.; Seo, Y.-T.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. J.; Ripmeester, J. A. Nature 2005, 434, 743. (13) Xu, M.; Elmatad, Y.; Sebastianelli, F.; Moskowitz, J. W.; Bacˇic´, Z. J. Phys. Chem. B 2006, 110, 24806. (14) Xu, M.; Sebastianelli, F.; Bacˇic´, Z. J. Phys. Chem. A 2007, 111, 12763. (15) Xu, M.; Sebastianelli, F.; Bacˇic´, Z. J. Chem. Phys. 2008, 128, 244715. (16) Sebastianelli, F.; Xu, M.; Bacˇic´, Z. J. Chem. Phys. 2008, 129, 244706. (17) Xu, M.; Sebastianelli, F.; Bacˇic´, Z. J. Phys. Chem. A 2009, 113, 7601. (18) Mak, T. C. W.; McMullan, R. K. J. Chem. Phys. 1965, 42, 2732. (19) Ulivi, L.; Celli, M.; Gianassi, A.; Ramirez-Cuesta, A. J.; Bull, D. J.; Zoppi, M. Phys. ReV. B 2007, 76, 161401. (20) Ulivi, L.; Celli, M.; Gianassi, A.; Ramirez-Cuesta, A. J.; Zoppi, M. J. Phys.: Condens. Matter 2008, 20, 104242. (21) Gianassi, A.; Celli, M.; Ulivi, L.; Zoppi, M. J. Chem. Phys. 2008, 129, 084705. (22) Strobel, T. A.; Sloan, E. D.; Koh, C. A. J. Chem. Phys. 2009, 130, 014506. (23) Strobel, T. A.; Taylor, C. J.; Hester, K. C.; Dec, S. F.; Koh, C. A.; Miller, K. T.; Sloan, E. D., Jr. J. Phys. Chem. B 2006, 110, 17121. (24) Chandler, D.; Wolynes, P. G. J. Chem. Phys. 1981, 74, 4078. (25) Hall, R. W.; Berne, B. J. J. Chem. Phys. 1984, 81, 3641.

20782

J. Phys. Chem. C, Vol. 114, No. 48, 2010

(26) Tuckerman, M. E.; Martyna, G. J.; Klein, M. L.; Berne, B. J. J. Chem. Phys. 1993, 99, 2796. (27) Ceperley, D. M.; Pollock, E. L. Phys. ReV. B 1984, 30, 2555. (28) Martyna, G. J.; Klein, M. L.; Tuckerman, M. E. J. Chem. Phys. 1992, 97, 2635. (29) Tuckerman. M. E. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: Oxford, 2010. (30) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6168. (31) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. J. Chem. Phys. 1996, 104, 5579.

Witt et al. (32) McDonald, S.; Ojama¨e, L.; Singer, S. J. J. Phys. Chem. A 1998, 102, 2824. (33) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2005, 123, 024507. (34) Diep, P.; Johnson, J. K. J. Chem. Phys. 2000, 112, 4465. (35) Tuckerman, M. E.; Yarne, D. A.; Samuelson, S. O.; Hughes, A. L.; Martyna, G. J. Comput. Phys. Commun. 2000, 128, 333.

JP107021T