Mechanism of Slow Crystal Growth of Tetrahydrofuran Clathrate

---

Subscriber access provided by GAZI UNIV

Article

Mechanism of Slow Crystal Growth of Tetrahydrofuran Clathrate Hydrate Takuma Yagasaki, Masakazu Matsumoto, and Hideki Tanaka J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b10293 • Publication Date (Web): 01 Feb 2016 Downloaded from http://pubs.acs.org on February 2, 2016

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 C 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 34

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

Mechanism of Slow Crystal Growth of Tetrahydrofuran Clathrate Hydrate

Takuma Yagasaki,† Masakazu Matsumoto,† and Hideki Tanaka†, ‡ †

Department of Chemistry, Faculty of Science, Okayama University, Okayama, 700-8530,

Japan ‡

Research Center of New Functional Materials for Energy Production, Storage and

Transport, Okayama, 700-8530, Japan

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 2 of 34

Abstract Tetrahydrofuran (THF) clathrate hydrate has been frequently used in experimental studies instead of gas hydrates because it forms at a temperature higher than the ice point under ambient pressure.

In this paper, we compare the crystal growth rates of THF

hydrate and ice using molecular dynamics simulations.

It is demonstrated that the crystal

growth of THF hydrate is much slower than that of ice.

The growth rates of THF hydrate

significantly deviate from a standard kinetic model known as the Wilson-Frenkel model whereas it reproduces the temperature dependence of the growth rate of ice.

The slow

crystal growth and the deviation from the Wilson-Frenkel model are attributed to the trapping of THF molecules in open small cages at the hydrate surface.

We calculate the

free energy profile of a THF molecule transferring from the bulk solution phase to the hydrate surface using the umbrella sampling technique.

It is shown that a THF molecule

trapped in an open small cage needs to cross one or two free energy barriers to escape from the surface region.

We also refer to the similarity between the mechanism of slow growth

of THF hydrate and the effect of kinetic hydrate inhibitors.

2 ACS Paragon Plus Environment

Page 3 of 34

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

Introduction Clathrate hydrates are crystalline inclusion compounds in which water molecules form polyhedral cavities.1

The hydrate cavities can incorporate various small nonpolar

molecules such as methane and hydrogen. clathrate hydrate.

Some weakly polar molecules also form

A typical example is tetrahydrofuran (THF).

THF is a liquid which is

miscible with water in any proportion, and forms structure II hydrate.

A unit cell of

structure II hydrates consists of sixteen small (512) cages and eight large (51264) cages. THF molecule occupies only a large cage because of its large molecular size.

A

Unlike

usual gas hydrates, THF hydrate forms at a temperature higher than the ice point under ambient pressure.

This property allows one to use a simpler experimental setup than the

case of gas hydrates which demands high pressure vessels, and therefore THF hydrate has been extensively investigated as an analogue to gas hydrates.2-5

For example, low dosage

hydrate inhibitors (LDHI), which are used to avoid plugging of oil and gas pipelines caused by clathrate hydrates, are often tested with THF hydrate in laboratory before field trials.6,7 Another important property is that the presence of THF significantly reduces the formation pressure of clathrate hydrates of various gas molecules.8,9

This is quite useful for the

purposes of gas storage and separation using clathrate hydrates.10-16 Molecular dynamics (MD) simulations have provided a wealth of information on the microscopic mechanism of hydrate formation.17-49

It was reported that the growth rate of

methane hydrate can be faster than that of ice when methane is highly supersaturated in the aqueous phase.20

It is also known that nucleation of gas hydrates in supersaturated

aqueous solutions of gases can be observed in all-atom MD simulations,17-19,22,25-34 whereas

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 4 of 34

pure liquid water hardly crystallizes into ice.50-53 The fast formation of methane hydrate compared with ice can be attributed to the fact that the first hydration shell of a hydrophobic molecule is highly structured.

A question is raised as to whether the

presence of THF in the aqueous phase also enhances crystal formation in a way similar to gas molecules. hydrates.

THF hydrate has been used in many experimental studies in place of gas

However, it is not sure whether the microscopic mechanism of formation of

THF hydrate is indeed similar to that of gas hydrates. Liang and Kusalik reported that formation of H2S hydrate is faster than methane hydrate in MD simulations owing to its relatively high solubility in water.23,24

This fact

infers that the formation of THF hydrate is even faster because THF is completely miscible with water.

If THF hydrate rapidly forms in MD simulations as expected, it would be

possible to examine nucleation of clathrate hydrate at a low computational cost by using THF instead of gas molecules. There are only a limited number of MD studies which examine the growth process of THF hydrate. Nada showed that THF molecules are trapped in open small cages on the hydrate surface.45

The author suggested that the arrangement of THF molecules at the

small cage sites might accelerate the formation of empty small cages in THF hydrate because water molecules located between guest molecules prefer to form stable four-coordinated structures.

Wu et al. carried out MD simulations of THF hydrate in

contact with an aqueous THF solution for several THF concentrations.46

They

demonstrated that the growth rate of THF hydrate is maximized at a certain THF concentration for a given temperature.

It was observed that the concentration of THF

4 ACS Paragon Plus Environment

Page 5 of 34

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

molecules near the growing hydrate surface is higher than the average value when THF hydrate is in contact with its melt.

The authors suggested that desorption of THF

molecules trapped in wrong sites on the hydrate surface is the rate limiting step of the crystal growth.

The effects of the THF molecules trapped on the hydrate surface in these

two MD studies are completely different from each other.

Therefore, the conflict on the

role of THF molecule should be settled by other methods than direct observations of crystal growth, such as the free energy calculation.

If the attractive forces between the incorrectly

trapped THF molecules and the hydrate surface are strong compared with the thermal energy, kT, the crystal growth of THF hydrate would be very slow.

However, if the

attractive force is weak, the trapped THF molecules would help the formation of empty cages as suggested by Nada, and THF hydrate could grow faster than ice and gas hydrates. Experimentally, in many cases, crystal growth is governed by the rate of removal of latent heat from the freezing interface.

Some simulation studies focused on the effects of

such macroscopic heat transfer on formation and dissociation of clathrate hydrates.25,30,54-56 The rate in the absence of heat transfer is also an important property to understand the intrinsic kinetics of crystal formation and dissociation.57-60

It is easy to evaluate the

