New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Thermodynamic stability of hydrogen clathrates

Edited by Russell J. Hemley, Carnegie Institution of Washington, Washington, DC, and approved October 13, 2003 (received for review February 14, 2003)
Abstract
The stability of the recently characterized type II hydrogen clathrate [Mao, W. L., Mao, H.K., Goncharov, A. F., Struzhkin, V. V., Guo, Q., et al. (2002) Science 297, 2247–2249] with respect to hydrogen occupancy is examined with a statistical mechanical model in conjunction with firstprinciples quantum chemistry calculations. It is found that the stability of the clathrate is mainly caused by dispersive interactions between H_{2} molecules and the water forming the cage walls. Theoretical analysis shows that both individual hydrogen molecules and nH_{2} guest clusters undergo essentially free rotations inside the clathrate cages. Calculations at the experimental conditions – 2,000 bar (1 bar = 100 kPa) and 250 K confirm multiple occupancy of the clathrate cages with average occupations of 2.00 and 3.96 H_{2} molecules per D5^{12} (small) and H5^{12}6^{4} (large) cage, respectively. The H_{2}–H_{2}O interactions also are responsible for the experimentally observed softening of the H—H stretching modes. The clathrate is found to be thermodynamically stable at 25 bar and 150 K.
Clathrate hydrates are a class of inclusion compounds in which guests (noble gases or small organic molecules) occupy, fully or partially, cages in the host framework made up of Hbonded water molecules (1). Clathrate hydrate research is of fundamental and practical importance and involves a broad variety of scientific disciplines. The behavior of clathrate hydrates under pressure can provide valuable information on water–water interactions and interactions of water with a wide range of guests. Methane hydrate is the most abundant natural form of clathrate. An estimate of the global reserve of natural gas in the hydrate form buried in the permafrost and sediments underneath the continental shelf is significantly larger than that from traditional fossil fuels and will be a valuable future energy resource (2). On the negative side, the blockage of natural gas pipeline by solid hydrocarbon hydrates is a potentially hazardous and expensive problem that has not been fully resolved (3). Methane hydrates also represent a potential source of climate instability. As warming proceeds downward toward the seafloor and reaches the limits of hydrate stability, the hydrate will decompose, and some of the methane gas will escape to the atmosphere and increase the greenhouse effect. There are recent reports suggesting that a cause of ancient global warming and mass extinction of many forms of life 183 million years ago may be traced to sudden eruption of oceanic methane hydrate (4, 5). Clathrate hydrates also may be the most abundant form of volatile materials in the solar system. The possible existence of gas hydrate is crucial to the modeling of bodies in the solar system. The identification of several highpressure forms of methane hydrates has helped to reconcile the origin of the presence of large amount of methane in the atmosphere of Saturn's moon Titan (6). Solid clathrate hydrate also has been suggested to exist in comets in order to explain the large difference between the latent heats of vaporization of the various ices in the Whipple's model (7). Recently, the discovery of new species of centipedelike worms living on and within mounds of methane hydrate on the floor of Gulf of Mexico challenges the conventional view that deep sea bottom is a monotonous habitat and indicates that methane hydrate may play a role in marine ecosystem (8). More recently, the synthesis of H_{2} hydrate, a hydrate with the smallest guest with multiple occupancy, questioned the conventional theory for the prediction of the stability of clathrate hydrate (9) and raised the prospect of using clathrate hydrate as an efficient H_{2} storage medium.
Under normal conditions, clathrate hydrates are known to have three distinct crystalline structures (10). Both structure I and II clathrate hydrates have a cubic structure and a guest:host water ratio of ≈1:6 (10). Usually, very small guest (rare gas) favors the formation of the type II structure (10). For a very large guest molecule, the hexagonal structure H with larger cavities is often formed (11). It was a general perception that hydrogen molecules are too small to fit into the hydrate cages, and no stable clathrate hydrate structure can be formed (10). Not until recently, under highpressure conditions, hydrogen clathrate with the type II structure was synthesized and characterized (9). More surprisingly, the hydrogen clathrate was found to have an unusually high H_{2}/H_{2}O ratio (1:2) by highpressure xray diffraction, Raman, and infrared spectroscopy (9). The chemical composition of these clathrates can only be accounted for if multiple (up to quadruple) occupancies of the clathrate cages are assumed. Moreover, the new clathrate remains apparently stable even at the normal atmospheric pressure as long as the temperature is below 150 K. The experimental observations clearly indicate that hydrogen incorporation in these structures is associated with a low freeenergy process, much lower than would be expected for the mechanical encapsulation of the hydrogen gas. Additionally, the experimentally observed H—H stretching modes in the new clathrate are softened as compared to the free molecule (9). This observation is in contrast to the behavior of hydrogen gas under similar pressures (<1 GPa) (12). Multiple occupancy of the clathrate cavities is a rare phenomenon. The only existing example is nitrogen clathrate synthesized under pressurized nitrogen atmosphere (13). The purpose of this article is to investigate the stability of hydrogen clathrate, in particular, the effects of occupancy in the empty cavities, by using statistical mechanical theory with firstprinciples quantum chemistry calculations.
Theoretical Background
Hydrogen clathrate forms a type II structure (14) with 136 water molecules and 24 cages per unit cell. Sixteen of the cages are pentagonal dodecahedra (D5^{12}). The remaining eight cages are 16hedra (H5^{12}6^{4}; see refs. 1 and 10 for nomenclature and further information). To avoid possible confusion with the hydrogen guest, we will refer to the D5^{12} and H5^{12}6^{4} as S (for small) and L (for large) cages, respectively. The structures of the cages are shown in Fig. 1. Experimentally, it was found that the new clathrate contains 61 ± 7 (9) hydrogen molecules, distributed among the cages, per unit cell. Placing two guest molecules in each S cage and four in each L cage would give an approximate 1:2 H_{2}/H_{2}O clathrate composition.
The structure and stability of several multipleoccupancy ice clathrates have been examined with molecular dynamics methods (15, 16). Here, we adopt a statistical thermodynamics approach to investigate the thermodynamics of H_{2}/H_{2}O clathrate formation. Encapsulation of successive hydrogen molecules is treated as a series of chemical reactions between gaseous hydrogen with an (initially empty) ice cage C: [1] [2] [3] [4] [5] ..., etc.
At moderate pressures, reactions 1–5 are characterized by the corresponding equilibrium constants: where x_{n} is the fraction of the cages of a given type with nH_{2} molecules inside them (∑_{i}_{=0} x_{i} = 1), and p_{H2} is the hydrogen gas pressure. The equilibrium constants K_{p,n} can be estimated with standard statistical thermodynamics (17) by using the energetics obtained from firstprinciples calculations of the potential energy surfaces.
Computational Details
Density functional theory (DFT) (18–20) calculations were performed with the Amsterdam density functional (adf) program package (refs. 21–24 and www.scm.com) by using Cartesian space numerical integration (25) and analytical gradients for the geometry optimization (26). An uncontracted doubleζ basis of Slatertype orbitals was employed for the 2s, 2p shells of oxygen and 1s orbitals of hydrogen atoms. The oxygen 1s shell was treated with the frozencore approximation (22). Closedshell spinrestricted calculations were performed employing the Perdew–Burke–Ernzerhof (27) Generalized Gradient Approximation (GGA) functional with revised exchange parameterization of Zhang and Yang (revPBE) (28), applied selfconsistently (29). Unless otherwise stated, guest H_{2} molecules were fully optimized in a fixed model water cage. All closedshell ab initio secondorder MöllerPlesset (MP2) calculations were performed by using gamess program (30) at optimized DFT geometries. These calculations used split valence Gaussian 3–21G (31) atom basis sets, augmented by a single set of 2p polarization functions (ζ = 1.1) on hydrogen atoms. Core 1s electrons on oxygen were not correlated.
The small cage was modeled with a 20molecule water cluster (Fig. 1). In this structure, 8 oxygen atoms are arranged in a perfect cube (O_{h}) 3.83 Å from the cage center. The remaining 12 oxygen atoms are placed 3.96 Å from the center of the cage and transform according to the T_{h} subgroup of O_{h}. For the large cage, a 28molecule model cluster was constructed in a distorted T_{d} symmetry. Oxygen atoms in this cluster are located between 4.32 and 4.84 Å from the cage center. For both model cages, hydrogen atoms, with a fixed O—H bond distance of 1.0 Å, were placed so as to maximize the hydrogenbonded network. These structural models correspond closely to the average xray structure of type II clathrates (14).
Structure and Stability of the Model H_{2}/H_{2}O Complexes
For the small clathrate cage, we considered complexes encapsulated with one, two, and three hydrogen molecules. For the large cage, up to five hydrogen molecules were placed inside the cage. As expected for weakly bound clusters, the calculated energy minima for fully optimized H_{2} clusters are rather soft with many nearzero harmonic vibrational frequencies. For several structures (nH_{2}@L, n = 1, 2, 3, and 5), we obtained one or more imaginary vibrational modes corresponding to rotations of the hydrogen molecules. Attempts to remove the imaginary frequencies by following the corresponding normal mode did not lead to any significant decrease in the total energy. In most cases, the total energy is insensitive to the orientation of the individual hydrogen molecules. Furthermore, the intramolecular H—H distances of the clusters are very similar to the free molecule with the exception of 2H_{2}@S, where the H—H bond length is slightly elongated by 0.004 Å. The calculated reaction enthalpies at 0 K are collected in Table 1.
For H_{2}@S and H_{2}@L complexes, density functional calculations predict a shallow minimum near the center of the respective cages. This very soft minimum is likely an artifact of the approximate exchangecorrelation functional, which is unable to properly describe the dispersive interactions (see below). In both complexes, the hydrogen stretching vibration is softened by 87 cm^{–1} (H_{2}@S) and 21 cm^{–1} (H_{2}@L) with respect to that of the free molecule. The translational displacement of the hydrogen molecules (radial displacement) from the center of both cages is associated with very low harmonic vibrational frequencies (120 cm^{–1} in H_{2}@S and 70 cm^{–1} in H_{2}@L).
For 2H_{2}@C (C = S, L) clusters, the distance between center of mass of the H_{2} molecules in these clusters is found to be 2.58 Å (S) and 2.80 Å (L). DFT calculations predict both complexes to be less stable than the separate hydrogen molecules and the cage by 4.0 and 3.4 kcal/mol for the small and large cages, respectively. The hydrogen stretching frequency is again softened by 190 to 15 cm^{–1} in the small cage and by ≈10 cm^{–1} in the large cage. The very large spread in calculated stretching frequencies in the 2H_{2}@S cluster indicates the sensitivity of frequency shifts to small changes in the cluster geometry. Because we model each nH_{2}@C cluster with a single representative structure, only qualitative conclusions for the frequency shifts can be drawn. A direct comparison with experimental IR and Raman spectra would require a considerably more extensive molecular dynamics simulation.
When a threemolecule cluster is placed inside a cage, the positions of the hydrogen molecules relax to form an approximate equilateral triangle. In the S cage, the H_{2}—H_{2} separations are between 2.47 and 2.50 Å. In the L cage, the center of mass H_{2}—H_{2} distances increase to 2.82–2.86 Å. Hydrogen stretch frequency shifts range from –26 to –4 cm^{–1} in the S cage and from –19 to +14 cm^{–1} in the L cage. Both 3H_{2}@S and 3H_{2}@L are predicted to have positive formation enthalpies of +15 and +3 kcal/mol, respectively.
Because of its smaller size, the S cage cannot accommodate more than three hydrogen molecules under realistic conditions. The calculated DFT formation enthalpy for a model 4H_{2}@S structure is in excess of +30 kcal/mol. In contrast, in the L cage four hydrogen molecules arrange to form a slightly distorted tetrahedron with H_{2}—H_{2} distances between 2.90 and 3.13Å. Once again, the hydrogen stretch modes are predicted to be softened by –26 to –3 cm^{–1}. The calculated DFT formation energy of +5 kcal/mol is somewhat endothermic.
Finally, we considered a cluster of five hydrogen molecules inside the L cage. For this structure, geometry optimization converges to a distorted tetrahedron of H_{2} molecules (center of mass distances between 3.74 and 4.23 Å). The fifth hydrogen molecule is situated approximately at the cage center and is 2.42–2.53 Å from the remaining H_{2} molecules at the corners of the tetrahedron. The hydrogen stretch mode in this structure is broadened compared to the free H_{2} molecule with shifts between –33 and +37 cm^{–1}. The DFT formation energy of this structure is +14 kcal/mol.
Except for the lowoccupancy complexes (H_{2}@S, H_{2}@L, and 2H_{2}@L) the calculated reaction enthalpies (Table 1) for hydrogen molecule(s) encapsulation are too endothermic to account for the clathrate formation at temperature T = 250 K (RT = 0.5 kcal/mol, where R is the universal gas constant) and 2,000 bar (1 bar = 100 kPa). The inconsistency is not entirely unexpected because all the popular approximate density functionals are incapable of describing dispersive interactions (32, 33). Such interactions are expected to be the dominant attractive contribution for the neutral, nonpolar hydrogen molecule. To correct for this deficiency, we performed singlepoint ab initio MP2 calculations of the hydrogen clusters at optimized DFT geometries. The inclusion of dispersive interactions increases stability of the clathrate structures, particularly at higher hydrogen molecule occupancies (Table 1). For multiply occupied S and L cages, the difference between DFT and MP2 energies approaches 7 and 4 kcal/mol per hydrogen molecule, respectively. In comparison, the difference between DFT and MP2 reaction energies is much smaller for lowoccupancy complexes (1H_{2}@S, 1H_{2}@L, and 2H_{2}@L).
This apparently counterintuitive behavior can be explained in terms of the r^{–6} dependence of the dispersive interaction energy. By using the standard Lennard–Jones 6–12 parameters for van der Waals interactions (see, e.g., ref. 34), the H_{2}–H_{2}O interaction potential is expected to be attractive for R_{0} >2.8 Å with an energy minimum at ≈3.2 Å. Because the water cages' radii are significantly larger than R_{0} (S ≈ 3.8 Å and L ≈ 4.0 Å), displacement of H_{2} molecules from the cage center in multipleoccupancy structures is expected to enhance the dispersive interactions. In view of the very low calculated DFT harmonic frequencies (<130 cm^{–1}) for the radial translation of a single H_{2} molecule in both model cages, it is likely that the equilibrium position of the hydrogen will move offcenter once dispersive interactions are taken into account. Unfortunately, ab initio MP2 full geometry optimization of complexes of this size would have been prohibitively expensive. Because the cage models are approximately spherically symmetric, to a first approximation, the interaction potential may be examined through a series of singlepoint MP2 calculations with the H_{2} molecule displaced from the center of the cage. The results of these calculations are compared with the DFT interaction energies evaluated at the same geometries in Fig. 2.
For the S cage, the MP2 energy minimum is located at an offcenter displacement of R_{min} = 1.1 Å. The corresponding H_{2} incorporation energy of –2.9 kcal/mol is 1.1 kcal/mol lower than the –1.8 kcal/mol calculated at the cage center. For the L cage, the MP2 energy minimum is 2.0 Å offcenter with a binding energy of –2.5 kcal/mol. From the MP2 potential energy curves, one can deduce the harmonic vibrational frequency for hydrogen radial translation of 250 cm^{–1} in the small cage and 160 cm^{–1} in the large cage.
Notably, the shape of the repulsive part of the MP2 potential energy curves (R > R_{min}) closely resembles their DFT counterpart (Fig. 2). From geometrical considerations, the average distances between H_{2} molecules (H_{2}–H_{2}) and between H_{2} and water (H_{2}–H_{2}O) forming the cage walls are closer in multiply occupied cages. Therefore, apart from a constant energy shift, DFT potential energy surfaces can be expected to provide a reasonable description of clusters at higher hydrogen occupancies.
It is important to examine the origin of the stabilization of the hydrogen clusters inside the ice cages. The relative importance of the H_{2}–H_{2} and H_{2}–H_{2}O interactions can be estimated from the MP2 formation energies of the “naked” nH_{2} cluster evaluated at the optimized geometry inside the model cages (Table 1). It is found that a “magic” size exists (n = 2 for the S cage and n = 4 for the L cage) up to which the repulsive H_{2}–H_{2} interactions seem to be relatively unimportant and contribute <0.5 kcal/mol per H_{2} molecule to the reaction energy. At higher cage occupancies, H_{2}–H_{2} repulsion increases sharply. In contrast, the stability gained from attractive H_{2}–water cage wall interactions reaches saturation at the magic occupancy. This trend is the consequence of the total area of H_{2}–wall contacts, which reaches a maximum at the magic cluster size and remains almost constant for larger clusters (see Table 1). Therefore, the net result is that the hydrogen encapsulation energy shows a pronounced minimum at the magic occupancies.
If a similar analysis is applied to H_{2} harmonic stretch frequencies, it is clear that the repulsive intracluster interactions lead to an increase in the stretch frequencies, of 9 to 60 cm^{–1}, depending on the geometry. The only exception is in the 2H_{2}/2H_{2}@S cluster, where one of the vibrational modes is softened by ≈90 cm^{–1} as compared with 190 cm^{–1} in 2H_{2}@S. In this case, the contribution to the frequency shift arises mainly from the H—H bond elongation. Both energy and harmonic frequency decomposition suggest that at the experimental conditions of 2,000 bar and 250 K, H_{2}O–H_{2} interactions dominate the structure and properties of the clathrate, whereas the H_{2}–H_{2} interactions are of secondary importance. Consequently, caution is advised in applying the results of highpressure studies on hydrogen gas to the H_{2}/H_{2}O clathrate.
Estimation of Thermodynamic Parameters
To assess the stability of the clathrate at experimental conditions, finite temperature and zeropoint motion effects must be taken into account. By using the MP2 reaction energies at 0 K (Table 1) and DFT harmonic vibrational frequencies, standard statistical thermodynamics techniques (17) can be used to estimate the correction of zeropoint vibrational energy and thermal contributions to reaction enthalpies and reaction entropies. The results for the energetics (H_{T}E_{0}, S_{T}) at 150 K, are shown in Fig. 3. Similar behavior is found at 250 K, with thermal contributions to the entropy and enthalpy increasing by 3–5 cal/mol·K and 0.5–1.0 kcal/mol, respectively. The calculated compositions of the clathrate are collected in Table 2.
At the experimental conditions where the clathrate were synthesized (T = 250 K, P = 2,000 bar), the harmonic approximation leads to average cage occupancies of 1.8 (S cage) and 1.6 (L cage). These guest occupancies are insufficient to prevent the clathrate structure from collapsing (10). At T = 150 K and P = 25 bar, the average occupancy of the L cage drops to almost zero (Table 2). Clearly, the harmonic vibrational approximation does not lead to an adequate description of the thermodynamics of the H_{2}/H_{2}O clathrate cages. The reason for this failure becomes clear on examination of the individual contributions to the calculated free energies of the hydrogen incorporation reactions. For example, the calculated reaction energy for 4L (4H_{2}+L = 4H_{2}@L) is exothermic at 250 K (Δ_{r}E_{250} = –4.9 kcal/mol). This energy gain is associated with a large decrease in entropy (TΔS_{250} = 20.7 kcal/mol). The entropic contribution is, in turn, dominated by the loss of translational entropy of free H_{2} (27.2 cal/mol per K) with no contributions of similar magnitude appearing for the nH_{2} cluster inside the cage. Therefore, it appears likely that the harmonic vibrational approximation failed to describe highamplitude motions of the nH_{2} clusters inside the model cages and led to underestimation of the entropies. A more realistic treatment of the vibrational entropy clearly is required.
For a spherically symmetric rigid host cage and freely rotating (or monoatomic) guest, the classical partition function can be reduced to a onedimensional integral of the interaction potential. This leads to a simple expression for the first equilibrium constant (35, 36): where ν(r) is the guestcage interaction energy at displacement r from the center of the cage. In the temperature range where H_{2}/H_{2}O clathrate is stable (T < 250K), the assumption of classical motion breaks down for the light H_{2} molecule leading to overestimation of the equilibrium constants (S.P. and S. Yurchenko, unpublished data). Any attempt to predict the thermodynamics of the model clathrate cages with classical molecular dynamics techniques will suffer from the same problem (37).
To attain a computationally affordable semiquantitative description of the thermal contributions to the free energy of the H_{2} incorporation reaction, we propose the following scheme. We will refer to this procedure as the “improved” approximation.
For each nH_{2} cluster inside a rigid model cage the 6n degrees of freedom are split into three categories.

