The formation of crystals has proven to be one of the most challenging phase transformations to quantitatively model—let alone to actually understand—be it by means of the latest experimental technique or the full arsenal of enhanced sampling approaches at our disposal. One of the most crucial quantities involved with the crystallization process is the nucleation rate, a single elusive number that is supposed to quantify the average probability for a nucleus of critical size to occur within a certain volume and time span. A substantial amount of effort has been devoted to attempt a connection between the crystal nucleation rates computed by means of atomistic simulations and their experimentally measured counterparts. Sadly, this endeavor almost invariably fails to some extent, with the venerable classical nucleation theory typically blamed as the main culprit. Here, we review some of the recent advances in the field, focusing on a number of perhaps more subtle details that are sometimes overlooked when computing nucleation rates. We believe it is important for the community to be aware of the full impact of aspects, such as finite size effects and slow dynamics, that often introduce inconspicuous and yet non-negligible sources of uncertainty into our simulations. In fact, it is key to obtain robust and reproducible trends to be leveraged so as to shed new light on the kinetics of a process, that of crystal nucleation, which is involved into countless practical applications, from the formulation of pharmaceutical drugs to the manufacturing of nano-electronic devices.

Understanding crystal nucleation is one of the fundamental ambitions of physical chemistry.1–4 Far from being a theoretical curiosity, the kinetics of crystallization has a crucial impact on many natural phenomena, such as the formation of ice5 or the process of biomineralization,6 and on a diverse range of practical applications, such as the formulation of pharmaceutical drugs7 or the design of novel nanostructures for data storage.8 

If we were to ask the reader to pick a single observable to characterize the nucleation process, the so-called crystal nucleation rate J would probably turn out to be a very popular choice. J is a scalar quantity that represents the average probability, per unit time and per unit volume, for a critical crystalline nucleus to occur within the supercooled liquid or the supersaturated solution of interest. The apparently straightforward nature of J makes it ideal to compare the crystallization kinetics of the same system in different conditions or indeed between different systems as well. Crucially, J can be both measured experimentally and estimated by means of computer simulations, thus providing, in principle, a much sought after connection between reality and modeling. However, long-standing major inconsistencies between experiments and simulations still persist as of today, despite the ever-growing capabilities of molecular simulations.3,4

When dealing with estimates of crystal nucleation rates obtained by means of atomistic simulations, the community has identified several outstanding issues through the years. Among the usual suspects, we can find the intrinsic limitations of the models we use to describe the interactions between atoms, or particles, or molecules.9–11 In addition, the rather dated theoretical framework provided by classical nucleation theory (CNT), while proving to be remarkably accurate in many cases, has been repeatedly put into question within the last few years.12,13

Here, we are going to focus on some (seven) aspects of interest for molecular simulations aimed at computing the crystal nucleation rate. Some of these issues, such as the potentially slow dynamics of the system,14 are often overlooked. Some others, for instance, the existence of finite size effects,15 are very well known—and yet, the extent of their impact on the estimate of J is still largely unexplored. Understanding the contributing factors to the uncertainty characterizing J is crucial, not only to keep working toward the perhaps overly ambitious goal of reconciling simulations with experiments but also to allow meaningful insight to be extracted from our simulations in the first place. Indeed, these are really exciting times for the community, as a number of recent works have provided an unprecedented level of detail into specific aspects (negatively) affecting the calculation of J. From the necessity to account for the microscopic kinetics of the system to the need of thoroughly investigating the adequacy of the order parameters we use, we take stock of the state-of-the-art and offer an opinionated perspective as to potential future developments within the field.

This paper is organized as follows: after having set the boundaries of our discussion and briefly reviewed the computational tools out our disposal as well as the above mentioned usual suspects (Sec. II), we focus on seven different aspects that are detrimental to the calculation of J by means of enhanced sampling atomistic simulations (Sec. III). Finally, in Sec. IV, we assess the challenges ahead and put forward several possibilities we hope the community will decide to explore in the near future.

In order to discuss the fundamental aspects of calculating nucleation rates via atomistic simulations, we have to remove as many layers of complexity as possible from both the system under investigation and the conditions in which we are working on it. We start by focusing on supercooling or supersaturation regimes that lie well away from the spinodal limit where the free energy barrier associated with the formation of a critical nucleus is vanishingly small. Conversely, we have to avoid, by necessity, very mild supercooling or supersaturation as well; in that case, the size of the critical nucleus is usually too large to be even taken into account by means of atomistic simulations. Then, we shall consider almost exclusively homogeneous nucleation, in principle, the simplest scenario available to us. However, homogeneous nucleation is rarely observed in reality, with heterogeneous nucleation being much more common. Indeed, it is often very challenging to experimentally measure crystal nucleation rates without any influence in terms of impurities, which almost always manage to facilitate the crystallization process.3,4 Note, however, that much of what is discussed in the following is equally relevant to studies of heterogeneous nucleation. In addition, we will largely avoid the emergence of two- or even multi-step nucleation processes (biomineralization being a prominent example) as well as confinement effects. Excellent reviews on these topics and, broadly speaking, on the subject of crystal nucleation as a whole can be found in, e.g., Refs. 2, 3, 13, 17, and 18.

At this stage, it is useful to highlight the distinction between “realistic” and “model” systems. Simulating the former usually represents an attempt to get as close as possible to the experimental reality by utilizing the most accurate interatomic potentials available. Note that leveraging ab initio simulations is, aside for very rare exceptions,8 simply not feasible, as both the time- and length-scales involved in the crystal nucleation process (particularly taking into account its stochastic nature) lie far beyond the reach of these accurate, but rather computationally expensive approaches. Machine learning (ML)-based interatomic potentials19 represent an especially intriguing avenue to strike the balance between the computational efficiency of classical potentials and the accuracy of first principles calculations; however, we argue that, aside for exceptional cases where the kinetics of crystallization is unusually fast (see, e.g., Ref. 20), these novel potential are still too computationally costly to be used to compute crystal nucleation rates with sufficient accuracy. Even when choosing the most accurate classical potentials at our disposal, the computational cost involved with the calculation of nucleation rates increases drastically, usually imposing a tough compromise between the quality of the simulations and the quantity of the results. In fact, the majority of rate calculations for realistic systems (notable examples would be Refs. 10, 21, and 22) focus on a single supercooling or supersaturation condition. Given that reaching an agreement between experimentally measured and simulated nucleation rates is, as we shall see shortly, very challenging, this is obviously sub-optimal.

On the other hand, computing crystal nucleation rates for “model” systems aims at employing simple and/or computational inexpensive interatomic potentials, so as to explore different, e.g., supercooling conditions and thus extract useful trends to be compared with the experimental data. Clearly, this often comes at the cost of sacrificing some accuracy in terms of the description of the actual system. However, it is important to point out that, in some cases, simple potentials can describe experimentally relevant systems to a very reasonable level of accuracy. Representative examples of such systems that we will often reference in this work are the Lennard-Jones (LJ) system and colloidal particles. These have all been extensively investigated through the years and thus offer the perfect opportunity to discuss several aspect related to the calculation of nucleation rates via atomistic simulations.

One issue that unfortunately unifies the vast majority of computational studies, not necessarily limited to those dealing with crystal nucleation rates, is the reproducibility of the results. This is especially evident in the case of LJ systems, where minimal modifications to details such as the truncation/shift (which vary wildly within the literature) lead to substantially different melting curves and nucleation rates. In fact, as we shall discuss in greater detail in Sec. II B, the estimate of J is incredibly susceptible to each and every feature of the interatomic potential of choice. In this context, it is comforting to witness the emergence of more and more open access databases where the actual input files used to obtain a particular computational results are reported in full (see, e.g., Ref. 23), thus greatly facilitating cross-validation.

In contrast, one key issue that stubbornly refuses to yield to the efforts of the community is the long-standing discrepancy between experimental and simulated crystal nucleation rates. Perhaps the starkest reminder of how much work remains to be done in the field is given by the current state of affairs with respect to colloidal systems. These qualify as “simple” systems, given that they can, in principle, be modeled via very inexpensive models, such as hard or soft spheres. Crucially, it is also possible to measure crystal nucleation rates of colloidal systems via relatively straightforward experimental techniques such as confocal microscopy.24 As can be seen in Fig. 1, however, while experiments and simulations agree quite well at strong supersaturation, discrepancies as large as ten orders of magnitude can be found at weak supersaturation, thus indicating that, perhaps, our simple models cannot, in fact, be used to simulate these systems accurately enough. This evidence is quite astonishing, as we have by now reached a point where different approaches to the calculation of J are consistent with each other. In particular, the simulation results reported in Fig. 1 for weak supersaturation (obtained by Fiorucci et al. in 202016) have been validated by three different methodologies. The fact that the discrepancy between experiments and simulations still persists is thus suggestive of the possibility that there must be some fundamental aspect of our simulations that we are still missing—and that has nothing to do with the choice of the method used. Potential culprits are the polydispersity and the deformability of the particles or effects related to sedimentation and hydrodynamics, albeit the latter has been very recently ruled out.16 

FIG. 1.

The inconsistencies between experimental and simulated nucleation rates persists. Nucleation rate J σ eff 5 / D 0 as a function of the volume fraction ϕ of a system of “nearly” hard spheres. This is often considered as a “model” system, given the relatively straightforward interactions between the constituent particles. Experimental results are reported in green, while the outcomes of a variety of simulations are reported in blue or red. Note the substantial discrepancy (∼10 orders of magnitude) between experiments and simulations in the weak supersaturation regime. The data points in red have been obtained in Ref. 16 by means of three different approaches and are impressively consistent with each other, which is indicative of the fact that said discrepancy has nothing to do with the methodology employed. Reprinted from Fiorucci et al., J. Chem. Phys. 152, 064903 (2020) with the permission of AIP Publishing LLC.

FIG. 1.

The inconsistencies between experimental and simulated nucleation rates persists. Nucleation rate J σ eff 5 / D 0 as a function of the volume fraction ϕ of a system of “nearly” hard spheres. This is often considered as a “model” system, given the relatively straightforward interactions between the constituent particles. Experimental results are reported in green, while the outcomes of a variety of simulations are reported in blue or red. Note the substantial discrepancy (∼10 orders of magnitude) between experiments and simulations in the weak supersaturation regime. The data points in red have been obtained in Ref. 16 by means of three different approaches and are impressively consistent with each other, which is indicative of the fact that said discrepancy has nothing to do with the methodology employed. Reprinted from Fiorucci et al., J. Chem. Phys. 152, 064903 (2020) with the permission of AIP Publishing LLC.

Close modal

What options are available to the computational scientist in order to calculate crystal nucleation rate by means of atomistic simulations? To start with, the so-called “brute force” simulations constitute the most straightforward option: the idea is to bring the simulated liquid (or solution) into a supercooled (or supersaturated) state and simply run a set of unbiased simulations in the hope that the system will crystallize on a reasonably short timescale. At that point, methodologies such as the mean first-passage time (MFPT) can be used to extract the nucleation rate.25,26 Brute force simulations are regarded as useful benchmarks to compare the results of enhanced sampling simulations with. However, not very many systems will, in fact, crystallize quickly enough, under experimentally relevant conditions, for this method to be applied. The reality is that, in most cases, enforcing spectacularly strong supercooling or supersaturation conditions is a necessity—and at that point, one should start questioning whether the nucleation events we are seeing do, in fact, conform to the definition of nucleation as a rare event, an aspect we will address in more detail in Sec. III A. More often than not, brute force simulations aimed at calculating crystal nucleation rates are used for model systems such as LJ liquids or colloids.27,28 However, particularly impressive results have been obtained of late for metallic alloys as well,29,30 in some cases via simulations involving millions of atoms.31 

In order to compute crystal nucleation rates for realistic systems or indeed for model systems within a wider range of supercooling or supersaturation, enhanced sampling techniques represent the tools of the trade. These can be classified as free energy-based or path sampling-based techniques. The former seeks to obtain information about the thermodynamics of the process, and as such, they do not offer immediate access to the kinetics of nucleation. However, it is possible, in principle, to use the outcomes of, e.g., umbrella sampling32–35 or metadynamics simulations36–39 (most prominently, the free energy profile or surface) as the starting point to compute the nucleation rate. Popular options in this context include the family of Bennett–Chandler methods for the calculation of rate constants40 and the approach recently pioneered by Tiwary and Parrinello.41,42

In contrast, path sampling-based techniques, such as transition path sampling (TPS),43–45 transition interface sampling (TIS),46–48 and forward flux sampling (FFS),10,21,22,49,50 all allow for the direct calculation of nucleation rates. Interestingly, until quite recently, not many examples of any of these methodologies being applied to the calculation of J could be found in the literature, most likely due to the high computational cost associated with converging these algorithms. However, the ever-growing capabilities of high performance computing facilities, together with the trivially parallel nature of the simulations often needed to sample the crystallization paths, have massively boosted the feasibility of path-sampling methods in the last few years.

