Pathways through Equilibrated States with Coexisting Phases for Gas

Subscriber access provided by TUFTS UNIV

Article

Pathways through Equilibrated States with Coexisting Phases for Gas Hydrate Formation Edyta Ma#olepsza, and Tom Keyes J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b06832 • Publication Date (Web): 01 Dec 2015 Downloaded from http://pubs.acs.org on December 9, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Pathways through Equilibrated States with Coexisting Phases for Gas Hydrate Formation Edyta Malolepsza and Tom Keyes∗ Department of Chemistry, Boston University, Boston, MA 02215-2521 E-mail: [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Under ambient conditions water freezes to either hexagonal ice or a hexagonal/cubic composite ice. The presence of hydrophobic guest molecules introduces a competing pathway: gas hydrate formation, with the guests in clathrate cages. Here the pathways of the phase transitions are sought as sequences of states with coexisting phases, using a generalized replica exchange algorithm designed to sample them in equilibrium, avoiding nonequilibrium processes. For a dilute solution of methane in water under 200 atm, initializing the simulation with the full set of replicas leads to methane trapped in hexagonal/cubic ice, while gradually adding replicas with decreasing enthalpy produces the initial steps of hydrate growth. Once a small amount of hydrate is formed, water rearranges to form empty cages, eventually transforming the remainder of the system to metastable β ice, a scaffolding for hydrates. It is suggested that configurations with empty cages are reaction intermediates in hydrate formation when more guest molecules are available. Free energy profiles show that methane acts as a catalyst reducing the barrier for β ice versus hexagonal/cubic ice formation.

Introduction “Water is the most extraordinary substance! Practically all its properties are anomalous...” wrote Albert Szent-Gyorgyi. 1 Water exhibits 72 anomalies, 2 solid structures that strongly depend on temperature and pressure (17 known crystalline forms and possibly still more to be found - ice XVI was experimentally obtained in 2014 3 ), plus amorphous forms of various density. Properties of mixtures of water with other substances add to this list. Nonpolar, hydrophobic compounds do not mix easily with polar water, and solutions usually remain as two-phase systems at ambient temperature and pressure, with low solubility of the minor component in each phase. The solubility of hydrophobic gases, e.g. methane, may be increased by increasing the pressure. For methane vapor in equilibrium with dissolved methane at 0◦ C, one estimate 4 is that the mole fraction, X, is ≈ 4×10−5 at 1 atm and ≈ 2


Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4×10−3 at 200 atm. However, as the methane pressure is increased near the water freezing temperature, something quite different from ordinary hydrophobic solvation can occur. Methane hydrate forms, with β ice providing a crystalline array of water cages, and methane molecules fill them to stabilize the system. We follow the terminology of van der Waals and Plateeuw, 5 who denote the metastable form of the hydrate matrix by “β”. The cages are almost round and are formed by flat pentagons and hexagons of water molecules, arrangements that largely preserve the hydrogen bond network despite changing the structure from the usual arrangement in the liquid and in hexagonal ice. Methane hydrate is important as an alternative energy source 6,7 or potential hazard for climate stability. 8–12 Full occupancy of β ice matrix cages leads to a ratio of 0.174 to 0.176 methanes per water for symmetry sI and sII, respectively, or X ≈ 0.15, corresponding to an increase in solubility over ordinary solvation of ≈ 37× at 200 atm and ≈ 3700× at 1 atm. Hydrate formation in nature does not take place with dissolved methane in equilibrium with the gas. 13 Concentrations on the sea floor may be much lower or, locally, higher. Regardless, hydrate formation starting with a dilute methane concentration very much less than the stoichiometric concentration is an important case. The process will progress until the methane in a region is consumed, and wait for more to enter, naturally pausing for convenient investigation. Computer simulation of first order phase transitions, be they liquid water freezing or hydrate formation, is challenging to conventional methods. We aim to test and demonstrate the power of a new method specifically targeted at first order transitions, the generalized replica exchange method (gREM). 14 While there have now been several applications of the gREM, 14–20 this is the first to mixtures, and the third of our isobaric molecular dynamics version. 19,20 In the gREM, enthalpy, H, and not temperature, is the control variable. To study concentration-limited growth we perform computer simulations on a system with X ≈ 0.008, with too few methanes to allow formation of a crystalline hydrate. Depending

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

on conditions, we find pathways leading either to formation of hexagonal and cubic ice or to partial formation of methane hydrate. In the former case we find nucleation and growth of ices similar to that in pure water, which is not surprising since the system is predominantly water. In the latter case, however, we observe an extraordinary influence of the dilute methanes, with a spontaneous expansion of β ice cages even when there are no more free guest molecules to stabilize them. Almost all the waters are converted to β ice by relatively very few methanes. The goals of this paper are twofold. In addition to developing the gREM, we seek to study the states of a dilute solution of methane, which cannot form a stochiometric hydrate, as a function of enthalpy, bearing in mind recent findings about the formation of ice with hexagonal and cubic layers stacked upon each other 20–28 in the zero concentration limit.

Methods MD gREM Hydrate formation and water freezing are first-order phase transitions. Constant-temperature algorithms are not well suited to such processes, as they cannot sample the states with coexisting phases, which are essential elements of the transition pathway when T ∼ Teq , where Teq is the equilibrium transition temperature. Deep supercooling or a bias potential is required to drive the transition, distorting the results by unknown amounts. In other words, temperature is not a good control variable. 29,30 The difficulties, discussed in detail elsewhere, 14–20 arise as follows. In an isothermal-isobaric ensemble at temperature T , the extrema, H ⋆ , of the enthalpy distribution obey the equation, T = TS (H ⋆ ), where TS (H) = 1/(dS(H)/dH) is the statistical temperature, and S is the configurational entropy. Away from a phase transition there is a single solution corresponding to the most probable enthalpy. At a transition, however, 4


