Methane Clathrate Hydrate Nucleation Mechanism by Advanced

Article pubs.acs.org/JPCC

Methane Clathrate Hydrate Nucleation Mechanism by Advanced Molecular Simulations Marco Lauricella,† Simone Meloni,*,†,‡ Niall J. English,§ Baron Peters,∥ and Giovanni Ciccotti†,⊥ †

School of Physics, University College Dublin, Room 302 UCD-EMSC, Belfield, Dublin 4, Ireland Laboratory of Computational Chemistry and Biochemistry, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland § SEC Strategic Research Cluster and Centre for Synthesis and Chemical Biology, School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland ∥ Department of Chemical Engineering, University of California, Santa Barbara, California 93106, United States ⊥ Dipartimento di Fisica and CNISM, Università La Sapienza, P. le A. Moro 5, 00185 Rome, Italy ‡

S Supporting Information *

ABSTRACT: The nucleation mechanisms of methane hydrates are studied using welltempered metadynamics and restrained molecular dynamics. The collective variables we used to follow the process are the methane−methane and methane−water coordination numbers, from which we computed the corresponding Landau free energy surface. This surface is characterized by two minima, corresponding to the two-phase methane bubble/ water solution and clathrate crystal, and a transition state. The clathrate crystal is of type II, while in the simulation conditions (T = 273 K and P = 500 atm) the most stable phase should be type I. We constructed the steepest ascent/descent path connecting the twophase methane bubble/water solution to the clathrate state and passing through the transition state. We interpret this path as the nucleation path, which shows four phases. First, the concentration of solvated methane increases in the aqueous domain via diffusion through the methane−water interface. Second, units of methane molecules solvated in water meet to form an unstructured cluster. Third, the water content of the nucleus decreases to a value compatible with the type II methane clathrate hydrate composition. Finally, a reordering process of solvated methane and water molecules occurs in a manner consistent with the “blob” hypothesis (Jacobson, L. C.; Hujo, W.; Molinero, V. J. Am. Chem. Soc. 2010, 132, 11806−11811).

1. INTRODUCTION Clathrate hydrates are (nonstoichiometric) crystalline inclusion compounds in which a water host lattice encages small guest atoms or molecules in its cavities. Their thermodynamic stability is driven by low temperatures, high pressures, weak interactions between the guest molecules and the host lattice, and relatively strong hydrogen bonds that stabilize the host lattice of water molecules.1,2 Vast reserves of natural gas hydrates exist in the permafrost and deep ocean regions. On a volumetric basis, their methane content is comparable to that of compressed natural gas. Thus, hydrates represent a possible future energy resource3−6 and also a potentially safe and economical medium for natural gas transport.6−9 An improved understanding of hydrate nucleation and growth kinetics would facilitate many technological applications, e.g., the design of inhibitors to prevent hydrate formation in pipelines1,2,10 or promoters to drive hydrate formation in gas transport applications. Several hydrate structures have been created in the laboratory. At most natural gas compositions, temperatures, and pressures, the type I structure (sI) hydrate is stable.6,11,12 In sI clathrate, water molecules form (small) dodecahedral cages composed of 12 pentagonal edge-sharing rings (“512” cages) © 2014 American Chemical Society

and (slightly larger) tetrakaidecahedral cages comprised of hexagonal and pentagonal rings (“51262” cages). The crystal structure contains these cages in a 2/6 ratio. An image of sI clathrate is shown in Figure 1. In sI, around 85−95% of the cages are occupied, giving a ratio between methane and water molecules slightly lower than the stoichiometric 8/46. Another relevant structure of clathrate hydrates is the type II structure (sII; see Figure 1). In sII clathrate, water molecules form 512 and 51264 (hexakaidecahedral) cages. sII clathrate is typically formed by large guest molecules, larger than methane (see, for example, ref 13). The guest/water ratio, 24/136, is slightly higher than that of sI. Three mechanisms have been proposed for homogeneous methane hydrate nucleation, but lingering questions surround each of them. According to the “labile cluster mechanism”,14 enclathrated guest molecules diffuse in the liquid phase, encounter one another, and assemble stable clusters that grow into the hydrate crystal. However, Guo et al. later showed that enclathrated methane molecules in solution are exceedReceived: May 28, 2014 Revised: September 11, 2014 Published: September 11, 2014 22847

dx.doi.org/10.1021/jp5052479 | J. Phys. Chem. C 2014, 118, 22847−22857

The Journal of Physical Chemistry C

Article

Figure 1. Type I and II clathrate hydrate crystals. Oxygen atoms of water molecules are at the vertexes of 512 (green), 51262 (orange), and 51264 (red) polyhedra. The crystallographic positions of the carbon atoms of methane molecules are represented by spheres, which are colored according to the hosting cage. The crystallographic unit cell is denoted by blue edges. Cages in the unit cell are incomplete; thus, we represented the water molecules of periodic images needed to complete them.

the growth front in which methane is highly concentrated.29 The model shows that homogeneous nucleation in this boundary layer can actually become quite fast. Understanding the mechanism of nucleation (both heterogeneous and homogeneous) can significantly benefit from the use of special techniques, such as metadynamics,30−32 temperature-accelerated molecular dynamics33/Monte Carlo34 simulations. In this work we use a combination of two rare event techniques: metadynamics and restrained molecular dynamics (RMD). We will first reconstruct the Landau free energy landscape of selected collective variables (CVs) and identify the most likely nucleation path. These CVs might be insufficient to fully characterize the nucleation mechanism. Thus, we analyze in detail the characteristics of the nucleus along the nucleation path according to an extended set of observables that might be important for nucleation. For this we use RMD, which enables us to collect copious data in the critical bottleneck region without “waiting” for thousands of spontaneous crystallization events. We “measure” the size of the nucleus and the degree of order in the nascent host lattice and around each guest molecule as the nucleus forms. In particular, the nucleus size is typically considered a key ingredient in nucleation. We show that this observable is well mapped by our CVs. Thus, even though the nucleus size is not explicitly included in our CVs, it is implicitly taken into account. The paper is organized as follows: In section 2, we describe the computational methods used in this work. In section 3, we present and discuss the results of our simulations. Finally, in section 4, we draw conclusions.

