Diffusion in Solids - Handwritten Notes | MCDB 220C, Papers of Biology

Download Diffusion in Solids - Handwritten Notes | MCDB 220C and more Papers Biology in PDF only on Docsity!

For, methane production, hydrate mining, and CO 2 sequestration applications, the time scale for methane diffusion is an important variable. Experiments on hydrates are difficult for several reasons.1,9-^12 The foremost difficulty in measuring diffusion rates within a hydrate crystal is that large single crystals of methane hydrate are difficult to grow and handle in the laboratory and that grain size distributions in polycrystalline hydrates are generally sensitive to conditions of growth.13,14^ Indirect estimates of methane diffusion rates have been made by fitting models of hydrate growth.^9 -^11 However, apparent mass transfer rates in polycrystalline materials are sensitive to factors like sample grinding and grain morphology.^15 Grinding could affect both grain sizes and the nature of grain boundaries, so it is difficult to separate contributions to the apparent rate of mass transfer through polycrystalline hydrates. The apparent rate of mass transfer into a polycrystalline sample could be limited by diffusion along grain boundaries, by diffusion into the hydrate lattice, or by crystal growth kinetics in a mechanism where the CH 4 hydrate decomposes and a CO 2 hydrate is reconstructed.^15 Understanding the mass transfer and exchange mechanisms would resolve critical questions that surround potential carbon sequestration and methane recovery strategies. For example, Ripmeester and co-workers observe that CH 4 remains in the hydrate even as CH 4 becomes infinitely dilute in the atmo- sphere.^15 Thermodynamics predicts an ever-increasing chemical potential driving force (∆μCH4,hydrate ) ∆μ^0 CH4,hydrate + RT ln P CH4) that pushes CH 4 out of the hydrate as the surrounding atmosphere becomes CH 4 depleted. Thus, mass transfer limita- tions may be important in hydrates. This paper investigates one of the possible diffusion mech- anisms that might explain observations like that of Ripmeester and co-workers.^15 We use an equilibrium path sampling method based on Bolas,^16 reactive flux simulations,^17 and kinetic Monte Carlo^18 to obtain the self-diffusion constant of the methane molecule in a natural gas hydrate. Our hypothesized mechanism for diffusion is thermally activated hopping from an occupied cage to an adjacent vacant cage. As proposed previously for the CO 2 diffusion mechanism in hydrates,^19 we hypothesize that a water vacancy between the donor and acceptor cages facilitates methane hopping. Our hypothesis implies a diffusion rate that depends on temperature and methane vacancy concentration. One consequence of concentration dependent diffusivity is the possibility of singular diffusion where wetting-front behavior could arrest the mass transfer process in a nonequilibrium state.^20 We emphasize that further simulations and experiments are also needed to investigate alternative diffusion mechanisms like Bjerrum defect assisted diffusion21,22^ and crystal decomposition-

reconstruction mechanisms like those involved in hydrate formation from ice.^9 -^11 Methodology Force Field Parameters and Simulation Details. Previous simulations have used ab initio calculations^23 -^26 and empirical force fields8,19,27-^35 to model hydrates. TIP4P water^36 describes the known phases of ice,37,38^ and recent studies show that it accurately describes the host lattice of the hydrate.8,28,36^ The TIP4P potential was used for the internal degrees of freedom of water and for water-water interactions. Long range electrostatics were treated using particle mesh Ewald summation with 32 k-points.^39 Internal degrees of freedom in CH 4 and CH 4 - CH 4 interactions were modeled using the OPLS potential.^40 The CH 4 and H 2 O interactions were specially developed for methane hydrates by Anderson et al.24,28^ All simulations were performed in CHARMM on a cubic system of 8 unit cells (2 × 2 × 2) of a structure I methane hydrate. Water orientations were selected with the algorithm of Buch et al.^41 Simulations were performed at two temperatures, 225 and 250 K, using the Nose-Hoover thermostat.^42 To approximate a constant pressure of 40 atm, the average box length for a fully occupied structure I methane hydrate (CH 4 · 6H 2 O) was computed for a 1 ns trajectory with an extended Lagrangian barostat.^43 We performed simulations to calibrate the size of a cubic box as a function of temperature. The hydrate remains stable for a nanosecond at temperatures as high as 325 K at 40 atm, but these high temperature structures are probably metastable. We chose the temperatures 225 and 250 K which are slightly below the melting temperature of TIP4P water based on the work of Vega et al.^38 At P ) 40 atm the box size is 23.33 ( 0.05 Å at T ) 225 K and 23.39 ( 0.05 Å at T ) 250 K, in reasonable agreement with the experimental lattice constants.^1 Subsequent calculations to compute rate constants were per- formed with NVT simulations using a cubic simulation box at the mean length from the NPT simulations at the corresponding temperature. Thermal expansion was incorporated in the simulation box size, ensemble effects generally diminish with system size,^44 and solids have small compressibility, so the free energy barriers computed within the NVT ensemble were interpreted as ap- proximate^44 Gibbs free energies.