Page 4 of 29

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TS (H) exhibits an S-loop due to surface effects in a finite simulation box, and T = TS (H ⋆ ) has three solutions, corresponding to the high and low enthalpy phases and a barrier, for T in the vicinity of the equilibrium transition temperature, Teq . 29 Cooling from T > Teq the system remains stuck in the high-enthalpy minimum even as T falls below Teq , until a nonequilibrium transition occurs at a temperature unrelated to Teq . The states with coexisting phases followed in the equilibrium pathway are not visited. Water has proved particularly difficult to freeze in simulations. 31 Similar considerations hold for approaching the transition from below. Conventional replica exchange is of limited help, because with constant-T replicas accepted exchanges are rare in the transition region. Thus, we employ an enhanced sampling algorithm based on generalized ensembles, the gREM, 14 specifically the recent isobaric MD version 20 that was implemented by our group in LAMMPS. 32,33 The gREM has a unique capability to sample states with coexisting phases. In the gREM, instead of a constant temperature, each replica α has a generalized temperature function, Tα (H), Tα (H) = λα + ηα (H − H0 ), that depends on enthalpy, H, where H0 is a reference enthalpy. Now the most probable enthalpy, Hα⋆ , obeys Tα (Hα⋆ ) = TS (Hα⋆ ). The slope, ηα , is negative and steep enough to ensure a single intersection of Tα (H) and TS (H) in the transition region, despite the characteristic backbending or S-loop in TS . 14–18,20 The result is that the replica robustly samples configurations with coexisting phases, that form the equilibrium transition pathway, within a unimodal enthalpy distribution. The parameter, λα , determines the intersection point, Hα⋆ , the center of the distribution, and the replicas are chosen to span the coexistence enthalpy range. What is the physical interpretation of the effective temperature? Because T = TS (H ⋆ ) at constant T , and Tα (Hα⋆ ) = TS (H ⋆ ) in the generalized ensemble, a replica with Tα (Hα⋆ ) samples configurations with enthalpies also found in an isothermal-isobaric simulation with T = Tα (Hα⋆ ). That is, the “effective temperature” really does indicate the analogous isothermal5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

isobaric temperature for configurational sampling. However, there is a crucial distinction: when the constant-T simulation would sample an enthalpy distribution with multiple peaks, a gREM replica samples from only one of the peaks, allowing dissection of the coexistence enthalpy range. While the gREM achieves enhanced configurational sampling, a second, constant “kinetic temperature”, T0 , governs ordinary Boltzmann sampling of the velocities. In addition to equilibrium sampling, the gREM allows us to find near-equilibrium pathways for hydrate formation by following the statistical temperature versus enthalpy through the coexisting states. The effective replica temperature is Tα (Hα⋆ ), and in some cases this corresponds to deep supercooling or superheating, so one might be tempted to say that our approach is no different from the conventional one. The point is that the deeply supercooled or superheated states are in equilibrium in the gREM ensemble, or (see below) at least in quasi-equilibrium within a distinct enthalpy landscape metabasin containing a subset of states with the selected enthalpy, and not unstable or metastable with respect to another phase with quite different enthalpy as in a constant-T algorithm. The pathways we find in the present approach follow decreasing enthalpy as determined by the effective temperature functions in the replicas. The route constructed from H ⋆ from neighbouring replicas might be called “walking along the statistical temperature”. We do not choose a temperature to drive hydrate formation, rather, it is naturally selected by the system. Upon passing through the most deeply supercooled states, nucleation and growth proceed with a rising temperature, due to the the S-loop.

Simulation details Simulations with the LAMMPS isobaric MD gREM were performed under 200 atm with four systems of different sizes, but constant concentration X ≈ 0.008 (122-123 waters per methane), denoted 17Met, 27Met, 39Met, 64Met. Details are shown in Table 1. The starting configurations were prepared as follows: a piece of methane clathrate of appropri6


Page 6 of 29

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ate size immersed in liquid water was first melted in the NpT ensemble at 200 atm and 700 K, and then the methane/water solution was cooled and briefly equilibrated at 160 K. At the chosen pressure and concentration, solvated methane is metastable, and almost homogeneously distributed in water (with local density fluctuations). Box sizes varied across the studied enthalpy ranges due to changes in density upon freezing. Water and methane were modeled as single-bead particles interacting via a Stillinger-Weber potential 34 using parameters determined by Molinero. 35,36 While the “mW” potential does not seek to provide the most accurate atomistic description, it reproduces several properties of water remarkably well, 35 provides a rich model for hydrate formation, 36–39 is computationally convenient, and is an appealing model for the first application of the isobaric MD gREM to hydrates. Atomistic potentials will be considered in future work. The gREM is implemented by modifying a standard isothermal-isobaric algorithm in LAMMPS, as explained in Refs. 19 and. 20 The thermostat and barostat were of the NoseHoover type with non-Hamiltonian equations of motion 40–43 with Tdamp and Pdamp set to 100 time steps. We used periodic boundary conditions. The modifications are as follows: 1. The temperature parameter now becomes the kinetic temperature, T0 , which controls the velocity distribution in an ordinary Boltzmann distribution, and was set to 330 K. 2. The forces are scaled by the factor, sα (H) = (T0 /Tα (H)). 3. The kinetic pressure in the barostat is scaled by the inverse factor, s−1 α (H). Statistical temperature curves were calculated by combining results from all replicas and applying the statistical temperature weighted histogram analysis method (ST-WHAM), 44 a Table 1: Characteristics of studied systems with X ≈ 0.008. 17Met 27Met 39Met 64Met 17 methanes 27 methanes 39 methanes 64 methanes 2084 waters 3310 waters 4781 waters 7845 waters 2101 molecules 3337 molecules 4820 molecules 7909 molecules 0.008 methanes per 1 water or 122-123 waters per methane

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

