The Gibbs free energy of homogeneous nucleation: from atomistic nuclei to the planar limit

TL;DR

This paper presents a method to improve the fitting of macroscopic nucleation models by combining atomistic simulations of nuclei and planar interfaces, reducing errors and enhancing prediction reliability.

Contribution

It introduces a combined fitting approach using both atomistic and planar interface simulations to better determine nucleation model parameters.

Findings

Fitting parameters from atomistic free energy profiles is highly sensitive to errors.

Combining planar interface and nucleus simulations reduces fitting errors.

The new method yields more reliable predictions for large nuclei.

Abstract

In this paper we discuss how the information contained in atomistic simulations of homogeneous nucleation should be used when fitting the parameters in macroscopic nucleation models. We show how the number of solid and liquid atoms in such simulations can be determined unambiguously by using a Gibbs dividing surface and how the free energy as a function of the number of solid atoms in the nucleus can thus be extracted. We then show that the parameters of a model based on classical nucleation theory can be fit using the information contained in these free-energy profiles but that the parameters in such models are highly correlated. This correlation is unfortunate as it ensures that small errors in the computed free energy surface can give rise to large errors in the extrapolated properties of the fitted model. To resolve this problem we thus propose a method for fitting macroscopic nucleation models that uses simulations of planar interfaces and simulations of three-dimensional nuclei in tandem. We show that when the parameters of the macroscopic model are fitted in this way the numerical errors for the final fitted model are smaller and that the extrapolated predictions for large nuclei are thus more reliable.

Figures

Tables

Table 1. Table 1: A comparison of the predictions for μ sl subscript 𝜇 sl \mu_{\text{sl}} and γ ∞ Φ superscript subscript 𝛾 Φ \gamma_{\infty}^{\Phi} from different models. For each number, the value in the bracket indicates the statistical uncertainty in the last digit.

	${\mu }_{\text{sl}}$				${\gamma }_{\infty }^{\Phi }$
model	0.56	0.58	0.60	0.6185	0.56	0.58	0.60	0.6185
CNT($\Phi$)+Tol	-0.100(9)	-0.061(5)	-0.033(4)	0.001(5)	0.42(4)	0.38(2)	0.39(2)	0.37(2)
Planar	-0.0991(2)	-0.0648(2)	-0.0305(2)	0.0012(2)	0.418(4)	0.393(3)	0.377(2)	0.360(2)
CNT($\Phi$)	-0.083(2)	-0.049(1)	-0.016(1)	0.017(1)	0.346(5)	0.330(2)	0.317(1)	0.302(3)

Equations

$$G(n)=\mu n+\gamma \Omega v^{\frac{2}{3}}{n}^{\frac{2}{3}},$$

$${\rho }_{\mathrm{sl}}\left(\Phi |n_{\text{s}}(\Phi ),N-n_{\text{s}}(\Phi )\right)\equiv {\rho }_{\text{ref}}\left(\Phi |n_{\text{s}}(\Phi ),N-n_{\text{s}}(\Phi )\right).$$

$$\begin{gathered} e^{-\beta \tilde{G}(\Phi )}= \\ {\int }_{n_{\text{cut}}}^{n^{\star }}\mathrm{d}{n}_{\text{s}}(\Phi ){\rho }_{\mathrm{sl}}\left(\Phi |n_{\text{s}}(\Phi ),N-n_{\text{s}}(\Phi )\right)N_s{e}^{-\beta G(n_{\text{s}}(\Phi ))}, \end{gathered}$$

$$G(n_{\text{s}}(\Phi ))={\mu }_{\text{sl}}{n}_{\text{s}}(\Phi )+{\gamma }_{\text{sl}}^{\Phi }A(n_{\text{s}}(\Phi )),$$

$$\delta {n_{\text{s}}}^{\Phi ,\Theta }=(n_{\text{s}}(\Theta )-n_{\text{s}}(\Phi ))/A$$

$$\Delta {\gamma }_{\infty }^{\Phi ,\Theta }={\gamma }_{\infty }^{\Theta }-{\gamma }_{\infty }^{\Phi }=-{\mu }_{\text{sl}}\delta {n_{\text{s}}}^{\Phi ,\Theta },$$