Perhaps unsurprisingly, ML has started to help in the context of enhanced sampling methods as well—not so much as to drive the nucleation process itself but to aid and complement the structural analysis of the newborn nuclei as well as the nucleation paths. Recent examples include the work of Desgranges and Delhommelle,51 where ML is used to infer the free energy of a LJ system over a wide range of densities and temperatures; the work of Bonati and Parrinello,52 in which ML is used to train a model capable of describing the crystal nucleation of silicon, starting from a dataset obtained by metadynamics simulations; and the recent work of Adorf et al.,53 where ML has been harnessed to analyze the nucleation pathways of colloidal systems.

Finally, seeding techniques54–56 can also be used as a starting point to extract crystal nucleation rates. These approaches relies on sampling the evolution of sets of crystalline nuclei of a given size that are inserted into the supercooled liquid or supersaturated solution. Particularly famous examples include the nucleation of ice57 and NaCl.9 

The inconsistencies between experimental and simulated crystal nucleation rates, such as those we have discussed in the case of colloidal systems in Sec. II, are often thought to be due to the usage of CNT and/or a particular force field. The two are not unrelated either, in that several thermodynamic parameters needed to estimate J via CNT, such as the chemical potential difference between the liquid and the crystal, Δμ, and crystal-liquid interfacial free energy γ, vary wildly according to the choice of a given force field. Several attempts to improve upon the value of Δμ, provided the force field in question, can be found in the recent literature. For instance, Wang et al.58 used a second-order Gibbs–Thompson equation to obtain a reliable description of the temperature dependence of Δμ in the case of ice nucleation. In other studies, the concept of an “effective” γ has been adopted in, e.g., Ref. 9 to describe the nucleation of NaCl from solution, instead of the often-used value of γ referring to a macroscopic flat interface. In fact, evidence from the field of ice nucleation confirms that values of γ for finite nuclei extracted from simulations via seeding59 or umbrella sampling60 are indeed quantitatively different from those computed from planar interfaces.61–63 Taking into account the Tolman correction64 contributes (partially) to bridging the difference, noting that the Tolman length encapsulates only the first term in an expansion of a difference, which may not be small.

Indeed, a major point of contention with respect to CNT is its usage at the microscopic scale, where the distinction between solid and liquid phases is not clear cut and the actual value of the thermodynamic parameters involved is bound to differ from their macroscopic counterparts. Even when applying this correction, however, discrepancies of about 15 orders of magnitude with respect to J can still be found in Ref. 9. The Joung–Cheatham force field65 for NaCl and the extended simple point charge model (SPC/E) used for water might be partially responsible for that discrepancy still, the authors suggest—but it has to be said that the choice is superior to that adopted in previous works where the GROMOS85 force field, which overestimates the stability of the NaCl wurtzite structure, was used.66 In Ref. 9, it is acknowledged that no current force field for NaCl will perfectly capture the chemical potential difference between ions in solution and ions encapsulated in a bulk crystal and hence reproduce exactly the experimental solubility. However, it is argued that provided one compares between simulation and experiment at the same value of this chemical potential difference (but different absolute concentration), one should be able to draw meaningful conclusions. This requires the crystal solubility within the chosen force field to be known accurately. In the case of the Joung–Cheatham model, various estimates of this solubility have been made,67–69 eventually reaching a consensus70 subsequently confirmed by increasingly advanced methods71 to be approximately half of the experimental value. This indicates that force fields tractable to nucleation simulations are far from sufficiently accurate in absolute terms. Later work72 re-examined the results of Ref. 9 in light of the revised estimates of model solubility, resulting in much improved agreement with experimental nucleation rates when compared at the equivalent chemical potential difference. That is, the authors found that reliable results in terms of J can be obtained when the driving force for nucleation, Δμ, and the solubility are consistent with each other and accurate with respect to the particular model/force field used. This means that even if J is in agreement with the experiments, the actual solubility of the force field is, in this case, not consistent with the experimental value at that particular supersaturation. This work also demonstrated a very strong sensitivity to the choice of order parameter used to quantify the nucleus size. We return to this question in Sec. III F.

To illustrate how sensitive the calculation of J is to the inaccuracies of the force field, we refer to the work of Haji-Akbari and Debenedetti,10 where the ice nucleation rate has been computed via FFS using the TIP4P/ice water model at the very strong supercooling of 42 K. This choice might appear rather extreme but it is not unusual: when dealing with realistic and/or computationally expensive force fields, resorting to very strong supercooling or supersaturation is often the only way to gain any insight into the nucleation process, even when using state-of-the-art enhanced sampling methods. This is important because in this scenario, one has to extrapolate the value of J at milder supercooling in order to be able to compare the computational result to the experimental reality. Once again, then, CNT comes into play and with it multiple potential sources of inaccuracies. In the case of Ref. 10, the authors observe a discrepancy of about eight orders of magnitude between simulated and experimental nucleation rates—a difference that can be explained, they argue, by noting that the TIP4P/ice water model yields a value of Δμ that, corresponds to an enthalpy difference about 20% smaller than the experimental value. Similar arguments can be found in other recent works. For instance, Arjun et al.73 put forward the usage of a specific water model as the reason why the free energy barriers relative to the nucleation of methane hydrates, obtained via TPS, are found to be lower with respect to previous results.

Indeed, the reliability of water models is a major issue when dealing with crystal nucleation from solution. We have already mentioned the case of NaCl, but perhaps the most prominent scenario in this context is that of biomineralization, where the water model has to capture the complex interaction with the mineral under investigation. In their recent review, Demichelis et al.11 argued that an accurate description of the solubility of the mineral, a quantity often neglected in the first parameterizations of force fields for simulations of biomineralization, is paramount to obtain robust results. The additional layer of complexity in that field is the scarcity of experimental data, particularly in terms of clustering and speciation (under conditions accessible by both experiments and simulations) available for the computational scientists to build their force field upon. Electronic structure calculations have been used to fill that gap, which, in turn, highlighted the absolute need for polarizable force fields when dealing with the crystal nucleation of biominerals.11 Unfortunately, the additional computational expense of such models makes quantitative calculation of J intractable at this time.

Another example, yet again from the ice nucleation field, is given by the recent work of Shi and Tanaka,74 where the authors have found that in the case of the TIP5P model, the treatment of electrostatic interactions has a huge impact on the nucleation process, to the point where results massively favored a specific (ferroelectric) crystalline phase. The fact that truncating electrostatic or indeed even non-bonded interactions can have a substantial effect on the property of the force field—and thus, indirectly, on the estimate of J—is well known (see, e.g., Ref. 75) but not necessarily discussed in a transparent, reproducible manner in all computational studies so far. In the same work,74 Shi and Tanaka also pointed out that while CNT assumes the liquid to be a perfectly homogeneous phase, this might not be the case.

The fact that supercooled liquids and supersaturated solutions exhibit both structural and dynamical heterogeneities is well established, but only in the last few years, we have started to witness the first attempts to make a connection between the pre-ordering of the liquid and crystal nucleation. As an example, Menon et al.76 recently investigated the crystal nucleation of molybdenum via TIS simulations, finding that the emergence of crystal-like precursors plays a role in polymorph selection. As it concerns dynamical properties, Fitzner et al.77 recently established a link between the dynamical heterogeneity of supercooled liquid water and the occurrence of ice nucleation.

In summary, we hope we have convinced the reader that the accuracy of the force fields we use is intertwined with the reliability of our CNT predictions. However, we argue that this picture is far from complete and that there exist several additional issues that might come into play when computing J via molecular simulations. In Sec. III, we will highlight seven such aspects, presenting some relevant examples from the recent literature with the aim to bring together the efforts of the community toward an ever-increasingly accurate picture of the crystal nucleation process from a microscopic perspective.

The central idea at the heart of the CNT kinetics was developed by Becker and Döring in 1935:78,79 the time evolution of the distribution of crystalline clusters is treated via a formalism equivalent to that of chemical rate equations—under a number of assumptions. For instance, they assumed that the nuclei either shrink or grow via losing or gaining a single atom, particle, or molecule to the existing nucleus. Most relevant to this section, though, is the assumption of a quasi-stationary distribution (QSD) of nuclei,80 which is supposed to not get depleted in time by the formation of critical nuclei—once critical clusters form, they are removed from the distribution and new smaller nuclei are added, so as to achieve a steady state. Provided that the free energy barrier associated with the formation of critical nucleus is “high enough,” the supercooled liquid (or the supersaturated solution) is to be found in a metastable state with respect to the crystalline phase within a timescale much longer than its relaxation time, thus allowing us to compute J as a steady-state crystal nucleation rate.

In this scenario, we are dealing with a rare event characterized by a survival probability Pliq(t*) for the supercooled liquid or supersaturated solution that decreases exponentially with time t as, e.g., Pliq(t*) = exp(−Jt*). In other words, the nucleation times (to be observed across a sufficiently large ensemble of either experiment or simulations, given the stochastic nature of the nucleation process) are distributed according to Poisson statistics. In the context of unbiased simulations of crystal nucleation, verifying this condition is relatively straightforward, as illustrated in Fig. 2 where we report the time evolution of 220–820 (see the  Appendix for further details) molecular dynamics (MD) trajectories of a LJ system at different supercooling and the corresponding survival probability for the liquid phase. At very strong supercooling [Fig. 2(a)], the free energy barrier is so low that the liquid is basically unstable, as opposed to metastable, with respect to the liquid phase. As such, there is no incubation (or waiting) time, as the system immediately proceeds to crystallize on a timescale comparable to that of its relaxation time. The resulting survival probability is thus a step function which tell us that, in these conditions, we are not dealing with a rare event in the first place. At milder supercooling [Fig. 2(b)], we reach a “butter zone” (i.e., a temperature regime ideally suited to extract, in this case J) where we can sample P(t) across a sufficiently large timescale, long tails included. This is an ideal setting to accurately compute J. Finally, at very mild supercooling [Fig. 2(c)], we are obviously dealing with even rarer events—however, the timescale associated with the nucleation process is such that we can only sample a very small portion of the relevant P(t), thus robbing us from the possibility of investigating the actual decay of this function and assessing whether it is truly consistent with a Poisson distribution.

FIG. 2.

Crystal nucleation of a LJ liquid at different supercooling. The freezing curves (left panels, shifted so that the nucleation times lie at Δ t* = 0) and survival probabilities Pliq(t*) (right panels, as a function of time elapsed, t*) for the same LJ system at different (reduced) temperatures T*. (a) Strong supercooling (T* = 0.68): the liquid is unstable with respect to the crystalline phase. (b) “Butter zone” (T* = 0.765): it is possible to sample the exponential decay of Pliq(t*) across a sufficiently long timescale. (c) Mild supercooling (T* = 0.8): it is possible to accumulate statistics with respect to the nucleation times only within a very limited timescale compared to the decay of Pliq(t*). We have fitted the survival probabilities according to the CNT-like expression Pliq(t*) = exp[−(Jt*)α], where α is a fitting parameter accounting for the possibility of non-exponential nucleation kinetics. The uncertainty associated with the unconstrained fit is shown as shaded regions. The details of the LJ simulations are given in the  Appendix.

FIG. 2.

Crystal nucleation of a LJ liquid at different supercooling. The freezing curves (left panels, shifted so that the nucleation times lie at Δ t* = 0) and survival probabilities Pliq(t*) (right panels, as a function of time elapsed, t*) for the same LJ system at different (reduced) temperatures T*. (a) Strong supercooling (T* = 0.68): the liquid is unstable with respect to the crystalline phase. (b) “Butter zone” (T* = 0.765): it is possible to sample the exponential decay of Pliq(t*) across a sufficiently long timescale. (c) Mild supercooling (T* = 0.8): it is possible to accumulate statistics with respect to the nucleation times only within a very limited timescale compared to the decay of Pliq(t*). We have fitted the survival probabilities according to the CNT-like expression Pliq(t*) = exp[−(Jt*)α], where α is a fitting parameter accounting for the possibility of non-exponential nucleation kinetics. The uncertainty associated with the unconstrained fit is shown as shaded regions. The details of the LJ simulations are given in the  Appendix.

Close modal

It is important to point out that this “butter zone” might simply not exist for some systems, particularly slow-diffusing ones for which the balance between supercooling (or supersaturation) and atomic/molecular mobility results in timescales too long to be probed by means of unbiased simulations even at strong supercooling. The opposite problem, that of “not rare enough” events, is also encountered experimentally, and it might coincide with the onset of spinodal decomposition. An intriguing take on this matter is given by Sear,81 who provocatively puts forward the idea that “the nucleation rate may not exist,” given that, by looking at the experimental data, the survival probabilities measured for a variety of systems are far from being exponential.

This is especially true when considering heterogeneous nucleation, where structural or even chemical changes in the impurity involved with the nucleation process often result in time-dependence nucleation regimes as opposed to the steady state J we are after. Fortunately, several studies aimed at rationalizing the experimental data available to us can be found in the recent literature. As an example, Maggioni and Mazzotti82 looked into an extensive set of literature data on the crystallization of p-aminobenzoic acid in three different solvents, developing a statistical analysis that allows us to quantify the uncertainty of experimental nucleation data—a crucial aspect to facilitate the comparison with simulated nucleation rates.

