Molecular Simulation of the Phase Diagram of Methane Hydrate: Free

Subscriber access provided by UNIV OF WESTERN ONTARIO

Article

Molecular Simulation of the Phase Diagram of Methane Hydrate: Free Energy Calculations, Direct Coexistence Method, and Hyper Parallel Tempering Dongliang Jin, and benoit coasne Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.7b02238 • Publication Date (Web): 09 Aug 2017 Downloaded from http://pubs.acs.org on August 14, 2017

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.

Langmuir 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 50

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

Langmuir

Molecular Simulation of the Phase Diagram of Methane Hydrate: Free Energy Calculations, Direct Coexistence Method, and Hyper Parallel Tempering Dongliang Jin and Benoit Coasne∗ Laboratoire Interdisciplinaire de Physique (LIPhy), CNRS and Universit´e Grenoble Alpes, F-38000 Grenoble, France E-mail: [email protected]

1

ACS Paragon Plus Environment

Langmuir

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 Different molecular simulation strategies are used to assess the stability of methane hydrate under various temperature and pressure conditions. First, using two water molecular models, free energy calculations consisting of the Einstein molecule approach in combination with semi-grand Monte Carlo simulations are used to determine the pressure–temperature phase diagram of methane hydrate. With these calculations, we also estimate the chemical potentials of water and methane and methane occupancy at coexistence. Second, we also consider two other advanced molecular simulation techniques that allow probing the phase diagram of methane hydrate: the direct coexistence method in the Grand Canonical ensemble and the hyper parallel tempering Monte Carlo method. These two direct techniques are found to provide stability conditions that are consistent with the pressure–temperature phase diagram obtained using rigorous free energy calculations. The phase diagram obtained in this work, which is found to be consistent with previous simulation studies, is close to its experimental counterpart provided the TIP4P/Ice model is used to describe the water molecule.

2


Page 2 of 50

Page 3 of 50

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

Langmuir

1.

INTRODUCTION

Methane hydrate is a non stoichiometric crystalline structure made up of water molecules forming hydrogen–bonded cages around a methane molecule. 1,2 Three primary crystalline hydrate structures have been identified: sI, 3 sII, 4 and sH. 5 These three structures differ from each other in their number and type of water cages (made up of five- and six-member rings of water molecules). Under typical environmental conditions where hydrates are encountered on Earth, methane hydrate exists as structure sI. 2,6 In this crystalline structure, 46 water molecules form two small pentagonal dodecahedral cages (512 ) and six tetracaidecahedral cages (512 62 ) so that a maximum of 8 methane molecules can be encapsulated. 6 Abundant methane hydrate resources on Earth, especially in deep seafloors and in the permafrost, 2,7,8 are important both for energy and environmental applications. 9–13 In particular, in the context of climate change and global warming, even a small temperature increase could induce the melting of methane hydrates and, therefore, the release of large amounts of methane into the atmosphere (methane leads to a far larger greenhouse gas effect than carbon dioxide). 14 Moreover, methane hydrate formation in oil and gas pipelines is known to be detrimental as it hinders flow. Finally, hydrates including methane hydrates are also thought to be a key ingredient in the geochemistry of planets, comets, etc. where the coexistence of water and gases leads to hydrate formation depending on temperature and pressure. 15,16 From a fundamental point of view, methane hydrate and other gas hydrates are model systems to gain insights into the complex thermodynamics and dynamics of non stoichiometric structures including the large family of clathrates. For instance, many porous materials such as zeolites and metal organic frameworks are synthesized by crystallizing cages around an organic template, therefore sharing some important features with hydrates. In addition, owing to their non stoichiometric nature, gas hydrates can be considered as prototypical examples of confined solids which also possess varying compositions 3


Langmuir

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

with temperature and pressure (even though their bulk counterpart exhibits constant stoichiometry). 17,18 As a result, owing to its importance for both fundamental and practical sciences, methane hydrate receives increasing attention with significant effort devoted to better understanding their physical and physicochemical properties. 2,19–28 In parallel to experimental studies devoted to determining stability conditions for methane hydrate, 2,29–32 many theoretical works have focused on the description of its pressure–temperature phase diagram. 19–22,33 Theoretical approaches relying on molecular simulation have been also used to derive appropriate order parameters to distinguish hydrates from cubic ice, hexagonal ice, and water, 26,34–45 and to investigate nucleation/dissociation of methane hydrate. 19,21,22,24,25,42,45–47 The phase diagram of methane hydrate is often predicted using the well-established semi empirical statistical mechanics theory of van der Waals and Platteeuw. 19,21,22,48 From a molecular simulation viewpoint, free energy strategies such as the Einstein crystal or Einstein molecule approaches have been used to determine the free energy and methane occupancy of bulk methane hydrate. 21,22,49–52 These calculations make it possible to determine the phase diagram of methane hydrate in equilibrium with liquid water (or ice at low temperature) and methane vapor using molecular simulation. As another possible strategy, other authors have estimated methane hydrate stability limits using the direct coexistence method; 6,53 With this technique, methane hydrate, liquid water, and methane vapor are initially placed into a simulation box, and several isobaric-isothermal molecular dynamics runs are performed to determine the stable phase at given pressure and temperature conditions. While this method allows in theory to estimate the liquid–hydrate–vapor (L–H–V) boundary conditions, it often requires very long computational times and tedious averaging over initial states. In this paper, we consider different molecular simulation strategies to assess the thermodynamics of bulk methane hydrate. First, for two different water models – namely TIP4P/2005 and TIP4P/Ice –, free energy calculations based on the Einstein molecule approach developed by Vega and coworkers 50,52 are used to determine the pressure–

4


Page 4 of 50

Page 5 of 50

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

Langmuir

temperature phase diagram of methane hydrate (in all simulations, methane is treated using a coarse-grained model known as the united-atom model). More precisely, we determine the stability conditions for 3 pressures: 1, 10 and 100 atm. For each pressure, in addition to determining the temperature range where methane hydrate is stable, we also estimate the methane occupancy of the hydrate and discuss the non negligible effect of the approximation used to treat methane vapor (exact equation of state as probed using molecular simulation versus thermodynamic integration from a supercritical ideal gas). While free energy calculations obviously constitute the most rigorous scheme to determine the phase diagram of such complex phases, we also consider in a second step less demanding strategies. First, we consider the direct coexistence method in which one generates an initial configuration where both the liquid and methane hydrate coexist to determine usnig molecular simulation the final, stable phase for many temperature and pressure conditions. While the direct coexistence method has already been used to investigate the thermodynamic stability of methane hydrate, 6,19 here a novel version is proposed; both water and methane are treated in the Grand Canonical ensemble using Monte Carlo simulations to account for large variations in the number of molecules upon melting and formation of the hydrate. Second, we also consider hyper parallel tempering molecular simulations in which several replicas of the system, taken at different temperatures and chemical potentials, are considered in parallel (following the work by De Pablo and coworkers, these simulations are referred to hyper parallel tempering rather than parallel tempering as the system is treated in the Grand Canonical ensemble). 54–56 While this method has been already used for simulating solid–liquid phase diagrams of confined mixtures, 17,57 it is the first time that such a hyper parallel tempering strategy is considered for methane hydrate. A recent paper, published after initial submission of the present article, has reported a very similar study on the free energy and phase diagram of methane and carbon dioxide hydrates. 58 In this paper, these authors also used the Einstein crystal approach and

5


Langmuir

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