$$G(n_{\text{s}}(V))={\mu }_{\text{sl}}{n}_{\text{s}}(V)+{\gamma }_{\text{sl}}^V\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(V).$$

$${\gamma }_{\text{sl}}^V(R)={\gamma }_{\infty }^V(1-2\delta /R+\mathcal{O}(1/R^2)).$$

$$G(n_{\text{s}}(V))\eqsim {\mu }_{\text{sl}}{n}_{\text{s}}(V)+{\gamma }_{\infty }^V\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(V)(1-ϵ{n_{\text{s}}}^{-\frac{1}{3}}(V)),$$

$$n_{\text{s}}(V)-n_{\text{s}}(\Phi )\eqsim \delta {n_{\text{s}}}^{\Phi ,V}\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(V),$$

$$G(n_{\text{s}}(\Phi ))\eqsim {\mu }_{\text{sl}}{n}_{\text{s}}(\Phi )+{\gamma }_{\infty }^{\Phi }\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(\Phi )(1+\zeta {n_{\text{s}}}^{-\frac{1}{3}}(\Phi )),$$

Full text

nucleation, atomistic simulation, classical nucleation theory, phase-field models, metadynamics

††preprint: APS/123-QED

I Introduction

Nucleation is an important component of many technological and natural processes, including metal casting, the assembly of microtubules in cells and the formation of water droplets and ice crystals in the atmosphere  Oxtoby (1992); Yi and Rutledge (2012); Sosso et al. (2016). In spite of its generic importance, however, our understanding of nucleation kinetics is somewhat unsatisfactory, even for the simplest cases of homogeneous nucleation. The main reason for this is that the experimental determination of nucleation rates is challenging. Many of the available analytical techniques can only be used to perform ex situ analysis, which cannot fully capture the dynamical processes that take place during nucleation.

A considerable number of studies in the last two decades have thus used atomistic simulation to study homogeneous nucleation ten Wolde et al. (1996); ten Wolde and Frenkel (1999); Moroni et al. (2005); Trudu et al. (2006); Lechner et al. (2011); Prestipino et al. (2012); Yi and Rutledge (2012); McCarty et al. (2016). Much of this computational effort has concentrated on investigating how the free energy changes with cluster size $G(n)$ because the homogeneous nucleation barrier $G^{\star}$ is the maximum of $G(n)$ and because the nucleation rate can be determined from $G^{\star}$ by using transition state theory or the Bennett-Chandler approach ten Wolde et al. (1996); Moroni et al. (2005); Jungblut and Dellago (2016). Researchers have thus been able to perform atomistic simulations in which solid critical nuclei containing hundreds of atoms have been observed to form from the melt ten Wolde et al. (1996); ten Wolde and Frenkel (1999); Moroni et al. (2005); Lechner et al. (2011); Prestipino et al. (2012).

When analyzing atomistic simulations of nucleation, macroscopic models that provide analytical expressions for the nucleation free energy profile $G(n)$ can be extremely useful as they can be used to extrapolate the $G(n)$ obtained for the sub-critical nuclei which form in atomistic simulations so as to obtain the free energies of larger nuclei. This is often essential when predicting $G^{\star}$, and when interpolating between results obtained at different temperatures. One simple model which is commonly used to rationalize nucleation is classical nucleation theory (CNT). This theory assumes that $G(n)$ can be expressed as the sum of a bulk and a surface term, i.e.

[TABLE]

Here $\mu$ is the chemical potential difference between the stable and the metastable phases, $\gamma$ is the area specific interfacial free energy, $v$ is the molar volume of the bulk stable phase, and $\Omega$ is a geometrical constant (e.g. $(36\pi)^{\frac{1}{3}}$ for a spherical nucleus).

CNT has been shown to give rise to significant systematic errors for all sizes of the nucleus Prestipino et al. (2012). A large part of problem is that the number of atoms $n$ in the cluster is not a well-defined quantity Moroni et al. (2005); Lechner et al. (2011); Prestipino et al. (2012); Cheng et al. (2015). Different definitions of $n$ inevitably affect the calculated free energy profile $G(n)$, which forces one to ask which definition of $n$ is most appropriate within the CNT framework. Another part of the problem is that the CNT expression in Eqn. (1) is generally used as a fitting model for $G(n)$, because the various physical quantities that enter this expression - the chemical potential $\mu$ and the interfacial free energy $\gamma$ - cannot be calculated directly from homogeneous nucleation simulations. As a consequence, when the results from simulations deviate from the predictions of CNT, it is not clear which assumption within the model is the main culprit. It is thus difficult to add additional correction terms to the model to compensate for the missing elements.