intrinsic rate using isothermal MD simulations because the coupling between the system and the thermal bath can be strong.

In contrast, experimental growth rates are affected by

the heat transfer because the rate of agitation is quite slow compared with the timescale of MD simulations.

It was experimentally shown that crystal growth is slower for THF

hydrate than for ice at the same subcooling.4,61-63

However, it is unclear how significant

the macroscopic heat transfer is in contributing to the measured rates and whether the

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 6 of 34

intrinsic growth rate of THF hydrate is slower than ice. In this paper, we perform crystal growth simulations of THF hydrate and ice Ih, and measure the intrinsic growth rates.

The early MD studies did not quantitatively examine

the effects of the diffusion coefficients of THF and water on the crystal growth rate of THF hydrate.45,46

Therefore, we analyze the growth rates using a Wilson-Frenkel type of

kinetic model which accounts for the temperature dependence of the diffusion coefficients.59,60

The comparison between THF hydrate and ice reveals the large effect of

the adsorption of THF molecules in open small cages at the hydrate surface.

We calculate

the free energy profile for transferring a THF molecule from the aqueous phase to the hydrate surface in order to examine binding sites and the adsorption affinity of THF on the surface.

We demonstrate that THF molecules inhibit the crystal growth in a way similar to

a type of LDHIs.

Computational details We perform MD simulations of the crystal growth of THF hydrate and ice. molecules are described by the TIP4P/Ice model.64

Water

This model well reproduces the

melting temperature of ice Ih and the three phase equilibrium temperature of methane hydrate.41

The structure and potential parameters for THF are taken from the OPLS

united-atom force field.65

In an OPLS THF molecule, the charge of the oxygen atom is

-0.5e and those of the two neighboring carbons are 0.25 e.

We present a snapshot of a MD

simulation of a mixture of OPLS THF and TIP4P/Ice water at 300 K in Figure 1a.

6 ACS Paragon Plus Environment

The

Page 7 of 34

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

system consists of 8704 water and 512 THF molecules.

This composition is the same as

that of THF hydrate.

THF molecules are randomly mixed with water molecules in the

initial configuration.

However, as shown in Figure 1a, THF molecules aggregate in the

solution within the simulation time of 10 ns. The OPLS model was parameterized to reproduce properties of pure THF solution.

The dipole moment of each THF molecule in

aqueous solutions should be larger than that of the original OPLS model because the dielectric constant of water is much higher than that of THF.

Figure 1.

Final configurations of the simulations of the water-THF mixture using (a) the

original OPLS force field and (b) the modified force field.

THF molecules are shown by

red pentagons and the hydrogen bond network of water is represented by gray lines.

(c)

Solvation free energy of THF in liquid water against the scaling factor for the partial charges, α.

The horizontal dashed line indicates the experimental value of -14.5 kJ mol-1.

We evaluate the solvation free energy of THF in liquid water using the acceptance ratio method with the soft core potential.66-68

The system is composed of one THF

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

molecule and 511 water molecules in this simulation.

Page 8 of 34

The calculated hydration free

energy of -9.1 kJ mol-1 is much higher than the experimental value of -14.5 kJ mol-1.69

To

improve the force field, we scale the partial charges of atoms in THF with a scaling factor α. Figure 1c shows the α dependence of the hydration free energy of THF. free energy decreases with increasing α.

The hydration

It is found that the hydration free energy agrees

with the experimental value when α = 1.136.

We perform a MD simulation with this

scaling factor for 20 ns starting from the configuration shown in Figure 1a. As shown in Figure 1b, the modified THF molecules mix with TIP4P/Ice water.

We also check the

self-diffusion coefficients of THF (DTHF) and water (Dwater) in the aqueous THF solution. We obtain DTHF/Dwater = 0.44 for the original OPLS force field.

This value is very close to

the experimental value of 0.43,70,71 indicating that the diffusion dynamics of THF in water is well reproduced by the OPLS model.

We find that the effect of the change in α on the

diffusion dynamics is very small: DTHF/Dwater is 0.42 for the modified force field.

All the

following simulations are performed with the modified force field. The melting point and the crystal growth rate of THF hydrate are determined by the direct coexistence method.41,72,73 hydrate/solution system.

Figure 2a presents the initial configuration of the

The hydrate slab is a 4 × 4 × 4 unit cell replica of structure II

THF hydrate (512 THF and 8704 water molecules).

The numbers of THF and water

molecules in the aqueous solution are the same as those in the hydrate slab. molecules are randomly located in the aqueous phase in the initial configuration. of the simulation cell is roughly 69 × 69 × 136 Å3.

THF

The size

The growth rate of ice in liquid water

is also determined by the direct coexistence method.

The ice/water system consists of

8 ACS Paragon Plus Environment

Page 9 of 34

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

19440 water molecules.

The cell size is approximately 68 × 71 × 129 Å3.

The half of

the system is ice Ih and the remaining is liquid water in the initial configuration. simulations are performed using the GROMACS 4.6 package.74,75 interactions are treated by the particle mesh Ewald method.76,77

MD

Long range coulomb

The pressure is kept at 1

bar using the Parrinello-Rahman method,78,79 and the temperature is controlled by the Nose-Hoover method.80,81

The time constant for the coupling between the system and the

thermostat is set to 0.6 ps, which ensures rapid removal of latent heat from the growing surface.

The simulations of the hydrate/solution system and the ice/water system are

carried out in the temperature ranges of 244 ≤ T ≤ 280 K and 228 ≤ T ≤ 274 K, respectively. MD simulations are performed for three different proton-disordered configurations at each temperature.

A simulation is continued for 700 ns or until complete disappearance of

either solid or liquid part.

Figure 2. Snapshots of a crystal growth simulation of THF hydrate in its melt at (a) t = 0 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

ns and (b) t = 700 ns.

The temperature is 268 K.

represented by gray lines.

Page 10 of 34

The hydrogen bond network of water is

Purple and black pentagons are THF molecules in the hydrate

and in the aqueous solution, respectively.

Blue spheres show occupied small cages.

A

purple sphere in panel b is an unoccupied large cage in the hydrate slab.

We categorize THF molecules into two types.

First, the F3 parameter, which is low

for highly tetrahedral structures,82 is evaluated for water molecules in the first solvation shell of a THF molecule (r < 6 Å).

If the number of surrounding water molecules with a