compared their results with experimental 59,60 and simulation 6,19,61–63 data obtained in previous works. While this study is similar to that reported in our paper, there are several differences which make it very interesting to compare the different sets of data/strategies. First, in the free energy calculations, while Waage et al. used an equation of state for methane vapor (PR-EOS), we performed molecular simulations in both the NPT and µVT ensembles to obtain the phase diagram of methane vapor (on the other hand, while we used the Gibbs-Duhem integration to estimate the chemical potential for liquid water, Waage et al. used thermodynamic integration). Second, while almost others in Ref. ( 6,19,63 ) considered the direct coexistence method in the isobaric-isothermal ensemble, we extended this technique to the Grand Canonical ensemble (for the reasons discussed above). Moreover, as discussed earlier, we also considered the Hyper Parallel Tempering technique as an additional direct simulation technique to probe the phase diagram of gas hydrate. Overall, as will be discussed later, the data reported in the present work were found to be fully consistent with those obtained by Waage et al. The remainder of this paper is organized as follows. In Section 2, we present the computational details: molecular models and interaction potentials, Monte Carlo algorithm to generate methane hydrate with structure sI, Molecular Dynamics simulations, and Monte Carlo simulations in different thermodynamic ensembles (Canonical, Semi-Grand Canonical, and Grand Canonical ensembles). In Section 3, we present general considerations regarding the liquid–hydrate–vapor phase equilibrium. In Section 4, free energy calculations of methane hydrate are first presented to determine the phase diagram of methane hydrate for the two water models selected in this work. In this part, we also determine the chemical potential for each species as well as methane occupancy for the different pressure/temperature coexistence conditions. In Section 4, we also present the stability conditions obtained using the direct coexistence method and the hyper parallel tempering method. The results obtained using the different methods above are compared with experimental data as well as data obtained in previous theoretical works. Section 5

6


Page 6 of 50

Page 7 of 50

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

Langmuir

presents some concluding remarks.

2.

COMPUTATIONAL DETAILS

2.1.

Models and interaction potentials

Methane was modeled as a single Lennard-Jones (LJ) sphere with the parameters taken from the OPLS-UA forcefield (UA stands for united-atom). 64,65 Water was modeled using the TIP4P model which consists of a rigid model containing 4 sites: an LJ site located on the oxygen atom, two sites corresponding to the hydrogen atoms, and a fourth site M corresponding to the negative charge of the oxygen atom located at a distance dOM from the oxygen atom toward the hydrogen atoms along the H–O–H angle bisector. Two versions of the TIP4P water model, 66 namely TIP4P/2005 67 and TIP4P/Ice 68 models, were used to describe the water molecules in methane hydrate. In both water models, the water ˚ and an H–O–H angle of 104.52◦ . The LJ molecule has an O–H bond length of 0.9572 A potential parameters for methane and water as well as the atomic charges and distance dOM for the two water models are given in Table 1. The TIP4P/2005 model reproduces qualitatively the liquid/solid coexistence but with a shift in temperature (20–30 K) and in pressure (100 MPa). 66,69 In contrast, the TIP4P/Ice model accurately reproduces the liquid/solid phase diagram but with some limitations in the coexistence lines for some dense ice forms (like Ice VII and Ice VIII). 66 Table 1: Interaction potential parameters corresponding to the OPLS-UA model for methane and the TIP4P/2005 and TIP4P/Ice models for water. For the two water models, we also indicate the melting temperature Tm as predicted from molecular modeling. Model TIP4P/Ice TIP4P/2005 CH4

ǫ/k B (K) 106.1 93.2 147.5

˚ σ (A) 3.1668 3.1589 3.7300

q H (e) 0.5879 0.5564 –

qO (e) -1.1758 -1.1128 –

˚ dOM (A) 0.1577 0.1546 –

Tm (K) 272.2 252.1 –

In our molecular simulation of methane hydrate, the total potential energy is the 7


Langmuir

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 8 of 50

sum of all intermolecular interactions (note that there is no intramolecular contribution since we only consider rigid water models and a united atom model for methane). The intermolecular potential between two atoms i and j includes a short-range repulsive contribution and an attractive dispersion contribution described using the LJ potential: uijLJ (r )

= 4ε ij

"

σij r

12

−

σij r

6 #

(1)

where r is the distance between the two atoms i and j while ε ij and σij are the corresponding LJ parameters. In our calculations, the LJ interaction potentials are truncated with a cut-off (rc = 0.85 nm). While like-atom LJ parameters are given in Table 1, the parameters between unlike atoms are obtained using the Lorentz–Berthelot mixing rules, ε ij = (ε ii ε jj )1/2 and σij = (σii + σjj )/2. 70,71 In addition to the repulsion/dispersion interactions, the intermolecular potential includes the the Coulomb potential between i and j separated by a distance r: uijC =

1 qi q j 4πε o r

(2)

where qi and q j are the atomic charges on atoms i and j (see Table 1). The Coulomb energy was computed using the Ewald summation technique to correct for the small size of the simulation box . 72–75

2.2.

Molecular structure of methane hydrate

Figure 1 shows a molecular configuration of methane hydrate corresponding to 2 × 2

× 2 unit cells of the sI structure (the unit cell has a length of 1.1877 nm). In this section, we describe the strategy used to generate such a molecular configuration of methane hydrate from the experimental crystallographic data. For methane hydrate, three criteria should be verified (more details can be found in Section S1 and Figure S1 of the Supporting Information): (1) proton disorder, (2) ice rules also known as Bernal–Fowler rules, and (3) zero dipole moment. To build a molecular structure obeying these criteria, we followed 8


Page 9 of 50

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

Langmuir

the stochastic strategy proposed by Buch et al. 76 1. A cubic box with dimensions L x = Ly = Lz = 2.3754 nm, corresponding to 2 × 2 × 2 unit cells, is constructed by first placing the oxygen atoms according to the experimental X-ray crystallographic data. 77 2. In order to comply with the ice rule, each pair of nearest neighbor oxygens must share a hydrogen atom which belongs either to the first or second oxygen atom. In what follows, the two oxygen atoms in each O-O pair are labelled O1 and O2 . Initially, a hydrogen atom is randomly assigned either to O1 or O2 for each O-O pair. The distance from the selected oxygen atom to this hydrogen atom is set according to the chemical O–H bond length of the TIP4P water model, dOH = 0.09578 nm. Due to the random assignment of the hydrogen atoms, the initial structure obtained according to this strategy is unrealistic; oxygen atoms are coordinated to Nc = 0, 1, 2, 3 or 4 hydrogen atoms (obviously, coordination numbers Nc 6= 2 are not physical). 3. The following stochastic/Monte Carlo approach is then performed to relax these nonphysical coordination numbers and reach realistic configurations where Nc = 2 for all oxygen atoms. We randomly choose a O-O pair. If the hydrogen atom is bonded to O1 (O2 ), attempt is made to transfer the hydrogen atom to O2 (O1 ). This move is accepted or rejected based on the change in the absolute difference in coordination numbers ∆Nc = | NcO1 − NcO2 |. More precisely, the move is accepted if the change in the absolute difference in coordination numbers ∆(∆Nc ) < 0 (because this leads overall to configurations with oxygen atoms having the same coordination numbers i.e. Nc = 2). The move is accepted with a probability 0.5 if ∆(∆Nc ) = 0. In contrast, the move is rejected if ∆(∆Nc ) > 0. Such moves are attempted until each oxygen atom is linked to two hydrogen atoms (in practice, 20000 moves are performed as it is found sufficient to reach physical configurations for the system size considered in this work). 9


Langmuir

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

Page 10 of 50

Page 11 of 50

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

Langmuir

and pressures P (details of the free energy calculations will be discussed later in this paper). All MD simulations were performed using LAMMPS; 78 The Velocity-Verlet algorithm 79 was used to integrate the equation of motion with a total time of at least 1 ns and a timestep of 1 fs. The Molecular dynamics simulations were performed in the isobaric–isothermal ensemble (constant number of particles N, temperature T, and pressure P) in which the temperature and pressure were controlled using Nose-Hoover thermostat/barostat with a typical relaxation time of 2 ps. 80,81 Monte Carlo simulation in the canonical ensemble (CMC) were used in our free energy calculations to determine (1) the free energy change between the non-interacting and the interacting Einstein molecules ∆A1 and (2) the free energy change from the Einstein molecule to the methane hydrate ∆A2 (again, details of the free energy calculations will be given later). In these canonical simulations (constant number of particles N, temperature T, and volume V), MC moves include rotations for the water molecules and translations for the water and methane molecules. In the framework of the Metropolis algorithm, each move from an old (o) to a new (n) microscopic states was accepted or rejected according to the acceptance probability Pacc = min{1, pnNVT /poNVT } where p NVT for a given configuration corresponds to the density of states in the canonical ensemble: VN p NVT (s ) ∝ exp N! N