(12) Bach-Verges, M.; Kitchin, S. J.; Harris, K. D. M.; Zugic, M.; Koh, C. A. J. Phys. Chem. B 2001 , 105 , 2699. (13) Udachin, K. A.; Lu, H.; Enright, G. D.; Ratcliffe, C. I.; Ripmeester, J. A.; Chapman, N. R.; Spence, G. Angew. Chem., Int. Ed. 2007 , 46 ,

(14) Klapp, S. A.; Klein, H.; Kuhs, W. F. Geophys. Res. Lett. 2007 , 34 , L13608. (15) Lee, H.; Seo, Y.; Seo, Y.-T.; Moudrakovski, I. L.; Ripmeester, J. A. Angew. Chem., Int. Ed. 2003 , 42 , 5048. (16) Radhakrishnan, R.; Schlick, T. J. Chem. Phys. 2004 , 121 , 2436. (17) Chandler, D. J. Chem. Phys. 1978 , 68 , 2959. (18) Gillespie, D. T. J. Comp. Phys. 1978 , 28 , 395. (19) Demurov, A.; Radhakrishnan, R.; Trout, B. L. J. Chem. Phys. 2002 , 116 , 702. (20) Shampine, L. F. Quart. App. Math. 1973 , Oct , 287. (21) Davidson D. W. and Ripmeester, J. A. In Inclusion Compounds, Vol. 3 ; Atwood, J. L., Davies, J. E. D., MacNicol, D. D., Eds.; London Academic: Oxford, 1984; p 69.

(22) Kirschgen, T. M.; Zeidler, M. D.; Geil, B.; Fujara, F. Phys. Chem. Chem. Phys. 2003 , 5 , 5247–5252. (23) Anderson, B. J. Fluid Phase Eqilib. 2007 , 254 , 144. (24) Anderson, B. J.; Tester, J. W.; Trout, B. L. J. Phys. Chem. B 2004 , 108 , 18705. (25) Klauda, J. B.; Sandler, S. I. J. Phys. Chem. B 2002 , 106 , 5722. (26) Alavi, S.; Ripmeester, J. A. Angew. Chem., Int. Ed. 2007 , 46 , 6102. (27) Rodger, P. M. J. Phys. Chem. 1990 , 94 , 6080–. (28) Anderson, B. J.; Tester, J. W.; Borghi, G. P.; Trout, B. L. J. Am. Chem. Soc. 2005 , 127 , 17852. (29) Radhakrishnan, R.; Trout, B. L. J. Chem. Phys. 2002 , 117 , 1786. (30) Baez, L. A.; Clancy, P. Ann. N.Y. Acad. Sci. 1994 , 715 , 177. (31) English, N. J.; MacElroy, J. M. D. J. Comput. Chem. 2003 , 24 , 1569. (32) Guo, G. J.; Zhang, Y. G.; Liu, H. J. Phys. Chem. C 2007 , 111 , 2595. (33) Tse, J. S.; Klein, M. L.; McDonald, I. R. J. Chem. Phys. 1984 , 81 ,

(34) Handa, Y. P.; Tse, J. S.; Klug, D. D. J. Chem. Phys. 1991 , 94 , 623. (35) Shpakov, V. P.; Tse, J. S.; Tulk, C. A.; Kvamme, B.; Belosludov, V. R. Chem. Phys. Lett. 1998 , 282 , 107. (36) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983 , 79 , 926. (37) Sanz, E.; Vega, C.; Abascal, J. L. F.; MacDowell, L. G. Phys. Re V_. Lett._ 2004 , 92 , 255701. (38) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005 , 122 , 114507. (39) Darden, T.; York, D.; Pederson, L. J. Chem. Phys. 1993 , 98 , 10089. (40) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984 , 106 , 6638. (41) Buch, V.; Sandler, P.; Sadlej, J. J. Phys. Chem. B 1998 , 102 , 8641. (42) Hoover, W. G. Phys. Re V_. A_ 1985 , 31 , 1695. (43) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994 , 101 ,

