Molecular-Dynamics Study of Propane-Hydrate Dissociation: Non

6 days ago - Non-equilibrium molecular dynamics (NEMD) simulations have been performed to investigate both thermal- and electric-field- driven break-u...
0 downloads 9 Views 7MB Size
Subscriber access provided by UNIV OF MISSOURI ST LOUIS

C: Physical Processes in Nanomaterials and Nanostructures

Molecular-Dynamics Study of Propane-Hydrate Dissociation: Non-Equilibrium Analysis in Externally-Applied Electric Fields Mohammad Reza Ghaani, and Niall Joseph English J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b12238 • Publication Date (Web): 19 Mar 2018 Downloaded from http://pubs.acs.org on March 21, 2018

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 37

The Journal of Physical Chemistry 1

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

Molecular-Dynamics Study of Propane-Hydrate Dissociation: NonEquilibrium Analysis in Externally-Applied Electric Fields Mohammad Reza Ghaania and Niall J. English *a a) School of Chemical and Bioprocess Engineering, University College Dublin, Belfield, Dublin 4, Ireland. Abstract Non-equilibrium molecular dynamics (NEMD) simulations have been performed to investigate both thermal- and electric-field- driven break-up of planar propane hydrate interfaces with liquid water at 250-300 K and in the 0-0.7 V/nm field-intensity range. The melting temperatures of each interface were estimated, and dissociation rates were observed to be strongly dependent on temperature, with higher dissociation rates at larger over-temperatures vis-à-vis melting. It was found that externally-applied electric fields below a certain intensity threshold do not lead to any marked structural distortion or dissociation effect on pre-existing bulk clathrates. However, field strengths higher than 0.7 V/nm led to statistically-significant differences in the observed initial dissociation temperature and rates. The activation energy in constant electric field was calculated based on the Arrhenius equation. The parameters of this equation, both in terms of kinetic and thermodynamics components (A and Ea), change significantly, accelerating the dissociation process.

Corresponding author: * [email protected], Tel: +353-1-7161646; Fax: +353-1-7161177.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 2 of 37 2

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 Clathrate hydrates are inclusion compounds, wherein guest gas molecules, like hydrogen, are trapped in a host water-lattice framework in a variety of polymorphs 1,2, with three more commonly-known hydrate structure types: (s)I, II and H. For instance, the sII clathrate comprises ‘large’ cages of type 51264 (i.e., having twelve pentagonal faces and six hexagonal ones) and small-cage 512 dodecahedra, each composed of five water pentamers, with each cage holding various numbers of guest molecules 1,2. Propane clathrate hydrate deserves further study, not least because of their high potential as a lower-pressure hosting medium for other gases. Because propane hydrate’s sII structure is stable near ambient conditions, neither cryogenic temperatures nor high pressures are necessary. The sII structure has a unit cell made up of 136 water molecules forming 16 dodecahedral 512 cages and 8 hexakaidecahedral 51264 cages. Propane only fits into the 51264 cages leaving the 512 cages available to store another gas, e.g., H2, for economic gas storage 3. 512 dodecahedral cages are made up of 20 water molecules forming 12 pentagonal faces; 51264 cages are made from 28 water molecules forming 12 pentagonal faces and 4 hexagonal faces. Molecular simulation has proved a rather invaluable tool in revealing much about equilibrium properties of clathrate hydrates 2, including structure scattering’ picture of phonons and guest motion in cavities characteristics

12,13

4–8

, a ‘resonant-

9–11

, and dynamical and energetic properties

, hydrogen-bonding

7,14

. Such molecular-

simulations studies on dynamical properties have often focused on vibrational properties of the host lattice and on the cavities, with ‘overlapping’ between these acoustic and optic modes, often interpreted in the context of resonant scattering

15

.

Such an understanding of equilibrium dynamical properties and phonon scattering ACS Paragon Plus Environment

Page 3 of 37