quick and non-iterative alternative to the WHAM method, 45 which converges slowly with increasing system size.

Algorithms to identify types of ice and hydrate cages Ice types. In our earlier studies 20 we proposed an algorithm to assign each water molecule to one of four types: liquid, hexagonal ice, cubic ice, or other. This algorithm takes directly into account the local arrangement of neighbours in the first shell of each water molecule using Steinhardt order parameters, 46 and indirectly in the second shell by using averaged order parameters introduced by Lechner and Dellago. 47 Hydrate cages. However, the algorithm does not distinguish water molecules forming hydrate cages as a new species, all are assigned as liquid. Therefore, following Ohmine’s approach, 48 we have developed another algorithm to identify cages in terms of closed polyhedrons made from water molecules. The algorithm is presented below. Capital letters correspond to panels in Fig. 1, where the algorithm is presented graphically. 1. Build database of pairs from molecules within 3.5 ˚ A. 2. Build database of triplets from pairs that share one molecule. 3. From analysis of molecules shared between triplets build databases of triangles, squares, pentagons, and hexagons, providing that all dihedral angles in each polygon are smaller than 55◦ to obtain flat polygons that are the building blocks of hydrate cages. 4. Make a list of polygon pairs that share an edge. 5. For each polygon pair (gray in A) sharing an edge (red in A) (a) Make a list of single (not shared, green in B) and double (already shared, red in B) edges. (b) Make a list of vertices on double edges.

8


Page 8 of 29

Page 9 of 29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


Results and discussion Statistical temperature TS (H) The statistical temperature curves provide both a quantitative and an intuitive picture of near-equilibrium pathways. Fig. 2 presents TS (H) for the four studied systems. We have found two distinct pathways, and each gave rise to a distinct TS (H). 64Met Temperature [K]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

39Met

27Met

17Met

300

300

270

270

240

240

210

210

-3900

-3750

-3600

-2300

-2200 -1610

-1540

-1000

-960

Enthalpy [eV]

Figure 2: Statistical temperature curves for four system sizes, all with X ≈ 0.008. Black lines correspond to regions where methane hydrate was observed, red lines to mixtures of hexagonal and cubic ice plus methane in liquid water, and blue lines to methane solvated in liquid water. The first pathway occurs when the simulation is started with a set of replicas covering at once the entire interesting range of enthalpy. Then, in replicas with Tα (Hα⋆ ) ≈ TS (Hα⋆ ) < Teq , the system is exposed to conditions of lower generalized temperature than when it was initially equilibrated - it experiences a “generalized quench” - leading to formation of mixtures of hexagonal and cubic ices (Sec. III. C), just as we observed during pure water freezing. 20 Ice forms before the rearrangements needed to make the hydrate are possible. Corresponding statistical temperature curves are shown as red lines in Fig. 2. The blue line represents methane solvated in liquid water. We have also performed simulations with gradual addition of the replicas. Initially some replicas were equilibrated such that Tα (Hα⋆ ) > Teq . We subsequently add a few replicas with lower λ (and therefore lower Tα (Hα⋆ )), each initialized with the final structure from the previously lowest replica, equilibrate, and repeat, until the entire interesting enthalpy region 10


Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


is covered. Despite being an MD algorithm containing kinetics (system evolves according to Newtonian equation of motions with scaled forces), gREM is designed for equilibrium sampling. In the presence of high barriers, the sampling can turn out to be quasi-equilibrium restricted to a specific metabasin on the enthalpy landscape. Here kinetics does enter, insofar as the system cannot escape from the metabasin in the available simulation time. In the three larger systems, the second approach allows the rearrangements needed for hydrate formation without interruption by ice formation, and we find a striking propagation of empty cages (Sec. III.D). Corresponding statistical temperature curves are shown as black lines in Fig. 2. For 17Met, methane hydrate in low generalized temperature replicas spontaneously transformed to ice. Therefore for this system we have a single TS , regardless of initialization, where the red section corresponds to an ice mixture, and the black to the presence of methane hydrate, but without the empty cage expansion. With two distinct pathways, equilibrium must be understood as occurring within metabasins that do not allow exchange on the time scale of our simulations. Nevertheless, the system is closer to equilibrium, sampling coexisting states in a time-independent fashion, than constant-T simulations requiring deep supercooling followed by a strongly driven transition. Furthermore, rather than a difficulty, we consider it desirable to be able to sample distinct, physically meaningful pathways.