An interesting aspect that unifies the majority of the modeling of nucleation data is the assumption that nucleation can be considered as a Markovian process—that is, there is no history dependence. However, it is becoming apparent that this assumption might not always hold. For instance, Jungblut and Dellago83 found that the nucleation dynamics of a particular LJ system shows non-Markovian aspects due to the lack of a good enough reaction coordinate (RC), an aspect we will discuss in greater detail in Sec. III F and that results in underestimating the nucleation rate if using MFPT methods. Indeed, the recent work of Kuhnhold et al.,14 also investigating the crystal nucleation of a LJ melt, highlighted the emergence of non-Markovian dynamics and went as far as saying that CNT has to be considered as the limiting case of a more general theoretical framework that includes history-dependent aspects. In the words of the authors, crystal nucleation may be “neither Markovian nor diffusive.” It is important to point out that the term “non-Markovian” can refer to two distinct issues, both of which will be discussed in this work. First, the occurrence of nucleation events itself might or might not be a Markovian process. But then, there is the question of whether the description of the crystallization process, by means of the evolution of the order parameters we shall discuss in Sec. III F, is Markovian or not.

At this stage, we should ask ourselves how would enhanced sampling fare in extracting accurate nucleation rates, given that even unbiased simulations struggle in some cases to recover reliable estimates of J. A central argument in this context is quantifying the uncertainty associated with any given methodology. The discussion contained in Ref. 84 highlights once more the importance of sampling a QSD: in particular, the work of Binder et al.85 on parallel replica dynamics offers some practical considerations that allow the computational scientist to analyze the obtained nucleation times. If the latter are exponentially distributed, we might be working in conditions close enough to a QSD—but not necessarily. Cases where the P(t) is exponential but the Markov process did not reach the QSD are still a possibility. Similar arguments apply to the case of extracting nucleation rates starting from the “infrequent metadynamics” framework of Tiwary and Parrinello.86 In particular, Salvalaglio et al. devised a simple and computationally inexpensive methodology, based on the Kolmogorov–Smirnov test, to assess the reliability of the kinetics of nucleation obtained from metadynamics simulations. This framework identifies several problematic cases where, e.g., the applied bias is strong enough for the system to move away from the QSD.41,87,88

The fact that finite size effects can have a substantial impact on the estimate of J should hopefully come as no surprise. A frequently encountered rule of thumb consist in working with a simulation box larger than two times the extent the critical nucleus size in any given dimension. However, the reality is that, without a thorough investigation of finite size effects using, e.g., increasingly large simulation boxes, the extent of this effect is impossible to quantity a priori. This is a real issue, as in many cases, one can rarely afford the computer time needed to repeat their simulations with larger simulation boxes—particularly when computing the nucleation rate for realistic systems. In addition, the impact of finite size effects is related to the supercooling or supersaturation conditions as well. The case of, e.g., mild supercooling is very well-known: as the critical nucleus size becomes (exponentially) larger as the supercooling decreases, larger boxes are, by definition, a necessity. However, at strong supercooling, the average density of pre-critical nuclei is much higher than what is observed at mild supercooling, thus posing the problem of the interaction between different smaller nuclei (as opposed to the very rare/large critical nuclei observed at mild supercooling). Finite size effects are also strongly connected to the issue of solute depletion we will discuss in Sec. III C.

Attempts to quantify the severity of finite size effects in simulations of crystal nucleation date back to the seminal work of Honeycutt and Andersen on LJ systems in 1986,89 where the density of critical nuclei was investigated. The actual influence of finite size effects on the calculation of J has been quantified much later, though—the work of Huitema et al. in 2000 being an excellent example.90 By then, the typical system size was in the region of 104 LJ particles. These days, however, in many cases, we can afford to investigate much larger systems. For example, Ouyang et al. recently probed the finite size effects on the crystal growth rate of LJ systems using models containing up to 106 LJ particles,28 taking advantage of graphics processing unit (GPU)-accelerated molecular dynamics simulations.

Perhaps one of the most striking results within the recent literature concerning finite size effects in the context of crystal nucleation rates is the work of Mahata and Asle Zaeem,91 which focuses on the crystal nucleation of elemental metals described via the second nearest-neighbor modified embedded atom method (2NN-MEAM) interatomic potential. As illustrated in Fig. 3, the authors considered models containing up to 107 atoms, providing reliable estimates of J as a function of system size. While many would probably argue that discrepancies barely spanning a single order of magnitude are to be considered as small in the context of J estimates, it is intriguing to observe that finite size effects are still very much present even in systems containing millions of atoms. In addition, the results reported in Fig. 3 show a non-monotonic dependence of J with respect to the system size. This is an intriguing finding, in that it is rather common to assume that the presence of finite size effects tends to consistently overestimate the nucleation rate. However, it appears that the picture is more nuanced—and certainly largely unexplored still. Indeed, the recent work of Hussain and Haji-Akbari,15 which thoroughly explored the impact of finite size effects in the context of simulations of heterogeneous ice nucleation, shows a positive linear correlation between the (log10 of) the nucleation rate and the system size. The same work provides a number of practical guidelines, particularly in terms of quantifying the spurious interactions between pre-critical nuclei due to small simulation boxes and how to differentiate such interactions from those originating from very strong supercooling regimes. These concepts can and should be considered when dealing with simulations of homogeneous nucleation as well.

FIG. 3.

The impact of finite size effect on the estimate of crystal nucleation rates. Nucleation rate for atomistic models of Al, obtained via mean the first-passage time (MFPT) method, for different system sizes. These simulations took advantage of the computational efficiency of the second nearest-neighbor modified embedded atom method (2NN-MEAM) interatomic potential to explore finite size effects arising when dealing with up to several millions of atoms. Adapted from Ref. 91.

FIG. 3.

The impact of finite size effect on the estimate of crystal nucleation rates. Nucleation rate for atomistic models of Al, obtained via mean the first-passage time (MFPT) method, for different system sizes. These simulations took advantage of the computational efficiency of the second nearest-neighbor modified embedded atom method (2NN-MEAM) interatomic potential to explore finite size effects arising when dealing with up to several millions of atoms. Adapted from Ref. 91.

Close modal

An interesting aspect of finite size effects is that the latter are almost always thought of in terms of structural correlations within the supercooled liquid or supersaturated solution. However, recent evidence suggests that dynamical correlations might also play a role (see the excellent review of Zanotto and Montazerian92). For instance, the work of Fitzner et al.77 established a correlation between dynamical heterogeneities within supercooled liquid water and the emergence of ice nuclei, and similar findings have been reported for a diverse portfolio of systems, from LJ liquids93 to oxides94,95 and metallic glass formers as well.96 Dynamical correlations can span much larger length-scales than the structural ones, thus prompting the question of whether up-to-now undetected finite size effects in terms of dynamical properties might be present in our simulations of crystal nucleation—a largely unexplored possibility.

Finally, finite size effects are important when computing the solubility of, e.g., ions in solution,55 which, in turn, is key to calculate J if leveraging CNT via, for instance, seeded approaches. The case of NaCl is especially relevant and will thus be discussed in Sec. III C, where we will tackle the emergence of solute depletion effects—an occurrence very much related to the size of the simulation box and thus to the possibility of finite size effects as well.

The molecular simulator is usually interested in quantifying the nucleation rate by studying a small volume of the system, representative of the bulk metastable parent phase. Simulations should ideally capture the coupling of this volume to its environment (i.e., the rest of the bulk system) via exchange of both particles and heat. Most simulations of nucleation from solution have been conducted with a constant number of solute particles. The nucleation of a solid phase when working with relatively small simulation boxes leads to a change in the effective supersaturation of the simulation due to the fact that the crystalline phase nucleates at the expenses of the solution phase, thus depleting the latter to an extent proportional to the number of particles/atoms/molecules contained in the simulation box. Real systems are large enough to compensate the resulting change in terms of chemical potential, but in atomistic simulations, this is a serious issue, which is less evident (albeit still present and, we argue, rather unexplored up until now) for crystal nucleation from a supercooled liquid (where as opposed to “depletion,” the emphasis is on the density difference between the liquid and the crystal) but quite spectacularly evident for crystal nucleation from a supersaturated solution. As an example, in Fig. 4(a), it can be seen that the free energy surface relative to the formation of crystalline NaCl from solution simply lacks the basin corresponding to the crystalline phase when performing simulations within the NVT ensemble. It is important to note, however, that solute depletion is not a prerogative of simulations and it can occur in actual nucleation processes as well. Depleted regions around the growing crystalline nuclei are also present in reality: they might be replenished quickly enough to conform to the CNT assumption of a constant background supersaturation, but they might also be replenished very slowly (depending on, e.g., the diffusivity of the solute)—a situation closer to simulations but not consistent with CNT.

FIG. 4.

Simulating the crystal nucleation at constant chemical potential. Free energy surfaces, obtained by means of metadynamics simulations, relative to the nucleation of NaCl crystals from solution. The two order parameters SO and SH provide a measure of the crystalline order and the hydration number of Na+ ions, respectively. (a) NVT ensemble. (b) Constant chemical potential via the CμMD method discussed in Ref. 101. (c) Same as (b) but at higher supersaturation. Note that the conventional NVT setup fails to identify the basin corresponding to the crystalline phase. Reprinted with permission from Karmakar et al., J. Chem. Theory Comput. 15, 6923–6930 (2019). Copyright 2019, American Chemical Society.

FIG. 4.

Simulating the crystal nucleation at constant chemical potential. Free energy surfaces, obtained by means of metadynamics simulations, relative to the nucleation of NaCl crystals from solution. The two order parameters SO and SH provide a measure of the crystalline order and the hydration number of Na+ ions, respectively. (a) NVT ensemble. (b) Constant chemical potential via the CμMD method discussed in Ref. 101. (c) Same as (b) but at higher supersaturation. Note that the conventional NVT setup fails to identify the basin corresponding to the crystalline phase. Reprinted with permission from Karmakar et al., J. Chem. Theory Comput. 15, 6923–6930 (2019). Copyright 2019, American Chemical Society.

Close modal

A number of different approaches to circumvent this long-standing issue have been proposed.4 As straightforward as it might sound, simply increasing the size of the simulation box might be a viable option. For instance, Zimmermann et al.9 provided practical guidelines as to how big the simulation box should be in order to avoid solute depletion effects when using seeded MD simulations to extract crystal nucleation rates. In particular, the box should be large enough to (i) contain enough solute to form a critical nucleus and (ii) allow, once a critical nucleus has formed, the solution to remain in a supersaturated state. In practice, simulations with fixed numbers of solute particles must be much larger than these lower bounds to be representative of nucleation in a system, which can exchange solute with its surrounding bulk environment to replenish concentration in the vicinity of the nucleus. An alternative is to use analytical corrections to the free energy surfaces (to be used as staring points in order to estimate J) obtained via “closed” simulations (e.g., NVT or NPT ensembles), as exemplified by the seminal work of Agarwal and Peters.17 In 2016, Salvalaglio et al. computed J (relative to the condensation of liquid droplets from vapor as opposed to crystal nucleation) for small systems using infrequent metadynamics86 and applying a bespoke analytical correction.42 Their results have been very recently validated by Bal.88 

More sophisticated methodologies often seek to modify the computational setup itself. Recent advances in this context include the work of Liu et al.,97 which harnesses the so-called string method in collective variables (SMCV98) to work in the osmotic ensemble (Nliquidμcrystal, P, T). This is an intriguing thermodynamic ensemble, which, applied to crystal nucleation, seeks to conserve the number of particles in the liquid phase Nliquid as well as the chemical potential of the crystalline phase μcrystal, in addition to temperature and pressure. However, it has to be said that, to our knowledge, the SMCV method has never been used to compute J directly so far—albeit it can be used to compare the relative nucleation rate of two different crystal polymorphs98 (a very useful feature).

Another approach is to build on multiscale methods, such as the adaptive resolution simulation (AdResS) scheme pioneered by Wang et al.99 In this case, the system is divided into smoothly connected regions where atoms/particles/molecules move from an atomistic to a coarse-grained description. In doing so, the balance of the different degrees of freedom involved with the atomistic and coarse-grained regions allows one to work in what is effectively a grand canonical ensemble (μ, V, T). The latest development of this approach can be found in Ref. 100, but its potential application to the actual calculation of J remains to be explored.

Crystal nucleation rates have been instead recently computed, once more in the case of NaCl, by Karmakar et al.101 using a variant of the constant chemical potential molecular dynamics (CμMD) method developed by Perego et al.102 In this case, the underlying framework is that of well-tempered metadynamics,103 which allows us to obtained free energy surfaces, such as those illustrated in Fig. 4. In stark contrast with the NVT scenario depicted in Fig. 4, the latest modification of the CμMD method correctly identifies the free energy basing corresponding to the emerging crystalline phase, as reported in Figs. 4(b) and 4(c)—where the supersaturation has been purposefully increased so as to highlight the stability of the crystal. It is also worth noting that a modification of the original CμMD approach, which should be able to lower the computational cost of the methodology, has recently been put forward by Chen and Ren.104 

