Niall J.English*,Mohammad Reza Ghaani

School of Chemical and Bioprocess Engineering,University College Dublin,Belfield,Dublin 4,Ireland

Keywords:Molecular dynamics Clathrate hydrates Crystallisation Thermostatting Radial temperature profile Fluctuation-dissipation

ABSTRACT Molecular-dynamics(MD)simulations have been performed for the growth of a spherical methane-hydrate nano-crystallite,surrounded by a supersaturated water-methane liquid phase,using both a hybrid and globalsystem thermostatting approach.It was found that hybrid thermostatting led to more sluggish growth and the establishment of a radial temperature profile about the spherical hydrate crystallite,in which the growing crystal phase is at a higher temperature than the surrounding liquid phase in the interfacial region,owing to latent-heat dissipation.In addition,Onsager's-hypothesis fluctuation-dissipation analysis of fluctuations in the number of crystal-state water molecules at the interface shows slower growth.

1.Introduction

Clathrate hydrates are crystalline inclusion compounds wherein a host lattice composed of water molecules encages small guest atoms or molecules in cavities[1,2].Methane hydrates are the most abundant in-situ hydrate,occurring in Nature in marine-sediment milieux and in the permafrost,typically adopting the sI-polymorph structure[1,2].In recent years,global interest in methane clathrates has increased substantially due to their potential as a new source of natural gas[3-5].This worldwide interest has motivated a wide range of moleculardynamics simulations,often gauging methane-hydrate extraction methods,usually from marine-sediment environments[6].

More widely,molecular simulation has proven rather invaluable in elucidating a great deal about underlying equilibrium properties of gas hydrates,including structure[6-10],phonon‘resonant-scattering'by guest molecules in cages[11-13],hydrogen-bonding behaviour[14,15],and dynamical [7-23]and energetic properties [10,16].Such a greater understanding of equilibrium phonon scattering and timedependent properties offered by molecular dynamics(MD)has resulted in a great deal of progress,especially vis-à-vis elucidation of underlying thermal-conduction mechanisms in clathrates[24-28],as well as with“cage-hopping”-based guest diffusion — of relevance to gas-storage and energy-applications[29-33].

Despite this striking recent progress in simulation-facilitated understanding of clathrate-equilibrium properties,hydrate-formation kinetics are understood rather poorly[1,2].There have been a number of molecular-simulation investigations into hydrate nucleation and crystallisation[2].The interest in both the inhibition and formation of natural gas hydrates lies at the heart of a great many technical problems in the oil and gas processing industries[1].For instance,the injection of chemicals into gas streams,such as methanol and various glycols,is inhibitory towards hydrate nucleation,agglomeration and growth[1].The exact molecular-scale mechanisms of hydrate growth remain to be determined fully,although there is more clarity than for nucleation[2].Mass transfer of guest molecules is certainly one key issue,together with limited timescales amenable to MD[34],but the(mis-)use of artificial thermostats at the interface in MD is a further compounding problem in the case of many studies.The presence of the latent heat of freezing is an important factor at any hydrate interface,and the temperature distribution thereat is not necessarily expected to be uniform；application of(global),system-wide thermostats leads to an artificially isothermal interface with potentially too-rapid conduction of heat,and removal of physically important temperature gradients.

In MD simulation of hydrate growth,the previous works of,for instance,English and MacElroy[34]and Bagherzadeh et al.[35]have shown the shortcomings of treating the hydrate-liquid interface in an isothermal manner,and they have applied local thermostats.In further important work,Vatamanu and Kusalik have sought to address this with the use of ‘local' thermostats,further away from the interface,but leaving the subtle details of temperature gradients and microcanonical dynamics undisturbed near the interface itself[36-38].In a recent and insightful study,Liang,Kusalik and co-workers have carried out a comprehensive analysis on the nucleation of H2S hydrate employing both NVT/NPT and NVE approaches[39,40],and concluded that more crystalline nuclei are formed under NVE conditions,probably attributable to artificially fast nucleation rates under globallythermostatted(NPT or NVT)approaches.