ingly rare and survive for only picoseconds before breaking apart.15,16 Radhakrishnan and Trout proposed the “local structuring” hypothesis,17 wherein well-solvated guest molecules meet, forming a cluster, with the guests in a configuration consistent with the clathrate crystal, while water molecules are still disordered. If the cluster contains a number of guest molecules larger than a given threshold value (critical nucleus), the water molecules can rearrange to form a proper hydrate framework, leading to hydrate crystallization. This mechanism was suggested in the context of CO2 hydrate formation. Jacobson et al. proposed an alternative mechanism that bears some resemblance to the local structuring hypothesis. According to their two-step “blob hypothesis”,18 the first step is formation of a disordered and highly concentrated cluster of solvated methane molecules, or “blob”. In the second step, if the blob is sufficiently large, the locally concentrated guests and water molecules can assemble in an ordered hydrate nucleus during the relatively long time required for the blob to dissolve by diffusion.19 Several molecular dynamics (MD) studies on homogeneous and heterogeneous crystallization have found that nucleation initially yields a product that is not fully consistent with the common bulk crystal structures and compositions. Depending on the simulation conditions and force fields, the nascent nucleation products have included amorphous structures,20−22 sII clathrate,23 mixed sI and sII clathrate,24,25 and structures assembled from atypical cages.25,26 Nucleation in each of the aforementioned studies was extremely fast because the simulations were performed at very high driving forces. Knott et al. showed that at natural pressure and temperature conditions the nucleation rate is essentially zero,27 thus confirming the predictions by Vatamanu and Kusalik that nucleation in natural environments must be heterogeneous.28 A recent model shows that, for the low temperatures where hydrates are stable, nonequilibrium freeze concentration driven by ice growth can induce hydrate nucleation. Specifically, growing ice rejects methane, thus creating a boundary layer at

2. THEORY In metadynamics we introduce a set of CVs, {θi(r)}i=1,m, where r is the 3N-dimensional vector of the atomic positions. The atoms evolve driven by the physical potential, U(r), plus a history-dependent bias potential, V({θi(r)}i=1,m,t). This historydependent potential is constructed by summing up Gaussian functions deposited along the trajectory: 22848



Article

lt

V ({θi(r)}i = 1, m , t ) =

methane−methane coordination number of the ith CH4 molecule. NMM is the average methane−methane coordination number up to a cutoff distance rMM cut = 5.5 Å, corresponding to the first minimum of the pair correlation function in fluid methane. CMM is obtained by dividing NMM by the methane− methane coordination number in the methane fluid phase. With this choice, CMM takes a value close to 1 in the biphasic water/ liquid methane state and a value near 0 in the all-clathrate phase where all methane−methane distances exceed rMM cut . We use a common smooth approximation to the coordination number, ni,M = ∑j≠iΘ(rMM cut − |rij|), where Θ(x) = exp(λx)/[exp(λx) + 1], and λ is a parameter controlling the smoothness of the approximation.38,39 The purpose of CMM is to controllably drive methane into the aqueous solution, thus allowing exploration of the effect of methane concentration fluctuations on the nucleation mechanism. In fact, as mentioned above, in a recent work29 it was proposed that driven or spontaneous fluctuations of supersaturation (methane concentration higher than the equilibrium value) can enhance nucleation. Indeed, increased methane concentration in solution is observed in correspondence with nucleation events in brute force MD simulations.22,25,40 The second CV, CMW, measures the coordination number of methane to water molecules, relative to the coordination number of methane to water in a perfect clathrate. The purpose of this second CV is to enforce the proper degree of solvation of methane, thus inducing nucleation. The definition of this CV is analogous to the previous one, with the i index of the sum now running over the nW water molecules. The parameter rMW cut is again set to 5.5 Å, corresponding to the position of the first minimum in the methane−water pair correlation function in an sI all-clathrate system at our simulation conditions. CMW takes value near 1 in a pure clathrate system and close to zero in a biphasic water/fluid methane sample (CMW ≠ 0 in this case because methane molecules at the CH4/H2O interface have a certain degree of water coordination). While it is evident that the collective variables we used are able to distinguish between the initial (two-phase system) and final (all clathrate) states, it is not clear that they can distinguish between the two possible clathrate phases, sI and sII. To show that this is the case, in Figure 2 we report the probability density function to observe values zMM and zMW of the collective variables CMM(r) and CMW(r), p(zMM,zMW), for the sI and sII methane clathrate hydrate phases. sI clathrate is characterized by a higher value of CMM and a lower value of CMW with respect to sII. The distributions p(zMM,zMW) in the two phases do not overlap, proving that CMM and CMW are a good set of CVs for the problem at hand. Since we can represent all the relevant states of the system with the selected CVs, our simulations are not biased toward the formation of a specific phase of clathrate hydrate. In the Supporting Information we further discuss the choice of cutoff values used in the definition of CMM(r) and CMW(r) and their effect on the modeling of the clathrate state. For the modeling of molecular interactions, we use the Jacobson and Molinero potential,41,42 in which methane and water molecules are represented as point particles. The water− water potential is of the Stillinger−Weber type, which favors a tetrahedral coordination consistent with that of water. Methane−methane and methane−water interactions are described by the two-body term of a Stillinger−Weber-type potential. The advantage offered by this potential is two-fold. First, it contains only short-range terms, which makes it

∑ wk k=1 m

⎡ |θ (r) − θ (r(k Δt ))|2 ⎤ i i ⎥ 2σ 2 ⎣ ⎦

∏ exp⎢− i=1

where the index k runs over the lt Gaussians “deposited” up to time t, and Δt is the deposition time interval. Here we use the well-tempered version of metadynamics,32 in which the weight of the Gaussian biasing terms, wk, is also history-dependent: wk = w0 exp[−V({θ(rk−1)}i=1,m,tk−1)/Vref)]. Vref is a parameter controlling the rate at which the weight decreases with the value of the biasing potential, and w0 is the initial (and highest) value of the weight. How we obtained a suitable value for Vref is explained below. Within well-tempered metadynamics, the free energy is computed via the relation F({Z i } i=1,m ) = limt→∞[−[(kBT + Vref)/Vref]V({Zi}i=1,m,t)] (see ref 32), where kBT is the thermal energy and zi is a realization of the collective variable θi. Well-tempered metadynamics offers several advantages over the standard method. First, the parameter Vref limits the exploration of z-space (essentially) to the region with free energy on the order of Vref + kBT. This implies that the dynamics continues to visit the relevant region of the z-space, depositing Gaussians with a smaller weight as the simulation goes [V({θi(r)}i=1,m,t) grows with t]. This results in a smoother and more accurate reconstruction of the free energy landscape. To ensure that the system visits the relevant part of the space, Vref must be set to a value higher than the relevant barrier of the process under investigation. In rare events this is typically much higher than kBT, which implies that the term (kBT + Vref)/Vref in the relation above is ∼1 and the free energy is equal to the biasing potential like in standard Metadynamics. A problem with well-tempered metadynamics is that the value of the relevant barrier to which Vref is set is, usually, not known in advance. Here, we first ran a short standard metadynamics simulation with a large value of w and then set Vref to twice the value of the barrier estimated in this simulation. RMD is used to compute conditional ensemble averages of selected observables at given values {zi}i=1,m of the CVs considered. In RMD the system is driven by the biased potential U(r) + ∑im= 1(κ/2)|θi(r) − zi|2. In the limit of large κ, this dynamics samples the conditional probability density P(r| {θi(r) = zi}i,m) = exp[−βU(r)]∏mi=1δ(θi(r) − zi)/8 z , with β = 1/k B T the inverse temperature and 8 z = ∫ dr exp[−βU(r)]∏mi=1δ(θi(r) − zi). Conditional ensemble averages can, therefore, be computed from configurations taken along RMD trajectories. A possible CV to study nucleation is the size of the growing cluster, measured, for example, by the number of molecules forming it.35 However, already for simple (Lennard-Jones), single-component liquids the identification of particles belonging to the nucleus is not straightforward.36,37 We use a different set of CVs, which are able to control/induce nucleation. Then we interpret our results in terms of the nucleus size (see section 3), which allows us to make a comparison with previous studies based on the classical nucleation theory. In this work we use two CVs. The first, CMM, is the mean methane−methane relative coordination number. CMM is n computed from NMM = (1/nM)∑i = M 1ni,M(r), with nM the number of methane molecules in the sample and ni,M(r) the 22849