The Journal of Physical Chemistry 3

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

offered by molecular dynamics (MD) has indeed led to progress in recent years in enhancing our understanding of thermal-conduction processes in clathrate hydrates 13,14,16,17

, together with “hopping”-mediated guest-diffusion processes of interest of

energy- and gas- storage applications 18–22. However, as ref. 2 outlines in the case for molecular simulation, the application of external electric fields to gas hydrates is of much technological interest, as well as in terms of elucidation of fundamental mechanistic understanding underpinning any such field effects. For instance, for inhibition of hydrate formation in natural-gas pipelines, aside from the use of either thermodynamic (e.g., methanol or ethylene glycol) or kinetic inhibitor additives injected into the gas stream, or localised application of higher temperature 1,23, external electric and electromagnetic (e/m) fields may serve to disrupt already-formed hydrates (e.g., as pipeline plugs or possibly in situ ‘natural’ hydrates) 1,2. This is certainly more attractive than downstream ‘cleaning’ - extraction of inhibitors – at great operational and economic cost. Indeed, Makogon has commented on the effect of electric fields on hydrate crystals 1, devising an expression for the shift in thermodynamic equilibrium conditions caused by a static electric or magnetic field, concluding that the electric-field intensity would have to be of the order of 107 V/m (i.e. 10-2 V.nm-1) to have any significant effect on equilibrium. Makogon also outlined a model for e/m field-induced melting of hydrate deposits from the porous rock of wells 1. Rojey has reserved a U.S. patent on a system to inhibit / prevent the formation of clathrate hydrates in pipelines by means of exposure to a field of e/m waves set up by a network of emitters spaced along the length of a pipeline

24

. Rojey stated briefly that the e/m fields, of appropriate frequency and

intensity, should prevent or minimise the organisation of water molecules to form a hydrogen-bonded crystalline lattice in a fluid containing an aqueous phase and hydrocarbons, disrupting the possibility of hydrate growth ACS Paragon Plus Environment

24

. External electric fields

The Journal of Physical Chemistry

Page 4 of 37 4

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

have the potential to act as an alternative for gas-pipeline flow-assurance to more established practices of using additives 1. Although electric-field penetration depths in liquid water are of the order of ~1 cm, in gas, this extends for up to meters

24

; this,

combined with repeated specular reflection along the pipe walls (in a ‘zig-zag’ fashion) serves to keep the amplitude relatively high over hundreds of metres (with Rojey’s patent specifying field emitters every few hundred metres of pipeline). Unfortunately, however, there has been essentially no developmental work of electricfield methods for gas-hydrate applications in the past 20 years reported in the open literature, with little such industrial R&D occurring. This relative neglect of this intriguing flow-assurance proposition serves as a motivation to us to study further electric-field effects on hydrates – especially for propane as a prototypical sII hydrateformer, given that its presence in natural-gas pipeline streams typically leads to the adoption of sII as being the primary polymorph type in any pipeline plugs (as opposed to sI for pure methane flows), with propane needing to occupy the larger 51264 cages in sII. Bearing in mind the scope and promise for molecular simulation in providing insights into hydrate microscopic processes, especially in electric fields 2, English and MacElroy have carried out non-equilibrium MD to assess how external e/m fields disrupt the formation of spherical methane-hydrate nano-crystallites

25

. They found

that the shifting dipolar alignment inside the crystallite weakens the structural hydrogen-bonding arrangements, thereby leading to break-up (or inhibition of formation in the first place). Naturally, this is dependent on field intensity and frequency, with intensities of the order of 0.1-0.5 V·nm-1 required to observe tangible results of inhibition within feasible MD simulation times. The ‘optimal’ frequencies for disruption were in the 20-200 GHz range (overlapping with the general hydrogen bond lifetimes). This confirms, to some extent, arguments advanced by ref. 25: in ACS Paragon Plus Environment

Page 5 of 37

The Journal of Physical Chemistry 5

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