Solute depletion represents one limitation of attempting to approximate an open system (coupled to its surrounding bulk parent phase) with a closed simulation. As well as exchange of particles with its surroundings, a real system can also exchange heat. This is relevant to simulation of nucleation from the melt as well as from solution.

In molecular dynamics simulations, exchange of heat with implicit surroundings is accomplished by augmenting equations of motion with stochastic and dissipative forces or by augmenting the system Hamiltonian with additional degrees of freedom, which accurately mimic the effect of coupling to a heat bath for the purposes of sampling. The resulting microscopic dynamics are however entirely fictitious. It is common to “thermostat” all degrees of freedom in this way.

Consider, for example, a growing nucleus. As the parent phase transforms into a crystal, latent heat is released at the interface. In a standard molecular dynamics simulation, this heat is immediately removed via the thermostat on a (typically) short timescale determined by the thermostat coupling parameter such that the temperature of the system remains both uniform and constant throughout the simulation. Similarly, a shrinking crystal will absorb heat from the parent phase, which is then rapidly replenished via the thermostat.

Thermostatted molecular dynamics simulations of this kind are entirely appropriate to the calculations of the free energy barrier to nucleation within CNT. The ensemble of nuclei which have size n is assumed within CNT to contain both shrinking and growing members such that the net flow of heat between nuclei and surroundings is zero, as is the case for an ensemble of configurations sampled via a simulation in which all degrees of freedom are coupled to a thermostat. This does, however, imply that techniques which infer the free energy gradient from the average dynamics of the nucleus size metric at each n54 must sample over “swarms” of trajectories which contain both growing and shrinking nuclei.

However, for “brute force” or rare event methods, it is not immediately clear that thermostatting all degrees of system is the appropriate choice. Consider again a growing nucleus releasing heat into its immediate surroundings. In a real/open system, this heat must be transported away from the nucleus before reaching the implicit heat bath represented by the surrounding bulk parent phase. This may occur on a timescale slower than the growth rate of the nucleus, meaning that subsequent additions to the nucleus take place at an elevated local kinetic temperature. This entirely physical effect would not be captured via a simulation in which a thermostat is applied globally. This has been discussed in the context of crystal growth at planar ice/water interfaces105 as an explanation for lower growth rates when compared to experiment but not explored quantitatively in simulations of nucleation to our knowledge. This is most likely to be most relevant at strong supercooling where rapid crystal growth is hindered by the limited speed at which heat can be transported away from the interface.

Within CNT, the rate of attachment to a critical nucleus appears in the exponential prefactor. The use of a global thermostat is therefore unlikely to have an impact on nucleation rates calculated from simulation in the case of freezing from the melt where uncertainties typically span several orders of magnitude.

In the case of nucleation from solution, the (additional) analogous problem is the transport of solute to/from the nucleus. Indeed simulating with a constant number of solute particles (as is common) is analogous to having no thermostat at all. A growing nucleus will deplete the surrounding region of solute, creating a region of effective lower supersaturation, which, in turn, will necessitate a larger critical nucleus. The crystal can grow only until the surrounding region is no longer supersaturated with respect to the growing phase and an equilibrium is reached, preventing the critical nucleus size from ever being reached, a much more serious problem. In fact, one could argue that these limitations in terms of transport affect simulations of nucleation in a very similar fashion to what happens when constraining the number of particles in the system (see Sec. III C).

As discussed above, implementing grand canonical molecular dynamics simulations is increasingly feasible but still unusual. With the sole exception of Ref. 102, there has been (to our knowledge) no study of how the location of the insertion/removal region can impact upon growth rates. Ideally, this should be optimized in terms of distance from the nucleus and insertion/removal frequency such that the diffusive kinetics of solute movement into and out of the simulated areas are accurately captured. For example, a solute/solvent force field that accurately captures the solute diffusion coefficient plus insertion sufficiently far from the nucleus such any memory of the insertion position/orientation of the solute is lost before it can be incorporated into the nucleus. As with the precise details of the thermostat implementation, this is unlikely to have a major impact on calculated nucleation rates due to the dominance of the exponential term.

In context of Monte Carlo simulations of nucleation in simple lattice models, there has been some investigation into how microscopic kinetics impacts the dynamics of the order parameter or collective variable used to characterize progress along the nucleation process, typically the cluster size. Kuipers and Barkema106 explored the usual assumption that a collective variable or order parameter (e.g., nucleus size) evolves according to Markovian dynamics. This appears to be a good approximation in the case where insertion or removal of lattice gas particles can occur at any position (including directly at the cluster boundary). However, a choice more appropriate to modeling solute precipitation on a lattice is to use Kawasaki-type local exchange moves107 which capture mass transport to/from the nucleus rather than instantaneous insertion/removal. Here, a simple Markovian model was found to be less than optimal.

Similarly, in Ref. 108, it is shown that the quality of fit to a Markovian model of nucleus size evolution is significantly worse when using diffusive microscopic kinetics instead of allowing particles within or at the boundary of a nucleus to be exchanged directly with a reservoir without including transport. This is clearly very relevant when dealing with seeded simulations. When inferring free energy gradients from nucleus size dynamics by fitting a Markovian model, this may lead to substantial sources of error. In the extreme case, nucleation which is diffusion limited implies no separation of timescales between the dynamics of solute particles and the dynamics of the nucleus size coordinate. In such cases, calculations of nucleation rates based on dynamics of that collective coordinate on a free energy landscape are inappropriate and would lead to potentially meaningless results.

We note here that even though a physical choice of microscopic kinetics may invalidate the assumption that the dynamics of reaction coordinate are Markovian, this does not necessitate that a Markovian model (such as that implied by CNT) is incapable of correctly capturing the nucleation rate, provided the thermodynamic assumptions of CNT are satisfied. In other contexts,109 it has been established that it must be possible to correctly predict rates from a Markovian model of the reaction coordinate dynamics (even if those dynamics are not Markovian), provided the reaction coordinate is chosen to be the committor. See Sec. III F for further discussion on optimization of the reaction coordinate, the committor, and its relationship to the nucleus size.

The known limitations of CNT are many and they have been reviewed extensively elsewhere (see, e.g., Refs. 3, 12, and 13 for some recent perspectives). However, it is important to stress that new issues related to CNT continue to appear as the community pushes the boundary of atomistic simulations to compute J for either more complex systems or employing novel enhanced sampling techniques. Two-step or even multi-step nucleation, a scenario which is not taken into account by the original formulation of CNT, is an excellent example. In fact, in the last few years, both experiments and simulations have reached a point where investigating complex crystallization processes, such as biomineralization and the formation of hydrates, is now a possibility. These systems are famous for their hotly debated mechanism of crystal nucleation and growth, which involves multiple steps not necessarily well described by even extended versions of CNT. Recent computational efforts in the field of methane hydrates that seems to validate at least, in part, the usage of CNT in that context include the work of Arjun and Bolhuis,48 who have used TIS to compute the homogeneous nucleation rate of methane hydrate at relatively low supersaturations. They find that the rate computed via TIS and the rate obtained via CNT (starting from the free energy barrier they have calculated) are, in fact, in good agreement, and they provide a thorough analysis of the uncertainties associated with their estimates as well. Various extensions to the original CNT framework have been put forward through the years in order to deal with multi-step nucleation: the one pioneered by Peters110 in 2011 has been recently used as the starting point to describe the homogeneous nucleation of methane hydrates111 via non-equilibrium MD.

Another issue with the potential of undermining one of the fundamental assumptions of CNT is the emergence of non-Markovian dynamics. Through the lens of CNT, crystal nucleation can be thought as the time evolution of the distribution of crystalline nuclei within the system, which is usually described via the Fokker–Planck equation as a Markovian (i.e., both stochastic and memory-less) process. This assumption has been repeatedly questioned: for instance, Kuipers and Barkema106 studied nucleation in the celebrated Izing model and concluded that when dealing with realistic diffusive dynamics, the Fokker–Plank equation cannot be used and non-Markovian effects should somehow be incorporated instead. The deviation of the nuclei from the spherical shape assumed by CNT (another point of contention) is also blamed as a potential source of uncertainty in the same work. A more recent example is the work of Kuhnhold et al.,14 where the authors argue that CNT can, in fact, be considered as a limiting case of a more general theory that contains memory and out-of-equilibrium effects: broadly speaking, nucleation cannot be considered, according to the authors, as either a Markovian or a diffusive process. Related to the issue of non-Markovian dynamics is the fact that the CNT expression for J refers to a steady state nucleation rate, while in many cases, we observe transient nucleation. This problem has been previously reviewed by, e.g., Sear.112 Here, we mention a recent work by Myint et al.,113 focusing on the crystallization of ice VII at high pressure, where some aspects of hydrodynamic have been incorporated in a theoretical treatment that allows us to relax the assumption of steady-state nucleation.

Very practical considerations about the reliability of CNT often arise when applying the latter to estimate J starting from the basic “ingredients,” that is, the kinetic prefactor and the free energy barrier ΔG(n∗) associated with the critical nucleus size n*. In turn, ΔG(n∗) can be obtained from the free energy difference Δμ between crystal and parent phase and the interfacial free energy γ between the two. The seeded MD approach we have discussed in Sec. II A provides the perfect platform for this discussion, as it usually relies on several estimates of these ingredients, in some cases computed independently and all of them affected by some degrees of uncertainty. An interesting as well as fundamental question concerns the definition of n, that is, the number of molecules within the crystalline nuclei. This deceptively simple question is not easily answered, as it depends on the choice of a specific order parameter—as we shall discuss in detail in Sec. III F. In a recent work, Cheng et al.114 argued that the usage of the solid–liquid Gibbs dividing surface is the one choice consistent with the CNT formalism and put forward a methodology, based on combining simulations of planar interfaces and three-dimensional nuclei, to alleviate the uncertainty associated with the estimate of ΔG(n∗) and, thus, of J. Another investigation that shed new light onto the error propagation involved with the calculation of J in the context of seeded MD is the work of Lifanov et al.108 on nucleation in a multi-species lattice model. In addition to some consideration of non-Markovian dynamics, the authors show that even when taking great care in computing Δμ, seeded MD cannot provide accurate estimates of ΔG(n∗) in the low supersaturation regime. To provide some context, the fact that ΔG(n∗) enters as an exponential in the expression of J implies that a 43% error with respect to ΔG(n∗) can lead to an uncertainty of some ten orders of magnitudes in terms of the nucleation rate.

The last aspect we believe it is worth discussing in the context of CNT limitations is the assumption that the nuclei are supposed to be perfectly spherical for every value of n. Note that while this is a very common assumption, it is perfectly possible to rewrite the relevant CNT expressions assuming nuclei shapes other than spherical—albeit the situation becomes complicated if one wants to allow for the occurrence of any arbitrary shape. Experimental measurements, most prominently on colloidal systems, have demonstrated that crystalline nuclei (particularly around the critical size and/or within the early stages of the crystallization process) are rarely spherical (see, e.g., Refs. 24 and 115). Simulations have been key to show that deviations from the spherical shape are commonly encountered across a much wider range of systems, particularly during the early stages of the nucleation process as well as strong supercooling, where the critical nuclei are smaller and thus affected to a greater extent by potential deviation from the spherical assumption.116,117 This is relevant for seeded simulations of crystal nucleation as well, as the crystalline seeds are usually built as spherical particles.59 Several modifications to CNT have been proposed throughout the years to account for the deviation of the nuclei from a spherical shape—see, e.g., the work of Prestipino et al.118,119 A more recent example is provided by the work of Lutsko,120 which entirely removes the need for the assumption about the spherical shape in the first place. We also note here that in the context of nucleation in the Ising model, well-established corrections for fluctuations in nucleus shape as well as effective “surface energies” (a concept we have already encountered in Sec. II B) for small nuclei have been explicitly tested against FFS and umbrella sampling calculations and found to be highly accurate.121 

Every route toward the estimate of J involves the choice of an order parameter, which is a mathematical object that allows us to identify which atoms, molecules, or particles belong to the crystalline nuclei. In this work, we are going to adopt the nomenclature of Peters,122 according to which a collective variable is any function of the full phase space coordinates, an order parameter is a collective variable that distinguishes the typical reactant and product configurations, and the reaction coordinate (RC hereafter) is a special order parameter that accurately quantifies dynamical progress from reactant to product. It is important to stress that when dealing with complex processes, such as crystal nucleation, the “true” RC might very well be a high-dimensional set of variables that we have no hope to identify correctly in their entirety. What we can do instead is to coarse-grain the manifold of said variables into a handful of tractable order parameters that we treat as our RC.

