Molecular Dynamics Study of Kinetic Hydrate ... - ACS Publications

Subscriber access provided by UNIV OF DURHAM

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Molecular Dynamics Study of Kinetic Hydrate Inhibitors: The Optimal Inhibitor Size and Effect of Guest Species Takuma Yagasaki, Masakazu Matsumoto, and Hideki Tanaka J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b09834 • Publication Date (Web): 31 Dec 2018 Downloaded from http://pubs.acs.org on January 1, 2019

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 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 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.

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

Molecular Dynamics Study of Kinetic Hydrate Inhibitors: The Optimal Inhibitor Size and Effect of Guest Species

Takuma Yagasaki, Masakazu Matsumoto, and Hideki Tanaka* Research Institute for Interdisciplinary Science, Okayama University, Okayama, 700-8530, Japan

1 Environment ACS Paragon Plus

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

Abstract We propose a model for slowing down of clathrate hydrate formation caused by kinetic hydrate inhibitors (KHIs) on the basis of the Gibbs-Thomson effect. The residence time of inhibitor molecules bound to the hydrate surface and the intrinsic growth rate of the clathrate hydrate without KHIs are key ingredients of the model. The binding free energies of the monomer, dimer, tetramer, and octamer of a KHI, polyvinylcaprolactam (PVCap), are calculated using molecular dynamics simulations to estimate the residence times which are far beyond the feasible simulation time. Our model accounts for the kinetic inhibition mechanism while reproducing experimental data of the size dependence of PVCap very well. We demonstrate that this model explains why blends of high and low molecular weight polymers show better performance than the KHI with a unimodal molecular weight distribution and why quaternary ammonium cations are good KHIs for tetrahydrofuran hydrate although they cannot inhibit formation of natural gas hydrates.


Page 2 of 40

Page 3 of 40 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


Introduction Blockage of gas and oil pipelines caused by formation of clathrate hydrates is a serious industrial problem.1-3 Clathrate hydrates are crystalline inclusion compounds consisting of water and guest molecules and there are various hydrate structures.3-8 Natural gas, which is mainly composed of methane, forms either CS-I or CS-II clathrate hydrate in the presence of water under high pressure.3 CS-I clathrate hydrate consists of the 512 and 51262 cages and CS-II clathrate hydrate consists of the 512 and 51264 cages. There are two categories of chemical methods to prevent formation of hydrate plugs: thermodynamic inhibitors (TIs) and low dosage hydrate inhibitors (LDHIs). Alcohols, glycols, and salts are TIs that can lower the dissociation temperature of clathrate hydrates as a result of reduction of the chemical potential of water in the liquid phase.3, 9-10 LDHIs are further categorized into kinetic hydrate inhibitors (KHIs) and antiagglomerants (AAs).1-3, 11-12 Both KHIs and AAs bind to hydrate surfaces. KHIs retard formation of clathrate hydrates and AAs prevent hydrate particles from aggregating. Usually, LDHIs are dosed at very low concentrations, 0.01 - 5 wt%, while TIs are dosed at concentrations of 20 - 50 wt%.12 Therefore, LDHIs are economically and environmentally more favorable. Molecular dynamics (MD) simulations have been used to examine adsorption of molecules on hydrate surfaces.13-28 Most KHIs have amide groups. Adsorption of KHIs on hydrate surfaces had been attributed to hydrogen bonding between the amide oxygen and water molecules on the hydrate surface.11-12 However, in a previous study, we found that binding of a standard KHI, polyvinylcaprolactam (PVCap), to the interface between clathrate hydrate and liquid water is not caused by the amide hydrogen bonds but driven by the entropic attraction between cavities at the hydrate surface and the hydrophobic parts of PVCap.29 This entropic attraction originates from the fact that the cavities at the hydrate surface are much larger than those in bulk liquid water.30-



33

Page 4 of 40

The same conclusion was reached by Bertolazzo et al. who examined the mode and strength of

binding of amphiphilic molecules to the hydrate surface using the coarse grained mW water model.24 Adsorption on the hydrate surface through hydrophobic groups has also been found in recent MD simulations of AAs25-27 and type I antifreeze protein.28 Recently, we investigated growth of clathrate hydrate in the presence of PVCap oligomers using MD simulations.34 It was found that the oligomers bind strongly to the growing hydrate surface and the surface becomes concave due to the trapped oligomers. The growth rate of the concave surface is much lower than that of the flat surface. This indicates that the mechanism of KHIs is explained from the Gibbs-Thomson effect likewise that of antifreeze proteins.35-44 Figure 1 shows the hydrate growth process in the presence of KHIs schematically. The clathrate hydrate is growing at a temperature of 𝑇 = 𝑇0𝑚 ―∆𝑇, where 𝑇0𝑚 is the equilibrium dissociation temperature of bulk clathrate hydrate and ∆𝑇 is the supercooling. The curvature of the surface, 1/r, increases as time evolves because of the inhibitor molecules bound to the surface. The 𝛼

dissociation temperature of the curved surface is expressed as 𝑇′𝑚 = 𝑇0𝑚 ― 𝑟 , where  is a positive constant ( ~ 73 K nm for methane hydrate).37, 45 This equation is known as the GibbsThomson equation. The effective supercooling for the surface curved by the KHI molecules, ∆𝑇′ 𝛼

= 𝑇′𝑚 ― 𝑇 = ∆𝑇 ― 𝑟 , decreases with increasing curvature. When the curvature reaches 1/r* =  𝑇/, the effective supercooling ∆𝑇′ becomes zero and the crystal stops growing. Such complete inhibition occurs when the distance between the adsorbed molecules is short enough and the temperature is not very low (the adsorbed molecules can inhibit the crystal growth when the distance between them is shorter than 2𝑟 ∗ = 2𝛼/(𝑇0𝑚 ― 𝑇)).38


Page 5 of 40 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


Figure 1. Schematic of the hydrate growth process in the presence of KHIs. The blue area represents clathrate hydrate and red spheres are inhibitor molecules irreversibly bound to the hydrate/water interface.

The mechanism of KHIs has still not been fully understood. Most KHIs are water soluble polymers.1, 11, 46-56 It is known that monomers of KHIs exhibit little inhibition effect and that the performance of KHIs increases sharply with increasing polymer size, and then decreases gradually under a constant wt% of the inihibitor.1, 11, 55-56 Figure 2 shows the size dependence of the supercooling ability of PVCap for methane hydrate taken from an experimental study of Sloan et al.55 In their experiment, methane hydrate formed at T = 6 K in the absence of PVCap when the temperature was decreased at a constant rate of 2.5 K/h at 206 bar (formation of methane hydrate was detected by a rapid decrease in pressure). The effect of monomers was found to be negligible in concentrations as high as 5 wt%. Methane hydrate formed at T = 19 K when PVCap oligomers of Mn = 900 were dosed at a concentration of 1 wt%, indicating that PVCap of this size inhibits hydrate formation efficiently. The supercooling decreases with increasing polymer size for Mn > 900.



