Arguably, the main challenge of nucleation theory is to accurately evaluate the work of formation of a critical embryo in the new phase, which governs the nucleation rate. In Classical Nucleation Theory (CNT), this work of formation is estimated using the capillarity approximation, which relies on the value of the planar surface tension. This approximation has been blamed for the large discrepancies between predictions from CNT and experiments. In this work, we present a study of the free energy of formation of critical clusters of the Lennard-Jones fluid truncated and shifted at 2.5σ using Monte Carlo simulations, density gradient theory, and density functional theory. We find that density gradient theory and density functional theory accurately reproduce molecular simulation results for critical droplet sizes and their free energies. The capillarity approximation grossly overestimates the free energy of small droplets. The incorporation of curvature corrections up to the second order with the Helfrich expansion greatly remedies this and performs very well for most of the experimentally accessible regions. However, it is imprecise for the smallest droplets and largest metastabilities since it does not account for a vanishing nucleation barrier at the spinodal. To remedy this, we propose a scaling function that uses all relevant ingredients without adding fitting parameters. The scaling function reproduces accurately the free energy of the formation of critical droplets for the entire metastability range and all temperatures examined and deviates from density gradient theory by less than one kBT.

Crystallization, condensation, and cavitation are phase transitions that govern a myriad of practical applications and natural phenomena.1–5 They are initiated by thermal fluctuations that lead to the formation of a nano-sized embryo of the new phase. In order to grow and for the phase transition to proceed, the size of the embryo must be sufficiently large to overcome the free energy cost of forming a surface. The thermodynamic principles of such processes were explained in the seminal work of Gibbs.6 Together with the theory to describe the kinetics developed by Volmer and Weber,7 Farkas,8 Becker and Döring,9 and Zeldovich,10 among others, established the basis of Classical Nucleation Theory (CNT), which still today leads our current understanding of nucleation phenomena. According to CNT, the rate of formation of embryos in the new phase per unit volume (J), called the nucleation rate, is described by the following formula:
J=Kexp(W/kBT),
(1)
where W is the work of formation of the critical embryo, kB is Boltzmann’s constant, T is the temperature, and K is a kinetic pre-factor.

Since nucleation is an activated process, the pressure of a vapor can be increased far above the saturation pressure, and the temperature of a liquid can be decreased well below the melting temperature before condensation or crystallization occurs. In thermodynamic theory, there is a limit to the degree of metastability that can be achieved called the spinodal (Fig. 1), beyond which the system becomes intrinsically unstable and phase transition is a spontaneous process. For a pure fluid, the thermodynamic spinodal can be defined as the state of diverging isothermal compressibility.11 Upon approaching the spinodal, the nucleation barrier and the critical size are expected to go to zero, and the nucleation time becomes extremely short.1 Therefore, these states are mostly studied indirectly.12–14 

FIG. 1.

Phase diagrams and isotherms for the Lennard-Jones fluid truncated and shifted at 2.5σ, calculated from the PeTS equation of state, in Tρ space (left) and Pρ space (right). The global phase diagram is represented by black curves, namely, the binodal (full black curve), the spinodals (dashed black curves), and the critical point (black dot). The colored curves correspond to isotherms: the green curve corresponds to T = 0.900, the gray corresponds to T = 0.741, and the blue corresponds to T = 0.625. The continuous region on the isotherms corresponds to the metastable vapor and coexisting stable bulk liquid, while the dotted region consists of the metastable liquid and unstable states.

FIG. 1.

Phase diagrams and isotherms for the Lennard-Jones fluid truncated and shifted at 2.5σ, calculated from the PeTS equation of state, in Tρ space (left) and Pρ space (right). The global phase diagram is represented by black curves, namely, the binodal (full black curve), the spinodals (dashed black curves), and the critical point (black dot). The colored curves correspond to isotherms: the green curve corresponds to T = 0.900, the gray corresponds to T = 0.741, and the blue corresponds to T = 0.625. The continuous region on the isotherms corresponds to the metastable vapor and coexisting stable bulk liquid, while the dotted region consists of the metastable liquid and unstable states.

Close modal

150 years after the pioneering work of Gibbs,6 nucleation phenomena keep fascinating scientists. The theoretical developments struggle to keep up with the increasingly accurate results from simulations and experiments that measure the nucleation rates of different substances and types of phase transitions.15–25 Arguably, the main challenge of nucleation theory is to accurately estimate the work of formation for the critical embryo. Since the nucleation rate depends exponentially on this quantity, a small error will cause the predicted nucleation rate to deviate by orders of magnitude from experimental results or simulations.1 

For pure fluids, CNT is qualitatively correct.26 However, it does not account for the expected vanishing of the nucleation barrier at the spinodal, and the predicted rates show systematic deviations from experiments, with errors reaching 20 orders of magnitude for the simple case of argon condensation.27 CNT assumes that the incipient embryo of the new phase is spherical and that its properties are identical to those of the bulk macroscopic stable phase, including the value of the surface tension of the planar vapor–liquid surface. This is referred to as the capillarity approximation, and it has been the main suspect in the problems of CNT.5 

To overcome the limitations of CNT, several routes have been pursued to provide more accurate estimates of the work of formation. Density gradient theory (DGT) assumes that thermodynamic properties depend on density gradients in addition to the local thermodynamic state. It obtains density profiles that minimize a constrained functional for the total Helmholtz energy.28 Its computational simplicity and robustness make DGT attractive. However, it is only the first order approximation to the full classical density functional theory (DFT).29 The accuracy and predictive ability of DGT has, therefore, been questioned, in particular at the highest metastabilities.30 Using full DFT has been touted as a promising alternative.31,32 In previous studies, DFT has been used with great success to estimate a range of important interfacial properties, such as surface tensions33–36 and adsorption isotherms.37 The performance of DFT is limited by a lack of knowledge about the exact functional to extremize34 and knowledge about the real intermolecular potential, even for simple substances like argon38 or water.39 

For many applications, it is impractical and too time-demanding to use DGT or full DFT to calculate the work of formation.5,40 This motivates the development of simple corrections to CNT that maintain the inherent simplicity of the theory. Among the most successful are the curvature corrections of the planar surface tension. The second order curvature expansion of the surface tension is called the Helfrich expansion, where the coefficients of the first- and second-order corrections are given by the Tolman length and rigidity constants, respectively. The challenge of using this approach has been to evaluate the magnitude of the corrections. There have even been controversies about the sign of the leading-order correction, the Tolman length.41–48 

Recently, it was shown how the Tolman length and rigidity constants can be calculated using DGT49–52 or DFT.53,54 Incorporating the curvature corrections into nucleation theory has been shown to significantly improve agreement with experimental results for water.49 It has also solved some of the inconsistencies that arise when CNT is extended to multi-component systems.52 The ultimate test of the accuracy of this curvature-corrected CNT (c-CNT) should be performed using a model system where all the relevant ingredients and properties are as much as possible under control. In this context, the Lennard-Jones (LJ) fluid stands out because it is a simple potential that can be implemented in simulations and has been extensively studied in the literature.55 

In this work, we have performed a comprehensive study of droplets in the Lennard-Jones fluid truncated and shifted at 2.5σ, combining umbrella sampling Monte Carlo simulations, DGT, and DFT to test the accuracy of the curvature-corrected CNT. The focus has been to evaluate the work of formation of liquid droplets for a wide range of conditions from the binodal to the spinodal and at different temperatures. We show that DFT and DGT yield similar results that are in excellent agreement with the molecular simulation results. Remarkably, c-CNT reproduces very accurately the results from simulations, DGT, and DFT for all temperatures and moderate degrees of metastability.

Since the c-CNT is based on the Helfrich expansion, which is a Taylor expansion around the planar surface, it does not give a vanishing nucleation barrier at the spinodal, unless by chance. This hampers its precision at high metastabilities and for sufficiently small embryos. We mend this shortcoming by formulating a scaling expression for the nucleation barrier that incorporates both the coefficients of the Helfrich expansion and the requirement of vanishing nucleation barriers at the spinodals as corrections to CNT. We demonstrate that the modified theory gives an excellent representation of the work of formation that differs from the predictions from DGT by less than one kBT for the entire metastability range and all temperatures examined.

Consider a metastable vapor at temperature T and chemical potential μ. For a spherical critical cluster, the work of formation is given by the formally exact relation,6,56
W=16πγ33(ΔP)2,
(2)
where γ is the surface tension of the cluster with the surface of tension as a dividing surface. The Laplace pressure is defined as
ΔP=P(T,μ)Pv(T,μ),
(3)
i.e., the pressure difference between the (stable) bulk liquid and the (metastable) bulk vapor having the same intensive properties (Tμ). The calculation of the work of formation thus requires the evaluation of two main ingredients: the pressure difference ΔP and the surface tension γ.
ΔP is the difference in bulk thermodynamic properties. It can be related to the surface tension through the Laplace relation,
ΔP=2γR.
(4)
Here, R is the radius of the critical cluster at the surface of tension.51,57 In this work, the quantities R, γ (with no subscript) are always associated with the surface of tension.
For nucleation in an open system, the formation of a critical droplet will increase the number of particles in the volume compared to the homogeneous metastable vapor. One definition of the excess number of particles is58 
ΔN=WμT,
(5)
which is independent of the dividing surface.

For a given metastable phase (Tμ), Eqs. (2)(4) constitute three equations for four properties. The closure of this system of equations is obtained by specifying a relation between γ and R. The closures considered in this work are detailed in Secs. II B–II F.

To calculate ΔP in Eq. (2) for condensation from a metastable vapor, one needs an equation of state (EoS) that is accurate for the stable liquid and the metastable vapor. For the Lennard-Jones fluid truncated and shifted at 2.5σ, LJTS(2.5σ), the PeTS equation of state by Heier et al.59 satisfies this criterion. (Note that an EoS of similar accuracy was also developed for the even-shorter-ranged Lennard-Jones spline potential,60 but for this potential simulations of droplets are missing.) We have verified the accuracy of PeTS in predicting the binodal by performing Gibbs ensemble Monte Carlo simulations61,62 for the three temperatures studied here (see the supplementary material).

The classical nucleation theory (CNT) formula for the nucleation barrier is obtained from Eq. (2) by using the capillarity approximation, which assumes that the surface tension γ is the same as that of the bulk liquid with a planar interface γ0,
WCNT=16πγ033(ΔP)2.
(6)
Accordingly, the ratio of the exact work of formation of a spherical critical droplet to that of the CNT-approximation is
WWCNT=γγ03.
(7)
Consequently, W/WCNT → 1 in the planar limit (ΔP → 0), and this is often interpreted as the work of formation from CNT being exact in the limit of large droplets. At the same time,
WWCNT=16π3(ΔP)2(γ3γ03),
(8)
=16π3(ΔP)2((γ0+αΔP+)3γ03),
(9)
=16π33γ02αΔP+O(1),
(10)
where α = (∂γ/ΔP)T. Therefore, if the surface tension has a non-zero derivative with respect to ΔP in the planar limit, the absolute error of CNT for the work of formation can be substantial for large droplets. More generally, the difference in work of formation from two theories will diverge upon approaching the binodal if the two theories predict differing values of γ0 or ∂γ/P) at the binodal (or both). This must be kept in mind when comparing theories close to the binodal. But obviously, the relative error (WWCNT)/WCNT goes to zero at the binodal, as expected for the thermodynamic limit.
Density gradient theory (DGT) models the grand free energy at the state defined by (Tμ, V) by the following functional of the density profile:
Ω[ρ(r)]=κ2(ρ(r))2+aEoS(ρ(r))μρ(r)dr.
(11)
The influence parameter κ is a positive parameter, usually fitted to reproduce the planar surface tension, and aEoS(ρ) is the Helmholtz energy density (per volume) of a fluid with uniform density ρ at the specified temperature. The function aEoS(ρ) = −PEoS(ρ) + μEoS(ρ)ρ is calculated from the PeTS EoS. We used the constant influence parameter κ = 2.7334 (in reduced, LJ-units) suggested by Heier et al.,59 which was shown to reproduce planar surface tension simulations almost within their error bars.

The grand canonical functional is stationary for a system at a fixed temperature T, volume, and chemical potential μ. The necessary conditions for a stationary point of Eq. (11) are provided by the Euler–Lagrange equations, δΩρr/δρr=0. To solve the Euler–Lagrange equations for DGT, we used the approach described in Ref. 51.

Heier et al.59 derived a density functional theory (DFT) for the LJTS(2.5σ) fluid based on a weighted density average (WDA) approach.34 To arrive at a formulation that is consistent with Barker–Henderson perturbation theory63 as represented by the PeTS EoS, the functional used the White Bear fundamental measure theory (FMT) for the hard sphere contribution.64 Density functional theory models the grand free energy at the state defined by (Tμ, V) by the following functional of the density profile:
Ωρ(r)=Aρ(r)ρrμdr,
(12)
where Aρ(r) is the Helmholtz energy functional. For the LJTS(2.5σ) fluid, the Helmholtz energy has contributions from the hard-sphere reference (FMT), dispersion (disp), and an ideal gas (id) term,
βAρr=ΦFMTnFMTr+Φdispndispr+Φidρrdr,
(13)
where Φ is the reduced Helmholtz energy density, and nr is a weighted density. The weighted densities are calculated via convolutions with the density profiles. For the contribution from FMT, the weights are listed in Ref. 65. The weighted density for the dispersion term is given by
ndispr=drρrωdisprr,
(14)
where ωdisp is a weight function of a sphere with diameter φdBH,
ωdispr=ΘφdBH|r|4π3(φdBH)3.
(15)
Here, Θ is the Heaviside function, dBH is the temperature dependent Barker–Henderson diameter,63 and φ is a parameter that depends on the potential range. Using the WDA approach, the dispersion contribution is a function of the dispersion term from the EoS as follows:
Φdispndispr=ndisprFEoS,dispndispr,
(16)
where F is the specific reduced Helmholtz energy.

The profiles from DFT can be obtained by solving the Euler–Lagrange equations, δAρr/δρr=μ. This can be efficiently accomplished by using the convolution theorem to solve the convolution integrals in Fourier space with cosine and sine transforms, for both the planar and spherical geometries.65,66 In addition, convergence of the overall equation system can be accelerated by use of Anderson mixing.67 

Sauer and Gross34 regressed the φ parameter to surface tension data of n-alkanes and found that φ = 1.3862 gave good results for most fluids. For the short ranged LJTS(2.5σ) fluid, Heier et al. had to use a smaller value of φ = 1.21 to achieve agreement with the simulation results for the surface tension;68 this is also the value used in this work.

For both DFT and DGT, the bulk behavior is governed by the underlying EoS, i.e., PeTS. Since the work of formation of a droplet (bubble) in DGT will be zero at the vapor (liquid) spinodal of the EoS, the location of this spinodal may severely affect the prediction of free energies of small clusters. Although non-local DFT yields more realistic density profiles than DGT,55 they have been found to yield similar free energy predictions54 in the planar limit. In this work, we will also compare them for small droplets.

Of interest to the current work is the behavior of mean field theories, such as DGT and DFT, close to the spinodal. Wilemski and Li30 studied in detail the scaling behavior of density gradient theory and mean field density functional theory. Both of these theories always predict that the work of formation vanishes at the spinodal, whereas the excess number of particles [Eq. (5)] of the critical cluster diverges to infinity. The mean field theories predict that the equimolar radius Re and the equimolar surface tension γe both diverge to infinity upon approaching the spinodal. For the surface of tension, however, γ → 0, R → 0 at the spinodal, whereas ΔP goes to a finite value.

The Helfrich expansion is based on a quadratic expansion of the surface tension in terms of the curvature 1/R around the planar limit,
γ(R)=γ02γ0δ/R+ks/R2,
(17)
where δ is the Tolman length, and ks=2k+k̄ is the spherical rigidity, with k and k̄ being the bending and Gaussian rigidity, respectively. All these coefficients can be computed using either DGT or DFT.50,51,54 Such an expansion can be performed for an arbitrary choice of dividing surface, but in this work, we use the surface of tension.
The Helfrich expansion can be viewed as a second-order refinement of the zeroth-order capillarity approximation γ(R) = γ0. Incorporating this expansion into CNT, one obtains a curvature corrected expression for the work of formation of critical clusters,52,
Wc-CNT4πγ0RCNT2316δRCNT+4πks,
(18)
where
RCNT=2γ0ΔP,
(19)
is the CNT prediction for the radius of the critical cluster. The physical significance of this curvature corrected CNT is corroborated by its success in improving the predictions of classical nucleation theory for the experimental nucleation rates of water50 and water-alcohol mixtures.52 In this work, we shall perform a more stringent test of the accuracy of this formula for the LJTS(2.5σ) using the values of the coefficients calculated from DGT according to the procedure of Ref. 51.

Since the pioneering work of McGraw and Laaksonen on scaling relations,26 there is increasing evidence suggesting that the work of formation of critical clusters normalized by the CNT prediction, fW/WCNT, may be a universal function of the degree of metastability, depending only weakly on the temperature and on the details of the intermolecular potential.69–73 Accordingly, many expressions for this scaling function f have been proposed in the literature,69,72,74–79 where the degree of metastability is typically expressed either in terms of the difference in chemical potential Δμ or the Laplace pressure ΔP. This scaling function should fulfill at least two requirements. First, CNT is expected to be valid at the binodal when the size of the critical cluster goes to infinity, such that f(0) = 1. Second, the work of formation of the critical cluster should be zero at the spinodal, such that fPsp) = 0, where ΔPsp represents the difference in bulk pressures between the stable and metastable phases at the spinodal.

The most straightforward scaling function fulfilling these two requirements is the linear function fP) = 1 − ΔPPsp, proposed (in terms of Δμ) by Talanquer.69 Along the same lines, different linear and parabolic scaling functions have been proposed in the literature.72,74–79 In fact, the c-CNT expression, Eq. (18), can be rewritten as
fc-CNT(ΔP)Wc-CNTWCNT=13δγ0ΔP+3ks4γ03ΔP2.
(20)
However, c-CNT is derived using only the Helfrich expansion, and although it predicts a vanishing of the nucleation barrier, this does not necessarily agree with the location of the thermodynamic spinodal (i.e., at ΔPsp, where the isothermal compressibility diverges).

A very similar parabolic scaling function was recently proposed by Kashchiev79 based on the nucleation theorem and using the assumption that the equimolar critical cluster radius only differs from the CNT prediction by the Tolman length. The main difference is that the second order term is obtained by imposing that the free energy of the critical cluster must vanish at the vapor spinodal, instead of using the value of the spherical rigidity. In fact, Kaschiev’s expression can be recovered assuming that the scaling function is parabolic and using f(0) = 1, fPsp)= 0, and the value of Tolman length δ as conditions to evaluate the coefficients.

We can formulate a more accurate expression for the scaling function fP), using all the ingredients that we know, namely, (i) f(0) = 1, (ii) the slope at the binodal given by the Tolman length δ, (iii) the spherical rigidity ks, (iv) zero work of formation at the vapor spinodal f(ΔPspv)=0, and (v) at the liquid spinodal f(ΔPsp)=0. Instead of using ks to specify the second derivative at the binodal, it is more convenient to use it to enforce the location of the maximum of the work of formation. We have enforced that the maximum of f occurs at
ΔPm=2δγ02ks,
(21)
namely, the same maximum as obtained from Eq. (20). An accurate scaling function that satisfies these five conditions is a fourth-order polynomial in the Laplace pressure,
f(ΔP)=(ΔPΔPsp)(ΔPΔPspv)(1+bΔP+c(ΔP)2)ΔPspΔPspv,
(22)
b=3δγ0+1ΔPspv+1ΔPsp,
(23)
c=(3bΔPm+2)ΔPm+(ΔPsp+ΔPspv)(2bΔPm+1)ΔPspΔPspvb4(ΔPm)33(ΔPsp+ΔPspv)(ΔPm)2+2ΔPspΔPspvΔPm.
(24)

The first two factors in the numerator of Eq. (22) warrant that the work of formation vanishes at both spinodals. The coefficients of the parabolic third factor are obtained from the slope at the binodal given by the Tolman length and the location of the maximum of the work of formation at ΔPm.

In Sec. IV, we will show that Eq. (22) accurately reproduces results from simulations, DGT, and DFT. The proposed scaling function in Eq. (22) is a modification of the Helfrich expansion so that the surface tension goes to zero at the right location of the spinodals while maintaining the correct behavior in the binodal limit described by the Tolman length and the spherical rigidity.

In order to test the accuracy of c-CNT and of the proposed scaling relation, we have performed umbrella sampling Monte Carlo simulations to calculate the free energy of critical droplets of the LJTS(2.5σ) fluid at two different temperatures: T = 0.625 and T = 0.900. Monte Carlo simulations in the NPT ensemble were implemented in a cubic box with periodic boundary conditions, using a self-developed code based on the algorithms described in Refs. 80 and 81. Simulations were performed with typically N = 4000 particles for the largest supersaturations and N = 10 000 for small supersaturations. In the MC simulations, each trial move consisted either of an attempted particle displacement or a trial volume change, chosen at random. The maximum displacement for both types of MC moves was adjusted to give a 50% acceptance ratio. A MC step is defined as N + 1 trial moves.

Liquid clusters were identified using the ten Wolde–Frenkel cluster criterion,81 where a molecule is classified as liquid-like if it has at least five neighbors within a distance smaller than the Stillinger radius rc, corresponding to the first minimum in the radial distribution function of the liquid at the desired temperature.

The number of liquid-like molecules in the largest liquid cluster n was used as a reaction coordinate in the harmonic umbrella biasing potential,
Wb(n)=12kn(nn0)2,
(25)
where kn and n0 set the width and location of the window. A total of 15–40 windows located at intervals of 20–150 molecules with typical values of kn = 0.005 − 0.02 were simulated at each temperature and pressure. For each window, 500k MC steps were used for equilibration, followed by 500k MC steps used for production. To speed up the simulations, a staging scheme was implemented, where a short trajectory of 25 MC steps was performed without bias and then accepted with a probability given by expβΔWb(n), where ΔWb(n) is the difference in biasing potential before and after the trajectory. At each temperature and pressure, we also performed an unbiased simulation to obtain the cluster size distribution P(n) of all clusters. In each window, the probability distribution function of the largest cluster size P(n) was recorded. The weighted histogram analysis method described in Ref. 82 was then used to reconstruct the free energy from the umbrella sampling simulations, joining the different windows. In order to obtain the free energy of any cluster size, βΔG(n), from the probability distribution of the largest cluster P(n), we used the formula
βΔG(n)=lnP(n)P(nmin)lnP(nmin),
(26)
where nmin is the size that we used to match the two distributions (typically nmin = 5). Figure 5 in the supplementary material shows a typical example of the reconstructed free energy landscape for T = 0.900 and p = 0.037. The location and value of the maximum of the reconstructed free energy βΔG(n) provide the size ΔN and energy of formation of the critical cluster W.

For T = 0.900, simulations were performed in the pressure range of 0.034 25 ≤ P ≤ 0.039 and for T = 0.625 in the range of 0.0028 ≤ P ≤ 0.078. For pressures smaller than that, the critical cluster size was too large to be simulated. For pressures larger than those values, the system was too close to the spinodal, and many clusters were formed spontaneously, hampering a reliable reconstruction of the barrier. Therefore, the range of nucleation barriers accessible through our umbrella sampling simulations is in the range of 20–150 kBT. Details from the simulations, such as the work of formation and the number of particles, are included in the supplementary material.

Simulation results for a third intermediate temperature T = 0.741 were also included in our comparison. For that temperature, we took advantage of the results of ten Wolde and Frenkel,81 who also used umbrella sampling to calculate the free energy at P = 0.007 83, and then used thermodynamic integration to obtain the nucleation barrier at other pressures.

We next discuss the results. The numerical values of the simulation results are listed in the supplementary material. All properties are reported in reduced Lennard-Jones units. Conversion formulas for real fluids are listed in Table I of Ref. 83.

Figures 24 compare the number of particles and the work of formation predicted from all theories considered in this work against simulations for the temperatures T = 0.625, T = 0.741, and T = 0.900. The x axes start from the binodal and end at the vapor spinodal for each temperature.

FIG. 2.

Comparison of theories and simulations for T = 0.625 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. The simulation uncertainties, estimated following Ref. 82, are shown with error bars. Helfrich coefficients were computed with DGT.

FIG. 2.

Comparison of theories and simulations for T = 0.625 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. The simulation uncertainties, estimated following Ref. 82, are shown with error bars. Helfrich coefficients were computed with DGT.

Close modal
FIG. 3.

Comparison of theories with simulation data from ten Wolde and Frenkel81 for T = 0.741 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. Helfrich coefficients were computed with DGT.

FIG. 3.

Comparison of theories with simulation data from ten Wolde and Frenkel81 for T = 0.741 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. Helfrich coefficients were computed with DGT.

Close modal
FIG. 4.

Comparison of theories and simulations for T = 0.900 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. The simulation uncertainties, estimated following Ref. 82, are shown with error bars. Helfrich coefficients were computed with DGT.

FIG. 4.

Comparison of theories and simulations for T = 0.900 as a function of the Laplace pressure ΔP, where the x axes are delimited by the binodal on the left and the vapor spinodal on the right. The simulation uncertainties, estimated following Ref. 82, are shown with error bars. Helfrich coefficients were computed with DGT.

Close modal

All theories, namely, DGT, DFT, CNT, and c-CNT (labeled “Helfrich expansion” in Figs. 24), and the scaling function, reproduce the excess particle number in the critical cluster ΔN obtained in the simulations, mostly within the simulation uncertainty. Although the excess number of particles is closely related to the work of formation through Eq. (5), in the simulations, it is calculated directly from the number of molecules in the critical cluster identified by the ten Wolde-Frenkel cluster criterion. Therefore, the fact that a theory can reproduce both the work of formation and the excess number of particles represents a stringent validation as well as an excellent consistency check of the simulations. Mean field theories, such as DFT and DGT, predict a diverging excess number of particles upon approaching the spinodal.30 This is not shown in Figs. 24 [panel (a)] since this divergence only occurs very close to the spinodal. However, note that none of the theories predict zero excess particle number at the spinodal; this is most evident from Fig. 3, panel (a).

For the work of formation, DGT and the proposed spinodal scaling relation are in closest agreement with the molecular simulation results. This can be seen by comparing the work of formation ratio W/WCNT and the value of WWDGT [Figs. 24, panels (c) and (d)] to the simulations.

The excellent match of the simulations and DGT with a temperature-independent influence parameter is remarkable, considering that the uncertainty of the planar surface tension values for the LJTS(2.5σ) fluid is on the order of several percent (cf. the supplementary material). Intuitively, one might expect DFT to give results that are closer to simulations than DGT. After all, DGT is only the first approximation to DFT. It should be kept in mind, however, that the exact functional to minimize in DFT is unknown, and the functional presented in Sec. II D represents only an approximation. A fitting parameter is needed in both theories to reproduce surface tensions from molecular simulations.

Stephan et al.55 showed that although the surface tension of the planar vapor-interface is accurately represented by both DFT and DGT, the theories give density and pressure tensor profiles that deviate from molecular simulation results. One theory did not show any clear advantage over the other (see Figs. 2–6 in Ref. 55). These discrepancies do not seem to influence the curvature dependence of the surface tension to a large degree, since both theories reproduce the work of formation from umbrella sampling simulations with impressive accuracy.

The c-CNT that uses the Helfrich expansion is also remarkably accurate when measuring the deviation in units of kBT for the droplet sizes studied in this work [Figs. 24, panel (d)]. Still, the deviation from simulations is only sligthly outside the simulation error bars at the largest metastabilities, as seen in Fig. 2. Only the CNT approximation results in large deviations from the simulations.

With the exception of CNT, all theories predict that the surface tension increases with metastability at low, positive values of ΔP, reaches a maximum, and then decreases monotonically toward the spinodal. In practice, homogeneous nucleation only occurs for metastabilities where W < 100kBT; for higher barriers, the nucleation rate is generally too low to be observable. For nearly all droplets for which W < 100kBT, the theories predict that the surface tension is lower than the planar surface tension for the LJTS(2.5σ) fluid.

The theories that only use the properties of the surface tension associated with the binodal, namely, CNT and the Helfrich expansion, differ markedly at high supersaturations from the theories that incorporate the spinodal. CNT does not predict a spinodal at all, whereas the Helfrich expansion predicts spinodals at significantly higher metastabilities than the PeTS EoS. Both theories, therefore, predict a positive work of formation W at the PeTS spinodal, which is largest for the lowest temperatures [compare Figs. 24 (panel d) for the three temperatures]. For CNT, this overprediction persists to the metastabilities of the simulations, showing clearly that CNT overpredicts the free energy of small droplets and that this is most severe at low temperatures and/or high metastabilities.

Figures 24 [panel (c)] illustrate that DFT and DGT yield very similar shapes of the ratio W/WCNT=(γ/γ0)3, where γ0 was calculated from the corresponding theory. For the value of W [Figs. 24 panel (d)], DGT performs slightly better than DFT in reproducing the simulation results. One possible reason is the different values of γ0 from the two theories (cf. Table I), as well as the different Helfrich coefficients δ,k,k̄. That the curves WCNTWDGT in panel (d) diverge to infinity at the binodal follows from Eqs. (8)(10), since the Tolman length of DGT is negative. We also see that WDFTWDGT diverges when approaching the binodal, where the divergence is to ∞ or −∞ depending on the sign of γ0DFTγ0DGT.

TABLE I.

Planar surface tension and Helfrich coefficients for the LJTS(2.5σ) fluid.

γ0δkk̄2k+k̄
T = 0.625 
DGT 0.732 11 −0.081 26 −0.322 32 0.192 49 −0.452 16 
DFT 0.715 81 −0.072 81 −0.323 49 0.200 50 −0.446 47 
T = 0.700 
DGT 0.585 99 −0.078 29 −0.323 00 0.189 84 −0.456 16 
DFT 0.580 60 −0.070 32 −0.336 53 0.204 46 −0.468 59 
T = 0.741 
DGT 0.509 38 −0.076 35 −0.323 10 0.188 45 −0.457 75 
DFT 0.507 58 −0.068 40 −0.341 72 0.205 52 −0.477 92 
T = 0.800 
DGT 0.403 18 −0.073 18 −0.322 17 0.186 09 −0.458 24 
DFT 0.404 30 −0.065 01 −0.345 78 0.205 16 −0.486 39 
T = 0.900 
DGT 0.235 23 −0.067 36 −0.313 89 0.178 93 −0.448 85 
DFT 0.237 23 −0.058 28 −0.338 83 0.197 00 −0.480 67 
T = 1.000 
DGT 0.088 27 −0.064 58 −0.281 11 0.158 62 −0.403 59 
DFT 0.088 89 −0.054 40 −0.296 58 0.169 26 −0.423 90 
γ0δkk̄2k+k̄
T = 0.625 
DGT 0.732 11 −0.081 26 −0.322 32 0.192 49 −0.452 16 
DFT 0.715 81 −0.072 81 −0.323 49 0.200 50 −0.446 47 
T = 0.700 
DGT 0.585 99 −0.078 29 −0.323 00 0.189 84 −0.456 16 
DFT 0.580 60 −0.070 32 −0.336 53 0.204 46 −0.468 59 
T = 0.741 
DGT 0.509 38 −0.076 35 −0.323 10 0.188 45 −0.457 75 
DFT 0.507 58 −0.068 40 −0.341 72 0.205 52 −0.477 92 
T = 0.800 
DGT 0.403 18 −0.073 18 −0.322 17 0.186 09 −0.458 24 
DFT 0.404 30 −0.065 01 −0.345 78 0.205 16 −0.486 39 
T = 0.900 
DGT 0.235 23 −0.067 36 −0.313 89 0.178 93 −0.448 85 
DFT 0.237 23 −0.058 28 −0.338 83 0.197 00 −0.480 67 
T = 1.000 
DGT 0.088 27 −0.064 58 −0.281 11 0.158 62 −0.403 59 
DFT 0.088 89 −0.054 40 −0.296 58 0.169 26 −0.423 90 

Although we have not compared to molecular simulation results for bubbles, we have included a comparison between the theories and methods used in this work for bubbles in the supplementary material. This comparison shows that the performance of the theories discussed above applies to the bubble branch as well. Remarkably, the proposed scaling function Eq. (22) reproduces almost exactly the predictions of DGT both for bubbles and droplets in the entire range of metastabilities. The supplementary material also includes versions of Figs. 24, where the quantities are plotted as functions of the supersaturation ratio P/Psat, which is commonly used as a measure of metastability in nucleation theory.

We will next compare our estimates of the Tolman length to those of van Giessen and Blokhuis,84 who determined the Tolman length δ of the LJTS(2.5σ) fluid at T = 0.900 from bulk pressures of large liquid drops via molecular dynamics simulations. Using the planar surface tension γ0 as a boundary condition in their analysis, they found δ = −0.100σ or δ = −0.0845σ depending on how they calculated the bulk pressures. This deviates somewhat from the value determined by DGT and DFT (cf. Table I). In fact, the values of the Tolman length we find, −0.067 36 from DGT and −0.058 28 from DFT, are closer to those estimated by Wilhelmsen et al.,50 namely, δ = −0.074σ. Wilhelmsen et al. used another implementation of DGT with a less accurate equation of state.50 This shows that the coefficients in the Helfrich expansion are sensitive toward the underlying thermodynamic description of the bulk fluid.

For T = 0.625, the simulation results for the largest droplets yield values of the surface tension that are higher than the planar value (Fig. 2), which is consistent with a negative Tolman length. To our knowledge, this is the first time this has been shown by umbrella sampling simulations. The regime with surface tensions that are higher than the surface tension of the planar surface was only reachable at the lowest temperatures. This is because the critical cluster size corresponding to the maximum and the required number of particles in the simulation volume increase quickly with temperature.

Interestingly, the value of the spherical rigidity ks=2k+k̄ appears to be approximately independent of temperature for T ≤ 0.9, with ks ∼ −0.45 for DGT (cf. Table I). This is in line with the findings by Wilhelmsen et al.50 How the W/WCNT curve varies with temperatures is therefore strongly connected to the value of the Tolman length.

For T = 0.900, van Giessen and Blokhuis84 calculated the planar surface tension of the LJTS(2.5σ) fluid for increasingly large systems (up to box sizes of 50σ × 50σ × 100σ) and extrapolated to the thermodynamic limit. This yielded the value γ0=0.227774, which is 3.2% lower than the value from DGT, γ0DGT=0.23523, and 4% lower than the value from DFT, γ0DFT=0.23723. We have attempted to refit the influence parameter of DGT to the presumably more accurate value from van Giessen and Blokhuis84 and found that this removes the good agreement between DGT and the simulations shown in Fig. 4 (cf. Fig. 4 in the supplementary material). One possible reason for this is the larger influence of capillary waves on large surfaces. Using the formulation of Planková et al.,85 capillary waves are expected to decrease the surface tension γ according to
γbareγ0=1+38πTTc12.72,
(27)
where Tc is the critical temperature, and Tc = 1.085 for LJTS(2.5σ). The capillary wave contribution to the surface tension is the largest close to the critical point, with a limiting value of γbare/γ0=1.044 at the critical temperature. For T = 0.900, this ratio equals 1.0364, which is in reasonable agreement with γ0DGT/γ0=1.0327. This may support the conclusions of Planková et al.,85 namely, that the small length scales of critical droplets tend to inhibit capillary waves, and it is therefore the bare surface tension that is relevant. It also seems to support the findings of Checa et al.,86 which indicate that DFT only captures fluctuations at length scales of (10 ± 2)σ and that large capillary waves must be separately accounted for.

In this work, we have performed a careful analysis of the work of formation of critical clusters of the truncated and shifted Lennard-Jones fluid LJTS(2.5σ) at three different temperatures and a wide range of pressures from the binodal to the spinodal. We have performed Monte Carlo umbrella sampling simulations and used density gradient theory (DGT) and density functional theory (DFT) to study the work of formation and the size of critical clusters. The aim was to test the accuracy of classical nucleation theory (CNT) and curvature corrected CNT (c-CNT). c-CNT relies on the Helfrich expansion, which is a second order expansion of the surface tension in curvature around the planar limit, where the first order correction is the Tolman length and the second order correction is the rigidity constants.

We have found that both DGT and DFT reproduce the simulation results with remarkable accuracy. As expected, since CNT uses the capillarity approximation, it overpredicts the nucleation barrier at high supersaturations and fails to predict a vanishing nucleation barrier at the spinodal. Remarkably, c-CNT reproduces accurately the results from molecular simulations at all temperatures and for most of the experimentally accessible region. This provides a solid confirmation of the validity of c-CNT, using a model system like the LJTS(2.5σ), with a simple intermolecular potential, and for which the equation of state and relevant properties have been characterized accurately by simulations.

The Tolman lengths from both DGT and DFT were found to be negative. Moreover, the second-order correction of the surface tension, the spherical rigidity, was found to be nearly independent of temperature. At the lowest temperature investigated, the work of formation from umbrella sampling simulations revealed a region where the surface tension is higher than the value of the planar vapor–liquid surface. This strongly suggests that the Tolman length is negative, in line with the predictions from DGT and DFT.

The molecular simulations for T = 0.900 were shown to be in agreement with capillary wave theory, suggesting that the capillary wave contributions to the planar surface tension must be removed before using it to tune DGT/DFT when modeling droplets. This warrants further investigation of the size-dependence of the surface tension for other temperatures, e.g., T = 0.625 and T = 0.741, and for other fluids to investigate whether this finding is general.

The only regime where c-CNT deviates from simulations is for the smallest droplets at the highest metastabilities in the vicinity of the spinodal. In fact, c-CNT mispredicts the location of the spinodal, identified as the vanishing of the nucleation barrier. This is expected, since c-CNT is based on a curvature expansion around the planar limit, i.e., at the binodal.

In order to solve this limitation, we have formulated a new scaling function that, in addition to the leading order curvature corrections (i.e., Tolman length and rigidity constants), also uses the location of the vapor and liquid spinodals. The resulting fourth order polynomial reproduces with excellent accuracy the results from simulations, DGT and DFT, in the whole range of studied temperatures and degrees of metastability, from the binodal to the spinodal. It also reproduces, with remarkable accuracy, the DGT and DFT predictions for the bubble branch.

Our results suggest that both c-CNT and the new scaling function constitute accurate yet simple models to accurately predict nucleation rates and solve the limitations of CNT. Here, we have confirmed this for the work of formation of droplets. However, the new scaling function is expected to be equally accurate for bubbles and could also possibly be extended to crystallization. Therefore, the proper incorporation of curvature corrections and the location of the spinodals might open the door to quantitatively accurate predictions of nucleation rates for other substances and a wide range of problems of scientific and technological interest.

The supplementary material contains a comparison of the predictions of the different theories on the bubble branch; a comparison between simulations and theories for critical droplets for T = 0.900 when refitting the influence parameter of DGT to reproduce the planar surface tension from van Giessen and Blokhuis;84 a table with all the simulation results; a validation using Gibbs ensemble Monte Carlo simulations of the PeTS EoS for the three temperatures studied in this work; a comparison of the calculated surface tensions from DFT and DGT against the simulation results by Vrabec et al.;68 and an alternative representation of Figs. 2–4 in terms of the supersaturation.

This work received funding from the Norwegian Research Council, Project No. 328679. A.A., Ø.W., and M.H. acknowledge funding from the Research Council of Norway (RCN), the Center of Excellence Funding Scheme, Project No. 262644, PoreLab. D.R. acknowledges funding from the Spanish government through Grant Nos. PGC2018-098373-B-I00 and PID2021-126570NB-I00 (MINECO/FEDER, UE).

The authors have no conflicts to disclose.

Ailo Aasen: Conceptualization (equal); Data curation (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Øivind Wilhelmsen: Conceptualization (equal); Investigation (equal); Supervision (equal); Writing – review & editing (equal). Morten Hammer: Data curation (equal); Investigation (equal); Writing – review & editing (equal). David Reguera: Conceptualization (equal); Data curation (equal); Investigation (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.

1.
P.
Debenedetti
,
Metastable Liquids: Concepts and Principles
(
Princeton University Press
,
Princeton
,
1996
).
2.
A.
Laaksonen
,
V.
Talanquer
, and
D. W.
Oxtoby
,
Annu. Rev. Phys. Chem.
46
,
489
(
1995
).
3.
D.
Kashchiev
,
Nucleation: Basic Theory with Applications
(
Butterworth-Heinemann
,
Oxford
,
2000
).
4.
K. F.
Kelton
and
A. L.
Greer
,
Nucleation in Condensed Matter
(
Elsevier
,
Amsterdam
,
2010
).
5.
H.
Vehkamäki
,
Classical Nucleation Theory in Multicomponent Systems
(
Springer-Verlag
,
Berlin
,
2006
).
6.
J. W.
Gibbs
,
The Scientific Papers of J. Willard Gibbs
(
Ox Bow Press
,
London
,
1993
).
7.
M.
Volmer
and
A.
Weber
,
Z. Phys. Chem.
119U
,
277
(
1926
).
8.
9.
R.
Becker
and
W.
Döring
,
Ann. Phys.
416
,
719
(
1935
).
10.
Y. B.
Zeldovich
,
Zh. Teor. Eksp. Fiz.
12
,
525
(
1942
).
11.
P.
Aursand
,
M. A.
Gjennestad
,
E.
Aursand
,
M.
Hammer
, and
Ø.
Wilhelmsen
,
Fluid Phase Equilib.
436
,
98
(
2017
).
12.
A.
Linhart
,
C. C.
Chen
,
J.
Vrabec
, and
H.
Hasse
,
J. Chem. Phys.
122
,
144506
(
2005
).
13.
E. S.
Loscar
,
E. E.
Ferrero
,
T. S.
Grigera
, and
S. A.
Cannas
,
J. Chem. Phys.
131
,
024120
(
2009
).
14.
E. S.
Loscar
,
C. G.
Ferrara
, and
T. S.
Grigera
,
J. Chem. Phys.
144
,
134501
(
2016
).
15.
E. M.
Dunne
,
H.
Gordon
,
A.
Kürten
,
J.
Almeida
,
J.
Duplissy
,
C.
Williamson
,
I. K.
Ortega
,
K. J.
Pringle
,
A.
Adamov
,
U.
Baltensperger
,
P.
Barmet
,
F.
Benduhn
,
F.
Bianchi
,
M.
Breitenlechner
,
A.
Clarke
,
J.
Curtius
,
J.
Dommen
,
N. M.
Donahue
,
S.
Ehrhart
,
R. C.
Flagan
,
A.
Franchin
,
R.
Guida
,
J.
Hakala
,
A.
Hansel
,
M.
Heinritzi
,
T.
Jokinen
,
J.
Kangasluoma
,
J.
Kirkby
,
M.
Kulmala
,
A.
Kupc
,
M. J.
Lawler
,
K.
Lehtipalo
,
V.
Makhmutov
,
G.
Mann
,
S.
Mathot
,
J.
Merikanto
,
P.
Miettinen
,
A.
Nenes
,
A.
Onnela
,
A.
Rap
,
C. L. S.
Reddington
,
F.
Riccobono
,
N. A. D.
Richards
,
M. P.
Rissanen
,
L.
Rondo
,
N.
Sarnela
,
S.
Schobesberger
,
K.
Sengupta
,
M.
Simon
,
M.
Sipilä
,
J. N.
Smith
,
Y.
Stozkhov
,
A.
Tomé
,
J.
Tröstl
,
P. E.
Wagner
,
D.
Wimmer
,
P. M.
Winkler
,
D. R.
Worsnop
, and
K. S.
Carslaw
,
Science
354
,
1119
(
2016
).
16.
G. C.
Sosso
,
J.
Chen
,
S. J.
Cox
,
M.
Fitzner
,
P.
Pedevilla
,
A.
Zen
, and
A.
Michaelides
,
Chem. Rev.
116
,
7078
(
2016
).
17.
T.
Karmakar
,
P. M.
Piaggi
, and
M.
Parrinello
,
J. Chem. Theory Comput.
15
,
6923
(
2019
).
18.
G.
Menzl
,
M. A.
Gonzalez
,
P.
Geiger
,
F.
Caupin
,
J. L. F.
Abascal
,
C.
Valeriani
, and
C.
Dellago
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
13582
(
2016
).
19.
C.
Duan
,
R.
Karnik
,
M.-C.
Lu
, and
A.
Majumdar
,
Proc. Natl. Acad. Sci. U. S. A.
109
,
3688
(
2012
).
20.
P. M.
Winkler
,
G.
Steiner
,
A.
Vrtala
,
H.
Vehkamäki
,
M.
Noppel
,
K. E. J.
Lehtinen
,
G. P.
Reischl
,
P. E.
Wagner
, and
M.
Kulmala
,
Science
319
,
1374
(
2008
).
21.
J.
Zhou
,
Y.
Yang
,
Y.
Yang
,
D. S.
Kim
,
A.
Yuan
,
X.
Tian
,
C.
Ophus
,
F.
Sun
,
A. K.
Schmid
,
M.
Nathanson
,
H.
Heinz
,
Q.
An
,
H.
Zeng
,
P.
Ercius
, and
J.
Miao
,
Nature
570
,
500
(
2019
).
22.
A. E. S.
Van Driessche
,
N.
Van Gerven
,
P. H. H.
Bomans
,
R. R. M.
Joosten
,
H.
Friedrich
,
D.
Gil-Carton
,
N. A. J. M.
Sommerdijk
, and
M.
Sleutel
,
Nature
556
,
89
(
2018
).
23.
E. B.
Moore
and
V.
Molinero
,
Nature
479
,
506
(
2011
).
24.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
,
J. Chem. Phys.
144
,
034501
(
2016
).
25.
K. E.
Blow
,
D.
Quigley
, and
G. C.
Sosso
,
J. Chem. Phys.
155
,
040901
(
2021
).
26.
R.
McGraw
and
A.
Laaksonen
,
Phys. Rev. Lett.
76
,
2754
(
1996
).
27.
K.
Iland
,
J.
Wölk
,
R.
Strey
, and
D.
Kashchiev
,
J. Chem. Phys.
127
,
154506
(
2007
).
28.
E.
Magnanelli
,
Ø.
Wilhelmsen
,
D.
Bedeaux
, and
S.
Kjelstrup
,
Phys. Rev. E
90
,
032402
(
2014
).
29.
R.
Evans
, in
Fundamentals of Inhomogeneous Fluids
, edited by
D.
Henderson
(
Marcel Dekker
,
New York
,
1992
).
30.
G.
Wilemski
and
J.-S.
Li
,
J. Chem. Phys.
121
,
7821
(
2004
).
31.
D. W.
Oxtoby
and
R.
Evans
,
J. Chem. Phys.
89
,
7521
(
1988
).
33.
X. C.
Zeng
and
D. W.
Oxtoby
,
J. Chem. Phys.
94
,
4472
(
1991
).
34.
E.
Sauer
and
J.
Gross
,
Ind. Eng. Chem. Res.
56
,
4119
(
2017
).
35.
J. F.
Lutsko
,
J. Chem. Phys.
134
,
164501
(
2011
).
36.
J. C.
Barrett
,
J. Chem. Phys.
124
,
144705
(
2006
).
38.
B.
Jäger
,
R.
Hellmann
,
E.
Bich
, and
E.
Vogel
,
J. Chem. Phys.
135
,
084308
(
2011
).
39.
P.
Rehner
and
J.
Gross
,
J. Chem. Eng. Data
65
,
5698
(
2020
).
40.
Ø.
Wilhelmsen
and
A.
Aasen
,
Chem. Eng. Sci.
248
,
117176
(
2022
).
41.
Y. A.
Lei
,
T.
Bykov
,
S.
Yoo
, and
X. C.
Zeng
,
J. Am. Chem. Soc.
127
,
15346
(
2005
).
42.
I.
Sanchez-Burgos
,
P. M.
De Hijes
,
P.
Rosales-Pelaez
,
C.
Vega
, and
E.
Sanz
,
Phys. Rev. E
102
,
062609
(
2020
).
43.
N.
Bruot
and
F.
Caupin
,
Phys. Rev. Lett.
116
,
056102
(
2016
).
44.
A.
Tröster
,
M.
Oettel
,
B.
Block
,
P.
Virnau
, and
K.
Binder
,
J. Chem. Phys.
136
,
064709
(
2012
).
45.
K. M.
Bal
and
E. C.
Neyts
,
J. Chem. Phys.
157
,
184113
(
2022
).
46.
M. N.
Joswiak
,
R.
Do
,
M. F.
Doherty
, and
B.
Peters
,
J. Chem. Phys.
145
,
204703
(
2016
).
47.
P.
Rosales-Pelaez
,
I.
Sanchez-Burgos
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
,
Phys. Rev. E
101
,
022611
(
2020
).
48.
J.
Diemand
,
R.
Angélil
,
K. K.
Tanaka
, and
H.
Tanaka
,
Phys. Rev. E
90
,
052407
(
2014
).
49.
Ø.
Wilhelmsen
,
D.
Bedeaux
, and
D.
Reguera
,
J. Chem. Phys.
142
,
171103
(
2015
).
50.
Ø.
Wilhelmsen
,
D.
Bedeaux
, and
D.
Reguera
,
J. Chem. Phys.
142
,
064706
(
2015
).
51.
A.
Aasen
,
E. M.
Blokhuis
, and
Ø.
Wilhelmsen
,
J. Chem. Phys.
148
,
204702
(
2018
).
52.
A.
Aasen
,
D.
Reguera
, and
Ø.
Wilhelmsen
,
Phys. Rev. Lett.
124
,
045701
(
2020
).
53.
E. M.
Blokhuis
and
A. E.
van Giessen
,
J. Phys.: Condens. Matter
25
,
225003
(
2013
).
54.
P.
Rehner
,
A.
Aasen
, and
Ø.
Wilhelmsen
,
J. Chem. Phys.
151
,
244710
(
2019
).
55.
S.
Stephan
,
J.
Liu
,
K.
Langenbach
,
W. G.
Chapman
, and
H.
Hasse
,
J. Phys. Chem. C
122
,
24705
(
2018
).
56.
A.
Obeidat
,
J.-S.
Li
, and
G.
Wilemski
,
J. Chem. Phys.
121
,
9510
(
2004
).
57.
J.
Rowlinson
and
B.
Widom
,
Molecular Theory of Capillarity
(
Clarendon Press
,
Oxford
,
1984
).
58.
T. L.
Hill
,
Thermodynamics of Small Systems
(
Courier Dover Publications
,
New York
,
2002
).
59.
M.
Heier
,
S.
Stephan
,
J.
Liu
,
W. G.
Chapman
,
H.
Hasse
, and
K.
Langenbach
,
Mol. Phys.
116
,
2083
(
2018
).
60.
T.
Van Westen
,
M.
Hammer
,
B.
Hafskjold
,
A.
Aasen
,
J.
Gross
, and
Ø.
Wilhelmsen
,
J. Chem. Phys.
156
,
104504
(
2022
).
61.
A. Z.
Panagiotopoulos
,
Mol. Phys.
61
,
813
(
1987
).
62.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
2017
).
63.
J. A.
Barker
and
D.
Henderson
,
J. Chem. Phys.
47
,
4714
(
1967
).
64.
R.
Roth
,
R.
Evans
,
A.
Lang
, and
G.
Kahl
,
J. Phys.: Condens. Matter
14
,
12063
(
2002
).
65.
R.
Stierle
,
E.
Sauer
,
J.
Eller
,
M.
Theiss
,
P.
Rehner
,
P.
Ackermann
, and
J.
Gross
,
Fluid Phase Equilib.
504
,
112306
(
2020
).
66.
M.
Hammer
,
G.
Bauer
,
R.
Stierle
,
J.
Gross
, and
Ø.
Wilhelmsen
,
J. Chem. Phys.
158
,
104107
(
2023
).
67.
J.
Mairhofer
and
J.
Gross
,
Fluid Phase Equilib.
444
,
1
(
2017
).
68.
J.
Vrabec
,
G. K.
Kedia
,
G.
Fuchs
, and
H.
Hasse
,
Mol. Phys.
104
,
1509
(
2006
).
69.
V.
Talanquer
,
J. Chem. Phys.
106
,
9957
(
1997
).
70.
V. K.
Shen
and
P. G.
Debenedetti
,
J. Chem. Phys.
114
,
4149
(
2001
).
71.
M.
Müller
,
L. G.
MacDowell
,
P.
Virnau
, and
K.
Binder
,
J. Chem. Phys.
117
,
5480
(
2002
).
72.
I.
Kusaka
,
J. Chem. Phys.
118
,
5510
(
2003
).
73.
I. E.
Parra
and
J. C.
Graňa
,
J. Chem. Phys.
132
,
034702
(
2010
).
74.
I.
Kusaka
,
J. Chem. Phys.
119
,
1808
(
2003
).
75.
K.
Koga
and
X. C.
Zeng
,
J. Chem. Phys.
110
,
3466
(
1999
).
76.
J.
Barrett
,
J. Chem. Phys.
111
,
5938
(
1999
).
77.
J. C.
Barrett
,
J. Chem. Phys.
131
,
084711
(
2009
).
78.
D.
Kashchiev
,
J. Chem. Phys.
118
,
1837
(
2003
).
79.
D.
Kashchiev
,
J. Chem. Phys.
153
,
124509
(
2020
).
80.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic Press
,
San Diego, CA
,
2002
).
81.
P. R.
ten Wolde
and
D.
Frenkel
,
J. Chem. Phys.
109
,
9901
(
1998
).
82.
F.
Zhu
and
G.
Hummer
,
J. Comput. Chem.
33
,
453
(
2012
).
83.
B.
Hafskjold
,
D.
Bedeaux
,
S.
Kjelstrup
, and
O.
Wilhelmsen
,
Phys. Rev. E
104
,
014131
(
2021
).
84.
A. E.
van Giessen
and
E. M.
Blokhuis
,
J. Chem. Phys.
131
,
164705
(
2009
).
85.
B.
Planková
,
V.
Vinš
, and
J.
Hrubý
,
J. Chem. Phys.
147
,
164702
(
2017
).
86.
R.
Checa
,
E.
Chacón
, and
P.
Tarazona
,
Phys. Rev. E
70
,
061601
(
2004
).
Published open access through an agreement with University of Barcelona

Supplementary Material