In the present work,bearing in mind these cautious remarks on globally-applied thermostats,we extend the seminal,early work of English and MacElroy [34]to MD simulation of seeded methanehydrate growth,assessing in more detail the impact of global versus local thermostatting,as well as developing crystallite-size and fluctuation-dissipation analysis.It should be borne in mind that the evolution in local structure in the interfacial region from crystalline to liquid entails changes in local density corresponding to local minima in the free-energy surface；these structural,density and heat changes could also be studied on a local level using grand-canonical-ensemble(GC)analysis.This would highlight the local regions of heat generation and consumption.Although GC would be an interesting avenue to pursue,it was also desired to estimate dynamical properties of the system during crystallisation,so MD has been used here as the tool of choice.

In any event,the present study expands the work of Ref.[34],in that it seeks to characterise more quantitatively the local radial temperature profile at the crystallite/liquid interface,measure and express growth rates as a function of crystallite size,examine relative methane and water diffusivities and also assess crystallite-size fluctuation in the context of Onsager's hypothesis.Taken together,this expands the present study a great deal beyond Ref.[34]in terms of more quantitative rigour；indeed,Ref.[34]sought to garner a more qualitative understanding of hydrate-nanocrystallite growth.

2.Methodology

As in Ref.[34],various non-equilibrium two-phase systems were generated to initiate crystal-growth MD simulations,and more detail as to the preparation of the systems is available there.In brief,this was accomplished by seeding a methane-supersaturated aqueous liquid phase by approximately spherical methane-hydrate nano-crystallites with radii of around 0.9 and 1.25 nm.Most analysis in the present work will focus on the larger,～1.25 nm cases,owing to Ref.[34]concluding that these exceeded the critical radius by a larger margin and were more into the outright crystal-growth régime.The extent of the cubic simulation box was about 4.2 nm[34],which was found to be adequate in size(with respect to potential concerns of periodic replicas of crystallite images in neighbouring periodically-imaged boxes interfering with growth)[41]；at no point do the grown crystals in the present study,or indeed Ref.[34],envelop almost the entire simulation box,which would naturally be expected to add to artificially fast crystallisation kinetics[2,42].There were 1949 water and 339 methane molecules in the resulting 2288-molecule system containing the 1.25 nm-radius crystallite,and 1944 water and 331 methane molecules in the 2275-molecule system containing the core of radius 0.9 nm[34].The system-temperature/time protocol is described in Ref.[34]；in any event,three 2288-molecule systems were created,labelled as ‘1',‘2'and‘3',with increasing levels of initial NVT relaxation[34].Systems were relaxed briefly via NPT-based MD to 220 K and 60 MPa,with no significant changes in the size,structure or identity of the crystal cores noted [34].These nano-clusters were very similar for the three 1.25 nm crystallites[34].

This study utilised the rigid/polarisable TIP4P-FQ water potential[43]and a five-site rigid methane model,comprising a Lennard-Jones(LJ)12-6 interaction site on the carbon atom and partial charges on the carbon and hydrogen atoms[44].These models were chosen for the crystal-growth simulations because of their good performance visà-vis the systems using flexible and non-polarisable water models for the simulation of hydrates and water[23,34,45].The Lekner method was used to handle long-range electrostatic interactions[46-48]and was implemented using an interpolation scheme described previously[49].The cut-off radius for LJ interactions was set at 1.0 nm and the Lorentz-Berthelot mixing rules[50]were applied for interactions between different types of LJ sites.The Báez-Clancy hydrate-recognition criteria were applied,as outlined in Ref.[34].

As in Ref.[34],to probe further the suitability of global thermostatting in the NPT ensemble for modelling crystal growth,we adopted the hybrid microcanonical/canonical ensemble of Ref.[34],in which the molecules whose centres of mass lie in a central spherical region of the simulation box,of radius 1.95 nm,are subjected to NVE dynamics,whilst the kinetic energy of the molecules in the surrounding region is regulated by means of a thermostat (i.e.NVT dynamics is applied to the molecules therein,with a set temperature of 220 K).The radius of the NVE region was selected as 1.95 nm so that approximately half of the molecules would be in each area of the simulation box.During crystal growth at 220 K and 60 MPa,there is a release of the latent heat of crystallisation at the solid/liquid interface and the consequent conduction of thermal energy away from the phase boundary.For the“NVE/T”case [34],the box volumes were chosen so that the total system pressure would not deviate greatly from 60 MPa,to afford a comparison with NNPT results at 220 K and 60 MPa.