Many order parameters to identify crystalline nuclei build on the famous work of Steinhardt et al.123 and are briefly reviewed in Ref. 124. Topological order parameters, such as the permutation invariant vector (PIV, reviewed in Ref. 125), represent a valid alternative, and of course, machine learning has found its way into this particular facet of the field as well, not only to craft new order parameters (see, e.g., Refs. 126 and 127) but also to mine the extensive portfolio of the existing ones in the attempt to pinpoint the combination yielding the best accuracy.128 In a similar fashion, there have been several attempts to automate the choice of the RC in the last few years, albeit none—to our knowledge—applied to crystal nucleation just yet. As an example, Krivov129 recently suggested an adaptive optimization of a multivalued RC to be used in the context of protein folding, a long-standing problem involving an incredibly high-dimensional space. Another approach is that of the spectral gap optimization of order parameters (SGOOP) put forward by Smith et al.130 and applied to the study of the kinetics of chemical reactions. It is important to note that these examples are still very much relevant to our discussion because it is now clear that, in many cases, we need to step away from the assumption that a single RC, such as the size n of the largest crystalline nucleus within the system, would be sufficient to provide a reliable estimate of J.

In fact, crystal nucleation has been traditionally thought of as a rather one-dimensional problem, partially because n is the one and only RC involved in the CNT framework. However, we now have evidence that this is not the case. Concerning the nucleation of crystal from supersaturated solutions, it is understood that, at the very least, one needs to work in a two-dimensional RC space involving the crystallinity as well as the local density (see, e.g., Ref. 131). When it comes to nucleation from the supercooled liquid phase, Jungblut and Dellago132 showed in the case of a LJ system that the choice of n as the only RC can lead to substantial memory effects, which invalidate the assumption of nucleation as a Markovian process and result into vastly different nucleation rates (which they evaluate via TIS). Interestingly, these effects are connected to relative abundance of the two polymorphs, bcc and fcc, within the crystalline seed that have been used to initiate the nucleation events. Indeed, polymorphism is a very common occurrence that is very difficult to assess a priori, as even when the most stable crystal form is known, there is no guarantee that the nucleation process would not involve intermediates characterized by different crystalline structures—again, the LJ system provides a striking example in this sense (see e.g., Ref. 133), or even result in a long-lived metastable bulk phase.

An additional layer of complexity has been recently brought to light by Liang et al.,134 who have investigated the crystal nucleation molybdenum via TIS simulations. The result of the committor probability P ( p B | n s * ) they have computed for this system according to the size of the crystalline nuclei alone, n s * , is reported in Fig. 5(a). Note that the committor distribution fails to yield a well-defined narrow peak for P ( p B | n s * ) = 0.5, thus demonstrating that this order parameter does not provide an accurate enough approximation of the reaction coordinate. In fact, not even introducing explicitly an additional degree of freedom in the form of a crystallinity parameter ( Q 6 c l * , inspired by the work of Refs. 135 and 136), is sufficient to improve the committor in this case. The authors had to resort to a third order parameter, namely, a measure of the short-range order (SRO2) in order to univocally determine which nuclei are more or less likely to either dissolve into the liquid or proceed toward post-critical sizes. This is illustrated in Fig. 5(b), where the committor probability P ( p B | n s * , Q 6 c l * , SRO 2 ) is narrowly peaked around zero and one for nuclei characterized by low and high values of SRO2, respectively, thus indicating that only sufficiently compact nuclei are likely to cross the free energy barrier associated with the nucleation process.

FIG. 5.

Investigating the suitability of order parameters for crystal nucleation. (a) Committor distribution P ( p B | n s * ) , where n s * corresponds to the number of atoms within the largest crystalline nucleus. (b) Committor distributions P ( p B | n s * , Q 6 c l * , SRO 2 ) for two sets of configurations characterized by the same n s * , the same degree of crystallinity (quantified by the order parameter Q 6 c l * , inspired by the work of Refs. 135 and 136), but different short-range order SRO2. Reprinted from Liang et al., J. Chem. Phys. 152, 224504 (2020) with the permission of AIP Publishing LLC.

FIG. 5.

Investigating the suitability of order parameters for crystal nucleation. (a) Committor distribution P ( p B | n s * ) , where n s * corresponds to the number of atoms within the largest crystalline nucleus. (b) Committor distributions P ( p B | n s * , Q 6 c l * , SRO 2 ) for two sets of configurations characterized by the same n s * , the same degree of crystallinity (quantified by the order parameter Q 6 c l * , inspired by the work of Refs. 135 and 136), but different short-range order SRO2. Reprinted from Liang et al., J. Chem. Phys. 152, 224504 (2020) with the permission of AIP Publishing LLC.

Close modal

The freezing of water to ice offers an especially challenging testing ground when it comes to the choice of the reaction coordinate. For instance, the work of Prestipino,137 focusing on the thermodynamic of crystal nucleation for a monoatomic water model, shows that the quality of a reaction coordinate cannot be assessed simply on the basis of the free energy barrier height obtained. On the contrary, different choices of reaction coordinates led to comparable free energy barriers but significantly different fractions of hexagonal ice, cubic ice, and ice-0 as well within the post-critical nuclei. Interestingly, in this particular case, the committor probabilities for the different reaction coordinates used are quite similar to each other, which implies that none of them alone captures the nucleation process accurately enough. In fact, it has becoming apparent that in some cases, we really need to combine multiple order parameters to obtain a sufficiently accurate reaction coordinate. As an example, Niu et al.39 recently computed the ice nucleation rate for the TIP4P/ice model by bringing together a local order parameter with a linear combination of seven long-range descriptors. Hence, the relevance of the above-mentioned approaches aimed at learning, automating, and selecting multi-dimensional degrees of freedom, as a potential avenue to lessen the burden associated with the often employed trial-and-error approach when building the reaction coordinate.

Quantifying the uncertainty associated with the choice of a particular reaction coordinate is also key. In principle, the use of reactive flux calculations (e.g., TIS and FFS) compensates exactly for potential inaccuracy in terms of the reaction coordinate,122 while indirect estimates of J obtained via combining free energy profiles with path-sampling methodologies can be corrected, under certain assumptions, via, e.g., a committor analysis. The latter is a popular, albeit computationally extensive, option, and specific flavors are available to identify and mitigate the inadequacy of the reaction coordinate.122,138 Systems characterized by an inherently slow dynamics, which we will discuss in greater detail in Sec. III G, are especially challenging to deal with.109,139–141

An interesting aspect of relevance for the computational community is the availability of the implementation of the many order parameters we use to construct reaction coordinates—not just as standalone mathematical objects but also in conjunction with MD and MC packages as well so that they can be leveraged to deploy the enhanced sampling method of choice. In this context, several highly collaborative projects have been born with the aim of facilitating the access to a variety of order parameters for enhanced sampling simulations, such as the SSAGES142 suite or the PLUMED consortium the authors are part of Ref. 23.

Finally, we discuss an assumption implicit in many path-sampling techniques applied to the calculation of nucleation rates. Forward Flux Sampling (FFS), in particular, has become the principle workhorse of many rare event studies.10,143 Alternative techniques such as transition interface sampling (TIS) are also used in the context of nucleation44 but generally acknowledged to be somewhat more complex to implement. The aimless shooting approach to path sampling is a third example,144 used in a number of nucleation studies.122 The relationship between these methods has been very recently summarized in an informative review by Bolhuis and Swenson.45 

Both FFS and TIS make the implicit assumption that the interval between nucleation events is dominant in determining the nucleation rate. The time taken for the event to occur is neglected. A further practical consideration is that the time taken for the nucleation event to occur must be tractably accessible to simulation. For example, the path-sampling stage of an FFS calculation requires contiguous sampling of trajectories which return to the boundary of the quasi-stationary parent state from an initial configuration containing a critical nucleus. TIS requires sampling of contiguous trajectories which leave the parent state, form a critical nucleus, and then return to the boundary of the parent state. Many such trajectories must be sampled, in practice, limiting the application of these techniques to fast (but still rare) nucleation processes which occur on nano-second timescales or less.

These two limitations have the effect of limiting the window of supercooling/supersaturation where FFS/TIS calculations can be accurately applied. In freezing from the melt, as supercooling increases, the interval between nucleation events becomes smaller, while dynamics become slower. Eventually the assumption of fast nucleation breaks down, as does our ability to sample appropriate trajectories. Applying these methods in regimes where the dynamics of nucleation are slow is likely to be fruitless.

One approach to obtaining nucleation rates in this regime is to rely on CNT. The method of Knott et al.54 (discussed earlier) pioneered the approach of fitting a free energy barrier to gradients inferred from mean growth/shrinkage rates over swarms of short trajectories, precisely to circumvent the need to simulate long trajectories. Approaches which do not rely on CNT are also available. The complete nucleation kinetics are subsequently reconstructed from a Markov state model which uses these crossing rates. This relies on a loss of memory between non-adjacent states, i.e., that the dynamics of the reaction coordinate are Markovian on at least the broad scale corresponding to crossings between windows. This assumption may lead to correct rates even if the true dynamics are non-Markovian, providing a suitably optimal choice of reaction coordinate is made, i.e., the reaction coordinate is a good approximation of the committor.109 

A similar method is milestoning,145 which can be used to inform both Markovian and non-Markovian models for the dynamics of the RC from ensembles of trajectories which span part of the nucleation process. In the context of nucleation, this has been combined with MFEP methods to inform the placement of “milestones” and used to compute rates of crystallization from the melt.146 However, simulations of nucleation in the slow dynamics regime remain relatively unusual compared to studies in the regime where the process is assumed rapid compared to the waiting time between nucleation events.

Molecular simulations have been playing an important role in complementing the experimental insight into the kinetics of crystal nucleation for decades. Recent advancements in the field of enhanced sampling, together with the ever-increasing accuracy of the available force fields as well as the unprecedented capabilities of the current high performance computing facilities, all provide exciting new avenues for the computational scientist. At a time when computing crystal nucleation rates is becoming a possibility for many research groups across the globe, we believe that it is important to take stock of what we have learned in terms of the limitations of our simulations, particularly in terms of the several aspects that can contribute toward the uncertainty of our results. In this Perspective, we focused on atomistic simulations and discussed seven of such aspects in the context of the recent literature.

We started by questioning whether the nucleation events we need to investigate in order to extract the nucleation rate can, in fact, be considered as “rare enough” to justify the assumption of a steady-state nucleation rate. This is key to quantify the uncertainty associated with our estimate of J via both unbiased and enhanced sampling simulations, and we have highlighted a few methodologies that are now available to do just that. Then, we moved onto finite size effects, a long-standing issue that affects atomistic simulations as a whole. These are exciting times for the community, as we are now able to probe the emergence of these spurious effect in models containing up to several million atoms. Related to this aspect is the problem of solute depletion, another practical aspect of molecular simulations that we can now overcome by leveraging a diverse portfolio of computational techniques.

Only within the last few years, the community has started to assess the effect of thermostatting the entire simulation box while simulating crystal nucleation—a choice that might play a role in determining whether the nucleation process can be considered as Markovian or not. Indeed, the emergence of non-Markovian RC dynamics is a problematic aspect that clashes with some of the assumptions CNT builds upon—and that we have discussed in this Perspective together with some perhaps not entirely obvious limitations of this theory.

It is impossible not to touch on the choice of the reaction coordinate when it comes to computing crystal nucleation rates. Much has been said in the past, but here we focused on the growing evidence that in many cases, we need to acknowledge the complexity of the nucleation process and construct very specific reaction coordinates starting from a diverse portfolio of order parameters. This is where, we believe, the rise of machine learning can contribute to further the current state of affairs, as a prime tool to identify, select, and potentially even learn the most accurate reaction coordinate available to us. We also argue, though, that at this point in time, we are running out of excuses for avoiding to validate our reaction coordinates, as we now have both the theoretical tools and—in most cases—the computational power as well.

Finally, we offer a perhaps provocative take on the usage of path-sampling techniques, which are righteously becoming more and more popular as the tool of the trade to compute crystal nucleation rates but that still suffer from a number of limitations.

At this stage, the reader might think that our perspective is painting a rather bleak picture: substantial inconsistencies between experimental and simulated nucleation rate persists; investigating realistic systems at conditions relevant to the experimental reality is, more often than not, still computationally prohibitive; and, even by using state-of-the-art methodologies, our estimates of J are likely to be affected by a sizeable degree of uncertainty, accumulated through a number of diverse aspects, only seven of which we have discussed in some details here.

However, we believe that the field is in a very different, in fact, much better spot than it was only a decade ago: this is because the importance of quantifying uncertainty has now been acknowledged across the board, with several research groups going back to the basics and questioning the assumptions we have built our results upon. Of course, this is an often painful process that is bound to uncover some more problematic aspects related to our simulations, but we argue that we should embrace the challenge and work toward promoting the transparency as well as the reproducibility of our work.

We also feel that the gap between the community working on the so-called model systems and those computational scientists interested instead in more realistic system is quickly narrowing, which, in turn, should facilitate knowledge exchange in terms of both techniques and expertise. In addition, it is worth pointing out that many of the issues we have addressed in here are by no means exclusive to crystal nucleation. The occurrence of rare events is ubiquitous in the natural sciences,147 and quantifying the kinetics characterizing these processes is a fundamental open question with reverberations across, e.g., climate science148,149 as well as epidemiology.150,151

