In this work, we use large-scale molecular dynamics simulations coupled to free energy calculations to identify for the first time a limit of stability (spinodal) and a change in the nucleation mechanism in aqueous NaCl solutions. This is a system of considerable atmospheric, geological, and technical significance. We find that the supersaturated metastable NaCl solution reaches its limit of stability at sufficiently high salt concentrations, as indicated by the composition dependence of the salt chemical potential, indicating the transition to a phase separation by spinodal decomposition. However, the metastability limit of the NaCl solution does not correspond to spinodal decomposition with respect to crystallization. We find that beyond this spinodal, a liquid/amorphous separation occurs in the aqueous solution, whereby the ions first form disordered clusters. We term these clusters as “amorphous salt.” We also identify a transition from one- to two-step crystallization mechanism driven by a spinodal. In particular, crystallization from aqueous NaCl solution beyond the spinodal is a two-step process, in which the ions first phase-separate into disordered amorphous salt clusters, followed by the crystallization of ions in the amorphous salt phase. By contrast, in the aqueous NaCl solution at concentrations lower than the spinodal, crystallization occurs via a one-step process as the ions aggregate directly into crystalline nuclei. The change of mechanism with increasing supersaturation underscores the importance of an accurate determination of the driving force for phase separation. The study has broader implications on the mechanism for nucleation of crystals from solutions at high supersaturations.

Crystallization from solution is a fundamental process that is of interest in many fields, including but not limited to atmospheric sciences, geochemistry, and biology.1–4 Despite its importance, the nature of the initial, rate-determining and highly non-equilibrium nucleation process, entailing the formation of microscopic ordered precursors (nuclei) of the stable crystal phase, has not been fully revealed at the molecular level, partly due to the fact that the existing experimental techniques lack the spatiotemporal resolution to probe the short-lived nanometer-scale nuclei at the early stages of the nucleation process. Molecular simulations, on the other hand, do not suffer from the lack of temporal or spatial resolution, and the last decade has seen a rapid growth in the number of simulation studies for both homogeneous and heterogeneous nucleation in various systems, ranging from simple models [e.g., Lennard-Jones (LJ) particles]5,6 to realistic systems represented by molecular force fields (e.g., NaCl, CaCO3, and urea).7–18 Many important features of the microscopic mechanisms underlying the nucleation/crystallization of solute molecules or ions from solution have been illustrated. In spite of the rapid progress in understanding the underlying mechanisms using simulations, the accurate calculation of rates remains a challenging task as the calculation is usually subject to large uncertainties,19 entails demanding computations,12 and is often limited to a narrow range of supersaturations.

In general, simulations of nucleation/crystallization from solution are conducted at high solute concentrations, where fast nucleation may be observed within a reasonable amount of simulation time.15,20,21 Despite the invaluable insights provided by many prior simulation studies, a precise determination of the driving force for nucleation, i.e., the solution supersaturation, or more precisely, the difference between chemical potential of the solute in solution and crystal phases, is generally missing in most prior studies. The supersaturation of a simulated solution is often calculated using the experimental solute solubility;14,16 however, the actual solubility of the underlying molecular force field may differ significantly from the experimental value, leading to severe under-/overestimation of the driving force. For the few studies that estimated the solubility of the underlying molecular force fields,11,22 the chemical potentials of the solute relative to the crystal at different supersaturations were either not known or not calculated with sufficient accuracy. For example, the calculation of solute chemical potentials in Ref. 11 was based on the assumption of an ideal solution with unitary activity coefficients. Since the nucleation process is highly sensitive to the driving force, the interpretation of nucleation mechanisms and rates from prior simulations becomes difficult when the supersaturation or the chemical potential of the solution or crystal phases is not properly established.

The connection between crystallization and liquid-liquid phase separation has been explored previously, both experimentally6,23–31 and by simulations.17,20 In protein and colloidal solutions,6,26,27 the liquid-liquid phase transition becomes metastable with respect to the liquid-solid transition when attractive interactions are sufficiently short-ranged. Experiments using polymer blends23 that can undergo both liquid-liquid phase separation and crystallization as a function of temperature show an increase in the nucleation rate of crystallization because of the concentration fluctuations caused by spinodal decomposition of two liquids at higher temperatures. Multivalent salts that tend to form hydrated salt phases, such as MgSO4 and CaCO3, have been experimentally observed30,31 to form ion-rich liquid phases in supersaturated solutions; a recent study based on molecular dynamics (MD) simulations of such solutions17 proposed a possible mechanism for nucleation proceeding via the formation of “hydrated clusters.” The study hypothesized the existence of an underlying metastable liquid-liquid transition interfering with nucleation at specific temperatures but retained the option of crystallizing the solid directly from solution at all concentrations. Two key questions that remain unanswered thus far are (a) how a solution of a salt such as NaCl, which does not form hydrated crystals, loses thermodynamic stability at sufficiently high concentrations and (b) the implications of the loss of stability on the nucleation mechanism. In the current work, we provide unambiguous answers to these questions, by careful simulations of two model systems. Specifically, we study the nucleation in a model Lennard-Jones (LJ) mixture and of NaCl from supersaturated aqueous solutions. We obtain component chemical potentials that accurately represent the driving force for nucleation, and we identify the stability limit of the supersaturated NaCl solution, for the first time in simulation studies. Prior to reaching the spinodal, there are no liquid pre-nucleation clusters and crystal nucleation follows a single-step mechanism. At and beyond the thermodynamic stability limit, we observe a shift from one-step to a two-step nucleation mechanism, i.e., a liquid/amorphous phase separation producing clusters, followed by crystallization, rather than a barrier-free crystal/solution spinodal decomposition. We also calculate crystal nucleation rates at several supersaturations.

We use the SPC/E (extended simple point charge)32 model of water and the Joung-Cheatham33 (JC) NaCl force fields to model supersaturated aqueous NaCl solutions. The SPC/E and JC force fields have been shown to provide reasonable predictions for several solution thermodynamic and transport properties.34 The JC NaCl model has an equilibrium solubility of 3.7 mol/kg in SPC/E water at 298.15 K and 1 bar (experimental value 6.1 mol/kg), confirmed by both chemical potential35 and direct coexistence methods.36 The nucleation rate for this system has been recently studied12,37 using forward flux sampling (FFS) methods and from seeding simulations in conjunction with classical nucleation theory (CNT);19,38 it was demonstrated that at modest supersaturations, nucleation follows a classical one-step mechanism.

The difference in chemical potentials of ions between solution and crystal phases is the driving force for nucleation. We calculate the chemical potential of NaCl (ions) in solution from 6.0 to 20.0 mol/kg, following the approach developed by Mester and Panagiotopoulos.35,39 In particular, the chemical potential is estimated in molecular dynamics (MD) simulations in the isothermal-isobaric ensemble as the change in the Gibbs free energy due to the insertion of a pair of ions into the solution. The Gibbs free energy of insertion is obtained through a thermodynamic integration process by slowly switching on the interactions between the inserted pair of ions and the solution. Further details on the calculation of electrolyte chemical potentials can be found in Refs. 35 and 39. For all MD simulations associated with the chemical potential calculation, the simulation box has 500 water molecules, and the number of ion pairs ranges from 55 to 180, corresponding to different salt concentrations. Because a relatively small number of ions are used in these calculations of the chemical potential, we do not observe significant crystal nucleus formation during the production stage (around 600 ns) of our simulations. More details on the calculation of electrolyte chemical potentials are provided in the supplementary material.