Article

3. RESULTS AND DISCUSSION The CMM−CMW Landau free energy landscape (F(CMM,CMW)) obtained from well-tempered metadynamics is shown in Figure 3. F(CMM,CMW) presents a global minimum at low CMM/high

Figure 2. Color map of the probability density function to observe values zMM and zMW of the collective variables CMM(r) and CMW(r), p(zMM,zMW), in the sI (top left) and sII (bottom right) methane clathrate phases. p(zMM,zMW) shown in the figure is obtained by the union of histograms of CMM(r) and CMW(r) obtained from MD simulations in the sI and sII states. Red (black) means high (low) probability density. Data in the figure show that the distributions in the two states do not overlap, which means that the two collective variables are able to distinguish between systems in the sI and sII states. This implies that there is no bias in our metadynamics simulation toward either of the clathrate phases.

Figure 3. Landau free energy landscape of the collective variables CMM(r) and CMW(r). Free energy is reported in kBT units. The green line represents the NP, and the five green points denote the CMM− CMW states at which we analyze in detail the structure of the sample. The yellow lines around the NP represent the external surface of the reaction tube (see the text), which is rather narrow all along the nucleation path.

computationally very efficient. Second, since it is based on a united atoms representation of CH4 and H2O, it allows a longer time step, 5 fs against a typical value of 1 fs for force fields with explicit hydrogen. This potential has already been validated and successfully used in the investigation of some aspects of nucleation of clathrates of several guest species.43 The (initial) simulation box is a cubic cell with sides of 48 Å. This box is filled, randomly, with 512 and 2944 molecules of methane and water, respectively. The composition corresponds to a methane molar fraction of 0.14, which is the stoichiometric composition of the sI methane hydrate structure.2 The system is equilibrated with an NPT MD simulation44 at 250 K and 500 atm (identical to that of Walsh et al.25) for 25 ns. We remark that these conditions amount to ∼55−60 K of supercooling with respect to the melting temperature of sI (307 K) and sII (303 K) methane clathrates obtained with a similar potential (same functional form, slightly different parameters).18 At the end of the equilibration process, we observe a single (approximately spherical) “drop” of liquid methane in water. We, then, perform a further 1 ns NPT MD simulation to measure the mole fraction of solvated methane, which is found to be 0.016. To determine a suitable value of Vref, we perform a 30 ns “standard” (i.e., non-well-tempered), coarse (w = 3 kBT, σ = 0.02, Gaussian deposition frequency 1 ps) NVT metadynamics simulation at 250 K on a sample at the density corresponding to the average density measured in the NPT simulation. The value of the largest free energy barrier, corresponding to the melting of clathrate, is estimated to be approximately 500 kBT. We set Vref to twice this value. Well-tempered metadynamics simulation is then performed at the same temperature and density with σ = 0.01 and a Gaussian deposition frequency of 1 ps. This simulation was run for 2 μs, during which the biasing potential reaches a value of 2400 kBT and wk ≤ 0.3 kBT.

CMW corresponding to the clathrate state. At high CMM/low CMW, the substantial local minimum corresponds to a state consisting of a drop of fluid methane in a methane/water solution; we refer to this as the liquid state. A visual inspection has shown that the structure of the clathrate minimum is of type II (Figure 4C). To confirm this observationm we identify those 512 cages surrounded by the correct number of 51264 cages in sII clathrate and vice versa. This analysis shows that 80% of enclathrated CH4 molecules are contained in sII-like cages. The rest are contained in sI domains, in cages of a different type, and in the water/methane solution. The fact that we observe the formation of sII methane clathrate in thermodynamic conditions in which the stable phase is sI might appear surprising. However, Molinero and coworkers also observed the formation of sII clathrate in their simulations with the same potential used in the present work.20 We believe that the formation of the sII structure is kinetically favored over that of the sI structure; thus, the former is more likely to form in (finite duration) simulations if no CVs specifically accelerating the formation of the latter are used. This hypothesis is further discussed below. We compute the steepest ascent/descent path connecting the local and global minima of the free energy surface (green line in Figure 3) using the string method.45 Here the potential energy used in the original method is replaced by the Landau free energy computed by metadynamics. We interpret this path as the nucleation path (NP). Along the NP we compute the clathrate nucleation barrier, i.e., the free energy barrier that the solution must overcome to nucleate a clathrate crystal. This barrier amounts to 130 kBT. Experiments are usually performed 22850



Article

Figure 4. Snapshots of the system at points 3 (A), 4 (B), and 5 (C) of the NP. Cyan spheres represent CH4 molecules belonging to the methane bubble and blue ones those belonging to the nucleus. Red spheres represent water molecules belonging to the H2O/CH4 solution (i.e., not assigned to the cluster), orange ones those belonging to the cluster but not forming any cage. Red polyhedra denote the structure II part of the nucleus, and yellow polyhedra denote 51262 cages, i.e., the structure I part of the nucleus. The central part of the ordered region of the TS nucleus is of the sII type, with the external shell made by sI. When the nucleus grows, the sII region entraps the sI part.