A 1.5,15 fs timestep structure was used[49],in conjunction with a TIP4P-FQ charge mass of 15 × 10-5(ps · eV-1)2·kcal · mol-1.For NPT dynamics,τQand τBwere set at 100 and 500 fs,respectively.For the hybrid NVE/T ensemble,τQwas set to 50 fs.In all cases,τQwas set to 50 fs for the control of the TIP4P-FQ electronic temperature,with the set value at 1 K.

3.Results and Discussion

Because crystal formation was observed in ref.[34]at a temperature and pressure of 220 K and 60 MPa,these conditions were selected for the simulations of this study.Each simulation was of 0.5 ns duration,during which appreciable growth on the nano-cluster“scaffolds”was seen.The number of hydrate-like methane molecules as averaged over the three initial-1.25 nm crystallite systems(1,2,3)is shown in Fig.1 for both the NPT and NVE/T cases.For both cases,a straight-line regression fit has been included on the size-time data after an initial phase of enhanced crystal formation lasting about 30 ps.As found in Ref.[34],the crystal size (in terms of number of hydrate-like molecules),both in terms of the number of enclathrated water and methane molecules,exhibits rapid fluctuations about an average value for small time intervals.The initial transient régime of rapid growth,over the first 30-40 ps arises primarily as a result of rotational relaxation of water molecules at the hydrate-liquid interface.Thereafter,crystal formation continued more slowly,with significant fluctuations in crystal size,(system 1)and 0.9 nm crystallites,respectively.The substantial variations in the cluster size and structures illustrated in snapshots(cf.Fig.4 of Ref.[34])result from instantaneous fluctuations in hydrate-phase classification.The newly attached hydrate structure was deformed to some extent vis-à-vis the original nanocrystal template(cf.Fig.4 of Ref.[34]).The water molecules in the newly attached part formed mostly pentagons and some tetramers.The interface formed between the hydrate cluster and the surrounding liquid phase consisted mostly of partial hydrate cavities holding one or no methane molecules.The plateaux and occasional substantial rises in the hydrate size-time results suggest that growth proceeds by consolidation and completion of individual cages.In order to stabilise the cavities to the greatest extent,methane molecules must diffuse from the surrounding liquid phase to occupy and vibrate about the lattice position,even in the case of partial cages,and this is evident in“staged”rises in Fig.1.Naturally,this methane-diffusion and water re-ordering process is a highly statistical,entropy-driven process.

Fig.1.(a)Number of enclathrated methane molecules as a function of time for NNPT case.(b)Number of enclathrated methane molecules as a function of time for the NVE/T case.

Fig.2.Time-averaged radial temperature profile in the NVE/T simulation.

Study of the NVE/T methane-size/-time results for the three 1.25 nm crystallites in Fig.1b shows that crystal growth after the initial 30-40 ps reorientational-adjustment phase is somewhat more sluggish than for the NPT case (in Fig.1a).For NVE/T,the overall system temperature and pressure were found to remain close to 220 K and 60 MPa throughout the simulations.Rotational relaxation is facilitated in the initial phase of the run in the NVE/T ensemble vis-à-vis the NPT case,as demonstrated by the slightly more rapid initial increase in crystal size(cf.Fig.1b vs.a,and Fig.2 in Ref.[34]).However,the slower growth rate after this period may be rationalised in terms of the more realistic treatment of latent-heat release at the crystal surface followed by heat conduction in the NVE/T ensemble:crystallisation proceeds by(i)transport of molecules to the surface(diffusion),(ii)incorporation at the interface(surface kinetics),and(iii)the release of latent heat to be transported away from the surface(heat conduction).If all of these processes are sufficiently rapid,ideal growth laws may be realised.In reality,however,the slowest process governs the growth rate as a whole.In the NPT ensemble,the temperature in the whole system is regulated and the issue of heat conduction from the interface is obviated.In the NVE/T case,one would expect the crystallite and the crystal-liquid interface to be at a higher temperature than the surrounding liquid phase,as the latent heat is conducted away from the crystal surface.The radial temperature profile for the NVE/T runs is shown in Fig.2,time-averaged over the duration of the simulation.The non-uniform profile was reasonably steady over the run.As found in Ref.[34],the temperature of the crystallite itself was higher than the liquid,being about 223 K in the central area (up to a radius of 0.8 nm),declining at the interface(at a radius of approximately 1.1 to 1.4 nm)to a liquid phase value of 217 K.Beyond the radial distance of 1.95 nm,the temperature was uniform at the set value of 220 K.The fact that the interface was higher in temperature than the surrounding liquid interfacial layer serves to decrease the rate of molecular attachment to the surface,explaining the more sluggish crystal growth relative to the NPT ensemble.It was found that the deposited cages in the NVE/T case,although fewer than for the NPT ensemble,were more crystal-like in structure；this is in agreement with the greater level of crystallinity in(similarly more sluggish)hydrate formation under NVE conditions,as compared to NVT,encountered in Refs.[39 and 40].Simply put,there is more ample time with slower conductive latent-heat dissipation for locally-microcanonical dynamics to allow for more water re-ordering around guest molecules to result in greater local,cage-level crystallinity at the growing crystallites' interface with the liquid.

