Structures of the I-, II-and H-Methane Clathrates and the Ice− Methane

ARTICLE pubs.acs.org/JPCA

Structures of the I-, II- and H-Methane Clathrates and the Ice-Methane Clathrate Phase Transition from Quantum-Chemical Modeling with Force-Field Thermal Corrections Annika Lenz and Lars Ojam€ae† Department of Chemistry, IFM, Link€oping University, SE-58183 Link€oping, Sweden

bS Supporting Information ABSTRACT: Methane hydrates with the three clathrate structures I, II, and H are studied by quantum-chemical methods. Hybrid density-functional theory B3LYP computations using periodic boundary conditions are combined with force-field methods for the thermal energy effects to calculate energetic, thermodynamic, and structural properties. The pressure dependencies for the crystal structures, lattice energies, and guest molecule interactions are derived. The quantum-chemical geometry optimizations predict too small cell volumes as compared to experimental data, but by including zero-point energy and thermal energy effects, we find the cell volumes increase and the correct densities are obtained. The phase transition from MH-I to ice Ih and methane was computed and found to occur at about 9.7 MPa.

1. INTRODUCTION Methane hydrates are solid-state phases formed by a mixture of methane and water molecules. The water molecules form clathrate structures consisting of cages enclosing the guest molecules, i.e., the methane molecules in methane hydrate. Methane hydrates form at low temperatures (close to 0 °C) and high pressures (a few MPa and up). In nature, methane hydrates can be found in permafrost regions and on the deep ocean floor. In pipelines, the conditions favor the formation of hydrates where they can block the flow and result in damage.1-3 Methane hydrates can play important roles in both environmental and energy supply issues. Since the methane hydrates have been estimated to contain around 12% of the organic carbon on earth,4 they constitute a large unemployed energy resource. On the other extreme, methane hydrates may have a negative effect, in that methane is a very efficient greenhouse gas. Leakage of methane from the hydrates, caused either by human attempts to recover the methane or by potential global warming, could possibly increase the greenhouse effect and cause further melting of the methane hydrates. Fast temperature changes in the history have been attributed to similar scenarios.2 A third area of interest concerns the ability of clathrates to enclose gas molecules of volumes up to 160 times the volume of the clathrate itself, which makes them of interest for storage and transport of other gases like natural gas,4 carbon dioxide,5 and hydrogen.6,7 There are three main crystalline structures of clathrates, termed I, II, and H, which contain cage formations of different dimensions. The number of cages, their sizes, and the number of water molecules per unit cell are listed in Table 1. The cage sizes r 2011 American Chemical Society

are described by the number of faces of each size. For example 51262 means that a cage is formed by 12 five- and 2 six-membered rings. For the methane hydrate, structure I is the most common structure but also other structures occur.8 Structures II and H most common for other guest molecules can also be formed for methane hydrate at higher pressure conditions.8-13 Theoretical insight into the factors that determine the properties of clathrates, and in particular how the clathrate stability depends on the presence of guest molecules and on pressure and temperature, is useful when discussing environmental issues and considering potential applications for clathrate structures. Most theoretical studies of methane hydrates have been performed using empirical potential models (for example, refs 14-18). Quantum-chemical calculations have been performed for guest molecules in a free water cage by Hermida-Ramon et al.,19 and Ikeda and Terakura studied the transformation from MH-I to filled ice using quantum-chemical CPMD simulations.20 In the present work the pressure dependencies of the crystal structure, lattice energies, and guest molecule interactions are studied using quantum-chemical periodic hybrid density functional theory B3LYP calculations in combination with force-field methods for zero-point energy and thermal corrections. Thermodynamic quantities are computed and the MH-I to ice Ih transition was studied. Special Issue: Victoria Buch Memorial Received: November 29, 2010 Revised: January 28, 2011 Published: February 22, 2011 6169

dx.doi.org/10.1021/jp111328v | J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

ARTICLE

Table 1. Number of Water Molecules and Cages of the Various Sizes in One Unit Cell for the Clathrate Structures I, II, and H cage typea and no. of H2O in the cage wall no. of H2O

512 and 20

I/cubic

46

2

II/face centered cubic

136

16

H/hexagonal

34

3

structure/unit cell

435663 and 20

51262 and 24