highly tetrahedral structure (F3 < 0.3) is more than 23, the THF molecule is regarded as being in a large cage in the hydrate phase.

Otherwise, the THF molecule is categorized as

a molecule dissolved in the aqueous phase.

Figure 2 demonstrates that THF molecules are

well classified by this method.

Results and Discussion Figure 2b presents a snapshot of the hydrate/solution system at t = 700 ns in a simulation performed at T = 268 K.

The hydrate slab grows at this temperature, indicating

that the melting point of THF hydrate is higher than 268 K.

There is an unoccupied large

cage in the hydrate, which is shown as a purple sphere in Figure 2b. unoccupied large cages are also observed in other trajectories.

A small number of

Blue spheres represent

occupied small cages in the hydrate slab or at the hydrate surface.

We refer to a small

cage as one which is occupied by a THF molecule if the distance between the center of the

10 ACS Paragon Plus Environment

Page 11 of 34

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

cage and the center of mass of the THF molecule is less than 2.5 Å.

It is shown that all

small cages in the bulk region of the hydrate are empty whereas there are a large number of occupied open small cages at the growing surface. Figure 3a shows the potential energy of the hydrate/solution two-phase coexistence system for several temperatures.

The potential energy increases with time at T = 276 K

because of melting while it decreases at T =274 K and 272 K due to crystal growth.

Three

simulations are performed at every temperature and similar results are obtained in all cases. Therefore, we determine the melting point of THF hydrate at 1 bar as 275 K. temperature is very close to the experimental value of 277.6 K.

This melting

The corresponding results

for the ice/water two-phase coexistence system are shown in Figure 3b.

The potential

energy does not change at T = 270 K for a long time, and thus we determine this temperature as the melting point of ice.

This result is in agreement with early simulation

studies of Vega et al.64,73

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 34

Figure 3 (a) Evolution of the potential energy as a function of time for the THF hydrate/solution two-phase coexistence system.

(b) Evolution of the potential energy for

the ice/water system.

In order to obtain the growth rates of THF hydrate and ice, we evaluate the position of the solid/liquid interface at each time step.

It is possible to distinguish the solid region

from the aqueous solution using the F3 parameter because the tetrahedrality of water molecules is much higher in ice and clathrate hydrates than in aqueous solutions.

Figure

4a is a typical example of the distribution function of water molecules with a highly tetrahedral structure (F3 < 0.3) in the hydrate/solution system.

The figure shows that the

hydrate slab occupies a region of -40 Å < z < 40 Å in the simulation cell. We determine the position of the surface located at z > 0 by fitting the distribution function with

12 ACS Paragon Plus Environment

Page 13 of 34

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

f + ( z) =

1 (ρs + ρl ) − 1 (ρs − ρl ) tanh (z − z0 )  , 2 2  w 

(1)

where z0 is the surface position, w represents the width, and ρs and ρl are the densities of the water molecules with F3 < 0.3 in the solid and solution phases, respectively.

The position

of the surface at z < 0 is determined by the fitting with

f− ( z) =

1 (ρs + ρl ) + 1 (ρs − ρl ) tanh ( z − z0 )  . 2 2  w 

A result of the fitting is shown by the black dashed curve in Figure 4a. of z0 is plotted in Figure 4b.

(2)

The time evolution

We define the slope of z0(t) as the growth rate of the surface.

Figure 4 (a) Distribution function of water molecules with a highly tetrahedral structure in the hydrate/solution system at t = 300 ns in a simulation at T = 268 K (red curve). fitting curve of Eq. 1 is shown by the black dashed curve.

A

(b) Time evolution of the

surface position z0 of the hydrate slab obtained from a simulation performed at T = 268 K.

In Figure 5a, we show the growth rates of ice and THF hydrate against the degree of subcooling, ∆T = (Tm – T), where Tm is the melting temperature.

13 ACS Paragon Plus Environment

Six points are plotted at

The Journal of Physical Chemistry

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 34

each temperature because three simulations are performed and there are two surfaces in a system.

Figure 5a demonstrates that the crystal growth of THF hydrate is significantly

slower than ice.

At ∆T = 10 K, for example, the growth rate of THF hydrate is

approximately 20 times slower than that of ice.

Experimental data showed that the growth

rate of THF hydrate is one or two order magnitude slower than that of ice when they are compared at the same subcooling.4,61-63

The simulation results are in agreement with the

experimental results, although the experimental rates are affected by the macroscopic heat transfer effect and it is unclear that they can be directly compared with the simulation results.

Figure 5 (a) Growth rates of ice (blue triangles) and THF hydrate (red circles) against the degree of subcooling, ∆T = (Tm – T). K for THF hydrate.

The melting temperature Tm is 270 K for ice and 275

The inset maginifies the growth rates of THF hydrate.

The solid

curves are obtained from the fitting to the Wilson-Frenkel equation. (b) Diffusion 14 ACS Paragon Plus Environment

Page 15 of 34

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

coefficients of water molecules in pure water (blue triangles), water molecules in the aqueous THF solution (red squares), and THF molecules in the aqueous THF solution (red circles).

The thermodynamic driving force for crystal growth, i.e., the difference in the molar free energy between the solution and solid phases, becomes larger as the temperature decreases.

In contrast, molecular motions are suppressed at lower temperatures.

As a

result, the growth rate of a crystal first increases and then decreases as the temperature is lowered from the melting point. Figure 5a. behavior.59,60

Both ice and THF hydrate exhibit this trend as shown in

The Wilson-Frenkel theory is a simple model that can reproduce this In this model, the activation energy for the transfer of a molecule from the

solution phase to the solid surface is assumed to be the same as that for the molecular diffusion in the bulk solution phase.

The growth rate is expressed as

  − ∆H f (Tm − T )   , k (T ) = CD(T )1 − exp k TT  B m  

(3)

where D(T) , ∆Hf, and C are the diffusion coefficient, the enthalpy of fusion, and a constant, respectively.

The diffusion coefficients of water and THF are shown in Figure 5b.

Diffusion dynamics of the THF solution is several fold slower than that of pure water. This slower diffusion dynamics cannot solely explain the remarkably slow growth of THF hydrate because the growth of THF hydrate can be 20 times slower than that of ice. the growth rates with Equation (2) using D(T) and ∆Hf as inputs. ice are 4.6 and 5.1 kJ mol-1, respectively.

We fit

∆Hf of THF hydrate and

