Reaction coordinates and rate constants for liquid droplet nucleation: quantifying the interplay between driving force and memory

TL;DR

This study develops a reaction coordinate for liquid droplet nucleation that accounts for local density fluctuations and memory effects, improving rate constant calculations across different supersaturations.

Contribution

The paper introduces a new reaction coordinate combining multiple order parameters, capturing non-Markovian effects at lower supersaturations, enhancing nucleation rate predictions.

Findings

Reaction coordinate shifts from simple atom count to include density fluctuations at lower supersaturation.

Density fluctuations exhibit long memory, affecting the Markovian assumption in nucleation modeling.

Infrequent metadynamics with the new RC yields four orders of magnitude faster and more accurate nucleation rate estimates.

Abstract

In this work we revisit the classic problem of homogeneous nucleation of a liquid droplet in a supersaturated vapor phase. We consider this at different extents of the driving force, which here is the extent of supersaturation, and calculate a reaction coordinate (RC) for nucleation as the driving force is varied. The RC is constructed as a linear combination of three order parameters, where one accounts for the number of liquid-like atoms, and the other two for local density fluctuations. The RC is calculated from all-atom biased and unbiased molecular dynamics (MD) simulations using the spectral gap optimization approach "SGOOP" [P. Tiwary and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A. 113, 2839 (2016)]. Our key finding is that as the supersaturation decreases, the RC ceases to simply be the number of liquid-like atoms, and instead it becomes important to explicitly consider local density fluctuations that correlate with shape and density variations in the nucleus. All three order parameters are found to have similar barriers in their respective potentials of mean force, however, as the supersaturation decreases the density fluctuations decorrelate slower and thus carry longer memory. Thus at lower supersaturations density fluctuations are non-Markovian and can not be simply ignored from the RC by virtue of being noise. Finally, we use this optimized RC to calculate nucleation rates in the infrequent metadynamics framework, and show it leads to more accurate estimate of the nucleation rate with four orders of magnitude acceleration relative to unbiased MD.

Tables

Table 1. Table 1: The metadynamics parameters used at different supersaturation levels: S 𝑆 S represents the supersaturation level, L 𝐿 L is the size of simulation cubic box which we used to control the supersaturation. Gaussian bias kernels of starting height h ℎ h and width ω 𝜔 \omega were added every Δ ​ t Δ 𝑡 \Delta t , which was kept same for metadynamics irrespective of free energy or kinetics calculation. γ 𝛾 \gamma is the bias factor for well-tempered metadynamics. Valsson et al. ( 2016 )

WTmetaD parameters
Label	$S$	$L$ (nm)	$h$ (kJ/mol)	${\omega }_n$	${\omega }_{\chi }$	$\Delta t$ ($ps$)	$\gamma$
$S_1$	13.65	9.9	0.01	0.5	0.08	25	5
$S_2$	12.80	10.1	0.05	0.5	0.08	25	8
$S_3$	11.43	10.5	0.2	0.5	0.08	25	8
$S_4$	9.87	11.0	0.2	0.5	0.08	25	8
$S_5$	9.04	11.3	0.2	0.5	0.08	25	8

Table 2. Table 2: The table shows the characteristic nucleation times τ N subscript 𝜏 𝑁 \tau_{N} by fitting Eq. 4 and the corresponding nucleation rates J 𝐽 J . Results are shown as obtained from the simulations biasing along n 𝑛 n as well as biasing along the optimized RC χ 𝜒 \chi . The labels corresponds to the supersaturation levels denoted in Table 1 . α 𝛼 \alpha is the mean acceleration factor for every set of simulations. For the fitted characteristic nucleation times τ N subscript 𝜏 𝑁 \tau_{N} we have also provided in square brackets the corresponding p-value of the fit when used in Kolmogorov-Smirnov test of Ref. Salvalaglio et al., 2014 .

$S$	RC	${\tau }_N\left(s\right)$ $\left[p-value\right]$	$J\left(1/{\mathrm{cm}}^3/s\right)$	$\alpha$
$S_1$	n	$4.16\pm 0.45\times {10}^{-9}$ $\left[0.17\right]$	$2.48\pm 0.27\times {10}^{26}$	1.1
	$\chi$	$4.56\pm 0.31\times {10}^{-9}$ $\left[0.67\right]$	$2.26\pm 0.15\times {10}^{26}$	1.1
$S_2$	n	$8.66\pm 0.85\times {10}^{-9}$ $\left[0.47\right]$	$1.12\pm 0.11\times {10}^{26}$	2.0
	$\chi$	$7.88\pm 0.51\times {10}^{-9}$ $\left[0.37\right]$	$1.23\pm 0.08\times {10}^{26}$	1.6
$S_3$	n	$1.00\pm 0.15\times {10}^{-7}$ $\left[0.87\right]$	$8.64\pm 1.30\times {10}^{24}$	$4.0\times {10}^1$
	$\chi$	$0.47\pm 0.07\times {10}^{-7}$ $\left[0.35\right]$	$1.84\pm 0.27\times {10}^{25}$	$6.4\times {10}^1$