particular, that e/m fields greater than 1.5 V·nm-1 in root-mean-square (r.m.s.) led to hydrate break-up under the effect of electric fields, i.e., “electro-dissociation” of hydrate nano-crystals, with the most rapid results in the 50-100 GHz range. However more recent non-equilibrium MD (NEMD) studies on bulk methane clathrates demonstrated that external e/m fields do not lead to “electro-dissociation” even with intensities approaching 2 V·nm-1 (r.m.s.), at least on timescales of the order of a nanosecond 26. Luis et al. applied subsequently static electric fields to bulk methane hydrate, as opposed to spherical nanocrystals, and found that fields greater than 1.5 V·nm-1 in intensity led to break up at a pressure of 2 MPa and temperature of 248 K 27. Luis et al. then performed simulations of further electro-dissociation of bulk methane hydrate at 260 K and 8 MPa, an also at 285 K and 40 MPa, in fields of 1 – 5 V·nm-1 28. In both cases, fields greater than 1.5 V·nm-1 led to dissociation, with formation of an (empty) ice-like and separate methane phase. Removal of applied fields led to ‘melting’ of the ice-like structure with separate liquid water and methane gas phase

28

. Waldron and

English also carried out NEMD studies on bulk methane hydrate, reporting an “electro-dissociation” threshold in the 1.2-1.6 V/nm range

29

, observing similar

gas/liquid nanoscale phase segregation upon melting. In addition to electric-field effects on quadrupole guests, the natural question arises as to potential differences in polymorph field response, especially in terms of sI versus sII dissociation. Bearing in mind the clear industrial/pipeline pertinence of sII hydrates in typical natural-gas environments mentioned previously, it is somewhat anomalous that all MD works on hydrate break-up in external electric fields have considered exclusively sI methane hydrates

25,27-29

. Indeed, one may conjecture that

water molecules in the larger 51264 cavity present in sII hydrates may be expected to respond more readily to external fields than those in the smaller 512and 51262 cages: ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 6 of 37 6

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 larger internal cage volume weakens somewhat the intrinsic electric field due to highly Coulombic water-water interactions. The intrinsic electric-field vector at any point in space is given by the Coulombic (electrostatic) force per unit charge thereat:26,29 in practice, one may evaluate this simply at charge sites by considering the electrostatic force acting thereon (conveniently available within MD at any instant, in any event) and dividing by its charge value. If one takes the magnitude of this vector and histograms for both sI and sII polymorphs, distributions in the ~10-30 V/nm emerge.26,29 In any case, the mean of the intrinsic field’s magnitude distribution in the centre of the larger 51264 cavities is some ~2 V/nm lower (in sII) vis-à-vis typical cage-centre values in sI polymorphs, owing to the larger cage sizes. This has the effect that the intrinsic field’s magnitude for sII polymorphs’ 51264 cages is lessened somewhat vis-à-vis those of the external fields in comparison to sI clathrates (although still being 20-2,000 times larger than external-field intensities used in the present study). As such, this has the effect that one would conjecture quite reasonably, ceteris paribus, that sII clathrates may be more prone to electric-field effects than sI for the same external-field intensity. In addition, this might be expected to lower the electro-dissociation threshold in sII clathrates with more ready molecular alignment, or else lead to more rapid decomposition for a given external-field intensity (exceeding the threshold). In addition to study of the dissociation threshold per se, in the present study, we consider the effects of static electric fields on the kinetics of hydrate dissociation. In any event, beyond hydrate electro-dissociation per se, there have been a number of previous molecular-dynamics (MD) studies of (thermal-driven) hydrate dissociation

30–36

. Báez and Clancy

30

pioneered the field by performing MD

simulations of dissolution of fully-occupied spherical hydrate clusters containing about 245 and 400 water molecules in surrounding liquid phases composed of melted ACS Paragon Plus Environment

Page 7 of 37

