J. Phys. Chem. 1983, 87,4198-4203
4198
Molecular Dynamics Studies of I c e IC and the Structure I Clathrate Hydrate of Methanet John S. Tse, Mlchael L. Klein,' Dbision of Chemistry, National Research Council of Canade, Ottawa, Canada K1A OR6
and Ian R. McDonald Deparlment of mysicei Chemistry, University In Final Form: January 19, 1983)
of Cambridge, CambrMge, CB2 lEP, United Kingdom (Received: September 24, 1982;
Computer simulationsare used to compare and contrast the dynamical behavior of structure I clathrate hydrates with that of ice IC. The calculations are based on recently proposed pairwise additive intermolecular potentials for the water molecules. The phonon densities of states of ice ICand the hydrates are found to be broadly similar, not withstanding their different crystal structures. This fact explains the similarity of the observed heat capacity and infrared data. In the case of the methane hydrate the simulation predicts a distinct dynamical (localized-mode) behavior for methane molecules trapped in each type of cage.
I. Introduction The clathrate hydrates are a special class of inclusion compounds consisting of water molecules which form an ice-like hydrogen-bonded structure. The empty clathrate has not been observed in nature, but the structure can be stabilized by inclusion of a second component, the guest specie8.l Depending on the size of the guest molecule, two types of hydrate (structure I and 11)are found. Structures I hydrates are formed when the size of the guest is relatively small, as with xenon and methane. For bigger enclathrated molecules, such as propane, a host lattice with larger cavities (structure 11)is required in order to encage the guest species. Both structures have a cubic unit cell containing cavities of two different sizes.2 If the cavities are fully occupied, the composition of the clathrate hydrates would be M.5.75H20 for structure I. More details about their structures and properties are described in several recent review^.^-^ Previous work on clathrate hydrates includes the measurement of dissociation press u r e ~and ~ studies of the infrared spectra,6-8 dielectric rela~ation,~ and nuclear magnetic re~onance.~JO Recent confirmation of extensive accumulation of methane hydrate in permafrost regions and in sediments on the deep ocean floor has revitalized interest in this subject."J2 A t present, most physical properties of the clathrate hydrates can be adequately accounted for by assuming the hydrate to be a nonstoichiometric ice-like solid formed from a mixture of water and a low molecular weight gas or liquid. The main exception is the anomalous behavior of the thermal cond~ctivity.'~Over a temperature range from 100 to 270 K, the thermal conductivities of the clathrate hydrates of ethylene oxide (structure I) and tetrahydrofuran (structure 11)are a factor of five less than that of ice Ih and they also decrease as the temperature is lowered. This behavior would be understandable if the phonon-scattering mechanism in the hydrates were substantially different from that in ice Ih.14 Each of these crystals is based upon a hydrogen-bonded network. The short-range order is similar but the overall structures differ substantially.'"" A complete knowledge of the phonon density of states (frequency spectra) in both the ice and hydrate systems may provide clues in understanding the thermal conductivity anomaly. 'Issued as NRCC No. 21155. 0022-3654/83/2087-4198$01.50/0
TABLE I: Intermolecular Potential Parameters O ( R )= 4€[(U/R)'* - (u/R)6] 4 ~ kJ, mol-' interaction 0 , a 0-0 CH,-0 CH,-H CH,-CH,
3.16 3.33 3.25 3.64
2.60 3.05 1.69 5.46
The purpose of the present paper is to report results obtained by computer simulation for the spectra of empty and filled structure I hydrates and of ice IC. The specific guest molecule we have considered is methane. The calculations were made by the method of molecular dynamics (MD),18the forces between molecules being represented by pairwise additive potentials. In the next section, we describe the potential model we have used and give some details of the MD calculations. Section I11 presents the results and compares them with available experimental data. (1)Van der Waals, J. H.; Plateeuw, J. C. Adu. Chem. Phys. 1959,2, 1. ( 2 ) (a) von Stackelberg, M. Naturwissenschaften 1949,36, 327. (b) von Stackelberg, M.; Miiller, H. R. 2.Electrochem. 1954,58,25. (3)Davidson, D.W. In "Water-A Comprehensive Treaties"; Franks, F., Ed.; Plenum Press: New York, 1973;Vol. 5. (4)Davidson, D. W.; Ripmeester, J. A. J . Glaciol. 1978,21, 33. (5)Glew, D. N. Can. J. Chem. 1960,38, 208. (6)Klug, D.; Whalley, E. Can. J . Chem. 1973,51, 4062. (7) Bertie, J. E.;Jacobs, S. M. Can. J. Chem. 1977,55, 1777. (8)Bertie, J. E.;Jacobs, S. M. J. Chem. Phys. 1979,69, 4105. (9)Gough, S.R.; Davidson, D. W. In "Physics and Chemistry of Ice"; Whalley, E., Jones, S. J., Gold, L. W., Ed.; Royal Society of Canada: Ottawa, 1973. (10)Garg, S.K.; Davidson, D. W. In "Physics and Chemistry of Ice"; Whalley, E., Jones, S. J., Gold, L. W., Ed.; Royal Society of Canada: Ottawa, 1973. (11)Chersky, N.;Makogon, Y. J. Oil GQSZnt. 1970,10, 82. (12) Stoll, R. D.; Bryan, G. M. J . Geophys. Res. 1979,84,1629. (13) (a) Ross, R. G.; Anderson, P.; Backstrom, G. Nature (London) 1981,290,332.(b) Ross, R. G.; Anderson, P. Con. J. Chem. 1982,60,881. (c) Cook, J. B.; Laubitz, M. J. "Proceedings of the 17th International Thermal Conductivity Conference", Gaithersburg, Maryland, June, 1981. (14)Reissland, J. A. In "The Physics of Phonons"; Wiley: London, 1972. (15) Whalley, E.In "The Hydrogen Bond"; Schuster, P., Fundel, T., Sandorfy, C. Ed.; Elsevier: New York, 1976;Vol. 111. (16)McMullan, R. K.; Jeffrey, G. A. J . Chem. Phys. 1965,42,2725. (17)Mak, T. C.W.; McMullan, R. K. J . Chem. Phys. 1965,42,2732. (18)Berne, B. J., Ed. "Modem Theoretical Chemistry", Plenum Press: New York, 1977;Vol. 6.
0 1983 American Chemical Society
The Journal of Physical Chemistry, Vol. 87, No. 21, 1983 4199
MD Studies of Ice IC and Methane Hydrate
TABLE 11: Details of the Molecular Dynamics Calculationsa
r; 2.0 J
B
->> Y
1.0
--
0.0
-
L
!
-
-U,kJ
-
W
Z
-
-1.0
L 3.0
I 3.2
3.4
3.6
3.8
4.0
R
(a)
4.2
4.4
46
4.8
5.0
Figure 1. CH,-H20 Interaction energies for several directions of approach of the methane molecule. The labeling system is the Same as that used by Kolos et ai. in ref 25.
11. Computational Procedure Intermolecular Potential Model. The adequacy of several intermolecular potential models for water molecule in predicting the static properties of different polymorphic forms of ice has recently been evaluated by Morse and Rice.lg Of the models considered, the ab initio pair potential of Matsuoka, Clementi, and Yoshimine (MCY),20 was shown to reproduce best the structures of various ice lattices, but the predicted densities were about 20% too low. In this study, we use the simple point charge (SPC) model for water recently proposed by Berendsen et a1.;21 this was not included in the survey of Morse and Rice.19 The parameters of the model were determined by fitting the internal energy and pressure of liquid water to the experimental values at 300 K, and with the further restriction that the model should yield a good approximation to the measured oxygen-oxygen distribution. The SPC model is a three-center one. Point charges (40 = 4.821e1, qH = 0.4llel) are placed on the atomic sites and the oxygen centers interact through a Lennard-Jones potential, no account is taken of any short-range H-H or 0-H forces. The water molecule is assumed to be rigid, with a 0-H bond strength of 1.0 A and an angle HOH of 109O28’. Details of the potential parameters are given in Table I. We have also made some calculations for the MCY potential, but the results differ little from those for the SPC model and no further mention of them will be made here. Because the MCY potential involves four interaction centers, it is computationally more expensive to study. For simplicity and to minimize the computational efforts, the methane molecules have been treated as spherical particles, the interaction with water being described by a Lennard-Jones potential. The methane potential parameters were estimated from data on krypton;23the similarity of the lattice spacing and binding energy of solid krypton and methane justifies this ~implification.~~ The resulting potential functions are displayed in Figure 1. The plots may be compared with those for an ab initio Ar-H20 po(19)Morse, M. D.;Rice, S. A. J. Chem. Phys. 1982, 76, 650. (20) Matsuoka, 0.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, 64, 1351. (21)Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In “International Forces”; Pullman, B., Ed.; Reidal: Dordrecht, Holland, 1981. (22)Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comp. Phys. 1977., 23. - , 327. -(23)Horton, G. K. In ‘‘Rare Gas Solids”; Klein, M. L., Venables, J. A. Ed.; Academic Press: London, 1976;Vol. 1, p 88. (24)Righini, R.; Maki, K.; Klein, M. L. Chem. Phys. Lett. 1981, 80, 301. ~
system
T, K
p, kbar
mol-’
empty hydrate empty hydrate CH, h y d r a t e CH, h y d r a t e CH, h y d r a t e
110 20 5 89 145 214
-4.0 -1.9 -5.4 -4.0 -2.1
.51.6 48.8 47.6 46.3 44.5
ice IC ice IC
69 252
-3.0
53.8 48.2
0.7
No corrections have been m a d e f o r t r u n c a t i o n of t h e Lennard-Jones 12-6 interactions for distances greater t h a n one-half of t h e simulation b o x . T h e electrostatic interactions are evaluated w i t h o u t t r u n c a t i o n (Ewald m e t h o d ) .
TABLE 111: Positions of Peaks ( c m - I ) in Calculated Frequency Spectraa system
T,K
1
2
3
empty hydrate empty hydrate CH, h y d r a t e CH, h y d r a t e ice IC
110
47 44 53
71 65 69 65 67
89 79 88 88
205 145 214 69
47 44
a T h e labels 1, 2, a n d 3 refer t o t h e peaks arrowed in Figures 7 a n d 8.
tential.25 It is gratifying that both the profile and orientation dependence shown in the ab initio calculations were well-reproduced by our simple potential model. Molecular Dynamics Calculations. The technical details of the MD calculations are the same as those in our earlier work on liquid water.26 The intermolecular forces and energies were evaluated by the Ewald method,27and the equations of motion of the water molecules were solved numerically by the method of constraints developed by Ryckaert et a1.22 Simulations of structure I hydrates were carried out for a single cubic unit cell of length 12.03 A16 containing 46 water molecules, with periodic boundary conditions imposed. The initial positions of the oxygen atoms were determined from the X-ray structure of ethylene oxide hydrate at -25 OC.16 For the MD calculations on ice IC, a larger system was used, consisting of 64 water molecules in a cubic box of length 12.70 A. In this case, the initial oxygen coordinates were taken from an X-ray diffraction study.28 As the protons of the water molecules in both the clathrate hydrate and ice ICare rotationally disordered, a proton assignment is also required, and initial proton coordinates were therefore selected in a random fashion to conform to the Bernal-Fowler rules.29 In the initial configuration the protons were always placed along the 0-0 bonds. In general, this procedure leads to water molecules with HOH angles slightly different from that required by the potential models. However, the method of constraints22that we use to solve the equations of motions immediately “repairs” the water molecules at the first step of the calculation and subsequent configurations and therefore based upon water molecules with the correct geometry. Molecular dynamics calculations were performed for the empty hydrate, hydrate fully occupied with methane, and (25)Kolos, W.;Corongiu, G.; Clementi, E. Znt. J. Quantum Chem. 1980, 17, 775. (26)Impey, R. W.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1981, 74. 647. (27)Cowley, E.R.; Jaccuci, G.; Klein, M. L.; McDonald, I. R. Phys. Reu. B 1976, 14, 1758. (28)Honjo, G.; Shimaoka, K. Acta Crystallogr. 1957, 10,710. (29) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, I, 515. r
-
4200
The Journal of phvsical Chemistry, Vol. 87, No. 21, 1983
4.0
I
0.0
a I
0.0
.
10 .
,
,
2.0
1
Tse et
al.
l IL
3.0
5.0
4.0
6.0
Figure 2. goo(R) for (a) ice IC at 69 K, (b) empty hydrate at 110 K, and (c)methane hydrate at 145 K.
the proton-ordered ice ICat several temperatures without varying the unit cell dimensions. In each simulation the trajectories of the molecules were followed for a total of 30 ps, about one-third of which was allowed to equilibration. In all cases, no drift in the calculated energy and pressure was detected during the last 20 ps of the run, indicating that the clathrate structures remained at least metastable over the period of the simulations. A summary of the MD data (temperatures, pressures, and configurational internal energies) is given in Table 11. Since we did not evaluate the free energy no conclusion can be drawn concerning the relative stability of the crystals we have studied. However, we note that at the lowest temperatures the enthalpies ( E + p v ) of all three systems are rather similar. The negative pressures observed in the calculations (see Table 111) are always small with respect to the bulk modulus of icem ( 100 kbar) and hence do not lead to mechanical instabilities of the crystals. This result arises partly from inadequacies in the potential model but equally significant is the neglect of the zeropoint energy contribution to the pressure at low temperatures. It should also be noted that we make no allowance for the effect of temperature on the unit cell dimensions. The power spectrum f ( u ) of the atomic motions was obtained by Fourier transformation of the velocity autocorrelation functions $(t), for 0, H, and CH4, i.e.
Figure 3. g&R) for (a) ice IC at 69 K, (b) empty hydrate at 110 K, and (c) methane hydrate at 145 K.
1 I
8.0
I
t
0.01 0.0
I
'
'
'
1.0
'
J
2.0
',I
I
,"
3.0
'
'
'
50
4.0
'
'
6.0
I
RH,(& Flgure 4. gHH(R) for (a) ice IC at 69 K, (b) empty hydrate at 110 K, and (c) methane hydrate at 145 K.
N
f(u) =
x m $ ( t cos ) ut dt
with
W) = (v(O)*v(t))/ (V(O)*V(O) ) where v(t) is the velocity of a particle at time t. For a monoatomic harmonic crystal, the power spectrum f ( u ) is just the phonon density of states. In the present case the proton and oxygen power spectra will be dominated respectively by the librational and translational motions of the water molecules. The power spectrum associated with the motion of the methane molecules will, of course, only contain translational contributions. 111. Results and Discussion Atomic Pair Distribution Functions. The calculated atomic pair distribution functions goo(R),goH(R),and gHH(R) for ice IC, empty hydrate, and CH, hydrate are compared in Figures 2-4. The positions of the peaks are dictated by the respective crystal structures. The distribution functions have been calculated out to half the (30) Gold, L. W. Can. J. Phys. 1958, 36, 1265.
length of the molecular dynamics cube edge, i.e., R = 6.015 = 6.35 A for ice IC. The first maximum in goo(R) for ice IC,located at 2.77 A, corresponds to the nearest-neighbor 0-0 distance.28 The distribution of the first shell of oxygen atoms in the hydrates is more asymmetric; there are four nonequivalent 0-0 distances ranging from 2.766 to 2.844 A.16 The most distinctive feature in the results for the hydrates is the disappearance of the third coordination shell observed in ice ICat about 5.25 A. This occurs because of the formation of cavities in the hydrate. The functions goH(R)and gHH(R)are very similar in the three systems. From the distribution functions we conclude that the various structures are indeed stable over the length of our MD simulations and that the presence of the methane molecule in the methane hydrate has little direct influence on the overall structure. The calculated methane-water distribution functions are depicted in Figure 5. There is a striking similarity in gcb4(R) between the present calculation and a previous Monte-Carlo study of a solution of methane in water.31 Integration to the first minimum of the clathrate results yields a coordination number of 23 which corresponds to the mean number of water molecules making up a hydrate cavity, averaged over both large and small cages. What is surprising is that the Monte-Carlo distribution function yields almost exactly the same coordination number. The CH,-0 distances in the hydrate range from 3.2 to 5.7 8, with a peak at 4.1 8, compared with a range of 3.1 to about 6.0 8, in the
A for the clathrate hydrates and R
~
~~
(31)(a) Owicki, J. C.; Scheraga, H. A. J . Am. Chem. SOC.1977, 99, 7413. (b) Swaminathan, S.;Harrison, S. W.; Beveridge, D. L.J . Am. Chem. SOC.1976, 100, 5705.
The Journal of Physical Chemistry. Vol. 87, No. 21, 1983 4201
MD Studies of Ice IC and Methane Hydrate
4.0
-
3.0
-
20
-
1.0
-
a
b
theory ( 2 5 2 K )
Figure 6. Translational frequency spectra for ice. (a) Experimental result Ih from ref 38 at 261 K. (b) MD result for ice IC obtained by using the SPC model at 252 K. 0.01 0.0
*
I
1.0
"
'
1 ' '
I
4.0
'
I
5.0
'
R%) Flgure 5. Radial distribution function for (a) CH4-O and (b) CH4-H in methane hydrate at 145 K.
The profile of gCH,-&3) closely resembles that of the CH4-O distribution. The distance of closet approach is only about 0.2 A shorter than for CH4-0, suggesting that on average the cage contains no proton pointing radially inward. A similar conclusion was reached in a Monte-Carlo simulation of argon in liquid water.32 Frequency Spectra. As mentioned earlier, the power spectra of the proton and oxygen motions are dominated, respectively, by the librations and translations of the water molecules. In all our calculations, the translational and librational modes are separated by a well-defined frequency gap, around 400 cm-'. We therefore discuss each type of motion separately. The translational density of states for ice IChas not yet been determined experimentally. Careful comparative s t ~ d i e s ~ of ~ - the ~ ' infrared spectra of ice Ih and ice IC suggest that their frequency distribution functions should be similar. It is therefore reasonable to compare the overall theoretical spectrum for ice ICwith the inelastic neutron data for ice Ih,= and this is done in Figure 6. As the figure shows, the two band profiles are in fairly good agreement over the entire range of frequency. The discrete peak structure of the low-frequency region of the calculations is a direct consequence of the small number of unit cells (8) employed. However, we have been able to show by evaluation of the dynamical structure factor that the prominent peaks at 29, 44, and 67 cm-l are due to zone boundary transverse acoustic phonons propagating along the high symmetry directions. These peaks should persist even if we could perform calculations on a larger number of unit cells. The only change that we anticipate would be a filling up of the gaps in the low-frequency region. With this in mind, the calculated transverse acoustic band (32) Alagona, G.;Tani, A. J. Chem. Phys. 1980, 72, 580. (33) Prask, H.;Boutin, H.; Yip, S. J. Chem. Phys. 1968, 48, 3367. (34) Bertie, J. E.; Whalley, E. J . Chem. Phys. 1964, 40, 1637. (35) Bertie, J. E.; Whalley, E. J. Chem. Phys. 1967, 46, 1271. (36) Bertie, J. E.;Whalley, E. J. Chem. Phys. 1969, 50, 4501. (37) Bertie, J. E.; Jacobs, S. M. J . Chem. Phys., 1977, 67, 2445. (38) Renker, K.B.; Blanckenhagen, P. V. In 'Physics of Ice"; Riehl, N., Bullemer, B., Engelhardt, H., Ed.; Plenum Press: New York, 1969.
cm-1 Flgure 7. Translational frequency spectra for (a) ice IC at 69 K, (b) empty hydrate at 110 K, and (c) methane hydrate at 145 K. Peaks marked with an arrow are discussed in the text.
occurs at a lower frequency than found experimentally. As with the negative pressures discussed earlier, this may indicate an inadequacy of the potential model. The calculated translational band extends up to about 325 cm-1.39 Earlier lattice dynamics calculations employing empirical force constants39 did not find any translational phonons having frequencies above about 275 cm-'. This failure has been attributed to the neglect of long-range forces, presumably electrostatic in origin, in the force-field models.35 It is clear from the present calculations that a more realistic model of the molecular interactions is capable of removing this discrepancy. The MD calculations also predict correctly that the phonon density of states for the low-frequency region (below 100 cm-') should be more intense than at higher frequencies (100-35O cm-I). Earlier calculations employing empirical force fields predicted the opposite. The calculated spectrum of Figure 6 contains a number of other resolved peaks. These correspond to phonons at the Brillouin zone b o u n d a r i e ~ . We ~ ~ ,tentatively ~~ assign (39) Prask, H.J.; Trevino, S. F.; Gault, J. D.; Logan, K. W. J. Chem. Phys. 1972, 56, 3217.
4202
The Journal of Physical Chemistty, Vol. 87, No. 21, 1983
Tse et al.
1.0
0.8
t
,
v)
I
w
+
g
"'I
0.6
%
> t-
iij 0.4
z w
a 0.2
i
0.0 0
100
200
300 0-
100
200
300
cm 1 Figure 8. Translaticoal frequency spectra for (a) empty hybate at 205 K and (b) methane hydrate at 214 K. Peaks marked with an arrow are discussed in the text.
the band centered a t approximately 240 cm-' to longitudinal and transverse optical modes; those at 140 and 180 cm-' are due to the longitudinal acoustic modes. The translational density of states for the empty and full clathrate hydrates are shown in Figure 7. The main features are similar to those found in ice ICat a comparable temperature (alsoshown in Figure 7). The most noticeable difference on passing from ice to clathrates is that the intensity of the peak at about 30 cm-' (labeled A in Figure 7) was drastically reduced. This behavior may be an artifact of the MD calculations. Variations in temperature lead to substantial redistribution of the density of states in the hydrates, as can be seen by comparison of Figures 7 and 8. In particular, there is an enhancement of the density of states in both the empty and full clathrate in the 100-300-cm-' region. In addition, there is a tendency for the vibrations to shift to lower frequencies as the temperature increases (see Table III). The same effect has been observed in infrared spectra of ice IC and ice Ih taken a t different temperature^.^^^^^ The results of Table 111also show that between ice and clathrate there is a significant upward shift in frequency of the peak labeled 1in Figures 7 and 8, and a further shift on inclusion of the guest molecules. Enclathration also causes the peak labeled 2 to broaden (see Figure 8),so that peak 3 is almost lost, but the frequencies of peaks 2 and 3 are almost unaltered. The frequency spectrum of the encaged methane molecules at 145 K is shown in Figure 9a. Three peaks, at 35,54, and 78 cm-', can be identified. Since there are two types of cages of different sizes in the structure I hydrate, we can further partition the vibrational spectrum into the contribution from the small cages3 (mean radius 4.0A) and from the large cages3 (mean radius 4.3 A). The results are shown in Figure 9 parts c and b, respectively. This partitioning makes it possible to assign the highest energy vibration (78 cm-') to the motion of methane molecules in the almost spherical small cages;41the more hindered environment in the smaller cavity increases the frequency of the vibration. The lower vibrational frequencies and the splitting of the translational energy of
- -
(40) Klein, M. L. In 'Computer Modeling of Matter";Lykos,P., Ed.; American Chemical Society: Washington, DC, 1978; ACS Symposium Series, No. 86,p 94. (41) Burgiel, J. C.;Meyer, H.; Richards, P. L. J.Chem. Phys. 1965,43, 4291.
0.0
lL--LL. cm-1
Flgure 9. Frequency spectra for CH, in methane hydrate at 145 K: (a) total spectrum, (b) large cage, and (c) small cage. r
I
I
1.0
0.8
v)
k W
E
8 0.6 >.
5z
0.4
0.2
0.0
0
Flgure 10. Librational frequency spectra for (a) ice IC at 69 K, (b) empty hydrate at 110 K, and (c) methane hydrate at 145 K.
the methane molecules in the large cages can be rationalized in terms of the greater size and asymmetry of the cavity. The methane frequencies fall in the same range as the prominent peaks visible below 100 cm-' in Figures 7 and 8. This is suggestive of a coupling between the translational motion of the guest molecules and that of the host lattice. There are no experimental infrared data available for methane hydrates. Recently, Bertie et al.42reported a low-temperature infrared spectrum for the xenon hydrate between 50 and 350 cm-l. In their main features, the results are very similar to those for ice IC,though it is apparent from the reported frequencies that there is a slight shift to higher energy from ice Ic.37342This is in general accord with the resulta of the presence calculations. The calculated frequency spectrum for the protons are shown in Figure 10; they span the range between 450 and 1000 cm-'. The proton motions can be studied experimentally by incoherent neutron scattering techniques.3yB,a For ice Ih, the neutron time-of-flight cross section shows a broad asymmetric profile with total width of about 600 (42) Bertie, J. E.; Jacobs, S. M. J. Chem. Phys. 1982, 77, 3230. (43) Harling, 0. K. J. Chem. Phys. 1969, 50, 5279.
MD Studies of Ice
The Journal of Physical Chemisfry, Vol. 87, No. 21, 1983 4203
IC and Methane Hydrate
TABLE IV: Heat Capacity (C,) for Ice IC and Clathrate Hydrates (J K-'m o l - ' ) system empty hydrate empty hydrate CH, hydrate CH,hydrate ice IC ice IC
trans- librap tionO guest
totalb
exptC
0.7 6.9 2.4 7.4 0.0 9.8
17.1 28.2 20.9 28.8 11.5 31.8
17.2 28.5 21.2 29.6 11.4 34.0
T, K lation= 110 16.4 205 21.3 145 18.5 214 21.3 69 11.5 252 22.0
24.2 24.5
a Translational and librational contributions are relative t o 1 mol of water. Contribution from guest species has Estimated from C, of ice Ih; see text. been omitted.
cm-' and a maximum located at 570 cm-1.39 This agrees very well with the results of Figure 10. The calculated spectrum for ice IC shows a narrow band peaked at 560 cm-' and there are two lower intensity bands centered at 680 and 175 cm-'. The neutron experiments were not of sufficient resolution to detect any such bands, but the main features of our spectrum appear also in ~ a l c u l a t i o n for s~~ ice ICbased on a force-constant model due to Shawyer and Dean.47 The proton frequency spectra of the empty and full clathrate hydrates do not resemble ice ICas closely as their translational counterparts. The first peak in the spectrum is shifted from 560 to 516 and 534 cm-l for the empty and filled hydrate, respectively. Such shifts are presumably due to a weaker hydrogen bond network in the hydrates relative to ice. However, we also observe a change in the profile. In particular, the density of states in the region around 600 cm-' increases markedly in the hydrates relative to ice IC. At least part of this effect arises from temperature differences in the three simulations shown in Figure 10. Heat Capacity Calculations. One interesting application of the frequency distribution function is the calculation of the heat capacity. The heat capacity at constant volume (C,) can be evaluated by convoluting the calculated power spectra with an Einstein function. This approach provides an approximate way of allowing for quantum effects. There are two components: the center-of-mass translation and the molecular libration contributions. In our calculations, each of the phonon distribution functions is normalized to 3N, where N is the number of molecules ~~
~~~~~
(44) Giauque, W. F.; Stout, J. W. J. Am. Chem. SOC.1936,58, 1144. (45) Leadbetter, A. J. h o c . R. SOC.London, Ser. A 1965, 287, 403. (46) Leaist, D.; Murray, J. J.; Post,M. L.; Davidson, D. W. J.Phys. Chem. 1982,86,4175. (47) Shawyer, R. C.; Dean, P. Discuss. Faraday SOC.1969, 48, 102.
in the simulation box. The calculated heat capacities for ice IC, empty hydrate, and methane hydrate at several temperatures are given in Table IV. The experimental values were estimated from the tabulated C, data for ice Ih,44145and a recent measurement on ethylene oxide hydrate.& In view of the limitations inherent in the method, the agreement between the theoretical and the measured results may be regarded as satisfactory. The heat capacity of the methane molecules remains virtually unchanged with variation in temperature and is close to the ideal gas value of 24.9 J K-'mol-l. Thus the temperature dependence of the hydrate heat capacity is due almost entirely to the temperature dependence of the lattice heat capacity. The constancy of the heat capacity for the molecule entrapped in the hydrate has also been observed in experiment~.~~
IV. Conclusion In this paper, we have demonstrated that the SPC model for the interaction of water molecules is able to describe many of the observed dynamical properties of ice ICand clathrate hydrates. From the numerical results presented here, it appears that the phonon densities of states of ice and the hydrates are broadly similar. This explains the close similarity of their infrared spectra and heat capacities. However, there are differences in detail in the densities of states, and these should be observable experimentally. Although the interaction between the encaged molecules and the water lattice is small, the guest species induces frequency shifts and enhancement of the phonon density of selective vibrational modes. This kind of weak coupling does not have a large effect on thermodynamic properties of the hydrate, but it may radically alter the phonon-scattering mechanism, leading to very different thermal conductivities from those of the ices. Confirmation of this speculation must await further work. In the case of methane hydrate we have found the presence of low-frequency localized modes that are associated with the translational motion of the methane molecules in their respective cages. Not surprisingly, the higher frequency motion is associated with methane molecules in the small cages. These predicted translation modes should be observable in infrared studies but to date no low-frequency data have been reported. Acknowledgment. We thank Don Davidson for introducing us to the subject of hydrates and several discussions on their properties, and Roger Impey for technical assistance. Registry No. Ice, 7732-18-5.