$S_4$	n	$1.26\pm 0.27\times {10}^{-6}$ $\left[0.64\right]$	$5.96\pm 1.28\times {10}^{23}$	$3.7\times {10}^2$
	$\chi$	$0.69\pm 0.12\times {10}^{-6}$ $\left[0.80\right]$	$1.09\pm 0.19\times {10}^{24}$	$1.1\times {10}^3$
$S_5$	n	$1.58\pm 0.19\times {10}^{-5}$ $\left[0.59\right]$	$4.32\pm 0.52\times {10}^{22}$	$5.7\times {10}^3$
	$\chi$	$0.62\pm 0.15\times {10}^{-5}$ $\left[0.13\right]$	$1.10\pm 0.27\times {10}^{23}$	$7.8\times {10}^3$

Equations

$$c_i=\sum _{j\ne i}\frac{1-(r_{ij}/r_c{)}^6}{1-(r_{ij}/r_c{)}^{12}}$$

$$n=\sum _{i=1}^N\frac{1-(c_l/c_i{)}^6}{1-(c_l/c_i{)}^{12}}$$

$${\mu }_2^2=\frac{1}{N}\sum _{i=1}^N(c_i-\overline{c}{)}^2,{\mu }_3^3=\frac{1}{N}\sum _{i=1}^N(c_i-\overline{c}{)}^3$$

$$P(t)=e^{-t/\tau }$$

$$K_{mn}=\Lambda \sqrt{\frac{{\pi }_n}{{\pi }_m}}$$

$$\Lambda =\frac{⟨N⟩}{\sqrt{{\pi }_m{\pi }_n}}$$

$$\Lambda =\frac{D_{\chi }}{2d^2}$$

$$\alpha (t)=\frac{\tau }{{\tau }_M}=⟨e^{\beta V(t)}⟩$$

Full text

I Introduction

The nucleation of one phase from another is considered as the first step of several phase transitions in chemical physics, with relevance to diverse and important problems in science and technology.Bartels-Rausch (2013); Murray et al. (2010); Mazur (1970); Cox et al. (2007); Erdemir et al. (2007); Sloan (2003); Hammerschmidt (1934); Harper et al. (1997); Cohen et al. (2013) Through experiments, simulations and theory, this problem has been extensively studied over the decades.Sosso et al. (2016) In spite of so much attention being lavished upon this problem, it continues to be a difficult challenge. For instance, the typical size of the critical nucleus is so small that it becomes difficult to observe and analyze it through experiments, or model through continuum based theories. The tens to hundreds of atoms size of the nucleus thus makes it in principle ideal for probing through all-atom molecular dynamics (MD) simulations. However, this is easier said than done due to the inherent rare event nature of the problem, where one nucleation event can take seconds, hours or longer, making it far beyond the microsecond timescale available through the fastest supercomputers. This has led to the development of a plethora of sampling schemes that attempt to enhance the process of nucleation in a controllable manner.Torrie and Valleau (1977); Kumar et al. (1992); Laio and Parrinello (2002); Barducci et al. (2008); van Erp et al. (2003); Moroni et al. (2004); Allen et al. (2006, 2009) These various sampling methods need the pre-determination of slow degree or degrees of freedom relevant to the nucleation process being studied. This slow degree of freedom which is most informative of the underlying physics is referred to as the reaction coordinate (RC).Berezhkovskii and Szabo (2005a); Tiwary and Berne (2016a); Berezhkovskii and Szabo (2005b) In sampling methods such as metadynamics,Barducci et al. (2008); Valsson et al. (2016) where one gradually deposits a time-dependent bias to escape free energy minimum, the need to know a reasonably good RC beforehand is well-documented. In a different class of methods such as forward flux sampling (FFS),Sarupria and Debenedetti (2012) recent work has started to highlight how FFS can benefit from pre-knowledge of adequate slow order parameters or the RC.DeFever and Sarupria (2019) Finally in methods such as transition path sampling (TPS)Dellago et al. (1998); Bolhuis et al. (1998) and variants thereof,van Erp et al. (2003); Moroni et al. (2004) this dependence on pre-knowledge of RC is somewhat mitigated, but instead one becomes reliant on the accuracy of the initial path used in the sampling. In any case, one can say with confidence that any sampling scheme for the study of nucleation can only benefit from a prior sense of an approximate RC for nucleation, with the degree of benefit varying from scheme to scheme.