The Journal of Physical Chemistry 7

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

hydrate and pure water. They found that dissociation was essentially stochastic, taking place more in a stepwise fashion; the collapse of metastable distorted partial cavities at the hydrate interface, following diffusion of gas molecules therefrom, governed break-up and led to a series of rapid decays and plateaux in the nanocrystal size

30

.

English and co-workers studied the break-up of bulk methane-hydrate systems in static electric fields 29, decomposition of spherical methane-hydrate nanoclusters and planar systems by both perturbation by e/m fields

25

, and dissociation of both

nanoclusters and planar systems via conventional heating findings, grosso modo, to Báez and Clancy

31-33

, reaching similar

30

. However, they also found that the

diffusion of methane into the liquid phase (whether pure water, or locally supersaturated water-methane mixtures) is another important step which influences the overall break-up rate. Clarke and English

33

and Myshakin et al.

35

have carried

out MD simulations of the dissociation of 85, 95 and 100%-occupied planar gashydrate interfaces in contact with pure water at over-temperatures relative to the estimated melting point of up to around 30 K, and have concluded that the break-up rate displays an Arrhenius temperature dependence and is relatively sensitive to guest composition The goals of the current study are to study carefully the decomposition of “completed-cage” interface of propane hydrate (as a prototypical sII case study, owing its flow-assurance relevance) in contact with pure water at a variety of overtemperatures and in the presence of a static electric field, and to apply a previouslydeveloped model

2,33,37

to rationalise the observed hydrate decomposition; we will

identify the dissociation rates and identify the electric field- and temperaturedependence. Now, the usage of a completed-cage interface is also, in and of itself, somewhat novel: to date, the most popular version of a planar-hydrate/water interface has been associated with 001-direct surface cleavage, or “linearly-trimmed”. To have ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 8 of 37 8

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 more directly comparable result with experimental data, for which cages at the liquid interface would naturally be complete prior to any thermal- or field- driven dissociation, some extra water molecules were added to the model to close the open cavities from direct 001-direction cleavage. Although the previous work of English and MacElroy

25

and Bagherzadeh et al 38

have shown the shortcomings of treating the hydrate-liquid interface in an isothermal manner, it was desired to have a ‘proven basis’ of isothermal treatment in the simulations of the present study with which to compare to previous studies, to avoid the thorny complication of latent-heat-dissipation and these heat-transfer dynamics rendering phase-change kinetics more sluggish 2. It must be recognised, however, that care needs to be taken in the application and study of global-system thermostats, potentially resulting in artificially rapid system response; future work here needs to examine applications of boundary thermostats. However, dissipation of heat arising from mechanical work induced by the external electric fields, via torques acting on water molecules’ dipoles and friction from partial alignment, is sine qua non in 37

assessing athermal electric-field effects

: this affords another motivation in the

present work to employ a global-system thermostat, despite any concerns. SIMULATION METHODOLOGY Propane hydrate was described by standard classical force fields. The TIP4P/2005 model

39

was used for water. Four-site water models, of TIP4P-type geometry, have

been used in many simulation studies on hydrate structures.2 Huang et al. determined that TIP4P/2005 yields the correct energy sequence, along with reasonable structural parameters of various ice polymorphs, consistent with dispersion-corrected Density Functional Theory (DFT) calculations

40

. Conde et al. used different water models

(SPC, SPC/E, TIP4P, TIP4P/2005, and TIP4P/Ice) to predict the behaviour of

ACS Paragon Plus Environment

Page 9 of 37

The Journal of Physical Chemistry 9

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

methane-hydrate phase diagram, and found that the closest results to the experiment values can be obtained with the TIP4P/Ice model, followed by TIP4P/2005

41

.