at constant pressure; thus, we computed a correction, described in the Supporting Information, to estimate the constantpressure free energy barrier from the one we obtained from NVT metadynamics. At 500 atm and 250 K, the NPT barrier is 173 kBT. We note that the NP consists of three branches, indicated by three arrows in Figure 3. In the first branch, CMW increases and CMM decreases. This corresponds to the shrinking of the methane bubble (decrease of CMM), with an increase of the concentration of methane in the methane/water solution (increase of CMW). In the second branch, which contains the transition state (TS), CMM remains relatively constant while CMW further increases. This corresponds to the formation of a sizable water-rich CH4/H2O clathrate nucleus (see below). Finally, in the third branch, CMW further increases and CMM decreases. In this postcritical branch the methane molecules evaporating from the CH4 liquid bubble are absorbed by the clathrate nucleus. The NP shown in Figure 3 represents the most likely but not unique nucleation path. At finite T, nucleation trajectories do not coincide with the NP. Thus, it is interesting to investigate the width of the “reactive tube”,46 i.e., the tube around the NP containing most of the nucleation trajectories (the theoretical aspects beneath the calculation of the reactive tube are explained in the Supporting Information). In our 2D CV space, we identify the reactive tube with the union of segments orthogonal to the NP having a Landau free energy within kBT from the point at which they intersect the NP. The external surface of this tube is represented by the two yellow lines on the two sides of the NP in Figure 3. The reactive tube is very narrow in the initial and final parts of the ascending and descending branches of the NP, and slightly wider in the TS region. In this region the surface of the reactive tube is also less smooth. This is due to the roughness of the free energy surface, on one hand, and the limited numerical accuracy in the determination of the straight lines orthogonal to the path, on the other hand. The limited smoothness of the tube, however, does not change the qualitative and general conclusion that the paths followed by the system during nucleation are expected to pass very close to the NP. To better understand the mechanism of nucleation, we compute the value of a series of observables along the NP. This

requires the calculation of averages over the conditional probability density ρ(r|CMW=zMW,CMM=zMM), representing the probability density to be at the configuration r and satisfying the conditions CMW(r) = zMW and CMM(r) = zMM, with zMW and zMM two possible values of the relative methane−water and methane−methane coordination numbers. This is obtained by RMD.33,38,39 We considered a set of five {zMW,i,zMM,i}i=1,5 points along the NP, corresponding to the liquid (i = 1) and clathrate (i = 5) minima, the TS (i = 3), and the intermediate points between the two minima and the TS (i = 2 and i = 4; see Figure 3). Since our system might contain further metastabilities, the sampling of ρ(r|CMW=zMW,i,CMM=zMM,i) by RMD can be biased by the choice of the initial configuration.39 To check that this is not the case, at each of the five (zMW,i, zMM,i) points, we run three RMD simulations. These simulations start from configurations taken from three reactant-to-product/productto-reactant branches of the metadynamic trajectory. In particular, these configurations are those with (CMW(r),CMM(r)) the closest to the point (zMW,i,zMM,i). In the presence of further metastabilities, some of these configurations should belong to a different basin (see ref 39), and observables computed in RMD simulations started from different configurations should produce inconsistent results. This phenomenon was not observed in our simulations, thus confirming that our sampling of ρ(r|CMW=zMW,i,CMM=zMM,i) is unbiased. In our analysis we consider six observables, defined and described in detail later: (i) the number of methane molecules in the largest clathrate nucleus, nnM; (ii) the methane molar fraction in the same nucleus, χnM; (iii) a measure of the structure of the water network in the nucleus, F4;47 (iv) the number of water cages per methane molecule in the nucleus, nc, and of the most relevant types of cages, nxc (with x = 512, 51262, etc.; see below); (v) the methane−methane pair correlation function, gMM(r); (vi) the methane molar fraction in solution, χsm. To count the number of methane molecules in the largest clathrate nucleus, we must first give a precise definition of clathrate nuclei. We consider well-hydrated methane molecules, denoted by the symbol CH4h in the following. CH4h molecules are those methane molecules coordinated with 18 or more water molecules. The threshold to define a methane molecule to be well hydrated is not a critical parameter because the level 22851



Article

instead of the supposedly more stable sI: if the formation of sII is much faster than that of sI in a finite simulation one could miss exploring the basin of the more stable phase. It must be remarked that the comparison between the present and previous results of Knott et al. is not straightforward. First, the thermodynamic conditions are different, 300 K and 900 atm in ref 27 vs 273 K and 500 atm in the present work. Second, the free energy computed in the two cases refers to different CVs. In principle, one could be computed from the other using the co-area formula,49 but this requires the calculation of derivatives between the two sets of CVs, which is not readily available in our case. Finally, in ref 27 nuclei are crystalline, while in the present work we show that at the TS they are still amorphous (see below). However, we notice that the cluster at the fourth point of the NP, which is well ordered and supercritical, is still sizably smaller (∼230 vs 300) than the crystalline nucleus of Knott et al. Finally, in ref 27 the methane molar fraction in solution is at the equilibrium value at 300 K and 900 atm. In the present work it changes along the NP (see below). We now focus on the change of the methane/water composition of the nucleus along the nucleation process. In particular, we measure the methane molar fraction, χnM, defined as the number of CH4 molecules over the total number of molecules, water plus methane, in the nucleus. This requires the counting of water molecules in the nucleus, which we define as those molecules coordinated to the CH4h assigned to the nucleus in the previous step. In Figure 6 we report the

of hydration of these molecules in solution, or in a clathrate, is much higher than that at the surface of the CH4 bubble (methane molecules inside the bubble have a water coordination equal to zero). We, then, search for clusters of CH4h molecules, defined as sets of molecules at a distance lower than 9 Å from (at least) another molecule of the set.34 For characterizing the nucleation process, we focus on the largest cluster, i.e., the one containing more CH4h molecules. The conditional distribution of nnM at the five points along the N P , p ( n Mn ; i ) ( i = 1 , . . . , 5 i s s h o r t h a n d f o r (CMM=zMM,i,CMW=zMW,i)), is reported in Figure 5. As a first

Figure 5. Distribution of the number of methane molecules in the nucleus at the five points along the NP described in the text, p(nM n ;i) (see also Figure 3). The color coding introduced here, blue for the first point, cyan for the second, orange for the third, red for the fourth, and violet for the fifth, is used throughout the paper. The peaks are narrow and well separated, and their positions grow along the path. These facts suggest that there is a map between CMM−CMW and nnM values.