In order to follow the progress of nucleation from aqueous solution, it is necessary to distinguish ions that are in the solution and in the crystalline phases. For each ion (i), the Steinhardt bond-orientational order parameter40q8 is calculated as

q8(i)=l=88|q8l(i)|2
(1)

and

q8l(i)=1NBj=0NBY8lθrij,ϕrij,
(2)

where Y8l are the spherical harmonics and θ(rij) and ϕ(rij) are the polar and azimuthal angles associated with the vector (rij) that connects the central ion (i) and one of its neighbor ions (j). The summation in Eq. (2) is over the 12 nearest neighbors of the ion i (NB = 12). We consider an ion to be in the crystalline phase if its q8 order parameter is larger than 0.45 and two crystalline ions that are separated by a distance less than 0.35 nm are considered in the same crystalline nucleus. A similar strategy has been used by Lanaro and Patey13 to follow the formation of NaCl nuclei in aqueous solutions.

When the free energy barrier for nucleation is low (of the order of a few kBT), spontaneous nucleation can be observed from unbiased MD simulations within hundreds of nanoseconds. For such spontaneous nucleation events, the nucleation free energy profile can be extracted from the mean first passage time (MFPT) as41,42

βG*(n*)=ln[B(n)]an*dxB(x)+C
(3)

and

B(x)=1Pst(x)axPst(x)dxτ(x)τ(b),
(4)

where τ(n*) is the MFPT collected from MD simulations as the average time required for the largest nucleus in a trajectory to reach a size of n* for the first time. Pst(n*) is the steady state probability that a configuration has a largest crystalline nucleus of size n*. a is the boundary of the solution domain, which is taken as the position where Pst shows a maximum for a solution before experiencing any crystallization. b is the boundary of the crystal domain, chosen to be 45 in this work, since a critical crystalline nucleus is generally much smaller than 45 for highly supersaturated solutions.

The above free energy profile from MFPT is expressed as a function of the size of the largest crystalline nucleus. However, the free energy barrier associated with a crystallization process corresponds to the reversible work needed to assemble a crystalline nucleus of size n, and it is related to the probability that a system has crystalline nuclei of size n (N(n)), rather than the probability of observing the largest nucleus in a system to have size n (N*(n)).43 Thus, the free energy profile of nucleation is expressed as

βG(n)=lnN(n)N(0),
(5)

where N(0) refers to the average number of ions in the solution phase. The above defined G(n) and G*(n*) profiles differ for small n (or n*) and are practically the same for large n (or n*).44 Therefore, following the approach of Lundrigan and Saika-Voivod,44 we first obtain the nucleus size distribution (N(n)) from a short MD simulation (<10 ns) and calculate G(n) using Eq. (5) up to a small n (<10) as crystallization has not occurred. G*(n) is then obtained from Eq. (3) and patched onto G(n) for large values of n. The nucleation rate J is determined, as a fitting parameter, from the MFPT (τ) using the following expression:44 

τ(n*)=12JV1+erfc(n*nc),
(6)

with nc (the size of a critical crystalline nucleus) and c also considered as fitting parameters. V is the system volume. Strictly speaking, Eq. (6) becomes less valid when the supersaturation is sufficiently high. However, a satisfactory fitting of τ (see Fig. S1 of the supplementary material for details) suggests that Eq. (6) is still valid for the supersaturations considered in this work.

We also obtained the free energy profile of nucleation using umbrella sampling from hybrid-Monte Carlo (MC) simulations.45 For the umbrella sampling simulations, the size of the largest crystalline nucleus is chosen as the reaction coordinate in the bias potential, expressed as ϕ(n,n0)=k0(nn0)2/2. Following Frenkel and Saika-Voivod et al.,46,47 we perform parallel simulations at different umbrella windows (n0) in order to drive the nucleation. At given n0, the probability that the largest nucleus in the system is of size n, defined as Pmax(n), is calculated, and a histogram of nucleus size distribution in the unbiased ensemble (N(n)) is computed from the distribution in the biased ensemble (N¯(n)) as N(n)=exp[βϕ(n,n0)]N¯(n)biased. We then trim N(n) by discarding histogram entries for which Pmax(n) is less than 0.05 to ensure good sampling statistics. With N(n) thus obtained, the nucleation free energy G(n) is determined in each umbrella sampling window up to an additive constant, and the free energy profile is then constructed by minimizing the difference between overlapping portions of G(n) in each window, as done in Refs. 46 and 47. More details for the umbrella sampling in hybrid-MC simulations are given in the supplementary material.

The chemical potential of NaCl in solution, together with the chemical potential of the crystal, defines the driving force for nucleation. The chemical potential of NaCl in solution up to a salt concentration of 20.0 mol/kg is given in Fig. 1, and the equilibrium solubility of the rock-salt crystal is also shown in the figure as a vertical dashed line. As shown in Fig. 1, the electrolyte chemical potential as a function of salt concentration reaches a maximum (−366.6 ± 0.3 kJ/mol) around 15.0 mol/kg and plateaus beyond this molality, within simulation uncertainty, up to at least 20 mol/kg (no calculations were performed beyond this salt concentration). Within numerical accuracy of our simulations, the derivative of chemical potential with respect to concentration is zero near 15.0 mol/kg, indicating that the system reaches a thermodynamic stability limit (spinodal) approximately at 15.0 mol/kg. At the spinodal, a new phase must necessarily emerge from the aqueous solution. It is worth mentioning that at and beyond the spinodal with the new phase coexisting with the aqueous solution phase, the chemical potential reported in Fig. 1, as the Gibbs free energy associated with the insertion of a pair of ions, refers to the average chemical potential of NaCl in the system. The constancy of chemical potentials between 15.0 and 20.0 mol/kg suggests that the system experiences a barrier-less phase separation as the mother phase (aqueous solution) is inherently unstable beyond the spinodal. However, it has been shown in the binary LJ fluid mixtures that the location of the spinodal identified as the composition where the chemical potential reaches its maximum depends on the simulation system size.48–50 In fact, the spinodal, although a useful concept, is really a mean-field approximation. In practice, the width of the region across which the system evolves from metastability to instability is finite.51 For finite fluid systems such as the binary LJ mixtures in Refs. 48–50 and the aqueous NaCl solutions studied here, the phase separation process evolves sharply (at a size-dependent spinodal) from nucleation/growth to spinodal decomposition as the free energy barrier of phase separation decreases with increasing supersaturation, and, for a finite system, vanishes at the spinodal. Because we focus primarily on the change of nucleation mechanism with solution supersaturation, and due to the very different types of interactions between the ionic aqueous solutions and short-ranged LJ mixtures, we do not investigate further the precise location of the spinodal and its possible system size dependence. Similar to the maximum chemical potential shown in Fig. 1, Moučka et al. identified in their osmotic ensemble Monte Carlo (OEMC) simulations a limiting chemical potential in aqueous NaCl52,53 and several other alkali-halide salts solutions52 when the salt concentration is sufficiently high and the solubility predicted by the underlying force field is low. The existence of a limiting chemical potential was found to be a generic feature across different force fields.54 Beyond the salt concentration that corresponds to the limiting chemical potential, a precipitate of ions that has rock-salt (FCC) crystalline structure was observed in the supersaturated NaCl solution (see Fig. 5 of Ref. 52). While the observed precipitation phenomenon was connected to the homogeneous crystallization of ions, the detailed mechanism associated with such precipitation remains unknown.