As shown by the blue solid curve in Figure 5a,

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 34

the growth rates of ice can be approximated by the Wilson-Frenkel equation. agreement with early MD studies.83,84

This is in

In contrast, the fitting curve significantly deviates

from the growth rates of THF hydrate (D(T) of THF is used for the fitting of the growth rates because it is lower than D(T) of water)85.

This result demonstrates that the growth

rate of THF hydrate is governed by a factor (or factors) other than diffusion coefficients. Figure 2b shows that many THF molecules are located in open small cages at the growing hydrate surface.

Nada inferred that the trapped THF at the surface may help the

formation of empty small cages during the growth process.45

The results of the present

study suggest that this effect is insignificant for the kinetics of crystal growth of THF hydrate.

Small cages in the bulk region of the hydrate cannot incorporate THF molecules,

and the THF molecules trapped in the open small cages at the surface must be excluded as the hydrate slab grows.

Wu et al. suggested that this desorption process is the rate

limiting step of the crystal growth of THF hydrate.46

They reported that the local

concentration of THF near the hydrate surface is higher than the average value and the local concentration slowly decreases with the growth of THF hydrate.

However, in their study,

the adsorption affinity of THF molecules on the hydrate surface was not examined.

The

precise positions of the trapped THF molecules were also not shown. To examine the adsorption mechanism of THF in details, we calculate the free energy profile of a THF molecule dissolved in the aqueous solution approaching the hydrate surface using the umbrella sampling technique.

The system employed to compute the free

energy profile consists of a 2 × 2 × 2 unit cell replica of THF hydrate (1088 water molecules and 64 THF molecules) and an aqueous THF solution which is composed of

16 ACS Paragon Plus Environment

Page 17 of 34

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

2176 water and 128 THF molecules.

Figure 6a presents a snapshot of the simulation cell.

The z coordinate of the center of mass of all guest molecules in the hydrate slab is defined as z = 0.

A THF molecule in the aqueous phase is selected as a probe and fixed around a

given z value using a harmonic potential with a force constant of 1000 kJ mol-1 nm-2. No constrains are imposed along the x and y directions.

The simulation is performed for

every 0.2 Å in the z direction, and the free energy profile of the probe THF molecule is calculated from the weighted histogram analysis method.86,87 must be carried out for equilibrium conditions.

Free energy calculations

If the hydrate slab grows or dissociates

during the simulation, the obtained free energy profile would change as time evolves. difficult to suppress the hydrate growth at a temperature lower than the melting point.

It is In

contrast, it is easy to prevent the hydrate slab from dissociation by locating guest molecules at the fixed positions in the hydrate.

Therefore, the simulations are performed at a

temperature slightly higher than the melting point with fixing all THF molecules in the hydrate slab at the lattice points using a harmonic potential of a force constant of 10000 kJ mol-1 nm-2.

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 34

Figure 6 (a) Snapshot of the simulation cell for the calculation of the free energy profile. (b) Free energy profile for transferring a THF molecule (red). energy profile around z = 20 Å.

The inset magnifies the free

The solid blue curve is ∆Ep(z) = Ep(z) – Ep(∞), where Ep

is the interaction energy between the probe THF molecule and all other molecules.

(c)

Distribution function of unconstrained THF molecules in the aqueous phase.

The red curve in Figure 6b shows the free energy profile of the probe THF molecule at T = 278 K.

Two minima are found at z = 19.4 and 18.6 Å.

The outermost layer of

occupied large cages in the hydrate slab is at z = 15.1 Å and the interval between the layers of large cages is roughly 4.3 Å.

The minimum at z = 19.4 Å thus corresponds to the layer

of open large cages at the hydrate surface.

Figure 7a presents a typical snapshot of the

probe THF molecule fixed at z = 19.4 Å.

The probe THF molecule is in an open large 18

ACS Paragon Plus Environment

Page 19 of 34

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

cage at the hydrate surface. fixed at z = 18.6 Å. hydrate surface.

Figure 7b presents a snapshot of the probe THF molecule

The probe THF molecule is incorporated in an open small cage at the

As pointed by the blue and red arrows in Figure 7a, unconstrained THF

molecules are located in open large and open small cages at the surface, indicating that these are indeed stable adsorption sites for THF. probe molecule on the xy plane.

Figure 8 presents the trajectory of the

The probe molecule bound at z = 18.6 Å diffuses on the

xy plane but it is frequently trapped in the small cage sites during the simulation time of 80 ns (a similar result is obtained for the probe molecule bound at z = 19.6 and 14.4 Å). Figure 6b shows that the free energy of a THF molecule at z = 18.6 Å is almost same as that at z = 19.4 Å.

That is, THF molecules favor the open small cages as well as the open large

cages at the hydrate surface.

The THF molecules in the open small cages must be

removed when the hydrate grows. energy barrier at z = 21.2 Å.

However, this process is unfavorable due to the free

The height of the barrier at z = 21.2 Å is 2.5 kJ mol-1, which

is comparable to the thermal energy at 278 K of 2.3 kJ mol-1.

One may anticipate that the

free energy barrier at z = 21.2 Å also hinders the transfer of a THF molecule from the solution to the hydrate surface. The barrier is, however, lower than the thermal energy for this transfer (from right to left in Figure 6b), and thus THF molecules in the solution can easily approach the surface.

The local concentration of THF molecules is higher near the

surface than in the bulk region due to the negative value of the free energy at z ~ 19 Å. Figure 6c shows the distribution function of unconstrained THF molecules in the aqueous phase.

It is evident that the concentration of THF molecules is indeed higher in the

surface region (z ~ 19 Å) than in the bulk region (z > 25 Å).

19 ACS Paragon Plus Environment

This result is in agreement

The Journal of Physical Chemistry

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 20 of 34

with the early study by Wu et al.46

Figure 7 Snapshots of the probe THF molecule fixed at (a) z = 19.4 Å, (b) 18.6 Å, and (c) 14.4 Å.

The probe THF molecule is colored red.

THF molecules in the hydrate and

those in the solution are represented by purple and black pentagons.

Panels (a) and (b) are

images looking down onto the (110) surface whereas panel (c) is that onto the (1 10) surface.

The black pentagons indicated by the blue arrow and the red arrow in panel (a)

are an unconstrained THF molecule located in a large cage and that in a small cage at the hydrate surface.