(44) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications , 2nd ed.; Academic Press: San Diego, 2002.

Path Sampling Study of Methane Diffusion in Hydrates A R T I C L E S

Finally, we note that lattice distortions associated with the transition state along a hopping pathway could be affected by their own periodic images. Larger scale simulations, e.g., a 27 unit cell simulation (3 × 3 × 3), would be needed to determine whether finite size effects are important. Spherical Bipolars for Reaction Coordinates and for Detecting Water Vacancies. If methane diffuses by hopping between adjacent cages of the hydrate lattice, then there are three pathways to consider between adjacent cages in structure I.^1 These are as follows: “L6L” from large donor cage to a large acceptor cage through a water vacancy in a six-membered water ring “S5L” from small donor cage to large acceptor cage through a water vacancy in a five-membered water ring “L5L” from large donor cage to large acceptor cage through a water vacancy in a five-membered water ring Figure 1 shows configurations along the S5L pathway with the methane in the donor (small) cage, between the donor and acceptor cages, and in the acceptor (large) cage. Bipolar coordinates ( q , θ) as shown in Figure 2 were used to describe the position of the methane molecule relative to the centers of the donor and acceptor cages. The reaction coordinate is q , and θ is a measure of deviation from the straight line path between the initial and final cage centers. These coordinates are convenient because they describe all three pathways (S5L, L5L, and L6L) between cages in the structure I hydrate lattice. For each configuration sampled, the donor and acceptor cage centers are needed to compute q and θ. The cage centers were found

using ghost atoms that avoid all atoms except for the fi V e atoms of the hopping, or “tagged”, methane molecule. The ghost potential V( r ), defined in Figure 3, is truncated at the minimum to leave a repulsive 1/ r^2 core that avoids any atom within 5 Å (except for the tagged methane). Because the cages are 3-4 Å in radius, the surrounding water molecules “push” the ghost atoms to the center of the acceptor and donor cages. The ghost particles do not alter the dynamics or energetics of the actual methane hydrate system. They are a fictitious device for obtaining the acceptor and donor cage centers as foci for the bipolar coordinates. For each config- uration sampled, the positions r g1 and r g2 of two ghost atoms were optimized using the potential in eq 1a.

V ( r g 1 , r g 2 , | rN ) ) ∑

untagged k

V(| r k - r g 1 |) +

untagged k

V(| r k - r g 2 |) +

(| r g 1 - r g 2 | - r EQ )^2 (1a)

and

V( r ) ) {5 ⁄^ r

(^2) - 2 ⁄ r + 1 ⁄ 5 r < 5 0 r g 5

(1b)

The last term in the ghost potential restrains the distance between the acceptor and donor cage centers near their equilibrium distances. The restraint helps prevent one ghost atom from hopping into a neighboring cage during the simulations. The equilibrium distances r (^) EQ in the restraints were the average distances from an unbiased simulation between two unrestrained ghost atoms in adjacent acceptor-donor cages with a water vacancy between the cages. The equilibrium distances for the L6L, L5L, and S5L cages are shown in Table 1. Table 1 also includes two values of q ( q min and q max ) and one value of θ (θmin ) for each hopping pathway. These values specify a region between donor and acceptor cages as shown in Figure 4. Our hypothesis requires a water vacancy between the donor and acceptor cages, so one water molecule must be missing from the region of Figure 4. Note that the region is actually a volume of rotation about the donor-acceptor axis. For the L6L hop, the region should contain 6 - 1 ) 5 water molecules, and for L5L or S5L hops, the region should contain 5 - 1 ) 4 water molecules. Testing the water vacancy hypothesis requires a sampling method that can include a constraint on the number of water molecules in the water vacancy region. Furthermore, q and θ are iteratively defined variables that are obtained through a numerical optimization of ghost atom positions. Importance sampling by simple molecular dynamics would be convenient, but forces corresponding to the constraints and the q and θ variables are not easily computed. To