Figure 2. Supercooling ability of PVCap plotted against the size of the polymer. The concentration of PVCap is 1 wt%. The values are taken from an experimental study of Sloan et al.55 The structure of PVCap is also shown. The molecular weight of the monomer unit is 137.

A possible mechanism for the little effect of monomers is that clathrate hydrate can grow by incorporating small monomers in hydrate cages with some distortion of the crystal lattice. There might be a contribution from this mechanism to the size effect. However, it cannot explain why AAs, which also bind to hydrate surfaces, do not inhibit growth of natural gas hydrates although they have long alkyl chains that are too bulky to be included in hydrate cages even with the distortion of the lattice.11 There is another interesting size effect: blends of high and low molecular weight polymers of a KHI exhibit better performance than the KHI with a unimodal molecular weight distribution.1, 11 In fact, the blends of PVCap have been used commercially. AAs are quaternary ammonium or phosphonium cations.11 One or two alkyl chain is very long and the remaining chains are butyl or phenyl groups. It has been reported that quaternary ammonium salts without long chains, such as tetraphenylammonium bromide (TPAB), also do not inhibit formation of natural gas hydrates.11 However, these cations inhibit efficiently growth


Page 6 of 40

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


of clathrate hydrate of tetrahydrofuran (THF).11, 48 This indicates that the performance of a KHI is dependent not only on its size but also on the guest molecule. Several questions arise naturally as to the above observations on the efficiency of KHIs. Why do small oligomers inhibit hydrate growth efficiently more than monomers and large polymers? Why is the performance of KHIs with a bimodal molecular weight distribution better than those with a unimodal distribution? Why do quaternary ammonium salts significantly retard growth of THF hydrate while they cannot inhibit formation of natural gas hydrate? In this paper, we address all these questions and demonstrate that the experimental observations can be explained by a single model which considers the residence time of inhibitors bound to the hydrate surface and the intrinsic growth rates of clathrate hydrates on the basis of the finding that the delay is due to the Gibbs-Thomson effect. The intrinsic growth rates can be taken from experimental results. In contrast, it is difficult to measure the residence time and there are no available experimental data. Therefore, we perform MD simulations of aqueous solutions of the monomer, dimer, tetramer, and octamer of PVCap coexisting with a slab of methane hydrate to calculate the residence time. MD simulations of TPAB are also performed.

Computational details An empty CS-I hydrate structure consisting of 26496 water molecules is generated using the GenIce tool,57 and 1536 methane molecules are inserted in cages of the empty hydrate so that 1/3 of the system becomes fully occupied hydrate. The size of the simulation cell is 9.3  9.8  9.8 nm3. Next, an isothermal isochoric MD simulation is performed for 50 ps at 1000 K with fixing coordinates of atoms in the fully occupied hydrate region to make a hydrate/water coexistence



system. Then, an isothermal isobaric MD simulation is performed for 100 ps at 276 K and 200 bar without any constraint to relax the system. PVCap molecules are inserted in the liquid water region. Several water molecules are removed to avoid overlapping in this process. The numbers of solute molecules are 48, 24, 12, and 6 for MD simulations of the monomer, dimer, tetramer, and octamer, respectively: the concentration of PVCap in the aqueous solution is 2 wt% for all systems. PVCap molecules are placed at z = 0 as shown in Figure 3a. We also perform MD simulations of TPAB. The number of TPAB is set to 24. The cations are placed at z = 0 while the anions are randomly inserted in the solution. The resultant structures are used as initial states of product runs, each of which is as long as 800 ns. We perform three independent MD simulations for each system with different initial atomic velocities that are determined randomly.

Figure 3. Snapshots along an MD simulation. A slab of CS-I methane hydrate is coexistent with an aqueous solution of PVCap. Water molecules and methane molecules are represented by small dots and white spheres, respectively. Six PVCap octamers are placed at z = 0 in the initial configuration. The {110} surface is exposed to the solution because this surface is more stable than any other surface of CS-I clathrate hydrate.58 The images looking down on the hydrate surfaces are given in Supporting Information.


Page 8 of 40

Page 9 of 40 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


A methane molecule is described by a single-site Lennard-Jones (LJ) particle. The LJ parameters of methane are taken from a text book ( = 0.3882 nm and  = 1.139 kJ mol–1).59 The TIP4P/Ice model is used for water.60 This combination of force fields well reproduces the experimental water/methane/hydrate three-phase coexistence temperature for a wide range of pressure.61 The OPLS-AA model is used for PVCap and TPAB.62-63 It was found that PVCap oligomers aggregate in aqueous solutions with the original OPLS charges.34 To avoid the aggregation, we use the point charges determined from an electronic structure calculation.34 All MD simulations are performed using the GROMACS 4.6 package.64-65 The particle mesh Ewald method is employed with a real space cutoff length of 9 Å.66-67 The temperature and pressure are maintained at 276 K and 200 bar using the Nosé-Hoover method68-69 and the Parrinello-Rahman method, respectively.70-71 This thermodynamic condition is similar to that of the experimental study of Sloan et al.55 The hydrate slab dissociates partly in the MD simulations although the temperature of 276 K is much lower than the three-phase coexistence temperature of ~290 K because the concentration of methane in the aqueous solution is zero in the initial state.61 To prevent this partial dissociation, the position of each methane molecule is constrained using a harmonic potential with a force constant of 5000 kJ mol–1 nm–2.

Results and Discussion Figure 3 shows snapshots along an MD simulation of octamers of PVCap. Six octamers are placed in the middle of the liquid region at t = 0. They diffuse from the initial positions as time evolves and bind to the hydrate surface eventually. Figure 4a shows the number density profiles of carbon in PVCap for this trajectory. The number density profiles of methane and water



molecules are also plotted. There are two hydrate/water interfaces at z = –2.5 and 2.5 nm. The local density of PVCap in the vicinity of the interface increases with time during the initial 200 ns of the simulation because of the adsorption. The number density profile at t = 200 ns (magenta) is almost the same as the number density averaged over 200 < t < 800 ns (black dashed) because the bound octamers are never released from the surface. The time-averaged density profile is asymmetric simply because the number of PVCap octamers, 6, is too small. We perform three MD simulations and find that the local density at the left surface is higher than that at the right surface in one trajectory and the opposite result is obtained for the other two trajectories. The profile would be symmetric if we have more trajectories and average over them. Figure 4b shows the corresponding result for an MD simulation of monomers. The monomers also bind to the hydrate surface. We again find that the number density profile of PVCap carbon at t = 200 ns (magenta) is similar to the time-averaged profile (black dashed). This is not because of irreversible adsorption but arises from the fact that the residence time is much shorter for the monomers and the system is almost completely equilibrated at t = 200 ns.


Page 10 of 40

Page 11 of 40 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


Figure 4. Number density profiles of carbon atoms of PVCap at t = 1 ns (green), 10 ns (yellow), 50 ns (red), and 200 ns (magenta) for (a) octamers and (b) monomers. The dashed curves are the number density profiles of water (grey), methane (blue) and PVCap carbon (black) averaged over the time interval between t = 200 ns and 800 ns. The number density profiles of PVCap carbon are scaled by a factor of 8.