An alternative method for probing nucleation involves investigating planar interfaces, which resemble the surfaces of a large nucleus. From the simulations of planar interfaces, the chemical potential as well as the interfacial free energies for different crystallographic directions can be evaluated Hoyt et al. (2001a); Davidchack and Laird (2000, 2003); Angioletti-Uberti et al. (2010a); Pedersen et al. (2013); Cheng et al. (2015). Various techniques have been developed for computing the interfacial free energy or the stiffness of the planar interfaces at the coexistence temperature Hoyt et al. (2001a); Davidchack and Laird (2000, 2003); Angioletti-Uberti et al. (2010a); Pedersen et al. (2013). Using the values of $\gamma$ obtained from such simulations within the framework of CNT is problematic though as nucleation occurs in undercooled conditions. Recently, however, the metadynamics method has been extended so that it also works in out-of-equilibrium conditions away from the melting point Cheng et al. (2015). This method thus opens the door to a systematic exploration into how the quantities obtained from simulations of planar interfaces can be used to predict the free energy associated with a three-dimensional nucleus under realistic conditions.

In the present study, we apply our recently introduced thermodynamic framework for defining the nucleation free energy profile $G(n)$ in an unambiguous way Cheng and Ceriotti (2017). We then fit the computed $G(n)$ to the expressions that come from macroscopic nucleation models, using the specific interfacial free energy $\gamma$ and the chemical potential difference $\mu$ as fitting parameters. We also show that a higher-order correction – generally referred to as the Tolman term – naturally emerges from our formalism. Finally, in the last part of the paper, we compare the predicted values for $\gamma$ and $\mu$ that we obtain from simulations of three-dimensional nuclei against those obtained from the simulations of planar interfaces.

II Theory

II.1 Gibbs dividing surface in atomistic systems

In order to link the atomistic descriptions and the macroscopic representations of a solid-liquid system, we define a Gibbs dividing surface using an arbitrarily-chosen extensive quantity $\Phi$ that characterises the state of all the atoms within the system. This extensive quantity might be the volume $V$, the entropy $S$, Piaggi et al. (2016) or a global collective variable (CV) $\Phi=\sum_{i}\phi(i)$ that is constructed by summing the order parameter values $\phi(i)$ for each of the atoms in the microstate Cheng et al. (2015); Cheng and Ceriotti (2017). We define the position of the $\Phi$-based Gibbs dividing surface by requiring that there be zero surface excess for the extensive quantity $\Phi$. In other words, the real system and a reference system in which the presence of the interface has no effect on the properties of the two bulk phases are setup so they have the same values for the extensive quantity Gibbs et al. (1928); Tolman (1948). Thermal fluctuations ensure that the instantaneous value of the extensive quantity $\Phi$ for this reference system is not fixed, but follows a distribution $\rho_{\text{ref}}\left(\Phi\middle|{n_{\text{s}}}(\Phi),N-{n_{\text{s}}}(\Phi)\right)$ with a finite width that can be calculated by performing simulations of the two bulk phases Yukalov (2011); Cheng and Ceriotti (2017). To have zero excess for $\Phi$ the real system and the reference system, both of which have ${n_{\text{s}}}(\Phi)$ solid atoms and $N-{n_{\text{s}}}(\Phi)$ liquid atoms, should have the same distribution for $\Phi$ Cheng and Ceriotti (2017), i.e.

[TABLE]