An approximate heat balance for the crystallite in NVE/T system is as follows:

where Ncrysis the number of molecules in the crystallite,Cv,crysis the heat capacity of the crystal,λ is the latent heat,˙Ncrysis the crystal formation rate and qcondis a term accounting for the rate of heat loss through conduction from the surface during growth.Integrating Eq.(1)gives

The growth of the smaller crystallite cores(with an initial approximate radius of 0.9 nm)showed similar qualitative trends to those in ref.[34]；the completion of the initial partial cavities results in a hydrate cluster after 0.5 ns of similar size to the initial 1.25 nm radius nanocrystals,albeit with less order in the newly attached part.For the given size of the smaller crystal,the growth rate appears to be consistent with that of the larger crystal(vide infra).This implies that the mirror images of the crystal in the surrounding infinite periodic array of replica simulation boxes [41,42]do not influence the growth of the larger crystallite greatly.

Fig.3.a.Oxygen-oxygen radial distribution function at various stages of the NPT crystallisation.b.Oxygen-carbon RDF at various stages of the NPT crystallisation.

Site-site radial distribution functions(RDFs)are depicted in Fig.3 at various stages throughout crystallisation,averaged for the systems with the 1.25 nm-radius crystallites.The RDFs were sampled over 4 ps windows,and the reported times of 2,6,58 and 500 ps correspond to the centre of these sampling periods.It may be seen that the peak structure in the water-water RDFs(cf.Fig.3a)becomes more pronounced as time progresses,i.e.,an increase occurs in the heights of the principal and secondary peaks along with a concomitant increase in the depths of the troughs.The increase in the number of hydrogen bonds present in the system,associated with the incorporation of water molecules onto the original nanocrystal template is manifested in the O-O RDF.Similar trends are observed in the O-C RDF(cf.Fig.3b).The initial rapid crystal‘growth',arising from reorientation of water molecules in the liquid phase at the interface,is readily apparent in the RDFs,in terms of the greater disparities of RDF results at earlier times,e.g.2 ps versus 6 ps.Reference to the RDFs in Refs.[49 and 23]for equilibrium liquid water and methane hydrate structures with TIP4P-FQ,respectively,shows that the RDFs in Fig.3 resemble those of the liquid state more closely.This is as expected,since 80%-90%of the molecules in the system are in the liquid state；as such,the RDFs serve only as a qualitative indication of the changes in system structure occurring as a result of hydrate formation.