We define immobile (solid-like) water molecules and classify whether a lactam ring is trapped in an open hydrate cage using the number of immobile water molecules around the ring. A trajectory is divided into time bins with a width of 2 ns and the following quantity is calculated for every water molecule in each bin 𝛿2𝑖 = 〈{𝐫𝑖(𝑡 + ∆𝑡) ― 𝐫𝑖(𝑡)}2〉 ,


(1)


where ri is the coordinate vector of the oxygen atom of i-th water molecule.34, 72 We set t = 0.1 ns. A water molecule is classified as an immobile molecule when the 𝛿2𝑖 value is smaller than 2.5 Å2. We plot the number of the immobile water molecules, 𝑁immobile, against time in Figure 5a. 𝑁immobile increases with time for 0 < t < 200 ns. This indicates that the hydrate surface is somehow solidified because of the adsorption of PVCap. We define that a lactam ring of PVCap is trapped in an open cage at the surface when the number of immobile water molecules within 0.65 nm from the center of the ring is more than 8 (the threshold of 0.65 nm corresponds to the first minimum of the radial distribution function). Figure 5b shows the total number of the total trapped lactam rings, 𝑁total trapped. As expected, the time evolution of 𝑁trapped is similar to that of

𝑁immobile.

Figure 5. (a) Number of immobile water molecules and (b) number of lactam rings trapped in open cages at the hydrate surface in the monomer (red) and octamer (black) solutions.


Page 12 of 40

Page 13 of 40 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 number of trapped rings for each oligomer, denoted as 𝑁(1) trapped, can be larger for larger oligomers. Figure 6 presents the probability distribution of 𝑁(1) trapped. The probability that two lactam rings in a dimer are trapped is only 26% and 42% of dimers bind to the hydrate surface by trapping of only one ring in one open cage. The remaining 32% of dimers stay in the liquid region. The maximum value of 𝑁(1) trapped for the tetramer is three because it is structurally impossible that all four rings are trapped simultaneously in four different open cages on a flat hydrate surface. Similarly, the maximum 𝑁(1) trapped value is five for the octamer. The mode value is 3. We present a snapshot of an octamer bound to the surface with 𝑁(1) trapped = 3 in Figure S2.

Figure 6. Probability distribution of the number of trapped lactam rings per PVCap molecule. The latter half (400 < t < 800 ns) of three trajectories, where 𝑁total trapped fluctuates around a value as shown in Figure 5b, are used to obtain the probability distribution. The probability distribution of the number of trapped phenyl groups of TPA+ is also shown.

As mentioned above, the residence time is much shorter for the monomer than for the octamer. We estimate the residence time using the bond correlation function given by



𝐶(𝑡) =

〈ℎ(𝑡)ℎ(0)〉 , 〈ℎ〉

(2)

where h(t) is a binary operator.73-74 We define that h = 1 when the distance between the center of a lactam ring and the center of an open cage at the surface is less than 0.293 nm and h = 0 otherwise (the threshold value of 0.293 nm is the inner radius of the 51262 cage which is defined as rinner = rcage - rwater, where rcage is the radius of the cage and rwater is the van der Waals radius of water).3 The center of each cage is determined using a method given in our previous paper.34 The red curve in Figure 7 is the bond correlation function of the monomer of PVCap. The residence time, , is 75 ns when it is defined as the time at which C(t) reaches e–1. We also calculate C(t) of the dimer, tetramer, and octamer and find that much longer simulations are required to determine  for them.

Figure 7. Bond correlation functions of the monomer (red), dimer (yellow), tetramer (green), and octamer (blue) of PVCap. Each correlation function is calculated from the latter half (400 < t < 800 ns) of three trajectories. The horizontal dashed line indicates e–1. The dotted gray curve shows the bond correlation function of TPA+.


Page 14 of 40

Page 15 of 40 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


Instead of estimating the residence time from direct simulations, we calculate it assuming the following relation with the aid of the binding free energy Gbind: 𝜏 ―1 = 𝐴exp( ― 𝐺bind/𝑘𝐵𝑇),

(3)

where kB is the Boltzmann constant and A is a pre-exponential factor. The binding free energy is calculated using the umbrella sampling method. First, we select a PVCap molecule trapped on the hydrate surface in the MD simulations. The selected PVCap molecule is fixed at z = zi using a harmonic potential with a force constant of 1000 kJ mol–1 nm–2. The width of a window is zi+1 – zi = 0.02 nm. We perform an MD simulation for 40 ns for each window. The free energy profile for transferring the selected molecule is obtained from the last 20 ns of the simulations using the weighted histogram analysis method.75-76 Figure 8a shows the free energy profiles of five different monomers. The binding free energy is defined as Gbind = Gbulk – Gmin, where Gbulk is the free energy of the molecule in bulk water and Gmin is the minimum value of the free energy. The binding free energy of the monomer is approximately 10 kJ mol–1. Figure 8b shows the free energy profiles of five octamers. The range of the binding free energies, 20 < Gbind < 50 kJ mol–1, is very wide reflecting the distribution of the probability of 𝑁(1) trapped shown in Figure 6. We plot (1) Gbind against 𝑁(1) trapped in Figure 8c (𝑁trapped is calculated from the MD simulation at the

minimum position of the free energy profile). There is a linear relation between Gbind and –1 𝑁(1) trapped with a slope of 10.1 kJ mol . We denote this slope as the binding free energy per 1 trapped hydrophobic group, 𝐺(1) bind. The value of 10.1 kJ mol is slightly lower than the binding

free energy of the monomer for CS-II methane hydrate of 11.3 kJ mol–1 calculated in our previous study because the 51262 cage of CS-I hydrate is smaller than the 51264 cage of CS-II hydrate.29



Figure 8. Free energy profiles for transferring PVCap along the z axis for (a) monomers and (b) octamers. The distance between the center of mass of PVCap and the center of the hydrate slab is defined as z’. (c) Binding free energy, Gbind, plotted against the number of trapped rings per molecule, 𝑁(1) trapped. The slope of the dashed line, i.e., the binding free energy per trapped hydrophobic group, is 10.1 kJ mol–1 for PVCap. The Gbind values of TPA+ also follow the linear relation (the slope for TPA+ may be different from that for PVCap but it is unclear because of the small number of samples).

We expect from the linear relation that the average of the binding free energy is approximated (1) (1) as 𝐺bind(𝑛) = 𝐺(1) bind𝑁trapped(𝑛), where n is the degree of polymerization and 𝑁trapped(𝑛) is the

average of the number of trapped lactam rings per bound PVCap molecule calculated from the distribution shown in Figure 6. We obtain A = 1.1 ns–1 for the monomer from  = 75 ns (Figure 7) and 𝐺bind(1) = 10.1 kJ mol1. The residence times of the dimer, tetramer, and octamer are 16 Environment ACS Paragon Plus

Page 16 of 40

Page 17 of 40 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


calculated from 𝐺bind(𝑛) assuming that A is independent of n. The average binding free energy and the corresponding residence time are shown in Figure 9.