FIG. 1.

Electrolyte chemical potential as a function of salt concentration at 298 K and 1 bar. The curved dashed line is a guide to the eye, and the vertical dashed line indicates the equilibrium solubility for the SPC/E + JC force field combination. The uncertainty of the simulation data is comparable to or smaller than the symbol size.

FIG. 1.

Electrolyte chemical potential as a function of salt concentration at 298 K and 1 bar. The curved dashed line is a guide to the eye, and the vertical dashed line indicates the equilibrium solubility for the SPC/E + JC force field combination. The uncertainty of the simulation data is comparable to or smaller than the symbol size.

Close modal

Since the aqueous NaCl solution reaches its stability limit approximately at 15.0 mol/kg, one would normally expect that crystallization of ions from the unstable aqueous solution at and beyond the spinodal is a barrier-free spinodal decomposition process. In order to understand the nucleation mechanism close to the spinodal, we calculated the free energy profile (G(n)) for crystallization of ions at 15.0 mol/kg using both the MFPT from MD simulations and umbrella sampling from hybrid-MC simulations, as described above. As shown in Fig. 2, the free energy profile (G(n)) at 15.0 mol/kg has a barrier of approximately 12 kBT, and both the MFPT and umbrella sampling methods give essentially identical results. From umbrella sampling simulations with systems that have 500 and 1000 ions, we obtain the same free energy barriers within simulation uncertainties, indicating that the system size effect on the free energy barrier of nucleation is negligible. If crystallization proceeds by spinodal decomposition, the free energy barrier would vanish. However, this is clearly not the case at the spinodal of the aqueous NaCl solution.

FIG. 2.

Free energy of formation crystalline nuclei of size n, G(n), at 15 mol/kg, 298 K, and 1 bar. The uncertainty of the simulation data is comparable to or smaller than the symbol size.

FIG. 2.

Free energy of formation crystalline nuclei of size n, G(n), at 15 mol/kg, 298 K, and 1 bar. The uncertainty of the simulation data is comparable to or smaller than the symbol size.

Close modal

Given that the system is not experiencing spinodal decomposition toward crystallization, it is interesting to investigate if classical nucleation theory (CNT), which is only valid at low supersaturations, is still applicable at the spinodal for the solution. From CNT, the nucleation rate is calculated as

JCNT=ρf+Zexp(ΔGkBT),
(7)

where ρ is the number density of NaCl (around 6.5 × 1027 m−3 at 15.0 mol/kg). f+ is the attachment rate of ions to a critical nucleus, which can be estimated as the slope of the mean squared evolution of the largest nucleus size, n*, with time (t), as [n*(t)n*(t=0)]2/2t. n*(t = 0) is close to the critical nucleus size. 100 independent MD simulations are initialized from systems that have largest crystalline nuclei around 20 particles (close to the size of critical crystalline nucleus at 15.0 mol/kg), and the attachment rate is estimated as around 4.6 ns−1 (see the supplementary material for more details). The Zeldovich factor (Z) is associated with the curvature of the free energy profile (G(n)) at the barrier and is given by Z=G(n)/(2πkBT). Using the G(n) obtained from umbrella sampling, the Zeldovich factor is calculated as 0.11 ± 0.03. Considering a nucleation barrier of 12 kBT, the nucleation rate J is calculated from Eq. (7) as 2 ± 1 × 1031 m−3 s−1. Without resorting to CNT, the nucleation rate at 15.0 mol/kg is obtained directly from MD simulations using Eq. (6) from MFPT as 4 ± 1 × 1031 m−3 s−1. The rate from CNT agrees very well with the result from the MFPT, within numerical uncertainties. This indicates that CNT is still valid for calculation of the nucleation rates at the spinodal for the prediction of crystal nucleation rate and again confirms that crystallization at the spinodal has a relatively large (12 kBT) barrier. It is important to note that while CNT may be valid to calculate nucleation rates, there is no evidence at this point to suggest that the crystallization of ions at and beyond spinodal is necessarily a CNT-like single-step process. For a multi-step nucleation mechanism, the nucleation rate may be well predicted by CNT, provided that the rate is determined only by the one step that has a significant higher free energy barrier than other steps.

We extended the analysis to a broader range of salt concentrations. Figure 3 shows the nucleation free energy barriers estimated either from MFPT or umbrella sampling, size of critical nuclei estimated from the position of maximum in G(n), and nucleation rates from fitting MFPT with Eq. (6), at different salt concentrations. In addition, the self-diffusion coefficients of the ions, estimated from the mean squared displacement using the Einstein relation, are shown in Fig. 3. As expected, the free energy barrier for nucleation of ions into the crystal decreases as the salt concentration increases. The size of critical nucleus (nc) decreases with salt concentration and seems to converge to a constant (around 20); however, the estimation of nc is subject to large uncertainty for systems close to the spinodal due to the relatively flat free energy profiles (shown in the supplementary material) near the barrier. It is noted that the number of ions in critical nuclei is much smaller than the total number of ions in the system, which suggests that the solution supersaturation does not change significantly due to nucleation. The nucleation rates (J) increase with salt concentration, while the diffusion coefficients of the ions decrease with concentration. If the crystallization of ions took place via spinodal decomposition, which is diffusion controlled, the nucleation rates would decrease with concentration since the diffusion of ions becomes slower as concentration increases. Therefore, the increasing nucleation rates in combination with decreasing diffusion coefficients rules out the possibility of spinodal decomposition-type crystallization. Also, as the system approaches the spinodal, the free energy barriers and nucleation rates vary smoothly with increasing concentration, without any detectable sudden change in their slopes.

FIG. 3.

(a) Nucleation free energy barriers (G(nc)) at different salt concentrations obtained from MFPT and umbrella sampling (at 13.2 and 15.0 mol/kg). (b) Size of the critical nucleus (nc) at different salt concentrations. Data at 10.0 and 12.0 mol/kg are analyzed from data collected during FFS calculation in Ref. 12. (c) Nucleation rates (J) at different salt concentrations obtained from MFPT using Eq. (6). (d) Self-diffusion coefficient of ions (DselfNa+ and DselfCl) at different salt concentrations. The dashed lines are guides to the eye.

FIG. 3.

(a) Nucleation free energy barriers (G(nc)) at different salt concentrations obtained from MFPT and umbrella sampling (at 13.2 and 15.0 mol/kg). (b) Size of the critical nucleus (nc) at different salt concentrations. Data at 10.0 and 12.0 mol/kg are analyzed from data collected during FFS calculation in Ref. 12. (c) Nucleation rates (J) at different salt concentrations obtained from MFPT using Eq. (6). (d) Self-diffusion coefficient of ions (DselfNa+ and DselfCl) at different salt concentrations. The dashed lines are guides to the eye.