20 ACS Paragon Plus Environment

Page 21 of 34

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

Figure 8 Trajectory of the probe molecule on the xy plane. at z = 18.6 Å.

The probe molecule is bound

The black and purple dots indicate the small and large cage sites on the

surface.

Figure 6b also shows that the free energy profile has a deep minimum at z = 14.4 Å. A snapshot of the probe THF molecule fixed at z = 14.4 Å is presented in Figure 7c. There is a layer of small cages at z = 13 Å.

The THF molecule is incorporated in one of

small cages in this layer in Figure 7c (the location of the minimum in the free energy profile, z = 14.4 Å, dose not coincide with the position of the center of the small cage, z = 13 Å, simply because the size of a THF molecule is larger than the inner diameter of a small cage).

The small cages at z = 13 Å are not completely closed because the large

cages at z = 19.4 Å are not fully occupied. these small cages.

Therefore, THF molecules can be included in

Once a THF molecule is trapped in an open small cage at z = 13 Å, the

THF molecule needs to overcome the free energy barrier at z = 16.6 Å in addition to that at z = 21.2 Å for the growth of hydrate.

The height of the free energy barrier at z = 16.6 Å is 21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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.2 kJ mol-1. quite slow.

Page 22 of 34

This is clearly higher than the thermal energy, indicating that this process is The present analysis demonstrates that the slow growth of THF hydrate is

indeed due to the trapping of THF molecules in open small cages at the hydrate surface. Recently, we investigated the adsorption mechanism of molecules on the surface of clathrate hydrate.88

The mechanism shown in ref. 88 is common for various nonpolar and

weakly polar molecules, and thus it can explain the preference of THF molecules in open cages at the surface.

It is possible to consider the free energy profile of a hard sphere

solute transferring from the bulk liquid region to the hydrate surface. solute has no attractive interaction with surrounding molecules.

The hard sphere

The hydration free energy

of a hard sphere solute, which is purely entropic, is given by G hyd = − kT ln P ,

(4)

where P is the insertion probability defined as the ratio of the volume of cavities that can accommodate the hard sphere solute in a unit volume.89,90

The free energy profile thus

can be expressed as

∆G( z ) = Ghyd ( z ) − Ghyd (∞) = −kT ln where P (z) is the z dependent insertion probability.

P( z ) , P(∞)

(5)

Open cages at the hydrate surface are

larger than voids which spontaneously occur in the bulk aqueous phase.

In other words,

the insertion probability is larger at the hydrate surface than in the bulk aqueous phase. Thus, the free energy of the solute located in an open cage at the hydrate surface can be lower than that in the bulk aqueous phase even when the solute has no attractive interactions with surrounding molecules. (Note that not only the minima but also the barriers in the free energy profile are explained from Eq. (5). 22 ACS Paragon Plus Environment

The possibility of formation

Page 23 of 34

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

of voids which can include a THF molecules is low at z = 16.6 and 21.2 Å because hydrogen bonds of water tend to occur there to form cages based on the existing crystalline structure.) We demonstrated that this is the dominant contribution to the preference of nonpolar and weakly polar molecules for the hydrate surface in ref. 88.

Other factors,

such as hydrogen bonding and van der Waals interaction between the adsorbed molecule and water, are insignificant.

In Figure 6b, we plot the change in the interaction energy

between the probe THF molecule and other molecules.

The energy profile does not

change by the transfer of the THF molecule, showing that the energetic term is insignificant for the change in the free energy. One may think that the change in the hydrogen bonding between water molecules contributes to the free energy change.

Wu et al. reported the enhanced hydrogen bonding

between water molecules caused by the presence of THF molecules near the hydrate surface.

This indicates that the interaction energy between water molecules, Ew, decreases

as THF molecules approach the hydrate surface (∆Ew < 0).

However, the enhancement of

the hydrogen bonding results in a decrease in entropy at the same time (T∆Sw < 0).

As

shown by Yu and Karplus, the decrease in the water-water energy is identical to the change in the entropy arising from the water-water interactions (∆Ew = T∆Sw).91 Therefore, the change in the hydrogen bonding between water molecules makes no contribution to the free energy change.92,93 In the free energy calculation, a flat (001) surface is exposed to the solution.

The free

energy profiles of other crystallographic faces would be different from that of the (001) surface.

The presence of surface defects such as kink and step sites would also affect the

23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

free energy profile.

Page 24 of 34

The investigation of various surfaces is beyond the scope of this study

because it requires a large computational cost.

However, it is certain from Eq. (5) that

THF molecules are trapped on any surface of THF hydrate because there must be open small cages and inclusion of THF molecules in them causes entropic stabilization. We finally mention the similarity between the mechanism of the slow growth of THF hydrate and the inhibition mechanism of a type of LDHIs known as kinetic hydrate inhibitors (KHIs). surface.3,94

KHIs are water-soluble polymers which adsorb on the hydrate

In the previous paper, we showed that the free energy of a monomer of a

common KHI, PVCap, in an open cage at the hydrate surface is lower than that in bulk water by 10.5 kJ mol-1 mainly due to the mechanism mentioned above.88 Because KHIs cannot be included in a hydrate crystal, the growth of the hydrate is inhibited by the adsorption of the KHIs on the hydrate surface.

This is quite similar to the mechanism of

the slow growth of THF hydrate: the crystal growth of THF hydrate is inhibited because THF molecules prefer the open small cages at the hydrate surface although the small cages in the bulk hydrate phase cannot contain THF molecules.

Conclusions We have investigated the growth mechanism of THF hydrate using MD simulations. The force field of THF is modified to reproduce the hydration free energy, and it yields a melting temperature of THF hydrate very close to the experimental one. the crystal growth of THF hydrate is much slower than that of ice.

24 ACS Paragon Plus Environment

It is found that The temperature

Page 25 of 34

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

dependence of the growth rate of THF hydrate significantly deviates from the Wilson-Frenkel model whereas that of ice is well approximated.

Both the slow growth

and the deviation from the Wilson-Frenkel model are attributed to the adsorption of THF molecules in open small cages at the hydrate surface.

It is shown that a THF molecule

trapped in an open small cage must cross one or two free energy barriers to escape from the surface region. THF hydrate is often used to evaluate the efficiency of KHIs as an analogue to methane hydrate.6,7

There are almost no methane molecules in the aqueous phase near the