Figure 9. Binding free energy (left axis) and the corresponding residence time (right axis) plotted against the degree of polymerization, n (open squares). The residence time is calculated from Eq. (3) with A = 1.1 ns–1. The red and gray horizontal dashed lines indicate t* of methane hydrate (60 s) and that of THF hydrate (0.6 s), respectively. The yellow circle shows  of the octamer calculated with A = 0.55 ns–1.

We propose a model to explain the size dependence of KHIs in which the residence time, , plays an important role. The growth process of clathrate hydrate is delayed by the adsorbed molecules because of the Gibbs-Thomson effect. Slowing down is not prominent immediately after the adsorption event because the curvature of the surface is small at this stage (Figure 10b). There must be a time interval between the adsorption event and noticeable slowing down of the crystal growth. We denote this duration as 𝑡 ∗ and the displacement of the front of the growing surface during 𝑡 ∗ as 𝐿 ∗ (Figure 10c). If the binding free energy is not high enough, the molecules would not continue to be trapped during 𝑡 ∗ (Figure 10f). The concave regions caused by the bound molecules flatten very quickly after desorption because curved surfaces are



thermodynamically unstable.34, 37 The same or other molecules would bind to the surface again but they also desorb quickly (Figures 10g and 10h). Therefore, the effective supercooling, ∆𝑇′, never reaches zero and the hydrate crystal grows continuously. Such molecules cannot inhibit hydrate formation efficiently even when their concentration is high. We refer to molecules for which τ < 𝑡 ∗ as “weakly bound molecules” and molecules which satisfy τ > 𝑡 ∗ as “strongly bound molecules”. An important feature is that  increases exponentially with increasing n as shown in Figure 9. This suggests that properties which do not change exponentially with n, such as the pre-exponential factor A, are not important for the classification of molecules into weakly or strongly bound molecules. The yellow circle in Figure 9 is  of the octamer calculated assuming that the factor A for the octamer is half that for the monomer (the self-diffusion coefficient of the octamer, 0.054  10–9 m2 s–1, is about half that of the monomer, 0.126  10–9 m2 s–1). We see that the shift caused by the change in A is indeed fairly small.

Figure 10. (a-d) Inhibition of crystal growth caused by strongly bound molecules. The dashed horizontal lines in panels b-d represent the hydrate/water interface at t = 0. (e-f) Continuous hydrate growth in the presence of weakly bound molecules. Desorbed molecules are represented by gray spheres.


Page 18 of 40

Page 19 of 40 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


We estimate 𝑡 ∗ for the classification of bound molecules. The relation between t* and L* can 𝑡∗

be given by 𝐿 ∗ = ∫0 𝑘(𝑡)𝑑𝑡, where k(t) is the growth rate of the hydrate. Since t* is defined as the time until which the delay of the crystal growth is not noticeable, we simplify the relation as 𝐿 ∗ ~𝑘0𝑡 ∗ , where k0 is the intrinsic growth rate of the hydrate, i.e., the growth rate without KHIs. It was shown that L* is approximately 1 nm for CS-I ethylene oxide (EO) hydrate.34 We assume that L* is 1 nm for methane hydrate ignoring the difference in guest species. The experimental k0 value is roughly 10 mm min–1 for methane hydrate at high supercooling.77 This value is, however, not suitable for the present purpose because it was obtained for macroscopic crystals. The onset temperature measured in the experimental study is dominated by the effect of KHI on invisibly small nuclei which grow slowly because of the Gibbs-Thomson effect.55 We assume tentatively that the growth rate of the invisibly small particles is one order of magnitude lower than that of macroscopic crystals, i.e., k0 = 1 mm min–1. We obtain t* = 60 s from L*= 1 nm and k0 = 1 mm min–1. In Figure 9, the red horizontal dashed line indicates t* = 60 s. The residence time  is much shorter than t* at n = 1, indicating that monomers of PVCap are weakly bound to the hydrate surface and cannot inhibit growth of methane hydrate. It is possible to convert the probability distribution of the number of trapped lactam rings shown in Figure 6 to the probability distribution of the residence time  using 𝐺(1) bind and Eq. 3. Figure 11a shows the probability distribution of  for the octamer (blue bars). The residence time of the octamer ranges from ~100 ns to ~109 ns (1 second). The vertical dashed line indicates t* =



Page 20 of 40

60 s. To estimate the ratio of the number of strongly bound molecules to the total number of bound molecules, xstrong, we approximate the probability distribution by a Gaussian function,

𝑃(log 𝜏) =

1

exp 2𝜋𝜎2

(

(log 𝜏 ― log 𝜏)2 2𝜎2

)

,

(4)

where log 𝜏 and  are the average and the standard deviation, respectively. The ratio xstrong is found to be 0.64 for the octamer. As shown in Figure 6, both log 𝜏 and  increase with increasing n. Figure 11b presents the n dependence of xstrong obtained under the assumption that log 𝜏 and  increase linearly with n. The ratio xstrong increases sharply and then approaches 1 gradually with increasing n.

Figure 11. (a) Probability distribution of the residence time of PVCap. The dashed red line indicates t* = 60 s. The inset shows the corresponding result for TPA+. The red and gray lines


Page 21 of 40 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


indicate t* of methane hydrate (60 s) and THF hydrate (0.6 s), respectively. (b) Ratio of the number of strongly bound PVCap molecules to the total number of bound molecules, xstrong, plotted against the degree of polymerization.

Figure 11b suggests that the increase in n does not contribute to the improvement of the performance of the KHI for n > 20 where xstrong is almost 1. Rather, the performance is expected to decrease with increasing n due to a different mechanism, as explained schematically in Figure 12. In this figure, a black bar represents a monomer of PVCap bound to the hydrate surface. The hydrate grows as time evolves and the region around each molecule is curved. The growth rate is delayed only in this curved region, which is represented by a light blue eclipse. We denote the area curved by a monomer as a1. The corresponding area for a dimer is not 2a1 but 2a1  s because of overlapping of two eclipses. The overlapping area is wider for larger polymers when the number of monomer units is the same. Therefore, large polymers are less effective than small oligomers.

Figure 12. Schematic for the mechanism which dominates the size dependence in the large n region. Each black bar represents a monomer of PVCap bound to the hydrate surface and the



Page 22 of 40

surrounding light blue area is the surface region curved by the monomer. The deep blue region is overlapping of light blue areas.

The size dependence can be represented by a simple equation. When the size of all KHI polymers is n, the surface area covered by them is given by 𝑎(𝑛) = 𝑁m𝑎1 ― 𝑁p(𝑛 ― 1)s,

(5)

where Np is the number of strongly bound polymers and Nm is the number of monomer units in them. Because Np = Nm/n, we obtain

(

𝑎(𝑛) = 𝑁m𝑎1 1 ―

)

𝑠 (𝑛 ― 1) . 𝑎1 𝑛

(6)

Nm is given by 𝑁m = 𝑁total m 𝑥bound(𝑛)𝑥strong(𝑛),

(7)