Close modal

Figure 4 shows two snapshots taken at different times from an MD simulation of nucleation at the spinodal (15.0 mol/kg). The system has 7500 ions with a box length around 10.5 nm. Only ions that are labeled crystalline [q8 > 0.45, see Eq. (1)] are shown in the snapshots for clarity. Here, we chose to use a system that is much larger than those used in free energy and rates calculations in order to enable a more clear observation of interconnected domains (if they exist), which is a characteristic feature of spinodal decomposition. The snapshots show that at the early stages of nucleation (t = 50 ns), only one small crystalline nucleus emerges from the solution and as the system evolves, more and larger crystalline nuclei are present (e.g., at t = 150 ns); however, the nuclei are not interconnected. The nuclei show clear rock-salt (FCC) structure with an approximately spherical boundary, which is more consistent with CNT rather than with the Cahn-Hilliard mean-field theory for spinodal decomposition,55 in which the nuclei are ramified with diffuse boundaries. Again, for the crystallization of ions at the spinodal, the system shows characteristics of nucleation/growth instead of spinodal decomposition. Additionally, 3 movies showing the progression of nucleation events at concentrations at 13.8 (before the spinodal), 15.0 (at the spinodal), and 18.5 mol/kg (beyond the spinodal), respectively, are included in the supplementary material. Spinodal decomposition-type crystallization is not observed in these representative movies.

FIG. 4.

Snapshots of crystalline nuclei at 15.0 mol/kg. (a) Snapshot at 50 ns and (b) snapshot at 150 ns. Only crystalline ions with q8 > 0.45 [see Eqs. (1) and (2)] are shown in the snapshots. Blue particles: Na+. Cyan particles: Cl.

FIG. 4.

Snapshots of crystalline nuclei at 15.0 mol/kg. (a) Snapshot at 50 ns and (b) snapshot at 150 ns. Only crystalline ions with q8 > 0.45 [see Eqs. (1) and (2)] are shown in the snapshots. Blue particles: Na+. Cyan particles: Cl.

Close modal

The above analysis demonstrates that nucleation/crystallization of ions at the limit of stability does not proceed by spinodal decomposition into a crystal phase. However, the presence of a spinodal indicated by the composition dependence of the chemical potential (Fig. 1) certainly suggests that spinodal decomposition occurs in the solution, and a new phase must necessarily emerge from the solution, at and beyond 15 mol/kg. In order to understand this apparent contradiction, we first investigate a binary Lennard-Jones (LJ) system with a LJ solute (S) supersaturated in an LJ solvent (SV), in analogy with the supersaturated NaCl solution.

The interaction parameters of the binary LJ system were taken from Ref. 36 (see Table I in Ref. 36), and the mass of argon is used for both solute and solvent LJ particles. The simulation of the binary LJ mixture was conducted at 50 K and 1 bar, corresponding to a reduced temperature (kBT/ε) of 0.44 for the solute (S) and 0.81 for the solvent (SV), respectively, making the binary system exhibit a crystal/solution phase separation at equilibrium, similar to the NaCl/water case. The stable crystal phase in the binary LJ mixture consists of pure solute (S) particles because of the disparity in particle size. The equilibrium solute solubility is (xS =) 0.036(2), as reported in Ref. 36.

Following the approach for the calculation of electrolyte chemical potentials, the chemical potential of the solute (S) in the solvent (SV) was obtained from MD simulations at different solute concentrations. The details of these simulations are given in the supplementary material. As shown in Fig. 5, similar to the supersaturated aqueous NaCl solution, the binary LJ solution has a spinodal at xS (mole fraction of solute) approximately to 0.2, indicated by a zero derivative of chemical potential with respect to the solute concentration. The chemical potential remains almost a constant from the spinodal to at least xS = 0.4, suggesting that the system undergoes phase separation by spinodal decomposition.

FIG. 5.

(a) Chemical potential of the solute LJ component (S) at different mole fractions. Ideal gas contribution to the chemical potential of solute is set to 0 for a system with a solute number density of 10−27 m3. (b) Free energy of formation of crystalline solute (S) nuclei of size n at a mole fraction of 0.4. The uncertainty of the simulation data is comparable to or smaller than the symbol size. (c) Snapshot for the binary LJ system with solute mole fraction (xS) at 0.24, near the spinodal. (d) Snapshot for the binary LJ system with solute mole fraction (xS) at 0.4, beyond the spinodal. The binary mixture has a total of 1000 LJ particles, and both snapshots are taken from NPT-MD simulations at t = 10 ns. Blue particles: solute (S). Cyan particles: solvent (SV).

FIG. 5.

(a) Chemical potential of the solute LJ component (S) at different mole fractions. Ideal gas contribution to the chemical potential of solute is set to 0 for a system with a solute number density of 10−27 m3. (b) Free energy of formation of crystalline solute (S) nuclei of size n at a mole fraction of 0.4. The uncertainty of the simulation data is comparable to or smaller than the symbol size. (c) Snapshot for the binary LJ system with solute mole fraction (xS) at 0.24, near the spinodal. (d) Snapshot for the binary LJ system with solute mole fraction (xS) at 0.4, beyond the spinodal. The binary mixture has a total of 1000 LJ particles, and both snapshots are taken from NPT-MD simulations at t = 10 ns. Blue particles: solute (S). Cyan particles: solvent (SV).

Close modal

The free energy barrier for crystallization of the solute particle at xS = 0.4, which is beyond the spinodal, was obtained using umbrella sampling from hybrid-MC simulations, as done for the NaCl solution. Details of the hybrid-MC simulations for the binary LJ systems are provided in the supplementary material. We use the procedure proposed by Jungblut and Dellago56 to identify the crystalline LJ nucleus. In particular, two solute (S) particles are considered as neighbors if they are within 1.5 σ and the normalized scalar product of their complex q6l40 vectors is larger than 0.5. A solute particle that has more than 8 neighbors is considered as a crystalline particle, and two crystalline particles that are within 1.5 σ of each other belong to the same crystalline nucleus. The free energy profile for crystallization at xS = 0.4 is given in Fig. 5. Again, similar to the NaCl solution, the free energy profile of the LJ mixture shows a non-vanishing barrier around 14 kBT at xS = 0.4, which is beyond the spinodal. Thus, the crystallization of the solute in the LJ solvent does not proceed by spinodal decomposition, again similar to the NaCl solution.

However, as shown in the snapshots for systems at two different solute concentrations beyond the spinodal, the solute and solvent LJ particles experience liquid/amorphous phase separation, in which the solute particles aggregate into a metastable disordered phase without crystalline structure. The presence of the solute phase suggests that the spinodal corresponds to a liquid/amorphous phase separation instead of the solution/crystal spinodal decomposition, and the free energy barrier shown in Fig. 5 corresponds to the crystallization of the solute in the metastable solute amorphous phase. A similar picture for binary LJ solutions has been previously observed by Anwar and Boateng.5 