growing surface of methane hydrate due to the low solubility of methane in water.

In this

case, the efficiency of a KHI, which depends on the adsorption affinity of the KHI on the hydrate surface, would be determined by a matching between the surface structure and the molecular shape of the KHI to a large extent.

In contrast, there are many THF molecules

near and on the growing surface of THF hydrate.

There could be a competitive effect

between the adsorption of THF and that of KHIs on the hydrate surface because the adsorption mechanisms of them are the essentially same.

If this effect is significant, the

use of THF hydrate instead of methane hydrate may be inadequate for the evaluation of KHIs.

MD simulations of a KHI near the interface between THF hydrate and an aqueous

THF solution would provide microscopic insights into the legitimacy on the substitution of THF hydrate for methane hydrate.

Acknowledgements

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 present work was supported by a Grant-in-Aid by JSPS and by HPCI Strategic Programs for Innovative Research (SPIRE) and Computational Materials Science Initiative (CMSI), MEXT, Japan.

Calculations were performed on the K computer at the RIKEN

Advanced Institute for Computational Science (Project ID: hp150221) and the computers at Research Center for Computational Science, Okazaki, Japan.

26 ACS Paragon Plus Environment

Page 26 of 34

Page 27 of 34

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

References (1) Sloan, E. D.; Koh, C. A. Clathrate hydrates of natural gases; CPC Press: Boca Raton, 2008. (2) Makogon, T. Y.; Larsen, R.; Knight, C. A.; Sloan, E. D. Melt growth of tetrahydrofuran clathrate hydrate and its inhibition: Method and first results. J. Cryst. Growth 1997, 1997 179, 258-262. (3) Larsen, R.; Knight, C. A.; Sloan, E. D. Clathrate hydrate growth and inhibition. Fluid Phase

Equilib. 1998, 1998 150–151, 353-360. (4) Bollavaram, P.; Devarakonda, S.; Selim, M. S.; Sloan, E. D. Growth kinetics of single crystal sII hydrates: Elimination of mass and heat transfer effects. Ann. N.Y. Acad. Sci. 2000, 2000 912, 533-543. (5) Zeng, H.; Wilson, L. D.; Walker, V. K.; Ripmeester, J. A. Effect of antifreeze proteins on the nucleation, growth, and the memory effect during tetrahydrofuran clathrate hydrate formation.

J. Am. Chem. Soc. 2006, 2006 128, 2844-2850. (6) Kelland, M. A. History of the development of low dosage hydrate inhibitors. Energy Fuels 2006, 2006 20, 825-847. (7) Perrin, A.; Musa, O. M.; Steed, J. W. The chemistry of low dosage clathrate hydrate inhibitors. Chem. Soc. Rev. 2013, 2013 42, 1996-2015. (8) Strobel, T. A.; Koh, C. A.; Sloan, E. D. Thermodynamic predictions of various tetrahydrofuran and hydrogen clathrate hydrates. Fluid Phase Equilib. 2009, 2009 280, 61-67. (9) Hester, K. C.; Strobel, T. A.; Sloan, E. D.; Koh, C. A.; Huq, A.; Schultz, A. J. Molecular hydrogen occupancy in binary THF−H2 clathrate hydrates by high resolution neutron diffraction. J. Phys. Chem. B 2006, 2006 110, 14024-14027. (10) Sloan, E. D. Fundamental principles and applications of natural gas hydrates. Nature 2003, 2003 426, 353-359. (11) Mao, W. L.; Mao, H. K. Hydrogen storage in molecular compounds. Proc. Natl. Acad. Sci.

U.S.A. 2004, 2004 101, 708-710. (12) Mao, W. L.; Mao, H.-K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Hydrogen clusters in clathrate hydrate. Science 2002, 2002 297, 2247-2249. (13) Struzhkin, V. V.; Militzer, B.; Mao, W. L.; Mao, H. K.; Hemley, R. J. Hydrogen storage in molecular clathrates. Chem. Rev. 2007, 2007 107, 4133-4151. (14) 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.

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Science 2004, 2004 306, 469-471. (15) Kang, S.-P.; Lee, H. Recovery of CO2 from flue gas using gas hydrate:  Thermodynamic verification through phase equilibrium measurements. Environ. Sci. Technol. 2000, 2000 34, 4397-4400. (16) Linga, P.; Kumar, R.; Englezos, P. The clathrate hydrate process for post and pre-combustion capture of carbon dioxide. J. Hazard. Mater. 2007, 2007 149, 625-629. (17) Moon, C.; Taylor, P. C.; Rodger, P. M. Molecular dynamics study of gas hydrate formation. J.

Am. Chem. Soc. 2003, 2003 125, 4706-4707. (18) Hawtin, R. W.; Quigley, D.; Rodger, P. M. Gas hydrate nucleation and cage formation at a water/methane interface. Phys. Chem. Chem. Phys. 2008, 2008 10, 4853-4864. (19) Vatamanu, J.; Kusalik, P. G. Unusual crystalline and polycrystalline structures in methane hydrates. J. Am. Chem. Soc. 2006, 2006 128, 15588-15589. (20) Vatamanu, J.; Kusalik, P. G. Molecular insights into the heterogeneous crystal growth of sI methane hydrate. J. Phys. Chem. B 2006, 2006 110, 15896-15904. (21) Vatamanu, J.; Kusalik, P. G. Heterogeneous crystal growth of methane hydrate on its sII [001] crystallographic face. J. Phys. Chem. B 2008, 2008 112, 2399-2404. (22) Vatamanu, J.; Kusalik, P. G. Observation of two-step nucleation in methane hydrates.

Phys. Chem. Chem. Phys. 2010, 2010 12, 15065-15072. (23) Liang, S.; Kusalik, P. G. Crystal growth simulations of H2S hydrate. J. Phys. Chem. B 2010, 2010 114, 9563-9571. (24) Liang, S.; Kusalik, P. G. Exploring nucleation of H2S hydrates. Chem. Sci. 2011, 2011 2, 1286-1292. (25) Liang, S.; Kusalik, P. G. Nucleation of gas hydrates within constant energy systems. J.