In this work, we consider what is arguably the simplest of nucleation problems, namely that of the homogeneous nucleation of a liquid droplet in a supersaturated vapor phase at different supersaturation levelsChkonia et al. (2009). The system is modeled using Lennard-Jones interactions.Salvalaglio et al. (2016); Chkonia et al. (2009) Even in this simplest of problems, we find that the RC for homogeneous nucleation deviates significantly from standard assumptions made so far in theoretical and simulation approaches.Salvalaglio et al. (2016); Kuipers and Barkema (2009); Ford (2004) Our calculations of the RC are performed using a spectral gap based optimization method “SGOOP”, originally proposed by Tiwary and Berne Tiwary and Berne (2016a, 2017a); Tiwary and Berne (2016b); Smith et al. (2018), for the automatic construction of RC from different trial order parameters. We find that there exists a supsersaturation (or more generally speaking, the driving force for nucleation) dependent interplay between size, density and shape of the nucleus. This interplay leads to a non-trivial RC that goes far beyond a spherical, uniformly dense nucleus assumed in classical nucleation theory (CNT).Weber and Volmer (1926); Farkas (1927); Becker and Döring (1935); Zeldovich (1943); Maris (2006); Kalikmanov (2013) While we define RC more rigorously in the main text, here we summarize it as a low-dimensional variable permitting a Markovian description of the underlying high-dimensional dynamics.Berezhkovskii and Szabo (2011a, 2005a); Kuipers and Barkema (2009) Our key finding is that as the supersaturation decreases, the RC becomes composed of not just the number of atoms in the largest liquid-like cluster, but it also becomes helpful to consider the spatial fluctuations of the aforementioned quantity. These fluctuations display similar barriers as the number of liquid-like atoms, but have a longer autocorrelation time (or equivalently, slower diffusion). This diffusion anisotropy becomes stronger as the supersaturation decreases. In the spirit of works by Szabo, Peters, Hynes and others,van der Zwan and Hynes (1982); Berezhkovskii and Szabo (2005a); Hynes (1985); Peters et al. (2013) we find that the RC itself starts to align with the direction of slowest diffusion or longest memory, given that the free energy barriers in the directions of various individual order parameters are similar. Finally, we use the optimized RC as a biasing variable in infrequent metadynamics calculations,Tiwary and Parrinello (2013) which allow recovering unbiased kinetic rate constants from biased simulations. We find that considering this diffusion anisotropy adjusted RC in infrequent metadynamics leads to more accurate estimates of the nucleation rate across different supersaturations with orders of magnitude speed-up relative to unbiased MD.

We believe that our work demonstrates the potential of using methods such as SGOOP in unraveling the subtle aspects of the RC in complex nucleation problems. Such a RC first of all directly gives useful physical insight into the processes at play, but secondly, as we show here it also serves as useful descriptor for performing enhanced sampling simulations including metadynamics and beyond.

II Theory and Methods

II.1 Order parameters

In order to motivate this work and the various order parameters we consider here, we start with a brief description of CNT which has been a basic building block in the study of nucleation. In CNT, the first liquid droplet formed in the vapor is treated as spherical shaped and uniformly dense.Weber and Volmer (1926); Farkas (1927); Becker and Döring (1935); Zeldovich (1943); Maris (2006); Kalikmanov (2013) The nucleation process is then modeled by balancing the surface tension penalty with chemical potential benefit.Rao et al. (1978) This simple theory, though it captures qualitatively how nucleation happens, however fails to quantify the true nucleation rate in any practical sense. It is believed that CNT makes several oversimplified assumptions especially incorrectly assuming that the cluster is spherical and uniform.Wang et al. (2007); Kalikmanov et al. (2008) By using numerical and experimental tools, the lack of sphericity and uniform density has indeed been documented in crystal nucleation and in nucleation in more complex systems.ten Wolde and Frenkel (1997); Trudu et al. (2006); Moroni et al. (2005); Kwon et al. (2015); Zhou et al. (2019) However such simulations and experiments are expensive, and it has been hard to quantitatively probe such effects even in the simple gas system such as the one used in this work.

A popular order parameter that goes beyond the spherical nucleus approximation of CNT was introduced by Frenkel and ten Wolde.ten Wolde and Frenkel (1998) This order parameter $n$ equals the number of liquid-like atoms in the system in a way that it is still a continuous and differentiable functionten Wolde and Frenkel (1998) of atomic coordinates, a necessity for the biased simulations we perform later. In this definition, an atom is classified as liquid if it has more than 5 neighboring atoms. The number of the neighborhood atoms of the atom with label $i$, or equivalently the coordination number $c_{i}$, is calculated through the use of a switching function as follows:

[TABLE]

where the summation is carried over all atoms $j\neq i$. The distance $r_{ij}$ between atoms $i$ and $j$ needs to be less than a cut-off $r_{c}$ to be considered as neighbors. The number of liquid phase atoms $n$ is then calculated using a similar form with threshold value $c_{l}$ which we take to be 5 here in spirit of Ref. ten Wolde and Frenkel, 1998:

[TABLE]

The above defined $n$ thus captures the total number of liquid-like atoms in the systemten Wolde and Frenkel (1998). It is however oblivious to details such as the density of the clusters in which these atoms are present, if there are more than 1 clusters, the shape of these clusters and other nuances. In order to consider these, we propose including the second and third moments of the distribution of coordination numbers, defined as $\mu^{2}_{2}$ and $\mu^{3}_{3}$ respectively. Moroni et al. (2005); Tribello et al. (2017, 2010) These are explicitly calculated as:

[TABLE]

where $\overline{c}$ is the average of the coordination number. In Fig. 1 we show a representative unbiased MD trajectory in $(n,\mu_{2}^{2})$ space where it can be seen that at roughly the critical size of $n=30$ different $\mu^{2}_{2}$ can represent clusters with strikingly different profiles in terms of shape, density and compactness.