Two degrees of freedom for each individual hydrogen molecule are treated as a free rotor. At 250 K, this term contributes 0.5n kcal/mol to ΔH, and 2.7n cal/mol·K to ΔS. This assumption is justified by the observation that spacing between rotational energy levels in H_{2} (experimental value, B_{e} = 60.9 cm^{–1}; DFT, B_{e} = 60.1 cm^{–1}) is ≈0.2 kcal/mol, which is comparable to the orientational dependence of the H_{2}–cage interaction energy. Experimental Raman spectra also suggest essentially free rotation of the clathrated H_{2} molecules (9). Within this approximation, the H_{2} rotational contributions to reaction free energies cancel identically.

Three degrees of freedom per nH_{2} cluster are treated as a rigid body rotation with moments of inertia calculated at the equilibrium structure of the cluster relative to the center of the cage. For the 4H_{2}@L cluster, this term amounts to 0.8 kcal/mol for ΔH_{250} and 19.8 cal/mol·K for S_{250}. This assumption is justified by the approximately spherical symmetry of the cages and weakness of the individual H_{2}–H_{2}O interactions.

The remaining 4n3 degrees of freedom are treated as harmonic oscillators with frequencies taken from DFT harmonic vibrational analysis. The n highfrequency modes correspond to the H—H stretching mode. The remaining modes correspond to the breathing modes of the cluster. For 4H_{2}@L, this contribution adds 11.9 cal/mol per K to S_{250} and 30.4 kcal/mol to ΔH_{250}.
An additional justification is required for treating the rattling motion of the hydrogen molecule in a singly occupied cage as rigid body rotation around the cage center. For the temperature range of interest (150–250 K), the RT value amounts to 0.3–0.5 kcal/mol. This value should be compared with the barrier height at the cage center: ≈1.0 kcal/mol for the small cages and 2.0 kcal/mol for the large cages (see Fig. 1). At these temperatures the hydrogen molecule may be expected to be localized within the thin, spherical shell with R ≈ R_{min}.
Thermal contributions to formation energy and entropy calculated following this procedure (approximations 1–3) are shown on Fig. 3. For singly occupied cages, the improved values compare favorably with the results from the direct solution of Schroedinger equations for the nuclear motion of the guest (S.P. and S. Yurchenko, unpublished data). In general, the approach outlined above leads to significantly lower calculated H_{T}E_{0} contributions and higher entropies than the harmonic approximation. The bulk of the change in the enthalpic contribution is attributable to the absence of the zeropoint energy from the free rotation contributions. The increase in the calculated entropy compared to the harmonic oscillator approximation derives largely from the rotations of the cluster. Treating these degrees of freedom as free rotations is essential for achieving agreement between theory and experiment.
Clathrate compositions at several (p, T) values calculated with the improved approximation are collected in Table 2. The average occupancies of the S and L cages at two representative temperatures are depicted in Fig. 2. At 250 K, the S cage is already fully occupied at relatively low pressures (≈100 bar). However, the bigger L cage is still empty at this pressure. At pressures >1,000 bar, both cages become fully occupied. This result agrees very well with the experimental formation conditions (T = 250 K, P = 2,000 bar). In view of the relatively crude approximations made in the calculations, this agreement should be considered excellent.
At 150 K, calculations suggest that the S cage should be able to retain both hydrogen molecules at H_{2} pressures well below 1 bar, whereas the L cage is already fully occupied at P > 25 bar. Again, the calculated results agree well with experiment and suggest that type II H_{2}/H_{2}O clathrate is thermodynamically stable with respect to H_{2} loss under mild temperature and pressure conditions. In the experiment, by warming the clathrate from 78 K to 115 K (9), it was observed that the intensity of the infrared absorption band assigned to vibrations of hydrogen in the small cages decreases rapidly. This observation may infer that the hydrogen content in the small cage decreases faster than in the large cages with increasing temperature, which is in apparent contradiction with the theoretical prediction. However, it should be recognized that the vibrational frequency of hydrogen molecule inside the cavities is very sensitive to the local environment. It may not be appropriate to associate that higher frequency vibrations with the motions of hydrogen in the more spatially confined smaller S cavity. In fact, the converse is observed for the C—H stretching modes of methane in methane clathrate (38). Furthermore, because the H—H stretching frequency is very sensitive to the location of the molecule inside the cavity, we may expect a very pronounced thermal broadening effect, leading to an apparent intensity decrease.
The analysis described above suggested that the hydrogen molecule clusters are rigid rotors in a weakly perturbed environment. We can then estimate the quantum rotation (tunneling) splittings in the limit of a free rigid rotor from the respective rotational constants. For hydrogen dimer in the S cage (2H_{2}@5^{12}), the maximum energy splitting is ≈2 cm^{–1}, and for the hydrogen tetramer in the large cage (4H_{2}@5^{12}6^{4}), the splitting is ≈0.9 cm^{–1}. These splittings may be observable in microwave spectroscopy experiments.
Conclusions and Outlook
We have examined the stability of the type II H_{2}/H_{2}O clathrate with respect to the hydrogen occupancy by using DFT and ab initio MP2 calculations on model water cages. The experimentally observed composition and formation conditions can only be explained if two assumptions are made: (i) individual hydrogen guest molecules undergo free rotation and (ii) the collective motions of the cluster of hydrogen molecules inside the individual host cages must be treated as rigid body rotations. Under these assumptions, the clathrate reaches the observed ≈1:2 H_{2}/H_{2}O composition at T = 250 K, P > 1,000 bar.
The incorporation of H_{2} into the cages is calculated to be exothermic (Δ_{r}H_{250} = –1.8 kcal/mol), and the stability of the clathrate increases at lower temperatures. Below 150 K, the clathrate should be thermodynamically stable at nearambient pressure, making it a promising and easy way to handle stored hydrogen. The stabilization of the hydrogen guest clusters arises mainly from attractive dispersive interactions with the water on the cage walls. As the surface contact is maximized for similar sizes of the guest and the host cage, sites within the smaller D5^{12} cage are more stable (Δ_{r}H_{250} = –2.6 kcal/mol) than in the H5^{12}6^{4} cage (Δ_{r}H_{250} = –1.2 kcal/mol). Thus, double clathrates with a suitable choice of guest in the larger cage (e.g., tetrahydrofuran) should form under even milder conditions. If hydrogen loss from the D5^{12} cage remains the stabilitydetermining factor, such clathrates should be stable at hydrogen pressures of ≈20 bar (T = 250 K).
Interactions with the cage walls are responsible for the experimentally observed softening of the H—H stretching modes in the clathrate. If these interactions are ignored, the stretching frequencies actually increase, as expected for moderately compressed hydrogen gas (9).
It should be noted that a more definitive determination of the occupancy number can be obtained from structural refinement of experimental neutron diffraction patterns. Furthermore, microwave spectroscopy will be invaluable to characterize the tunneling splittings if the hydrogen molecule clusters were indeed behave as rigid rotors.
In summary, we have performed a theoretical investigation of the stability of hydrogen clathrate with respect to the guest occupancy. Although the occupancy may well be the dominating stability factor, other thermodynamic factors such as the stability of the ice framework at different cage occupations and kinetic (e.g., H_{2} diffusion rates) factors need to be examined. Moreover, we only considered a single representative structure for each occupancy of the host cage. The structures of the guest clusters appear to be highly fluxional, leading to broad complex features in the experimental IR and Raman spectra. Understanding these features would require a detailed quantum molecular dynamics simulation of the guest clusters (37). In this study, quantum effects have been included in the consideration of the rotational motions only. The quantum contributions to the translational motions were neglected. It is expected that the equilibrium constants will be somewhat reduced if this effect was included. Finally, it is clear that the discovery of the type II hydrogen clathrate is opening new exciting horizons in the field of clathrate research.
Acknowledgments
We thank Wendy Mao, David Mao, and Russell Hemley for communicating their results prior to publication and for many helpful discussions.
Footnotes