Phys. Chem. B 2013, 2013 117, 1403-1410. (26) Pirzadeh, P.; Kusalik, P. G. Molecular insights into clathrate hydrate nucleation at an ice-solution interface. J. Am. Chem. Soc. 2013, 2013 135, 7278-7287. (27) Lauricella, M.; Meloni, S.; Liang, S.; English, N. J.; Kusalik, P. G.; Ciccotti, G. Clathrate structure-type recognition: Application to hydrate nucleation and crystallisation. J. Chem. Phys. 2015, 2015 142, 244503. (28) Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Microsecond simulations of spontaneous methane hydrate nucleation and growth. Science 2009, 2009 326, 1095-1098. (29) Walsh, M. R.; Beckham, G. T.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Methane hydrate nucleation rates from molecular dynamics simulations: Effects of aqueous methane

28 ACS Paragon Plus Environment

Page 28 of 34

Page 29 of 34

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

concentration, interfacial curvature, and system size. J. Phys. Chem. C 2011 2011, 115, 21241-21248. (30) 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, 015 17, 8870-8876. (31) Sarupria, S.; Debenedetti, P. G. Homogeneous nucleation of methane hydrate in microsecond molecular dynamics simulations. J. Phys. Chem. Lett. 2012, 2012 3, 2942-2947. (32) Jiménez-Ángeles, F.; Firoozabadi, A. Nucleation of methane hydrates at moderate subcooling by molecular dynamics simulations. J. Phys. Chem. C 2014, 2014 118, 11310-11318. (33) Jiménez-Ángeles, F.; Firoozabadi, A. Enhanced hydrate nucleation near the limit of stability. J. Phys. Chem. C 2015, 2015 119, 8798-8804. (34) English, N. J.; Lauricella, M.; Meloni, S. Massively parallel molecular dynamics simulation of formation of clathrate-hydrate precursors at planar water-methane interfaces: Insights into heterogeneous nucleation. J. Chem. Phys. 2014, 2014 140, 204714. (35) English, N. J.; MacElroy, J. M. D. Theoretical studies of the kinetics of methane hydrate crystallization in external electromagnetic fields. J. Chem. Phys. 2004, 2004 120, 10247-10256. (36) Jacobson, L. C.; Hujo, W.; Molinero, V. Nucleation pathways of clathrate hydrates: Effect of guest size and solubility. J. Phys. Chem. B 2010, 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, 2010 132, 11806-11811. (38) Jacobson, L. C.; Molinero, V. Can amorphous nuclei grow crystalline clathrates? The size and crystallinity of critical clathrate nuclei. J. Am. Chem. Soc. 2011, 2011 133, 6458-6463. (39) Nguyen, A. H.; Jacobson, L. C.; Molinero, V. Structure of the clathrate/solution interface and mechanism of cross-nucleation of clathrate hydrates. J. Phys. Chem. C 2012, 2012 116, 19828-19838. (40) Knott, B. C.; Molinero, V.; Doherty, M. F.; Peters, B. Homogeneous nucleation of methane hydrates: Unrealistic under realistic conditions. J. Am. Chem. Soc. 2012, 2012 134, 19544-19547. (41) Conde, M. M.; Vega, C. Determining the three-phase coexistence line in methane hydrates using computer simulations. J. Chem. Phys. 2010, 2010 133, 064507. (42) 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, 2015 142, 044501. (43) Nada, H. Growth mechanism of a gas clathrate hydrate from a dilute aqueous gas solution: A molecular dynamics simulation of a three-phase system. J. Phys. Chem. B 2006, 2006 110,

29 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

16526-16534. (44) Tung, Y.-T.; Chen, L.-J.; Chen, Y.-P.; Lin, S.-T. Molecular dynamics study on the growth of structure I methane hydrate in aqueous solution of sodium chloride. J. Phys. Chem. B 2012, 2012 116, 14115-14125. (45) Nada, H. Anisotropy in growth kinetics of tetrahydrofuran clathrate hydrate: A molecular dynamics study. J. Phys. Chem. B 2009, 2009 113, 4790-4798. (46) Wu, J.-Y.; Chen, L.-J.; Chen, Y.-P.; Lin, S.-T. Molecular dynamics study on the equilibrium and kinetic properties of tetrahydrofuran clathrate hydrates. J. Phys. Chem. C 2015, 2015 119, 1400-1409. (47) Wu, J.-Y.; Chen, L.-J.; Chen, Y.-P.; Lin, S.-T. Molecular dynamics study on the growth mechanism of methane plus tetrahydrofuran mixed hydrates. J. Phys. Chem. C 2015, 2015 119, 19883-19890. (48) Alavi, S.; Ripmeester, J. A. Effect of small cage guests on hydrogen bonding of tetrahydrofuran in binary structure II clathrate hydrates. J. Chem. Phys. 2012, 2012 137, 054712. (49) Alavi, S.; Ripmeester, J. A.; Klug, D. D. Molecular-dynamics simulations of binary structure II hydrogen and tetrahydrofurane clathrates. J. Chem. Phys. 2006, 006 124, 14704. (50) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular dynamics simulation of the ice nucleation and growth process leading to water freezing. Nature 2002, 2002 416, 409-413. (51) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Spontaneous liquid-liquid phase separation of water. Phys. Rev. E 2014, 2014 89, 020301. (52) Espinosa, J. R.; Sanz, E.; Valeriani, C.; Vega, C. Homogeneous ice nucleation evaluated for several water models. J. Chem. Phys. 2014, 2014 141, 18C529. (53) English, N. J. Massively parallel molecular-dynamics simulation of ice crystallisation and melting: The roles of system size, ensemble, and electrostatics. J. Chem. Phys. 2014, 2014 141, 234501. (54) Bagherzadeh, S. A.; Englezos, P.; Alavi, S.; Ripmeester, J. A. Molecular simulation of non-equilibrium methane hydrate decomposition process. J. Chem. Thermodyn. 2012, 2012 44, 13-19. (55) Bagherzadeh, S. A.; Alavi, S.; Ripmeester, J. A.; Egnglezos, P. Evolution of methane during gas hydrate dissociation. Fluid Phase Equilib. 2013, 2013 358, 114-120. (56) Alavi, S.; Ripmeester, J. A. Nonequilibrium adiabatic molecular dynamics simulations of methane clathrate hydrate decomposition. J. Chem. Phys. 2010, 2010 132, 144703. (57) Kim, H. C.; Bishnoi, P. R.; Heidemann, R. A.; Rizvi, S. S. H. Kinetics of methane hydrate decomposition. Chem. Eng. Sci. 1987, 1987 42, 1645-1653.