The distribution of the coordination numbers captures the variations in density of the liquid phase and thus can be used to study the local properties of the liquid droplets. Another nice feature of these moments is that they can easily be generalized for multi-component systems to take into account the variation in density related to specific species.ten Wolde and Frenkel (1997); Kwon et al. (2015); Kathmann et al. (2004) Our RC is then expressed as a linear combination of these three order parameters $n$, $\mu^{2}_{2}$ and $\mu^{3}_{3}$.

II.2 Nucleation rate

The process of nucleation is inherently stochastic in nature and satisfies the law of rare events. In other words, different independent observations of nucleation should give a distribution of nucleation times adhering to a Poisson process.Resnick (2013) If we let $P(t)$ denote the survival probability of not observing any liquid droplets at and until time $t$, it will satisfy the following relation valid for all Poisson processes:

[TABLE]

where $\tau$ is the characteristic time for the first nucleation event, the inverse of which can be interpreted as the nucleation rate. In this work, we find the characteristic time by performing multiple independent simulations starting from the system in gaseous state with randomized velocities (other simulation details in Sec. II.3), and collect the statistics of transition times until the first nucleation event. The nucleation time is then obtained by performing a Poisson fit to these independent observations following the protocol outlined in Ref. Salvalaglio et al., 2016, 2014.

II.3 Simulation set-up

The simulations were performed under the constant number, volume, temperature (NVT) ensemble with N=512 argon atoms and average temperature fixed at 80.7 K. The volume of the simulation box was set in order to correspond to desired supersaturation levels $S$ detailed in Table 1. In order to compare our results with unbiased nucleation rates in Ref. Chkonia et al., 2009, our cubic box size ranged from 9.5 nm to 11.5 nm. The interaction between atoms were modeled through a Lennard-Jones potential with $\epsilon=0.99797$ kJ/mol and $\sigma=0.3405$ nm Chkonia et al. (2009). The potential was truncated with cutoff at 6.75 $\sigma$. The velocity rescale thermostat with time constant of 0.1 $ps$ was used to do temperature coupling. Bussi et al. (2007) All simulations were performed using GROMACS version 2016.5Lindahl et al. (2001) patched with PLUMED version 2.4.2.Massimiliano Bonomi (2019)

II.4 Reaction coordinate

In order to quantify how these various order parameters $n$, $\mu^{2}_{2}$, and $\mu^{3}_{3}$ matter for the process of nucleation, we intend to learn a RC $\chi$ as their linear combination. In addition to quantifying exactly how much these order parameters matter for driving nucleation, this RC will also serve as a crucial input for biased simulations to be performed later in this work. We first carefully define what exactly we mean by RC.

The RC for a given molecular system is traditionally defined as an abstract low-dimensional coordinate that best captures progress along relevant reaction pathway. While this intuitive notion can be formalized and quantified in several different ways, here we use the definition of RC as follows. For a given multidimensional complex system undergoing a certain dynamics, it is an optimal low-dimensional variable such that the multidimensional dynamics of the full system in terms of movement between different metastable states can be mapped into Markovian dynamics between various states viewed as a function of the RCBerezhkovskii and Szabo (2005a, 2011b). Thus an optimal RC is a low-dimensional mapping which best satisfies (i) thermodynamic truthfulness: demarcating between the various relevant metastable states present in the actual high-dimensional system, (ii) kinetic truthfulness: preserving pathways for moving between these different states, and (iii) timescale separation: displaying a clean-cut separation of timescales between the relaxation times in the various metastable states, and the time spent in the actual event of crossing from one state to another.

II.5 SGOOP

To find such a RC, here we use the method “Spectral gap optimization of order parameters (SGOOP)”. This method uses the principle of maximum caliber (“MaxCal”), which is similar to path entropy,Dixit et al. (2015); Jaynes (1980); Pressé et al. (2013); Dixit et al. (2018) to construct a transition probability matrix along any candidate RC, and then calculates its eigenvalues $\lambda_{0}=1>\lambda_{1}>\lambda_{2}>...$. Here $\lambda_{0}=1$ corresponds to stationary state, while the other eigenvalues carry information about the timescales of various dynamical processes. The best RC will then produce a transition matrix $K$ with a maximal timescale separation between visible slow and hidden fast processes. This timescale separation, also known as spectral gap, is quantified as the difference $\lambda_{n}-\lambda_{n+1}$, where $n$ is the number of discernible energy wells along the putative RC. SGOOP needs two key inputs: (i) an estimate of the stationary probability density $\pi$ along any putative RC, and (ii) some dynamical observables or constraints. With these inputs, the SGOOP transition matrix $K$ can be formulated as follows:

[TABLE]

where $\pi_{m}$ is the stationary probability along any putative, spatially-discretized RC $\chi$ with $m$ denoting the grid index, and $\Lambda$ is a dynamical observable we will revisit shortly. $K_{mn}$ gives the rate for moving from grid $m$ to grid $n$ in a small time interval. The input (i), namely the stationary density $\pi$ can come from unbiased MD at high enough supersaturations, or if the supersaturation is too low to permit unbiased MD, it can come through the use of preliminary metadynamics along a trial RC, followed by reweighting.Tiwary and Parrinello (2014) For input (ii), namely calculation of the dynamic observable needed to constrain the maximum caliber estimate of rate matrix, we run short unbiased MD runs which calculate the mean number of nearest neighbor transitions $\langle N\rangle$ along any putative RC. It is then easy to showTiwary and Berne (2017b) that the dynamical observable $\Lambda$ in Eq. 5 is given by

[TABLE]

EquivalentlyTiwary and Berne (2017b), if one was to completely by-pass the MaxCal framework, a very similar equation as Eq. 5 can be derived by comparing a master equation along $\chi$ with a discretized Smoluchowski equation along the same. Bicout and Szabo (1998) Then the prefactor $\Lambda$ becomes:

[TABLE]

where $D_{\chi}$ is the position-dependent diffusivity along the coordinate $\chi$ and $d$ is the grid spacing along $\chi$.

Eqs. 5–7 collectively show that the rate matrix $K$, hence the spectral gap, and consequently the optimized RC, depend not just on the free energy barriers that would be encapsulated in the stationary density $\pi$, or equivalently in the associated free energy, but that the dynamics of the system as captured in the diffusivity of the various order parameters can also play a significant role in the RC. As we will show later in Sec. III, we find this to be a very important point in the context of liquid droplet nucleation.

II.6 Metadynamics

For high enough supersaturation such as $S_{1}$, $S_{2}$, $S_{3}$ in Table 1, we can perform unbiased simulations directly in reasonable computer time, both for the calculation of nucleation kinetics and for feeding stationary density into SGOOP for constructing the RC. However, for lower supersaturations we need to apply enhanced sampling methods since nucleation becomes a rare event. In this work, we use well-tempered metadynamicsBarducci et al. (2008); Valsson et al. (2016) along a trial RC to obtain preliminary stationary density estimates, and infrequent metadynamicsTiwary and Parrinello (2013, 2014) to calculate the kinetics of nucleation.

In metadynamics,Laio and Parrinello (2002); Barducci et al. (2008); Valsson et al. (2016) the system is encouraged to visit new states by adding a history-dependent Gaussian bias $V(s,t)$ as a function of a biasing variable $s$. Here specifically we use the well-tempered variant of metadynamics in which the height of the Gaussian is tempered through a bias factor each time a point is revisited. The different parameters of the Gaussian bias are listed in Table 1 in which $h$ is the starting height, $\omega$ is the width, $\Delta t$ is the deposition interval, and $\gamma$ is the bias factor. The final output of a traditional metadynamics run is the free energy along the variable $s$ or along any other degree of freedom which can be expressed as function of the atomic coordinates of the system. In theory, a relation connecting the free energy with the deposited bias can be derived irrespective of the precise choice of biasing variable $s$, which will be asymptotically valid in the limit of long simulation time. In practice however, it helps if $s$ is as close to the true RC as possible.

More recently, a simple extension to well-tempered metadynamics was introduced which allows recovering not just static free energies but also unbiased kinetic information from metadynamics. This protocol has been dubbed “infrequent metadynamics”.Tiwary and Parrinello (2013, 2014) The key idea here is that as long as the bias deposition rate is infrequent enough compared to barrier-crossing timescales, in principle we should be able to reweight the biased timescales from well-tempered metadynamics directly to obtain unbiased kinetics through a simple acceleration factor:

[TABLE]

where $\tau$ is the unbiased transition time we seek to learn and $\tau_{M}$ is the biased transition time we actually observe in metadynamics. $V(t)$ is the net bias deposited until time $t$. The central assumption in infrequent metadynamics is that the biasing variable does a good job of timescale separation between time spent in the free energy basin and the time spent during barrier crossing. Thus the need to have a more accurate RC for biasing, that satisfies the criteria of Sec. II.4 becomes even more significant for infrequent metadynamics than in traditional metadynamics. For instance, as we will show in the Results section, infrequent metadynamics becomes more accurate if the biased variable includes all relevant slow modes with long autocorrelation time, and any hidden modes not considered in the biasing variable are as markovian (or quickly decorrelating) as possible. Here we learn such a 1-dimensional RC $\chi$ as a linear combination of our order parameters. Our bias potential then becomes $V(\chi,t)=V(w_{1}n+w_{2}\mu^{2}_{2}+w_{3}\mu^{3}_{3},t)$. The weights of the different order parameters $(w_{1},w_{2},w_{3})$, are determined with SGOOP Tiwary and Berne (2016c).

III Results

III.1 RC predicted from SGOOP

We first describe the RC, as introduced and defined in Sec. II.4, that we identify for the condensation of a liquid droplet across different supersaturation values. To learn this RC we have used SGOOP, Tiwary and Berne (2016a, 2017a) with the stationary probability density $\pi$ estimated through unbiased MD at high supersaturations $S_{1}$, $S_{2}$ and preliminary metadynamicsValsson et al. (2016) at low supersaturations $S_{4}$, $S_{4}$, $S_{5}$. The preliminary metadynamics runs were performed biasing $n$. All runs were complemented with short unbiased MD runs (50 ns) for obtaining dynamical constraints for MaxCal. For every supersaturation, SGOOP is initiated from a given choice of trials weights ($w_{1},w_{2},w_{3}$) for the RC $\chi$ expressed as $\chi=w_{1}n+w_{2}\mu^{2}_{2}+w_{3}\mu^{3}_{3}$.

A key question that immediately arises is whether at any given supersaturation $S$ there is a unique RC, or if there are multiple possible combinations of the weights ($w_{1},w_{2},w_{3}$) which meet equally well the criteria for an optimal RC described in Sec. II.4. Yet another question which we ask and answer is how transferable is the RC learnt at one supersaturation $S$ across different values of $S$. To answer the first question, we perform several exhaustive SGOOP trials to estimate the optimized RC, first in the 2-d $(n,\mu_{2}^{2})$ space where we do an explicit grid based search over the full space, and then in the 3-d $(n,\mu_{2}^{2},\mu_{3}^{3})$ space where we start SGOOP from different initial weights. In the latter case, the optimization over weights in SGOOP is performed using a basin hopping algorithm which is a global search algorithm with several stochastic jumps aiding the system from not getting trapped in local minima.

In the 2-d $\chi=w_{1}n+w_{2}\mu_{2}^{2}$ optimization, we do an explicit search among all possible RCs by rotating the putative RC $\chi\equiv\cos(\theta)n+\sin(\theta)\mu^{2}_{2}$ in the $(n,\mu_{2}^{2})$ space. Here as shown in Fig. 2(a) for $S$ = 11.43, we find that the spectral gap profile has a sharp peak when the RC is almost exclusively comprised of $\mu_{2}^{2}$, i.e. $\theta\approx\pi/2$ and $\mu_{2}^{2}$ has around 8 times higher weight in the RC than $n$. Such a $\mu_{2}^{2}$–heavy RC is obtained irrespective of any $S$ value, showing unequivocally that the second moment $\mu_{2}^{2}$ plays a more important role in the RC than $n$ itself. Fig. 2(a) also shows the variation of the kinetic pre-factor $\Lambda$ of Eq. 5 with RC choice, and we will revisit this profile in Sec. III.2. Next, we perform optimization in the full 3-d $(n,\mu_{2}^{2},\mu_{3}^{3})$ space. Here we find that there are many combinations $(w_{1},w_{2},w_{3})$ with similarly enhanced spectral gaps relative to the traditional choice of $\chi=n$, but with a common theme that the second moment and the third moment consistently show up in the optimized RC.

Thus to summarize so far: (a) RC optimization in $(n,\mu_{2}^{2})$ leads to RC predominantly comprised of $\mu_{2}^{2}$, (b) RC optimization in $(n,\mu_{2}^{2},\mu_{3}^{3})$ leads to a RC invariably with weights for all 3 variables, but with multiple local maxima in the spectral gap profile. In other words, the RC is quite degenerate, but considering $\mu_{2}^{2}$ and $\mu_{3}^{3}$ in the RC is important for a more accurate description of the nucleation process. In the SI we also show results from a full grid search over spectral gaps in the ($w_{1},w_{2},w_{3}$) space at at $S$ = 11.43 further illustrating the findings from SGOOP. Here among the first few largest local maxima from 3 different trajectories, we picked ($w_{1},w_{2},w_{3}$) = $(0.15,0.65,-0.15)$ for use in further calculations across all supersaturations $S$. In Fig. 2 (b), we plot the ratio between the spectral gap along RC = $0.15n+0.65\mu^{2}_{2}-0.15\mu^{3}_{3}$ and that along RC = $n$ at different $S$. As can be seen there, at all $S$ values the optimized RC gives higher spectral gaps than just $n$, and the improvement increases sharply with decreased supersaturation. That is, as the supersaturation decreases the importance of consider shape and density fluctuations in the nuclei become more and more important, which is one of the central findings of this paper. Furthermore, the optimized RC learnt at one supersaturation gives improved spectral gaps at other supersaturations, and hence the RC is transferable across supersaturations. Thus in Sec. III.3 we use the RC $\chi=0.15n+0.65\mu^{2}_{2}-0.15\mu^{3}_{3}$ at all supersaturations for enhanced sampling based calculations of the nucleation rate.

III.2 Understanding the RC

SGOOP optimizes the RC by finding a low-dimensional projection with highest gap between slow and fast processes. In most cases this amounts to selecting a projection with the highest barrier separating the metastable states. To understand if the RC learnt in Sec. III.1 can be attributed to simply barriers in the free energy profile, or if dynamical concerns such as the prefactor $\Lambda$ in Eqs. 5–7 play a role, we construct free energies along various 1-d and 2-d components (totaling 6 combinations) of $(n,\mu_{2}^{2},\mu_{3}^{3})$. These free energies were obtained by running metadynamics with same parameters defined in Table. 1 and bias potentials added along $n$. We then averaged over 10 independent metadynamics runs with each trajectory reweighted using the free estimator described in Ref. Tiwary and Parrinello, 2014.

From the various 1-d and 2-d free energy profiles shown in Fig. 3 (a)-(f) for $S$ = 11.43, it is hard to distinguish between the importance of the various order parameters $n$, $\mu_{2}^{2}$,and $\mu_{3}^{3}$. The 2-d profiles show that starting from the gas phase (red stars in Fig. 3 (a)-(c)), all three order parameters change in a very correlated manner until the barrier is reached and nucleation is essentially complete ($n$ >100). The 1-d free energies along the three order parameters (Fig. 3 (d)-(f)) show that the free energy barrier that needs to be overcome is also very similar for each of the 3 order parameters, though there are some systematic differences which we revisit shortly in Fig. 3 (g)-(i). Comparing the 1-d free energy along $\mu^{3}_{3}$ (Fig. 3 (f)) with the corresponding 2-d free energies (Fig. 3 (b)-(c)), we can see that unlike $n$ and $\mu_{2}^{2}$, the 1-d projection along $\mu_{3}^{3}$ does a very poor job of describing the pathway in higher dimension space, further justifying our choice of RC $\chi$ in the previous section with higher weight for $\mu_{2}^{2}$ than for $\mu_{3}^{3}$. In Fig. 3 (g)-(i), we show the free energies along three different RC choices, namely $n$, $\mu_{2}^{2}$ and the optimized $\chi=0.15n+0.65\mu^{2}_{2}-0.15\mu^{3}_{3}$, for three different $S$ values. As $S$ is decreased, invariably there is a small but consistent improvement in the barrier height when viewed as function of $\chi$ or $\mu^{2}_{2}$, relative to when viewed as function of $n$. However, firstly this difference is very small (0.25 $kJ$ or 0.1 $k_{B}T$), and secondly, it appears to get even smaller with decreasing $S$ (Fig. 3 (g)-(i), left to right). Thus, the free energy barrier can not be used to explain the behavior of spectral gap versus supersaturation shown in Fig. 2 (b). Here we showed that at all supersaturation levels we considered, the spectral gaps of the optimized RC are better than of $n$. It was also pointed out that as supersaturation decreases the spectral gap improvement increases. This tells us that the optimized RC works better at lower supersaturation, which is inconsistent with the change in free energy barriers along different order parameters with supersaturation.

Our next step is therefore explaining why SGOOP finds that $\mu^{2}_{2}$ has a role to play in $\chi$, and why the advantage in considering $\mu^{2}_{2}$ increases with decreasing supersaturation $S$. In Fig. 2(a), we provided a profile of how the prefactor $\Lambda$ varied with the RC choice and correspondingly with the spectral gap. It can be seen there that the prefactor $\Lambda$ has a strong inverse correlation with the spectral gap of $\chi$ – the maximum spectral gap coincides with minimum $\Lambda$. Thus the minuscule increase in barrier height with varying RC is compensated by the slowness of the dynamics along the RC, as captured by $\Lambda$ or the average number of first neighbor transitions in a unit time.

To gain further insight into this, we calculated time-autocorrelation functions along our three different order parameters (see Fig. 4) as higher autocorrelation time represents slower diffusivity. Our calculations show that $\mu^{2}_{2}$ and $\mu^{3}_{3}$ have longer autocorrelation times than n, and therefore lose memory slower than $n$.Kuipers and Barkema (2009); Ford (2004) Furthermore, the increase in autocorrelation times of the two order parameters $\mu^{2}_{2}$ and $\mu^{3}_{3}$ relative to $n$ becomes more and more pronounced as the supersaturation $S$ decreases (Fig. 4). This is in striking contrast to Fig. 3, where we found an opposite trend looking at the free energy barriers along these order parameters.

We therefore conclude this section with the observation that anisotropic diffusion in the space of order parameters becomes an important factor in determining the RC especially at lower supersaturations. The longer autocorrelation times are linked to less Markovian behavior, which means $\mu^{2}_{2}$ and $\mu^{3}_{3}$ carry longer memory than $n$.Kuipers and Barkema (2009); Ford (2004) Coupled with the finding that all three order parameters have similar barriers in their respective potential of mean force, this means that change in nuclei characteristics such as shape and density become slower as supersaturation $S$ decreases, and it becomes important to explicitly consider this in the construction of a Markovian RC.

III.3 Nucleation kinetics

Now that we have identified an optimized RC $\chi$ = $0.15n+0.65\mu^{2}_{2}-0.15\mu^{3}_{3}$ with improved spectral gap relative to the Frenkel-ten Wolde parameter $n$, we perform two sets of enhanced sampling simulations (specifically, infrequent metadynamics) using $n$ and $\chi$ as biasing variable respectively. We use Eq. 8 to reconstruct the unbiased timescale estimates from these biased runs. At high enough supersaturations we are able to run unbiased MD as well and together with the results of Reguera et alChkonia et al. (2009) these constitute a valuable set of results to benchmark our findings against. At each supersaturation level, we launched 40 independent metadynamics runs with 20 of them biasing $n$ and the other 20 biasing the optimized RC $\chi$. For each independent run, in order to be able to compare our results with previously published work Chkonia et al. (2009); Salvalaglio et al. (2016) we defined the nucleation event as when the coordination number $n$ reaches 30 for the first time. Every independent observation of such an event in terms of its metadynamics time was scaled by the acceleration factor (Eq. 8) to obtain an unbiased observation of the nucleation time. With these 20 independent estimates of the nucleation time, we can compute the characteristic time (Eq. 4) of observing the first nucleation event $\tau_{N}$ by fitting a Poisson distribution to the statistics, where $\tau_{N}$ is the expected value of the fitted Poisson distribution. The corresponding nucleation rates are then calculated through the formula $J=1/(\tau_{N}V)$ and $J=1/(t_{N}V)$ where V is the volume of system. The results are shown in Table 2 and in Fig. 5.

We find that the use of $n$ as a biasing variable in infrequent metadynamics does a remarkably decent job of obtaining nucleation rates (in agreement with the findings of Ref. Salvalaglio et al., 2016) even with very significant acceleration factors or computational boost relative to unbiased MD. There is nonetheless further improvement that can be obtained in the quality of the nucleation rate if the optimized RC is used instead of $n$, especially as the supersaturation is decreased. As can be seen from Table 2, the acceleration factor in metadynamics relative to unbiased MD increases steadily as $S$ decreases, reaching almost four orders of magnitude at the lowest $S$. All reweighted nucleation times, irrespective of whether they came from biasing $n$ or biasing $\chi$ give p-values above the recommended cut-off in the Kolmogorov-Smirnov test from Ref. Salvalaglio et al., 2014. At $S=11.43$, the use of $\chi$ as biasing variable instead of $n$ leads to much better agreement with the unbiased estimate of Reguera et al Chkonia et al. (2009), as can be seen in Fig. 5. In general, the characteristic times for nucleation from runs biasing the RC $\chi$ are significantly lower than that those from biasing $n$, and roughly speaking this difference increases as $S$ decreases. In addition to the explicit agreement with unbiased estimate of Reguera et al Chkonia et al. (2009) at $S=11.43$, the lower characteristic time (with similar p-values) can be seen as further evidence of the benefit of biasing $\chi$ instead of $n$. This is because in metadynamics the presence of missing slow degrees of freedom from explicit consideration in the biasing variable typically leads to hysteresis during free energy calculations, or overestimate of the accelerated time through Eq. 8, as pointed out in Ref. Salvalaglio et al., 2014 and Ref. Casasnovas et al., 2017.

IV Discussion

In this work, we used new tools Tiwary and Berne (2016a); Tiwary and Parrinello (2013) to revisit a classic problem in nucleation, namely that of the formation of liquid droplet from gaseous precursor as a function of varying driving force for nucleation, namely supersaturation. Our interest was in (a) constructing a Markovian reaction coordinate (RC) for this process, and (b) testing if there is any gain to be had through the use of a more Markovian RC in enhanced sampling calculations of nucleation kinetics. To answer these questions especially at low supersaturations where access to unbiased trajectories of nucleation is difficult (needed by many other RC optimization methods such as Ref. Peters and Trout, 2006 and Ref. Best and Hummer, 2005), we use the spectral gap optimization method from Ref. Tiwary and Berne, 2016a to construct optimized RC from input biased simulations. Our calculations demonstrate unequivocally that it is not sufficient to consider only the typical order parameter used to describe nucleation, namely the number of liquid like atoms in the system. By considering further variables that account for heterogeneity in the system, such as higher moments $\mu^{2}_{2}$ and third moment $\mu^{3}_{3}$ of the distribution of coordination numbers, we could obtain a much more Markovian RC. Interestingly these various order parameters have nearly identical free energy barriers, and they differ primarily only in associated diffusivities. The importance of these variables further increases with decreasing supersaturation as their associated autocorrelation time increases sharply. In other words, shape and density fluctuations in the nucleating clusters cease to stay rapidly equilibrating variables which can be entirely ignored from a Markovian low-dimensional description of nucleation. We conclude that diffusion anisotropy plays a more important role at lower $S$, which is supported by our analysis of autocorrelation functions and autocorrelation times. While previous work has demonstrated how infrequent metadynamics can predict nucleation time with only $n$ as the RC, we show in this work that the prediction of nucleation time can be further improved by biasing along an optimized RC. It will be interesting to see if the use of such a more Markovian RC makes improvement in the reliability and efficiency of other enhanced sampling methods such as forward flux sampling.

Finally, it should be mentioned that in this calculation we didn’t consider effects due to the finite size of the system, which can be done using the method proposed in Ref. Salvalaglio et al., 2016, but was not the main objective here. Similarly our findings might change with constant number, pressure and temperature (NPT) simulations. The present work can be redone taking these important nuances into account. Finally, strictly speaking ours was a model system with model parameters. This work is a proof of principle that ideas such as SGOOP for RC optimization are potentially useful for study of nucleation through enhanced sampling or otherwise. In future we will be extending this work to systems such as crystal nucleation, multiple polymorphs, systems with multiple pathways or multiple species, where there will be even more order parameters to be considered. All of these continue to be very difficult yet important problems for understanding nucleation pathways and rates, and we are hopeful our tools will allow us and others to systematically study these.

Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for partial support of this research (PRF 60512-DNI6). We also thank Deepthought2, MARCC and XSEDE (projects CHE180007P and CHE180027P) for computational resources used in this work. SGOOP code is available at https://github.com/zwsmith200/SGOOP/