The mean square displacement(MSD)curves,originated at the true zero simulation time,are presented in Fig.4 for methane molecules(with those of water shown in Fig.5 of Ref.[34]).The MSD indicates the extent of molecular mobility and diffusion,especially in the liquid phase and close to the crystal interface.Results are shown for all three 1.25 nm-radius-crystallite systems in the NPT ensemble,and for system 1 in the NVE/T ensemble,in addition to the curves for 0.9 nm-radius-crystallite system(NPT).The MSD values are a factor of 2 to 2.5 times larger for the methane molecules relative to the water molecules [34].This is attributable to the possibility of these molecules forming local regions in the liquid phase richer in methane,in addition to methane molecules diffusing to various partial cavities,whereas the water molecules must adopt specific orientational alignments to be incorporated into a(partial)cage structure,which requires a greater amount of time.Another reason that methane diffuses more quickly and independently is because of its weaker interactions with water relative to water-water interactions.The water molecules are hydrogen bonded,and these bonds would need to be adjusted in order for water molecules to diffuse independently,i.e.the activation energy EAfor diffusion of water molecules,characterised qualitatively by Eq.(4),is higher:

Fig.4.MSD of the methane molecules.

The slower diffusion of water molecules may be rationalised further in a qualitative manner,using the Wilke-Chang correlation for diffusion coefficients of solutes in dilute solutions:[51].

where the solute diffusivity is in m2·s-1,φ is an association factor for the solvent(2.6 for water),M and μ are the molecular masses and the viscosity of the solvent,respectively,T is the temperature and Vmis the molar volume of the solute at its boiling point in m3·kmol-1.Since the activation energy for the complete dissociation of a water molecule from the hydrogen bonds to its nearest neighbours is likely to be quite high,it is unrealistic to consider water molecules as diffusing independently.Using the molar-volume values of 0.0189 and 0.0296 m3·kmol-1reported by Gambill[52]for free water and methane molecules,respectively,one may tabulate theparameters for methane molecules and various water clusters,which gives an approximate indication of the extent of molecular mobility — cf.Table 1.Theratio is 2.24,which agrees well with the observed ratio of methane to water MSD slopes of 2 to 2.5(cf.Fig.4 vs.Fig.5 of Ref.[34]),corroborating the idea that water molecules tend to diffuse in hydrogen-bonded clusters,e.g.tetramers,pentamers,hexamers or larger groups.This is an important point,in that the prevailing‘picture'of liquid-phase diffusion is often predicated on the notion of isolated,individual-molecule events；however,given the importance of reorganisation of solvent coordination shells around solvated methane molecules (with coordination number of 20)to 24 for large 51262cages in methane (sI)hydrates [1],the diffusion of both methane and water requires local,collective processes of local structural rearrangement on a continual basis.Therefore,it should be less surprising that diffusion takes place featuring these more clustered formations of molecules,showcasing the importance of hydrogen bonds in collective rearrangement and diffusion.

Table 1 Wilke-Chang diffusivity parameters for methane and water clusters

After the initial period of reorientational relaxation and the greater molecular mobility associated with this realignment,the slopes of the MSD curves level out.The MSD curves for the three 1.25-nm-radius-crystallite systems in the NPT ensemble are very similar for water molecules,whilst the results for methane molecules are in somewhat less absolute agreement but exhibit similar long-time slopes.Clearly,the role of statistics,as demonstrated by the use of several simulations starting from similar conditions,is important,especially for the methane molecules.The lower rate of crystal growth exhibited by system 1 in the NVE/T ensemble manifests itself in smaller and slightly less steep MSD curves,as less molecules from the surrounding liquid are transported to the hydrate surface,due to the interface being at a higher temperature than the surrounding interfacial liquid layer(cf.Fig.2).The MSD curves for the 0.9-nm-radius-crystallite system are similar to those systems with the larger hydrate clusters,implying that the molecular mobility in the surrounding liquid phase makes the substantial contribution to the MSD.

The crystal-formation rates,beyond the initial transient phase of interfacial reorientational relaxation,as characterised by the regression lines in Fig.1 and Figs.1-3 of Ref.[34](vide ante),are specified in Table 1 of Ref.[34].The formation rate is consistent for the three 1.25nm radius nanocrystal systems in the NPT ensemble,being about 0.32 water molecules per picosecond and 0.046 methane molecules per picosecond.The water growth rate is about seven times greater than the methane formation rate:less methane is being incorporated into the crystal than may be inferred from the equilibrium water to methane ratio of 5.75:1.This leads to non-equilibrium defect structures in the newlyadded parts of the crystallites,with a greater occurrence of empty complete and partial cavities.As mentioned earlier,the growth rate for the NVE/T ensemble is lower than the corresponding rate in the NPT case,being about two thirds of the NPT value,due to the more realistic modelling of heat conduction from the crystal interface.The formation rate of 0.182 water mol·ps-1and 0.014 methane mol·ps-1for system featuring the nano-crystallite with an initial radius of 0.9nm is lower than that of the systems with the larger initial hydrate cluster.It shall be shown,however,that for the given size of the smaller crystal,the growth rate appears to be consistent with that of the larger crystal.