where 𝑁total is the total number of monomer units which is a constant when the performance is m compared at the same weight percent of the KHI. xbound(n) is the ratio of the number of bound molecules to that of all molecules. If the supercooling ability is simply proportional to a(n), we have

(

∆𝑇(𝑛) ― ∆𝑇pure ∝ 𝑥bound(𝑛)𝑥strong(𝑛) 1 ―

)

𝑠 (𝑛 ― 1) , 𝑎1 𝑛

(8)

where Tpure is the supercooling without additives. We fit the experimental data of Sloan et al. to Eq. (8) assuming that n is identified as the experimental weigh-average degree of polymerization.


Page 23 of 40 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 polydispersity index (PDI = Mw/Mn) of PVCap is 2.78 As shown in Figure 6, xbound is 0.59 for the monomer and 0.68 for the dimer. For larger polymers, 𝑥bound is approximated as 1. The best fit is shown as the black solid curve in Figure 13. There is excellent agreement between the model and the experimental data. The ratio s/a1 obtained from the fitting is 0.96. This fairly wide overlapping area is consistent with the observation in our previous study that a PVCap molecule concaves a wide range of the surface (the area of overlap, s, increases with increasing a1).34

Figure 13. Fit of Eq. (8) to the experimental data of Sloan et al.55 The experimental supercooling is shown by open squares. The black solid curve is the best fit of Eq. 8 when t* is set to 60 s. The dashed curves are obtained from the fitting with different t* values.

Figure 13 also shows the fitting curves obtained with several different t* values. The peak position, i.e., the optimal PVCap size, shifts toward larger n with increasing t*. However, even when t* is changed from 6 s to 600 s, the optimal size changes only from n = 7 to n = 13. The performance of the monomer is somewhat overestimated when t* is 6 s and the peak height is underestimated when t* is 600 s, implying that t* = 60 s is a good guess although it is based on the crude assumptions.



The present model can explain why blends of high and low molecular weight polymers show very good performance.1, 11 Tiny hydrate particles grow slowly because of the Gibbs-Thomson effect, meaning that t* is long for them. Our model suggests that large polymers can inhibit the growth of such particles because of large  values. In contrast, small polymers can inhibit growth of relatively large (but still invisibly small) particles efficiently. The performance of the blends is good probably because they can inhibit growth of both small and large hydrate particles. It has been reported that amphiphilic molecules promote formation of clathrate hydrates.79-81 A model proposed by Betrolazzo et al. suggests that the nucleation rate of clathrate hydrate increases due to the reduction of the hydrate/water interfacial free energy caused by sodium dodecyl sulfate bound to the surface of nascent nuclei.24 PVCap may also decrease the interfacial free energy because it is amphiphilic in character. Adsorbents must desorb from the surface so that the hydrate particles grow and reach the critical size. Therefore, strongly bound large polymers cannot act as nucleation promoters. In contrast, the hydrate particles can grow in the presence of monomers or very small oligomers which cannot bind to the surface strongly as schematically shown in Figures 10e-10h. The supercooling ability may become somewhat lower than that shown in Figure 13 for the small n region when the promotion effect is considered. We also perform MD simulations of the aqueous TPAB solution coexisting with methane hydrate. As shown in Figure 6, typically two of four phenyl groups are trapped in open hydrate cages (see also Figure S2). Figure 7 shows that the decay of the bond correlation function of TPA+ is faster than that of the PVCap tetramer. Correspondingly, as shown in Figure 9, the binding free energy of TPA+ is lower than that of the PVCap tetramer. The average residence time is 4 s for TPA+ if its pre-exponential factor A is the same as that for the PVCap monomer.


Page 24 of 40

Page 25 of 40 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


Experimental studies reported that TPA+ cannot inhibit formation of natural gas hydrates while it delays growth of THF hydrate efficiently.11,

48

This can be explained from the large

difference in t* between gas hydrates and THF hydrate. THF is completely miscible with liquid water. Crystal growth of THF hydrate is much faster than that of gas hydrates because of the absence of the slow transfer process of guest molecules from the gas phase to the growing hydrate surface in water. The experimental growth rate of macroscopic THF hydrate crystals is one order of magnitude higher than that of methane hydrate when they are compared at the same supercooling.77, 82 In addition, macroscopic crystals are used in experimental studies of KHIs for THF hydrate while the experiments for gas hydrates measure the effects of KHIs on microscopic hydrate crystals. Therefore, the t* value for THF hydrate is much shorter than that of methane hydrate. We assume that t* is 100 times shorter for THF hydrate. The gray horizontal dashed line in Figure 9 shows t* of THF hydrate. The average  value of TPA+, 4 s (gray square), is between the gray and red dashed lines, indicating that TPA+ is a KHI for THF hydrate while not for methane hydrate. The inset of Figure 11a shows the probability distribution of  for TPA+. It is found that 74% of TPA+ cations are categorized as strongly bound molecules for THF hydrate while the ratio becomes only 10% of for methane hydrate. The  value of TPA+ for CS-II THF hydrate might be different from that for CS-I methane hydrate. As mentioned above, the binding free energy of the PVCap monomer is higher for CS-II hydrate because of the large 51264 cage of CS-II. However, the difference is only ~1 kJ mol1.29 In addition, the number density of the small 512 cages is higher for CS-II. Therefore, it is expected that  is not so different from the present value when it is estimated with CS-II hydrate.



The {100} surface of CS-II hydrate was exposed to the aqueous solution in the MD simulation of ref. 29 while TPA+ adsorbs on the {111} surface which is the slowest growing surface of CSII THF hydrate. This is also not problematic because the surface affinity is insensitive to the crystallographic surface. PVCap oligomers can adsorb on various surfaces of clathrate hydrate without specific preferential surfaces as revealed in ref. 34. The matching of the spacing between hydrophobic groups of a KHI and the periodicity of the surface structure is not required for binding because of the fairly strong attractive interaction between each open cage and hydrophobic group. This behavior is contrast to that of antifreeze proteins (AFPs) each of which prefers one or a few specific surfaces.36, 83 An MD study of Bellucci et al. showed that the residence time of an AA, n-dodecyl-tri(nbutyl) ammonium cation, is ~ 0.9 s at 277 K when two butyl groups are trapped in open hydrate cages at the {111} surface of methane-propane CS-II hydrate.26 This value is much shorter than t* for methane hydrate, 60 s, suggesting that this AA is not a KHI for this clathrate hydrate. The residence time of tetrabutylammonium cation, TBA+, is expected to be almost the same as that of this AA. The value of 0.9 s is shorter than  of TPA+ calculated in this study, 4 s, but somewhat longer than t* for THF hydrate of 0.6 s. This is consistent with the experimental finding that TBA+ is a KHI for THF hydrate but it is less efficient than TPA+.48

Conclusions We have investigated the mechanism of KHIs using MD simulations. A model is proposed to explain the size (degree of polymerization) dependence of KHIs. It is assumed that a molecule can act as a kinetic inhibitor only when the residence time of the molecule at the hydrate surface


Page 26 of 40

Page 27 of 40 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