Figure 1. S5L hop showing the methane molecule (black) in the donor cage (S, top), between the cages (middle), and in the acceptor cage (L, bottom). Hydrogen atoms on the methane were omitted for clarity. Also note the missing water molecule between the two cages.

Figure 2. Coordinates q ) ln[ R D/ R A] and θ ) angle[DMA] describe the position of a methane molecule (M) with respect to ghost atoms A and D that occupy the initial (D) and final (A) cages. The plane between the two cages is the isosurface q ) 0. Contours of θ are shown in 30° intervals, and contours of q are shown in intervals of 1 / 2.

Figure 3. Ghost atom potential V( r ) used to obtain the initial and final cage centers as foci of the bipolar coordinate system.

Table 1. Equilibrium Distance between Ghost Atoms and the Values of q and θ That Define a Region That Must Contain a Water Vacancy water vacancy region r EQ/Å q min q max θmin L6L 5.68 - 0.32 0.32 85 ° L5L 7.14 - 0.32 0.32 95 ° S5L 6.59 - 0.32 0.16 90 °

A R T I C L E S Peters et al.

energy reaches a minimum and then rises again at large values of q because of the diminishing volume between q and q + dq at large values of | q |. Note that Alavi and Ripmeester^26 found that H 2 easily passes through five- and six-membered water rings even when there is no water vacancy in the structure. A comparison of these results suggests that guest molecule size has a large effect on the diffusion mechanism and on the diffusion rate. Concentration of Water Vacancies. Our hypothesis is that water vacancies facilitate diffusion. To test that hypothesis, we consider two pathways by which methane can hop from a donor cage to a vacant “acceptor” cage in the hydrate lattice: (1) The methane plows through an intact clathrate cage. This mechanism corresponds to the free energy of Figures 5 and 6. (2) The methane “waits” until a water vacancy appears between the donor and acceptor cages and then slips through the defect site with a reduced barrier. For the second pathway, the rate of passage through a water vacancy defect must be multiplied by the fraction of time that the water vacancy is present in the donor-acceptor channel. The free energy to form a water vacancy is related to the fraction of time that a water vacancy defect resides in the donor-acceptor channel. Thermodynamic integration,^44

∆ F ) (^) ∫ 0

1 dλ〈∂ E ∂λ 〉λ

was used to compute the reversible work to create a water vacancy by scaling the partial charges and Lennard-Jones interaction parameters of one TIP4P water molecule by a factor λ. Demurov et al.^19 computed the concentration of water vacancies in a type-I hydrate using SPC/E water and by transferring the water molecule to an ideal gas phase. Our results are expected to differ from theirs because we have used TIP4P water and transferred the water to the stable phase of TIP4P

water at the temperature and pressure of the simulation. The end states of thermodynamic integration are depicted in Figure 7. From the work of Vega et al.^38 the melting point of TIP4P water at 40 atm can be estimated as 231 K ( 5 K. Thus, at 250 K the TIP4P water molecule was transferred from the hydrate phase to a TIP4P water bath at 250 K and 40 atm. At 225 K where TIP4P is stable as ice, the calculation is more complicated. The calculation is not intended to compare the free energy to form a water V acancy in the hydrate with the free energy to fill a water V acancy in ice. Therefore, the TIP4P water molecule was transferred from a hydrate at 225 K to a water bath at 232 K. The result was corrected by the free energy change to freeze (∆ G ) 0) and cool the system with one additional water from 232 to 225 K. ( S ) S o^ + Cp ln[ T / T o] with d G ) - S d T ) For the correction, experimental values at 1 atm were used for the heat capacity and for S o.^56 The difference between two Poynting corrections for pressurization and depres- surization to and from 40 atm at 225 and 232 K were ignored. The estimated free energy correction was nearly insignificant, 0.286 J/mol. Results are shown in Figure 8. At 250 K and 40 atm thermodynamic integration predicts that approximately 1 water molecule in the hydrate lattice per 74 is vacant. Further within the region of hydrate stability at 225 K and 40 atm, approximately 1 water in 211 is vacant. The 4.3 kT and 5.3 kT free energies required to create a vacancy in the hydrate at 250 and 225 K must be included in comparing the barriers to methane diffusion via the water vacancy assisted and unassisted mechanisms. Free Energy Barriers with a Water Vacancy. Figure 9 shows (54) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; the free energy barriers for methane hopping with a water Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005 , 26 , 1781. (55) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992 , 13 , 1011.