Smirnov carried out melting and superheating of sI methane hydrate, determining that the TIP4P/2005 model gives superior agreement with the experimental coexistence line than the TIP4P/Ice model over a wider pressure range 42. In previous work,43 we determined that the melting point of propane hydrate was reproduced well by the TIP5P/2005 potential (vide infra). Thus, it is not at all unreasonable to construct a P-T phase diagram of water and gas hydrates based upon the TIP4P/2005 model. Propane molecules were parameterised by General Amber Force Field (GAFF) Atomic point-charges in propane were calculated by RESP procedure

45,46

44

.

based on

the HF/6-31G(d) electrostatic potential, consistent with Amber force fields. Interactions between molecules i and j, Vij, were calculated using pairwise potentials as the sum of Coulombic and Lennard-Jones terms:    =    +  ( ) =





×

  

+ 4 

 







!

−  " 

(1)

Here, qi and qj are the charges of atoms i and j, rij is the distance between these atoms, σ is the van der Waals diameter whilst ε is the depth of the VLJ potential well; σ and ε are given in Table 1. The parameters for pair interactions of atoms were determined using the Lorenz–Berthelot rule:47 # =

 $

(2)



 = %  

(3)

Table 1. Partial charges and van der Waals potential parameters used in calculations Atom

σ, Å

ε, kJ/mol

qatom, e

0.00E+00 7.75E-01 0.00E+00 Propane

0.5564 0.0000 -1.1128

Water H O M

0.000 3.159 0.000

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 10 of 37 10

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

C (CH3) H (CH3) C (CH2) H (CH2)

3.816 2.974 3.816 2.974

4.58E-01 6.57E-02 4.58E-01 6.57E-02

-0.24226 0.036563 0.50563 -0.12025

The bond-stretch and bond-angle intramolecular parameters for propane are listed in Table 2, with dihedral interactions set to nought, in a manner consistent with GAFF force-fields.

Table 2- Bond-stretch and bond-angle intramolecular interactions for propane

Bond Components C-C C-H Angles Components C-C-C C-C-H H-C-H

Bond-Length (Å) Force term (kJ/mol/Å) 1.5350 126847.3 1.0920 141160.1 Angles (Degree) 110.630 110.050 108.350

Force term (kJ/mol/deg.) 264.534 194.058 165.015

The cut-off radius for Lennard-Jones interaction parameters was 6.8 Å. The smooth particle mesh Ewald method 47-49 was used to handle long-range electrostatics. Although English and MacElroy 50,51, Jiang et al 52 and Burnham and English 53 have demonstrated previously that polarisable water models tend to be quantitatively superior, especially for the treatment of hydrates, this has been gauged, in the main, from agreement with experimental dynamical properties, such as vibrational and librational spectra (i.e., density of states) gleaned from neutron-scattering experiments. One important lacuna in the polarisable-potential literature as regards hydrates lies in critiquing the fidelity, or otherwise, of predicting thermodynamic coexistence accurately; in this respect, there is much to learn from other polarisable models such as rigorous work on, inter alia, Gibbs-Duhem analysis of the GPCM and BK3 of Nezbeda and co-workers in the context of pure-water phases.54-56 Indeed, although polarisable models have not been used in the current study, their future

ACS Paragon Plus Environment

Page 11 of 37

The Journal of Physical Chemistry 11

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

deployment towards propane hydrate would be welcome, bearing in mind gauging their thermodynamic as well as dynamical-property suitability. In any event, the reasonable prediction of the propane-hydrate melting point at the current pressure of interest,43 is reassuring for the present work (vide infra). Aside from the relative lack of use of polarisable potentials per se, there have been few attempts made to date in hydrate simulation to use specifically optimised guest-water interaction parameters, other than by Sun and Duan

57

and Anderson and co-workers

58,59

to use water and

guest potential models and interaction parameters which have been parameterised specifically for water-guest or hydrate systems. The velocity Verlet scheme was used for MD

47

with a time step of 2 fs. For

extended-system dynamics in the NVT ensemble, a Nosé-Hoover relaxation time of 0.02 ps was applied

47