exceeds a threshold value determined by the intrinsic growth rate. The ratio of the number of KHI molecules whose residence time is longer than the threshold increases with size because the number of trapped hydrophobic groups per molecule increases. However, the increase in the polymer size also results in the reduction of the coverage area per monomer unit. As a result of these two opposing effects, the performance of a KHI is maximized at a certain polymer size. This model, together with the residence time estimated from free energy calculations and the intrinsic growth rate taken from experimental data, reproduces the size dependence of the supercooling ability of PVCap for methane hydrate very well. Our model also suggests that small polymers inhibit growth of large particles efficiently while large polymers can inhibit growth of small hydrate particles and thus the performance of blends of high and low molecular weight polymers is better than that of the KHI with a unimodal molecular weight distribution. MD simulations are also performed for a quaternary ammonium cation, TPA+. This cation cannot inhibit growth of methane hydrate because of its short residence time. However, TPA+ is a good kinetic inhibitor for THF hydrate because of the very high growth rate of THF hydrate arising from the high solubility of the guest in water. THF hydrate has often been used as an analogue to natural gas hydrates because experiments can be performed under ambient pressure.11 The present study suggests that growth of THF hydrate is delayed even with weakly bound molecules that cannot inhibit gas hydrate formation. Therefore, the assessment of KHIs with THF hydrate may not be a useful way to find efficient KHIs for natural gas hydrates. In this study, all guest molecules are fixed in hydrate cages and the hydrate surface is completely flat during the MD simulations. Taking surface roughness and curvature into consideration may affect the residence time. However, this effect is expected to be insignificant



Page 28 of 40

because of the cancelation between strong binding to open cages at concave regions of the surface and weak binding to convex regions. The delay of crystal growth caused by the Gibbs-Thomson effect has been observed for various crystals such as ice and calcite.39-40 In many cases, irreversible adsorption is assumed to model the effect of adsorbent. In contrast, our model takes account of the finite residence time and its size dependence as well as the intrinsic growth rate of the crystal. This treatment may provide new insight into the mechanism of growth inhibition for crystals other than clathrate hydrates. It is known that the thermal hysteresis of AFPs is higher for larger proteins although the range of size examined has been limited, and there are attempts to synthesize high performance AFPs by polymerization (or oligomerization).39,

84-87

For example, the thermal

hysteresis of the dimer of AFP III is twice that of the monomer and the trimer and tetramer exhibit even higher activity.85-86 Our model suggests that the improvement of the activity due to polymerization levels off at a certain size and further polymerization causes reduction of the efficiency if the mechanism of the AFP is similar to that of KHIs. Such a behavior was found for -helical TmAFP.84

Supporting Information The supporting information contains images looking down to the hydrate surfaces trapping KHI molecules and enlarged snapshots of a PVCap octamer and TPA+ on the surface.

Author Information 28 Environment ACS Paragon Plus

Page 29 of 40 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


Corresponding Author Phone: +81-(0)86-251-7769 E-mail: [email protected]

Acknowledgments The present work was supported by a grant of MORINO FOUNDATION FOR MOLECULAR SCIENCE and MEXT as "Priority Issue on Post-Kcomputer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use) using computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp180204). Calculations were also performed on the computers at Research Center for Computational Science, Okazaki, Japan.

References 1.

Kelland, M. A., Production Chemicals for the Oil and Gas Industry; CPC Press: Boca Raton, 2014.

2.

Creek, J. L., Efficient Hydrate Plug Prevention. Energy Fuels 2012, 26, 4112-4116.

3.

Sloan, E. D.; Koh, C. A., Clathrate Hydrates of Natural Gases; CPC Press: Boca Raton, 2008.

4.

Ripmeester, J. A.; Tse, J. S.; Ratcliffe, C. I.; Powell, B. M., A New Clathrate Hydrate Structure. Nature 1987, 325, 135-136.



5.

Jeffrey, G. A., Hydrate Includion Compounds. J. Inclusion Phenom. 1984, 1, 211-222.

6.

Huang, Y.; Zhu, C.; Wang, L.; Cao, X.; Su, Y.; Jiang, X.; Meng, S.; Zhao, J.; Zeng, X. C., A New Phase Diagram of Water under Negative Pressure: The Rise of the LowestDensity Clathrate S-III. Sci. Adv. 2016, 2, e1501010.

7.

Huang, Y.; Zhu, C.; Wang, L.; Zhao, J.; Zeng, X. C., Prediction of a New Ice Clathrate with Record Low Density: A Potential Candidate as Ice XIX in Guest-Free Form. Chem. Phys. Lett. 2017, 671, 186-191.

8.

Matsui, T.; Hirata, M.; Yagasaki, T.; Matsumoto, M.; Tanaka, H., Communication: Hypothetical Ultralow-Density Ice Polymorphs. J. Chem. Phys. 2017, 147, 091101.

9.

Yagasaki, T.; Matsumoto, M.; Andoh, Y.; Okazaki, S.; Tanaka, H., Dissociation of Methane Hydrate in Aqueous Nacl Solutions. J. Phys. Chem. B 2014, 118, 11797-11804.

10.

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, 17, 32347-32357.

11.

Kelland, M. A., History of the Development of Low Dosage Hydrate Inhibitors. Energy Fuels 2006, 20, 825-847.

12.

Perrin, A.; Musa, O. M.; Steed, J. W., The Chemistry of Low Dosage Clathrate Hydrate Inhibitors. Chem. Soc. Rev. 2013, 42, 1996-2015.

13.

Carver, T. J.; Drew, M. G. B.; Rodger, P. M., Inhibition of Crystal Growth in Methane Hydrate. J. Chem. Soc., Faraday Trans. 1995, 91, 3449-3460.

14.

Carver, T. J.; Drew, M. G. B.; Rodger, P. M., Characterisation of the {111} Growth Planes of a Type II Gas Hydrate and Study of the Mechanism of Kinetic Inhibition by Poly(Vinylpyrrolidone). J. Chem. Soc., Faraday Trans. 1996, 92, 5029-5033.


Page 30 of 40

Page 31 of 40 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


15.

Carver, T. J.; Drew, M. G. B.; Rodger, P. M., Configuration-Biased Monte Carlo Simulations of Poly(Vinylpyrrolidone) at a Gas Hydrate Crystal Surface. Ann. N.Y. Acad. Sci. 2000, 912, 658-668.

16.

Storr, M. T.; Taylor, P. C.; Monfort, J.-P.; Rodger, P. M., Kinetic Inhibitor of Hydrate Crystallization. J. Am. Chem. Soc. 2004, 126, 1569-1576.

17.

Freer, E. M.; Sloan, E. D., An Engineering Approach to Kinetic Inhibitor Design Using Molecular Dynamics Simulations. Ann. N.Y. Acad. Sci. 2000, 912, 651-657.

18.

Davenport, J. R.; Musa, O. M.; Paterson, M. J.; Piepenbrock, M.-O. M.; Fucke, K.; Steed, J. W., A Simple Chemical Model for Clathrate Hydrate Inhibition by Polyvinylcaprolactam. Chem. Commun. 2011, 47, 9891-9893.