(56) Lide, D. R. Handbook of Chemistry and Physics , 88th ed.; CRC Press: Boca Raton, FL, 2007.

Figure 5. Free energy landscape for L6L hop with no water vacancy. q -Isosurfaces are shown decreasing and increasing in intervals of 1 / 2 from q ) 0, the plane of symmetry. θ-Isosurfaces are shown in intervals of 30° decreasing from θ ) 180 °, the straight line between the donor and acceptor cages. The distance between the cage centers is 5.68 Å.

Figure 6. F ( q )/ kT for the L6L hop with no water vacancy at 250 K. The dotted portion is the symmetric image of the solid q > 0 curve.

Figure 7. A vacancy in the host lattice of the hydrate was created by thermodynamic integration. The interaction parameters of one TIP4P water molecule were scaled from 1 to 0 in the hydrate lattice and from 0 to 1 in a liquid TIP4P water simulation at 250 K and 40 atm.

Figure 8. Free energy to create a water vacancy in the hydrate lattice as a function of the coupling parameter λ. The TIP4P molecule was transferred to the equilibrium 40 atm phase of TIP4P water (or ice-h) at 250 and 225 K. The small correction on the last data point at T ) 225 K is explained in the text.

A R T I C L E S Peters et al.

vacancy at 250 K in the regions listed in Table 1. The figure again includes the Jacobian factor of eq 2 for transformation to Cartesian coordinates. The L6L pathway has a barrier of just 12.6 kT , 17 kT lower than the barrier for hopping without a water vacancy. With the 4.3 kT units of free energy to create a water vacancy, the barrier for the hypothesized mechanism assisted by a water vacancy is still much lower than the barrier for the unassisted pathway. The L5L and L6L barriers were computed by using symmetry across the q ) 0 isosurface. The S5L barrier is asymmetric, but it is similar to the L5L pathway because both have a broad, shallow intermediate at the top of the barrier. The intermediate reflects configurations where the hopping methane molecule

occupies the vacancy created by a missing water molecule or a neighboring cage to the side of the donor-acceptor pathway. The intermediate configurations effectively create a doubly occupied cage^57 adjacent to both the donor and acceptor cages that is slightly stable. Figure 1 showed the configurations along the S5L pathway. Interestingly, the barriers are higher for the S5L and L5L pathways than for the L6L pathways. Transition State Theory. The rate constants for individual hops can be obtained from the computed free energy barriers using transition state theory.^58

kTST )

2 xA

〈δ[ q - q * ]〉〈| q˙ |〉 q ) q * (4a)

where

〈| q˙ |〉 q ) q * )

〈δ[ q - q * ]| q˙ |〉 〈δ[ q - q * ]〉

(4b)

and x A is the equilibrium fraction of time the system spends on the reactant side of the barrier. 〈 | q˙ |〉 q ) q * can be estimated by averaging the Jacobian scale factor^59 for the curvilinear spherical bipolar^60 reaction coordinate, hq ) ( r EQ/2)/[cosh( q ) - (θ)]:

〈| q˙ |〉 q ) q * ) 〈|∇ xDA q ( x , D , A ) · d( x , D , A ) ⁄ d t |〉 q ) q /

≈ 〈|∇ x q · x |〉 q ) q / ) 

2 kT π m

〈|∇ x q |〉 q ) q /

8 kT π mr EQ^2

(cosh( q /) - 〈cos(θ)〉 q ) q /)

Here x , D , and A are Cartesian vector positions of the methane molecule and the two ghost particles, respectively. The approximation assumes the ghost atoms, D and A , are static and that flux across q ) q * comes only from motion of the methane molecule itself, i.e. from d x /d t. r EQ values are given in Table 1. Numerical integration of cosh( q ) - cos(θ) weighted by a histogram of θ-values on the surface q ) q * gave approximately the value cosh( q ) - cos(θ) where θ is the value of θ at the saddle point in the free energy landscape. Invoking this approximation and using units of time
h ) h / kT , length λ T ) h /(2π mkT )1/2^ , and mass m ) 16 amu for methane as a united atom gives