Based on the behavior of the binary LJ solution, we expect the supersaturated NaCl solution to also exhibit a liquid/amorphous phase separation at and beyond the spinodal. Here, we choose to analyze an MD trajectory for the simulation of a NaCl solution at a salt concentration of 18.5 mol/kg, which is beyond the spinodal. Ions that have more than Ncut ion neighbors within a radius of rcut are considered to be in an amorphous phase regardless of the local structure of their neighboring ions, and two amorphous (not solution) ions that are within rcut are considered in the same cluster. We set rcut as 0.45 nm, corresponding roughly to the position of the second peak of the cation/anion pair correlation function. Ncut is set to 8, as in the binary LJ mixture. Figure 6 shows the four largest clusters identified using such criteria, as well as their hydration water molecules surrounding these clusters. A water molecule is considered as hydration water if it is within 0.4 nm of any ion in a cluster. The snapshots are from a system with 9000 ion pairs and 27 000 water molecules, taken at 10 ns (a) and 70 ns (b) of a MD simulation. The two largest crystalline nuclei are also shown in the snapshots.

FIG. 6.

Snapshots from a MD simulation of aqueous NaCl solution at a concentration of 18.5 mol/kg. (a) Snapshot at t = 10 ns and (b) snapshot at t = 70 ns. Only ions with more than 8 ion neighbors (and their hydration water) are shown in the snapshots. Blue and cyan particles: “amorphous salt” ions; red and orange particles: crystalline ions. Small red and white atoms represent the oxygen and hydrogen atoms of water.

FIG. 6.

Snapshots from a MD simulation of aqueous NaCl solution at a concentration of 18.5 mol/kg. (a) Snapshot at t = 10 ns and (b) snapshot at t = 70 ns. Only ions with more than 8 ion neighbors (and their hydration water) are shown in the snapshots. Blue and cyan particles: “amorphous salt” ions; red and orange particles: crystalline ions. Small red and white atoms represent the oxygen and hydrogen atoms of water.

Close modal

As shown in Fig. 6, the clusters are disordered with dimensions on the order of few nanometers. While the above criteria to identify amorphous ions clusters are ad-hoc, we note that water molecules only surround, but do not penetrate, these relatively large clusters, indicating that the amorphous clusters correspond to a more dense phase compared to the aqueous solution. At t = 10 ns, which corresponds to the early stages of crystallization, ions form the disordered liquid amorphous clusters, with only a small number of ions (a total of less than 20) forming crystalline nuclei (not shown in the snapshots). Unlike a crystalline nucleus (see Fig. 4), the amorphous clusters are ramified and have diffuse boundaries, which is consistent with the Cahn-Hilliard mean field theory.55 As the system evolves, at t = 70 ns, part of the ions inside the amorphous salt clusters have already rearranged into crystalline nuclei, shown as red and orange particles. Similar to the disordered pre-nucleation NaCl clusters here, recent optical Kerr-effect spectroscopy experiments confirmed the presence of large metastable ion clusters in supersaturated sodium thiosulfate solutions.57 

Based on the above arguments, the spinodal observed from the electrolyte chemical potential vs. concentration curve (Fig. 1) identifies the stability limit of the aqueous NaCl solution, with respect to concentration fluctuations, not crystallization, and corresponds to the spontaneous formation of amorphous salt clusters, i.e., a solution/amorphous salt phase separation, rather than the solution/crystal spinodal decomposition. Beyond the spinodal, the system loses its thermodynamic stability as a metastable solution and starts to form a microscopic amorphous salt phase. Since the amorphous salt clusters have not lost their thermodynamic stability with respect to crystal, the free energy barrier for crystallization is non-vanishing. The NaCl crystal nucleation is thus a two-step process at concentrations beyond the spinodal: a solution/amorphous (liquid/amorphous) phase separation with a vanishing (or low) free energy barrier followed by rearrangement of amorphous ions into crystalline nuclei (liquid/solid phase separation) with a non-vanishing free energy barrier.

Despite the fact that different crystallization processes involve liquid/liquid (or liquid/amorphous) phase separation, the effects of liquid/liquid phase separation on crystallization processes could well be quite varied. In protein and colloid solutions, the bulk liquid/liquid phase separation that occurs near the fluid/fluid critical point induces large density fluctuations, and the free-energy barrier for crystal nucleation is strongly reduced.6 In polymer blends, it has been observed that crystallization may be enhanced by liquid/liquid phase separation as such phase separation generates phase boundaries or precursors, which subsequently induce heterogeneous crystal nucleation.58,59 In contrast to the protein/colloid solutions and polymer blends, the crystallization of NaCl ions evolves smoothly from single-step to two-step as the nucleation free energy barrier and rates (see Fig. 3) do not change abruptly as the system crosses the spinodal.

As discussed above, the nucleation of ions from solution into a crystal at and beyond the spinodal is a two-step process: ions first form disordered amorphous salt clusters, and ions in the amorphous clusters rearrange into crystalline nuclei with a free energy barrier. The two-step (or non-classical) nucleation/crystallization mechanism has been observed, often at large supersaturations, in several systems for crystallization from solution,26,60,61 in which the density/concentration fluctuation first leads to the formation of dense phases/precursors, and then solute molecules in the dense phases/precursors transform into an ordered crystal phase with structure fluctuations in a second step. On the other hand, for a single-step nucleation mechanism, concentration and structure fluctuations occur simultaneously, leading directly to the formation of crystalline nuclei. Several experiments and simulation studies have shown that the single-step nucleation pathway is less favored than the indirect, two-step mechanism for nucleation/crystallization from solutions.62 However, there are interesting scenarios where two- and single-step crystallization mechanisms coexist,63 or a system prefers single-/two-step pathway under different conditions. For example, despite many experimental and simulation studies that show a two-(or multi-) step mechanism, crystallization of CaCO3 on organothiol self-assembled monolayer follows a direct one-step pathway without the presence of amorphous precursors prior to nucleation.64 Glucose isomerase, which is known to exhibit characteristics of multi-step nucleation at a concentration of 100 mg/ml, crystallizes in solution following a classical one-step mechanism at a concentration of 0.1 mg/ml.62 The nucleation of urea follows a two-step mechanism in water and acetonitrile solutions, while a single-step pathway is more favored in methanol and ethanol.22 