System size dependence and pathways Hydrate growth is collective, with a barrier imposing a surface free energy penalty on the formation of small clusters of cages. Here we are deliberately studying systems with low methane concentration, with system sizes such that the total number of methanes, ranging from 17 to 64, is small. Finite-size effects must be considered. It is not surprising that the smallest system, 17Met, does not show significant growth of hydrate structures, specifically of empty cages. The hydrate which is formed initially transform into ice with decreasing H. 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Note that we do no standard quenches with T as the control variable: in our approach, one follows the pathway by the statistical temperature. For the three larger systems, 64Met, 39Met and 27Met, hydrate forms before ice with decreasing temperature (black curve deviates from blue of solvated methane at higher T than does red). While this is a normal scenario for hydrate formation, again, only small hydrate crystallites, and empty cages can form. The presence of empty cages not stabilized by guests, plus surface effects, explain why an ice mixture, plus methane, was found to be more stable than a hydrate-containing system at low H and TS (H) (left hand side of TS curves) having (i) lower enthalpy for a given temperature and (ii) higher temperature as the limit of stability (maximum of TS (H)). With a large enough system, or higher concentration of methane, and water in excess the low-enthalpy state should be hydrate crystal in ice, with hydrate in liquid water found at higher H, and solvated methane in liquid water above the nucleation point. Moving from 27Met to 64Met, with decreasing H the initial rise of the black curves steepens, and proceeds to a higher TS value. The rise in TS corresponds to a fall in the enthalpy of basins visited on the enthalpy landscape, so the stability of the hydrate structures is increasing with system size. Correspondingly the entropy penalty for hydrate formation is decreasing, as a larger temperature enters the denominator of dH/TS (H), making the pathways closer, thermodynamically. Note also that the kinetics of mW water, and the dissipation of latent heat, 50 are too fast, possibly causing ice to form too fast and favoring the red curves. The slower kinetics of atomistic water could further decrease the difference between the two pathways. Further studies are indicated for future work. The hydrate pathway leads to metastable product states. However as they are stabilized with more guests at higher concentration or in larger systems, if formation of empty cages is fast compared to diffusion of methane to occupy them, empty cages will persist as reaction intermediates. We suggest that, when guest availability is limited, configurations with sub-

12


Page 12 of 29

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


stantial numbers of empty cages can act as reaction intermediates - not as stable product states. The TS (H) presented in Fig. 2 represent one of a family of possible curves. As shown earlier, 20 in simulations differing by the random seed used to generate initial velocities, or by different ηα and λα in the replicas, nuclei of slightly different shape and ratio of hexagonal to cubic ices will form, leading to different structures of layered ice, and therefore to different TS (H). We expect similar behaviour for hydrate formation and therefore anticipate some differences in TS (H) depending on starting conditions, with an overall broad similarity. Perfectly reproducible statistical temperatures can be obtained by a melting process when starting from the perfect structure of a solid. 20 Nevertheless, the red curves of all four systems show more or less well defined plateaus at T ≈ 270K, corresponding to a coexisting slab of ice and liquid methane solution. Changing the enthalpy in these configurations by melting or freezing does not change the surface area and therefore does not change the surface effects influencing the shape of TS , hence the plateau, and therefore the plateau temperature gives a quick estimate of the equilibrium temperature which is roughly correct.

Ice formed in the presence of methane A mixture of hexagonal and cubic ice forms in all studied systems if initialized with all the replicas. Representative structures are shown in Fig. 3. Just as for pure water, 20–28 in the case of dilute methane solutions we see formation of layered ice with hexagonal and cubic surfaces stacked on each other. The ratio of hexagonal to cubic ice depends on initial conditions in the simulations, but in general is about 1:1. For Fig. 3 we choose snapshots with some liquid still present to show that indeed watermethane solutions form a pure water solid phase upon freezing, while methane (cyan balls) stays in the remaining liquid until all liquid is gone and methane ends up as a bubble trapped in ice, as shown on the right panel of Fig. 4. When freezing, a slab of ice is growing parallel to 13




Page 14 of 29

Page 15 of 29




The β ice topology follows the same pattern seen in cubic, hexagonal, and β ice melting, 17,20 or in Lennard-Jones fluids, 30 with sphere, cylinder, and slab forming with increasing amount. Table 2: Cages in sI and sII hydrates. cage type number of waters forming a cage number of cages in primary cell of sI number of cages in primary cell of sII

512 20 2 16

512 62 24 6 -

512 64 28 8

The relevant hydrate structures are sI and sII. 6 They contain three types of cages, see Table 2. A cage of type 512 consists of twelve water pentagons, 512 62 means twelve pentagons and two hexagons, while 512 64 means twelve pentagons and four hexagons. The listed cages require 20, 24, or 28 waters, respectively. sI and sII differ not only by cage types, but also by their geometrical arrangement; simple cubic for sI and faced centered cubic for sII. We applied our cage identification algorithm to determine numbers of cage types versus enthalpy along the statistical temperature; results are presented in Fig. 6. The numbers of cages found, compared to the numbers of methanes in the systems, establishes that they are predominantly empty, with approximately 10 times more unoccupied than occupied ones. 64Met

39Met

27Met

600

# cages

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 29

20 24 26 28

400 200 0

-3800

-3700

-2320

-2240

-1600

-1550

Enthalpy [eV]