51264 and 28

51268 and 36

6 8 2

1

5 6 means that the cage is formed by x five- and y six-membered rings.

a x y

2. METHOD The crystallographic oxygen framework for the three structures (I, II, and H) were taken from refs 21-23. Regarding the positions of the H atoms, a proton-ordered H-bond network was created for each of the clathrate unit cells. For each H-bond in these structures there are basically two different relative orientations for the two water molecules participating in the bond (“C” and “D” pairs using the terminology in ref 24). The H atoms were inserted manually in such a way that approximately the most probable ratio between the two types of H-bonded pairs was obtained. The energetics of H-bond topologies has earlier been discussed for ice24-26 and water clusters.27-30 To form the methane hydrates, the clathrates were filled with one methane molecule in each cage. Their structures are shown in Figure 1. The lowestenergy proton-ordered form of ice Ih24,31 is also included in the study. Quantum-chemical electronic structure calculations were performed using the periodic HF and DFT quantum-chemistry program Crystal0332 employing the B3LYP33,34 hybrid density functional and the 6-31G(d,p)35 basis set. Geometry optimizations of the crystallographic structures were performed for a series of cell volumes to find each energy-minimum cell structure and to evaluate the dependence of the cell energy and structure on the pressure. The cell axis length ratios and angles were held constant when the cell volume was varied. Since the Crystal03 program does not compute analytical second derivatives, and the computational cost to obtain these numerically would be prohibitive, zero point energy (ZPE) and thermal energy corrections36 were instead taken from force-field computations. The vibrational thermal energy corrections are taken from a normal-mode analysis using the module Discovery in the Materials Studio program37 with the SPC38 and CVFF39-41 analytical model potentials for the water and methane molecules, respectively. A flexible model was used for water where the intramolecular part was taken from the CVFF force field.41 These potentials were found to give a very similar energy-volume dependence as that obtained from the B3LYP/6-31G(d,p) calculations. The vibrational thermal energy and entropy corrections are computed at the temperature 273 K for the series of cell volumes used in the B3LYP calculations. The energies are presented as the energy difference between the energy of the clathrate unit cell and the sum of energies of the corresponding number of free water and methane molecules. The energies are presented relative to that of one water molecule: E = [E(unit cell) - E(free monomers)]/n(H2O). The computed interaction energies using B3LYP/6-31G(d,p) are considerably larger in magnitude than the experimental energies, due to limitations in both the functional and the basis sets. To partly compensate for this, the electronic energy minima are scaled in such a way that the computed electronic energy for

Figure 1. I-, II-, and H-clathrate structures with one methane molecule in each cage. For MH-II the primitive cell with 34 water molecules is shown.

ice Ih at zero pressure and zero temperature agrees with the experimentally estimated energy (sans ZPE) for ice Ih at zero 6170

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

ARTICLE

pressure and zero temperature, -58.82 kJ/mol.42 EðscaledÞ ¼ Emin ðcalcÞ 3 F þ ðEðcalcÞ - Emin ðcalcÞÞ

ð1Þ

The scaling factor F in eq 1 was 0.6871 for the B3LYP/6-31G(d,p) energies and 0.9650 for the SPC/CVFF energies respectively. The pressure can be calculated from the volume dependence of the Helmholtz free energy, A, which here is the sum of the electronic interaction energy E, the vibrational zero-point energy, and thermal corrections Evib and the entropy contribution -TS:36 A ¼ U - TS ¼ E þ Evib - TS

 DA DðU - TSÞ p¼¼DV T DV T

 DU DS ¼þT DV T DV T


 DE DEvib DS ¼þT DV T DV T DV T

ð2Þ

ð3Þ

From the pressure volume dependence the isothermal compressibility, κ, can be obtained:36
 1 DV k¼ð4Þ V DP T The vibrational thermal energy correction including the zeropoint energy and the entropy are given by36
 hνi 1 1 þ Evib ¼ nR ð5Þ 2 expðhνi =kTÞ - 1 i k

∑

S ¼ nR

∑i




 hνi =kT - lnð1 - expð - hνi =kTÞÞ expðhνi =kTÞ - 1 ð6Þ

where νi are the vibrational frequencies, k is Boltzmann’s constant, and h is Planck’s constant. The Gibbs free energy, G, is be obtained from the Helmholtz free energy and the pressure:36 G ¼ A þ pV