19.

Hawtin, R. W.; Rodger, P. M., Polydispersity in Oligomeric Low Dosage Gas Hydrate Inhibitors. J. Mater. Chem. 2006, 16, 1934-1942.

20.

Moon, C.; Hawtin, R. W.; Rodger, P. M., Nucleation and Control of Clathrate Hydrates: Insights from Simulation. Faraday Discuss. 2007, 136, 367-382.

21.

Kvamme, B.; Kuznetsova, T.; Aasoldsen, K., Molecular Dynamics Simulations for Selection of Kinetic Hydrate Inhibitors. J. Mol. Graphics Modell. 2005, 23, 524-536.

22.

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, 127, 17852-17862.

23.

Gómez Gualdrón, D. A.; Balbuena, P. B., Classical Molecular Dynamics of Clathrate−Methane−Water−Kinetic Inhibitor Composite Systems. J. Phys. Chem. C 2007, 111, 15554-15564.



24.

Bertolazzo, A. A.; Naullage, P. M.; Peters, B.; Molinero, V., The Clathrate-Water Interface Is Oleophilic. J. Phys. Chem. Lett. 2018, 9, 3224-3231.

25.

Jiménez-Ángeles, F.; Firoozabadi, A., Hydrophobic Hydration and the Effect of NaCl Salt in the Adsorption of Hydrocarbons and Surfactants on Clathrate Hydrates. ACS Cent. Sci. 2018, 4, 820-831.

26.

Bellucci, M. A.; Walsh, M. R.; Trout, B. L., Molecular Dynamics Analysis of AntiAgglomerant Surface Adsorption in Natural Gas Hydrates. J. Phys. Chem. C 2018, 122, 2673-2683.

27.

Mehrabian, H.; Bellucci, M. A.; Walsh, M. R.; Trout, B. L., Effect of Salt on Antiagglomerant Surface Adsorption in Natural Gas Hydrates. J. Phys. Chem. C 2018, 122, 12839-12849.

28.

Bagherzadeh, S. A.; Alavi, S.; Ripmeester, J. A.; Englezos, P., Why Ice-Binding Type I Antifreeze Protein Acts as a Gas Hydrate Crystal Inhibitor. Phys. Chem. Chem. Phys. 2015, 17, 9984-9990.

29.

Yagasaki, T.; Matsumoto, M.; Tanaka, H., Adsorption Mechanism of Inhibitor and Guest Molecules on the Surface of Gas Hydrates. J. Am. Chem. Soc. 2015, 137, 12079-12085.

30.

Widom, B., Potential-Distribution Theory and the Statistical Mechanics of Fluids. J. Phys. Chem. 1982, 86, 869-872.

31.

Widom, B.; Bhimalapuram, P.; Koga, K., The Hydrophobic Effect. Phys. Chem. Chem. Phys. 2003, 5, 3085-3093.

32.

Pohorille, A.; Pratt, L. R., Cavities in Molecular Liquids and the Theory of Hydrophobic Solubilities. J. Am. Chem. Soc. 1990, 112, 5066-5074.


Page 32 of 40

Page 33 of 40 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


33.

Chandler, D., Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640-647.

34.

Yagasaki, T.; Matsumoto, M.; Tanaka, H., Adsorption of Kinetic Hydrate Inhibitors on Growing Surfaces: A Molecular Dynamics Study. J. Phys. Chem. B 2018, 122, 33963406.

35.

Raymond, J. A.; DeVries, A. L., Adsorption Inhibition as a Mechanism of Freezing Resistance in Polar Fishes. Proc. Natl. Acad. Sci. U.S.A. 1977, 74, 2589-2593.

36.

Nada, H.; Furukawa, Y., Antifreeze Proteins: Computer Simulation Studies on the Mechanism of Ice Growth Inhibition. Polym. J. 2012, 44, 690-698.

37.

Naullage, P. M.; Qiu, Y.; Molinero, V., What Controls the Limit of Supercooling and Superheating of Pinned Ice Surfaces? J. Phys. Chem. Lett. 2018, 9, 1712-1720.

38.

Kristiansen, E.; Zachariassen, K. E., The Mechanism by Which Fish Antifreeze Proteins Cause Thermal Hysteresis. Cryobiology 2005, 51, 262-280.

39.

Shtukenberg, A. G.; Ward, M. D.; Kahr, B., Crystal Growth with Macromolecular Additives. Chem. Rev. 2017, 117, 14042-14090.

40.

Bar Dolev, M.; Braslavsky, I.; Davies, P. L., Ice-Binding Proteins and Their Function. Annu. Rev. Biochem 2016, 85, 515-542.

41.

Mochizuki, K.; Molinero, V., Antifreeze Glycoproteins Bind Reversibly to Ice via Hydrophobic Groups. J. Am. Chem. Soc. 2018, 140, 4803-4811.

42.

Knight, C. A.; Devries, A. L., Melting Inhibition and Superheating of Ice by an Antifreeze Glycopeptide. Science 1989, 245, 505-507.

43.

Knight, C. A.; Cheng, C. C.; DeVries, A. L., Adsorption of Alpha-Helical Antifreeze Peptides on Specific Ice Crystal Surface Planes. Biophys. J. 1991, 59, 409-418.



44.

Celik, Y.; Graham, L. A.; Mok, Y. F.; Bar, M.; Davies, P. L.; Braslavsky, I., Superheating of Ice Crystals in Antifreeze Protein Solutions. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 5423-5428.

45.

Jacobson, L. C.; Molinero, V., Can Amorphous Nuclei Grow Crystalline Clathrates? The Size and Crystallinity of Critical Clathrate Nuclei. J. Am. Chem. Soc. 2011, 133, 64586463.

46.

O’Reilly, R.; Ieong, N. S.; Chua, P. C.; Kelland, M. A., Missing Poly(N-Vinyl Lactam) Kinetic Hydrate Inhibitor: High-Pressure Kinetic Hydrate Inhibition of Structure II Gas Hydrates with Poly(N-Vinyl Piperidone) and Other Poly(N-Vinyl Lactam) Homopolymers. Energy Fuels 2011, 25, 4595-4599.

47.

Kelland, M. A.; Kvæstad, A. H.; Astad, E. L., Tetrahydrofuran Hydrate Crystal Growth Inhibition by Trialkylamine Oxides and Synergism with the Gas Kinetic Hydrate Inhibitor Poly(N-Vinyl Caprolactam). Energy Fuels 2012, 26, 4454-4464.

48.

Chua, P. C.; Kelland, M. A., Tetra(Iso-Hexyl)Ammonium Bromide—the Most Powerful Quaternary Ammonium-Based Tetrahydrofuran Crystal Growth Inhibitor and Synergist with Polyvinylcaprolactam Kinetic Gas Hydrate Inhibitor. Energy Fuels 2012, 26, 11601168.

49.

Reyes, F. T.; Malins, E. L.; Becer, C. R.; Kelland, M. A., Non-Amide Kinetic Hydrate Inhibitors: Performance of a Series of Polymers of Isopropenyloxazoline on Structure II Gas Hydrates. Energy Fuels 2013, 27, 3154-3160.