It is useful and instructive to translate the crystal-formation rates defined above into radial growth rates for the crystallite.Making the approximation that the growing crystal is spherical (as confirmed by visualisation)[34],then the growth rate R may be estimated also fromR ＝where the radius of the crystal,r,is known as a function of time.Using the linear regression fits shown in Fig.1,and Figs.1-3 of Ref.[34],produces t hevalues listed in Table 1 of Ref.[34](where this is～0.32 and 0.045 water and methane molecules per picosecond,respectively,for 1.25nm crystallites under NPT conditions,and～0.22 and 0.029 mol·ps-1respectively for NVE/T)[34].The crystal radius may be found from the volume of the hydrate crystal,Vcrys,by r ＝The crystal volume may be estimated by

where NS(W)and NS(M)refer to the number of water and methane molecules in the hydrate cluster,andare the volumes per water and methane molecule,respectively,of the hydrate phase.The smoothed results for NS(W)and NS(M),depicted by the regression curves in Fig.1,and also in Figs.1-3 of Ref.[34],may be used to define the time evolution of Vcrysusing Eq.(6),once reasonable estimates forandare available.This allows then the calculation of the crystal radius as a function of time,from which the growth rate R may be found by numerical differentiation of r.Alternatively,one may differentiate r analytically:

Eq.(7)shows a non-linear relationship between R and NS(W)and NS(M),i.e.a declining radial growth rate as the crystal increases in size,as one would expect for the growth of a spherical crystal[53].

Prior to the use of the procedure outlined above,it is necessary to obtain reliable estimates ofAn NPT equilibrium simulation of eight hydrate unit cells(consisting of 368 water and 64 methane molecules),analogous to those of Ref.[23],was carried out for 0.1 ns at identical temperature and pressure conditions to the present crystalformation studies(220 K and 60 MPa).It was found that the isotropic box length was 2.324 nm.The methane molecular radius was found to be 0.218 nm in the clathrate state by von Stackelberg and Müller from X-ray diffraction experiments [54].Usingwhere rM=0.218nm,gives=0.0434 nm3.Using the box volume from the clathrate simulation at 220 K and 60 MPa and subtracting the total volume occupied by the methane molecules(i.e.64×0.0434 nm3)results in a specific volume per water molecule of 0.02656 nm3.Although it may appear,at first glance,somewhat unorthodox to use an experimental estimate forto then infer a(primarily)simulation-based one forit is permissible,and even necessary,to adopt this approach,because it is difficult to establish exactly the de-facto methane,nearspherical volume in clathrate cages directly from molecular simulation.The Lennard-Jones parameter for the carbon atom in methane may serve as a rough proxy for this,but it is much more rigorous to use a“tried-and-tested”experimental estimate forin the sI hydrate cages,allowing a reasonable volume partitioning then forin a manner that is consistent,and pragmatic,with the molecularsimulation approach.

The approximate smoothed crystal radii versus time are presented in Fig.5,determined from Eq.(6)in conjunction with the smoothed number of water and methane molecules represented by the regression lines in Fig.1 and Figs.1-3 of Ref.[34].The additional period of NVT melting of the surrounding liquid layer for systems 2 and 3 is apparent in the larger initial crystal sizes vis-à-vis system 1 for the NPT crystal formation simulations with the 1.25-nm radius crystals,as the more relaxed liquid layer surrounding the crystal surface is more ordered and hence more liable to an increased rate of reorientational alignment.Once this initial,～30 ps period is complete,the formation rates for the three systems become consistent,as noted previously in Ref.[34].The size of the clathrate cluster for system 1 in the NVE/T ensemble is initially greater than its counterpart in the NPT ensemble,as the more reasonable treatment of heat conduction seems to allow for an enhanced rate of interfacial rotational realignment,but the size increases less quickly thereafter,as commented on earlier in Ref.[34].Similar previous trends are noted in the size progression of the smaller crystallite.