Figure 6: Number of cages in 64Met, 39Met, and 27Met, as a function of enthalpy. Right to left is the direction of cage formation. In each case the most abundant cage was the 20-water 512 , followed by three times fewer 28-water 512 64 cages of sII, and smaller numbers of the 24-water 512 62 of sI (in perfect sI crystal 512 62 is 3× more abundant than 512 ) and 26-water 512 63 cages that do not enter any crystal but were proposed 39,53–55 as a bridge between 512 62 and 512 64 in early growth. 16


Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Methane ultimately forms sI, so we confirm 36,37,53,54 that initially elements of both structures coexist. In our study of early growth, sII is dominant. Visual analysis of the hydrate structure in Fig. 5 corroborates this conclusion - cages are arranged in a rhomboidal pattern in contrast to the cubic one for sI. 17 This makes sense with empty cage expansion, since the empty sII lattice is more stable than empty sI. 39 Periodic boundary conditions stabilize periodic structures in finite systems, 50,56 so it it worth asking if this influences our findings of empty cage expansion. Note that hexagonal, cubic, and β ice are all periodic, so there is no reason to think that β ice is favored relative to the other crystal forms. In 17Met, the smallest system and presumably most susceptible to this effect, we do not see empty cage expansion; rather, β ice rapidly transforms to a mixture of hexagonal and cubic ice. Furthermore, with gREM and mW, we have accurately described the β ice/liquid transition 17 with 2944 waters, and the cubic and hexagonal ice/liquid transitions, 20 with 3337 and 3168 atoms, respectively, obtaining transition temperatures and the surface tensions. It does not seem that such success would be possible if the crystal was unduly stabilized vis-a-vis the liquid, and 64Met is twice as large. Finally, even if the boundary conditions distort the results when the box is nearly full of metastable β ice, the effect is weaker at earlier stages of the reaction pathway, and initial-intermediate stage propagation of empty cages is how we envision them serving as reaction intermediates in the formation of thermodynamically stable hydrates.

Stability of ice mixture versus hydrate/β ice For a more quantitative comparison of the two competing pathways we calculated the Gibbs free energy profiles versus enthalpy, Fig. 7, at the equilibrium temperature, Teq , for ice mixture formation, determined from a Maxwell construction of equal areas 17,20 on 1/TS (H). Entropy is obtained to within an additive constant as the integral of 1/TS . In order to have 17



the same value of this constant for both hydrate and ice mixture, and to obtain positive values of entropy across the coexistence region, we start the integration at high enthalpy, where both systems are identical solvated methane and the TS curves overlap, and shift the entropy curves up by 2 eV. The entropies are also identical at high enthalpy. Fig. 7 expresses the correct free energy difference of hydrate forming versus ice forming, but, because of the constant entering the entropy, not the absolute values. -3945 hydrate/β ice growth

-3950

fo

rm

at

io

n

-3955 -3960

ice