30 ACS Paragon Plus Environment

Page 30 of 34

Page 31 of 34

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

(58) Clarke, M.; Bishnoi, P. R. Determination of the activation energy and intrinsic rate constant of methane gas hydrate decomposition. Can. J. Chem. Eng. 2001, 2001 79, 143-147. (59) Kirkpatrick, R. Crystal growth from the melt: A review. Am. Mineral. 1975, 1975 60, 798-814. (60) Jackson, K. The interface kinetics of crystal growth processes. Interface Sci. 2002, 2002 10, 159-169. (61) Barduhn, A. J.; Huang, J. S. Why does ice grow faster from salt water than from fresh water? Desalination 1987, 1987 67, 99-106. (62) Fernandez, R.; Barduhn, A. J. The growth rate of ice crystals. Desalination 1967, 1967 3, 330-342. (63) Iida, T.; Mori, H.; Mochizuki, T.; Mori, Y. H. Formation and dissociation of clathrate hydrate in stoichiometric tetrahydrofuran–water mixture subjected to one-dimensional cooling or heating. Chem. Eng. Sci. 2001, 2001 56, 4747-4758. (64) Abascal, J. L. F.; Sanz, E.; García Fernández, R.; Vega, C. A potential model for the study of ices and amorphous water: TIP4P/Ice. J. Chem. Phys. 2005, 2005 122, 234511. (65) Briggs, J. M.; Matsui, T.; Jorgensen, W. L. Monte carlo simulations of liquid alkyl ethers with the opls potential functions. J. Comput. Chem. 1990, 1990 11, 958-971. (66) Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. J.

Comput. Phys. 1976, 1976 22, 245-268. (67) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 1994 222, 529-539. (68) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Effects of thermodynamic inhibitors on the dissociation of methane hydrate: A molecular dynamics study. Phys. Chem. Chem. Phys. 2015, 2015

17, 32347-32357. (69) Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. Group contributions to the thermodynamic properties of non-ionic organic solutes in dilute aqueous solution. J. Solution Chem. 1981, 1981 10, 563-595. (70) Trappeniers, N. J.; Gerritsma, C. J.; Oosting, P. H. The self-diffusion coefficient of water, at 25°C, by means of spin-echo technique. Phys. Lett. 1965, 1965 18, 256-257. (71) Leaist, D. G.; MacEwan, K.; Stefan, A.; Zamari, M. Binary mutual diffusion coefficients of aqueous cyclic ethers at 25 °C. Tetrahydrofuran, 1,3-dioxolane, 1,4-dioxane, 1,3-dioxane, tetrahydropyran, and trioxane. J. Chem. Eng. Data 2000, 2000 45, 815-818. (72) Ladd, A. J. C.; Woodcock, L. V. Triple-point coexistence properties of the Lennard-Jones

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

system. Chem. Phys. Lett. 1977, 1977 51, 155-159. (73) García Fernández, R.; Abascal, J. L. F.; Vega, C. The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface. J. Chem. Phys. 2006, 2006 124, 144506. (74) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 2008 4, 435-447. (75) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 2005 26, 1701-1718. (76) Darden, T.; York, D.; Pedersen, L. Particle mesh ewald: An n⋅log(n) method for ewald sums in large systems. J. Chem. Phys. 1993, 1993 98, 10089-10092. (77) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh ewald method. J. Chem. Phys. 1995, 1995 103, 8577-8593. (78) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol.

Phys. 1983, 1983 50, 1055-1076. (79) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 1981 52, 7182-7190. (80) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol.

Phys. 1984, 1984 52, 255-268. (81) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 1985 31, 1695-1697. (82) Baez, L. A.; Clancy, P. Computer-simulation of the crystal-growth and dissolution of natural-gas hydrates. Ann. N.Y. Acad. Sci. 1994, 1994 715, 177-186. (83) Rozmanov, D.; Kusalik, P. G. Temperature dependence of crystal growth of hexagonal ice (Ih). Phys. Chem. Chem. Phys. 2011, 2011 13, 15501-15511. (84) Weiss, V. C.; Rullich, M.; Köhler, C.; Frauenheim, T. Kinetic aspects of the thermostatted growth of ice from supercooled water in simulations. J. Chem. Phys. 2011, 2011 135, 034701. (85) Tang, C.; Harrowell, P. Anomalously slow crystal growth of the glass-forming alloy CuZr.

Nature Mater. 2013, 2013 12, 507-511. (86) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. J.

Comput. Chem. 1992, 1992 13, 1011-1021. (87) Hub, J. S.; de Groot, B. L.; van der Spoel, D. G_wham—a free weighted histogram analysis

32 ACS Paragon Plus Environment

Page 32 of 34

Page 33 of 34

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

implementation including robust error and autocorrelation estimates. J. Chem. Theory Comput. 2010, 2010 6, 3713-3720. (88) Yagasaki, T.; Matsumoto, M.; Tanaka, H. Adsorption mechanism of inhibitor and guest molecules on the surface of gas hydrates. J. Am. Chem. Soc. 2015, 2015 137, 12079-12085. (89) Pohorille, A.; Pratt, L. R. Cavities in molecular liquids and the theory of hydrophobic solubilities. J. Am. Chem. Soc. 1990, 1990 112, 5066-5074. (90) Chandler, D. Interfaces and the driving force of hydrophobic assembly. Nature 2005, 2005 437, 640-647. (91) Yu, H. A.; Karplus, M. A thermodynamic analysis of solvation. J. Chem. Phys. 1988, 1988 89, 2366-2379. (92) Yagasaki, T.; Iwahashi, K.; Saito, S.; Ohmine, I. A theoretical study on anomalous temperature dependence of pKw of water. J. Chem. Phys. 2005, 2005 122, 144504. (93) Yagasaki, T.; Saito, S.; Ohmine, I. Effects of nonadditive interactions on ion solvation at the water/vapor interface: A molecular dynamics study. J. Phys. Chem. A 2010, 2010 114, 12573-12584. (94) Anderson, B. J.; Tester, J. W.; Borghi, G. P.; Trout, B. L. Properties of inhibitors of methane hydrate formation via molecular dynamics simulations. J. Am. Chem. Soc. 2005, 2005 127, 17852-17862.

33 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

TOC graphic

34 ACS Paragon Plus Environment

Page 34 of 34