Fig.5.Approximate smoothed crystal radius.

Fig.6.Approximate smoothed crystal growth rate(nm·ps-1)for crystallisation(notation same as Fig.5).

The radial formation rates R(in nm·ps-1),beyond the initial transient period are shown in Fig.6,calculated by differentiation of the crystal radii[cf.Eq.(7)].The growth rates decrease from 0.0005 nm·ps-1to less than 0.0004 nm·ps-1for the three systems containing the 1.25nm radius crystal in the NPT ensemble.The radial rate for the system containing the smaller clathrate crystallite agrees well with that of the value for the larger cluster.The NVE/T rate is about two thirds the NPT rate for system 1,as observed previously.The radial growth rates decline with time.This is hardly surprising,in that the growth rate of the crystal,as characterised by Eq.(7),declines as the size of the crystal increases with time.The formation rate is depicted in Fig.7 as a function of the smoothed crystal radius.It can be seen that the three systems in the NPT ensemble are reasonably consistent,and that the NVE/T rate is lower for a similar size range.However,the rate for the smaller crystallite appears to be consistent with the rate for the larger crystallite,when one considers that the initial size of the 0.9nm cluster is probably not much greater than the critical size at the given temperature and pressure.This suggests that the mirror images of the nano-crystallite in the surrounding infinite periodic array of replica simulation boxes do not influence the growth of the larger crystal significantly[41,42].

Fig.7.Smoothed crystal growth rate (nm·ps-1)versus approximate crystal radius(notation same as Fig.5).

It is apparent from Fig.1 and Figs.1-3 in Ref.[34]that substantial instantaneous fluctuations in the number of water and methane molecules in the clathrate phase take place about the underlying rates of formation,as characterised by the regression lines.It is interesting to investigate the dynamics of these fluctuations.The fluctuation dynamics may be related to Onsager's hypothesis,which states that fluctuations about the equilibrium state decay on average according to macroscopic laws.Since the crystal-formation rates are relatively slow,one may attempt to model the fluctuations about the underlying slow nonequilibrium phenomenon within the same framework.From Fig.1 and Figs.1-3 of Ref.[34],the instantaneous number of solid water or methane molecules shall be denoted as NS(i)(t),where i=W or M.The corresponding number on the underlying regression line shall be designated asmeaning that the instantaneous fluctuation,ΔNS(i)(t),would be ΔNS(i)(t)=NS(i)(t)-.The normalised autocorrelation function of ΔNS(i)(t)is

The decay of the ACF,C(t),may be treated as a negative exponential,i.e.C（t）＝A exp（）beyond a certain short-time transient,according to the Onsager hypothesis.

The autocorrelation functions (ACFs)of ΔNS(i)(t)are shown in Figs.8a and 9a for water and methane,respectively,up to 6 ps on a log-linear plot.The initial period of interfacial relaxation was excluded in their compilation,but the calculation was based on crystal number fluctuations up to 0.5 ns.In the case of Fig.8a,for water,the ACF decays more quickly,so the ACF axis is plotted up to 0.1,rather than 1.Linearregression fits of the form A exp（）have been included with the original curves beyond 2 ps.In the case of the three systems containing the 1.25nm-radius crystal in the NPT ensemble,the effects of the additional period of NVT melting in the system preparation [34]appearevident on the ACF.Although the background growth rates are identical,the actual ACF values are larger for system 3,followed by system 2,and then by system 1,indicating that the ΔNS(i)(t)fluctuations are most correlated for the most relaxed system.Also,the actual noise,or extent of fluctuation,in the C(t)curves decreases with a greater degree of overall system relaxation,e.g.there is less noise for the system 3 curve than for system 1.The ACF curves occur at lower values still,and have a more substantial degree of noise,for system 1 in the NVE/T ensemble and for the system containing the smaller crystallite,indicating that the interfacial structure in these systems appear less relaxed.

Fig.8.a:Normalised ACF of crystal molecule number fluctuation for water(based on crystal number —time data for 30 to 500 ps).b:Normalised ACF of crystal molecule number fluctuation for water(based on crystal number—time data for 350 to 500 ps).