remark, we notice that all the distributions are monomodal and narrow and, apart from the first two points, well separated. The width of these distributions probably could be further reduced by using an even higher value of the restraint constant κ in the RMD simulations. These facts suggest that there is a map between points along the NP in the CMM−CMW CV space and the size of the clathrate nucleus. Thus, even though we did not include the nucleus size among our CVs, this key observable of nucleation theories is, indeed, implicitly taken into account in our simulations. We compare our results with those of a previous computational work.27 We focus on this work because it is based on a potential similar to the one used in the present work. We compare the size of the largest cluster at the TS with results reported by Knott et al.27 Knott et al. computed the size of the critical nucleus for the water/methane-to-type I clathrate crystallization by following the fate of initially perfectly crystalline and spherical nuclei immersed in a water/methane solution at the target methane molar fraction. If a nucleus grows, it is supercritical if it shrinks, it is subcritical, and if it does not change in size, it is critical. From the size of the critical nucleus, measured to be ∼300 cages at 273 K and 900 atm,48 and from the estimation of other parameters, Knott et al. computed the value of the nucleation barrier to be 300 kBT. The size of the nucleus at the TS in the present work (250 K and 500 atm) is ∼65 cages and ∼110 CH4 molecules (not all hosted in well-formed water cages) and the barrier ∼173 kBT. Both quantities are sizably smaller than the corresponding values of ref 27. This fact might explain why we obtain sII

Figure 6. Distribution of the methane molar fraction in the nucleus along the NP. At the beginning of nucleation, the distribution is wide, indicating large fluctuations of the water/methane composition of the nucleus. However, already at the transition state the distribution is quite sharp. A narrow distribution indicates that adding/removing water from a nucleus is energetically expensive, which supports our hypothesis of why formation of sII clathrate is kinetically favored.

probability distribution of χnM along the NP. We notice that the fraction of methane increases with the size of the nucleus, and it is initially much lower than that in an ideal methane clathrate. This can be explained as follows. The methane/water solutions contain CH4h(H2O)n units, with n ≈ 18, to be compared with 46/8 ≈ 5.75 and 136/24 ≈ 5.6 in structure I and II clathrates, respectively. When these units collide with, and stick to, a nucleus, they produce a new, larger, nucleus with a lower methane molar fraction. To form a proper clathrate structure (I or II) the extra water must be removed. Our results indicate that the initially amorphous cluster50 evolves into an ordered clathrate-like nucleus by expelling part of the water and 22852



Article

superimposed, indicating that the ordering of water in the clathrate is complete. Our results confirm previous findings that the clathrate nucleus is initially disordered.20−22 However, all these studies, including the present one, have been performed at high supercooling. An interesting point we plan to address in the future is how the degree of supercooling affects the degree of disorder in subcritical and critical nuclei: is the nucleus more or less amorphous in conditions of lower supercooling? Another important aspect concerning the structure of water molecules in the nucleus is their arrangement in cages. As we mentioned in the Introduction, in the sI clathrate CH4 molecules are hosted in 512 and 51262 H2O cages, while in sII they are hosted in 512 and 51264 cages. It has been proposed25,53 that along the nucleation other cages can be present as well, namely, 51263, 4151062, 4151063, and 4151064. Cages can be identified following an algorithm suggested by Molinero and co-workers.54 In Figure 8 we report the distribution of the total

arranging the remaining water in the proper geometry. This two-step mechanism is consistent with previous results.22 We speculate that the change of the water content in the nucleus along the crystallization explains the proposed faster kinetics of formation of sII over sI. The water content of the nucleus is initially higher than that in sI and sII clathrates. During the process the water content in the nucleus decreases, reaching first the composition of sII, which contains more water. If the energetic cost to further reduce the water content is higher than that of transforming it into an sII nucleus, the second process will prevail. To analyze how the framework formed by water molecules in the nucleus changes along the NP, we measure the F4 order parameter47 of H2O molecules belonging to the nucleus at the five points of Figure 3. F4(r) is related to the Oi−Oj−Ok−Ol dihedral angles of four-term chains of hydrogen-bonded H2O molecules: F4(r) = (1/nd)∑i∑j≠i∑k≠i,j∑l≠i,j,k cos(3θi,j,k,l), nd denoting the number of dihedral angles in the nucleus and θi,j,k,l the value of the dihedral angle formed by the molecules in the subscript. The O−O−O−O dihedral angles in liquid water take all possible values, with a ∼25% difference between the most (60°, 180°, and 300°) and least (0, 120°, and 240°) probable angles.51 Thus, F4(r) in bulk water is zero. In clathrates, where water molecules form (essentially) planar five- and six-term rings, most of the dihedral angles take values close to 0°. However, the nucleus also contains O−O−O−O chains formed by water molecules belonging to different faces of a cage, or different cages, whose dihedral angles take values different from 0° (the value depends on the types of faces and cages). For bulk methane clathrate we observe a value of F4 ≈ 0.8, slightly different from values reported in the literature. This is because we use a slightly different definition of F4 (see ref 52). F4 allows liquid-like, disordered (low F4) and clathrate-like, ordered (high F4) water to be distinguished. In Figure 7 we report the F4

Figure 8. Distribution of the total number of water cages per methane molecule along the NP. As in the case of F4, the position of the peak of the nc distribution also changes suddenly along the NP, further supporting the hypothesis that water ordering is an induced event.

number of cages per methane molecule along the NP, p(nc,i), while in Figure 9 the distribution of cages of types 512, 51262, and 51264, p(nxc ,i) (x = 512, 51262, 51264), is shown. In the inset of Figure 9 we also report the ratio between the numbers of cages of types 51262 (51264) and 512 along the NP. The numbers of cages of other types are small at all points along the NP, and their trend is not discussed. nc is lower than 1 all along the NP, indicating that not all the methane molecules in the nucleus are hosted in well-formed cages. In the early stage of nucleation (first and second points along the NP), water molecules form no cages. At the transition state, ∼60% of the methane molecules are hosted in cages. Finally, in the second half of the NP, ∼90% of the CH4 molecules are in cages. As for F4, this value does not change significantly between the fourth and fifth points. There are two principal reasons why not all water molecules form cages. First, as discussed above, the crystal structure of the clathrate minimum is sII, which has a methane/water composition inconsistent with that of our sample (our sample, with a 0.14 CH4 mole fraction, is consistent with the sI structure). Second, even a small misalignment of the “lattice” directions of the growing nucleus with respect to the edges of the simulation box prevents the formation of a perfect crystal, which results in a

Figure 7. F4 at the five points along the NP. The distributions are initially very broad, becoming sharper along the path. The sudden change of the position of the F4 peaks between points 2 and 3 and points 3 and 4 indicates that the ordering of water occurs suddenly, which supports the assumption that it is triggered by another event: the ordering of CH4 molecules.