In Ref. Cheng and Ceriotti (2017), we argued that Eqn. (2) can be used to link a free energy profile $\tilde{G}(\Phi)$, which is expressed as a function of an extensive quantity $\Phi$ for the whole system, with the nucleation free energy for a single solid cluster which is expressed as a function of ${n_{\text{s}}}(\Phi)$. This connection can be made under the reasonable assumption that the probability of observing a single large sub-critical cluster inside an undercooled bulk liquid can be expressed as $N_{s}\exp(-\beta G({n_{\text{s}}}(\Phi)))$ when $n_{\text{cut}}\leq{n_{\text{s}}}(\Phi)\leq n^{\star}$. In these expressions $N_{s}$ is the number of nucleation sites, which in this work we take to be simply the number of atoms, while $G({n_{\text{s}}}(\Phi))$ is the free energy excess associated with the solid cluster. Using the law of total probability, the probability distribution for $\Phi$ in such systems follows

[TABLE]

This expression is valid for values of $\Phi$ that satisfy $\rho_{\mathrm{sl}}\left(\Phi\middle|{n_{\text{s}}}(\Phi),N-{n_{\text{s}}}(\Phi)\right)\approx 0$ for all ${n_{\text{s}}}(\Phi)<n_{\text{cut}}$. To be clear $\tilde{G}(\Phi)$ can be directly computed from atomistic simulations of nucleation. Values of $G({n_{\text{s}}}(\Phi))$ for ${n_{\text{s}}}(\Phi)>n_{\text{cut}}$ can then be found by inverting Eqn.(3) numerically Cheng and Ceriotti (2017).

II.2 A general formulation of classical nucleation theory

With a $\Phi$-based dividing surface, the free energy of the solid-liquid system relative to that of a bulk liquid can be naturally decomposed into a bulk and a surface term:

[TABLE]

Here $\mu_{\text{sl}}=\mu_{\text{s}}-\mu_{\text{l}}$ is the difference between the per-atom chemical potentials of the solid and liquid phases. The free energy excess associated with the $\Phi$-based dividing surface is then the product of a specific energy term $\gamma_{\text{sl}}^{\Phi}$ and an extensive area term $A({n_{\text{s}}}(\Phi))$.

For planar interfaces, the surface area $A({n_{\text{s}}}(\Phi))$ in Eqn. (4) is fixed by the boundary conditions so the specific planar interfacial free energy $\gamma_{\text{sl}}^{\Phi}=\gamma_{\infty}^{\Phi}$ is a constant that depends on the crystallographic direction of the planar interface. Notice that if the extensive quantity $\Theta$ is used to define the dividing surface instead of $\Phi$ the composition of the reference system changes. This is important as a change in the composition of the reference state will affect the value obtained for the interfacial free energy. We quantify this change in composition using:

[TABLE]

and further note that this quantity should be constant as replacing $\Phi$ with $\Theta$ shifts the location of the dividing surface by a fixed amount Cheng et al. (2015). Furthermore, when $\Phi$ is replaced by $\Theta$ the resulting change in the planar interfacial free energy

[TABLE]

is also a constant, because the total free energy in Eqn. (4) should be unaffected by the change in the extensive quantity Cheng et al. (2015).

For the curved interfaces around a three-dimensional nucleus, the surface area $A({n_{\text{s}}}(\Phi))$ is not equally well-defined. Classical nucleation theory assumes that the nucleus has the same density as the bulk solid, which is what ensures that $A=\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}$. Critically, however, this bulk density assumption is only valid when the equimolar dividing surface is used to define the reference state, which implies that only quantities determined using a dividing surface with no excess volume should be inserted into Eqn. (4). This equation should thus read:

[TABLE]

Even when the reference state is defined in a manner consistent with this assumption about the bulk density of the nucleus, however, one still has to include a term that incorporates the surface excess free-energies’ dependence on curvature into $\gamma_{\text{sl}}^{V}$ Tolman (1948); Block et al. (2010). This dependence of $\gamma_{\text{sl}}^{V}$ on the effective radius of the nucleus $R=(3/4\pi)^{\frac{1}{3}}v_{\text{s}}^{\frac{1}{3}}{n_{\text{s}}}^{\frac{1}{3}}(V)$ can be written as

[TABLE]

In this expression $\delta$ is the Tolman length - a quantity that measures the difference between the locations of the equimolar dividing surface and the surface of tension at the planar limit. $\gamma_{\infty}^{V}$, meanwhile, is the specific interfacial free energy associated with the equimolar surface in the planar limit. Taking only the leading terms in Eqn. (8), Eqn. (7) can be rewritten as