h 〈| q˙ |〉 q ) q * ≈ 4(cosh( q * ) - cos(θ * ))λ T ⁄ r EQ (6) Table 2 summarizes the transition state theory parameters for eqs 4a and 6. Rate constants are omitted from Table 2 because the relationship between equilibrium cage occupancy and hopping kinetics reveals a missing component of the rate constants. The missing component is subtle, so we take some care to explain its origin. The probability to create a water vacancy depends mostly on the four strong hydrogen bonds that must be broken. However, the water molecule being removed also interacts with nearby methane molecules. Each water molecule has four methane neighbors in a fully occupied hydrate with full methane loading. In the situation relevant for hopping, one of the four methanes is missing, i.e., the vacant acceptor cage. Hopping via the water vacancy assisted mechanism can be

(57) van Klaveren, E. P.; Michels, J. P. J.; Schouten, J. A.; Klug, D. D.; Tse, J. S. J. Chem. Phys. 2001 , 114 , 5745. (58) Hanggi, P.; Talkner, P.; Borkovec, M. Re V_. Mod. Phys._ 1990 , 62 , 251. (59) Schenter, G. K.; Garrett, B. C.; Truhlar, D. G. J. Chem. Phys. 2003 , 119 , 5828. (60) Korn G. A. and Korn, T. M. Mathematical Handbook for Scientists and Engineers ; McGraw-Hill: New York, 1961.

Figure 9. Free energy landscapes at 250 K for the L6L (a), S5L (b), and L5L (c) hops. All plots are colored according to the same scale. q -Isosurfaces are shown decreasing and increasing in intervals of 1 / 2 from the q ) 0 plane of symmetry. θ-Isosurfaces are shown in intervals of 30° decreasing from the straight line path between donor and acceptor centers, θ ) 180 °. F ( q )/ kT is shown to the right of each two-dimensional free energy surface.

Table 2. Transition State Theory Parameters at 250 and 225 K a

250 K L6L L5L L5S S5L F D/ kT 0.1 0.2 0.1 0. F / kT 12.7 16.6 15.3 15. ∆ F / kT 12.6 16.4 15.2 14. q * 0.00 0.60 0.42 0. θ 135 ° 125 ° 120 ° 120 °
h 〈| q˙ |〉 q * 0.33 0.27 0.27 0.

225 K L6L L5L L5S S5L F D/ kT 0.1 0.2 0.4 0. F / kT 12.8 18.6 16.4 16. ∆ F / kT 12.7 18.4 16.0 16. q * 0.00 0.60 - 0.90 - 0. θ 135 ° 120 ° 115 ° 115 °
h 〈| q˙ |〉 q * 0.35 0.27 0.33 0. a (^) The subscript D denotes the free energy of the donor basin: -
F D ) ln ∫ d q exp[-
F ( q )] where integration is over q < q . q * is the value of q at the maximum in F ( q ), and θ is the minimum on q ) q * in the free energy landscapes of Figure 9.
h 〈| q˙ |〉 q * is the dimensionless absolute average velocity along q at q ) q *.

Path Sampling Study of Methane Diffusion in Hydrates A R T I C L E S

The Gillespie algorithm^18 simulates diffusion by updating a list of possible hopping “reactions” after each hopping event. When a hop occurs, the list of vacant cages and their neighbor occupancies is updated to give a new list of possible reactions. The diffusion rate is then calculated from the Einstein diffusion relation: 〈 R^2 〉 ) 6 Dt. Simulations were conducted for 8 × 8 × 8 unit cells and 1 000 000 hopping events at various concentrations. Periodic boundary conditions were used, but net displacements were monitored without periodic boundary conditions to ensure that 〈 R^2 〉 did not include finite size effects. Initial conditions were prepared with uniformly distributed vacancy patterns to ensure fast equilibra- tion with no long-lived gradients in the vacancy distribution. The systems were equilibrated for 1 000 000 hopping events before taking data. The results are shown in Figure 12. Diffusion rates are reported relative to the Einstein estimate:19,44^ D ) kL^2506 LKX (6 Å)^2 where X is the fraction of methane vacancies and 6 Å is the approximate distance for an L6L hop. The rates are nearly an order