−U ( s n ) kB T

(3)

where s N is the set of coordinates of the N molecules in a given microscopic configuration and U (s N ) is the corresponding intermolecular potential energy. Semi-Grand Monte Carlo (SGMC) simulations were performed to determine the numH inside the methane hydrate as a function of their chemical ber of methane molecules Nm

potential µH m at given T and P (here, the subscript m refers to methane while the superscript H refers to the hydrate phase). In this hybrid ensemble, methane is treated at constant chemical potential µH m and temperature T while water is treated at constant number of molecules

11


Langmuir

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 12 of 50

NwH and temperature T. On the other hand, the volume V is allowed to fluctuate since the system is at constant pressure P. For each T and P, we start from an equilibrium configuration obtained using isobaric-isothermal MD simulations. MC moves in SGMC simulations include rotations and translations for water and translations, insertions, and deletions for methane. Moreover, volume changes are also attempted. In the framework of the Metropolis algorithm, moves from an old (o) to a new (n) microscopic states are accepted or rejected according to the acceptance probability Pacc = min{1, pnµm Nw PT /poµm Nw PT } where pµm Nw PT for a given configuration corresponds to the density of states in the semi-grand canonical ensemble: VN pµm Nw PT (s ) ∝ exp N! N

− PV kB T

exp

Nm µm kB T

−U ( s n ) exp kB T

(4)

As in the case of canonical Monte Carlo simulations, s N is the set of coordinates for the N molecules in the microscopic configuration while U (s N ) is the corresponding intermolecular potential energy. V and Nm are the volume and number of methane molecules in the configuration. Grand Canonical Monte Carlo simulations (GCMC) were used in the direct coexistence method and the hyper parallel tempering technique. In the grand canonical ensemble, the system has a constant volume and methane and water are at constant chemical potentials µm , µw and temperature T. Monte Carlo moves in the grand canonical ensemble include rotations, translations, insertions, and deletions for both water and methane. In this ensemble, moves from an old (o) to a new (n) microscopic states are accepted or rejected using a Metropolis scheme with an acceptance probability Pacc = min{1, pnµm µw VT /poµm µw VT } where pµm µw VT for a given configuration corresponds to the density of states in the grand canonical ensemble: VN exp pµm µw VT (s ) ∝ N! N

Nm µm + Nw µw kB T

12

−U ( s n ) exp kB T


(5)

Page 13 of 50

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

Langmuir

s N is the set of coordinates of the N molecules in the microscopic configuration while U (s N ) is the corresponding intermolecular potential energy. Nw and Nm are the numbers of water and methane molecules in the configuration.

3. 3.1.

LIQUID–HYDRATE–VAPOR EQUILIBRIUM Phase coexistence condition

Methane hydrate (H) is a binary mixture of water (w) and methane (m) that coexists with liquid water (L) (or ice at sufficient low T) and methane vapor (V) in specific temperature T and pressure P ranges (i.e., for a given P, there exists a T at which the three phases L–H–V coexist – the hydrate phase being stable at low T/high P). At P and T where the three phases coexist, the chemical potentials µiΦ for each species (i = w, m) in all phases (Φ = L, H, V) are equal. µiΦ at given T and P varies with methane and water mole fractions (xm and xw , respectively) so that L–H–V equilibrium depends also on xm and xw . 29,30 Since xw = 1 − xm for a binary system, the L–H–V equilibrium condition can be expressed using xm only: V µLw ( xm , T, P) = µH w ( xm , T, P ) = µw ( xm , T, P )

µLm ( xm , T, P)

=

µH m ( xm , T, P )

=

(6)

µV m ( xm , T, P )

Such L–H–V equilibrium can be recast as 2 two-phase coexistence conditions: (1) liquid water–methane hydrate (L–H) and (2) methane hydrate–methane vapor (H–V): µLw ( xm , T, P) = µH w ( xm , T, P ) µH m ( xm , T, P )

=

(7)

µV m ( xm , T, P )

As indicated by the experimental Henry constant (xm ∼0.003–0.001 for methane in liquid water at 100 bar for T ranging between 275 and 310 K), 82,83 the solubility of methane 13


Langmuir

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 14 of 50

in liquid water is very low so that the effect of methane on the chemical potential of water in the liquid phase can be neglected, i.e. µLw ( xm ∼ 0, T, P) ∼ µLw ( xm = 0, T, P). 20,21 Similarly, the chemical potential of methane in the vapor can be approximated by that of V pure methane vapor, i.e. µV m ( xm ∼ 1, T, P ) ∼ µm ( xm = 1, T, P ). With these approximations,

the L–H–V coexistence conditions defined in Eqs. (7) become µLw ( xm = 0, T, P) = µH w ( xm , T, P )

(8)

V µH m ( xm , T, P ) = µm ( xm = 1, T, P )

The description above shows that determining phase coexistence requires to estimate H L the four following chemical potentials: µH m ( xm , T, P ), µw ( xm , T, P ), µw ( xm = 0, T, P ), and

µV m ( xm = 1, T, P ).

3.2.

Estimation of the different chemical potentials

In the previous section, it was shown that the following chemical potentials are required H L to estimate rigorously L–H–V phase coexistence: µH m ( xm , T, P ), µw ( xm , T, P ), µw ( xm =

0, T, P), and µV m ( xm = 1, T, P ). In the next paragraph, we show that the two chemical potentials for pure phases, µLw ( xm = 0, T, P) and µV m ( xm = 1, T, P ), can be estimated in H a straightforward way. In contrast, µH m ( xm , T, P ) and µw ( xm , T, P ) will be estimated in a

second step using free energy calculations. L µV m ( xm = 1, T, P ) and µw ( xm = 0, T, P ). The chemical potential of methane in the vapor

phase µV m ( T, P ) was computed using its equation of state determined as follows. At a given T, isobaric-isothermal MD simulations are performed to determine the density of methane as a function of pressure, i.e. ρm ( T, P). In parallel, GCMC simulations are performed to determine the relation between the chemical potential and density of methane vapor, V ρ m ( µV m , T ). By inverting these two relationships, one obtains µm ( T, P ) as a function of T

and P (Table S1 in the Supporting Information displays µV m ( T, P ) for the various T and P

14


Page 15 of 50

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

Langmuir

considered in this work). The chemical potential of pure liquid water (µLw ( xm = 0, T, P)) at given T and P can be estimated using the Gibbs-Duhem equation: L NwL dµLw = −Sw dT + VwL dP

(9)

L , N L , and V L are the entropy, number of water molecules, and volume of the where Sw w w

liquid phase. If we assume that the density ρLw = NwL /VwL of liquid water is constant (incompressible liquid), integration of the Gibbs-Duhem equation at constant temperature T = T0 leads to: µLw ( xm = 0, T0 , P) = µLw ( xm = 0, T0 , P0 ) +

P − P0 ρw ( T0 , P0 )

(10)

It is convenient to take the L–V phase coexistence of water (T0 , P0 ) as a reference state since it is well-known for the different water models considered in this work. 66 In particular, for the temperature and pressure ranges considered here, water vapor along the L–V coexistence line can be treated as an ideal gas so that the chemical potential at coexistence L is readily obtained from the bulk saturating vapor pressure µV w ( T0 , P0 ) = µw ( T0 , P0 ) = √ k B T0 ln P0 Λ3 /k B T0 (Λ = h/ 2πmk B T is the thermal wavelength with h Plank’s constant

and m the molecular mass of water). Table S2 in the Supporting Information shows the chemical potential of water as a function of T and P (both the data for TIP4P/2005 and TIP4P/Ice are shown). H µH m ( xm , T, P ) and µw ( xm , T, P ). While the chemical potentials for pure phases (L and H V) are rather easy to assess, µH m ( xm , T, P ) and µw ( xm , T, P ) must be computed using a