In summary, we hope we have managed with this Perspective not only to highlight some problematic aspects but also to showcase some of the excellent work that has recently been done to overcome these challenges and break new ground. The field is moving forward at a very fast pace, leveraging concepts righteously taken from other disciplines—it is our duty to keep an eye on the accuracy of our simulations as we embark in the next chapter of this long-standing endeavor.

We gratefully acknowledge the use of Athena at HPC Midlands+, which was funded by the EPSRC through Grant No. EP/P020232/1, through the HPC Midlands+ consortium, as well as the high-performance computing facilities provided by the Scientific Computing Research Technology Platform at the University of Warwick. We would also like to acknowledge the EPSRC program (Grant No. EP/R018820/1; “Crystallisation in the Real World: Delivering Control through Theory and Experiment”) as well as the EPSRC Centre for Doctoral Training in Modeling of Heterogeneous Systems (EPSRC Grant No. EP/S022848/1). We gratefully acknowledge Dr. Matteo Salvalaglio for insightful discussions and for reading an earlier version of the manuscript.

The data we have used to put together Fig. 2 are available from the corresponding author upon reasonable request.

Details about the LJ calculations reported in Fig. 2. Simulations were performed in the NPT ensemble using the LAMMPS MD package (version released June 5, 2019). The Lennard-Jones parameter was implemented with a cutoff of 3.5 σ and an associated tail correction. A periodic cubic box of 4000 liquid atoms was produced by melting a face centered cubic crystal. A small amount of linear quench was then applied to lower the system to the desired temperature before the simulation began. An isotropic pressure of p* = 5.68 was applied through the simulation using a chain of 5 thermostats attached to a Hoover barostat of damping parameter 0.5t* with an applied Martyna–Tobias–Klein correction.

Thermostatting was performed using a chain of five Nosé–Hoover thermostats of damping coefficient 0.1t*. The size of the cluster was calculated using the ten Wolde q6 order parameter135 as implemented in the molecular dynamics package LAMMPS [we are aware of the issue with q6(n) in LAMMPS, but errors from this source should be small by the time the critical size is reached]. The system was observed over 5 000 000 time steps of length t* = 0.002, with q6(n) being recorded every 100 time steps.

Simulations (258 for T* = 0.68; 224 for T* = 0.765; and 823 for T* = 0.8) were analyzed in Python by fitting a sigmoid curve to the value of q6(n). If this never exceeded 2000, the simulation was determined not to have crystallized. Shifted curves were created by removing the time associated with the midpoint of the sigmoid curve from the time steps of the simulations.

The survival probabilities were found by dividing the output of simulations into five approximately equal-sized sets. The location of the midpoints of the sigmoid curves were binned into 15 bins ranging from t* = 0 to the maximum value of t* at which a midpoint was recorded. The plotted values are the mean of the values for these five sets, with the error bar corresponding to the error in the mean (standard deviation divided by the square root of the number of sets). Exponentials of the form exp J t * α , with α either allowed to vary or constrained to 1, were fitted using scipy.optimize.curve_fit. Bounds were imposed to ensure that a reasonable fit was given, with details of these bounds available in the analysis code. For T* = 0.68, all errors on points were 0, so errors were neglected. For T* = 0.8, only the error on the initial point was 0, and for T* = 0.765, errors on the first and last point were 0; in these cases, fitting was performed on data points with associated errors, and data points with 0 error were omitted.