Gibbs free energy [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 29

-3965 -3900

-3800

-3700

-3600

-3500

Enthalpy [eV]

Figure 7: Gibbs free energy profile calculated as H − Teq (S(H)+2) for 64Met with Teq = 243 K < 273 K due to finite size surface effects, presence of methane and absence of thermodynamically stable pure hexagonal ice. Black and red lines correspond to profiles for hydrate/β ice and hexagonal/cubic ice mixture formation, respectively. Green dashed line shows equal depth of minima for the ice mixture. With decreasing enthalpy, the red line shows a double well potential and free energy barrier for the pathway forming an ice mixture. The black free energy profile for hydrate formation and empty cage expansion, after the initial overlap with the red line for ice, exhibits a lower barrier height and continues to a shallow, high “products” well, reflecting the greater stability of the hexagonal/cubic ice mixture versus β ice. It is as if methane acts as a catalyst, lowering the barrier by inducing β ice formation which acts, as mentioned earlier, as a reaction intermediate. In the finite-size, dilute solutions of methane simulated, the most stable low-enthalpy states are mixtures of hexagonal and cubic ices, plus methane. However, at intermediate enthalpies in the coexistence range, near the barrier top for ice formation, due to the propagation of empty cages, the pathway for forming hydrate and β ice has a lower free energy. 18


Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Thus we expect that the presence of a methane reservoir and the possibility of growing a large piece of hydrate with stabilized, occupied cages will eventually lead to a situation where hydrate is overall the most stable. Several other authors 57–61 have observed formation of small numbers of of empty cages in their simulations, reaching 3 to 5% of all formed cages. They discussed the role of empty cages in the nonequilibrium hydrate growth pathway, and the need for guests and full cages to nucleate and stabilize them. On the whole these studies use strong driving forces. On the other hand, as in our simulations, Ref. 38 finds large numbers of empty cages generated after a critical nucleus has formed, in a simulation at 210 K and 500 atm. The explanation is simply that, under such strong driving conditions, the β ice lattices are stable, and grow too fast for guests to fill the cages, even though additional guests are present. In our study a very large multiplier effect for empty cage expansion induced by occupied cages occurs along a pathway through the coexisting states, with quasi-equilibrium sampling in the replicas versus the usual nonequilibrium process, and with filling of the cages impossible due to low concentration of guest molecules. The stages of empty cage expansion follow the black section of TS , right to left, in Fig. 2. The desirability of considering low concentrations for comparison with real conditions, as it is herein, was recognized in a massively parallel MD study of mW with up to seven million atoms, 62 along with the resulting difficulty 63 of observing extremely slow nucleation rates. With strong driving forces at 225 K and 500 bar, it was found that hydrate formation via heterogeneous nucleation at the water-methane interface became faster with increasing system size. In part this was ascribed to the presence of longer-wavelength fluctuations. We do not address kinetics, our goal is to stay close to equilibrium, and our system is homogeneous and does not possess a fluctuating interface (although of course other fluctuations are important for homogeneous nucleation), but this is consistent with our finding that the hydrate pathway is more favored in larger systems; also, equilibration is faster. Other non-conventional methods, forward flux sampling, 64 and well-tempered metady-

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

namics and restrained molecular dynamics, 65 have been recently applied to hydrate growth with the mW model. This is natural for the reasons given in the Introduction, namely, that simulating a first-order transition without excessive supercooling is challenging. The cited methods depend upon a choice of order parameter or biased potential to guide the simulation through the transition region; results depend upon the choice. While orderparameter dependent pathways are valuable and intuitive, we note that gREM is based on the strategy of sampling states with coexisting phases utilizing the behavior of the statistical temperature - the “caloric curve” - and the behavior of ensembles with enthalpy-dependent generalized temperatures. Analysis 17–20 is based upon the entropic approach 29,30,66,67 in which the pathway is characterized by the fundamental thermodynamics of the entropy as a function of enthalpy (or the inverse entropy derivative, TS (H)). No choice of order parameter is required. Thus the entropic perspective is applied to hydrates.

Conclusions Having constructed an isobaric MD version of an enhanced sampling algorithm, gREM, designed to study first-order phase transitions by sampling in equilibrium states with coexisting phases along the transition pathway, avoiding strongly nonequilibrium processes, we directed our attention to the liquid↔solid transition of water in the presence of a low concentration of a hydrophobic guest molecule, methane. In our approach the pathway is naturally expressed in terms of progress along the statistical temperature function, TS (H). We observed two competing pathways for solidifying at 200 atm, each resulting in a distinct TS (H). When the system is initialized with all the replicas, the system transforms into a hexagonal/cubic ice mixture with trapped methane. If the replicas are added gradually, water molecules around methane can reorganize into clathrate cages, and then induce further changes in water structure leading to the formation of empty cages - β ice. As such, at low concentrations small numbers of methanes can induce many more empty cages. If they

20


Page 20 of 29

Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


persist until more methane enters to fill and stabilize them, they will induce more empty cage formation, etc. Thus the relative rates of empty cage growth and guest diffusion to fill the cages will be a significant factor in hydrate growth. Calculated free energy profiles for both processes corroborate this conclusion - with a small amount of methane the most stable form is a hexagonal/cubic ice mixture, but even at such low concentration methane makes the system ready to form methane hydrate and β ice by reducing the barrier below the barrier for ice formation. Methane behaves as a catalyst. The product state of the hydrate pathway is metastable and not observable in equilibrium. Beta ice is not the stable form for any conditions of interest, as shown in Ref. 39 and must eventually change to mixed cubic/hexagonal ice and then hexagonal ice. However, once β ice growth begins with only the higher-enthalpy replicas present, the simulations proceed to sample the metabasin of clathrate as lower-H replicas are added. Despite the enhanced sampling of the gREM, transitions between different crystal structures are too much to expect. This picture is consistent with Fig. 7. When cage formation is not interrupted by ice formation, the barrier is lowered. Empty cage propagation occurs near the barrier top, where one one expects to see intermediates, but the pathway proceeds to a state which is metastable with respect to layered ice plus trapped methane. In experiment, eventually, stable ice would form or more methane would enter to stabilize hydrate growth. Our suggestion is that the ease with which empty cages propagate before such events occur indicates that there may be near-equilibrium pathways in which they play a major role as reaction intermediates. Large numbers of empty cages were not found in Ref. 63 studying conditions of low driving force. However, clathrate seeds much larger than any cluster in our systems were introduced, indicating a later stage of growth. gREM studies of larger systems are indicated to probe those later growth stages and product states. The present application of gREM to hydrates has already revealed considerable information about the statistical temperatures, quasi-equilibrium pathways, barriers and reaction intermediates.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acknowledgement The research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences under Award DESC0008810. We thank Profs. William Klein and Harvey Gould and Dr. Jaegil Kim for helpful discussions.

References (1) Szent-Gyorgyi, A. The Living State: With Observations on Cancer ; Academic Press, Inc., New York, 1972. (2) http://www1.lsbu.ac.uk/water/water anomalies.html (accessed Nov 9, 2015). (3) Falenty, A.; Hansen, T. C.; Kuhs, W. F. Formation and Properties of Ice XVI Obtained by Emptying a Type sII Clathrate Hydrate. Nature 2014, 516, 231–233. (4) Duan, Z.; Moller, N.; Greenberg, J.; Weare, J. H. The Prediction of Methane Solubility in Natural Waters to High Ionic Strength from 0 to 250◦ C and from 0 to 1600 bar. Geochim. Cosmochim. Acta 1992, 56, 1451–1460. (5) van der Waals, J. H.; Platteeuw, J. C. Clathrate Solutions. Adv. Chem. Phys. 1959, 2, 1–57. (6) Sloan, E. D.; Koh, C. Ed. Clathrate Hydrates of Natural Gases; CRC Press/TaylorFrancis, Boca Raton, FL, 2007. (7) Koh, C.; Sum, A.; Sloan, E. D. Gas Hydrates: Unlocking the Energy from Icy Cages. J. Appl. Phys. 2009, 106, 061101. (8) Kennett, J.; Cannariato, K.; Hendy, I.; Behl, R. Methane Hydrates in Quaternary Climate Change: The Clathrate Gun Hypothesis. Spec. Publ. Ser. American Geophysical Union 2003, 54, 1–217. 22


Page 22 of 29

Page 23 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9) Winguth, A.; Shellito, C.; Shields, C.; Winguth, C. Climate Response at the Paleocene– Eocene Thermal Maximum to Greenhouse Gas Forcing—A Model Study with CCSM3. J. Clim. 2010, 23, 2562–2584. (10) Phrampus, B.; Hornbach, M. Recent Changes to the Gulf Stream Causing Widespread Gas Hydrate Destabilization. Nature 2012, 490, 527–530. (11) Shakhova, N.; Semiletov, I.; Leifer, I.; Sergienko, V.; Salyuk, A.; Kosmach, D.; Chernykh, D.; Stubbs, C.; Nicolsky, D.; et al., V. T. Ebullition and Storm-Induced Methane Release from the East Siberian Arctic Shelf. Nat. Geosci. 2014, 7, 64–70. (12) Kennedy, M.; Mrofka, D.; von der Borch, C. Snowball Earth Termination by Destabilization of Equatorial Permafrost Methane Clathrate. Nature 2008, 453, 642–645. (13) E. D. Sloan, C. K.; Sum, A. Gas Hydrate Stability and Sampling: the Future as Related to the Phase Diagram. Energies 2010, 3, 1991–2000. (14) Kim, J.; Keyes, T.; Straub, J. E. Generalized Replica Exchange Method. J. Chem. Phys. 2010, 132, 224107. (15) Lu, Q.; Kim, J.; Straub, J. E. Order Parameter Free Enhanced Sampling of the VaporLiquid Transition Using the Generalized Replica Exchange Method. J. Chem. Phys. 2013, 138, 104119. (16) Lu, Q.; Kim, J.; Straub, J. E. Exploring the Solid-Liquid Phase Change of an Adapted Dzugutov Model Using Generalized Replica Exchange Method. J. Phys. Chem. B 2012, 116, 8654–8661. (17) Malolepsza, E.; Kim, J.; Keyes, T. Entropic Description of Gas Hydrate Ice/Liquid Equilibrium via Enhanced Sampling of Coexisting Phase. Phys. Rev. Lett. 2015, 114, 170601.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(18) Malolepsza, E.; Kim, J.; Keyes, T. Generalized Ensemble Method Applied to Study Systems with Strong First Order Transitions. J. Phys.: Conf. Ser. 2015, 640, 012003. (19) Secor, M.; Malolepsza, E.; Keyes, T. Isobaric Molecular Dynamics Version of the Generalized Replica Exchange Method: Liquid-Vapor Equilibrium. J. Phys. Chem. B 2015, 119, 13379–13384. (20) Malolepsza, E.; Keyes, T. Water Freezing versus Ice Melting. J. Chem. Theor. Comp. 2015, in press, DOI: 10.1021/acs.jctc.5b00637. (21) Li, T.; Donadio, D.; Russo, G.; Galli, G. Homogeneous Ice Nucleation from Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 19807–19813. (22) Pirzadeh, P.; Kusalik, P. G. On Understanding Stacking Fault Formation in Ice. J. Am. Chem. Soc. 2011, 133, 704–707. (23) Moore, E. B.; Molinero, V. Is It Cubic? Ice Crystallization From Deeply Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 20008–20016. (24) Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of Ice Crystallized From Supercooled Water. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1041–1045. (25) Kuhs, W. F.; Sippel, C.; Falenty, A.; Hansen, T. C. Extent and Relevance of Stacking Disorder in “Ice Ic ”. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 21259–21264. (26) Quigley, D. Communication: Thermodynamics of Stacking Disorder in Ice Nuclei. J. Chem. Phys. 2014, 141, 121101. (27) Carr, T. H. G.; Shephard, J. J.; Salzmann, C. G. Spectroscopic Signature of Stacking Disorder in Ice I. J. Phys. Chem. Lett. 2014, 5, 2469–2473. (28) Malkin, T. L.; Murray, B. J.; Salzmann, C. G.; Molinero, V.; Pickering, S. J.; Whale, T. F. Stacking Disorder in Ice I. Phys. Chem. Chem. Phys. 2015, 17, 60–76. 24