ð7Þ

For the reference monomers, translational and rotational thermal energy corrections must be included. Interaction energies are calculated relative to noninteracting monomers in the gas phase at 1 atm and 273.15 K. The translational energy correction at temperature T is obtained from36 !
 2πMkT 3=2 Ve ð8Þ Atrans ¼ - nRT ln h2 nNA where M is the mass of a the molecule and V is the volume for an ideal gas at 1 atm. The rotational energy correction is given by36 0 !1=2 1 1=2 3 π T A ð9Þ Arot ¼ - nRT ln@ σ ΘA ΘB ΘC where σ is the rotational symmetry number (σ = 2 for H2O and

Figure 2. Electronic energies, ΔE, obtained from the quantum-chemical calculations for MH-I (O), MH-II (/), MH-H (0) (black lines), the empty clathrate structures (red lines), and ice Ih (Δ, blue line) versus the volume per water molecule.

σ = 12 for CH4) and ΘA,B,C = p2/2kIA,B,C are the rotational temperatures where IA,B,C are the moments of inertia. The methane molecules in the cavities of the clathrates are only weakly bound. A normal-mode analysis consequently revealed that many near-zero frequency vibrations were present, which were consistent with the three translational and three rotational modes that each methane molecule would exhibit if it was unbound. We therefore modeled the methane molecules to possess free rotational and translational motions in the cavities, and the corresponding low-frequency methane normal modes were thus removed. The translational and rotational free energies of methane in a cavity are instead given by eqs 8 and 9, where the volume was taken as the cavity volume (estimated from the smallest wall water oxygen-cavity center distance for each cavity type subtracted by the oxygen van der Waals radius). This is probably a fair approximation for the structures at the relatively moderate pressures present in this study, whereas it might be more adequate to model the CH4 intermolecular motions as vibrations at elevated pressures. To model the equilibrium between methane hydrates and ordinary ice Ih, the methane molecules are for simplicity assumed to form a separate phase in coexistence with ice. For methane the ideal gas law at high pressures will fail since it predicts volumes smaller than the volume of the methane molecule. Instead, we make the approximation that a hard sphere equation of state43 is valid, where the total volume includes the volume of the methane molecule: V ¼ nVCH4 þ

nRT p

ð10Þ

The diameter of a methane molecule was crudely estimated from the energy minimum distance between two methane molecules (using B3LYP/6-31G(d,p)). The volume of one molecule is then obtained from the volume of a unit cell built from close-packed spheres of that diameter divided by the number of spheres (VCH4 = 4.26  10-5 m3/ mol).

3. RESULTS AND DISCUSSION Structures and Energetics for Geometry-Optimized Geometries. The electronic energies for the three methane hydrates 6171

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

ARTICLE

Table 2. Interaction Energies (Electronic, ΔE, and Helmholtz Free Energy, ΔA, at 273.15 K) for the Geometry Optimized (Coordinates and Cell Volumes with Cell Axis Ratios Fixed) Electronic Energy-Minimum Structures for the I, II, and H Methane Hydrates and Empty Clathrates and for Ice Iha ΔE type

empty clathrate

ΔA methane hydrate

empty clathrate

methane hydrate

structure I

-58.73 (-59.21)

-60.19 (-63.04)

-2.89 (-3.37)

-0.83 (-3.68)

structure II

-58.80 (-59.44)

-59.87 (-63.00)

-2.73 (-3.37)

-0.33 (-2.88)

structure H

-58.53 (-58.89)

-59.63 (-62.47)

-2.44 (-2.80)

ice Ih

-59.30 (-59.32)

0.35 (-2.49) -3.44 (-3.58)

a

The energies are given per water molecule and relative noninteracting monomers in the gas phase at 1 atm and 273.15 K (kJ/mol). SPC/CVFF energies are given in brackets.

Figure 4. Radii for the cages in MH-I (black), MH-II (blue), and MH-H (green) versus the volume per water molecule. Figure 3. Electronic interaction energy for one free methane molecule absorbed in an empty clathrate structure for MH-I (O), MH-II (/), and MH-H (0) versus the volume per water molecule.