of magnitude slower than the Einstein estimate because of the hydrate lattice topology. The fastest L6L pathways are arranged in parallel (nonintersecting) lines in the hydrate. The diffusion along L6L pathways alone is therefore a one-dimensional diffusion process. To change directions, a methane molecule must hop via one of the slower S5L or L5L pathways. At typical methane vacancy concentrations found in the hydrate, our calculations predict a diffusion constant on the order of 10 -^15 m^2 /s, a few orders of magnitude faster than the (approximate) rate of methane diffusion in an ice-h crystal^10 and a few orders of

Figure 11. Reactive flux correlation function κ( t ) for the L6L (a), S5L (b), and L5L (c) hops. Separate contributions to the reactive flux from trajectories that were initially crossing the transition state surface in the forward direction (κ+) and backward (κ-) directions are also shown. The evolving density of a swarm of trajectories initiated from the transition state surface is shown next to each correlation function.

Table 3. Summarized Transmission Coefficient Data at 250 K a

250 K L6L S5L L5L n trajectories 10000 2590 7500 t plateau 1 ps 16 ps 5 ps q * 0.00 0.45 0.
h 〈| q˙ |〉 q ) q * 0.28 0.49 0. κ( t plateau) 0.52 0.29 0. 〈 H [ q ( t plateau)]〉 q (0)) q * 0.424 0.722 0. a (^) The table also shows numerical estimates of the initial

dimensionless absolute velocity along the reaction coordinate and the fraction of trajectories from the initial swarm that commit to the product state.

Table 4. Activation Free Energy ∆ F / kT , κ,
h 〈| q˙ |〉 q * , and Rate Constants for Each Pathway a 250 K L6L L5L L5S S5L ∆ F / kT 12.6 16.4 15.2 - 16.1 13.5- 14. κ 0.52 0.24 0.29 0.
h 〈| q˙ |〉 q * 0.28 0.48 0.49 0. k /(ms-^1 ) 101.0 1.50 6.4 - 2.6 35.0- 14.

225 K L6L L5L L5S S5L ∆ F / kT 12.7 18.4 16.0 - 17.9 14.3- 16. κ 0.52* 0.24* 0.29* 0.29*
h 〈| q˙ |〉 q * 0.30 0.51 0.52 0. k /(ms-^1 ) 31.8 0.07 0.99 - 0.14 5.3- 0. a (^) The bold L5S free energy barrier is uncorrected and paired with a corrected (non-bold) S5L free energy barrier. Similarly, the bold S5L barrier is uncorrected and paired with a corrected L5S value. The corrections enforce methane’s preference for large cages over small cages as described in the text. The correction is 2 kT at 225 K and 0.9 kT at 250 K.

Path Sampling Study of Methane Diffusion in Hydrates A R T I C L E S

magnitude slower than the rate of methane diffusion through polycrystalline ice9,10^ and hydrate layers.^11 This finding supports claims that the diffusion rates in polycrystalline ice and hydrates are dominated by diffusion along grain boundaries.^9 -^11 Interest- ingly, Klapp et al.^14 demonstrated a method to determine grain size distributions in natural gas hydrates. Comparison of methane diffusion rates through hydrate samples with different grain sizes may help determine the rate of diffusion in a pure crystal. Direct measurements of the methane diffusion rate may also be possible using new single crystal growth techniques.^65 The self-diffusion constant is approximately proportional to the guest vacancy concentration over the range 0.02 < X < 0.1, which corresponds to measured hydration numbers.^66 -^70 (The hydration number “ n ” in the expression CH 4 · n H 2 O can be converted to a guest vacancy concentration X using the expression n ) 5.75/(1 - X ).) For small X , most of the vacant cages are small cages because of methane’s affinity for large cages. For X near 0.25, methane vacancies will also become common in the large cages. This would facilitate diffusion by opening the L6L pathway, but hydrates with X g 0.25 are likely to be unstable. The relationship between the diffusion constant, lattice topology, and guest molecule affinity for large cages may have important implications for larger molecules with a stronger preference for large cages. For example, CO 2 hydrates have mostly unoccupied small cages.^71 The ready availability of the S5L pathway may lead to a different relationship between diffusivity and guest occupancy in CO 2 hydrates. These effects have not been verified, but they are interesting questions for future work.

Conclusions

This paper uses an equilibrium path sampling approach^16 with reactive flux^64 and kinetic Monte Carlo^18 simulations to estimate