Page 24 of 29

Page 25 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(29) Bachmann, M. Thermodynamics and Statistical Mechanics of Macromolecular Systems; Cambridge University Press, Cambridge, UK, 2014. (30) Binder, K.; Block, B.; Virnau, P. Beyond the Van Der Waals Loop: What Can Be Learned from Simulating Lennard-Jones Fluids Inside the Region of Phase Coexistence. Am. J. Phys. 2012, 80, 1099–1109. (31) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular Dynamics Simulation of the Ice Nucleation and Growth Process Leading to Water Freezing. Nature 2002, 416, 409–413. (32) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1–19. (33) LAMMPS. http://lammps.sandia.gov (accessed Nov 9, 2015). (34) Stillinger, F. H.; Weber, T. A. Computer Simulation of Local Order in Condensed Phases of Silicon. Phys. Rev. B 1985, 31, 5262–5271. (35) Molinero, V.; Moore, E. B. Water Modeled As an Intermediate Element between Carbon and Silicon. J. Phys. Chem. B 2009, 113, 4008–4016. (36) 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. (37) Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 2010, 132, 11806–11811. (38) Jacobson, L. C.; Matsumoto, M.; Molinero, V. Order Parameters for the Multistep Crystallization of Clathrate Hydrates. J. Chem. Phys. 2011, 135, 074501. (39) Jacobson, L. C.; Hujo, W.; Molinero, V. Thermodynamic Stability and Growth of Guest-Free Clathrate Hydrates: A Low-Density Crystal Phase of Water. J. Phys. Chem. B 2009, 113, 10298–10307. 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(40) Martyna, G. L.; Tobias, D. J.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177–4189. (41) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. (42) Tuckerman, M. E.; Alejandre, J.; Lopez-Rendon, R.; Jochim, A. L.; Martyna, G. J. A Liouville-Operator Derived Measure-Preserving Integrator for Molecular Dynamics Simulations in the Isothermal–Isobaric Ensemble. J. Phys. A: Math. Gen. 2006, 39, 5629–5651. (43) Shinoda, W.; Shiga, M.; Mikami, M. Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation under Constant Stress. Phys. Rev. B 2004, 69, 134103. (44) Kim, J.; Keyes, T.; Straub, J. E. Communication: Iteration-Free, Weighted Histogram Analysis Method in Terms of Intensive Variables. J. Chem. Phys. 2011, 135, 161103. (45) Ferrenberg, A. M.; Swendsen, R. H. Optimized Monte Carlo Data Analysis. Phys. Rev. Lett. 1989, 63, 1195–1198. (46) Steinhardt, P.; Nelson, D. R.; Ronchetti, M. Bond-Orientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, 784–805. (47) Lechner, W.; Dellago, C. Accurate Determination of Crystal Structures Based on Averaged Local Bond Order Parameters. J. Chem. Phys. 2008, 129, 114707. (48) Matsumoto, M.; Baba, A.; Ohmine, I. Topological Building Blocks of Hydrogen Bond Network in Water. J. Chem. Phys. 2007, 127, 134504. (49) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 2011, 32, 2319–2327.