(Figure 1), the empty clathrates, and ice Ih as a function of volume per water molecule are shown in Figure 2. As expected, ice Ih is seen to be a denser structure than the clathrates at energy minimum. The methane-filled clathrates are seen to have larger energy-minimum volumes than the empty clathrates. When the clathrate structures are filled with methane, their electronic interaction energies decrease due to the attractive watermethane dipole-induced dipole forces.44 For the geometry-optimized structures electronic energies and Helmholtz free energies (at 273 K) relative the free monomers are shown in Table 2. For the empty clathrates structure II is lowest in electronic energy while structure I is lowest in Helmholtz free energy. For the methane hydrates MH-I is lowest in electronic and Helmholtz free energy. Structure H is higher in energy than structure I and II for both empty clathrates and methane hydrates when the electronic energy as well as when the Gibbs free energy is considered. Ice Ih is as expected the lowest in electronic energy of the pure water structures. However, that the structures are geometry-optimized does not imply that the pressures on the crystals are zero, as we will see below. The geometry-optimized (coordinates and cell volumes with the same fixed axis ratios as in the B3LYP calculations) energies obtained using the SPC/CVFF potentials are also listed in Table 2. The analytical-potential calculations manage to predict

Figure 5. Volume-pressure diagram for empty clathrates, methane hydrates, and ice Ih, where the pressure is obtained from the negative free energy (quantum-chemical electronic energies plus thermal and entropic corrections from analytical-potential calculations) derivative with respect to volume.

the same relative ordering of the geometry-optimized methane hydrate and clathrate energies as the quantum-chemical calculations does. Contrary to the quantum-chemical results the SPC 6172

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

ARTICLE

Table 3. Interaction Energies at Zero Pressure for the I, II, and H Structures of Methane Hydrate and Empty Clathrates and for Ice Iha S-I

S-II

S-H

MH-I

MH-II

MH-H

ice Ih

ΔE

-57.21 (-57.50)

-57.56 (-57.77)

-57.34 (-57.86)

-58.93 (-61.77)

-58.51 (-60.71)

-58.13 (-61.49)

-59.51 (-58.23)

ΔH ΔG

-44.18 (-44.62) -3.44 (-4.09)

-44.46 (-44.93) -3.40 (-4.19)

-44.30 (-44.79) -3.30 (-3.84)

-45.01 (-47.81) -1.70 (-4.42)

-44.71 (-47.14) -1.15 (-3.93)

-44.45 (-47.53) -0.77 (-3.37)

-45.00 (-45.04) -3.87 (-4.25)

V

22.65

22.64

23.24

22.40

23.09

23.76

19.31

Electronic energy, ΔE, enthalpy, ΔH and Gibbs free energy, ΔG (which equals ΔA at 0 Pa), are given at 273.15 K and 0 Pa. The molar volumes are given per mol water molecules (cm3). The energies are given per water molecule and relative noninteracting monomers in the gas phase at 1 atm and 273.15 K (kJ/mol). SPC/CVFF energies are given in brackets. a

Table 4. Cell Sizes and Densities for the B3LYP/6-31G(d,p) Geometry-Optimized (Electronic Energy-Minimum) Structures and for the Predicted Structures at the Pressure 0 Pa and at the Pressures Employed in Selected Experiments Where the Cell Sizes/ Densities Are Knowna type

cell size (Å)

density (g/cm3)

S-I

geometry-optimized

11.59  11.59  11.59

0.884 (0.88)

90°, 90°, 90°

0 Pa, 273 K

11.9  11.9  11.9

0.795 (0.80)

0.1 MPa, 273 Kc geometry-optimized

12.0  12.0  12.0b 11.81  11.81  11.81

0.81b 0.873 (0.87)

primitive cell

0 Pa, 273 K

12.1  12.1  12.1

0.796 (0.80)

60°, 60°, 60°

experimental 12.2  12.2  12.2b

0.80b

experimental S-II

0.1 MPa, 273 Kc S-H

geometry-optimized

11.84  11.84  9.83

0.852 (0.85)

90°, 90°, 120°

0 Pa, 273 K

12.2  12.2  10.1

0.775 (0.78)

MH-I

geometry-optimized

11.62  11.62  11.62

1.013 (1.02)

90°, 90°, 90°

0 Pa, 273 K 0.5 GPa, 273 K