↵* To whom correspondence should be addressed. Email: john.tse{at}nrc.ca.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: DFT, density functional theory; MP2, secondorder MöllerPlesset.
 Received February 14, 2003.
 Copyright © 2003, The National Academy of Sciences
References
 ↵
Jeffrey, G. A. (1984) in Inclusion Compounds, eds. Atwood, J. L., Davies, J. E. D. & MacNicol, D. D. (Academic, New York).
 ↵
Suess, E., Bohrmann, G., Greinert, J. & Lausch, E. (1999) Sci. Am. 281 (5), 52–59.
 ↵
Sloan, E. D. (2000) Hydrate Engineering Monograph (Soc. Plastic Engineer., Richardson, TX), Vol. 21.
 ↵
Katz, M. E., Pak, D. K., Dickens, G. R. & Miller, K. G. (1999) Science 286, 1531–1533.pmid:10567252
 ↵
 ↵
 ↵
Delsemme, A. H. & Seings, P. (1952) Ann. Astrophys. 15, 1–6.
 ↵
 ↵
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. (2002) Science 297, 2247–2249.pmid:12351785
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Irikura, K. K. (1998) in Computational Thermochemistry: Prediction and Estimation of Molecular Thermodynamics, eds. Irikura, K. K. & Frurip, D. J. (Am. Chem. Soc., Washington, DC).
 ↵