more complex formalism which requires to combine SGMC simulations and free energy calculations. Let us consider a methane hydrate made up of Nm methane molecules and Nw water molecules at given T and P. For this system, an infinitely small change in the internal energy dU writes: 15


Langmuir

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 50

H dU = TdS − PdV + µH m dNm + µw dNw

(11)

where V and S are the volume and entropy of the methane hydrate, respectively. Legendre transformation of U with respect to S, V, Nw and Nm leads to: H U = TS − PV + µH m Nm + µw Nw

(12)

By comparing Eq. (11) with the derivative of Eq. (12), we obtain: H Nw dµH w = − SdT + VdP − Nm dµm

(13)

which is the Gibbs–Duhem equation for a binary mixture. Considering that Nw is constant in methane hydrate (owing to its crystalline structure), we can integrate Eq. (13) at constant T and P to obtain the change ∆µH w in the chemical potential for water between the zerooccupancy and occupied methane hydrate (i.e., as the methane mole fraction increases from 0 to xm ): ∆µH w

=

H µH w ( xm ) − µw ( xm

1 = 0) = − Nw

Z µH ( x m ) m

µH m ( x m =0)

Nm dµm

(14)

While Nm can be determined as a function of µH m using SGMC simulations as described in Section 2, the later equation shows that determining the chemical potential of water µH w in the hydrate phase requires to estimate the same chemical potential in the zero-occupancy H hydrate phase µH w ( xm = 0). The determination of µw ( xm = 0) is not straightforward and

requires free energy calculations that are reported in the next section.

16


Page 17 of 50

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

Langmuir

4.

PHASE DIAGRAM OF METHANE HYDRATE

4.1.

Free energy approach

Einstein molecule method. In Section 3, it was shown that determining the condition for L–H–V phase coexistence requires to estimate the chemical potentials for water in the liquid and hydrate phases and for methane in the vapor and hydrate phases: µLw ( xm = V H 0, T, P), µH w ( xm , T, P ), µm ( xm = 1, T, P ) and µm ( xm , T, P ). While the estimation of the

chemical potentials for the pure phases µLw ( xm = 0, T, P) and µV m ( xm = 1, T, P ) and for methane in the hydrate phase µH m ( xm , T, P ) does not raise important technical issues, the estimation of the the chemical potential for water in the hydrate phase µH w ( xm , T, P ) is not straightforward. However, as shown at the end of Section 3, µH w ( xm , T, P ) can be estimated from its value in the zero-occupancy hydrate µH w ( xm = 0, T, P ) (See Eq. (14)). By noting that the chemical potential is defined as the Gibbs free energy per water H molecule µH w ( xm = 0, T, P ) = Gw ( xm = 0) /Nw , the chemical potential of water in the zero-

occupancy methane hydrate can be estimated from the Helmholtz free energy AH w ( x m = 0):

µH w ( x m = 0) =

H ( x = 0) AH ( xm = 0) + PV Gw m = w Nw Nw

(15)

where the contribution PV is determined using molecular dynamics in the isobaricisothermal ensemble (NPT). In this section, we estimate AH w ( xm = 0) using free energy calculations based on the Einstein molecule approach developed by Vega and coworkers. 33,52 This technique, which derives from the Einstein crystal approach, consists of estimating AH w ( xm = 0) along a reversible thermodynamic path linking the real solid to an Einstein molecule; the Einstein molecule is an ideal crystalline structure without any intermolecular interactions in which each molecule is attached to its reference lattice position and orientation by a harmonic potential. The canonical partition function and free energy of this reference state are known analytically. For technical reasons, it is convenient to compute the partition function of the 17


Langmuir

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 50

Einstein molecule with one of its molecules at a fixed reference position (we note that the position of this reference molecule is constant but molecular rotation is allowed). Figure 2 shows the thermodynamic path used in the Einstein molecule approach to determine the free energy of the zero-occupancy methane hydrate. Throughout the manuscript, the superscript * indicates that the system has one of its water molecules at a fixed position (this molecule is shown by the big pink ’+’ sign in Figure 2). The reversible integration path considered in the Einstein molecule approach consists of 4 steps which transform the ideal Einstein molecule into the zero-occupancy methane hydrate: 1. We start from the non-interacting Einstein molecule (A) whose free energy AA is known analytically; AA = −k B T ln QA where QA is the canonical partition function of the non-interacting Einstein molecule. The first step in the Einstein molecule approach consists of fixing the position of one of its water molecules to form a constrained, non–interacting Einstein molecule (A*). The free energy change corre sponding to this transformation is simply ∆AA→A* = AA* − AA = k B T ln V/Λ3

where V is the volume of the Einstein molecule and Λ the thermal wavelength of the water molecule;

2. The constrained, non-interacting Einstein molecule (A*) is transformed into the corresponding interacting Einstein molecule (B*) by adding the intermolecular potential energy between water molecules (which includes the Lennard-Jones and Coulomb potentials as described in Section 2). In this step, both the non-interacting and interacting Einstein molecules have one of their water molecules at a fixed position so that both of these structures are referred to as ”constrained”. The free energy difference along this step, ∆A1 = AB* − AA* , is determined using a perturbation treatment described below. 3. The constrained interacting Einstein molecule (B*) is transformed into the corresponding constrained, zero-occupancy methane hydrate (C*) by gradually switching 18


Page 19 of 50

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

Langmuir

off the harmonic potentials UA that attach the water molecules to their reference lattice position in the Einstein molecule. The free energy difference in this step, ∆A2 = AC* − AB* , is determined by thermodynamic integration as described below; 4. The zero-occupancy methane hydrate (C) is obtained from the constrained, zerooccupancy methane hydrate (C*) by releasing the constraint over the fixed water molecule. The free energy change for this step simply writes ∆AC*→C = AC − AC* = −k B T ln V/Λ3 .

The thermodynamic path above allows writing the free energy of the zero-occupancy methane hydrate as AC = AA + ( AA* − AA ) + ( AB* − AA* ) + ( AC* − AB* ) + ( AC − AC* )

= AA + k B T ln

V V + ∆A1 + ∆A2 − k B T ln 3 = AA + ∆A1 + ∆A2 3 Λ Λ

(16)

where we used that constraining (step 1) and unconstraining (step 4) the position of one reference water molecule in the thermodynamic path cancel out. While these free energy calculations should not depend on a specific choice for the Einstein molecule (provided a reasonable configuration is used), we followed here the annealing approach suggested by Noya and coworkers. 50 First, the Einstein molecule is selected with a volume identical to that of real methane hydrate as obtained using isobaric–isothermal MD simulations at P = 1, 10, and 100 atm. Then, a simulated annealing strategy (canonical ensemble) is used to determine the final configuration; the temperature is decreased from T = 180 K to 1 K with temperature steps of 10 K. Eq. (16) shows that only the three following contributions must be calculated to determine the free energy of the zero-occupancy methane hydrate: AA , ∆A1 and ∆A2 . In the rest of this subsection, we determine these three contributions before gathering all the data to estimate the free energy of the zero-occupancy methane H H hydrate AC = AH w ( xm = 0) and the chemical potentials µw and µm in the real (i.e. methane

occupied) methane hydrate. 19


Langmuir


Page 20 of 50

Page 21 of 50

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

Langmuir

Free Energy AA of the non-interacting Einstein molecule. Water molecules in the non-interacting Einstein molecule (A) are attached to their reference lattice position and orientation through harmonic potentials so that its potential energy writes: N

UA (Ri , φa,i , φb,i ) = λ T

N

(0)