[TABLE]

where $\epsilon=(32\pi/3)^{\frac{1}{3}}v_{\text{s}}^{-\frac{1}{3}}\delta$ is a constant.

At the atomic scale, due to the magnitude of fluctuations, it is difficult to distinguish between the different phases based only on the difference in molar volume. Hence, it is often not convenient to use this quantity to determine the number of solid particles in a simulation. One can use a different extensive quantity to determine the location of the diving surface as it is possible to convert the free-energy profile obtained using one dividing surface to the result that would have been obtained for a different dividing surface using Eq. (5) Cheng et al. (2015). It would be useful, however, to derive a CNT expression based on an arbitrary choice of reference thermodynamic variable $\Phi$. To derive such an expression, we first assumed that the difference between the location of the $\Phi$-based dividing surface and the location of the equimolar dividing surface is much smaller than the effective radius of the nucleus $R$, as this ensures that the difference between ${n_{\text{s}}}(\Phi)$ and ${n_{\text{s}}}(V)$ can be approximated using

[TABLE]

where $\delta{n_{\text{s}}}^{\Phi,V}$ is the constant that appears in Eqn. (5) and  (6) and where $\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(V)$ is a measure of the surface area of the nucleus. We can substitute ${n_{\text{s}}}(V)$ for ${n_{\text{s}}}(\Phi)$ in Eqn. (9) using the above approximation, and can also use the relation $\gamma_{\infty}^{V}-\gamma_{\infty}^{\Phi}=-\mu_{\text{sl}}\delta{n_{\text{s}}}^{\Phi,V}$. After dropping the higher order correction terms, we obtain

[TABLE]

where $\zeta=\frac{2}{3}\delta{n_{\text{s}}}^{\Phi,V}\Omega v_{\text{s}}^{\frac{2}{3}}-\epsilon$ is a constant. This $\zeta{n_{\text{s}}}^{-\frac{1}{3}}(\Phi)$ term in Eqn. (11), which has the same mathematical form as the Tolman correction, has two components: the $\delta{n_{\text{s}}}^{\Phi,V}$ term stems from the difference between the $\Phi$-based dividing surface and the equimolar surface, while the other $\epsilon$ term is determined from the Tolman length and has a value that is independent of the choice of the dividing surface. Hereafter, we will use CNT($\Phi$)+Tol to denote the expression in Eqn. (11) in order to highlight its dependence on the chosen extensive quantity. An awareness of this dependence on the choice of $\Phi$ is crucial when comparing studies performed with different protocols. Eqs. (9) and (11) are both valid expressions for the free energy of a nucleus that are consistent with Tolman-corrected CNT. These two expressions differ, however, when they come to defining the size of the nucleus and the value and interpretation of the planar-interface excess free energy and the coefficient for the finite-size correction.

III Simulation methods

We simulated the processes of homogeneous nucleation for a simple but realistic Lennard-Jones system Davidchack and Laird (2003); Angioletti-Uberti et al. (2010b); Benjamin and Horbach (2014). The NPT ensemble was employed throughout with the Nose-Hoover thermostat and isotropic barostat. The time step was set equal to 0.004 Lennard-Jones time units and a supercell containing 23328 atoms was used throughout.

Each independent simulation at each temperature was run for approximately $6\times 10^{6}$ steps. To accelerate the sampling so as to obtain reversible formation of a solid nucleus in a viable amount of simulation time, we performed biased sampling using the well-tempered metadynamics protocol with adaptive Gaussians  Barducci et al. (2008); Branduardi et al. (2012), and the collective variable $\Phi=\sum_{i}S(\kappa(i))$ that was employed in Ref. Cheng et al. (2015). In essence, $S(\kappa(i))$ is a local structural fingerprint for atom $i$, which takes a value of $1$ for a perfect fcc crystal and [math] for a homogeneous liquid. Fast implementation of this complex simulation setup was made possible by the flexibility of the PLUMED code Tribello et al. (2014) in combination with LAMMPS Plimpton (1995). See the Supplemental Material for sample input files SI . In addition to the simulations of homogeneous nucleation, we also simulated planar interfaces along the $\langle{100}\rangle$, $\langle{111}\rangle$, and $\langle{110}\rangle$ crystallographic directions of the fcc lattice. For the planar interface simulations, the setups are identical to the ones described in Ref. Cheng et al. (2015).