distribution at the five points along the NP. We notice that the distribution is initially very broad and centered around zero. This indicates that water in the cluster is liquid-like. At the TS, the F4 distribution becomes more peaked and the position of its maximum moves toward higher values, indicating that water has a more ordered, clathate-like configuration. Finally, at the last two points, the distribution is very sharp and essentially 22853



Article

nucleus. Then, when the nucleus grows, they remain trapped within the sII crystallite (Figure 4B,C). A question discussed in the literature is whether it is the ordering of water that induces the ordering of methane (or vice versa) or if the ordering of the two species occurs at the same time (see refs 24 and 55). To address this question, we analyze the methane−methane pair correlation function, gMM(r), among the methane molecules belonging to the cluster along the NP (Figure 10). gMM(r) at the first and second points is

Figure 9. Distribution of the most relevant water cage type along the 12 2

12

12 4

12

NP. In the inset we report n5c 6 /n5c (▲) and n5c 6 /n5c (■) at points 12 2

12

12 4

12

3−5 of the NP. We recall that n5c 6 /n5c = 1/3 in sI and n5c 6 /n5c = 1/2 in sII clathrates. The cage composition of the nucleus differs from these reference values all along the NP. This is because it is made of sI and sII regions and disordered domains.

Figure 10. Methane−methane pair correlation function at the five states along the NP. In the case of small nuclei (points 1 and 2 along the NP), due to the poorer statistics, g(r) is “noisy”. The quality rapidly improves along the NP. It is worth remarking that, from the third point on, the second peak of g(r) shows only small changes, indicating that the methane molecules in the cluster achieved their final configuration already at the TS.

disordered structure of the water molecules belonging to the external shell of the clathrate nucleus. It is also interesting to follow the distribution of the number of water cages of a given type per methane molecule, nxc , along the NP. We focus our analysis on points 3−5 of the NP because the small number (or complete absence) of cages in the first two points prevents this analysis from being performed. At the 12 2 12 12 4 12 TS, n5c 6 /n5c ≈ 0.1 and n5c 6 /n5c ≈ 0.4. Both take values below 12 2

12

12 4

12

12 2

12

quite noisy. This is because the number of methane molecules in the nucleus is small (see Figure 5). Nevertheless, the presence of a peak at r ≈ 3.2 Å is quite evident, indicating that in small nuclei there are methane molecules in direct contact. When the nucleus grows, the number of these pairs decreases. At the TS, the amount of methane molecules in direct contact with each other is already very small and mostly concentrated in the external shell of the nucleus, where the cages around CH4 are not well formed. At the same time, the clathrate-like peak at 6.8 Å becomes higher and sharper. Between points 4 and 5, we still observe some small change in the first peak, indicating a further reduction of CH4 molecules in direct contact. On the contrary, we observe no changes in the clathrate-like peak, suggesting that the ordering of CH4 molecules in the nucleus is essentially completed already at the TS. If we consider that gMM(r) does not change after the TS, while F4 and nc at the TS and at the next points are significantly different (compare Figures 7, 8, and 10), we must conclude that it is the methane ordering, which occurs first, triggering the formation of the water framework. This observation is consistent with the blob hypothesis of Jacobson et al.18 and previous simulation results on a sample with an ad hoc configuration of CO2 molecules.55 Lastly, we discuss the effect of fluctuation of the composition in the methane/water solution on clathrate nucleation. Recently, Poon and Peters29 discussed freeze concentration as a route to large supersaturations at clathrate nucleation conditions. Walsh et al.,25 in their brute force MD simulations, observed nucleation in conditions of large supersaturation, χsm = 3.9 × 10−2, to be compared with an equilibrium mole fraction χsm = 2.5 × 10−3. English et al. have observed solute enrichment,

the ideal composition of sI (n5c 6 /n5c = 1/3) and sII (n5c 6 /n5c = 1/2) structures. However, this cage composition suggests that the nucleus is mainly the sII type with, possibly, a small part of 12 4 12 sI. n5c 6 /n5c grows toward the ideal sII value at point 4 and, then, slightly decreases at point 5. On the contrary, n5c 6 /n5c steadily grows along the NP. To measure more precisely the amount of sI and sII 12 2 12 12 4 12 clathrates in the nucleus, the n5c 6 /n5c and n5c 6 /n5c ratios are not sufficient. In a clathrate crystallite, 512 and 51262 cages, for sI structures, and 512 and 51264 cages, for sII structures, have a precise spatial distribution. Thus, we performed an analysis analogous to that performed above to identify the type of structure of the global minimum of the free energy. In Figure 4 we report snapshots of the nucleus in which we have highlighted sI and sII clathrate subdomains in points 3−5 of the NP. To make this analysis more quantitative, we counted the number of CH4 molecules enclathrated in the sII subdomain over the total number of methane molecules belonging to the nucleus. We find that at the TS ∼50% of the methane molecules are enclathrated in the sII subdomain, while at points 4 and 5 this number grows to ∼80%. We repeated a similar analysis for 512 and 51262 cages, but we found that only a small number of them are in an arrangement consistent with sI clathrate. When the nucleus is small (at the TS; see Figure 4A), these cages lie in the external shell of the 22854



Article

As a first remark, we observe the formation of the sII methane clathrate, rather than the expected sI structure. We propose that the sII structure is kinetically favored over sI, and thus, in the absence of CVs explicitly enhancing nucleation of the sI structure, the structure with a lower barrier prevails in our 2 μs metadynamics simulation. The nucleation barrier is estimated to be 173 kBT, so most likely the formation of clathrate hydrates in nature occurs via heterogeneous rather than homogeneous nucleation. The nucleus is initially disordered and becomes more ordered when it approaches the transition state. This result is consistent with previous studies.18,22 We also found that the methane−water composition of the nucleus changes with its size. When the nucleus is small, its water content is well above that of the clathrate phases. The water content decreases to the composition of sII while the nucleus grows. More water need to be expelled to create an sI nucleus, potentially explaining the selection of sII over sI. Also the level of ordering of water, measured by both F4 and the number of water cages per methane molecule, grows with the size of the nucleus. However, a well-ordered water framework is obtained only toward the end of the nucleation, at points 4 and 5 along the NP in our simulations. On the contrary, methane molecules are in a good clathrate-like configuration already at the TS, indicating that it is the ordering of the guest species that induces the ordering of the water framework in the amorphous nucleus.