11.9  11.9  11.9 11.8  11.8  11.8

0.929 (0.93) 0.976

12.0  12.0  12.0

0.92b

experimental quenched from 0.5 GPa, 273 Kd MH-II

geometry-optimized

11.88  11.88  11.88

0.993 (1.01)

primitive cell

0 Pa, 273 K

12.2  12.2  12.2

0.903 (0.92)

0.25 GPa, 273 K

12.2  12.2  12.2

0.930 0.92b 0.969 (0.98)

60°, 60°, 60°

experimental MH-H

0.25 GPa, 298 K e geometry-optimized

12.1  12.1  12.1b 11.91  11.91  9.89

90°, 90°, 120°

0 Pa, 273 K

12.2  12.2  10.1

0.896 (0.90)

0.6-0.8 GPa, 273 K

12.1  12.1  10.1

0.931-0.945

0.6 GPa, 298 Ke

12.0  12.0  10.0

0.95b

0.8 GPa, 298 Kf

12.0  12.0  10.0

0.95b

ice Ih

geometry-optimized

4.38  7.59  7.16

1.004 (0.99)

90°, 90°, 90°

0 Pa, 273 K experimental

4.50  7.80  7.35

0.933 (0.91)

0 Pa, 269 Kg

4.52  7.83  7.38b

0.917

experimental

a

Densities using SPC/CVFF energies are given in brackets. b Recalculated from the cell dimensions in the experimental references. c Reference 50. d Reference 51. e Reference 8. f Reference 10. g Reference 52.

potential gives the empty clathrate structure II to be slightly more stable than ice Ih (at least for the current set of fixed axis ratios). The electronic interaction energies, ΔE, for a methane molecule interacting with the surrounding clathrate structures are shown in Figure 3 as a function of volume. The variation of the radii (as obtained from the closest O-O distance of water molecules on opposite sides of the wall) of the methane clathrate cavities as a function of the clathrate molar volume is seen in

Figure 4. From Figure 3 it is seen that it is energetically more favorable to fill structure I with methane than to fill structure II or H. The highest energy gain (i.e., the largest absorption energy) is seen for larger volumes, where favorable attraction dominates over any remaining unfavorable steric repulsion forces. For structure I the interaction energies at volumes above 21 cm3/ mol of H2O are almost constant. The magnitude of these methane interaction energies could possibly be somewhat 6173

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

Figure 6. Gibbs free energies, ΔG, at 273.15 K for ice Ih plus the Gibbs free energy of a nonideal methane gas and for methane hydrate MH-I as a function of pressure, showing the phase transition corresponding to the decomposition of methane hydrate into ice Ih and methane gas. The inset shows the corresponding diagram for MH-II and MH-H.

underestimated due to the failure of the DFT functional to account for dispersion interactions (e.g., in the methane dimer the interaction energy obtained by B3LYP is ∼0.1 kJ/mol CH4, whereas that from an ab initio second-order perturbation method45,46 is ∼0.5 kJ/mol). The here adopted energy scaling procedure, however, probably at least partly compensates for this shortcoming. Volume-Pressure Relationship and Properties at Zero Pressure. A calculated pressure (eq 3) versus volume diagram is shown in Figure 5. An interesting note is that paradoxically the computed pressures (obtained from the free energy derivative) at the electronic energy-minima volumes (from Figure 2) are not zero but found to be about 500-1300 MPa. It is apparent that the cell volume corresponding to zero pressure increases when vibrational zero-point and thermal energies are added. The volume is seen to be compressed more rapidly for the empty clathrate structures than for the methane hydrates when the pressure is increased. The calculated isothermal compressibility at zero pressure is 0.118 GPa-1 for MH-I, 0.128 GPa-1 for MHII, 0.164 GPa-1 for MH-H, and 0.1040 GPa-1 for ice Ih. In ref 47 the isothermal compressibilities were estimated from experimental data to be 0.128 GPa-1 for MH-I, 0.133 GPa-1 for a methane/ethane hydrate structure II, and 0.116 GPa-1 for ice Ih (extrapolated to zero pressure). The electronic energies, enthalpies, and free energies at zero pressure for the different clathrate structures are given in Table 3. The lowest electronic energy structure at zero pressure is structure II for the empty clathrates and structure I for the methane hydrates, which was also found for the geometryoptimized structures in Table 2. For the methane hydrates the energetic stability decreases in the sequence I, II, and H, whereas for the empty clathrates the sequence is II, H, and I, where thus structure I is the least stable at zero pressure. These relative orderings are valid also for the enthalpies. But when Gibbs free energies are considered, structure I becomes the most stable structure, followed by II and lastly H, both for the empty clathrates and for the methane hydrates.