IV Results and Discussions

At each of the temperatures $0.56$, $0.58$, $0.60$ and $0.6185$, a total of 12 independent metadynamics runs were performed. For each simulation run we first computed the free energy profile $\tilde{G}(\Phi)$ with respect to the extensive quantity $\Phi$ of the system. We then extracted the nucleation free energy for a single solid cluster as a function of ${n_{\text{s}}}(\Phi)$ using the framework introduced above, which is thoroughly described in Ref. Cheng and Ceriotti (2017). Each free energy profile $G({n_{\text{s}}}(\Phi))$ is plotted as a thin blue curve in Figure 1. For each $G({n_{\text{s}}}(\Phi))$, we then performed a CNT($\Phi$)+Tol fit using Eqn. (11) with $\mu_{\text{sl}}$, $\gamma_{\infty}^{\Phi}\Omega$, and $\zeta$ as fitting parameters. Each of these fitted curves is shown in red in Figure 1.

The red curves in Figure 1 start to diverge at larger sizes even though they almost overlap at small sizes. In order to understand the origin of this divergence we show the two parameters ($\mu_{\text{sl}}$, $\gamma_{\infty}\Omega$) obtained from CNT+Tol($\Phi$) fits for each free energy profile at $T=0.60$ in Figure 2. The two parameters $\mu_{\text{sl}}$ and $\gamma_{\infty}\Omega$ are clustered around a straight regression line in Figure 2, which suggests that there is a very large correlation between them. The consequence of this correlation is that a tiny change in the segment of data used for fitting can make the two parameters vary collectively along the regression line by a significant amount. Thus, small uncertainties in the computed free energy profile for the small sub-critical nuclei that form in simulations can propagate and amplify when these curves are extrapolated to large nuclei.

We wanted to determine whether or not the values obtained for $\mu_{\text{sl}}$ and $\gamma_{\infty}\Omega$ from the CNT($\Phi$)+Tol model are physically meaningful. To this end, we also show results obtained from simulations of the planar interfaces in Figure 2. The value of $\mu_{\text{sl}}$ for an interface perpendicular to any lattice direction can be computed by performing biased simulations of that interface Pedersen et al. (2013); Cheng et al. (2015). We computed these planar interfacial free energies in our previous work Cheng et al. (2015) and used a deterministic framework to locate the Gibbs dividing surface. In the Supplementary Material, we show that the values of the interfacial free energies do not change significantly when the probabilistic framework in Eqn. (2) is applied so for this reason we chose to adopt the values for $\gamma^{\Phi}$ and $\gamma^{V}$ reported in Ref. Cheng et al. (2015). To estimate $\gamma_{\infty}^{\Phi}\Omega$, we first obtained the values of $\gamma_{100}^{\Phi}$, $\gamma_{111}^{\Phi}$ and $\gamma_{110}^{\Phi}$ for planar interfaces with normal vectors parallel to the specified lattice directions. We then used a common assumption; namely, that for this Lennard-Jones system the interfacial free energy surface $\gamma^{\Phi}(\vec{n})$ can be expanded using cubic harmonics to the third order Hoyt et al. (2001b), whose coefficients can be parametrized using the values of $\gamma_{100}^{\Phi}$, $\gamma_{111}^{\Phi}$ and $\gamma_{110}^{\Phi}$. As shown in the Supplementary Material, the shape of the equimolar surface of the nucleus $R(\vec{n})$ can then be reproduced by performing a Wulff construction using the $\gamma^{V}$ values Adams et al. (1994); SI . As the difference between the locations of the $\Phi$-based dividing surface and the equimolar dividing surface was assumed to be insignificant compared to the radius of the nucleus, the shape, $R(\vec{n})$, that emerges from the Wulff construction can be considered to be very close to the shape of the $\Phi$-based dividing surface of the nucleus. From this shape, we can, therefore, estimate the geometrical constant $\Omega$, and an effective value of $\gamma_{\infty}^{\Phi}$ for the surface of the whole nucleus by computing the surface integral $\iint_{R}\gamma^{\Phi}(\vec{n})\,dA/\iint_{R}\,dA$. At all the temperatures we considered, the shapes $R(\vec{n})$ that we obtained were close to spherical. In fact, the corresponding $\Omega$ values were within 0.5% of the geometrical constant for a sphere.