∑ (Ri − Ri )2 + λR ∑ [(sin2 φa,i + (

i =1

i =1

φb,i 2 )] π

(17)

where the sum runs over each molecule i of the N molecules in the system. The first term in Eq. (17) corresponds to harmonic potentials acting on each molecule position Ri (0)

with an equilibrium position defined as the reference position Ri . Similarly, the second term in Eq. (17) corresponds to harmonic potentials acting on each molecule orientation (0)

defined by two vectors a and b with equilibrium vectors ai

(0)

and bi

corresponding to the

reference molecule orientation. As shown in Figure S3 of the Supporting Information, the two orientation vectors can be chosen as a = (l1 − l2 )/|l1 − l2 | and b = (l1 + l2 )/|l1 + l2 | where l1 and l2 are the vectors along the O–H bonds in the water molecule. For each water (0)

(0)

molecule i, φa,i = cos(ai ·ai ) and φb,i = cos(bi ·bi ). Following previous works, 21,22,49,52 2

˚ = λ T /k B T = 25000 the spring constants in UA (Ri , φa,i , φb,i ) were selected as λ R /k B T A (note that when reasonable choices are made for these parameters, AA is independent of these values as harmonic oscillators only depend on temperature). The Helmholtz free energy AA of the non-interacting Einstein molecule, which can be computed from its canonical partition function QA , subdivides into a translation AA,T and a rotation AA,R contributions: A A ln QA AA =− = A,T + A,R Nk B T N Nk B T Nk B T

(18)

where all free energy contributions are normalized to the total thermal energy Nk B T. As shown in Section S4 of the Supporting Information, these two contributions can be

21


Langmuir

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 22 of 50

expressed as: 2 AA,T Λ λT 3 1 NΛ3 1 + 1− ln = ln Nk B T N V 2 N k B Tπ 2 Λ λT 3 ln ∼ N →∞ 2 k B Tπ AA,R = − ln Nk B T

Z

1 0

Z 1 λR λR 2 2 (1 − x ) dx y dy exp − exp − kB T kB T 0

(19)

(20)

Calculations based on these expressions, including numerical integration of Eq. (20), can be found in Section S4 of the Supporting Information and lead to A A,T /( Nk B T ) = 29.43, A E,R /( Nk B T ) = 16.01. These values are fully consistent with those reported by Vega and coworkers for hexagonal ice.

52

Free energy difference ∆A1 . The free energy change ∆A1 = AB* − AA* between the non-interacting and interacting Einstein molecules is estimated through a perturbation approach. We write that the potential energy in the interacting Einstein molecule UB∗ is the sum of the non-interacting Einstein molecule U A∗ and the intermolecular potential energy U, i.e. UB∗ = U A∗ + U. For large λ R and λ T , U Tm . In general, such direct coexistence simulations are conducted in the isobaric-isothermal ensemble (NPT) because phase transitions occurs at constant T and P. As a result, all direct coexistence method strategies applied to methane hydrate have been carried out so far in this ensemble. 6,19,85 However, for binary compounds such as methane hydrate, such coexistence simulations can be performed in the Grand Canonical ensemble where the system volume V, temperature T, and chemical potentials for water µw and methane µm are constant. In the present work, we adopted this strategy which has not been considered previously to the best of our knowledge. Considering such an open ensemble in which the numbers of water and methane molecules fluctuate present several advantages over constant number of molecules ensemble (such as NVT or NPT simulations). First, this allows considering small system sizes since the number of methane molecules will adjust upon methane hydrate formation even though the initial number of methane molecules is small. In contrast, with constant N simulations, one has to simulate a large domain of methane molecules that acts as a methane source to fill the water cages upon methane 31


Langmuir

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 32 of 50

hydrate formation. Moreover, by considering an ensemble where density will change through molecule numbers fluctuations, one avoids difficulties due to inefficient/limited sampling in volume changes. Finally, in GCMC simulations, molecule insertion/deletions are attempted randomly, homogeneously throughout the simulation box so that difficulties inherent to slow diffusion between the methane hydrate and liquid/fluid phases are overcome. For such complex systems, DCM should be used with caution because of the initial coexisting system can be chosen in different ways. According to Gibbs’ phase rule, in the temperature/pressure range where methane hydrate is stable, it coexists with the liquid (water-rich) and vapor (methane-rich) phases. As a result, initial phase coexistence in DCM can be chosen as a system made of two of these three phases or three phases. In the present work, we chose to consider phase coexistence between the liquid phase and methane hydrate; while this corresponds to an approximation, the use of the Grand Canonical ensemble ensures that three-phase coexistence is simulated in fact; because the system is in equilibrium with an infinite reservoir of bulk molecules at chemical potentials corresponding to those of the water-rich liquid and methane-rich vapor, DCM simulations in this specific ensemble are equivalent to simulating a system with threephase coexistence. In order to prepare the initial system (i.e. methane hydrate coexisting with liquid water), several strategies are possible. In this work, we started from a methane hydrate phase having the following dimensions: L x = Ly = Lz = 2.3754 nm. Periodic boundary conditions were applied in each direction to avoid finite size effects. We started from a hydrate phase equilibrated at low T (we recall that the pressure was set to 100 atm). Then, molecules located in the region z < 0 were frozen while the rest of the simulation box was equilibrated at high temperature T to melt the hydrate located in the region z > 0. In so doing, one obtains a coexisting system made of methane hydrate in equilibrium with the liquid phase (Figure 7(a)). Obviously, this system is maintained at coexistence condition in an unphysical fashion and, depending on the temperature used in subsequent

32


Page 33 of 50

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

Langmuir

GCMC simulations, the system will melt or form hydrate (unless in the very unlikely event that the chosen temperature and chemical potentials exactly correspond to phase coexistence). From a practical point of view, for P = 100 atm, we performed M = 18 simulations with temperatures in the range T = 180–350 K (the temperature interval is 10 K). Our DCM simulations in the Grand Canonical ensemble at a given pressure and temperature require to specify chemical potentials for water and methane. In the present work, as described in Section 3.2., the chemical potential for water in the liquid phase was chosen equal to that of pure liquid water while the chemical potential for methane in the vapor phase was chosen equal to that or pure methane vapor. Figure 7(b) shows the methane xm and water xw mole fractions in the system in the course of the GCMC simulation (i.e. the number of MC moves performed with one MC move corresponding to a molecule translation, rotation, insertion or creation). Results for different temperatures are shown: T = 260, 270, 280, 290, and 300 K. On the one hand, at high temperature, T ≥ 290 K, the system melts as evidenced by the decrease in the methane mole fraction xm . As expected, xm ∼ 0 (xw ∼ 1) in the liquid phase, which further justifies our choice in the L–H–V equilibrium condition to assume that µLw ( xm , T, P) ∼ µLw ( xm = 0, T, P). On the other hand, at low temperature, T ≤ 280 K, the methane mole fraction increases (while xw decreases) upon methane hydrate formation. While melting does not suffer from ambiguity since all methane hydrate is transformed into liquid, it should be emphasized that hydrate formation was found to be inefficient; due to the low probability to nucleate hydrate cages (inherent to their very small entropy), it was observed that formation of the hydrate is incomplete. As a result, despite the coexistence with an already formed hydrate, many long GCMC runs (about 7-8 ×108 MC moves for a system size of the order of ∼ 102 − 103 molecules) were not sufficient to lead to perfect methane hydrates. Despite this drawback of the direct coexistence method, the results above show that the equilibrium temperature for hydrate/liquid coexistence is comprised between 280 K and 290 K. As can be seen in Figure 6, this coexistence temperature is in

33


Langmuir


Page 34 of 50

Page 35 of 50

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

Langmuir

For each replica, conventional GCMC moves are performed: molecule translation, rotation, deletion and insertion. Moreover, trial swap moves between configuration α (energy U α , α methane molecules) in replica (1) and configuration β (energy Nwα water molecules and Nm β

β

U β , Nw water molecules and Nm methane molecules) in replica (2) are attempted. The swap move is accepted or rejected according to the following Metropolis probability:

  T 3( Nmβ + Nwβ − Nmα − Nmα )/2 2 Pacc (α1 , β 2 → α2 , β 1 ) = min 1,  T1 1 1 β α U −U − exp k B T2 k B T1 " ! #) µ2i µ1i β α ∏ exp k B T1 − k B T2 Ni − Ni i =m,w

(26)

In this work, the different replicas were considered at temperatures and chemical potentials corresponding to a pressure P = 100 atm. The temperature of the different replicas ranges from 283 to 298 K with a temperature difference between two successive replicas of ∆T = 1 K. In theory, hyper parallel tempering should provide a rigorous description of methane hydrate formation/dissociation as a function of temperature provided that both configurations corresponding to the liquid phase and the methane hydrate phase are considered in the initial replicas; for long enough simulations, swapping between the liquid and solid phases at different temperatures should lead to an accurate estimate of the phase transition temperature Tm with liquid configurations for T > Tm and methane hydrate configurations for T > Tm . However, in practice, very low swapping probabilities were observed between liquid and methane hydrate configurations due to the large differences in water and methane molecule numbers in these two states (as can be seen in the acceptance probability in the equation above, the difference in the number of molecules is an important parameter). We found that this issue can be overcome by considering in the initial replicas composite configurations corresponding

35


Langmuir

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 36 of 50

to mixtures of the liquid and hydrate phases (in the spirit of the mixture considered as the initial configuration for the direct coexistence method). As shown in Figure 8(a), in addition to pure liquid and hydrate configurations, several configurations corresponding to methane hydrate regions coexisting with the liquid phase were considered (these mixtures correspond to different hydrate volume fractions ranging from 0.25 to 0.75). The total number of methane and water molecules in each replica is of the order of ∼ 102 − 103 . Equilibration was reached after 9 × 108 Monte Carlo steps and water and methane mole fractions were averaged over another 1 × 108 Monte Carlo steps. Figure 8(b) shows the methane xm and water xw mole fractions as a function of temperature T once equilibrium has been reached. The sharp decrease (increase) at Tm = 289.5 K in xm (xw ) indicates melting of the methane hydrate. Such a transition temperature for P = 100 atm is consistent with the values obtained using free energy calculations and the direct coexistence method. These results show that such a hyper parallel tempering technique improves the sampling of phase space and allows determining accurately the melting temperature of complex, non stoichiometric systems such as methane hydrates (by preventing the system from being trapped in local metastable states).

5.

CONCLUSION

Using different molecular simulation strategies, we determined the pressure–temperature phase diagram for bulk methane hydrate. For two different water models, we first determined the liquid–hydrate–vapor phase coexistence using rigorous free energy calculations based on the Einstein molecule approach. Our data, which are consistent with previous molecular simulation works, shows that the different thermodynamic approximations such as the description of methane vapor are important. Overall, in agreement with previous studies, it is shown that the choice of the water model is a key problem and that TIP4P/Ice, which was specifically developed to reproduce crystalline phases of water, reproduces

36


Page 37 of 50

Langmuir


Langmuir

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 38 of 50

molecules upon hydrate dissociation/formation. This allows considering calculations with the direct coexistence method using system sizes that remains small (otherwise, large methane regions in the system have to be considered to act as methane molecules source/sink upon melting/crystallization of the hydrate). In addition to the direct coexistence method extended to the grand canonical ensemble, we also considered hyper parallel tempering which consists of considering several replicas of the system at different temperatures and chemical potentials – the system being therefore treated in the grand canonical ensemble to allow for large changes in its composition upon hydrate formation/dissociation. Despite the reduced accuracy/robustness compared to more rigorous approaches based on free energy techniques, both the direct coexistence method and hyper parallel tempering technique were found to lead to reasonable predictions for phase coexistence. However, while the results reported in this work shows that these two direct techniques can be used to estimate stability conditions for methane hydrate, we emphasize that several refinements and ”tricks” were needed to lead to sufficient sampling of the phase space and accurate phase coexistence predictions. First, as mentioned above, both the direct coexistence and hyper parallel techniques were used with water and methane treated in the Grand Canonical ensemble; we found that this was needed to efficiently sample large molecule number fluctuations upon hydrate formation/dissociation. Moreover, in the case of hyper parallel tempering, we also found that the initial replicas (i.e. at different temperatures and chemical potentials) must include composite systems where both the hydrate and liquid phases coexist. Such coexisting states allow sufficient swapping along the hyper parallel tempering simulation between the low and high temperature replicas. Otherwise, considering the Metropolis acceptance probability in this hyper grand canonical ensemble given in Eq. (27), the large difference in the numbers of water and methane molecules between the liquid and hydrate phases lead to very low swapping probabilities (too low to allow efficient sampling). As a result, while our data show that accurate hydrate

38


Page 39 of 50

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

Langmuir

stability conditions can be estimated in principle using hyper parallel tempering, the latter drawback constitutes an important limitation to this technique. Finally, in addition to being more robust than the direct coexistence and hyper parallel tempering methods, free energy calculations provide accurate estimates for the chemical potentials for water and methane in the hydrate phase, including their values at phase coexistence (in contrast, with the two direct techniques, one has to estimate in an approximate fashion the chemical potentials that lead to phase equilibrium). This is a key asset of the free energy technique over direct methods since such chemical potentials at phase coexistence will be used in subsequent work on the stability of methane hydrate confined in porous media (which are in equilibrium with an external methane and water mixture or hydrate imposing its chemical potentials at constant temperature).

Acknowledgement We are pleased to contribute with this paper to the special issue of Langmuir honoring Professor Keith Gubbins on the occasion of his 80th birthday. Among his numerous works, Keith Gubbins has made significant contributions to the field of statistical mechanics and molecular simulation of liquids with important original studies on the phase diagram of confined liquids and solids. D. Jin thanks the China Scholarship Council for financial support (No. 201506450015). The authors wish to thank F. Calvo for very interesting discussions. We are also grateful to J. L. Barrat, J. M. Herri, and L. Gelb for helpful comments.

Supporting Information Available The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.xxxxxxx. Details about the procedure to generate methane

39


Langmuir

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 40 of 50

hydrate configurations, free energy calculations and additional data such as chemical potentials of methane and water as a function of the number of methane in the hydrate. An xyz configuration of the 2 × 2 × 2 unit cells of methane hydrate generated according to the strategy described in our paper is also available (the box dimensions are L x = Ly = Lz = 2.3754 nm). In this configuration, methane is described as a united atom and water has a O–H bond length of 0.1 nm and a bond angle of ∼109◦ (note that this configuration must be relaxed with specific bond constraints to reach a water molecule conformation corresponding to a given water model such as TIP4P, SPC, etc.).

References (1) Davy, H. The Bakerian Lecture: On Some of the Combinations of Oxymuriatic Gas and Oxygene, and on the Chemical Relations of These Principles, to Inflammable Bodies. Philos. Trans. R. Soc. London, Ser. B 1800, 1, 385–388. (2) Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; Chemical Industries; CRC Press, 2007. (3) McMullan, R. K.; Jeffrey, G. A. Polyhedral Clathrate Hydrates. IX. Structure of Ethylene Oxide Hydrate. J. Chem. Phys. 1965, 42, 2725. (4) Mak, T. C. W.; McMullan, R. K. Polyhedral Clathrate Hydrates. X. Structure of the Double Hydrate of Tetrahydrofuran and Hydrogen Sulfide. J. Chem. Phys. 1965, 42, 2732–2737. (5) Ripmeester, J. A.; Tse, J. S.; Ratcliffe, C. I.; Powell, B. M. A New Clathrate Hydrate Structure. Nature (London) 1987, 325, 135. (6) Michalis, V. K.; Costandy, J.; Tsimpanogiannis, I. N.; Stubos, A. K.; Economou, I. G. Prediction of the Phase Equilibria of Methane Hydrates Using the Direct Phase Coexistence Methodology. J. Chem. Phys. 2015, 142, 044501. 40


Page 41 of 50

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

Langmuir

(7) Kvenvolden, K. A. Methane hydrate – A Major Reservoir of Carbon in the Shallow Geosphere? Chem. Geol. 1988, 71, 41. (8) MacDonald, G. J. The Future of Methane as an Energy Resource. Annu. Rev. Energy 1990, 15, 53. (9) Udachin, K. A.; Enright, G. D.; Ratcliffe, C. I.; Ripmeester, J. A. Structure, Stoichiometry, and Morphology of Bromine Hydrate. J. Am. Chem. Soc. 1997, 119, 11481–11486. (10) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N.; Sloan, E. D. Stable Low-Pressure Hydrogen Clusters Stored in a Binary Clathrate Hydrate. Science 2004, 306, 469. ¨ (11) Schuth, F. Technology: Hydrogen and Hydrates. Nature (London) 2005, 434, 712. (12) Lee, H.; Lee, J.; Kim, D. Y.; Park, J.; Seo, Y.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A. Tuning Clathrate Hydrates for Hydrogen Storage. Nature (London) 2005, 434, 743. (13) Strobel, T. A.; Sloan, E. D.; Koh, C. A. Raman Spectroscopic Studies of Hydrogen Clathrate Hydrates. J. Chem. Phys. 2009, 130, 014506. (14) Henriet, J. P.; Mienert, J. Gas Hydrates: Relevance to World Margin Stability and Climatic Change; Geological Society of London, 1998. (15) Kieffer, S. W.; Lu, X. L.; Bethke, C. M.; Spencer, J. R.; Marshak, S.; Navrotsky, A. A Clathrate Reservoir Hypothesis for Enceladus’ South Polar Plume. Science 2006, 314, 1764–1766. (16) Hersant, F.; Gautier, D.; Lunine, J. I. Enrichment in Volatiles in the Giant Planets of the Solar System. Planet. Space Sci. 2004, 52, 623–641.

41


Langmuir

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 42 of 50

(17) Coasne, B.; Czwartos, J.; Gubbins, K. E.; Hung, F. R.; Sliwinska-Bartkowiak, M. Freezing and Melting of Binary Mixtures Confined in a Nanopore. Mol. Phys. 2004, 102, 2149–2163. (18) Czwartos, J.; Coasne, B.; Gubbins, K. E.; Hung, F. R.; Sliwinska-Bartkowiak, M. Freezing and Melting of Azeotropic Mixtures Confined in Nanopores: Experiment and Molecular Simulation. Mol. Phys. 2005, 103, 3103–3113. (19) Conde, M. M.; Vega, C. Determining the Three-phase Coexistence Line in Methane Hydrates Using Computer Dimulations. J. Chem. Phys. 2010, 133, 064507. (20) Docherty, H.; Galindo, A.; Vega, C.; Sanz, E. A Potential Model for Methane in Water Describing Correctly the Solubility of the Gas and the Properties of the Methane Hydrate. J. Chem. Phys. 2006, 125, 074510. (21) Jensen, L.; Thomsen, K.; von Solms, N.; Wierzchowski, S.; Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Calculation of Liquid Water–Hydrate–Methane Vapor Phase Equilibria from Molecular Simulations. J. Phys. Chem. B 2010, 114, 5775– 5782. (22) Wierzchowski, S. J.; Monson, P. A. Calculation of Free Energies and Chemical Potentials for Gas Hydrates Using Monte Carlo Simulations. J. Phys. Chem. B 2007, 111, 7274–7282. (23) English, N. J.; Gorman, P. D.; MacElroy, J. M. D. Mechanisms for Thermal Conduction in Hydrogen Hydrate. J. Chem. Phys. 2012, 136, 044501. (24) Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 2010, 132, 11806–11811, PMID: 20669949. (25) Knott, B. C.; Molinero, V.; Doherty, M. F.; Peters, B. Homogeneous Nucleation of

42


Page 43 of 50

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

Langmuir

Methane Hydrates: Unrealistic under Realistic Conditions. J. Am. Chem. Soc. 2012, 134, 19544–19547, PMID: 23148735. (26) Nguyen, A. H.; Molinero, V. Identification of Clathrate Hydrates, Hexagonal Ice, Cubic Ice, and Liquid Water in Simulations: the CHILL+ Algorithm. J. Phys. Chem. B 2015, 119, 9369–9376, PMID: 25389702. (27) Pefoute, E.; Prager, M.; Russina, M.; Desmedt, A. Quasi-elastic Neutron Scattering Investigation of the Guest Molecule Dynamics in the Bromomethane Clathrate Hydrate. Fluid Phase Equilibria 2016, 413, 116 – 122, Special Issue: Gas Hydrates and Semiclathrate Hydrates. (28) Desmedt, A.; Bedouret, L.; Pefoute, E.; Pouvreau, M.; Say-Liang-Fat, S.; Alvarez, M. Energy Landscape of Clathrate Hydrates. Eur. Phys. J. Spec. Top. 2012, 213, 103–127. (29) Sloan, E. D. Fundamental Principles and Applications of Natural Gas Hydrates. Nature 2003, 426, 353–363. (30) Huo, Z.; Hester, K.; Sloan, E. D.; Miller, K. T. Methane Hydrate Nonstoichiometry and Phase Diagram. AIChE J. 2003, 49, 1300–1306. (31) Wu, R.; Kozielski, K. A.; Hartley, P. G.; May, E. F.; Boxall, J.; Maeda, N. Probability Distributions of Gas Hydrate Formation. AIChE J. 2013, 59, 2640–2646. (32) Ueno, H.; Akiba, H.; Akatsu, S.; Ohmura, R. Crystal Growth of Clathrate Hydrates Formed with Methane + Carbon Dioxide Mixed Gas at the Gas/Liquid Interface and in Liquid Water. New J. Chem. 2015, 39, 8254–8262. (33) Conde, M. M.; Torr, J. P.; Miqueu, C. Revisiting the Thermodynamic Modelling of Type I Gas–Hydroquinone Clathrates. Phys. Chem. Chem. Phys. 2016, 18, 10018–10027. (34) Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V. Freezing, Melting

43


Langmuir

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 44 of 50

and Structure of Ice in a Hydrophilic Nanopore. Phys. Chem. Chem. Phys. 2010, 12, 4124–4134. (35) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bond-orientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, 784–805. (36) Radhakrishnan, R.; Trout, B. L. Nucleation of Crystalline Phases of Water in Homogeneous and Inhomogeneous Environments. Phys. Rev. Lett. 2003, 90, 158301. (37) Radhakrishnan, R.; Trout, B. L. Nucleation of Hexagonal Ice in Liquid Water. J. Am. Chem. Soc. 2003, 125, 7743–7747. (38) Errington, J. R.; Debenedetti, P. G. Relationship between Structural Order and the Anomalies of Liquid Water. Nature 2001, 409, 318–321. (39) Radhakrishnan, R.; Trout, B. L. A New Approach for Studying Nucleation Phenomena Using Molecular Simulations: Application to CO2 Hydrate Clathrates. J. Chem. Phys. 2002, 117, 1786. (40) Moon, C.; Taylor, P. C.; Rodger, P. M. Molecular Dynamics Study of Gas Hydrate Formation. J. Am. Chem. Soc. 2003, 125, 4706–4707. (41) Chakraborty, S. N.; Grzelak, E. M.; Barnes, B. C.; Wu, D. T.; Sum, A. K. Voronoi Tessellation Analysis of Clathrate Hydrates. J. Phys. Chem. C 2012, 116, 20040–20046. (42) Barnes, B. C.; Beckham, G. T.; Wu, D. T.; Sum, A. K. Two-component Order Parameter for Quantifying Clathrate Hydrate Nucleation and Growth. J. Chem. Phys. 2014, 140, 164506. (43) Jacobson, L. C.; Matsumoto, M.; Molinero, V. Order Parameters for the Multistep Crystallization of Clathrate Hydrates. J. Chem. Phys. 2011, 135, 074501. (44) Moon, C.; Hawtin, R. W.; Rodger, P. M. Nucleation and Control of Clathrate Hydrates: Insights from Simulation. Faraday Discuss. 2007, 136, 367–382. 44


Page 45 of 50

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

Langmuir

(45) 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–4864. (46) Chakraborty, S. N.; Gelb, L. D. A Monte Carlo Simulation Study of Methane Clathrate Hydrates Confined in Slit-Shaped Pores. J. Phys. Chem. B 2012, 116, 2183–2197. (47) Gasser, U.; Weeks, E. R.; Schofield, A.; Pusey, P. N.; Weitz, D. A. Real-Space Imaging of Nucleation and Growth in Colloidal Crystallization. Science 2001, 292, 258–262. (48) van der Waals, J. H.; Platteeuw, J. C. Clathrate Solutions. Adv. Chem. Phys. 1959, 2, 1–57. (49) Vega, C.; Noya, E. G. Revisiting the Frenkel-Ladd Method to Compute the Free Energy of Solids: The Einstein Molecule Approach. J. Chem. Phys. 2007, 127, 154113. (50) Noya, E. G.; Conde, M. M.; Vega, C. Computing the Free Energy of Molecular Solids by the Einstein Molecule Approach: Ices XIII and XIV, Hard-dumbbells and a Patchy Model Proteins. J. Chem. Phys. 2008, 129, 104704. (51) Frenkel, D.; Ladd, A. J. C. New Monte Carlo Method to Compute the Free Energy of Arbitrary Solids. Application to the FCC and HCP Phases of Hard Spheres. J. Chem. Phys. 1984, 81, 3188–3193. (52) Vega, C.; Sanz, E.; Abascal, J. L. F.; Noya, E. G. Determination of Phase Diagrams via Computer Simulation: Methodology and Applications to Water, Electrolytes and Proteins. J. Phys.: Condens. Matter 2008, 20, 153101. (53) Conde, M. M.; Gonzalez, M. A.; Abascal, J. L. F.; Vega, C. Determining the Phase Diagram of Water from Direct Coexistence Simulations: The Phase Diagram of the TIP4P/2005 Model Revisited. J. Chem. Phys. 2013, 139, 154505. (54) Yan, Q.; de Pablo, J. J. Hyper-parallel Tempering Monte Carlo: Application to the

45


Langmuir

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 46 of 50

Lennard-Jones Fluid and the Restricted Primitive Model. J. Chem. Phys. 1999, 111, 9509–9516. (55) de Pablo, J. J.; Laso, M.; Suter, U. W. Simulation of Polyethylene Above and Below the Melting Point. J. Chem. Phys. 1992, 96, 2395–2403. (56) Yan, Q.; de Pablo, J. J. Hyperparallel Tempering Monte Carlo Simulation of Polymeric Systems. J. Chem. Phys. 2000, 113, 1276–1282. (57) Coasne, B. Freezing of Mixtures Confined in a Slit Nanopore. Adsorption 2005, 11, 301–306. (58) Waage, M. H.; Vlugt, T. J. H.; Kjelstrup, S. Phase Diagram of Methane and Carbon Dioxide Hydrates Computed by Monte Carlo Simulations. J. Phys. Chem. B 2017, 121, 7336–7350. (59) Ohgaki, K.; Hamanaka, T. Phase-Behavior of CO2 Hydrate–Liquid CO2 –H2 O System at High Pressure. Kagaku Kogaku Ronbunshu 1995, 21, 800–803. (60) Nakano, S.; Moritoki, M.; Ohgaki, K. High-Pressure Phase Equilibrium and Raman Microprobe Spectroscopic Studies on the CO2 Hydrate System. J. Chem. Eng. Data 1998, 43, 807–810. (61) Ravipati, S.; Punnathanam, S. N. Calculation of three-phase methaneethane binary clathrate hydrate phase equilibrium from Monte Carlo molecular simulations. Fluid Phase Equilibria 2014, 376, 193–201. (62) Costandy, J.; Michalis, V. K.; Tsimpanogiannis, I. N.; Stubos, A. K.; Economou, I. G. The role of intermolecular interactions in the prediction of the phase equilibria of carbon dioxide hydrates. J. Chem. Phys. 2015, 143, 094506. (63) Mguez, J. M.; Conde, M. M.; Torr, J.-P.; Blas, F. J.; Pieiro, M. M.; Vega, C. Molecular

46


Page 47 of 50

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

Langmuir

dynamics simulation of CO2 hydrates: Prediction of three phase coexistence line. J. Chem. Phys. 2015, 142, 124505. (64) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. (65) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638–6646. (66) Vega, C.; Abascal, J. L. F.; Nezbeda, I. Vapor–Liquid Equilibria from the Triple Point up to the Critical Point for the New Generation of TIP4P-like models: TIP4P/Ew, TIP4P/2005, and TIP4P/ICE. J. Chem. Phys. 2006, 125, 034503. (67) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (68) Abascal, J. L. F.; Sanz, E.; Fernndez, R. G.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/ICE. J. Chem. Phys. 2005, 122, 234511. (69) Aragones, J. L.; Conde, M. M.; Noya, E. G.; Vega, C. The Phase Diagram of Water at High Pressures as Obtained by Computer Simulations of the TIP4P/2005 Model: the Appearance of a Plastic Crystal Phase. Phys. Chem. Chem. Phys. 2009, 11, 543–555. (70) Berthelot, D. Sur le mlange des gaz. Comptes rendus hebdomadaires des sances de l’Acadmie des Sciences 1898, 126, 1703–1855. (71) Lorentz, H. A. Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Annalen der Physik 1881, 248, 127–136. (72) Kolafa, J.; Perram, J. W. Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems. Mol. Simul. 1992, 9, 351–368.

47


Langmuir

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 48 of 50

(73) Di Pierro, M.; Elber, R.; Leimkuhler, B. A Stochastic Algorithm for the IsobaricIsothermal Ensemble with Ewald Summations for All Long Range Forces. J. Chem. Theory Comput. 2015, 11, 5624–5637. (74) Ewald, P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik 1921, 369, 253–287. (75) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Computational Science; Academic Press, 2002. (76) Buch, V.; Sandler, P.; Sadlej, J. Simulations of H2 O Solid, Liquid, and Clusters, with an Emphasis on Ferroelectric Ordering Transition in Hexagonal Ice. J. Phys. Chem. 1998, 102, 8641–8653. (77) Kirchner, M. T.; Boese, R.; Billups, W. E.; Norman, L. R. Gas Hydrate Single-Crystal Structure Analyses. J. Am. Chem. Soc. 2004, 126, 9407–9412. (78) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (79) Verlet, L. Computer Experiments on Classical Fluids: I. Thermodynamical Properties of Lennard–Jones Molecules. Phys. Rev. 1967, 159, 98–103. (80) Hoover, W. G. Canonical Dynamics: Equilibrium Phase–Space Distributions. Phys. Rev. A 1985, 31, 1695–1697. (81) Nos´e, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. (82) Harvey, A. H. Semiempirical Correlation for Henry’s Constants over Large Temperature Ranges. AIChE J. 1996, 42, 1491–1494. (83) Harvey, A. H.; Sengers, J. M. H. L. Correlation of Aqueous Henry’s Constants from 0◦ C to the Critical Point. AIChE J. 1990, 36, 539–546. 48


Page 49 of 50

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

Langmuir

(84) Nagle, J. F. Lattice Statistics of Hydrogen Bonded Crystals. I. The Residual Entropy of Ice. J. of Math. Phys. 1966, 7, 1484–1491. (85) Tung, Y.-T.; Chen, L.-J.; Chen, Y.-P.; Lin, S.-T. The Growth of Structure I Methane Hydrate from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 10804– 10813, PMID: 20669917.

49


Langmuir


Page 50 of 50

Molecular Simulation of the Phase Diagram of Methane Hydrate: Free

Recommend Documents