Since the nucleation of NaCl crystals follows a two-step process at and beyond the spinodal, it would be interesting to investigate if the two-step mechanism also holds before the system reaches its spinodal. Here, we analyze trajectories collected from unbiased MD simulations that include nucleation events, for systems at different salt concentrations. Figure 7 shows the distributions of sizes of largest amorphous salt clusters (Ncluster) versus the sizes of largest crystalline nuclei (Ncrystal) for NaCl solutions before reaching the spinodal (at concentrations of 10.0 and 12.0 mol/kg) and after the spinodal (at 16.0 mol/kg). Each symbol in Fig. 7 corresponds to a configuration collected during MD simulations. The data at 10.0 mol/kg and 12.0 mol/kg are based on configurations collected at milestones of forward flux sampling (FFS) calculations of nucleation rates in Ref. 12. The data at 16.0 mol/kg are obtained using configurations sampled from an unbiased MD simulation (with 1000 ions pairs) of 120 ns, which produces a crystalline nucleus of 45 ions. Crystalline nuclei that are very small (less than 3-4 ions) and do not reside in the largest disordered clusters are excluded from the analysis. It is noted that the configurations collected at different FFS milestones are joined by stochastic events due to the randomization of particle velocities required by the FFS algorithm, while the configurations collected for system at 16.0 mol/kg belong to a deterministic MD trajectory. We do not expect such a difference in sampling methods to affect Ncluster and Ncrystal as these quantities convey structural information. The distribution of sizes of the largest amorphous salt clusters (Ncluster) quantifies concentration fluctuations, while the distribution of sizes of the largest crystalline nuclei (Ncrystal) measures structural fluctuations. As shown in Fig. 7, for solutions before reaching the spinodal (10.0 and 12.0 mol/kg), Ncluster correlates linearly with Ncrystal, which indicates a simultaneous fluctuation of concentration and structure, consistent with a single-step nucleation mechanism. In fact, at low supersaturations away from the spinodal, the largest amorphous salt cluster, identified using the criteria discussed above, is essentially the ordered largest crystalline nucleus with minor differences associated with the identification of ions at the interface of the cluster/nucleus. As nucleation proceeds, the ions nucleate into a crystalline nucleus simultaneously as they aggregate. After the system passes the spinodal (at 16.0 mol/kg), Ncluster has a much wider distribution at early stages of nucleation (small Ncrystal), and the concentration fluctuations are not coupled with structure fluctuations, which indicates a two-step nucleation mechanism.

FIG. 7.

Distribution for sizes of the largest ion clusters (Ncluster) and sizes of the largest crystalline nuclei (Ncrystal) at salt concentrations of 10.0, 12.0, and 16.0 mol/kg. Solid lines indicate the evolution of Ncluster and Ncrystal with time.

FIG. 7.

Distribution for sizes of the largest ion clusters (Ncluster) and sizes of the largest crystalline nuclei (Ncrystal) at salt concentrations of 10.0, 12.0, and 16.0 mol/kg. Solid lines indicate the evolution of Ncluster and Ncrystal with time.

Close modal

The above analysis has shown that there is a transition from the single to two-step crystallization mechanism in the supersaturated NaCl solution, and such transition is driven by the microscopic solution/amorphous salt spinodal. It is possible that the two-step nucleation is a mechanism that is only associated with systems that are close to, or beyond, their stability limits, and particles nucleate via the single-step mechanism for systems well before reaching the spinodal. The general belief that classical nucleation theory is only strictly applicable at low supersaturations may be interpreted more precisely as suggesting that the theory is only strictly valid at supersaturations well before reaching the solution spinodal. Such hypothesis may possibly explain the preference of single/two-step crystallization pathway under different solution conditions. Regardless of the validity of the proposed hypothesis, the spinodal driven transition of crystallization mechanisms suggests that an accurate determination of crystallization driving forces from the calculation of solution/crystal chemical potentials is of special importance in order to interpret the nucleation mechanism from simulations.

In this work, we have examined the nucleation of ions from aqueous NaCl solutions using molecular simulations. For the SPC/E water and JC ions force fields, the chemical potential of the electrolyte starts to plateau at salt concentrations around 15.0 mol/kg, suggesting the presence of a spinodal. The spinodal indicates that the aqueous NaCl solution reaches its limit of stability with respect to concentration fluctuations, and a phase separation through spinodal decomposition occurs in the solution. However, the nucleation of ions into a rock-salt crystal at the spinodal shows no signatures of solution-crystal spinodal decomposition: the free energy barrier to crystallization does not vanish, the critical crystalline nuclei remains finite and have an approximately spherical shape, and the nucleation rates increase with ion concentration, in spite of the corresponding slowing down of ion diffusive motion. Therefore, the spinodal in the system is not associated with crystal/solution phase separation. It is found that the solution at the spinodal experiences a microscopic liquid/amorphous phase separation similar to that previously observed in a binary Lennard-Jones solution.5 The ions form a metastable amorphous salt phase in transient coexistence with the solution, and the nucleation of ions occurs in the amorphous salt clusters with a free energy barrier. While the nucleation is a two-step process beyond the spinodal, a single-step nucleation to the crystalline phase is observed for solutions before reaching the spinodal. The transition from one-step to two-step crystallization mechanism, driven by the solution spinodal, suggests that an accurate knowledge of solution and crystal chemical potentials is of importance in order to obtain a proper interpretation of mechanisms from simulations of nucleation. We hypothesize that the two-step mechanism observed in NaCl aqueous solutions and in other systems with different solute/solvent interactions (e.g., binary LJ system) may be a general feature for systems in the vicinity of a stability limit.

Note added in proof. After the original submission of this work, a highly relevant study by Poole and co-workers65 has appeared on the role of fluctuations and metastability in two-step nucleation.

See supplementary material for details of molecular dynamics and hybrid Monte Carlo simulations, mean first passage time for estimation of nucleation rates, crystallization free energy profile at salt concentrations of 13.2, 13.8, 14.6, and 15.0 mol/kg, attachment rates of critical nuclei at 15.0 mol/kg, and movies showing progression of nucleation events at 13.8, 15.0, and 18.5 mol/kg.

Financial support for this work was provided by the Office of Basic Energy Sciences, U.S. Department of Energy, under Award No. DE-SC0002128. Additional support was provided by the National Oceanic and Atmospheric Administration (Cooperative Institute for Climate Science Award No. AWD 1004131) and the U.S. National Science Foundation under Award No. CBET-1402166.