Figure 2 suggests that the values of $\mu_{\text{sl}}$ and $\gamma_{\infty}^{\Phi}\Omega$ from the CNT($\Phi$)+Tol fits are consistent with the planar limit results. However, there is a large spread in the parameter values because of the strong correlation between these two parameters in the fitting. Table 1 shows all the values of $\mu_{\text{sl}}(T)$ together with the values of $\gamma_{\infty}(T)$ that are predicted both from the CNT($\Phi$)+Tol model and the planar limit results and further illustrates the very good agreement between the two models at all temperatures considered. We think this good agreement is a strong indication that the CNT($\Phi$)+Tol model is a good macroscopic model for the free energies of atomistic nuclei. Indeed, when a different model is used, the values of the two parameters $\mu_{\text{sl}}$ and $\gamma_{\infty}^{\Phi}\Omega$ may not agree with planar limit results. For example, we also performed a so-called CNT($\Phi$) fit on each computed $G({n_{\text{s}}}(\Phi))$, using a conventional CNT formulation $G({n_{\text{s}}}(\Phi))=\mu_{\text{sl}}{n_{\text{s}}}(\Phi)+\gamma_{\infty}^{\Phi}\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(\Phi)$ SI . The values of $\mu_{\text{sl}}$ and $\gamma_{\infty}^{\Phi}\Omega$ from the CNT($\Phi$) fits are indicated in Figure 2 and Table 1. Once again there are very strong inter-correlations between the parameters. What is even worse, however, is that the fitted values are no longer consistent with the planar interface results. This example highlights the perils associated with using the wrong model to fit the free energy profile for nucleation: although one might obtain a good fit for the data points, the values of the fitted parameters of the model are physically unrealistic and the model will most likely have limited predictive power in scenarios that are outside the range of the existing data set.

The comparison between the values of $\mu_{\text{sl}}$ and $\gamma_{\infty}^{\Phi}\Omega$ from the CNT($\Phi$)+Tol and the planar limit results suggests that one should just use the results obtained at the planar interface limit in Eqn. (11) so as to have a model that has $\zeta$ as the only fitting parameter. This model would have much greater statistical accuracy as there are no correlations between the fitting parameters. We performed such fittings and compared the values of the parameter $\zeta$ with the values obtained from the regular CNT($\Phi$)+Tol fits using all three parameters. The bottom panel of Figure 3, shows that this new approach results in much smaller uncertainties in the values of $\zeta$. Therefore, performing simulations of planar interfaces as well as simulations of three dimensional nuclei is worthwhile when studying homogeneous nucleation, especially given how computationally inexpensive such simulations are. We would recommend first computing the values of $\mu_{\text{sl}}$ and $\gamma_{\infty}^{\Phi}\Omega$ from simulations of planar interfaces. This allows one to determine the leading terms $\mu_{\text{sl}}{n_{\text{s}}}(\Phi)+\gamma_{\infty}^{\Phi}\Omega v_{\text{s}}^{\frac{2}{3}}{n_{\text{s}}}^{\frac{2}{3}}(\Phi)$ in a CNT-type model. Comparing these leading terms with the actual nucleation free energy profiles computed from the simulations of three-dimensional nuclei, then allows one to extract the other higher-order correction terms in the CNT-type model with much higher statistical accuracy.