50.

Reyes, F. T.; Guo, L.; Hedgepeth, J. W.; Zhang, D.; Kelland, M. A., First Investigation of the Kinetic Hydrate Inhibitor Performance of Poly(N-Alkylglycine)s. Energy Fuels 2014, 28, 6889-6896.


Page 34 of 40

Page 35 of 40 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


51.

Reyes, F. T.; Kelland, M. A.; Kumar, N.; Jia, L., First Investigation of the Kinetic Hydrate Inhibition of a Series of Poly(Β-Peptoid)S on Structure II Gas Hydrate, Including the Comparison of Block and Random Copolymers. Energy Fuels 2015, 29, 695-701.

52.

Kelland, M. A.; Abrahamsen, E.; Ajiro, H.; Akashi, M., Kinetic Hydrate Inhibition with N-Alkyl-N-Vinylformamide Polymers: Comparison of Polymers to N-Propyl and Isopropyl Groups. Energy Fuels 2015, 29, 4941-4946.

53.

Sa, J. H.; Kwak, G. H.; Lee, B. R.; Park, D. H.; Han, K.; Lee, K. H., Hydrophobic Amino Acids as a New Class of Kinetic Inhibitors for Gas Hydrate Formation. Sci. Rep. 2013, 3, 2428.

54.

Qin, H.-B.; Sun, Z.-F.; Wang, X.-Q.; Yang, J.-L.; Sun, C.-Y.; Liu, B.; Yang, L.-Y.; Chen, G.-J., Synthesis and Evaluation of Two New Kinetic Hydrate Inhibitors. Energy Fuels 2015, 29, 7135-7141.

55.

Sloan, E. D.; Subramanian, S.; Matthews, P. N.; Lederhos, J. P.; Khokhar, A. A., Quantifying Hydrate Formation and Kinetic Inhibition. Ind. Eng. Chem. Res. 1998, 37, 3124-3132.

56.

Magnusson, C. D.; Kelland, M. A., Nonpolymeric Kinetic Hydrate Inhibitors: Alkylated Ethyleneamine Oxides. Energy Fuels 2015, 29, 6347-6354.

57.

Matsumoto, M.; Yagasaki, T.; Tanaka, H., GenIce: Hydrogen-Disordered Ice Generator. J. Comput. Chem. 2018, 39, 61-64.

58.

Smelik, E.; King, H. E., Crystal-Growth Studies of Natural Gas Clathrate Hydrates Using a Pressurized Optical Cell. Am. Mineral. 1997, 82, 88-98.



59.

Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B., Molecular Theory of Gases and Liquids; Wiley: New York, 1954.

60.

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, 122, 234511.

61.

Tanaka, H.; Yagasaki, T.; Matsumoto, M., On the Thermodynamic Stability of Clathrate Hydrates Vi: Complete Phase Diagram. J. Phys. Chem. B 2018, 122, 297-308.

62.

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.

63.

Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L., Evaluation and Reparametrization of the Opls-Aa Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001, 105, 6474-6487.

64.

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, 4, 435-447.

65.

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, 26, 1701-1718.

66.

Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald: An N⋅Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092.

67.

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G., A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577-8593.

68.

Nosé, S., A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255-268.


Page 36 of 40

Page 37 of 40 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


69.

Hoover, W. G., Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695-1697.

70.

Nosé, S.; Klein, M. L., Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 1983, 50, 1055-1076.

71.

Parrinello, M.; Rahman, A., Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182-7190.

72.

Vatamanu, J.; Kusalik, P. G., Molecular Dynamics Methodology to Investigate SteadyState Heterogeneous Crystal Growth. J. Chem. Phys. 2007, 126, 124703.

73.

Rapaport, D. C., Hydrogen Bonds in Water Network Organization and Lifetimes. Mol. Phys. 1983, 50, 1151-1162.

74.

Luzar, A.; Chandler, D., Hydrogen-Bond Kinetics in Liquid Water. Nature 1996, 379, 55-57.

75.

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, 13, 1011-1021.

76.

Hub, J. S.; de Groot, B. L.; van der Spoel, D., G_Wham—a Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. J. Chem. Theory Comput. 2010, 6, 3713-3720.

77.

Tanaka, R.; Sakemoto, R.; Ohmura, R., Crystal Growth of Clathrate Hydrates Formed at the Interface of Liquid Water and Gaseous Methane, Ethane, or Propane: Variations in Crystal Morphology. Cryst. Growth Des. 2009, 9, 2529-2536.



78.

O'Reilly, R.; Ieong, N. S.; Chua, P. C.; Kelland, M. A., Crystal Growth Inhibition of Tetrahydrofuran Hydrate with Poly(N-Vinyl Piperidone) and Other Poly(N-Vinyl Lactam) Homopolymers. Chem. Eng. Sci. 2011, 66, 6555-6560.

79.

Lo, C.; Zhang, J. S.; Somasundaran, P.; Lu, S.; Couzis, A.; Lee, J. W., Adsorption of Surfactants on Two Different Hydrates. Langmuir 2008, 24, 12723-12726.

80.

Zhang, J. S.; Lee, S.; Lee, J. W., Kinetics of Methane Hydrate Formation from SDS Solution. Ind. Eng. Chem. Res. 2007, 46, 6353-6359.

81.

Karaaslan, U.; Parlaktuna, M., Surfactants as Hydrate Promoters? Energy Fuels 2000, 14, 1103-1107.

82.

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, 912, 533-543.

83.

Olijve, L. L. C.; Meister, K.; DeVries, A. L.; Duman, J. G.; Guo, S.; Bakker, H. J.; Voets, I. K., Blocking Rapid Ice Crystal Growth through Nonbasal Plane Adsorption of Antifreeze Proteins. Proc. Natl. Acad. Sci. U.S.A. 2016, 113, 3740-3745.

84.

Marshall, C. B.; Daley, M. E.; Sykes, B. D.; Davies, P. L., Enhancing the Activity of a ΒHelical Antifreeze Protein by the Engineered Addition of Coils. Biochemistry 2004, 43, 11637-11646.

85.

Baardsnes, J.; Kuiper, M. J.; Davies, P. L., Antifreeze Protein Dimer: When Two IceBinding Faces Are Better Than One. J. Biol. Chem. 2003, 278, 38942-38947.

86.

Nishimiya, Y.; Ohgiya, S.; Tsuda, S., Artificial Multimers of the Type III Antifreeze Protein. Effects on Thermal Hysteresis and Ice Crystal Morphology. J. Biol. Chem. 2003, 278, 32307-32312.


Page 38 of 40

Page 39 of 40 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


87.

Stevens, C. A.; Drori, R.; Zalis, S.; Braslavsky, I.; Davies, P. L., Dendrimer-Linked Antifreeze Proteins Have Superior Activity and Thermal Recovery. Bioconjug. Chem. 2015, 26, 1908-1915.



TOC Graphic


Page 40 of 40

Molecular Dynamics Study of Kinetic Hydrate ... - ACS Publications

Recommend Documents