Analysis and input scripts can be found on GitHub (https://github.com/keb721/7DeadlySins).

1.
J. J.
De Yoreo
and
P. G.
Vekilov
, “
Principles of crystal nucleation and growth
,”
Rev. Mineral. Geochem.
54
,
57
93
(
2003
).
2.
P.
Yi
and
G. C.
Rutledge
, “
Molecular origins of homogeneous crystal nucleation
,”
Annu. Rev. Chem. Biomol. Eng.
3
,
157
182
(
2012
).
3.
S.
Karthika
,
T. K.
Radhakrishnan
, and
P.
Kalaichelvi
, “
A review of classical and nonclassical nucleation theories
,”
Cryst. Growth Des.
16
,
6663
6681
(
2016
).
4.
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
).
5.
T.
Bartels-Rausch
, “
Chemistry: Ten things we need to know about ice and snow
,”
Nature
494
,
27
29
(
2013
).
6.
S.
Weiner
and
L.
Addadi
, “
Crystallization pathways in biomineralization
,”
Annu. Rev. Mater. Res.
41
,
21
40
(
2011
).
7.
A. Y.
Lee
,
D.
Erdemir
, and
A. S.
Myerson
, “
Crystal polymorphism in chemical process development
,”
Annu. Rev. Chem. Biomol. Eng.
2
,
259
280
(
2011
).
8.
G. C.
Sosso
,
G.
Miceli
,
S.
Caravati
,
F.
Giberti
,
J.
Behler
, and
M.
Bernasconi
, “
Fast crystallization of the phase change compound GeTe by large-scale molecular dynamics simulations
,”
J. Phys. Chem. Lett.
4
,
4241
4246
(
2013
).
9.
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
).
10.
A.
Haji-Akbari
and
P. G.
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
).
11.
R.
Demichelis
,
A.
Schuitemaker
,
N. A.
Garcia
,
K. B.
Koziara
,
M.
De La Pierre
,
P.
Raiteri
, and
J. D.
Gale
, “
Simulation of crystallization of biominerals
,”
Annu. Rev. Mater. Res.
48
,
327
352
(
2018
).
12.
P. J. M.
Smeets
,
A. R.
Finney
,
W. J. E. M.
Habraken
,
F.
Nudelman
,
H.
Friedrich
,
J.
Laven
,
J. J.
De Yoreo
,
P. M.
Rodger
, and
N. A. J. M.
Sommerdijk
, “
A classical view on nonclassical nucleation
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
E7882
E7890
(
2017
).
13.
P. G.
Vekilov
, “
Nonclassical nucleation
,” in
Crystallization via Nonclassical Pathways Volume 1: Nucleation, Assembly, Observation & Application
, ACS Symposium Series Vol. 1358 (
American Chemical Society
,
2020
), Sec. 2, pp.
19
46
.
14.
A.
Kuhnhold
,
H.
Meyer
,
G.
Amati
,
P.
Pelagejcev
, and
T.
Schilling
, “
Derivation of an exact, nonequilibrium framework for nucleation: Nucleation is a priori neither diffusive nor Markovian
,”
Phys. Rev. E
100
,
052140
(
2019
).
15.
S.
Hussain
and
A.
Haji-Akbari
, “
How to quantify and avoid finite size effects in computational studies of crystal nucleation: The case of heterogeneous ice nucleation
,”
J. Chem. Phys.
154
,
014108
(
2021
).
16.
G.
Fiorucci
,
G. M.
Coli
,
J. T.
Padding
, and
M.
Dijkstra
, “
The effect of hydrodynamics on the crystal nucleation of nearly hard spheres
,”
J. Chem. Phys.
152
,
064903
(
2020
).
17.
V.
Agarwal
and
B.
Peters
, “
Nucleation near the eutectic point in a Potts-lattice gas model
,”
J. Chem. Phys.
140
,
084111
(
2014
).
18.
D.
Gebauer
,
P.
Raiteri
,
J. D.
Gale
, and
H.
Cölfen
, “
On classical and non-classical views on nucleation
,”
Am. J. Sci.
318
,
969
988
(
2018
).
19.
J.
Behler
, “
Perspective: Machine learning potentials for atomistic simulations
,”
J. Chem. Phys.
145
,
170901
(
2016
).
20.
G. C.
Sosso
and
M.
Bernasconi
, “
Harnessing machine learning potentials to understand the functional properties of phase-change materials
,”
MRS Bull.
44
,
705
709
(
2019
).
21.
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
).
22.
G. C.
Sosso
,
T. F.
Whale
,
M. A.
Holden
,
P.
Pedevilla
,
B. J.
Murray
, and
A.
Michaelides
, “
Unravelling the origins of ice nucleation on organic crystals
,”
Chem. Sci.
9
,
8077
8088
(
2018
).
23.
M.
Bonomi
,
G.
Bussi
,
C.
Camilloni
,
G. A.
Tribello
,
P.
Banáš
,
A.
Barducci
,
M.
Bernetti
,
P. G.
Bolhuis
,
S.
Bottaro
,
D.
Branduardi
,
R.
Capelli
,
P.
Carloni
,
M.
Ceriotti
,
A.
Cesari
,
H.
Chen
,
W.
Chen
,
F.
Colizzi
,
S.
De
,
M.
De La Pierre
,
D.
Donadio
,
V.
Drobot
,
B.
Ensing
,
A. L.
Ferguson
,
M.
Filizola
,
J. S.
Fraser
,
H.
Fu
,
P.
Gasparotto
,
F. L.
Gervasio
,
F.
Giberti
,
A.
Gil-Ley
,
T.
Giorgino
,
G. T.
Heller
,
G. M.
Hocky
,
M.
Iannuzzi
,
M.
Invernizzi
,
K. E.
Jelfs
,
A.
Jussupow
,
E.
Kirilin
,
A.
Laio
,
V.
Limongelli
,
K.
Lindorff-Larsen
,
T.
Löhr
,
F.
Marinelli
,
L.
Martin-Samos
,
M.
Masetti
,
R.
Meyer
,
A.
Michaelides
,
C.
Molteni
,
T.
Morishita
,
M.
Nava
,
C.
Paissoni
,
E.
Papaleo
,
M.
Parrinello
,
J.
Pfaendtner
,
P.
Piaggi
,
G.
Piccini
,
A.
Pietropaolo
,
F.
Pietrucci
,
S.
Pipolo
,
D.
Provasi
,
D.
Quigley
,
P.
Raiteri
,
S.
Raniolo
,
J.
Rydzewski
,
M.
Salvalaglio
,
G. C.
Sosso
,
V.
Spiwok
,
J.
Šponer
,
D. W. H.
Swenson
,
P.
Tiwary
,
O.
Valsson
,
M.
Vendruscolo
,
G. A.
Voth
,
A.
White
, and
PLUMED consortium
, “
Promoting transparency and reproducibility in enhanced molecular simulations
,”
Nat. Methods
16
,
670
673
(
2019
).
24.
Z.
Wang
,
F.
Wang
,
Y.
Peng
, and
Y.
Han
, “
Direct observation of liquid nucleus growth in homogeneous melting of colloidal crystals
,”
Nat. Commun.
6
(
1
),
6942
(
2015
).
25.
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
).
26.
D. A.
Nicholson
and
G. C.
Rutledge
, “
Analysis of nucleation using mean first-passage time data from molecular dynamics simulation
,”
J. Chem. Phys.
144
,
134105
(
2016
).
27.
V. G.
Baidakov
and
K. R.
Protsenko
, “
Spontaneous crystallization of a supercooled Lennard-Jones liquid: Molecular dynamics simulation
,”
J. Phys. Chem. B
123
,
8103
8112
(
2019
).
28.
W.
Ouyang
,
B.
Sun
,
Z.
Sun
, and
S.
Xu
, “
Entire crystallization process of Lennard-Jones liquids: A large-scale molecular dynamics study
,”
J. Chem. Phys.
152
,
054903
(
2020
).
29.
S.
An
,
Y.
Li
,
J.
Li
,
S.
Zhao
,
B.
Liu
, and
P.
Guan
, “
The linear relationship between diffusivity and crystallization kinetics in a deeply supercooled liquid Ni50Ti50 alloy
,”
Acta Mater.
152
,
1
6
(
2018
).
30.
A.
Mahata
,
M. A.
Zaeem
, and
M. I.
Baskes
, “
Understanding homogeneous nucleation in solidification of aluminum by molecular dynamics simulations
,”
Modell. Simul. Mater. Sci. Eng.
26
,
025007
(
2018
).
31.
Y.
Shibuta
,
S.
Sakane
,
E.
Miyoshi
,
S.
Okita
,
T.
Takaki
, and
M.
Ohno
, “
Heterogeneity in homogeneous nucleation from billion-atom molecular dynamics simulation of solidification of pure metal
,”
Nat. Commun.
8
(
1
),
10
(
2017
).
32.
J.-M.
Leyssale
,
J.
Delhommelle
, and
C.
Millot
, “
A molecular dynamics study of homogeneous crystal nucleation in liquid nitrogen
,”
Chem. Phys. Lett.
375
,
612
618
(
2003
).
33.
S.
Punnathanam
and
P. A.
Monson
, “
Crystal nucleation in binary hard sphere mixtures: A Monte Carlo simulation study
,”
J. Chem. Phys.
125
,
024508
(
2006
).
34.
P.
Yi
and
G. C.
Rutledge
, “
Molecular simulation of crystal nucleation in n-octane melts
,”
J. Chem. Phys.
131
,
134902
(
2009
).
35.
L.
Filion
,
M.
Hermes
,
R.
Ni
, and
M.
Dijkstra
, “
Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: A comparison of simulation techniques
,”
J. Chem. Phys.
133
,
244115
(
2010
).
36.
D.
Quigley
and
P. M.
Rodger
, “
Metadynamics simulations of ice nucleation and growth
,”
J. Chem. Phys.
128
,
154518
(
2008
).
37.
D.
Quigley
and
P. M.
Rodger
, “
A metadynamics-based approach to sampling crystallisation events
,”
Mol. Simul.
35
,
613
623
(
2009
).
38.
F.
Giberti
,
M.
Salvalaglio
, and
M.
Parrinello
, “
Metadynamics studies of crystal nucleation
,”
IUCrJ
2
(
2
),
256
266
(
2015
).
39.
H.
Niu
,
Y. I.
Yang
, and
M.
Parrinello
, “
Temperature dependence of homogeneous nucleation in ice
,”
Phys. Rev. Lett.
122
,
245501
(
2019
).
40.
G.
Menzl
,
A.
Singraber
, and
C.
Dellago
, “
S-shooting: A Bennett–Chandler-like method for the computation of rate constants from committor trajectories
,”
Faraday Discuss.
195
,
345
364
(
2016
).
41.
M.
Salvalaglio
,
P.
Tiwary
, and
M.
Parrinello
, “
Assessing the reliability of the dynamics reconstructed from metadynamics
,”
J. Chem. Theory Comput.
10
,
1420
1425
(
2014
).
42.
M.
Salvalaglio
,
P.
Tiwary
,
G. M.
Maggioni
,
M.
Mazzotti
, and
M.
Parrinello
, “
Overcoming time scale and finite size limitations to compute nucleation rates from small scale well tempered metadynamics simulations
,”
J. Chem. Phys.
145
,
211925
(
2016
).
43.
D.
Moroni
,
P. R.
ten Wolde
, and
P. G.
Bolhuis
, “
Interplay between structure and size in a critical crystal nucleus
,”
Phys. Rev. Lett.
94
,
235703
(
2005
).
44.
A.
Arjun
and
P. G.
Bolhuis
, “
Molecular understanding of homogeneous nucleation of CO2 hydrates using transition path sampling
,”
J. Phys. Chem. B
125
,
338
349
(
2021
).
45.
P. G.
Bolhuis
and
D. W. H.
Swenson
, “
Transition path sampling as Markov chain Monte Carlo of trajectories: Recent algorithms, software, applications, and future outlook
,”
Adv. Theory Simul.
4
,
2000237
(
2021
).
46.
T. S.
van Erp
and
P. G.
Bolhuis
, “
Elaborating transition interface sampling methods
,”
J. Comput. Phys.
205
,
157
181
(
2005
).
47.
C.
Valeriani
,
E.
Sanz
, and
D.
Frenkel
, “
Rate of homogeneous crystal nucleation in molten NaCl
,”
J. Chem. Phys.
122
,
194501
(
2005
).
48.
A.
Arjun
and
P. G.
Bolhuis
, “
Rate prediction for homogeneous nucleation of methane hydrate at moderate supersaturation using transition interface sampling
,”
J. Phys. Chem. B
124
,
8099
8109
(
2020
).
49.
J. T.
Berryman
and
T.
Schilling
, “
Sampling rare events in nonequilibrium and nonstationary systems
,”
J. Chem. Phys.
133
,
244101
(
2010
).
50.
N. B.
Becker
,
R. J.
Allen
, and
P. R.
ten Wolde
, “
Non-stationary forward flux sampling
,”
J. Chem. Phys.
136
,
174118
(
2012
).
51.
C.
Desgranges
and
J.
Delhommelle
, “
Crystal nucleation along an entropic pathway: Teaching liquids how to transition
,”
Phys. Rev. E
98
,
063307
(
2018
).
52.
L.
Bonati
and
M.
Parrinello
, “
Silicon liquid structure and crystal nucleation from ab initio deep metadynamics
,”
Phys. Rev. Lett.
121
,
265701
(
2018
).
53.
C. S.
Adorf
,
T. C.
Moore
,
Y. J. U.
Melle
, and
S. C.
Glotzer
, “
Analysis of self-assembly pathways with unsupervised machine learning algorithms
,”
J. Phys. Chem. B
124
,
69
78
(
2020
).
54.
B. C.
Knott
,
V.
Molinero
,
M. F.
Doherty
, and
B.
Peters
, “
Homogeneous nucleation of methane hydrates: Unrealistic under realistic conditions
,”
J. Am. Chem. Soc.
134
,
19544
19547
(
2012
).
55.
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
).
56.
Y.
Sun
,
H.
Song
,
F.
Zhang
,
L.
Yang
,
Z.
Ye
,
M. I.
Mendelev
,
C.-Z.
Wang
, and
K.-M.
Ho
, “
Overcoming the time limitation in molecular dynamics simulation of crystal nucleation: A persistent-embryo approach
,”
Phys. Rev. Lett.
120
,
085703
(
2018
).
57.
A.
Goswami
,
I. S.
Dalal
, and
J. K.
Singh
, “
Seeding method for ice nucleation under shear
,”
J. Chem. Phys.
153
,
094502
(
2020
).
58.
C.
Wang
,
H.
Wang
,
J.
Wu
, and
Z.
Zhang
, “
Classical nucleation theory of ice nucleation: Second-order correction of thermodynamic parameters
,”
J. Chem. Phys.
154
,
234503
(
2021
).
59.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
, “
Seeding approach to crystal nucleation
,”
J. Chem. Phys.
144
,
034501
(
2016
).
60.
A.
Reinhardt
and
J. P. K.
Doye
, “
Free energy landscapes for homogeneous nucleation of ice for a monatomic water model
,”
J. Chem. Phys.
136
,
054501
(
2012
).
61.
D. T.
Limmer
and
D.
Chandler
, “
Phase diagram of supercooled water confined to hydrophilic nanopores
,”
J. Chem. Phys.
137
,
044509
(
2012
).
62.
J. R.
Espinosa
,
C.
Vega
, and
E.
Sanz
, “
Ice–water interfacial free energy for the TIP4P, TIP4P/2005, TIP4P/ice, and mW models as obtained from the mold integration technique
,”
J. Phys. Chem. C
120
,
8068
8075
(
2016
).
63.
M.
Ambler
,
B.
Vorselaars
,
M. P.
Allen
, and
D.
Quigley
, “
Solid–liquid interfacial free energy of ice Ih, ice Ic, and ice 0 within a mono-atomic model of water via the capillary wave method
,”
J. Chem. Phys.
146
,
074701
(
2017
).
64.
M. N.
Joswiak
,
N.
Duff
,
M. F.
Doherty
, and
B.
Peters
, “
Size-dependent surface free energy and tolman-corrected droplet nucleation of TIP4P/2005 water
,”
J. Phys. Chem. Lett.
4
,
4267
4272
(
2013
).
65.
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
).
66.
F.
Giberti
,
G. A.
Tribello
, and
M.
Parrinello
, “
Transient polymorphism in NaCl
,”
J. Chem. Theory Comput.
9
,
2526
2530
(
2013
).
67.
F.
Moučka
,
M.
Lísal
,
J.
Škvor
,
J.
Jirsák
,
I.
Nezbeda
, and
W. R.
Smith
, “
Molecular simulation of aqueous electrolyte solubility. 2. Osmotic ensemble Monte Carlo methodology for free energy and solubility calculations and application to NaCl
,”
J. Phys. Chem. B
115
,
7849
7861
(
2011
).
68.
J. L.
Aragones
,
E.
Sanz
, and
C.
Vega
, “
Solubility of NaCl in water by molecular simulation revisited
,”
J. Chem. Phys.
136
,
244508
(
2012
).
69.
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
).
70.
A. L.
Benavides
,
J. L.
Aragones
, and
C.
Vega
, “
Consensus on the solubility of NaCl in water from computer simulations using the chemical potential route
,”
J. Chem. Phys.
144
,
124504
(
2016
).
71.
S.
Boothroyd
,
A.
Kerridge
,
A.
Broo
,
D.
Buttar
, and
J.
Anwar
, “
Solubility prediction from first principles: A density of states approach
,”
Phys. Chem. Chem. Phys.
20
,
20981
20987
(
2018
).
72.
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
).
73.
Arjun
,
T. A.
Berendsen
, and
P. G.
Bolhuis
, “
Unbiased atomistic insight in the competing nucleation mechanisms of methane hydrates
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
19305
19310
(
2019
).
74.
R.
Shi
and
H.
Tanaka
, “
Homogeneous nucleation of ferroelectric ice crystal driven by spontaneous dipolar ordering in supercooled TIP5P water
,”
J. Chem. Phys.
151
,
024501
(
2019
).
75.
M.
Fitzner
,
L.
Joly
,
M.
Ma
,
G. C.
Sosso
,
A.
Zen
, and
A.
Michaelides
, “
Communication: Truncated non-bonded potentials can yield unphysical behavior in molecular dynamics simulations of interfaces
,”
J. Chem. Phys.
147
,
121102
(
2017
).
76.
S.
Menon
,
G.
Díaz Leines
,
R.
Drautz
, and
J.
Rogal
, “
Role of pre-ordered liquid in the selection mechanism of crystal polymorphs during nucleation
,”
J. Chem. Phys.
153
,
104508
(
2020
).
77.
M.
Fitzner
,
G. C.
Sosso
,
S. J.
Cox
, and
A.
Michaelides
, “
Ice is born in low-mobility regions of supercooled liquid water
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
2009
2014
(
2019
).
78.
R.
Becker
and
W.
Döring
, “
Kinetische behandlung der keimbildung in übersättigten dämpfen
,”
Ann. Phys.
416
,
719
752
(
1935
).
79.
K. A.
Jackson
, “
Nucleation
,” in
Kinetic Processes
(
Wiley-VCH Verlag GmbH & Co. KGaA
,
Weinheim, FRG
,
2005
), pp.
175
201
.
80.
J. F.
Lutsko
and
M. A.
Durán-Olivencia
, “
Classical nucleation theory from a dynamical approach to nucleation
,”
J. Chem. Phys.
138
,
244908
(
2013
).
81.
R.
Sear
, “
What do crystals nucleate on? What is the microscopic mechanism? How can we model nucleation?
,”
MRS Bull.
41
,
363
368
(
2016
).
82.
G. M.
Maggioni
and
M.
Mazzotti
, “
Stochasticity in primary nucleation: Measuring and modeling detection times
,”
Cryst. Growth Des.
17
,
3625
3635
(
2017
).
83.
S.
Jungblut
and
C.
Dellago
, “
Caveats of mean first-passage time methods applied to the crystallization transition: Effects of non-Markovianity
,”
J. Chem. Phys.
142
,
064103
(
2015
).
84.
A.
Gonzalo
, “
New methods: General discussion
,”
Faraday Discuss.
195
,
521
556
(
2016
).
85.
A.
Binder
,
T.
Lelièvre
, and
G.
Simpson
, “
A generalized parallel replica dynamics
,”
J. Comput. Phys.
284
,
595
616
(
2015
).
86.
P.
Tiwary
and
M.
Parrinello
, “
From metadynamics to dynamics
,”
Phys. Rev. Lett.
111
,
230602
(
2013
).
87.
S. A.
Khan
,
B. M.
Dickson
, and
B.
Peters
, “
How fluxional reactants limit the accuracy/efficiency of infrequent metadynamics
,”
J. Chem. Phys.
153
,
054125
(
2020
).
88.
K. M.
Bal
, “
Nucleation rates from small scale atomistic simulations and transition state theory
,” arXiv:2101.09234 [physics] (
2021
).
89.
J. D.
Honeycutt
and
H. C.
Andersen
, “
Small system size artifacts in the molecular dynamics simulation of homogeneous crystal nucleation in supercooled atomic liquids
,”
J. Phys. Chem.
90
,
1585
1589
(
1986
).
90.
H. E. A.
Huitema
,
J. P.
van der Eerden
,
J. J. M.
Janssen
, and
H.
Human
, “
Thermodynamics and kinetics of homogeneous crystal nucleation studied by computer simulation
,”
Phys. Rev. B
62
,
14690
14702
(
2000
).
91.
A.
Mahata
and
M.
Asle Zaeem
, “
Size effect in molecular dynamics simulation of nucleation process during solidification of pure metals: Investigating modified embedded atom method interatomic potentials
,”
Modell. Simul. Mater. Sci. Eng.
27
,
085015
(
2019
).
92.
E. D.
Zanotto
and
M.
Montazerian
, “
Dominant effect of heterogeneous dynamics on homogenous crystal nucleation in supercooled liquids
,”
Front. Phys.
8
,
20
(
2020
).
93.
U. R.
Pedersen
,
I.
Douglass
, and
P.
Harrowell
, “
How a supercooled liquid borrows structure from the crystal
,”
J. Chem. Phys.
154
,
054503
(
2021
).
94.
P. K.
Gupta
,
D. R.
Cassar
, and
E. D.
Zanotto
, “
Role of dynamic heterogeneities in crystal nucleation kinetics in an oxide supercooled liquid
,”
J. Chem. Phys.
145
,
211920
(
2016
).
95.
A. S.
Abyzov
,
V. M.
Fokin
,
N. S.
Yuritsyn
,
A. M.
Rodrigues
, and
J. W. P.
Schmelzer
, “
The effect of heterogeneous structure of glass-forming liquids on crystal nucleation
,”
J. Non-Cryst. Solids
462
,
32
40
(
2017
).
96.
F.
Puosi
and
A.
Pasturel
, “
Nucleation kinetics in a supercooled metallic glass former
,”
Acta Mater.
174
,
387
397
(
2019
).
97.
C.
Liu
,
G. P. F.
Wood
, and
E. E.
Santiso
, “
Modelling nucleation from solution with the string method in the osmotic ensemble
,”
Mol. Phys.
116
,
2998
3007
(
2018
).
98.
L.
Maragliano
,
A.
Fischer
,
E.
Vanden-Eijnden
, and
G.
Ciccotti
, “
String method in collective variables: Minimum free energy paths and isocommittor surfaces
,”
J. Chem. Phys.
125
,
024106
(
2006
).
99.
H.
Wang
,
C.
Schütte
, and
L.
Delle Site
, “
Adaptive resolution simulation (AdResS): A smooth thermodynamic and structural transition from atomistic to coarse grained resolution and vice versa in a grand canonical fashion
,”
J. Chem. Theory Comput.
8
,
2878
2887
(
2012
).
100.
L.
Delle Site
,
C.
Krekeler
,
J.
Whittaker
,
A.
Agarwal
,
R.
Klein
, and
F.
Höfling
, “
Molecular dynamics of open systems: Construction of a mean-field particle reservoir
,”
Adv. Theory Simul.
2
,
1900014
(
2019
).
101.
T.
Karmakar
,
P. M.
Piaggi
, and
M.
Parrinello
, “
Molecular dynamics simulations of crystal nucleation from solution at constant chemical potential
,”
J. Chem. Theory Comput.
15
,
6923
6930
(
2019
).
102.
C.
Perego
,
M.
Salvalaglio
, and
M.
Parrinello
, “
Molecular dynamics simulations of solutions at constant chemical potential
,”
J. Chem. Phys.
142
,
144113
(
2015
).
103.
A.
Barducci
,
G.
Bussi
, and
M.
Parrinello
, “
Well-tempered metadynamics: A smoothly converging and tunable free-energy method
,”
Phys. Rev. Lett.
100
,
020603
(
2008
).
104.
W.
Chen
and
Y.
Ren
, “
Molecular dynamics simulations of polymerisation and crystallisation at constant chemical potential
,”
Mol. Simul.
46
,
823
828
(
2020
).
105.
D.
Rozmanov
and
P. G.
Kusalik
, “
Temperature dependence of crystal growth of hexagonal ice (Ih)
,”
Phys. Chem. Chem. Phys.
13
,
15501
15511
(
2011
).
106.
J.
Kuipers
and
G. T.
Barkema
, “
Limitations of a Fokker-Planck description of nucleation
,”
Phys. Rev. E
82
,
011128
(
2010
).
107.
K.
Kawasaki
, “
Diffusion constants near the critical point for time-dependent Ising models. I
,”
Phys. Rev.
145
,
224
230
(
1966
).
108.
Y.
Lifanov
,
B.
Vorselaars
, and
D.
Quigley
, “
Nucleation barrier reconstruction via the seeding method in a lattice model with competing nucleation pathways
,”
J. Chem. Phys.
145
,
211912
(
2016
).
109.
S. V.
Krivov
, “
On reaction coordinate optimality
,”
J. Chem. Theory Comput.
9
,
135
146
(
2013
).
110.
B.
Peters
, “
Supersaturation rates and schedules: Nucleation kinetics from isothermal metastable zone widths
,”
J. Cryst. Growth
317
,
79
83
(
2011
).
111.
M.
Lauricella
,
G.
Ciccotti
,
N. J.
English
,
B.
Peters
, and
S.
Meloni
, “
Mechanisms and nucleation rate of methane hydrate by dynamical nonequilibrium molecular dynamics
,”
J. Phys. Chem. C
121
,
24223
24234
(
2017
).
112.
R. P.
Sear
, “
Quantitative studies of crystal nucleation at constant supersaturation: Experimental data and models
,”
CrystEngComm
16
,
6506
6522
(
2014
).
113.
P. C.
Myint
,
A. A.
Chernov
,
B.
Sadigh
,
L. X.
Benedict
,
B. M.
Hall
,
S.
Hamel
, and
J. L.
Belof
, “
Nanosecond freezing of water at high pressures: Nucleation and growth near the metastability limit
,”
Phys. Rev. Lett.
121
,
155701
(
2018
).
114.
B.
Cheng
,
G. A.
Tribello
, and
M.
Ceriotti
, “
The Gibbs free energy of homogeneous nucleation: From atomistic nuclei to the planar limit
,”
J. Chem. Phys.
147
,
104707
(
2017
).
115.
U.
Gasser
,
E. R.
Weeks
,
A.
Schofield
,
P. N.
Pusey
, and
D. A.
Weitz
, “
Real-space imaging of nucleation and growth in colloidal crystallization
,”
Science
292
,
258
262
(
2001
).
116.
X.-M.
Bai
and
M.
Li
, “
Test of classical nucleation theory via molecular-dynamics simulation
,”
J. Chem. Phys.
122
,
224510
(
2005
).
117.
G.
Díaz Leines
,
R.
Drautz
, and
J.
Rogal
, “
Atomistic insight into the non-classical nucleation mechanism during solidification in Ni
,”
J. Chem. Phys.
146
,
154702
(
2017
).
118.
S.
Prestipino
,
A.
Laio
, and
E.
Tosatti
, “
Systematic improvement of classical nucleation theory
,”
Phys. Rev. Lett.
108
,
225701
(
2012
).
119.
S.
Prestipino
,
A.
Laio
, and
E.
Tosatti
, “
Shape and area fluctuation effects on nucleation theory
,”
J. Chem. Phys.
140
,
094501
(
2014
).
120.
J. F.
Lutsko
, “
Systematically extending classical nucleation theory
,”
New J. Phys.
20
,
103015
(
2018
).
121.
S.
Ryu
and
W.
Cai
, “
Numerical tests of nucleation theories for the Ising models
,”
Phys. Rev. E
82
,
011603
(
2010
).
122.
B.
Peters
, “
Reaction coordinates and mechanistic hypothesis tests
,”
Annu. Rev. Phys. Chem.
67
,
669
690
(
2016
).
123.
P. J.
Steinhardt
,
D. R.
Nelson
, and
M.
Ronchetti
, “
Bond-orientational order in liquids and glasses
,”
Phys. Rev. B
28
,
784
805
(
1983
).
124.
G. A.
Tribello
,
F.
Giberti
,
G. C.
Sosso
,
M.
Salvalaglio
, and
M.
Parrinello
, “
Analyzing and driving cluster formation in atomistic simulations
,”
J. Chem. Theory Comput.
13
,
1317
1327
(
2017
).
125.
F.
Pietrucci
, “
Novel enhanced sampling strategies for transitions between ordered and disordered structures
,” in
Handbook of Materials Modeling
(
Springer
,
Cham
,
2020
), pp.
597
619
.
126.
M.
Fulford
,
M.
Salvalaglio
, and
C.
Molteni
, “
DeepIce: A deep neural network approach to identify ice and water molecules
,”
J. Chem. Inf. Model.
59
,
2141
(
2019
).
127.
Q.
Kim
,
J.-H.
Ko
,
S.
Kim
, and
W.
Jhe
, “
GCIceNet: A graph convolutional network for accurate classification of water phases
,”
Phys. Chem. Chem. Phys.
22
,
26340
26350
(
2020
).
128.
H.
Doi
,
K. Z.
Takahashi
, and
T.
Aoyagi
, “
Mining of effective local order parameters for classifying crystal structures: A machine learning study
,”
J. Chem. Phys.
152
,
214501
(
2020
).
129.
S. V.
Krivov
, “
Protein folding free energy landscape along the committor—The optimal folding coordinate
,”
J. Chem. Theory Comput.
14
,
3418
3427
(
2018
).
130.
Z.
Smith
,
D.
Pramanik
,
S.-T.
Tsai
, and
P.
Tiwary
, “
Multi-dimensional spectral gap optimization of order parameters (SGOOP) through conditional probability factorization
,”
J. Chem. Phys.
149
,
234105
(
2018
).
131.
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
).
132.
S.
Jungblut
and
C.
Dellago
, “
On the reaction coordinate for seeded crystallisation
,”
Mol. Phys.
113
,
2735
2741
(
2015
).
133.
C.
Desgranges
and
J.
Delhommelle
, “
Controlling polymorphism during the crystallization of an atomic fluid
,”
Phys. Rev. Lett.
98
,
235502
(
2007
).
134.
Y.
Liang
,
G.
Díaz Leines
,
R.
Drautz
, and
J.
Rogal
, “
Identification of a multi-dimensional reaction coordinate for crystal nucleation in Ni3Al
,”
J. Chem. Phys.
152
,
224504
(
2020
).
135.
P. R.
ten Wolde
,
M. J.
Ruiz-Montero
, and
D.
Frenkel
, “
Numerical evidence for bcc ordering at the surface of a critical fcc nucleus
,”
Phys. Rev. Lett.
75
,
2714
2717
(
1995
).
136.
S.
Jungblut
and
C.
Dellago
, “
Pathways to self-organization: Crystallization via nucleation and growth
,”
Eur. Phys. J. E
39
(
8
),
77
(
2016
).
137.
S.
Prestipino
, “
The barrier to ice nucleation in monatomic water
,”
J. Chem. Phys.
148
,
124505
(
2018
).
138.
R.
Yappert
,
K.
Kamat
, and
B.
Peters
, “
The overdamped transmission coefficient: Recovering the true mean first passage time from an inaccurate reaction coordinate
,”
J. Chem. Phys.
151
,
184108
(
2019
).
139.
G. T.
Beckham
and
B.
Peters
, “
Optimizing nucleus size metrics for liquid–solid nucleation from transition paths of near-nanosecond duration
,”
J. Phys. Chem. Lett.
2
,
1133
1138
(
2011
).
140.
W.
Lechner
,
C.
Dellago
, and
P. G.
Bolhuis
, “
Reaction coordinates for the crystal nucleation of colloidal suspensions extracted from the reweighted path ensemble
,”
J. Chem. Phys.
135
,
154110
(
2011
).
141.
A. M.
Berezhkovskii
and
A.
Szabo
, “
Diffusion along the splitting/commitment probability reaction coordinate
,”
J. Phys. Chem. B
117
,
13115
13119
(
2013
).
142.
H.
Sidky
,
Y. J.
Colón
,
J.
Helfferich
,
B. J.
Sikora
,
C.
Bezik
,
W.
Chu
,
F.
Giberti
,
A. Z.
Guo
,
X.
Jiang
,
J.
Lequieu
,
J.
Li
,
J.
Moller
,
M. J.
Quevillon
,
M.
Rahimi
,
H.
Ramezani-Dakhel
,
V. S.
Rathee
,
D. R.
Reid
,
E.
Sevgen
,
V.
Thapar
,
M. A.
Webb
,
J. K.
Whitmer
, and
J. J.
de Pablo
, “
SSAGES: Software suite for advanced general ensemble simulations
,”
J. Chem. Phys.
148
,
044104
(
2018
).
143.
E. E.
Borrero
and
F. A.
Escobedo
, “
Simulating the kinetics and thermodynamics of transitions via forward flux/umbrella sampling
,”
J. Phys. Chem. B
113
,
6434
6445
(
2009
).
144.
B.
Peters
and
B. L.
Trout
, “
Obtaining reaction coordinates by likelihood maximization
,”
J. Chem. Phys.
125
,
054108
(
2006
).
145.
A. K.
Faradjian
and
R.
Elber
, “
Computing time scales from reaction coordinates by milestoning
,”
J. Chem. Phys.
120
,
10880
10889
(
2004
).
146.
E. E.
Santiso
and
B. L.
Trout
, “
A general method for molecular modeling of nucleation from the melt
,”
J. Chem. Phys.
143
,
174109
(
2015
).
147.
N.
Malik
and
U.
Ozturk
, “
Rare events in complex systems: Understanding and prediction
,”
Chaos
30
,
090401
(
2020
).
148.
F.
Ragone
,
J.
Wouters
, and
F.
Bouchet
, “
Computation of extreme heat waves in climate models using a large deviation algorithm
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
24
29
(
2018
).
149.
R. J.
Webber
,
D. A.
Plotkin
,
M. E.
O’Neill
,
D. S.
Abbot
, and
J.
Weare
, “
Practical rare event sampling for extreme mesoscale weather
,”
Chaos
29
,
053109
(
2019
).
150.
L.
Billings
and
E.
Forgoston
, “
Seasonal forcing in stochastic epidemiology models
,”
Ric. Mat.
67
,
27
47
(
2018
).
151.
T.
Grafke
and
E.
Vanden-Eijnden
, “
Numerical computation of rare events via large deviation theory
,”
Chaos
29
,
063118
(
2019
).