Fig.9.a:Normalised ACF of crystal molecule number fluctuation for methane(based on crystal number —time data for 30 to 500 ps).b:Normalised ACF of crystal molecule number fluctuation for methane(based on crystal number—time data for 350 to 500 ps).

The time constants τ for the exponential decay are specified in Table 2.The results for τ suggest that as the system becomes morerelaxed,the decay rate of C(t)declines,i.e.τ becomes larger.The τ values are larger for methane than water.There does not appear to be any consistent connection between the actual magnitude of C(t),noise or decay dynamics and the underlying,non-equilibrium clathrate formation rate,at least when one considers the whole simulation period(apart from the initial period of interfacial relaxation)for the definition of the ACF.

Table 2 Decay-time constants for crystal-number fluctuation ACFs

The finding that the ACF curves of Figs.8a and 9a seem to be dependent on the state of initial interfacial relaxation suggests that the systems have a long memory of the initial non-equilibrium structure,which persists in the fluctuation dynamics for perhaps several hundred picoseconds.To investigate this,the time range of the number-time data of Fig.1 and Figs.1-3 of Ref.[34],from which the ACFs were extracted,was subdivided into three separate segments with a length of 0.15-0.16 ns,and C(t)defined for each segment up to 6 ps,with A exp（-t/τ）fits applied beyond 2 ps.One might expect the resultant τ values,as specified below,to be more consistent in the later segments,and this is found to be the case in Table 3 where these values are tabulated.The results for τ in the last segment(i.e.,0.35-0.5 ns)are in better agreement for the three different initial geometries and for the NVE/T case for both water and methane.The ACFs for the last segment are shown in Figs.8b and 9b for water and methane,respectively.They are much closer in actual magnitude and in decay dynamics(as characterised by τ)than their analogues defined for the whole of the number-time data.This suggests that the systems have largely‘lost memory'of the initial non-equilibrium state by this stage,and that the time required for this effect on the fluctuation dynamics is in excess of 0.3 ns.

Table 3 Decay time constants for ACFs from segmented number-time data

4.Conclusions

For the hydrate nano-crystallites,the radial growth rates declined with time and crystal size,as expected.In the case of the system containing the smaller crystallite,the radial hydrate formation rate was consistent with the value for the larger cluster,especially when one considers that the initial size of the 0.9nm cluster may not be much greater than the critical size at the given temperature and pressure.

The formation rate was more sluggish in the hybrid NVE/T ensemble,being about two thirds of the rate of the corresponding system in the NPT ensemble.This may be rationalised in terms of heating at the interface as crystallisation occurs；indeed,we are not aware of a radial-temperature profile(Fig.2)having been presented in the literature before for molecular simulation of crystallisation.In addition,the finding of more crystal-like cages is significant and important for the more sluggish (and realistic)crystal growth experienced under locally-microcanonical dynamics at the crystallisation front,which agrees with systems under full (global)NVE dynamics vis-à-vis globally-thermostatted ones[39,40].Although often neglected by the vast majority of molecular-simulation practitioners attempting to study crystallisation (or melting)processes,latent-heat dissipation must ideally be handled with care in both time and space,rather than being“eliminated”by artificial and unrealistic global thermostatting to be an instantaneous process[2].In that respect,the early description of a radial-temperature profile in crystallisation by English and MacElroy was important [34],together with important advances in local themostatting of crystallisation professes in refs.[35-38].It is to be hoped that more sophisticated“hybrid”thermostatting schemes shall be developed to allow more physically reasonable latent-heat dissipation in the locale of (diffuse)crystal-liquid interfaces.Indeed,more rigorous thermostatting for more general systems in which there is a transfer of latent heat is also worth tackling,such as mineral-/fluid-adsorption problems and accurate estimation of thermodynamic and transport properties.

Acknowledgments

We thank J.M.D.MacElroy and Peter Kusalik for interesting conversations,as well as the anonymous referees whom provided insightful and helpful comments.MRG thanks the Irish Research Council for his Government-of-Ireland postdoctoral fellowship,under grant no.GOIPD/2016/365.