ARTICLE

The difference in enthalpy between that the empty clathrate and ice Ih is found to be 820 J/mol for S-I and 540 J/mol for S-II. Earlier studies48-50 have found enthalpy differences in the range from 714 to 1398 J/mol for S-I and from 538 to 1025 J/mol for S-II. The differences in Gibbs free energy for S-I and S-II relative ice Ih are in the present study computed to be 430 and 480 J/mol, respectively. Data in the literature vary from 537 to 1297 J/mol for S-I and from 414 to 955 J/mol for S-II.48-50 We note that the clathrate I and II - ice Ih enthalpy and free energy differences obtained in the present study are within about 100 J from the values found in the molecular dynamics simulation study by Jacobson et al.48 In Table 4 the cell geometries and the corresponding densities are listed for the electronic energy-minimum structures and for the zero-pressure (free energy-minimum) structures. The cell lengths for the quantum-chemical electronic energy-minimum structures are too small compared to the experimental cell parameters, leading to too high densities. Inclusion of the zero point energy and thermal energy corrections results in an expansion of the unit cell. The calculated densities at the pressures used in some experiments are also given. The computed densities with the pressure and thermal energy corrections included show a good agreement with experiments. A final note of more methodological nature: as seen in Table 4, the computed densities based on SPC/CVFF energies (plus SPC/CVFF thermal corrections) agree well with those obtained when B3LYP energies (plus SPC/CVFF thermal corrections) are used. As noted in the Methods, the obtained energy-volume relationships from B3LYP and SPC/CVFF were very similar (even though the two methods predict different energies for the energy-minimum structures as shown in Table 2). This is also reflected in the calculated pressures being similar: for example, the calculated pressure for the geometry-optimized structure of MH-I is 974 MPa when B3LYP energies are used and 989 MPa when purely SPC/CVFF energies are used. Ice-Methane Hydrate Phase Transition. The phase transition in which methane hydrate decomposes and releases methane is of immediate interest to study. The calculated phase transition pressure between ice Ih and MH-I occurs at 9.71 MPa, as seen from Figure 6 (if the ideal gas law or the van der Waals equation with empirical constants43 are used instead of the hard sphere equation of state the computed phase transition occurs at 9.87 and 9.47 MPa, respectively). This is of the same magnitude as the experimental pressure of 2.39 MPa1 at 270.9 K. Possible causes for the overestimation are discussed below. The purely hypothetical phase transitions from ice Ih and methane to MH-II and MH-H are found to be at 45.1 and 116 MPa, respectively. The pressure of the phase transition can be seen from Figure 6 to change by a substantial amount if the free energies are altered by for example a tenth of a kJ. This is really at the limit of our computational accuracy for the electronic energies and forcefield thermal corrections, even though the scaling may manage to compensate for some of the limitations of the functional (such as neglect of dispersion) and basis set effects. A phenomenon that is not included in these calculations is the occupancy of the cages. In the calculations all cages are filled with exactly one methane molecule per cage, while methane hydrates in nature are only partially occupied. The estimated1 fractional occupation of MH-I is 0.87 in small cages and 0.973 in large cages and 0.672 and 0.057 for the small and large cages in MH-II, respectively. A partial occupancy of the cages will not only change the energies but also result in an entropy contribution for the guest molecules due to 6174

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A the different possible distributions among the cages,18 which we have not taken into account here. We have also neglected the entropy due to H-bond proton disorder,25,26,28,53 which could affect the results. A proper inclusion of entropy stemming from partial filling of cages and H-bond proton disorder would probably result in a more realistic description of the methane hydrate-ice equilibrium. This will be the subject of future studies.