26


Page 26 of 29

Page 27 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(50) English, N. Massively Parallel Molecular-Dynamics Simulation of Ice Crystallisation and Melting: The Roles of System Size, Ensemble, and Electrostatics. J. Chem. Phys. 2014, 141, 234501. (51) Shultz, M. J.; Bisson, P. J.; Brumberg, A. Best Face Forward: Crystal-Face Competition at the Ice–Water Interface. J. Phys. Chem. B 2014, 118, 7972–7980. (52) Poon, G.; Peters, B. A Stochastic Model for Nucleation in the Boundary Layer during Solvent Freeze-Concentration. Cryst. Growth Des. 2013, 13, 4642–4647. (53) Walsh, M.; Koh, C.; Sloan, E. D.; Sum, A.; Wu, D. Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326, 1095–1098. (54) Benedetti, P. D.; Sarupia, S. Hydrate Molecular Ballet. Science 2009, 326, 1070–1071. (55) Vatamanu, J.; Kusalik, P. Observation of Two-Step Nucleation in Methane Hydrates. Phys. Chem. Chem. Phys. 2010, 12, 15065–15072. (56) English, N.; Tse, J. Massively Parallel Molecular Dynamics Simulation of Formation of Ice-Crystallite Precursors in Supercooled Water: Incipient-Nucleation Behavior and Role of System Size. Phys. Rev. E 2015, 92, 032132. (57) Sarupria, S.; Debenedetti, P. G. Homogeneous Nucleation of Methane Hydrate in Microsecond Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2012, 3, 2942–2947. (58) Bai, D.; Chen, G.; Zhang, X.; Wang, W. Nucleation of the CO2 Hydrate from Three– Phase Contact Lines. Langmuir 2012, 28, 7730–7736. (59) Pirzadeh, P.; Kusalik, P. G. Molecular Insights into Clathrate Hydrate Nucleation at an Ice–Solution Interface. J. Am. Chem. Soc. 2013, 135, 7278–7287. ´ (60) Jim´enez-Angeles, F.; Firoozabadi, A. Nucleation of Methane Hydrates at Moderate Subcooling by Molecular Dynamics Simulations. J. Phys. Chem. C 2014, 118, 11310– 11318. 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(61) Zhang, Z.; Walsh, M. R.; Guo, G.-J. Microcanonical Molecular Simulations of Methane Hydrate Nucleation and Growth: Evidence that Direct Nucleation to sI Hydrate Is Among the Multiple Nucleation Pathways. Phys. Chem. Chem. Phys. 2015, 17, 8870– 8876. (62) English, N.; 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. (63) Knott, B.; Molinero, V.; Doherty, M.; Peters, B. Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions. J. Am. Chem. Soc. 2012, 134, 19544. (64) Bi, Y.; Li, T. Probing Methane Hydrate Nucleation through the Forward Flux Sampling Method. J. Phys. Chem. B 2014, 118, 13324–13332. (65) Lauricella, M.; Meloni, S.; English, N.; Peters, B.; Ciccotti, G. Methane Clathrate Hydrate Nucleation Mechanism by Advanced Molecular Simulations. J. Phys. Chem. C 2014, 118, 22847–22857. (66) Schnabel, S.; Seaton, D.; Landau, D.; Bachmann, M. Microcanonical Entropy Inflection Points: Key to Systematic Understanding of Transitions in Finite Systems. Phys. Rev. E 2011, 84, 011127. (67) MacDowell, L. G.; Virnau, P.; M¨ uller, M.; Binder, K. The Evaporation/Condensation Transition of Liquid Droplets. J. Chem. Phys. 2004, 120, 5293–5308.

28


Page 28 of 29

Page 29 of 29



Pathways through Equilibrated States with Coexisting Phases for Gas

Recommend Documents