the methane diffusion rate in structure I gas hydrates.^1 In a hydrate I crystal at 250 K with nearly all cages occupied by methane, we estimate D ≈ 7 × 10 -^15 X m^2 /s where X is the fraction of unoccupied cages. At 225 K we estimate D ≈ 1 × 10 -^15 X m^2 /s. The calculations support a diffusion mechanism where water vacancies in the hydrate lattice provide a geometric opening for methane to pass from an occupied “donor” cage to an adjacent “acceptor” cage.^19 The free energy landscape reveals differences between the three pathways by which methane hops between different cages in the hydrate structure. The pathway between large cages that are separated by six-membered water rings (“L6L”) has the smallest activation barrier and no stable intermediates along the pathway. Pathways between large cages separated by five-membered water rings (“L5L”) have a large barrier with a shallow intermediate state at the top of the barrier. Reactive flux simulations17,64^ show that the intermediate leads to a low transmission coefficient. Structures associated with the intermediate state have methane occupying the position vacated by the water molecule and methane creating a doubly occupied cage that is off-pathway from donor to acceptor. A similar high barrier with a shallow intermediate at the top of the barrier was found for the pathway between a large cage and a small cage in the hydrate. Kinetic Monte Carlo simulations based on the computed rate constants showed that self-diffusion was nearly an order of magnitude slower than the Einstein estimate. The discrepancy was attributed to the connectivity of different pathways in the structure I hydrate lattice. Specifically, the fastest pathways for methane hopping are arranged in parallel (nonintersecting) lines in the hydrate. CH 4 must hop via one of the slower pathways to change directions. Because of the lattice connectivity and guest molecule preferences for large cages,^3 the diffusion rate is coupled nonlin- early to the concentration of guest molecule vacancies. The nonlinearity is small for methane which can fit in both small and large cages, but it may be pronounced for larger guest molecules with a stronger preference for large cages. From a computational perspective, this paper demonstrates that EPS^16 can compute free energies for a broader class of coordinates than umbrella sampling with molecular dynamics.^50 Reactive flux correlation functions17,64^ provided dynamically correct rates and additional insight into the effects of a stable intermediate near the top of the activation barrier. Finally, kinetic Monte Carlo^18 simulations enabled an estimate of a technologically important diffusion constant that has been difficult to measure. We hope these calculations will stimulate efforts to measure diffusion rates within single hydrate crystals and further simulations to investigate other possible diffusion mechanisms.

Acknowledgment. B.P. and N.Z. are grateful to Berend Smit for mentorship, generous support, encouragement, and helpful suggestions. We also thank Brian Anderson for hydrate forcefield files. B.P. and N.Z. were supported at the CECAM by the EC through the Marie Curie EXT and EST Projects MEXTCT-2003-023311 and MEST- CT-2005-020491, respectively. We thank Enitechnologie for support in addition to the Singapore-MIT Alliance for funding.

JA802014M

(65) Kirchner, M. T.; Boese, R;.; Billups, W. E. ; Norman, L. R J. Am. Chem. Soc. 2004 , 126 , 9407. (66) Uchida, T.; Hirano, T.; Ebinuma, T.; Narita, H.; Gohara, K.; Mae, S.; Matsumoto, R AIChE J. 1999 , 45 , 2641. (67) Seo, Y.; Lee, H.; Ryu, B. H. Geophys. Res. Lett. 2002 , 29 , 1244. (68) Ripmeester, J. A.; Ratcliffe, C. I. J. Phys. Chem. 1988 , 92 , 337. (69) Gallaway, T. J.; Ruska, W.; Chappelear, P. S.; Kobayashi, R. Ind. Eng. Chem. Fundam. 1970 , 9 , 237. (70) De Roo, J. L.; Peters, C. J.; Lichtenthaler, R. N.; Diepen, G. A. M. AIChE J. 1983 , 29 , 651. (71) Udachin, K. A.; Ratcliffe, C. I.; Ripmeester, J. A. J. Phys. Chem. B 2001 , 105 , 4200.

Figure 12. CH 4 self-diffusion rates as a function of X , the fraction of hydrate cages that do not contain a methane molecule. The solid gray regions show the range of values that can result depending on how the correction η is applied. The normalization factor k L6L(6 Å)^2 is 3.6 × 10 -^14 m^2 /s.

A R T I C L E S Peters et al.