4. CONCLUSIONS The structures of methane hydrates have successfully been modeled by combined quantum-chemical and force-field computations. Cell sizes obtained directly from quantum-chemical geometry optimizations are smaller than the experimental cell sizes, but inclusion of zero-point energy and thermal energy corrections from force-field computations as well as pressure effects bring the computed densities in agreement with experiment. Modeling of the transition of ice and methane to methane hydrate was performed and found to occur at 9.7 MPa, which constitutes a slight overestimation compared to experimental data. ’ ASSOCIATED CONTENT

bS

Supporting Information. Crystallographic information files containing the geometry-optimized methane hydrate I, II, and H structures. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author †

E-mail: [email protected].

’ ACKNOWLEDGMENT Financial support from the Swedish Research Council (VR) and grants from the Swedish supercomputer centre, SNAC, are gratefully acknowledged. ’ REFERENCES (1) Sloan., E. D. Clathrate hydrates of natural gases; Dekker: New York, 1990. Electronic reproduction: CRC Press: Boca Raton, FL, 2008. (2) Sloan, E. D., Jr. Nature 2003, 426, 353. (3) Pellenbarg, R. E.; Max, M. D. J. Chem. Educ. 2001, 78, 896. (4) Gbaruko, B. C.; Igwe, J. C.; Gbaruko, P. N.; Nwokeoma, R. C. J. Petrol. Sci. Eng. 2007, 56, 192. (5) Alavi, S.; Woo, T. K. J. Chem. Phys. 2007, 126, 044703. (6) Strobel, T. A.; Koh, C. A.; Sloan, E. D. Fluid Phase Equilib. 2007, 261, 382. (7) Struzhkin, V. V.; Milizer, B.; Mao, W. L.; Mao, H.-K.; Hemley, R. J. Chem. Rev. 2007, 107, 4133. (8) Chou, I.-M.; Sharma, A.; Burruss, R. C.; Shu, J.; Mao, H.-K.; Hemley, R. J.; Gonchrov, A. F.; Stern, L. A.; Kirby, S. H. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 13484. (9) Sun, Q.; Zheng, T.-Y.; Ji, J.-Q.; Wu, X.-Y. J. Chem. Phys. 2005, 122, 024714. (10) Hirai, H.; Uchihara, Y.; Fujihisa, H.; Sakashita, M.; Katoh, E.; Aoki, K.; Nagashima, K.; Yamamoto, Y.; Yagi, T. J. Chem. Phys. 2001, 115, 7066. (11) Shimizu, H.; Kumazaki, T.; Kume, T.; Sasaki, S. J. Phys. Chem. B 2002, 106, 30.

ARTICLE

