Homogeneous SPC/E water nucleation in large molecular dynamics simulations
HHomogeneous SPC/E water nucleation in large molecular dynamics simulations
Raymond Ang´elil, J¨urg Diemand, Kyoko K. Tanaka, and Hidekazu Tanaka Institute for Computational Science, University of Zurich, 8057 Zurich, Switzerland Institute of Low Temperature Science, Hokkaido University, Sapporo 060-0819, Japan (Dated: August 28, 2018)We perform direct large molecular dynamics simulations of homogeneous SPC/E water nucleation,using up to ∼ · molecules. Our large system sizes allow us to measure extremely low andaccurate nucleation rates, down to ∼ cm − s − , helping close the gap between experimentallymeasured rates ∼ cm − s − . We are also able to precisely measure size distributions, stickingefficiencies, cluster temperatures, and cluster internal densities. We introduce a new functional formto implement the Yasuoka-Matsumoto nucleation rate measurement technique (threshold method).Comparison to nucleation models shows that classical nucleation theory over-estimates nucleationrates by a few orders of magnitude. The semi-phenomenological nucleation model does better,under-predicting rates by at worst, a factor of 24. Unlike what has been observed in Lennard-Jonessimulations, post-critical clusters have temperatures consistent with the run average temperature.Also, we observe that post-critical clusters have densities very slightly higher, ∼ J vs. S scaling relation using both experimental and simulationdata, finding remarkable consistency in over 30 orders of magnitude in the nucleation rate range,and 180 K in the temperature range. PACS numbers: 05.10.-a, 05.70.Fh, 05.70.Ln, 05.70.Np, 36.40.Ei, 36.40.Qv, 64.60.qe, 64.70.F, 64.60.Kw,64.10.+h, 68.35.Md, 83.10.Mj, 83.10.Rs, 83.10.Tv
I. INTRODUCTION
The vapor-to-liquid transition of water is a common phe-nomenon in nature, relevant to many areas of technologyand science. Attempts to predict the rate of homogeneouswater nucleation often fail because of the lack of under-standing of the properties of the tiny seeds of the interme-diate phase, which are not necessarily large enough to havereached the bulk liquid properties. The relevant proper-ties of the tiny clusters which affect predicted nucleationrates include surface tension, temperature, and density.Molecular dynamics simulation has proven to be a pow-erful test of thermodynamic analytical nucleation models,now that codes are efficient enough, and computers fastenough. Realistic, atmospheric nucleation rates are toolow to be possible in direct computer simulations, due tothe large number of molecules required. The lowest waternucleation rates performed in simulations and reported inthe literature are ∼ − cm − s − [1, 2], usually beyondthe spinodal limit. Laboratory water nucleation rates onthe other hand are far lower - usually < cm − s − , al-though a few experiments have managed to measure farhigher rates ∼ cm − s − [3–5]. Our simulations of ho-mogeneous SPC/E water nucleation, which we report on inthis paper, manage to close the gap considerably, resolvingnucleation rates down to ∼ cm − s − .Nucleation models, which seek to provide explanationsand predictions for nucleation rates, have a long history offalling short when compared to experimental results[3, 4, 6–14]. For the case of water, rate predictions from the classi-cal nucleation theory disagree with experimental measure-ments by factors of 10 − [5–7, 15–23]. These models alsohave difficulty when predicting rates measured in numeri-cal molecular dynamics nucleation simulation experiments.However, with molecular simulation, one can make mea- FIG. 1. A slice through the simulation T325f after 581 ns. Thecolor-map indicates the density, meaning that the white spotsrepresent large clusters. By the end of the simulation, the largestcluster in this run has 527 members. This simulation box is ∼ µ m × µ m × µ m, although only a thin slice into the z − direction is visible. surements more detailed and accurate than what’s possiblein laboratory experiments. Size distributions, nucleationrates, cluster densities, temperatures, and even cluster pres-sures, shapes, angular momenta, and surface tension mea-surements are possible. Understanding the properties of thetiny yet complex, many-body clusters which form is vitalfor the development of a complete and successful thermo-dynamic description of the phase transformation[24]. Sim-ulations allow us to identify the shortcomings in the as- a r X i v : . [ phy s i c s . a t m - c l u s ] J u l umptions made by existing nucleation models, and sug-gest ways they may be improved. Cluster properties arenoisy, necessitating large systems with many millions ofmolecules. This demands costly compute power, and onlyrecently have some of these direct measurement techniquesbecome possible[25–27].Direct vapor-to-liquid molecular dynamics simulation fora Lennard-Jones fluid has become a popular exercise duethe computational accessibility of the short-range, single-site potential[26, 28–32]. Water is significantly more de-manding. For the same system size, more complicatedmolecular interaction potentials like SPC/E[33–35] andTIP4P[2, 36, 37] necessitate a few orders of magnitude morecomputational power than a pure Lennard-Jones simula-tion. An exception is mW water, a comparatively simplemonoatomic single-site water model[38–44]. MD nucleationsimulations of mW water have been carried out, yet only onsmall systems with relatively high nucleation rates[44, 45].The monoatomic water model proposed by Zipoli et al.(2013)[46] offers similar advantages. However, we foundthat short-range potentials require extremely long equili-bration times to form the correct equilibrium abundanceof small clusters (dimers, trimers, etc.) in a supersatu-rated vapor, because interactions are rare, especially thethree body encounters required for dimer formation. Thisdrawback makes it computationally expensive to simulaterealistic, steady state vapor-to-liquid with such short-rangepotentials - despite their low cost per time-step - and wedo not use them in this work. Matsubara et al. (2007) [34]simulate homogeneous vapor-to-liquid nucleation using theSPC/E water model and include a Lennard-Jones carriergas, measuring rates down to 2 . · cm − s − . SPC/Esimulations by Tanaka et al. (2014)[1] manage to reachnucleation rates 3 · cm − s − . Both efforts addition-ally measure critical cluster sizes, formation energies, sizedistributions and sticking probabilities for systems in the T = 300 −
390 K, providing ample opportunity for modelcomparison and development.In this study, we continue in similar spirit, yet simulat-ing the SPC/E water vapor-to-liquid phase change in evenlarger computational volumes using longer time integra-tions. This allows for the measurement of lower nucleationrates than previously possible by a few orders of magnitude,and for the first time, measurements of naturally-formedSPC/E cluster density and temperature profiles. Our re-sults provide opportunities for the verification and calibra-tion of the standard assumptions which go into nucleationmodels, in a previously unexplored temperature and satu-ration regime.
II. SIMULATIONSA. Simulation code, setup and parameters
We use the molecular dynamics SPC/E[33] water model.SPC/E is a rigid 3-site model, which registers Coulombicinteractions, as well as polarization corrections to each site, and further adds a Lennard-Jones component to the oxygenatom potential.The Large-scale Atomic/Molecular Massively ParallelSimulator (or LAMMPS) computer program[47], developedat the Sandia National Laboratories and distributed underthe GPL license, was used to perform the SPC/E simula-tions. We have verified that our runs produce the sameresults as found in similar, yet smaller SPC/E numeri-cal nucleation experiments[1], which used an independentmolecular dynamics code. We cut the short range Lennard-Jones component to the force field off at 9 . t = 2 fs, commonfor SPC/E water simulations[1]. A typical simulation runsfor 72 hours on 1024 cores. Our largest simulation ran for1000 hours on 8192 cores on the Piz Daint supercomputerat Centro Svizzero di Calculo Scientifico (CSCS), perform-ing 3 · integration time-steps.The simulation box has periodic boundary conditions.Initially the molecules are given random non-overlappingpositions and random velocities. This is done at 1000 K,after which the ensemble is cooled and the box size ex-panded until the simulation reaches the target temperatureand pressure. The run continues in this state under NVTconditions, regulated by a Nose-Hoover thermostat[51–53]with temperature damping timescales of 1000 fs.At this stage the gas is allowed to equilibrate for a fixedamount of time - dependent on the run temperature (Re-fer to table I for the chosen equilibration timescales t e ateach temperature). During this phase the subcritical clus-ter equilibrium distribution forms. Around this stage webegin to make nucleation rate, size distribution, and clustergrowth rate measurements. For most runs, the nucleationrate is low enough that unnatural effects from the interven-tions due to the thermostat are minimal. Our nucleationrates are low enough that the latent heat of transformationin the simulations is extremely small, resulting in only afaint influence from the thermostat. Our largest run seesa total energy increase of ∼ .
1% over the steady-statephase, i.e. our simulations are very close to NVE (micro-canonical) ensembles.The first few columns of table II lists the runs which werecarried out, their target temperatures, box sizes, number ofmolecules, and their run times.2
ABLE I. Thermophysical quantities and parameters at eachtemperature. The vapor equilibrium pressure P v , the planarsurface tension γ , and the bulk liquid density ρ b for SPC/Ewater are determined from the fitting functions in Matsubara etal. (2007)[34]. η and ξ are nucleation model parameters[30]. t e the time over which we allow our simulations to equilibrate intothe steady-state before taking size distribution measurements.T P v γ ρ b r η ξ t e [K] [dyn/cm ] [dyn/cm] [g/cm ] [10 − cm] [ns]300 8 . · .
93 6.05 8.95 25325 4 . · .
94 5.29 7.47 20350 1 . · .
95 4.63 6.31 10375 4 . · .
97 4.02 5.38 6
B. Simulation Analysis
We use the simple Stillinger criterion[54] (also known asthe friends-of-friends method) to identify clusters. As thesimulation runs, the cluster size distribution is regularlycalculated and outputted, typically resulting in > ∼ × molecules. We perform this at the endof the simulation, well-within the steady-state nucleationregime. This calls for per-atom outputs of velocity and po-sition information. Sections VIII and IX detail how thedensity profile and temperature measurements respectivelyare made, and discuss the results. III. NUCLEATION RATES
We use a modified Yasuoka-Matsumoto method[35](threshold method) to measure nucleation rates. In thesteady-state nucleation regime, the time rate of increase ofthe number of clusters above a certain size N is the nu-cleation rate. However, the simulations must equilibrate -form the sub-critical size distribution - and properly popu-late it before they reach the steady state nucleation regime.How long the simulations take to transition into the steadystate regime is not known a priori . Nucleation rates esti-mated from the first nucleation event alone (e.g. mean firstpassage time or survival probability methods) can be ordersof magnitude smaller than the true steady state nucleationrates[57, 58]. Thus we use the following method: To thesize-threshold curves, we fit the following function, whichis able to capture the transition from the equilibration tothe steady-state phase, N ( > i ) = J · N ( > i, t ) + N ( > i, , (1)where N ( > i, t ) = (cid:20) π t − t t r (cid:21) · t − t π · V, (2)where J is the nucleation rate, t is the lag time, t r isthe relaxation timescale and V the volume of the simula-tion box. This function captures the system’s transitionfrom the initial equilibration phase to the intermediate re-laxation phase as the clusters begin to form, through tothe steady-state regime. We count clusters above a certainpost-critical size N , frequently throughout the simulation,and fit the count to this curve, allowing J , t , and t r tovary. Visual inspection of this approach is provided in theupper left panels of figures 2 for run T375c respectively.Our nucleation rate measurements are listed in II. Figure3 plots our simulations’ nucleation rates against supersatu-ration, and includes comparison to earlier results[1], whichused smaller simulations and were therefore restricted tolower nucleation rates. Estimates for the critical clustersizes using the first nucleation theorem, via i ∗ NT = (cid:18) ∂ ln J∂ ln S (cid:19) T − • High temperature (325 K, 350 K, 375 K): Runs atthese temperatures have nucleation rates in the range ∼ − cm − s − . These runs have generally lowerrors on the nucleation rates. For the higher nucle-ation rates, there is an error on the supersaturation,as the pressure drops significantly due to the largenumber of clusters forming quickly. • Low temperature (300 K): Here we measure nucle-ation rates in the range ∼ − cm − s − . Theseruns suffer from extremely long equilibration periods,3 ABLE II. Run temperature T , supersaturation S as calculated from the run monomer number density, box length L , moleculenumber N , runtime t end , nucleation rate measured from simulation J MD , critical cluster size i ∗ from the first nucleation theorem, J SP Semi-phenomenological model prediction, J MCNT modified classical nucleation theory prediction, i ∗ ∆ G critical size from the ∆ G reconstruction.Run ID T S L N t end J MD i ∗ α J SP J MCNT i ∗ ∆ G [K] [nm] [ · ] [ns] [cm − s − ] [cm − s − ] [cm − s − ]T300a 300 23 . ± .
89 4859 . . . ± . · - 1 .
31 4 . · . · . ± .
71 6581 . . ± . · - 1 .
11 1 . · . · . ± .
94 7591 . . ± . · - 0 .
59 4 . · . · . ± .
46 3307 . . . ± . · - 0 .
52 1 . · . · . ± .
33 3441 . . . ± . · ± .
59 1 . · . · . ± .
20 4863 . . . ± . · ± .
26 3 . · . · . ± .
06 6461 . . . ± . · ± .
25 2 . · . · . ± .
03 8167 . . . ± . · ± .
15 1 . · . · . ± . . . . ± . · ± .
16 8 . · . · -T350a 350 4 . ± .
15 2573 . . . ± . · - 0 .
38 8 . · . · . ± .
09 2680 . . . ± . · ± .
23 2 . · . · . ± .
001 2783 . . . ± . · ± .
21 2 . · . · . ± .
003 2893 . . . ± . · ± .
18 3 . · . · . ± .
001 4993 . . . ± . · ± .
06 3 . · . · -T375a 375 3 . ± .
06 2017 . . . ± . · ±
10 0 .
40 7 . · . · . ± .
005 2107 . . . ± . · ± .
26 1 . · . · . ± .
001 2158 . . . ± . · ± .
23 2 . · . · . ± .
001 2937 . . . ± . · ± .
25 7 . · . · - which continue while the initial large, stable clustersare already forming. In other words, the sub-criticaldistribution formation timescale t r is longer than thenucleation timescale 1 / ( J · V ). This leads to largeerrors in both the nucleation rate measurements andthe supersaturation measurements. IV. RATE COMPARISON WITH ANALYTICALMODELS
Nucleation models endeavor to describe the phase changeprocess in purely thermodynamic terms. The standard ap-proach tries to find the balance between the Gibbs freeenergy gain and cost due to the creation of volume andsurface. The classical nucleation theory (CNT)[55, 59–64]is the most basic of them all, and forms the basis uponwhich many appendages have since been added. In theCNT, the surface energy term in the Gibbs free energy issimply calculated using the planar surface tension, with noadditional corrections. The CNT nucleation rate is [65] J CNT = (cid:114) πγ m r (cid:18) p g k b T (cid:19) exp (cid:34) π r γ
27 ( k B T ) (log S ) (cid:35) , (4)where m is the molecular mass, γ the planar surface tensionat the run temperature, p g the monomer partial pressurein the simulation box (assuming an ideal gas) gas pressure and S the supersaturation S = p g p v , (5)where p v is the equilibrium vapor pressure at the run tem-perature. r is a characteristic molecular radius: r = (cid:18) ρ l π (cid:19) / , (6)where ρ l is the bulk liquid density at the run temperature.Table I includes the thermodynamic variables for SPC/Ewater, which we use in our analysis and comparison to nu-cleation models. The CNT predictions for the nucleationrates of SPC/E water at our runs’ supersaturations areshown as solid curves in figure 3. We find that the CNTpredicts too-high nucleation rates by factors of 10 − .Various authors [6, 8, 66] employ a 2-parameter, temper-ature dependent correction factor J corr = J CNT exp (cid:18) A + BT (cid:19) . (7)The Manka et al. (2010)[6] (see their figure 4) lami-nar flow diffusion chamber experiments find that the pa-rameter pair ( A, B ) = ( − . , − B/A = 235 . S ), but less strongly than in CNT. These4 c o un t > N N = 25, J = (3.38±0.04)e22N = 50, J = (2.90±0.03)e22 N = 90, J = (2.58±0.02)e22J-rate units in cm −3 s −1
10 20 30 40 50time [nano-sec ]0.5 1.0 1.5 2.0 2.5timesteps 1e710 -6 -5 -4 -3 -2 -1 i m e r - m o n o m e r r a t i o i=2i=3 i=4i=5 i=6i=7 i=8i=9 p r e ss u r e [ a t m ]
10 20 30 40 50time [nano-sec ]0.5 1.0 1.5 2.0 2.5timesteps 1e72.82.93.03.13.23.33.43.5 S a t u r a t i o n i α=0.183 FIG. 2. For the 768000 molecule simulation T375c, clockwise from the upper left panel: (1) Nucleation rate measurements, (2)largest-cluster growth curve, (3) number density and monomer partial pressure, and (4) i -mer concentrations, all over the entire runperiod. Nucleation rates are measured by counting the number of clusters above a specified threshold size, at periodic time intervals.The steady-state regime slope is the nucleation rate. (Refer to section II B and equation (1).) The dotted vertical lines indicate thelag times for each size. The cluster sticking probabilities α are calculated from the measured slopes, di/dt , using equation (10). Wecan see that this run took ∼
20 ns to equilibrate, and spent another ∼ −
15 ns in the lag phase before reaching the steady-stateregime (for the chosen threshold sizes of N = 50 and N = 90). The probability that a cluster-monomer encounter results in thecluster growing by one molecule is the sticking probability α = 0 . . corrected CNT predictions, when extrapolated to our su-persaturations (dashed curves in figure 3) under-predict ourmeasurements by 1-3 orders of magnitude. Using our dataat temperatures T = 325 , ,
375 K to determine the best-fit parameter pair, we find (
A, B ) = ( − . , − B/A = 297 . − ). Figure 4 shows the ratio betweenthe MCNT model predictions (red markers) and the directMD measurements. Refer to table II for the MCNT modelnucleation rate predictions.The Semi-Phenomenological model (SP)[56, 69–72] at-taches a ∼ /R (or ∼ i − / ) correction to the surface ten- sion, where R is the cluster size under the assumption ofsphericity. This radial dependence is functionally equiva-lent to that introduced by the Tolman length[11, 73, 74],although the motivation is different: the coefficient to thisterm is set by the second virial coefficient B [34] so thatthe dimer number density is correctly predicted. The nu-cleation rate predictions for the SP model relative to themeasured values are plotted with red markers in figure 4.The predictions at T = 300 K are somewhat accurate -within a factor of 5 of the measured values, although, asnoted in section III the measurements at these tempera-tures carry significant uncertainty in the nucleation rate.At the higher temperatures, the SP model under-predictsthe measured MD rates by factors of 4 −
80. Table II liststhe SP model nucleation rate predictions.
V. NUCLEATION RATE SCALING
In this section we examine the scaling of the nucleationrates. For the case of water, Hale (2005)[75] (and simi-larly for Lennard-Jones in Hale (2010)[76]) uses a scalingrelation[77] for experimentally measured nucleation ratesover the range J = 10 − cm − s − of5 S10 J [ c m − s − ]
76 ±7412 ±617 ±519 ±521 ±222 ±976 ±6722 ±616 ±426 ±141 ±223 ±826 ±635 ±337 ±7
Tanaka et al. 2014T = 300KT = 325KT = 350KT = 375KManka experimentsCNTbest shift
FIG. 3. The filled solid markers are the nucleation rates wemeasure from our simulations. The solid curves correspond tothe classical nucleation theory (4) predictions. The dashed curveincludes the CNT correction factor (7) used in Manka et al.(2011)[6], and the dotted one is our best-fit correction factor(the fit excludes the runs at T = 300 K). ln S ( T c /T − . . (8)Tanaka et al. (2014)[31] showed that this scaling relationworks well for large scale Lennard-Jones simulations andArgon laboratory experiments, albeit with an exponent of1.3 instead of 1.5. We confirm that the same scaling rela-tion (8) applies well to our SPC/E water nucleation ratemeasurements. However, we find that the combined nucle-ation rates from both SPC/E simulations and laboratoryexperiments with water are even better scaled byln S ( T c /T − . . (9)Figure 5 shows the nucleation rates as a function of (9).This empirical scaling relation seems to work well over asurprisingly wide nucleation rate range - from J = 10 − to J = 10 cm − s − for both MD simulations and experi-ments. The results from the MD simulations join smoothlywith the experiments with the scaling by ln S/ ( T c /T − . .Figure 5 also shows, using solid curves, the nucleation ratespredicted by the SP model for various temperatures. Thisscaling relation also works very well for the SP model.
300 35010 −2 J m ode l / J M D T [K]
MCNT newMCNT previousSP newSP previous
FIG. 4. A comparison between the measured molecular dynam-ics nucleation rates, and those predicted by the analytic SP andMCNT nucleation models. ‘new’ and ‘previous’ refer to thesimulations detailed in this paper, and those by Tanaka et al.(2014)[1] respectively. While the SP model is more reliable thanthe MCNT, the SP model shows a trend of rate under-predictionfor lower supersaturations. ln(S)/(T c /T−1) J [ c m − s − ] previous MDT=375KWolk & Strey 2001Miller et al. 1983Khan et al. 2003Manka et al. 2010Brus et al. 2009Kim et al. 2004T=350KT=325KT=300K FIG. 5. Nucleation rates as function of ln S/ ( T /T c − . for MDsimulations and experiments[1, 3, 5, 6, 8, 13, 14]. Our simula-tions are filled circles and the previous simulations ‘+’ symbols.The solid curves show the SP model for various temperatures(210, 315, 350, and 380 K). For thermodynamic quantities suchas the surface tension and the saturated vapor pressure, we usethose of the SPC/E at 315, 350, and 380 K for the comparisonwith the MD results, while real water at 210 K (pink curve) forthe comparison with the experiment. I. STICKING PROBABILITIES
The sticking probabilities α can be calculated from therate at which large, stable clusters grow. For each run,we observe the size of the largest cluster, and measureits growth rate di/dt over the second half of the simula-tion. Early on in the simulations, before stable clustershave formed, the largest designation jumps between clus-ters. However, the first stable cluster to form is likely to re-main the largest until the end of the simulation. The upperright panel of figure 2 shows our cluster growth rate mea-surements for run T375c. We find the cluster size i ( t ) to bestrongly cubic within the steady-state regime. The clustergrowth rate is therefore proportional to the surface area.This is consistent with what has been found in Lennard-Jones nucleation simulations[29, 30]. We may determine α [1, 30] from α = 34 πr v th n (1) (cid:18) − S (cid:19) − di / dt . (10)The supersaturation S dependence includes the effect ofthe evaporation of molecules from the clusters into thegas. We list the measured sticking probability results intable II. Sticking probability results for our low temper-ature T = 300 K runs are somewhat unreliable and canexceed unity due to the large number of dimers, trimers,and tetrames which also contribute to cluster growth. Eq(10) considers the accretion and evaporation of monomersonly. Our sticking probability measurements are consistentwith those measured at slightly higher supersaturations inTanaka et al. (2014) [1]. The upper panel of figure 6 plots α against S . The sticking probability is a necessary pre-requisite for performing the ∆ G landscape reconstructionprocedure for post-critical clusters (see section VII).While we are, due to computational constraints, unableto probe the low nucleation rates observed in laboratoryexperiments, it is possible to measure cluster growth ratesunder laboratory conditions. We have performed an addi-tional simulation from the end state of T325c, in which wemeasured a sticking probability α = 0 .
26. We target thetemperature and saturation conditions found in Brus et al.(2008)[22]. Using a Nose-Hoover thermostat we maintainthe temperature, and gently increase the box size until thesupersaturation S = 2 .
5, after which we continue runningfor 60 ns. Under these low pressure conditions, no new clus-ters nucleate (Brus et al. (2008)[22] report nucleation rates ∼ cm − s − ) due to our comparatively small and short-lived system. However, clusters which had previously nu-cleated and then grown under the original conditions, per-sist. Using the largest of these still-post critical clusters,we measure a decreased growth rate: an i / slope shal-lower by a factor of ∼
7. Including this, and the reduced(by a factor ∼
3) monomer number density into Eq. (10)gives a sticking probability for these laboratory-like growthrates of α = 0 .
15. In nucleation models the sticking effi-ciency is usually taken to be unity, entering linearly in thetransition growth rate (typically denoted R + ), as a prefac-tor to the ∆ G i exponent. We find that in the T = 325 K −1 −1 ln S (T/273K) ln S αα T=300K325350375
FIG. 6. Upper panel: sticking probability measurements forour simulations (solid markers) and from previous simulations(‘+’ symbols)[1]. Our simulations continue the expected trendof lower growth rates with decreasing gas pressure. The lowerpanel shows a temperature-dependent scaling relation which re-duces the results to a single curve. The blue donut markerindicates the sticking probability measured under experimental saturations (refer to the end of section VI). and S = 2 . VII. FREE ENERGY RECONSTRUCTION
In this section, we evaluate the formation free energyof a cluster ∆ G i ( S ) directly from our molecular dynamicssimulations, even for post-critical cluster sizes. We obtain∆ G i ( S ) from the equilibrium size distribution of the clus-ter. The equilibrium size distribution can be obtained us-ing the steady state size distribution, the accretion rate ofmolecule on a cluster, and the nucleation rate, all of whichcan be measured directly in from the MD simulations. Re-fer to Tanaka et al. (2014)[31] for a thorough explanationof the technique. The cluster size distributions are mea-sured in the MD simulations and time-averaged over thesteady state nucleation phase. In the accretion rate, weuse the value of the sticking probability obtained from MDsimulations. With the use of them, we reconstruct the fullequilibrium size distribution (at all sizes i , where we havegood abundance estimates, including i >> i ∗ ) and then the7ntire free energy function ∆ G i ( S = 1). We can further de-rive ∆ G i ( S = 1), which is a surface term corresponding tothe work required to form the vapor-liquid interface, bysubtracting the volume term from ∆ G i ( S ):∆ G i ( S = 1) = ∆ G i ( S ) + ( i −
1) ln S. (11)Figure 7 shows ∆ G ( S = 1) obtained from the MD resultsat 375 K, and various supersaturations. Since ∆ G ( S =1) is supersaturation independent the values from all runsshould overlap at all sizes. However, our simulation data isonly good enough for accurate abundance estimates belowa certain cluster size, which depends on the run properties.The highest nucleation rates run of these (T375a) produceda large number of clusters over the entire plotted size rangeand allows the most reliable reconstruction of ∆ G i ( S ). Theresults from the other runs are only accurate at smallersizes, where they overlap with (T375a). Figure 7 also showsthe surface energy ∆ G i ( S = 1) divided by that of the CNT,i.e., ∆ G i ( S = 1) / ( ηi / kT ). In the figure, we also show theresults of the SP model, given by∆ G i ( S = 1) / ( ηi / kT ) = 1+( ξ/η ) i − / − ( ξ/η ) i − / . (12)The simulation results deviate from the SP model at 375 K.In Figure 7, we can fit the reconstructed ∆ G i ( S = 1) with∆ G i ( S = 1) / ( ηi / kT ) = 1 + Ai − / − Ai − / , (13)using a fitting parameter A = 0 . A =0.9, 1.0, 1.0 and 1.5 at 375, 350, 325, and 300 K, re-spectively) and the MD simulations for two cases: one inwhich α = 1 and the other in which α is set to be valueobtained directly from simulation. In Figure 8, the predic-tions from the SP model are also shown for comparison.We find the new model agrees with the simulations withinone order of magnitude for all cases. At 375 K this is nosurprise, since this data was used to to determine the pa-rameters of our fitting function for the surface term (13).The good agreement at the other temperatures is encourag-ing and might motivate using (13) also to predict nucleationrates at different temperatures and supersaturations. VIII. CLUSTER DENSITIES
It has been shown[78] that for spherical clusters, liquid-vapor interface densities are well-approximated by ρ ( r ) = 12 (cid:20) ρ c + ρ g − ( ρ c − ρ g ) tanh (cid:18) r − Rd (cid:19)(cid:21) , (14)where ρ c is the number density within the cluster, ρ g the gasnumber density, R the interface position, and d its width.In each cluster’s center-of-mass frame, we bin the spheri-cal number density, using a bin size of 1 . . The numberdensity profiles for clusters of the same size are used tomake ensemble averages, to which equation (14) can befit. This method of measuring internal cluster densities
T=375K i ∆ G ( S = ) / k T i −1/3 ∆ G ( S = ) / ( η i / k T ) t375at375bt375ct375d t375at375bt375ct375d FIG. 7. Reconstructed Gibbs free energy curves shifted to S = 1(top panel). Runs at T = 375 K at different supersaturationswere used. To get to S = 1 the CNT volume term was sub-tracted, the resulting ∆ G i ( S = 1) can be interpreted as thesurface term. The bottom panel shows the ∆ G i ( S = 1) dividedby the surface term from CNT. The surface term from the SPmodel and a simple fitting function (Eq. (13)) are shown withdotted and solid lines. is robust only for clusters which are large enough to pos-sess a constant density core. Clusters with i < − ρ c against R for clusters in T325f. We observean over-density for clusters between 4 − . ∼
00 35010 −2 J m ode l / J M D T [K]
FIG. 8. Comparisons of nucleation rate from the MD simula-tions and several model predictions: the SP model (grey filledcircles), our new surface term fit (Eq. (13)) with α = 1 (filledcircles) and using the α values measured in the MD simulations(open circles). clusters with i = i ∗ implied larger surface areas, and there-fore larger-than-expected surface energies, resulting in anincreased free energy cost to form a critical cluster, whichlowered nucleation rates from model predictions. We sus-pect that nucleation rate predictions are more successful forSPC/E water than they are for Lennard-Jones because theassumption of small clusters possessing the bulk density for i = i ∗ is more realistic for the case of SPC/E water. IX. TEMPERATURES
We define the temperature of an ensemble of atoms fromtheir mean kinetic energy kT ≡ (cid:104) E kinetic (cid:105) = 13 N N (cid:88) i =1 mv i . (15)Using full per-particle velocity information outputted at theend of the simulation, we are able to investigate the clustersize dependence of temperature. We find that sub-criticalclusters are at the run average temperature, as observedin similar Lennard-Jones simulations[26]. However, con-trary to what has been observed in Lennard-Jones nucle-ation simulations, post-critical SPC/E water clusters pos-sess temperatures consistent with the run average temper-ature. Figure 9 plots the ensemble average (at each clustersize i ) of their temperatures against the density profile in-terface midpoint R (i.e., the cluster radius).The latent heat from condensation has been efficientlydissipated back into the gas, leaving the post-critical clus- T − [ K ] run temperature t325f FIG. 9. Cluster temperatures for T325f. We observe clustertemperatures consistent with the run target temperature. Thereis no signal of residual latent heat for post-critical clusters, con-trary to what has been observed in Lennard-Jones nucleationsimulations[26]. We suspect this may be due to the long-rangenature of the SPC/E interaction potential, which enables effi-cient energy exchange between members of the cluster and thegas. ρ c [ n m ] bulk liquid t325f FIG. 10. Cluster densities for T325f, plotted against interfacepositions, both determined from density profile fits to (14). Weobserve an over-density for clusters between 4 − . < . ters in thermodynamic equilibrium with their surroundings.This finding is consistent with the post-critical clusters den-sity profile measurements, which finds their densities at theexpected bulk density. Had the clusters significant latentheat retention, their densities would be correspondinglylower. We conjecture that the efficient kinetic energy ex-change is effected by the long-range Coulombic interactions- even molecules deep within a cluster may exchange energyand angular momentum with members of the gas - resultingin kinetic energy equipartition on shorter timescales thancluster growth rates. A molecule impinging on a cluster im-parts heat onto into the cluster-system, yet the heat doesnot linger. This may have implications for non-isothermalnucleation models[79], which include latent heat retentionin the thermodynamic description of growing droplets.9 . CONCLUSIONS We have performed molecular dynamics simulations ofSPC/E water, and significantly closed the nucleation-rategap between simulation and experiment, measuring nucle-ation rates low as ∼ − cm − s − . This is an hithertounexplored saturation and temperature region for waternucleation experiments and water nucleation simulation.Nucleation rate results in this new regime will provide mod-els with further testing comparison opportunities, to com-plement the already-existing lower nucleation rates fromexperiment, and higher nucleation rates from other simu-lations. We summarize our most significant contributionsbelow. • We introduce a new functional form, Eq. (1) in or-der to implement the Yasuoka-Matsumoto nucleationrate measurement. This modified version smoothlycaptures the system’s transition between the lagphase, relaxation phase, and onto the steady-stateregime. • As expected, the CNT over-estimates nucleation ratesby a few orders of magnitude. The empirical CNTcorrection factor (7) [66], when using the Mankaet al.[6] best-fit parameter values (calibrated in thelow nucleation rate, low saturation regime) under-estimates our rates by a few orders of magnitude.When fitting their proposed correction function toour results, we find that the slopes are not steepenough. We conclude that this empirical and purelytemperature-dependent correction factor to the CNTis not rich enough to reproduce the qualitative be-havior we observe in our regime. • The MCNT nucleation model continues to over-predict nucleation rates, by factors of up to 10 . TheSP model on the other hand, does somewhat bet-ter, under-predicting rates at worst by a factor of24. Despite these failings, we note that these modelpredictions are significantly more accurate than thecorresponding predictions for the Lennard-Jones fluidvapor-to-liquid nucleation[26, 29–32]. • Performing a cluster growth rate measurement sim-ulation under laboratory conditions (those found inBrus et al. (2008)[22]) of T = 325 K and S = 2 .
5, wemeasure a sticking probability of α = 0 .
15. This sug-gests that in this regime, nucleation rate predictionsfrom models should lowered by a factor of seven. • We find the cluster size i ( t ) to be strongly cubicwithin the steady-state regime. The cluster growthrate is therefore proportional to the surface area, aresult new to water nucleation. This is consistentwith what has been found in Lennard-Jones nucle-ation simulations[29, 30]. • Unlike Lennard-Jones nucleation simulations, we findthat post-critical clusters have temperatures consis-tent with the simulation average temperature: Grow-ing clusters are in thermal equilibrium with their sur-roundings. Latent heat is not retained as the clustersgrow, it is efficiently dissipated back into the gas.We suspect this efficiency is due to the long-rangeCoulombic interactions, not present in the Lennard-Jones case. This could have an impact on nucleationmodels which include non-isothermal processes intothe thermophysical modeling of cluster properties[79]. • Post critical clusters have densities consistent withwhat is expected from the bulk liquid. There isa possible indication of an over-density for clustersaround the critical size. This would imply a lower-than-expected surface area, which lowers the totalsurface energy, decreasing the free energy cost to forma critically-sized cluster and would result in higher-than-expected nucleation rates. • The scaling relation ln S/ ( T /T c − . is remarkablysuccessful in reducing the 3-parameter T vs. S vs. J surface into a 2-parameter curve. It accurately linksnucleation rates from simulation and experiment fromover 30 orders of magnitude in the nucleation raterange, and a temperature range of 180 K. [1] K. K. Tanaka, A. Kawano, and H. Tanaka, The Journal ofChemical Physics , 114302 (2014).[2] Duska, Michal, Nemec, Tomas, Hruby, Jan, Vins, Vaclav,and Plankova, Barbora, EPJ Web of Conferences , 02013(2015).[3] Y. J. Kim, B. E. Wyslouzil, G. Wilemski, J. W¨olk, andR. Strey, The Journal of Physical Chemistry A , 4365(2004), http://dx.doi.org/10.1021/jp037030j.[4] M. A. L. J. Fransen, J. Hruby, D. M. J. Smeulders, andM. E. H. van Dongen, The Journal of Chemical Physics , 164307 (2015).[5] A. Khan, C. H. Heath, U. M. Dieregsweiler, B. E. Wys-louzil, and R. Strey, The Journal of Chemical Physics ,3138 (2003). [6] A. A. Manka, D. Brus, A.-P. Hyvrinen, H. Lihavainen,J. Wolk, and R. Strey, The Journal of Chemical Physics , 244505 (2010).[7] V. B. Mikheev, P. M. Irving, N. S. Laulainen, S. E. Barlow,and V. V. Pervukhin, The Journal of Chemical Physics ,10772 (2002).[8] J. W¨olk and R. Strey, The Journal of Physical Chemistry B , 11683 (2001), http://dx.doi.org/10.1021/jp0115805.[9] Y. Viisanen, R. Strey, and H. Reiss, The Journal of Chem-ical Physics , 8205 (2000).[10] Y. Viisanen, R. Strey, and H. Reiss, The Journal of Chem-ical Physics , 4680 (1993).[11] V. Holten, D. G. Labetski, and M. E. H. van Dongen, TheJournal of Chemical Physics , 104505 (2005).
12] C. C. M. Luijten, K. J. Bosschaart, and M. E. H. vanDongen, The Journal of Chemical Physics , 8116 (1997).[13] R. C. Miller, R. J. Anderson, J. L. Kassner, and D. E.Hagen, The Journal of Chemical Physics , 3204 (1983).[14] D. Brus, V. Zdimal, and J. Smolik, The Journal of Chem-ical Physics , 174501 (2008).[15] R. C. Miller, R. J. Anderson, J. L. Kassner, and D. E.Hagen, The Journal of Chemical Physics , 3204 (1983).[16] Y. Viisanen, R. Strey, and H. Reiss, The Journal of Chem-ical Physics , 4680 (1993).[17] C. C. M. Luijten, K. J. Bosschaart, and M. E. H. vanDongen, The Journal of Chemical Physics , 8116 (1997).[18] M. P. Anisimov, P. K. Hopke, I. N. Shaimordanov, S. D.Shandakov, and L.-E. Magnusson, The Journal of ChemicalPhysics , 810 (2001).[19] J. W¨olk and R. Strey, The Journal of Physical Chemistry B , 11683 (2001), http://dx.doi.org/10.1021/jp0115805.[20] Y. J. Kim, B. E. Wyslouzil, G. Wilemski, J. W¨olk, andR. Strey, The Journal of Physical Chemistry A , 4365(2004), http://dx.doi.org/10.1021/jp037030j.[21] V. Holten, D. G. Labetski, and M. E. H. van Dongen, TheJournal of Chemical Physics , 104505 (2005).[22] D. Brus, V. ˇZd´ımal, and J. Smol´ık, The Journal of Chem-ical Physics , 174501 (2008).[23] D. Brus, V. ˇZd´ımal, and H. Uchtmann, The Journal ofChemical Physics , 074507 (2009).[24] M. Schweizer and L. M. C. Sagis, The Journal of ChemicalPhysics , 224102 (2014).[25] M. Horsch, J. Vrabec, and H. Hasse, Phys. Rev. E ,011603 (2008).[26] Angelil, Raymond and Diemand, J¨urg and Tanaka, KyokoK. and Tanaka, Hidekazu, The Journal of Chemical Physics , 074303 (2014).[27] M. Horsch, H. Hasse, A. K. Shchekin, A. Agarwal, S. Eck-elsbach, J. Vrabec, E. A. M¨uller, and G. Jackson, Phys.Rev. E , 031605 (2012).[28] Yasuoka, Kenji and Matsumoto, Mitsuhiro, The Journal ofChemical Physics , 8451 (1998).[29] K. K. Tanaka, H. Tanaka, T. Yamamoto, and K. Kawa-mura, The Journal of Chemical Physics , 204313 (2011).[30] Diemand, J¨urg and Angelil, Raymond and Tanaka, KyokoK. and Tanaka, Hidekazu, The Journal of Chemical Physics , 074309 (2013).[31] K. K. Tanaka, J. Diemand, R. Ang´elil, and H. Tanaka, TheJournal of Chemical Physics , 194310 (2014).[32] K. K. Tanaka, K. Kawamura, H. Tanaka, and K. Nakazawa,The Journal of Chemical Physics , 184514 (2005).[33] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma,The Journal of Physical Chemistry , 6269 (1987),http://dx.doi.org/10.1021/j100308a038.[34] H. Matsubara, T. Koishi, T. Ebisuzaki, and K. Yasuoka,The Journal of Chemical Physics , 214507 (2007).[35] K. Yasuoka and M. Matsumoto, The Journal of ChemicalPhysics , 8463 (1998).[36] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W.Impey, and M. L. Klein, The Journal of Chemical Physics , 926 (1983).[37] J. L. F. Abascal and C. Vega, The Journal of ChemicalPhysics , 234505 (2005).[38] V. Molinero and E. B. Moore, The Journal of PhysicalChemistry B , 4008 (2009).[39] M. H. Factorovich, V. Molinero, and D. A. Scherlis, Journalof the American Chemical Society , 4508 (2014).[40] V. Molinero, AIP Conference Proceedings , 83 (2013). [41] B. Song and V. Molinero, The Journal of Chemical Physics , 054511 (2013).[42] V. Holten, D. T. Limmer, V. Molinero, and M. A. Anisi-mov, The Journal of Chemical Physics , 174501 (2013).[43] M. H. Factorovich, V. Molinero, and D. A. Scherlis, TheJournal of Chemical Physics , 064111 (2014).[44] T. Li, D. Donadio, G. Russo, and G. Galli, Phys. Chem.Chem. Phys. , 19807 (2011).[45] A. V. Mokshin and B. N. Galimzyanov, The Journal ofPhysical Chemistry B , 11959 (2012), pMID: 22957738.[46] F. Zipoli, T. Laino, S. Stolz, E. Martin, C. Winkelmann,and A. Curioni, The Journal of Chemical Physics ,094501 (2013).[47] S. Plimpton, Journal of Computational Physics , 1(1995).[48] R. W. Hockney and J. W. Eastwood, Computer SimulationUsing Particles (Taylor & Francis, Inc., Bristol, PA, USA,1988).[49] E. Pollock and J. Glosli, Computer Physics Communica-tions , 93 (1996).[50] J. P. Ryckaert, G. Ciccotti, and J. Berendsen, The Journalof Computational Physics , 327 (1977).[51] S. Nos´e, The Journal of Chemical Physics , 511 (1984).[52] W. G. Hoover, Phys. Rev. A , 1695 (1985).[53] W. Shinoda, M. Shiga, and M. Mikami, Phys. Rev. B ,134103 (2004).[54] F. H. Stillinger, The Journal of Chemical Physics , 1486(1963).[55] V. I. Kalikmanov, Nucleation Theory , Lecture Notes inPhysics (Springer, Dordrecht, 2013).[56] D. W. Oxtoby, Journal of Physics: Condensed Matter ,7627 (1992).[57] V. A. Shneidman, The Journal of Chemical Physics ,051101 (2014).[58] A. V. Mokshin and B. N. Galimzyanov, The Journal ofChemical Physics , 024104 (2014).[59] M. Volmer and W. A., Z. Phys. Chem. (1926).[60] R. Becker and W. Dring, Annalen der Physik , 719(1935).[61] J. Zeldovich, J. Exp. Theor. Phys. (1942).[62] J. Feder, K. Russell, J. Lothe, andG. Pound, Advances in Physics , 111 (1966),http://dx.doi.org/10.1080/00018736600101264.[63] D. Kashchiev, Nucleation (Butterworth-Heinemann, 2000).[64] I. Ford, Physical Review E , 5615 (1997).[65] R. Becker and W. D¨oring, Annalen der Physik , 719(1935).[66] J. Wolk, R. Strey, C. H. Heath, and B. E. Wyslouzil, TheJournal of Chemical Physics , 4954 (2002).[67] Y. Viisanen, R. Strey, and H. Reiss, The Journal of Chem-ical Physics , 4680 (1993).[68] Y. J. Kim, B. E. Wyslouzil, G. Wilemski, J. Wlk, andR. Strey, The Journal of Physical Chemistry A , 4365(2004), http://dx.doi.org/10.1021/jp037030j.[69] A. Laaksonen, I. J. Ford, and M. Kulmala, Phys. Rev. E , 5517 (1994).[70] A. Dillmann and G. E. A. Meier, The Journal of ChemicalPhysics , 3872 (1991).[71] C. F. Delale and G. E. A. Meier, The Journal of ChemicalPhysics , 9850 (1993).[72] V. I. Kalikmanov and M. E. H. van Dongen, The Journalof Chemical Physics , 4250 (1995).[73] R. C. Tolman, The journal of chemical physics , 333(1949).