1.
D.
Erdemir
,
A. Y.
Lee
, and
A. S.
Myerson
, “
Nucleation of crystals from solution: Classical and two-step models
,”
Acc. Chem. Res.
42
,
621
629
(
2009
).
2.
G. C.
Sosso
,
J.
Chen
,
S. J.
Cox
,
M.
Fitzner
,
P.
Pedevilla
,
A.
Zen
, and
A.
Michaelides
, “
Crystal nucleation in liquids: Open questions and future challenges in molecular dynamics simulations
,”
Chem. Rev.
116
,
7078
7116
(
2016
).
3.
A.
Peckhaus
,
A.
Kiselev
,
R.
Wagner
,
D.
Duft
, and
T.
Leisner
, “
Temperature-dependent formation of NaCl dihydrate in levitated NaCl and sea salt aerosol particles
,”
J. Phys. Chem. B
145
,
244503
(
2016
).
4.
A.
van Driessche
,
M.
Kellermeier
,
L. G.
Benning
, and
D.
Gebauer
,
New Perspectives on Mineral Nucleation and Growth
(
Springer International Publishing
,
2017
).
5.
J.
Anwar
and
P. K.
Boateng
, “
Computer simulation of crystallization from solution
,”
J. Am. Chem. Soc.
120
,
9600
9604
(
1998
).
6.
P.
ten Wolde
and
D.
Frenkel
, “
Enhancement of protein crystal nucleation by critical density fluctuations
,”
Science
277
,
1975
1978
(
1997
).
7.
A.
Haji-Akbari
and
P.
Debenedetti
, “
Direct calculation of ice homogeneous nucleation rate for a molecular model of water
,”
Proc. Natl. Acad. Sci. U. S. A.
112
,
10582
10588
(
2015
).
8.
G. C.
Sosso
,
T.
Li
,
D.
Donadio
,
G. A.
Tribello
, and
A.
Michaelides
, “
Microscopic mechanism and kinetics of ice formation at complex interfaces: Zooming in on kaolinite
,”
J. Phys. Chem. Lett.
7
,
2350
2355
(
2016
).
9.
G. D.
Soria
,
J. R.
Espinosa
,
J.
Ramirez
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
, “
A simulation study of homogeneous ice nucleation in supercooled salty water
,”
J. Chem. Phys.
148
,
222811
(
2018
).
10.
K.
Mochizuki
,
Y.
Qiu
, and
V.
Molinero
, “
Promotion of homogeneous ice nucleation by soluble molecules
,”
J. Am. Chem. Soc.
139
,
17003
17006
(
2017
).
11.
M.
Salvalaglio
,
C.
Perego
,
F.
Giberti
,
M.
Mazzotti
, and
M.
Parrinello
, “
Molecular-dynamics simulations of urea nucleation from aqueous solution
,”
Proc. Natl. Acad. Sci. U. S. A.
112
,
E6
E14
(
2015
).
12.
H.
Jiang
,
A.
Haji-Akbari
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
Forward flux sampling calculation of homogeneous nucleation rates from aqueous NaCl solutions
,”
J. Chem. Phys.
148
,
044505
(
2018
).
13.
G.
Lanaro
and
G. N.
Patey
, “
Birth of NaCl crystals: Insights from molecular simulations
,”
J. Phys. Chem. B
120
,
9076
9087
(
2016
).
14.
M.
Mucha
and
P.
Jungwirth
, “
Salt crystallization from an evaporating aqueous solution by molecular dynamics simulations
,”
J. Phys. Chem. B
107
,
8271
8274
(
2003
).
15.
G.
Lanaro
and
G. N.
Patey
, “
The influence of ion hydration on nucleation and growth of LiF crystals in aqueous solution
,”
J. Chem. Phys.
148
,
024507
(
2018
).
16.
G. A.
Tribello
,
F.
Bruneval
,
C.
Liew
, and
M.
Parrinello
, “
A molecular dynamics study of the early stages of calcium carbonate growth
,”
J. Phys. Chem. B
113
,
11680
11687
(
2009
).
17.
A.
Wallace
,
L.
Hedges
,
A.
Fernandez-Martinez
,
P.
Raiteri
,
J.
Gale
,
G.
Waychunas
,
S.
Whitelam
,
J.
Banfield
, and
J.
De Yoreo
, “
Microscopic evidence for liquid-liquid separation in supersaturated CaCO3 solutions
,”
Science
341
,
885
889
(
2013
).
18.
J. D.
Shore
,
D.
Perchak
, and
Y.
Shnidman
, “
Simulations of the nucleation of AgBr from solution
,”
J. Chem. Phys.
113
,
6276
6284
(
2000
).
19.
N. E. R.
Zimmermann
,
B.
Vorselaars
,
D.
Quigley
, and
B.
Peters
, “
Nucleation of NaCl from aqueous solution: Critical sizes, ion-attachment kinetics, and rates
,”
J. Am. Chem. Soc.
137
,
13352
13361
(
2015
).
20.
R.
Demichelis
,
P.
Raiteri
,
J. D.
Gale
,
D.
Quigley
, and
D.
Gebauer
, “
Stable prenucleation mineral clusters are liquid-like ionic polymers
,”
Nat. Commun.
2
,
590
(
2011
).
21.
D.
Chakraborty
and
G. N.
Patey
, “
How crystals nucleate and grow in aqueous NaCl solution
,”
J. Phys. Chem. Lett.
4
,
573
578
(
2013
).
22.
M.
Salvalaglio
,
M.
Mazzotti
, and
M.
Parrinello
, “
Urea homogeneous nucleation mechanism is solvent dependent
,”
Faraday Discuss.
179
,
291
307
(
2015
).
23.
X.
Zhang
,
Z.
Wang
,
X.
Dong
,
D.
Wang
, and
C. C.
Han
, “
Interplay between two phase transitions: Crystallization and liquid-liquid phase separation in a polyolefin blend
,”
J. Chem. Phys.
125
,
024907
(
2006
).
24.
F.
Arai
,
H.
Takeshita
,
M.
Dobashi
,
K.
Takenaka
,
M.
Miya
, and
T.
Shiomi
,
“Effects of liquid–liquid phase separation on crystallization of poly(ethylene glycol) in blends with isotactic poly(methyl methacrylate)
,”
Polymer
53
,
851
856
(
2012
).
25.
N. D.
Loh
,
S.
Sen
,
M.
Bosman
,
S. F.
Tan
,
J.
Zhong
,
C. A.
Nijhuis
,
P.
Král
,
P.
Matsudaira
, and
U.
Mirsaidov
, “
Multistep nucleation of nanocrystals in aqueous solution
,”
Nat. Chem.
9
,
77
82
(
2016
).
26.
P. G.
Vekilov
, “
The two-step mechanism of nucleation of crystals in solution
,”
Nanoscale
2
,
2346
(
2010
).
27.
N.
Asherie
, “
Phase diagram of colloidal solutions
,”
Phys. Rev. Lett.
77
,
4832
4835
(
1996
).
28.
P. E.
Bonnett
,
K. J.
Carpenter
,
S.
Dawson
, and
R. J.
Davey
, “
Solution crystallisation via a submerged liquid–liquid phase boundary: Oiling out
,”
Chem. Commun.
0
,
698
699
(
2003
).
29.
M. A.
Bewernitz
,
D.
Gebauer
,
J.
Long
,
H.
Cölfen
, and
L. B.
Gower
, “
A metastable liquid precursor phase of calcium carbonate and its interactions with polyaspartate
,”
Faraday Discuss.
159
,
291
(
2012
).
30.
X.
Wang
,
I.-M.
Chou
,
W.
Hu
, and
R. C.
Burruss
, “
In situ observations of liquid–liquid phase separation in aqueous MgSO4 solutions: Geological and geochemical implications
,”
Geochim. Cosmochim. Acta
103
,
1
10
(
2013
).
31.
L. B.
Gower
and
D. J.
Odom
, “
Deposition of calcium carbonate films by a polymer-induced liquid-precursor (PILP) process
,”
J. Cryst. Growth
210
,
719
734
(
2000
).
32.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
, “
The missing term in effective pair potentials
,”
J. Phys. Chem.
91
,
6269
6271
(
1987
).
33.
I. S.
Joung
and
T. E.
Cheatham
, “
Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations
,”
J. Phys. Chem. B
112
,
9020
9041
(
2008
).
34.
G. A.
Orozco
,
O. A.
Moultos
,
H.
Jiang
,
I. G.
Economou
, and
A. Z.
Panagiotopoulos
, “
Molecular simulation of thermodynamic and transport properties for the H2O + NaCl system
,”
J. Chem. Phys.
141
,
234507
(
2014
).
35.
Z.
Mester
and
A. Z.
Panagiotopoulos
, “
Temperature-dependent solubilities and mean ionic activity coefficients of alkali halides in water from molecular dynamics simulations
,”
J. Chem. Phys.
143
,
044505
(
2015
).
36.
J. R.
Espinosa
,
J. M.
Young
,
H.
Jiang
,
D.
Gupta
,
C.
Vega
,
E.
Sanz
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
On the calculation of solubilities via direct coexistence simulations: Investigation of NaCl aqueous solutions and Lennard-Jones binary mixtures
,”
J. Chem. Phys.
145
,
154111
(
2016
).
37.
H.
Jiang
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
Communication: Nucleation rates of supersaturated aqueous NaCl using a polarizable force field
,”
J. Chem. Phys.
149
,
141102
(
2018
).
38.
N. E. R.
Zimmermann
,
B.
Vorselaars
,
J. R.
Espinosa
,
D.
Quigley
,
W. R.
Smith
,
E.
Sanz
,
C.
Vega
, and
B.
Peters
, “
NaCl nucleation from brine in seeded simulations: Sources of uncertainty in rate estimates
,”
J. Chem. Phys.
148
,
222838
(
2018
).
39.
Z.
Mester
and
A. Z.
Panagiotopoulos
, “
Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations
,”
J. Chem. Phys.
142
,
044507
(
2015
).
40.
P. J.
Steinhardt
, “
Bond-orientational order in liquids and glasses
,”
Phys. Rev. B
28
,
784
805
(
1983
).
41.
J.
Wedekind
,
R.
Strey
, and
D.
Reguera
, “
New method to analyze simulations of activated processes
,”
J. Chem. Phys.
126
,
134103
(
2007
).
42.
J.
Wedekind
and
D.
Reguera
, “
Kinetic reconstruction of the free-energy landscape
,”
J. Phys. Chem. B
112
,
11060
11063
(
2008
).
43.
S. M. A.
Malek
,
G. P.
Morrow
, and
I.
Saika-Voivod
, “
Crystallization of Lennard-Jones nanodroplets: From near melting to deeply supercooled
,”
J. Chem. Phys.
142
,
124506
(
2015
).
44.
S. E. M.
Lundrigan
and
I.
Saika-Voivod
, “
Test of classical nucleation theory and mean first-passage time formalism on crystallization in the Lennard-Jones liquid
,”
J. Chem. Phys.
131
,
104503
(
2009
).
45.
J.
Guo
,
A.
Haji-Akbari
, and
J. C.
Palmer
, “
Hybrid Monte Carlo with LAMMPS
,”
J. Theor. Comput. Chem.
17
,
1840002
(
2018
).
46.
S.
Auer
and
D.
Frenkel
, “
Numerical prediction of absolute crystallization rates in hard-sphere colloids
,”
J. Chem. Phys.
120
,
3015
3029
(
2004
).
47.
I.
Saika-Voivod
,
F.
Romano
, and
F.
Sciortino
, “
Nucleation barriers in tetrahedral liquids spanning glassy and crystallizing regimes
,”
J. Chem. Phys.
135
,
124506
(
2011
).
48.
L. G.
MacDowell
,
P.
Virnau
,
M.
Muller
, and
K.
Binder
, “
The evaporation/condensation transition of liquid droplets
,”
J. Chem. Phys.
120
,
5293
(
2004
).
49.
S. K.
Das
and
K.
Binder
, “
Thermodynamic properties of a symmetrical binary mixture in the coexistence region
,”
Phys. Rev. E.
84
,
061607
(
2011
).
50.
S. K.
Das
, “
Atomistic simulations of liquid-liquid coexistence in confinement: Comparison of thermodynamics and kinetics with bulk
,”
Mol. Simul.
41
,
382
401
(
2015
).
51.
D.
Stauffer
and
K.
Binder
, “
Statistical theory of nucleation, condensation and coagulation
,”
Adv. Phys.
25
,
343
396
(
1976
).
52.
F.
Moučka
,
M.
Lísal
, and
W. R.
Smith
, “
Molecular simulation of aqueous electrolyte solubility. 3. Alkali-halide salts and their mixtures in water and in hydrochloric acid
,”
J. Phys. Chem. B
116
,
5468
5478
(
2012
).
53.
F.
Moučka
,
I.
Nezbeda
, and
W. R.
Smith
, “
Molecular force fields for aqueous electrolytes: SPC/E-compatible charged LJ sphere models and their limitations
,”
J. Chem. Phys.
138
,
154102
(
2013
).
54.
I.
Nezbeda
,
F.
Moučka
, and
W. R.
Smith
, “
Recent progress in molecular simulation of aqueous electrolytes: Force fields, chemical potentials and solubility
,”
Mol. Phys.
114
,
1665
1690
(
2016
).
55.
P. G.
Debenedetti
,
Metastable Liquids: Concepts and Principles
(
Princeton University Press
,
Princeton, NJ, USA
,
1996
).
56.
S.
Jungblut
and
C.
Dellago
, “
Crystallization of a binary Lennard-Jones mixture
,”
J. Chem. Phys.
134
,
104501
(
2011
).
57.
J.
Reichenbach
and
K.
Wynne
, “
Frustration vs prenucleation: Understanding the surprising stability of supersaturated sodium thiosulfate solutions
,”
J. Phys. Chem. B
122
,
7590
7596
(
2018
).
58.
M. K.
Mitra
and
M.
Muthukumar
, “
Theory of spinodal decomposition assisted crystallization in binary mixtures
,”
J. Chem. Phys.
132
,
184908
(
2010
).
59.
L. R.
Gómez
, “
Amorphous precursors of crystallization during spinodal decomposition
,”
Phys. Rev. E
83
,
021501
(
2011
).
60.
A.
Myerson
and
B.
Trout
, “
Nucleation from solution
,”
Science
341
,
855
856
(
2013
).
61.
J.
De Yoreo
, “
More than one pathway
,”
Nat. Mater.
12
,
284
285
(
2013
).
62.
M.
Sleutel
,
J.
Lutsko
,
A. E.
Van Driessche
,
M. A.
Durán-Olivencia
, and
D.
Maes
, “
Observing classical nucleation theory at work by monitoring phase transitions with molecular precision
,”
Nat. Commun.
5
,
5598
(
2014
).
63.
S.
Karthika
,
T. K.
Radhakrishnan
, and
P.
Kalaichelvi
, “
A review of classical and nonclassical nucleation theories
,”
Cryst. Growth Des.
16
,
6663
6681
(
2016
).
64.
Q.
Hu
,
M. H.
Nielsen
,
C. L.
Freeman
,
L. M.
Hamm
,
J.
Tao
,
J. R. I.
Lee
,
T. Y. J.
Han
,
U.
Becker
,
J. H.
Harding
,
P. M.
Dove
, and
J. J.
De Yoreo
, “
The thermodynamics of calcite nucleation at organic interfaces: Classical vs. non-classical pathways
,”
Faraday Discuss.
159
,
509
(
2012
).
65.
D.
James
,
S.
Beairsto
,
C.
Hartt
,
O.
Zavalov
,
I.
Saika-Voivod
,
R. K.
Bowles
, and
P. H.
Poole
, “
Phase transitions in fluctuations and their role in two-step nucleation
,”
J. Chem. Phys.
150
,
074501
(
2019
).

Supplementary Material