anticipating the formation of precursors of nucleation in their nonequilibrium simulations of a planar water/methane interface.56 Finally, Vatamanu and Kusalik22 have shown that local fluctuations of solute concentration beyond the equilibrium value might enhance nucleation. Analogous effects have also been observed in other clathrates (see, for example, ref 40). Thus, it is interesting to follow the methane molar fraction in solution along the NP (Figure 11). χsM is computed as the ratio

Figure 11. Distribution of the methane molar fraction in solution at the five points along the NP. The distributions at points 4 and 5 are affected by an algorithmic artifact explained in the text. It is worth remarking that the methane content in the solution at the second point along the NP is significantly higher than the equilibrium value. This suggests that there are methane molar fraction fluctuations triggering clathrate nucleation.

■

ASSOCIATED CONTENT

S Supporting Information *

Validation of the collective variables, estimation of the NPT barrier from NVT metadynamics simulations, and description of the reactive tube. This material is available free of charge via the Internet at http://pubs.acs.org/.

between the numbers of methane and water molecules assigned neither to the nucleus nor to the methane bubble. It is worth recalling that a water molecule is assigned to the nucleus if it coordinates a CH4 molecule already assigned to the nucleus (see the definition of χnM). Thus, all water molecules at the boundary between the nucleus and aqueous solution are assigned to the nucleus. At the fourth and fifth points of the NP, when the nucleus almost completely occupies the simulation box, this protocol results in an artificially high χsM. We, thus, do not consider these points in our analysis. On the contrary, at the first three points, i.e., up to the TS (included), the nucleus is still quite small and the arbitrary assignment of nucleus/solution boundary water molecules to the former does not critically affect the estimate of χsM. Figure 11 shows that along the nucleation χsM first increases and then decreases. This suggests, consistently with refs 22, 25, 40, and 56, that a temporary enrichment of the solute concentration in methane/water solution beyond the equilibrium value is the event that triggers nucleation.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: simone.meloni@epfl.ch. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Valeria Molinero for providing her code for cage identification and interesting discussions. G.C., S.M., and M.L. acknowledge financial support from SFI Grant No. 08-IN.1I1869. G.C. and S.M. acknowledge financial support from the Istituto Italiano di Tecnologia under SEED Project Grant No. 259 SIMBEDDAdvanced Computational Methods for Biophysics, Drug Design and Energy Research. S.M. acknowledges financial support from MIUR-FIRB Grant No. RBFR10ZUUK. B.P. gratefully acknowledges the US National Science Foundation for support through CE-1125235. We thank Ireland’s High-Performance Computing Centre (ICHEC) and the IT service of the University College Dublin for the provision of computational resources.

4. CONCLUSIONS In this work, we have investigated the mechanism of (homogeneous) nucleation of methane clathrate hydrate by means of well-tempered metadynamics and restrained molecular dynamics. Nucleation is a very complex phenomenon, and the relevant degrees of freedom controlling the process are still debated.57 Thus, here we first reconstruct the Landau free energy of two collective variables able to produce the nucleation and then characterize the nucleation path connecting the liquid−transition−clathrate states according to a series of observables measuring the size, composition, and structure of the nucleus.

■

REFERENCES

(1) Makogon, Y. Hydrates of Hydrocarbons; PennWell Books: Tulsa, OK, 1997. (2) Sloan, E.; Koh, C. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press, Taylor & Francis: Boca Raton, FL, 2007. (3) MacDonald, G. J. The Future of Methane as an Energy Resource. Annu. Rev. Energy 1990, 15, 53−83.

22855



Article

(4) Kvenvolden, K. A. Methane Hydrates: A Major Reservoir of Carbon in the Shallow Geo-Sphere? Chem. Geol. 1988, 71, 41−51. (5) Moridis, G. J.; et al. Challenges, Uncertainties, and Issues Facing Gas Production from Gas-Hydrate Deposits. SPE Reservoir Eval. Eng. 2011, 14, 76−112. (6) Koh, C.; Sloan, E.; Sum, A.; Wu, D. Fundamentals and Applications of Gas Hydrates. Annu. Rev. Chem. Biomol. Eng. 2011, 2, 237−257. (7) Chatti, I.; Delahaye, A.; Fournaison, L.; Petitet, J. Benefits and Drawbacks of Clathrate Hydrates: A Review of Their Areas of Interest. Energy Convers. Manage. 2005, 46, 1333−1343. (8) Wen, Y.; Chen, Q.; Chen, Y.; Fan, S. Research Progress on Hydrate Self-Preservation Effect Applied to Storage and Transportation of Natural Gas. Adv. Mater. Res. 2013, 772, 795−801. (9) Nakoryakov, V.; Misyura, S. The Features of Self-Preservation for Hydrate Systems with Methane. Chem. Eng. Sci. 2013, 104, 1−4. (10) Hammerschmidt, E. G. Formation of Gas Hydrates in Natural Gas Transmission Lines. Ind. Eng. Chem. 1934, 26, 851−855. (11) Lundgaard, L.; M?llerup, J. Calculation of Phase Diagrams of Gas-Hydrates. Fluid Phase Equilib. 1992, 76, 141−149. (12) Kvenvolden, K. A. Gas HydratesGeological Perspective and Global Change. Rev. Geophys. 1993, 31, 173−187. (13) Sloan, E. D. Fundamental Principles and Applications of Natural Gas Hydrates. Nature 2003, 426, 353−363. (14) Sloan, E. D.; Fleyfel, F. A Molecular Mechanism for Gas Hydrate Nucleation from Ice. AIChE J. 1991, 37, 1281−1292. (15) Guo, G.; Zhang, Y.; Zhao, Y.; Refson, K.; Shan, G. Lifetimes of Cagelike Water Clusters Immersed in Bulk Liquid Water: A Molecular Dynamics Study on Gas Hydrate Nucleation Mechanisms. J. Chem. Phys. 2004, 121, 1542−1547. (16) Guo, G.; Zhang, Y.; Li, M.; Wu, C. Can the Dodecahedral Water Cluster Naturally Form in Methane Aqueous Solutions? A Molecular Dynamics Study on the Hydrate Nucleation Mechanisms. J. Chem. Phys. 2008, 128, 194504. (17) Radhakrishnan, R.; Trout, B. L. A New Approach for Studying Nucleation Phenomena Using Molecular Simulations: Application to CO Hydrate Clathrates. J. Chem. Phys. 2002, 117, 1786. (18) Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 2010, 132, 11806−11811. (19) Peters, B. On the Coupling between Slow Diffusion Transport and Barrier Crossing in Nucleation. J. Chem. Phys. 2011, 135, 044107. (20) Jacobson, L. C.; Molinero, V. Can Amorphous Nuclei Grow Crystalline Clathrates? The Size and Crystallinity of Critical Clathrate Nuclei. J. Am. Chem. Soc. 2011, 133, 6458. (21) Jacobson, L. C.; Matsumoto, M.; Molinero, V. Order Parameters for the Multistep Crystallization of Clathrate Hydrates. J. Chem. Phys. 2011, 135, 074501. (22) Vatamanu, J.; Kusalik, P. G. Observation of Two-Step Nucleation in Methane Hydrates. Phys. Chem. Chem. Phys. 2010, 12, 15065−15072. (23) Moon, C.; Taylor, P. C.; Rodger, P. M. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003, 125, 4706− 4707. (24) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Gas Hydrate Nucleation and Cage Formation at a Water/Methane Interface. Phys. Chem. Chem. Phys. 2008, 10, 4853. (25) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326, 1095−1098. (26) Sarupria, S.; Debenedetti, P. G. Homogeneous Nucleation of Methane Hydrate in Microsecond Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2012, 3, 2942. (27) Knott, B.; Molinero, V.; Doherty, M.; Peters, B. Homogeneous Nucleation of Methane Hydrates: Unrealistic at Realistic Conditions. J. Am. Chem. Soc. 2012, 134, 19544−19547. (28) Vatamanu, J.; Kusalik, P. G. Molecular Insights into the Heterogeneous Crystal Growth of sI Methane Hydrate. J. Phys. Chem. B 2006, 110, 15896−15904.