Parr, R. G. & Yang, W. (1989) Density Functional Theory of Atoms and Molecules (Oxford Univ. Press, Oxford).
 ↵
Salahub, D. R., Castro, W. & Proynov, E. Y., (1994) in Relativistic and Electron Correlation Effects in Molecules and Solids, ed. Malli, G. L. (Plenum, New York).
 ↵
Baerends, E. J., Autschbach, J. A., Bérces, A., Bo, C., Boerrigter, P. M., Cavallo, L., Chong, D. P., Deng, L., Dickson, R. M., Ellis, D. E., et al. (2002) adf2002.1 (Scientific Computing and Modelling NV, Theoretical Chemistry, Vrije Universiteit, Amsterdam).
 ↵

Fonseca Guerra, C., Visser, O., Snijders, J. G., te Velde, G. & Baerends, E. J. (1995) Methods and Techniques in Computational Chemistry METECC95, eds. Clementi, E. & Corongiu, G. (STEF, Cagliari, Italy), p. 305.
 ↵
te Velde, G., Bickelhaupt, F. M., Baerends, E. J., Fonseca Guerra, C., van Gisbergen, S. J. A., Snijders, J. G. & Ziegler, T. (2000) J. Comp. Chem. 22, 931–967.
 ↵
te Velde G. & Baerends, E. J. (1992) J. Comput. Phys. 9, 84–98.
 ↵
 ↵
 ↵
 ↵
Patchkovskii, S. (2002) SelfConsistent Implementation of Perdew–Burke–Ernzerhof and related GGAs inadfTechnical Report (Univ. of Calgary, Calgary, AB, Canada).
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