It is worth discussing the higher-order correction term $\zeta=\frac{2}{3}\delta{n_{\text{s}}}^{\Phi,V}\Omega v_{\text{s}}^{\frac{2}{3}}-\epsilon$ in Eqn. (11) a little further. This factor has two distinct components but we treat it as a single fitting parameter in the present study. To evaluate the first term separately we can exploit the fact that the difference between the number of solid-like atoms per area $\delta{n_{\text{s}}}^{\Phi,V}$ for different definitions of the dividing surface can be directly evaluated from simulations of planar interfaces as discussed in Ref. Cheng et al. (2015). We evaluated $\delta{n_{\text{s}}}^{\Phi,V}$ for planar interfaces with normal vectors parallel to the $\langle{100}\rangle$, $\langle{111}\rangle$ and $\langle{110}\rangle$ lattice directions at each temperature, and then approximated $\delta{n_{\text{s}}}^{\Phi,V}(\vec{n})$ using a cubic harmonic expansion, and thus estimated its averaged value for a three-dimensional nucleus. These estimated values for $\frac{2}{3}\delta{n_{\text{s}}}^{\Phi,V}\Omega v_{\text{s}}^{\frac{2}{3}}$ are shown in the bottom panel of Figure 3. The other term $\epsilon=(32\pi/3)^{\frac{1}{3}}v_{\text{s}}^{-\frac{1}{3}}\delta$ that enters $\zeta$ can be determined from the Tolman length $\delta$. This quantity can be evaluated from the pressure tensor at each position along the normal direction of a planar interface. However, because calculations of this sort are difficult to converge Lei et al. (2005); Blokhuis and Kuipers (2006); van Giessen and Blokhuis (2009), we chose to subtract the values we obtained for $\frac{2}{3}\delta{n_{\text{s}}}^{\Phi,V}\Omega v_{\text{s}}^{\frac{2}{3}}$ from the value we obtained for $\zeta$ and to extract the Tolman length from $\epsilon$. In this way, we predicted the Tolman length to be on the order of $-0.3$ Lennard-Jones length units at the temperatures considered.

Finally, we want to point out that the fact that the CNT($\Phi$)+Tol model succeeds in describing small nuclei that contain hundreds of atoms, while being consistent with simulations performed in the limit of a planar solid-liquid interface has theoretical implications. Even though the actual interface between the solid and the liquid is diffuse and fluctuating, the concept of the Gibbs dividing surface can be used to convert the atomistic descriptions of the solid-liquid system into a macroscopic representation. Moreover, finite-size effects are well-captured by simple corrections that take into account the deviation between the chosen dividing surface, the equimolar dividing surface and the surface of tension. In this case, these two corrections are of comparable size.

V Conclusions

In summary we have studied, by means of atomic-scale simulations, the various different terms that appear in classical nucleation theory. When analysing our simulations we define a Gibbs dividing surface from the value of a macroscopic order parameter when calculating the dependence of the free energy on the number of particles in the nucleus. Then, by juxtaposing explicit simulations of a three-dimensional nucleus with simulations of a planar solid-liquid interface, performed in equivalent thermodynamic conditions, we identify the effects that lead the results on small nuclei to deviate from the predictions of classical nucleation theory. We find that deviations occur because curvature-dependent corrections to the planar-interface surface energy (the so-called Tolman term) are required. These genuine departures from CNT should not be confused with the deviations that occur when an arbitrary order parameter is used to identify the nucleus, however. Such corrections are required simply because we require zero surface excess for the arbitrary extensive quantity $\Phi$ when calculating the surface excess free energy. As $\Phi$ is not the volume the value we obtain for this excess free energy differs from the value that would have been obtained had we found the true equimolar dividing surface. The leading order correction for this effect has the same functional form as the Tolman term. It should not, however, be considered as a genuine departure from CNT as it is simply an artifact in the analysis.

In the final parts of the manuscript we discussed the statistical efficiency of different approaches for finding the CNT parameters. We showed that the surface and bulk parameters in the CNT model are strongly correlated, which makes the fitted values of these parameters extraordinarily sensitive to small deviations in the computed free energy surfaces that are used to fit them. We thus concluded that fitting these terms using a free energy surface obtained from a simulation of a 3D nucleus is not optimal and that also using information from simulations of planar-interfaces is thus beneficial. In fact, given that planar-interface models converge faster with respect to both size and simulation time, and that such simulations can be more easily analyzed in terms of different order parameters, we suggest that they should always be used in a preliminary phase. Explicit simulations of 3D nucleation can then be used to identify genuine finite-size effects. We hope that our careful analysis will help resolve some of the ambiguities in atomic-scale studies of both homogeneous and heterogeneous nucleation, and that this work will lay the foundations for the modelling of more complex materials and for the study of different kinds of phase transitions.