. Although the use of system-wide thermostats in non-

equilibrium crystal-liquid systems (at temperatures above the crystallite’s melting point) would be expected perhaps to lead to artificially rapid dissociation

25,60–63

, it

was decided to employ this approach in the current study, as mentioned previously, as well as the essential matter of dissipation of external-field heating (via mechanical work acting on torques 37). In the hydrate phase, the starting coordinates of the oxygen atoms in the unit cells of sII hydrate were taken from x-ray diffraction data 64. The initial unit cell lengths were 17.32 Å and the initial orientations of the water molecules were selected in a random manner so as to conform to the Bernal-Fowler rules

65

and so that the total dipole

moment of the system would be vanishingly small by the Rahman-Stillinger procedure66. A 2×2×2 supercell was constructed by replication. The supercell was also surrounded in the + and – z-directions by relaxed liquid water ‘faces’ of equivalent cross-sectional area in the x-y plane, and of length of approximately 34 Å along the z-axis to construct a planar geometry crystallite. The constructed simulation ACS Paragon Plus Environment

The Journal of Physical Chemistry

Page 12 of 37 12

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

box is 34.02 × 34.04 × 107.44 Å, containing 15104 atoms. There were cavities at the interface which were initially open, as well as those which were complete. As mentioned previously, to avoid the crude and physically unrealistic practice of 001surface cleavage, we adopted a completed-cage configuration, in which specific extra water molecules added to the surface to complete and close all open cavities. Propane molecules were placed in all of the large cavities of each unit cell (100% large-cage occupancy). Prior to NVT production simulations under periodic boundary conditions (PBC), the hydrate framework was kept fixed during initial NVT MD runs for 200 ps, 280 K range to let the surrounding water molecules completely cover the hydrate structure. The simulations were run at different temperatures in order to determine (approximately) melting points (by the onset of dissociation), as well as carrying out outright non-equilibrium melting simulations at 250-300 K (over the respective melting temperatures), as a function of applied static electric field. NVT production simulations were carried out for circa 10 ns at the abovementioned range of temperatures. The geometric hydrate-ice-liquid distinction criteria of Báez and Clancy were employed to distinguish between the hydrate, ice lattices, and liquid-phase 30. Briefly, these involve the calculation of an angular order parameter, Fi, of each water molecule i to quantify the tetrahedral nature of bonding for nearest-neighbour water molecules, followed by the recognition of five-membered rings of water molecules present in hydrate structures but absent in liquid water and ice. The order parameter Fi is found by summation of bond angle cosines for the angles j-i-k:30 1 3 1 ∑-2$()*+,- ()*+,- + 0.11 & = ∑2



ACS Paragon Plus Environment

(4)

Page 13 of 37

The Journal of Physical Chemistry 13

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

where ni is the number of nearest neighbour for molecule i. For tetrahedral bonding in ice and hydrate crystals, there are six independent angles per molecule and the cosine of each angles ,- is close to -0.33. Therefore, the parameter Fi expresses the degree of deviation from perfect tetrahedral bonding: a value near zero indicates tetrahedral bonding, whilst larger values imply that the molecule is in a liquid-like region. This allows a preliminary classification of hydrate-, ice- and liquid-like molecules, which is refined further by grouping hydrate-like molecules into clusters and taking into account the identities of neighbouring water molecules.30 We ran non-equilibrium MD for various static-field intensities up to 1.4 V.nm-1 (given that the bulk-hydrate electro-dissociation threshold is ~1.5 V.nm-1)

27-29

, and it

would be expected to perhaps be lower for a planar-liquid interface. A particle of charge q in the presence of a static electric field (E) experiences a force

4 (5) = 67

(5)

Utilizing the velocity-Verlet algorithm, at each time step of particle position (r) and velocities (v), we have: (5 + ∆5) = (5 ) + ∆59 (5) +

9 (5 + ∆5) = 9 (5 ) +

∆: