# Extending and validating bubble nucleation rate predictions in a Lennard-
Jones fluid with enhanced sampling methods and transition state theory

Kristof M. Bal<EMAIL_ADDRESS>Department of Chemistry and NANOlab
Center of Excellence, University of Antwerp, Universiteitsplein 1, 2610
Antwerp, Belgium  Erik C. Neyts  Department of Chemistry and NANOlab Center of
Excellence, University of Antwerp, Universiteitsplein 1, 2610 Antwerp, Belgium

(26 September 2022)

###### Abstract

We calculate bubble nucleation rates in a Lennard-Jones fluid through explicit
molecular dynamics simulations. Our approach—based on a recent free energy
method (dubbed reweighted Jarzynski sampling), transition state theory, and a
simple recrossing correction—allows us to probe a fairly wide range of rates
in several superheated and cavitation regimes in a consistent manner. Rate
predictions from this approach bridge disparate independent literature studies
on the same model system. As such, we find that rate predictions based on
classical nucleation theory, direct brute force molecular dynamics
simulations, and seeding are consistent with our approach and one another.
Published rates derived from forward flux sampling simulations are, however,
found to be outliers. This study serves two purposes. First, we validate the
reliability of common modeling techniques and extrapolation approaches on a
paradigmatic problem in materials science and chemical physics. Second, we
further test our highly generic recipe for rate calculations, and establish
its applicability to nucleation processes.

kinetics, free energy barriers, nucleation, transition state theory, bubbles

##  I Introduction

Nucleation processes are important throughout nature and technology, and are
therefore a long-standing research area within the physical and chemical
sciences. Even so, many uncertainties remain around our theoretical
understanding of nucleation mechanisms. [1] Atomistic simulations are, in
principle, ideally suited to study nucleation with a high level of detail.
Accurate, physically meaningful simulation setups are however challenging to
construct with many, sometimes subtle, possible sources of error. [2]

The difficulties associated with nucleation simulations are nicely illustrated
by a seemingly simple and rather well-defined system: Homogeneous bubble
nucleation in a Lennard-Jones (LJ) fluid. Different modeling studies disagree
about the nucleation mechanism, [3, 4] the nucleation rate, [3, 5, 6] the
validity of classical nucleation theory (CNT) for the process, [3, 5] and the
assumptions underpinning existing CNT models. [7, 8] These literature studies
employed different methodologies, such as direct brute force molecular
dynamics (MD), [4, 6] forward flux sampling (FFS), [3, 5] or seeding [6] which
have been compared to CNT predictions [4, 7, 8] Due to different time and
length scale restrictions of these methods, different physical conditions were
probed, resulting in limited overlap between rate data. A rigorous cross-
validation of these methodologies, and their respective rate predictions, is
therefore mostly lacking.

Recently, we proposed a generic strategy based on transition state theory
(TST) to evaluate rates of diverse processes. [9, 10] A key advantage of this
strategy is that it can fully leverage the rich toolbox of enhanced sampling
methods for free energy calculation. [11] As a result, very wide time scale
ranges can be studied within a consistent simulation strategy. [10]

In this work, we revisit the bubble nucleation in LJ fluids using this
methodology, simultaneously verifying existing rate predictions, and closing
the gap between disparate literature conditions.

##  II Methodology

Rate estimates for bubble nucleation are obtained in a system- and process-
agnostic manner based on recent developments in the field of enhanced sampling
approaches. [9, 10, 12] The approach employs molecular dynamics simulations in
relatively small simulation cells to ultimately yield macroscopic nucleation
rates over a wide range of conditions and time scales. The methodology
summarized here has already been successfully applied to droplet nucleation
from supersaturated vapor in a Lennard-Jones system. [10, 12]

###  II.1 Free energy calculation in a finite simulation cell

Nucleation is a rare event and therefore, in general, difficult to observe in
molecular simulations. One way to overcome the time scale problem, is the
application of a bias potential $V$. Such a bias, if properly designed, allows
to sample the sample critical bubbles as well as metastable liquid states
within the same simulation, so that a nucleation free energy surface (FES) can
be constructed from the (reweighted) marginal probability density along the
nucleation path.

In the specific approach used here, dubbed reweighted Jarzynski sampling, [12]
the bias potential is generated from a small number of nonequilibrium
simulations in which the system is pushed from the metastable liquid towards
states beyond the the critical bubble; a putative nucleation free energy
surface $\widetilde{G}$ is learned from the associated nonequilibrium work
distribution by an approximation of the Jarzynski equality. [13] The final
free energy estimate $G$ is then obtained from longer sampling runs under
influence of the bias potential $V=-\widetilde{G}$.

For any collective variable $\chi(\mathbf{R})$ that is a function of the
system coordinates $\mathbf{R}$, the free energy $G(\chi)$ at a constant
temperature $T$ and pressure $p$ is defined as

$G(\chi)=-k_{B}T\ln P(\chi),$ (1)

in which $k_{B}$ is the Boltzmann constant and $P(\chi)$ the marginal
probability distribution of $\chi$ under the considered conditions. $P(\chi)$
can be sampled in a biased simulation by accumulating the histogram for $\chi$
while reweighting each sampled data point $\chi_{i}$. The weight $w_{i}$ of
each sampled data point $i$ in the histogram is given by the umbrella sampling
relation $w_{i}=e^{V_{i}/k_{B}T}$, with $V_{i}$ the instantaneous value of a
time-independent $V$. [14]

Reweighted Jarzynski sampling addresses the exploration–convergence conundrum
inherent to adaptive bias enhanced sampling methods. [15] When using imperfect
collective variables, a single sampling run can either be optimized to target
many transitions (exploration) or a stable free energy estimate (convergence).
This issue is here essentially sidestepped by splitting up the simulation in a
forced exploration phase with a nonequilibrium bias and a convergence phase
with a non-adaptive bias.

An advantage of the outlined reweighting strategy is that $V$ need not be a
function of $\chi$. That is, biasing and sampling can be performed on
different collective variables. This is precisely what we will take advantage
of in this work. Most nucleation studies rely on reaction coordinates that
implement some variant of the ten Wolde–Frenkel parameter $n$, [16] which
counts the number of atoms that are part of the nucleating phase. For bubble
nucleation, this means that $n$ should count the number of vapor atoms. We use
a continuous expression for $n$ based on geometric switching functions to
count the number of atoms with less than five neighbors closer than
$1.6\sigma$. [3, 5] This manybody expression for $n$ is however very expensive
to evaluate, especially if also atomic gradients are needed to apply a bias
potential $V(n)$. Therefore, we perform initial nonequilibrium simulations by
biasing the average molar volume $v_{m}$ in the system, yielding
$\widetilde{G}(v_{m})=-V(v_{m})$. The final nucleation FES $G(n)$ is then
calculated from a subset of the configurations sampled in the biased
simulation under influence of $V(v_{m})$.

It must be noted that our implementation of $n$ is evaluated for all atoms in
the system, and thus counts all vapor-like atoms. This means that it is not
strictly the parameter defined by ten Wolde and Frenkel—according to CNT only
the largest cluster will drive the nucleation process, and $n$ should only
count the atoms in this nucleus. As previously discussed, we only require $n$
to be a sufficiently good reaction coordinate that allows us to distinguish
pre- and postcritical states and parametrize a dividing surface $n=n^{*}$ (as
described in the next section). [10] Such a set-up is consistent with previous
enhanced sampling studies of droplet nucleation, where droplet growth was
analyzed and driven by a parameter $n$ that counts all liquid-like atoms,
defined as atoms with more than five neighbors. [17, 18, 10]

###  II.2 Rates from transition state theory

Transition state theory offers a theoretical framework for the calculation of
rare event rates. If we wish to employ the Eyring formulation of the TST rate,

$k^{\mathrm{TST}}=\frac{k_{B}T}{h}e^{-\Delta^{\ddagger}G/k_{B}T},$ (2)

the free energy barrier $\Delta^{\ddagger}G$ must be defined as [19, 9]

$\Delta^{\ddagger}G=G(n^{*})+k_{B}T\ln\frac{\langle|\nabla
n|\rangle^{-1}_{n=n^{*}}}{h}\sqrt{2\pi mk_{B}T}-G_{l},$ (3)

with

$G_{l}=-k_{B}T\ln\int_{n<n^{*}}\mathrm{d}n\,e^{-G(n)/k_{B}T}.$ (4)

In these equations, $h$ is the Planck constant, $m$ the mass of the particles,
and $\langle|\nabla n|\rangle^{-1}_{n=n^{*}}$ the average norm of the gradient
of $n$ with respect to all atomic coordinates at $n=n^{*}$. The TST rate
measures the total flux through the dividing surface $n=n^{*}$, which is
always an upper bound to the effective rate of interest. Therefore, the value
of $n^{*}$, i.e., the location of the transition state, can be found by
maximizing $\Delta^{\ddagger}G$ and, hence, minimizing the rate.

The key assumption underpinning TST is that the candidate reaction coordinate
(here $n$) can correctly parametrize the dividing surface. This idea is
connected to the notion that $n$ is the slowest degree of freedom relevant to
the transition. A proper choice of the reaction coordinate is therefore
critical. It is possible to optimize the definition of the reaction coordinate
itself by maximizing its time scale separation with other degrees of freedom.
[20] Such an optimized reaction coordinate, which also contains information on
the nucleus shape, was developed in the context of droplet nucleation. [18]
The droplet analogue of $n$, when plugged into our procedure, however still
turned out to yield accurate rate estimations, which is why we also employ it
here. [10]

Two factors still separate the TST rate $k^{\mathrm{TST}}$ from the target
macroscopic nucleation rate $J$:

  1. 1.

The FES $G(n)$ and rate $k$ are only defined for the specific simulation
model, i.e., only within a small periodic simulation cell with a few thousand
atoms;

  2. 2.

The TST rate contains all crossings of the dividing surface, even those that
do not result into effective state-to-state transitions, which means that
$k\leq k^{\mathrm{TST}}$.

###  II.3 Obtaining the final macroscopic rate

The relation between the per-cell nucleation rate $k$ and macroscopic per-
volume rate $J$ is given by

$J=\frac{k}{V},$ (5)

provided that the cell volume $V$ is large enough to avoid self-interaction of
the nucleating bubble across periodic images.

Before arriving at $J$, we must first obtain $k$. The relation between $k$ and
$k^{\mathrm{TST}}$ can be expressed as

$k=\kappa k^{\mathrm{TST}},$ (6)

in which the transmission coefficient $\kappa$ accounts for recrossings of the
TST dividing surface.

We have recently proposed a simple strategy for the determination of $\kappa$
that fits in the workflow of our free energy calculation. [10] In order to
ascertain the ability of an approximate reaction coordinate to properly
discriminate between the dividing surface, and the (meta)stable states on that
it separates, one can use committor analysis of the putative dividing surface
$n=n^{*}$. As part of this committor analysis we prepare a number of
equilibrated configurations confined at $n=n^{*}$ by using restraints and
monitor their evolution after restraints are lifted. The committor $p_{g}$ is
the fraction of trajectories that results in a successful nucleation event,
i.e., further growth of the gas phase.

Committor analysis thus serves as an a posteriori validation of the employed
order parameter. A necessary condition for $n=n^{*}$ to be a dividing surface
for bubble nucleation is observing $p_{g}\approx 0.5$.

If $n=n^{*}$ passes the committor test (so $p_{g}=0.5$), additional insights
can be extracted from the collection of committor trajectories. The average
number of crossings of the dividing surface $\langle
j_{\mathrm{cross}}\rangle$ is correlated with $\kappa$. If we assume that
recrossings are the only contributor to $\kappa$, we have, by definition:

$\kappa=\frac{1}{2\langle j_{\mathrm{cross}}\rangle}.$ (7)

We have previously speculated that this definition of $\kappa$ might also
compensate deficiencies in the chosen reaction coordinate $n$, provided $n$
does not deviate too much from the true reaction coordinate. [10]

##  III Computational details

All simulations were carried out with LAMMPS [21, 22] and the PLUMED plugin.
[23, 24] Machine learning algorithms were used as implemented in the scikit-
learn library. [25]

Pairwise interatomic interactions were described using a truncated force-
shifted Lennard-Jones (TFS–LJ) potential

$U(r)=\phi(r)-\phi(r_{c})-(r-r_{c})\left|\frac{d\phi(r)}{dr}\right|_{r=r_{c}},$
(8)

in which $\phi(r)$ is the standard Lennard-Jones potential of the interatomic
distance $r$:

$\phi(r)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right].$
(9)

We set $\sigma=\epsilon=k_{B}=m=1$, thus using reduced units throughout. The
cutoff distance was $r_{c}=2.5$. These choices allowed us to compare computed
rates to literature simulations based on an identical potential. More
specifically, we revisited literature work on the $T=0.7$ and $T=0.855$
isotherms, [3, 5, 4] and the $p=0.026$ isobar of the same TFS–LJ fluid. [5, 4,
6] Sampling of the isobaric–isothermal (NpT) ensemble was performed with an
isotropic Nosé–Hoover style thermo- and barostat, [26] integrated with a 0.005
time step. Initial thermalization used a Langevin thermostat. [27]

Three distinct types of simulation were carried out for each condition.

In the first step, a set of 10 steered MD (SMD) simulations were carried out
to obtain a nonequilibrium work distribution. A moving harmonic restraint was
used to bring $v_{m}$ from a low value (close to that of the liquid in
equilibrium) to high value (beyond the critical nucleus) over $2\times 10^{6}$
MD steps. The range of $v_{m}$ was determined for each set of simulations on a
trial-and-error basis. The nonequilibrium work $W$ was recorded at regular
time steps. For each trajectory $i$, the noisy collection of $(v_{m}(t),W(t))$
data mapped to a smooth $W_{i}(v_{m})$ curve using kernel ridge regression
with a regularization strength of $10^{-3}$. Using the cumulant expansion of
the Jarzynski equality, [28] the full set of $W_{i}$ curves was then finally
converted into $\widetilde{G}(v_{m})$, which was then fitted to an artificial
neural network (ANN) consisting of a single 12-neuron hidden layer.

The second step involved five longer equilibrium sampling runs using the ANN
bias $V(v_{m})=-\widetilde{G}(v_{m})$. Sampling efficiencies were improved by
limiting the simulation to a narrow range of $v_{m}$, based on the shape of
$\widetilde{G}(v_{m})$; harmonic restraints were used to prevent the system
from moving too far beyond the critical nucleus. The final free energy surface
$G(n)$ was then obtained from the reweighted histogram $P(n)$ of $n$ according
to Eq. (1). Sufficient sampling was found to be possible within $2\times
10^{7}$ MD steps per simulation. The histogram was accumulated using the
kernel density estimation functionality within PLUMED, and written to a grid.

The third and final step of the procedure was committor analysis—after
evaluation of the barrier Eq. (3), TST rate Eq. (2), and identification of the
approximate dividing surface $n=n^{*}$. SMD simulations were used to prepare
10 configurations at $n=n^{*}$ over $2\times 10^{6}$ MD steps each. The
subsequent committor trajectories were analyzed to verify that $p_{g}\approx
0.5$, and $\langle j_{\mathrm{cross}}\rangle$ was then also evaluated. In this
step the $n$ CV was biased. A multiple time stepping scheme [29] was used with
a stride of 10 to limit the number of expensive force evaluations on $n$.

Representative input decks that contain the full set of simulation parameters
and implementation details are available on PLUMED-NEST (www.plumed-nest.org),
the public repository of the PLUMED consortium [24], as plumID:22.025.[30]

##  IV Results and discussion

###  IV.1 Free energy surface and choice of CV

Figure 1:  Bubble nucleation FES at (a–b) $p=0.01438$ and (c–d) $p=0.026$ on
the $T=0.855$ isotherm of the TFS–LJ fluid, projected on $v_{m}$ and $n$,
respectively. Reduced LJ units are used throughout, meaning that $G$ is in
units of $\epsilon$.

We have used two CVs in our simulation: The inexpensive molar volume CV
$v_{m}$ for biasing, and the expensive manybody CV $n$ for evaluation of the
final barrier. The effectiveness of this choice hinges on the assumption that
a FES projected on $v_{m}$ sufficiently captures the barrier on $G(n)$. That
is, we assume that although $n$ is the best possible reaction coordinate and
the true barrier can be evaluated from $G(n)$, the barrier on $G(v_{m})$ only
misses the true barrier by a few $k_{B}T$. $V(v_{m})=-G(v_{m})$, while
imperfect, therefore ought to be a perfectly cromulent bias to allow for
frequent barrier crossing and good sampling.

We illustrate the applicability of our approach for two conditions on the
$T=0.855$ isotherm. They represent a low barrier ($p=0.01438$, Figure 1a–b)
and a higher barrier case ($p=0.026$, Figure 1c–d). In both systems, the local
maximum of $G(v_{m})$ is lower than that of $G(n)$. The relative valley-to-
peak height of the FES maximum is not, strictly speaking, the free energy
barrier, although it is usually quite close. [9] In any case, it might be
inferred that $n$ is a somewhat better reaction coordinate because it
discriminates more sharply between metastable states and is the order
parameter with the largest spectral gap. [20] In absolute terms, however, the
difference remains in the order of 1–2 $k_{B}T$. A bias based on $G(v_{m})$
can therefore still be successful.

Figure 2:  FES $G(n)$ and values of the inverted bias
$-V(v_{m})=\widetilde{G}(v_{m})$ recorded at different times $t$, plotted at
$n(t)$ for (a) $p=0.01438$ and (b) $p=0.026$.

The bias $V(v_{m})$ is based on the approximate $\widetilde{G}(v_{m})$, rather
than the true $G(v_{m})$. Both $n$ and $v_{m}$ are function of $\mathbf{R}$,
which is in turn a function of time. To illustrate how well the bias
$V(v_{m})$ compensates for the underlying FES $G(n)$, we can plot
$(n(t),-V(v_{m}(t)))$ points over $G(n)$. An example of such analysis is given
in Figure 2. The effective bias does not perfectly match $G(n)$, which is
especially visible in the transition state region. $v_{m}$ cannot sharply
discriminate configurations with critical bubbles quite as well as $n$. This
is because $v_{m}$ is a _global_ order parameter, whereas bubble formation
entails a _local_ density change. As a result, configurations not directly
around the top of $G(n)$ receive a bias meant for transition state
configurations, and are therefore subjected to only a small biasing force.
This diminishes the ability of the bias to facilitate transitions. However,
the limited difference in barrier height on $G(n)$ and $G(v_{m})$,
respectively, means that the bias $V(v_{m})$ does adequately compensate the
FES on average.

###  IV.2 Superheated isotherm

Figure 3:  Rate predictions from our approach and literature simulations, for the $T=0.855$ superheated isotherm of the TFS–LJ fluid. (a) Extrapolation of our data and MD rates using CNT relations. (b) Comparing explicit rate predictions from simulation to CNT predictions based on different literature corrections to the surface tension (see text for details). A first order correction employing a Tolman length $\delta_{T}=0.128$ yields rates that overlap with the second order line. Table 1:  TST nucleation barriers $\Delta^{\ddagger}G$, extensive per-cell TST rates $k^{\mathrm{TST}}$, transmission coefficients $\kappa$ and final rate estimates $J$ for bubble nucleation on the $T=0.855$ isotherm of the TFS–LJ fluid at different pressures $p$.222Reduced LJ units are used throughout. For comparison, predicted critical bubble radii $R_{s}^{*}$ and CNT barriers $G^{*}$ are reported for the “best fit” CNT estimate, i.e., based on a second order correction to the surface tension. $p$ | $\Delta^{\ddagger}G$ | $k^{\mathrm{TST}}$ | $\kappa$ | $J$ | $R_{s}^{*}$ | $G^{*}$
---|---|---|---|---|---|---
0.01438 | $2.11\pm 0.13$ | $1.15\times 10^{-2\pm 0.06}$ | $0.05\pm 0.02$ | $1.99\times 10^{-8\pm 0.2}$ | $6.23\pm 0.17$ | $13.96\pm 0.65$
0.01701 | $4.64\pm 0.19$ | $5.97\times 10^{-4\pm 0.1}$ | $0.04\pm 0.01$ | $7.11\times 10^{-10\pm 0.2}$ | $6.82\pm 0.18$ | $16.77\pm 0.78$
0.02004 | $8.58\pm 0.23$ | $5.99\times 10^{-6\pm 0.1}$ | $0.03\pm 0.01$ | $5.45\times 10^{-12\pm 0.2}$ | $7.64\pm 0.20$ | $21.15\pm 0.98$
0.02383 | $16.03\pm 0.34$ | $9.80\times 10^{-10\pm 0.2}$ | $0.02\pm 0.01$ | $6.21\times 10^{-16\pm 0.2}$ | $8.98\pm 0.24$ | $29.40\pm 1.37$
0.02600 | $23.31\pm 0.36$ | $1.96\times 10^{-13\pm 0.2}$ | $0.02\pm 0.01$ | $1.14\times 10^{-19\pm 0.2}$ | $9.98\pm 0.27$ | $36.40\pm 1.69$

We have probed five conditions in the superheated regime between $p=0.01438$
and $p=0.026$ on the $T=0.855$ isotherm, well below the coexistence pressure
$p_{\mathrm{eq}}=0.0461$. Nonequilibrium SMD simulation were carried out from
$v_{m}=1.7$ to $v_{m}=2.2$ to parametrize $V(v_{m})$ in cubic simulation cells
containing $N=17576$ atoms. The computed barriers and rates are summarized in
Table 1.

Our rate estimates can be compared with previously published data for the same
system. On one end, Diemand et al. [4] carried out MD simulations at low
pressure ($p\leq 0.02383$) in very large simulation cells ($N\approx 5\times
10^{8}$) in which several nucleation events could be observed directly. On the
other end, Wang et al. [3] and Meadley & Escobedo [5] used forward flux
sampling (FFS) to obtain rates at $p=0.026$ in small systems of $N=3375$ and
$N=8000$ atoms, respectively. We have simulated the same pressures $p$ as
these literature studies to allow for a direct one-to-one comparison. The
advantage of our approach is that it allows us the evaluate the nucleation
rate over the full considered pressure range, so that disparate methodologies
can now be compared to one consistent data set.

It can be seen from Figure 3a that our rate estimates closely match large-
scale MD data at low $p$, but are orders of magnitude lower than FFS
predictions at $p=0.026$. The FFS estimates of Wang et al. overshoot our
values the most: As pointed out by Meadley & Escobedo, these rates are likely
affected by finite size effects due to the small cell size that was employed
in the FFS simulations. Nevertheless, even Meadley & Escobedo’s simulations in
a larger cell produce a rate that is over two orders of magnitude higher than
our estimate.

Rates from low-pressure MD simulations can, in principle, be extrapolated to
higher pressures through CNT-derived relations. According to CNT, $\ln
J\sim(\Delta p)^{-2}$, in which $\Delta p$ is the pressure difference between
the bubble and the surrounding liquid. Assuming that $\Delta p$ is
proportional to $p_{\mathrm{eq}}-p$, it thus becomes possible to extrapolate
the MD results, as shown in Figure 3a. The extrapolated MD results match our
TST-based estimates within an order of magnitude over the whole considered
range. In addition, the TST-derived rates are very well-represented by a
linear fit of $\log J$ to $(p_{\mathrm{eq}}-p)^{-2}$.

The applicability of CNT for bubble nucleation in superheated liquids has been
controversial. [3, 4, 7, 8] Wang et al. analyzed temperature profiles in their
FFS simulations. They concluded that bubble formation is driven by the
occurrence of local hot spots in the liquid. [3] Because CNT (and also TST)
are based on the assumption of local thermal equilibrium, such a fact would
undermine the validity of the theory. Diemand et al., however, observed that
local hot spots never _precede_ bubble formation. Rather, hot spots are a
consequence of the excess kinetic energy that is produced by rapid expansion
of the bubble. [4] Hot spots therefore only occur _after_ a larger-than-
critical bubble is already formed. The good cross-agreement between MD data,
TST, and CNT corroborates this observation.

CNT in principle offers a path to nucleation rates that does not require any
explicit simulation of the nucleation process. Only a few bulk properties at
coexistence are needed: the vapor/liquid surface tension $\gamma$, pressure
$p_{\mathrm{eq}}$, liquid density $\rho_{l}$ and vapor density $\rho_{g}$.
Then, the CNT nucleation barrier $G^{*}$ is

$G^{*}=\frac{16\pi}{3}\frac{\gamma^{3}}{(\Delta p)^{2}},$ (10)

in which $\Delta p$ is the pressure difference between the pressure inside and
outside the bubble, which can be approximated though the Poynting correction
as $\Delta p=(p_{\mathrm{eq}}-p)\delta_{P}$ in which

$\delta_{P}\approx
1-\left(\frac{\rho_{g}}{\rho_{l}}\right)+\frac{1}{2}\left(\frac{\rho_{g}}{\rho_{l}}\right)^{2}.$
(11)

At $T=0.855$, $\rho_{l}=0.729$ and $\rho_{g}=0.0198$, so $\delta_{P}=0.870$.
The nucleation rate can then be expressed as

$J=\rho_{l}\left(\frac{2\gamma}{\pi m}\right)^{1/2}e^{-G^{*}/k_{B}T}.$ (12)

Both our TST-based approach and direct application of CNT rely on the concept
of a free energy barrier on which the rates depends exponentially,
supplemented by some prefactor. In addition, both approaches assume that the
system remains in local thermal equilibrium throughout the nucleation process,
i.e., that there exists a time scale separation between bubble growth (slow)
and atomic motion along other degrees of freedom (fast).

There are also two key differences. First of all, the relation Eq. (12) can be
derived from the explicit mechanistic assumption of a single bubble that grows
or shrinks through evaporation or condensation of single atoms. The TST
equation Eq. (2), in contrast, merely describes the total flux through the
dividing surface parametrized by the CV of choice; recrossings are explicitly
accounted for by the committor analysis. Second, the CNT barrier $G^{*}$ and
TST barrier $\Delta^{\ddagger}G$ have different meanings. The former is an
intensive property, while the latter is extensive, as discussed before. [31,
10] The extensive nature of $\Delta^{\ddagger}G$ therefore explains its lower
absolute value in our large simulation cells, as can be seen in Table 1.
$\Delta^{\ddagger}G$ is also reconstructed directly from an explicitly sampled
thermodynamic ensemble of the simulated system, rather than computed from
equilibrium properties only. A direct comparison between the two approaches is
therefore only possible on the basis of rates $J$.

Note that uncertainties in $\gamma$ (of about 20%) have a rather large impact
on $G^{*}$ and, therefore, result in CNT estimates of $J$ with an uncertainty
of around one order of magnitude—larger than the uncertainties on our TST-
based estimates (Table 1). [4]

A direct, zeroth order, application of CNT based on macroscopic input data
yields fair agreement with our data and MD results, in line with observations
of Diemand et al. (Figure 3b). However, discrepancies remain. The CNT rates
can be further improved by taking into account the curvature dependence of
surface tension $\gamma$. A first order correction is given by the Tolman
equation [32]

$\gamma(R_{s})=\frac{\gamma_{\infty}}{1+2\delta_{T}/R_{s}}.$ (13)

Here, $R_{s}$ is the radius of the critical bubble or droplet,
$\gamma_{\infty}$ the surface tension of a planar interface
($R_{s}\rightarrow\infty$), and $\delta_{T}$ the Tolman length. $\delta_{T}$
is defined as the difference between the Gibbs equimolar radius $R_{e}$ and
radius of the surface of tension $R_{s}$, although it is in practice an
empirical parameter in the context of rate calculation. Diemand et al. found a
value of $\delta_{T}=0.25$ to fit their MD data well. [4] Schmelzer & Baidakov
later questioned the validity of a first order correction, especially for
small $R_{s}$. [8] They argued that only a second order correction is able to
accurately describe a wider range of curvatures:

$\gamma(R_{s})=\frac{\gamma_{\infty}}{1+2\delta_{T}/R_{s}+l^{2}/R_{s}^{2}},$
(14)

where they proposed $\delta_{T}=0.128$ and $l^{2}=1.56$. To evaluate these
corrections in the context of the rate we must know $R_{s}^{*}$, the radius of
the critical nucleus, which CNT gives as $R_{s}^{*}=2\gamma(R_{s}^{*})/\Delta
p$. We therefore obtain $\gamma(R_{s}^{*})$ iteratively.

Application of the different corrections to $\gamma$ reveals their impact on
the rate estimate, as can be seen in Figure 3b. Setting
$\gamma=\gamma_{\infty}=0.0895$ results in an underestimation of the MD data
and our results. First order corrections bring the CNT estimate closer in
absolute terms, although they appear to overcorrect somewhat. Second order
corrections then shift the rate back down. This final data set has the closest
agreement to the other rates estimates (except FFS) over the full pressure
range. Yet, the simple first order Tolman correction still appears to be
adequate in this particular regime: If we set $\delta_{T}=0.128$ (i.e., using
the first order coefficient of the second order expansion) rate predictions
are almost indistinguishable from the second order correction. It appears
that, at least under the present conditions, the good performance of the
proposed second order correction can be attributed to its first order
component. For similar bubble sizes, Sanchez-Burgos et al. found a first order
$\delta_{T}\approx 0.15$ from the extrapolation
$\delta_{T}=\displaystyle\lim_{R_{s}\rightarrow\infty}(R_{e}-R_{s})$. [33]

###  IV.3 Cavitation isotherm

Figure 4:  Rate predictions from our approach and literature simulations, for the $T=0.700$ cavitation isotherm of the TFS–LJ fluid. Comparing explicit rate predictions from simulation to CNT predictions based on different literature corrections to the surface tension (see text for details). A first order correction employing a Tolman length $\delta_{T}=0.128$ yields rates that overlap with the second order line. Table 2:  TST nucleation barriers $\Delta^{\ddagger}G$, extensive per-cell TST rates $k^{\mathrm{TST}}$, transmission coefficients $\kappa$ and final rate estimates $J$ for bubble nucleation on the $T=0.700$ isotherm of the TFS–LJ fluid at different pressures $p$. Predicted critical bubble radii $R_{s}^{*}$ and CNT barriers $G^{*}$ are calculated as in Table 1. $p$ | $\Delta^{\ddagger}G$ | $k^{\mathrm{TST}}$ | $\kappa$ | $J$ | $R_{s}^{*}$ | $G^{*}$
---|---|---|---|---|---|---
$-$0.16627 | $6.32\pm 0.13$ | $1.35\times 10^{-5\pm 0.1}$ | $0.05\pm 0.02$ | $6.66\times 10^{-11\pm 0.2}$ | $3.53\pm 0.04$ | $16.02\pm 0.34$
$-$0.16222 | $7.01\pm 0.23$ | $4.97\times 10^{-6\pm 0.1}$ | $0.06\pm 0.03$ | $2.60\times 10^{-11\pm 0.3}$ | $3.62\pm 0.04$ | $16.86\pm 0.36$
$-$0.15808 | $7.78\pm 0.24$ | $1.67\times 10^{-6\pm 0.2}$ | $0.06\pm 0.02$ | $9.74\times 10^{-12\pm 0.2}$ | $3.72\pm 0.05$ | $17.79\pm 0.37$
$-$0.15000 | $9.66\pm 0.23$ | $1.13\times 10^{-7\pm 0.1}$ | $0.08\pm 0.03$ | $8.05\times 10^{-13\pm 0.2}$ | $3.92\pm 0.05$ | $19.81\pm 0.42$
$-$0.14724 | $10.94\pm 0.15$ | $1.82\times 10^{-8\pm 0.1}$ | $0.09\pm 0.04$ | $1.46\times 10^{-13\pm 0.2}$ | $3.99\pm 0.05$ | $20.57\pm 0.43$
$-$0.14000 | $13.07\pm 0.14$ | $8.66\times 10^{-10\pm 0.1}$ | $0.14\pm 0.05$ | $1.07\times 10^{-14\pm 0.2}$ | $4.19\pm 0.05$ | $22.79\pm 0.48$

We considered six conditions on the $T=0.700$ isotherm, where negative
pressures are required to induce nucleation via cavitation and the coexistence
pressure is $p_{\mathrm{eq}}=0.01186$. Nonequilibrium SMD simulation were
carried out at pressures between $p=-0.16627$ and $-0.14000$ from $v_{m}=1.4$
to $v_{m}=1.6$. Due to the smaller bubble sizes, cubic simulation cells with
$N=4096$ atoms were found to be sufficient. The computed barriers and rates
are summarized in Table 2.

MD data are available for the lowest pressures ($p\leq-0.15808)$, [4] and
match our predictions well. At higher pressures, only one FFS estimate is
available, at $p=-0.15$. [5] The FFS rate overestimates our prediction by an
order of magnitude, similar to the superheated case. Our rate predictions
appear to follow a CNT-style dependence on the pressure (Figure 4).

We verify the direct applicability of CNT-based formulas for rate calculation
in this system. Relevant parameters for this computation are
$\rho_{g}=0.0198$, $\rho_{l}=0.729$ and $\gamma_{\infty}=0.329$. Assuming a
curvature-independent $\gamma$ leads to a very large underestimation of the
rate, an issue exacerbated in this system by the smaller bubble size (and
larger curvatures) compared to the superheated system. A first order
correction using $\delta_{T}=0.25$ overcorrects, however. It therefore appears
that the second order correction is necessary, since it does lead to an
excellent agreement with our data and MD results. Yet, as in the superheated
example, a first order correction with a smaller $\delta_{T}=0.128$ performs
almost identically.

###  IV.4 Superheated isobar

Figure 5:  Rate predictions from our approach and literature simulations, for the $p=0.026$ superheated isobar of the TFS–LJ fluid. Comparing explicit rate predictions from simulation to CNT predictions based on literature second order corrections to the surface tension. Table 3:  TST nucleation barriers $\Delta^{\ddagger}G$, extensive per-cell TST rates $k^{\mathrm{TST}}$, transmission coefficients $\kappa$ and final rate estimates $J$ for bubble nucleation on the $p=0.026$ isobar of the TFS–LJ fluid at different temperatures $T$. Predicted critical bubble radii $R_{s}^{*}$ and CNT barriers $G^{*}$ are calculated as in Table 1, using empirical linear functions $\gamma_{\infty}(T)$ and $\Delta p(T)$.444The different calculation of $\gamma_{\infty}$ and $\Delta p$ explains the different numerical values and error bars for $T=0.855$ compared to Table 1. $T$ | $\Delta^{\ddagger}G$ | $k^{\mathrm{TST}}$ | $\kappa$ | $J$ | $R_{s}^{*}$ | $G^{*}$
---|---|---|---|---|---|---
0.8550 | $23.31\pm 0.36$ | $1.96\times 10^{-13\pm 0.2}$ | $0.02\pm 0.01$ | $1.14\times 10^{-19\pm 0.2}$ | $9.98\pm 0.61$ | $35.89\pm 3.79$
0.8580 | $15.21\pm 0.24$ | $2.74\times 10^{-9\pm 0.1}$ | $0.05\pm 0.02$ | $4.64\times 10^{-15\pm 0.2}$ | $9.03\pm 0.58$ | $27.72\pm 3.10$
0.8600 | $11.79\pm 0.23$ | $1.53\times 10^{-7\pm 0.1}$ | $0.02\pm 0.01$ | $9.00\times 10^{-14\pm 0.2}$ | $8.43\pm 0.57$ | $23.24\pm 2.71$
0.8640 | $5.65\pm 0.16$ | $1.99\times 10^{-4\pm 0.1}$ | $0.03\pm 0.01$ | $2.03\times 10^{-10\pm 0.2}$ | $7.34\pm 0.54$ | $16.14\pm 2.05$
0.8660 | $3.40\pm 0.10$ | $2.73\times 10^{-3\pm 0.1}$ | $0.05\pm 0.02$ | $4.86\times 10^{-9\pm 0.2}$ | $6.83\pm 0.52$ | $13.35\pm 1.77$
0.8675 | $1.80\pm 0.14$ | $1.73\times 10^{-2\pm 0.1}$ | $0.04\pm 0.02$ | $2.51\times 10^{-8\pm 0.2}$ | $6.46\pm 0.51$ | $11.54\pm 1.59$
0.8680 | $1.44\pm 0.20$ | $2.64\times 10^{-2\pm 0.1}$ | $0.06\pm 0.02$ | $5.39\times 10^{-8\pm 0.2}$ | $6.35\pm 0.51$ | $10.98\pm 1.53$

Rosales-Pelaez et al. studied nucleation on the $p=0.026$ isobar using direct
brute force MD simulations as well as a seeding approach. [6] We consider
temperatures within the same range, in casu between $T=0.855$ and $T=0.868$.
These conditions intersect with our work on the superheated isotherm for
$(p,T)=(0.026,0.855)$, for which FFS results are also available. [3, 5]
Likewise, initial nonequilibrium sampling was carried out from $v_{m}=1.7$ to
$v_{m}=2.2$ in cubic simulation cells containing $N=17576$ atoms. The computed
barriers and rates are summarized in Table 3.

Rosales-Pelaez et al. produced two sets of rate estimates from seeding, which
differed in how the critical bubble radius was defined. They noted that
neither set could simultaneously match their MD results at high temperatures,
and literature FFS data at $T=0.855$. As it turns out, seeding predictions
based on the Gibbs dividing surface (GDS)—which Rosales-Pelaez et al. found to
be consistent with MD simulations at low temperature—are further validated by
our data (Figure 5). While the FFS estimates at $T=0.855$ overshoot our rate
prediction by orders of magnitude, as noted earlier, the GDS-based rate
derived from seeding at $T=0.858$ agrees very closely with our prediction.

From a CNT point of view, the pressure-dependence of rate is primarily
impacted by the change in $\Delta p$. The effect of a change in temperature is
somewhat more complex to account for, since $\gamma$ and $\rho_{g}$ and
$\rho_{l}$ are also temperature-dependent. This means that extrapolation of
nucleation rates to lower $T$ is more difficult. We can use the empirical
relations $\gamma_{\infty}(T)=1.3557-1.4809T$ and $\Delta
p(T)=-0.2071+0.2624T$. [6] Compared to literature, we have shifted the
$\gamma_{\infty}(T)$ curve to match the value of $\gamma_{\infty}$ at
$T=0.855$ reported by Diemand et al. [4] and allow for consistency with our
simulations on the $T=0.855$ isotherm. If we further assume $\rho_{g}$ and
$\rho_{l}$ to be constant in the considered temperature range, and use the
second-order correction to $\gamma$, we can plot CNT-predicted rates in Figure
5. Although a linear fit of our $\log J$ data to $(\Delta p)^{2}$ might seem
appropriate within the studied range, the nonlinear behavior of the CNT rate
reveals that linear extrapolation far beyond the simulated temperature range
is not advisable. The uncertainties on the $\gamma_{\infty}(T)$ fit are also
fairly high, which affects the fidelity of CNT for rate predictions (compare
$R_{s}^{*}$ and $G^{*}$ in Tables 1 and 3). If accurate temperature-dependent
surface tensions and coexistence pressures are known, however, a
straightforward application of CNT accurately describes the nucleation rate.

###  IV.5 General remarks

Our results allow us to reconcile rate predictions from direct MD simulations,
seeding approaches, and CNT. By combining a modern, efficient free energy
method with a rate evaluation based on a generic implementation of TST, we are
able to span the full range of rates and conditions that were previously
reported within a single, consistent approach that is mostly free of
mechanistic assumptions. It therefore becomes easier to cross-validate other
rate computation approaches, and to investigate the applicability of
extrapolation techniques. We note that our approach was already validated for
droplet nucleation processes, highlighting its general applicability. [10]

One consistent outlier has been forward flux sampling (FFS), yielding bubble
nucleation rates that are consistently above the range established by other
methods by several orders of magnitude. Finite size effects can explain some
of the inconsistencies, [5] but not all; the FFS rate overestimation persists
even for the cavitation process when using larger cells than in our
simulations. [5] Large discrepancies, up to seven orders of magnitude, between
FFS and other methods have also been reported for ice nucleation. [34] Given
that FFS has become a workhorse tool for the computation of nucleation rates
in diverse systems, [35] it might be prudent to further validate the method
and investigate the origins of such large deviations. That being said, certain
deficiencies of FFS in the context of nucleation have already been addressed
in later developments of the method. [36]

Cross-method benchmarking, reproducibility, and accuracy of nucleation rate
prediction has seen a renewed interest, as recently discussed in an excellent
review. [2] We hope our benchmarking strategies and methodological insights
can also contribute to this discussion.

We note that the methodology as outlined in in Sec. II is not completely set
in stone. Its only fixed aspects are the reconstruction of the nucleation free
energy surface $G(\chi)$ along a suitable order parameter $\chi$—which can be
fed into the TST expressions Eqs. (2)–(4)—followed by committor analysis and
the recrossing correction of Eq. (7). There exists a wide variety of methods
that can reconstruct free energies, each with their own characteristics and
strengths. [11] The particular choice of free energy method in practice mostly
depends on the problem at hand, the user’s expertise, and code availability.
We eat our own cooking by using reweighted Jarzynski sampling, but are also
motivated by its good sampling efficiency and convergence in head-to-head
comparisons with established adaptive bias methods. [12] The specific choice
of reaction coordinate $\chi$ will be system-dependent and can take much more
sophisticated forms in studies of crystallization. [37, 38, 39] As noted
before, sampling efficiency and accuracy of rate is highly dependent on an
appropriate choice of this reaction coordinate. That being said, the choice of
$\chi$ is the _only_ system-specific aspect of the overall procedure. Our rate
calculation procedure is therefore as generic and process-agnostic as the free
energy methods on which it is based, equally applicable to nucleation as it is
to chemical reactions. [10]

Finally, we have confirmed the validity of CNT for the considered conditions.
First, relations derived from CNT allow for the extrapolation of explicit rate
predictions to conditions that are more difficult to sample. Second, if high-
quality equilibrium properties of the liquid are known, it is possible to
produce accurate rate estimates with CNT. An appropriate curvature correction
to the surface tension is however required. Although we find that a first
order Tolman correction is sufficient in our case, it is possible that second
order corrections are needed for a consistent description of even larger
curvatures, i.e., in the case of higher superheatings (or supercoolings in
case of droplet nucleation). [8] A recent study demonstrated that a consistent
simultaneous treatment of bubble and droplet nucleation is possible using the
same first order Tolman correction, although only systems with fairly large
bubble/droplet sizes were fitted. [33]

##  V Conclusions

We have revisited a paradigmatic nucleation problem—rate calculation of bubble
formation in a Lennard-Jones fluid—with a set of recently developed enhanced
sampling methodologies based on the Jarzynski equality, transition state
theory, and a simple recrossing correction. The employed methodology is highly
generic—the only system- or process-dependent aspect being the choice of an
approximate reaction coordinate. A very wide range of simulation conditions
can thus be consistently probed within an identical simulation setup. We were
therefore able to reconcile disparate literature data on the same system.

On one hand, we validate simulated nucleation rate predictions based on brute
force molecular dynamics and seeding approaches, while forward flux sampling
(FFS) appears to produce outliers. We also show that simple analytical
expressions derived from classical nucleation theory can produce satisfactory
rate predictions, provided that the curvature dependence of the surface
tension is empirically corrected. Another viable approach is to extend
explicit rate predictions to a wider range of conditions with extrapolation
schemes that exploit CNT trends.

On the other hand, our results confirm that our generic TST-based approach to
rate calculation is equally applicable to complex nucleation processes as it
is to simple chemical reactions. It may therefore be a useful additional tool
for rate calculation in future modeling studies, either to be used on its own
or for cross-validation purposes. We therefore hope that our insights can help
improve the accuracy and reproducibility of rate calculation also for other
types of (nucleation) processes.

###### Acknowledgements.

K.M.B. was funded as a junior postdoctoral fellow of the FWO (Research
Foundation – Flanders), Grant 12ZI420N. The computational resources and
services used in this work were provided by the HPC core facility CalcUA of
the Universiteit Antwerpen, and VSC (Flemish Supercomputer Center), funded by
the FWO and the Flemish Government.

## Author Declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request. Sample inputs to reproduce the
reported simulations are deposited on PLUMED-NEST (www.plumed-nest.org), the
public repository of the PLUMED consortium [24], as plumID:22.025.[30]

## References

  * Karthika, Radhakrishnan, and Kalaichelvi [2016] S. Karthika, T. K. Radhakrishnan, and P. Kalaichelvi, “A review of classical and nonclassical nucleation theories,” Cryst. Growth Des. 16, 6663–6681 (2016).
  * Blow, Quigley, and Sosso [2021] K. E. Blow, D. Quigley, and G. C. Sosso, “The seven deadly sins: When computing crystal nucleation rates, the devil is in the details,” J. Chem. Phys. 155, 040901 (2021).
  * Wang, Valeriani, and Frenkel [2009] Z.-J. Wang, C. Valeriani, and D. Frenkel, “Homogeneous bubble nucleation driven by local hot spots: A molecular dynamics study,” J. Phys. Chem. B 113, 3776–3784 (2009).
  * Diemand  _et al._ [2014] J. Diemand, R. Angélil, K. K. Tanaka, and H. Tanaka, “Direct simulations of homogeneous bubble nucleation: Agreement with classical nucleation theory and no local hot spots,” Phys. Rev. E 90, 052407 (2014).
  * Meadley and Escobedo [2012] S. L. Meadley and F. A. Escobedo, “Thermodynamics and kinetics of bubble nucleation: Simulation methodology,” J. Chem. Phys. 137, 074109 (2012).
  * Rosales-Pelaez  _et al._ [2019] P. Rosales-Pelaez, M. I. Garcia-Cid, C. Valeriani, C. Vega, and E. Sanz, “Seeding approach to bubble nucleation in superheated Lennard-Jones fluids,” Phys. Rev. E 100, 052609 (2019).
  * Tanaka  _et al._ [2015] K. K. Tanaka, H. Tanaka, R. Angélil, and J. Diemand, “Simple improvements to classical bubble nucleation models,” Phys. Rev. E 92, 022401 (2015).
  * Schmelzer and Baidakov [2016] J. W. P. Schmelzer and V. G. Baidakov, “Comment on “Simple improvements to classical bubble nucleation models”,” Phys. Rev. E 94, 026801 (2016).
  * Bal  _et al._ [2020] K. M. Bal, S. Fukuhara, Y. Shibuta, and E. C. Neyts, “Free energy barriers from biased molecular dynamics simulations,” J. Chem. Phys. 153, 114118 (2020).
  * Bal [2021a] K. M. Bal, “Nucleation rates from small scale atomistic simulations and transition state theory,” J. Chem. Phys. 155, 144111 (2021a).
  * Hénin  _et al._ [2022] J. Hénin, T. Lelièvre, M. R. Shirts, O. Valsson, and L. Delemotte, “Enhanced sampling methods for molecular dynamics simulations,”  (2022), arXiv:2202.04164 .
  * Bal [2021b] K. M. Bal, “Reweighted Jarzynski sampling: Acceleration of rare events and free energy calculation with a bias potential learned from nonequilibrium work,” J. Chem. Theory Comput. 17, 6766–6774 (2021b).
  * Jarzynski [1997] C. Jarzynski, “Nonequilibrium equality for free energy differences,” Phys. Rev. Lett. 78, 2690–2693 (1997).
  * Torrie and Valleau [1977] G. M. Torrie and J. P. Valleau, “Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling,” J. Comput. Phys. 23, 187–199 (1977).
  * Invernizzi and Parrinello [2022] M. Invernizzi and M. Parrinello, “Exploration vs convergence speed in adaptive-bias enhanced sampling,” J. Chem. Theory Comput. 18, 3988–3996 (2022).
  * ten Wolde and Frenkel [1998] P. R. ten Wolde and D. Frenkel, “Computer simulation study of gas–liquid nucleation in a Lennard-Jones system,” J. Chem. Phys. 109, 9901–9918 (1998).
  * Salvalaglio  _et al._ [2016] 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).
  * Tsai, Smith, and Tiwary [2019] S.-T. Tsai, Z. Smith, and P. Tiwary, “Reaction coordinates and rate constants for liquid droplet nucleation: Quantifying the interplay between driving force and memory,” J. Chem. Phys. 151, 154106 (2019).
  * Vanden-Eijnden and Tal [2005] E. Vanden-Eijnden and F. A. Tal, “Transition state theory: Variational formulation, dynamical corrections, and error estimates,” J. Chem. Phys. 123, 184103 (2005).
  * Tiwary and Berne [2016] P. Tiwary and B. J. Berne, “Spectral gap optimization of order parameters for sampling complex molecular systems,” Proc. Natl. Acad. Sci. 113, 2839–2844 (2016).
  * Plimpton [1995] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117, 1–19 (1995).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comput. Phys. Commun. 271, 108171 (2022).
  * Tribello  _et al._ [2014] G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi, “PLUMED 2: New feathers for an old bird,” Comput. Phys. Commun. 185, 604–613 (2014).
  * The PLUMED consortium [2019] The PLUMED consortium, “Promoting transparency and reproducibility in enhanced molecular simulations,” Nat. Methods 16, 670–673 (2019).
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” J. Mach. Learn. Res. 12, 2825–2830 (2011).
  * Martyna, Tobias, and Klein [1994] G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant pressure molecular dynamics algorithms,” J. Chem. Phys. 101, 4177–4189 (1994).
  * Bussi and Parrinello [2007] G. Bussi and M. Parrinello, “Accurate sampling using Langevin dynamics,” Phys. Rev. E 75, 056707 (2007).
  * Park and Schulten [2004] S. Park and K. Schulten, “Calculating potentials of mean force from steered molecular dynamics simulations,” J. Chem. Phys. 120, 5946–5961 (2004).
  * Ferrarotti  _et al._ [2015] M. J. Ferrarotti, S. Bottaro, A. Pérez-Villa, and G. Bussi, “Accurate multiple time step in biased molecular simulations,” J. Chem. Theory Comput. 11, 139–146 (2015).
  * Bal [2022] K. M. Bal, “Bubble nucleation rate predictions in a Lennard-Jones fluid,” https://www.plumed-nest.org/eggs/22/025 (2022), PLUMED-NEST, plumID:22.025.
  * Yi and Rutledge [2012] P. Yi and G. C. Rutledge, “Molecular origins of homogeneous crystal nucleation,” Annu. Rev. Chem. Biomol. Eng. 3, 157–182 (2012).
  * Tolman [1949] R. C. Tolman, “The effect of droplet size on surface tension,” J. Chem. Phys. 17, 333–337 (1949).
  * Sanchez-Burgos  _et al._ [2020] I. Sanchez-Burgos, P. Montero de Hijes, P. Rosales-Pelaez, C. Vega, and E. Sanz, “Equivalence between condensation and boiling in a Lennard-Jones fluid,” Phys. Rev. E 102, 062609 (2020).
  * Cheng, Dellago, and Ceriotti [2018] B. Cheng, C. Dellago, and M. Ceriotti, “Theoretical prediction of the homogeneous ice nucleation rate: disentangling thermodynamics and kinetics,” Phys. Chem. Chem. Phys. 20, 28732–28740 (2018).
  * Hussain and Haji-Akbari [2020] S. Hussain and A. Haji-Akbari, “Studying rare events using forward-flux sampling: Recent breakthroughs and future outlook,” J. Chem. Phys. 152, 060901 (2020).
  * Haji-Akbari [2018] A. Haji-Akbari, “Forward-flux sampling with jumpy order parameters,” J. Chem. Phys. 149, 072303 (2018).
  * Piaggi and Car [2020] P. M. Piaggi and R. Car, “Phase equilibrium of liquid water and hexagonal ice from enhanced sampling molecular dynamics simulations,” J. Chem. Phys. 152, 204116 (2020).
  * Zou, Tsai, and Tiwary [2021] Z. Zou, S.-T. Tsai, and P. Tiwary, “Toward automated sampling of polymorph nucleation and free energies with the SGOOP and metadynamics,” J. Phys. Chem. B 125, 13049–13056 (2021).
  * Karmakar  _et al._ [2021] T. Karmakar, M. Invernizzi, V. Rizzi, and M. Parrinello, “Collective variables for the study of crystallisation,” Mol. Phys. 40, e1893848 (2021).