(12) Hirai, H.; Kondo, T.; Hasegawa, M.; Yagi, T.; Yamamoto, Y.; Komai, T.; Nagashima, K.; Sakashita, M.; Fujihisa, H.; Aoki, K. J. Phys. Chem. B 2000, 104, 1429. (13) Loveday, J. S.; Nelmes, R. J.; Guthrie, M.; Belmonte, S. A.; Allan, D. R.; Klug, D. D.; Tse, J. S.; Handa, Y. P. Nature 2001, 410, 661. (14) Koyama, Y.; Tanaka, H.; Koga, K. J. Chem. Phys. 2005, 122, 074503. (15) Egorov, A. V.; Brodskaya, E. N.; Laaksonen, A. Comput. Mater. Sci. 2006, 36, 166. (16) English, N. J.; Macelroy, J. M. D. J. Comput. Chem. 2003, 24, 1569. (17) Okano, Y.; Yasuoka, K. J. Chem. Phys. 2006, 124, 024510. (18) Belosludov, R. V.; Subbotin, O. S.; Mizuseki, H.; Kawazoe, Y.; Belosludov, V. R. J. Chem. Phys. 2009, 131, 244510. (19) Hermida-Ramon, J. M.; Grana, A. M.; Estevez, C. M. Struct. Chem. 2007, 18, 649. (20) Ikeda, T.; Terakura, K. J. Chem. Phys. 2003, 119, 6784. (21) Gutt, C.; Asmussen, B.; Press, W.; Johnson, M. R.; Handa, Y. P.; Tse, J. S. J. Chem. Phys. 2000, 113, 4713. (22) Rawn, C. J.; Rondinone, A. J.; Chakoumakos, B. C.; Circone, S.; Stern, L. A.; Kirby, S. H.; Ishii, Y. Can. J. Phys. 2003, 81, 431. (23) Udachin, K. A.; Ratcliffe, C. I.; Enright, G. D.; Ripmeester, J. A. Supramol. Chem. 1997, 8, 173. (24) Hirsch, T. K.; Ojam€ae, L. J. Phys. Chem. B 2004, 108, 15856. (25) Knight, C.; Singer, S. J.; Kuo, J-.L.; Hirsch, T. K.; Ojam€ae, L.; Klein, M. L. Phys. Rev. E 2006, 73, 056113. (26) Singer, S. J.; Kuo, J-.L.; Hirsch, T. K.; Knight, C.; Ojam€ae, L.; Klein, M. L. Phys. Rev. Lett. 2005, 94, 135701. (27) Lenz, A.; Ojam€ae, L. Phys. Chem. Chem. Phys. 2005, 7, 1905. (28) McDonald, S.; Ojam€ae, L.; Singer, S. J. J. Phys. Chem. A 1998, 102, 2824. (29) Anick, D. J. J. Mol. Struct. (THEOCHEM) 2002, 587, 87. (30) Anick, D. J. J. Mol. Struct. (THEOCHEM) 2002, 587, 97. (31) Jackson, S. M.; Nield, V. M.; Whitworth, R. W.; Oguro, M.; Wilson, C. C. J. Phys. Chem. B 1997, 101, 6142. (32) Sanders, V. R.; Dovesi, R.; Roetti, C.; Orlando, R.; ZicovichWilson, C.; Harrison, N. M.; Doll, K.; Civeralleri, B.; Buch, I. J.; D’Arco, P.; Llunell, M. CRYSTAL03 User’s Manual; University of Torino: Torino, 2003. (33) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (34) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (35) Hehre, W. J.; Dichfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (36) McQuarrie, D. A.; Simon, J. D. Molecular Thermodynamics; University Science Books: Sausalito, CA, 1999. (37) Accelrys. MS Modeling Getting Started, manual Release 4.0; Accelrys Software Inc.: San Diego, 2006. (38) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1982; p 331. (39) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins Struct. Funct. Genet. 1988, 4, 31. (40) Maple, J.; Dinur, U.; Hagler, A. T. Proc. Natl. Acad. Sci. U. S. A. 1988, 85, 5350. (41) Lau, K. F.; Alper, H. E.; Thacher, T. S.; Stouch, T. R. J. Phys. Chem. 1994, 98, 8785. (42) Whalley, E. J. Chem. Phys. 1984, 81, 4087. (43) Mortimer, R. G. Physical Chemistry; Elsevier: London, 2008. (44) Docherty, H.; Galindo, A.; Vega, C.; Sanz, E. J. Chem. Phys. 2006, 125, 074510. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B. et al. Gaussian 03; Gaussian Inc.: Wallingford, CT, 2004. (46) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (47) Helgerud, M. B.; Waite, W. F.; Kirby, S. H.; Nur, A. J. Geophys. Res. 2009, 114, B02212. (48) Jacobson, L. C.; Hujo, W.; Molinero, V. J. Phys. Chem. B 2009, 113, 10298. 6175

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176

The Journal of Physical Chemistry A

ARTICLE

(49) Dharnawarhana, P. B.; Parrish, W. R.; Sloan, E. D. Ind. Eng. Chem. Fundam. 1980, 19, 410. (50) Bakker, R. J. In Gas hydrates: relevance to world margin stability and climate change; Henriet, J.-P., Mienert, J., Eds.; Geological Society Special Publication No. 137; The Geological Society: London, 1998; pp 75-82. (51) Ogienko, A. G.; Kurnosov, A. V.; Manakov, A. Y.; Larionov, E. G.; Ancharov, A. I.; Sheromov, M. A.; Nesterov, A. N. J. Phys. Chem. B 2006, 110, 2840. € (52) Nordling, C.; Osterman, J. Physics handbook for Science and Engineering; Studentlitteratur: Lund, Sweden, 2004. (53) Rick, S. W.; Freeman, D. L. J. Chem. Phys. 2010, 132, 054509.

6176

dx.doi.org/10.1021/jp111328v |J. Phys. Chem. A 2011, 115, 6169–6176