(29) Poon, G.; Peters, B. A Stochastic Model for Nucleation in the Boundary Layer during Solvent Freeze-Concentration. Cryst. Growth Des. 2013, 13, 4642−4647. (30) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (31) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics. J. Phys. Chem. B 2005, 109, 6714−6721. (32) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (33) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168−175. (34) Ciccotti, G.; Meloni, S. Temperature Accelerated Monte Carlo (TAMC): A Method for Sampling the Free Energy Surface of NonAnalytical Collective Variables. Phys. Chem. Chem. Phys. 2011, 13, 5952−5959. (35) Beckham, G.; Peters, B. Optimizing Nucleus Size Metrics for Liquid-Solid Nucleation from Transition Paths of Near-Nanosecond Duration. J. Phys. Chem. Lett. 2011, 2, 1133−1138. (36) Lechner, W.; Dellago, C.; Bolhuis, P. Reaction Coordinates for the Crystal Nucleation of Colloidal Suspensions Extracted from the Reweighted Path Ensemble. J. Chem. Phys. 2011, 135, 154110. (37) Jungblut, S.; Singraber, A.; Dellago, C. Optimising Reaction Coordinates for Crystallization by Tuning the Crystallinity Definition. Mol. Phys. 2013, 111, 3527−3533. (38) Orlandini, S.; Meloni, S.; Colombo, L. Order-Disorder Phase Change in Embedded Si Nanoparticles. Phys. Rev. B 2011, 83, 235303. (39) Orlandini, S.; Meloni, S.; Ciccotti, G. Combining Rare Events Techniques: Phase Change in Si Nanoparticles. J. Stat. Phys. 2011, 145, 812−830. (40) Liang, S.; Kusalik, P. G. Exploring Nucleation of H2S Hydrates. Chem. Sci. 2011, 2, 1286. (41) Molinero, V.; Moore, E. B. Water Modeled as an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2008, 113, 4008. (42) Jacobson, L. C.; Molinero, V. A Methane-Water Model for Coarse-Grained Simulations of Solutions and Clathrate Hydrates. J. Phys. Chem. B 2010, 114, 7302−7311. (43) Jacobson, L. C.; Hujo, W.; Molinero, V. Nucleation Pathways of Clathrate Hydrates: Effect of Guest Size and Solubility. J. Phys. Chem. B 2010, 114, 13796−13807. (44) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177. (45) Weinan, E.; Ren, W.; Vanden-Eijnden, E. Simplified and Improved String Method for Computing the Minimum Energy Paths in Barrier-Crossing Events. J. Chem. Phys. 2007, 126, 164103. (46) Vanden-Eijnden, E. Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology; Springer: Berlin, 2006; Vol. 1, pp 453−493. (47) Rodger, P. M.; Forester, T. R.; Smith, W. Simulations of the Methane Hydrate/Methane Gas Interface Near Hydrate Forming Conditions. Fluid Phase Equilib. 1996, 116, 326−332. (48) In ref 27, since the nucleus is created crystalline, the numbers of cages and methane in the nucleus coincide, and the size of the nucleus is measured by the former property. Thus, to make it fair, the comparison is made confronting the number of complete cages at the TS. (49) Hörmander, L. The Analysis of Linear Partial Differential Operators; Springer: Berlin, New York, 2003. (50) A possible doubt is that the nucleus is formed amorphous because of the finiteness (and small) size of the sample. However, at the TS the nucleus is significantly smaller than the simulation box (see Figure 4). Thus, we believe that finite size effects do not affect the qualitative conclusions that can be drawn from our simulations. (51) Mezei, M.; Speedy, R. J. Simulation Studies of the Dihedral Angle in Water. J. Phys. Chem. 1984, 88, 3180−3182. 22856



Article

(52) In its original definition, F4 is computed from H−O−O−H dihedral angles of two hydrogen-bonded molecules. In our atomistic model water molecules are represented by point particles; thus, we have to adapt F4 to this representation. This gives rise to the small difference from reference data found in the literature. (53) Walsh, M. R.; et al. The Cages, Dynamics, and Structuring of Incipient Methane Clathrate Hydrates. Phys. Chem. Chem. Phys. 2011, 13, 19951−19959. (54) Jacobson, L. C.; Hujo, W.; Molinero, V. Thermodynamic Stability and Growth of Guest-Free Clathrate Hydrates: A LowDensity Crystal Phase of Water. J. Phys. Chem. B 2009, 113, 10298− 10307. (55) Hirai, S.; Okazaki, K.; Tabe, Y.; Kawamura, K. CO2 ClathrateHydrate Formation and Its Mechanism by Molecular Dynamics Simulation. Energy Convers. Manage. 1997, 38, S301−S306. (56) English, N. J.; Lauricella, M.; Meloni, S. Massively Parallel Molecular Dynamics Simulation of Formation of Clathrate-Hydrate Precursors at Planar Water-Methane Interfaces: Insights Into Heterogeneous Nucleation. J. Chem. Phys. 2014, 140, 204714. (57) Agarwal, V.; Peters, B. Advances in Chemical Physics; Wiley: New York, 2014; Vol. 155; pp 97−160.

22857


Methane Clathrate Hydrate Nucleation Mechanism by Advanced

Recommend Documents