# Free energy calculations along entropic pathways: II. Droplet nucleation in
binary mixtures.

Caroline Desgranges and Jerome Delhommelle  Department of Chemistry,
University of North Dakota, Grand Forks ND 58202

###### Abstract

Using molecular simulation, we study the nucleation of liquid droplets from
binary mixtures and determine the free energy of nucleation along entropic
pathways. To this aim, we develop the $\mu_{1}\mu_{2}VT-S$ method, based on
the grand-canonical ensemble modeling the binary mixture, and use the entropy
of the system $S$ as the reaction coordinate to drive the formation of the
liquid droplet. This approach builds on the advantages of the grand-canonical
ensemble, which allows for the direct calculation of the entropy of the system
and lets the composition of the system free to vary throughout the nucleation
process. Starting from a metastable supersaturated vapor, we are able to form
a liquid droplet by gradually decreasing the value of $S$, through a series of
umbrella sampling simulations, until a liquid droplet of a critical size has
formed. The $\mu_{1}\mu_{2}VT-S$ method also allows us to calculate the free
energy barrier associated with the nucleation process, to shed light on the
relation between supersaturation and free energy of nucleation, and to analyze
the interplay between the size of the droplet and its composition during the
nucleation process.

##  I Introduction

The nucleation of liquid droplets is a ubiquitous phenomenon, central to many
applications in chemistry, physics and atmospheric sciences Yasuoka and
Matsumoto (1998); Oxtoby (1992); Shen and Debenedetti (1999); Weakliem and
Reiss (1993); Schenter et al. (1999); Zeng and Oxtoby (1991a); Yi et al.
(2002); Kinjo et al. (1999); Toxvaerd (2001); Ford (1996); Talanquer and
Oxtoby (1995a); Reiss et al. (1990); Kalikmanov and Van Dongen (1995); Horsch
et al. (2008); Neimark and Vishnyakov (2005); Oxtoby and Evans (1988); Lutsko
(2008); Wang et al. (2008); Ten Wolde et al. (1999); Gonzalez et al. (2015);
Loeffler et al. (2015); Sosso et al. (2016); Xu et al. (2015); Keasler and
Siepmann (2015); Wilhelmsen et al. (2015); van Meel et al. (2015); Hale (1986,
2005); Hale and Thomason (2010); Yuhara et al. (2015); Lauricella et al.
(2015); Singh and Müller-Plathe (2014); Ni and Dijkstra (2013); Reinhardt and
Doye (2014); Ten Wolde and Frenkel (1998); Chen et al. (2001); Oh and Zeng
(1999); Chen et al. (2002); Zhukhovitskii (1995); Nishi et al. (2015); Lupi et
al. (2016); Santiso and Trout (2015); Berryman et al. (2016); Zimmermann et
al. (2015); Lam et al. (2015); Kratzer and Arnold (2015); Bolhuis and Dellago
(2015); Lau et al. (2015); Toxvaerd (2016). Nucleation from a single component
system can be rationalized in terms of a change in a single intrinsic
variable, such as e.g. the chemical potential or pressure. In this case, the
supersaturation of the parent vapor phase is simply characterized by a given
value of the pressure Ten Wolde and Frenkel (1998); Tanaka et al. (2005);
Kraska (2006); Oh and Zeng (2000); Senger et al. (1999); Lau et al. (2015);
Toxvaerd (2016) (or equivalently of the chemical potential) that departs from
the pressure at coexistence. On the other hand, nucleation from a mixture
involves a myriad of pathways arising from the larger dimension of the system,
with as additional intrinsic variables, the mole fractions for each of the
components Kulmala and Laaksonen (1990); Zeng and Oxtoby (1991b); Oxtoby and
Kashchiev (1994); Napari and Laaksonen (1999); Jaecker-Voirol and Mirabel
(1988); Talanquer and Oxtoby (1995b); ten Wolde and Frenkel (1998); Laaksonen
and Oxtoby (1995); Yoo et al. (2001); Napari and Laaksonen (2000); Braun et
al. (2014); Shimizu and Tanaka (2015); Pinho et al. (2014); Gao et al. (2014);
Alekseechkin (2015); Watson et al. (2011); Desgranges and Delhommelle (2014a).
This also makes the definition of an appropriate reaction coordinate for the
system especially challenging since, for instance, the total number of
particles in the cluster is not the only significant variable, as the numbers
of particles of each type also need to be taken into account to fully
characterize the nucleation process. Here, we propose and implement a new
approach that uses the entropy $S$, which captures the interplay between size
increase and molecular selectivity during nucleation, as the reaction
coordinate for the process.

The aim of this work is to shed light on the entropic pathways followed during
the nucleation of two binary mixtures. The first example we consider is a
binary mixture of two highly miscible gases, specifically the $Ar-Kr$ mixture.
The second example involves a mixture of carbon dioxide with an alkane (here,
as an example we consider $C_{2}H_{6}$) that is of technological relevance for
the oil industry and for separation applications. We develop here the
$\mu_{1}\mu_{2}VT-S$ method, where $\mu_{1}$ and $\mu_{2}$ are the chemical
potentials for the two components of the mixture and $V$ and $T$ are the
volume and temperature of the system, to simulate the nucleation of a liquid
droplet for these two binary mixtures. In the first paper of the series
Desgranges and Delhommelle (2016a), we discussed how, in the case of a single
component system, the $\mu VT-S$ approach provided a direct connection with
classical nucleation theory McGraw and Laaksonen (1996), as any arbitrary
value of $\mu$ (or, equivalently, of the supersaturation $\Delta\mu$) could be
applied, and allowed the system to overcome the free energy barrier of
nucleation. In practice, this was achieved by using the entropy of the system,
which can be readily calculated in the grand-canonical ensemble, as the
reaction coordinate and by driving the system along an entropic pathway using
the umbrella sampling simulation technique. Extending this approach to the
case of mixtures is especially appealing since the grand-canonical ensemble
naturally allows the number of molecules of each component to vary as
nucleation proceeds. It is therefore very well suited to shed light on the
impact of the choice of a given supersaturation (i.e. $\Delta\mu_{1}$ and
$\Delta\mu_{2}$ for a binary mixture) on the free energy barrier of nucleation
as well as on the interplay between the size of the droplet and its
composition as nucleation takes place.

The paper is organized as follows. In the next section, we present the
simulation method as well as the molecular models used in this work. We
explain how we set up the $\mu_{1}\mu_{2}VT-S$ simulations, detailing how we
proceed with the choice of chemical potentials, supersaturations, range of
entropies to be sampled and the relation between the conditions of nucleation
and the composition of the bulk. We then discuss the results obtained during
the simulations of the nucleation process from supersaturated vapor phases of
$Ar-Kr$ and $CO_{2}-C_{2}H_{6}$. For all systems, we determine the free energy
profile of nucleation and show how the choice of the conditions of nucleation
impacts the height of the free energy barrier. We also focus on the analysis
of the nucleation mechanism and on unraveling the interplay between size and
composition during the formation of the liquid droplet, before finally drawing
the main conclusions from this work in the last section.

##  II Simulation Method

###  II.1 Spanning entropic pathways

We extend to the case of mixtures the simulation method developed in the first
paper in this series Desgranges and Delhommelle (2016a). The approach proposed
is termed as $\mu VT-S$ and consists in sampling configurations of the system
around a value of the entropy $S_{0}$ in the grand-canonical ($\mu VT$)
ensemble. As discussed in the first part of this series, the $\mu VT-S$
simulation method provides a direct connexion with classical nucleation theory
McGraw and Laaksonen (1996), since simulations of the nucleation process can
be carried out for any value of the supersaturation $\Delta\mu$. Another
advantage of this method is that it allows the calculation of the entropy of
the system $S$. In the case of binary mixtures, we therefore define the method
as $\mu_{1}\mu_{2}VT-S$ and evaluate the entropy during the simulations
through

$S={{U-N_{1}\mu_{1}-N_{2}\mu_{2}}\over{T(N_{1}+N_{2})}}$ (1)

where $\mu_{1}$ and $\mu_{2}$ are the chemical potentials for the binary
mixture, $N_{1}$ and $N_{2}$ are the number of atoms/molecules for each of the
two components and $U$ is the internal energy for the entire system given by,
in the case of a binary mixture of atoms,

$U=U_{pot}+{3\over 2}N_{1}k_{B}T+{3\over 2}N_{2}k_{B}T$ (2)

where $U_{pot}$ is the potential energy for the system. In the case of linear
molecules, we add a contribution of $k_{B}T$ per molecule to account for the
rotational degrees of freedom.

The $\mu_{1}\mu_{2}VT-S$ simulation relies on gradually decreasing the entropy
of the system through the application of a bias potential. This bias potential
is defined within the framework of the umbrella sampling technique as

$U_{bias}={1\over 2}k(S-S_{0})^{2}$ (3)

in which $S_{0}$ is the target value for the entropy, $S$ is the current value
of the entropy of the system and $k$ is a spring constant. This bias potential
is then added to the potential energy of the system, and the total potential
energy is used in the conventional Metropolis criteria for the acceptance of
the different types of Monte-Carlo ($MC$) steps. For the $Ar-Kr$ mixture, $MC$
steps include the insertion ($12.5$ % of the attempted $MC$ steps) or deletion
($12.5$ % of the attempted $MC$ steps) of atoms as well as the translation of
a single atom ($75$ % of the attempted $MC$ steps). In the case of the
$C_{2}H_{6}-CO_{2}$ mixture, we have the following rates: rotation ($37.5$ %
of the attempted $MC$ steps), translation ($37.5$ %), insertion ($12.5$ %) and
deletion ($12.5$ %).

Successive umbrella simulations with decreasing values for the target entropy
($S_{0}$) are carried out to achieve the formation of a liquid droplet of a
critical size. During each of these simulations, histograms for the number of
times a given entropy interval is visited are collected, allowing for the
calculation of the free energy profile associated with the nucleation process
Torrie and Valleau (1977); Allen and Tildesley (1987); Desgranges and
Delhommelle (2009, 2007); Desgranges and Delhommelle (2014a).

###  II.2 Simulation models

The simulation models used in this work for $Ar$ and $Kr$ are based on the
Lennard-Jones potential Vrabec et al. (2001), with the following parameters
for $Ar$: $\sigma_{Ar}$=$3.3952$ Åand $\epsilon_{Ar}/k_{B}$=$116.79~{}K$ and
for $Kr$: $\sigma_{Kr}$=$3.6274$ Åand $\epsilon_{Kr}/k_{B}$=$162.58$ K.

For the $CO_{2}-C_{2}H_{6}$, we use a force field that models the dispersion-
repulsion interactions through an exp-6 functional form Potoff et al. (1999);
Errington and Panagiotopoulos (1999a, b); Errington (2003).

$\begin{array}[]{lll}u(r)&={\epsilon\over{1-6/\alpha}}\left[{6\over\alpha}exp\left(\alpha\left[1-{r\over
r_{m}}\right]\right)-\left({r_{m}\over r}\right)^{6}\right]&(r>r_{max})\\\
&=\infty&(r<r_{max})\\\ \end{array}$ (4)

In Eq. LABEL:exp6, $r_{m}$ is the distance for which the potential reaches a
minimum, $r_{max}$ is the smallest positive distance for which $du(r)/dr=0$
and $\epsilon$ and $\alpha$ are two potential parameters. As discussed by
Errington et al. Errington and Panagiotopoulos (1999a), it is convenient to
discuss the potential parameters in terms of $\sigma$ (i.e. the distance for
which $u(r)=0$, obtained numerically by solving the equation $u(\sigma)=0$)
rather than in terms of $r_{m}$. In the case of $C_{2}H_{6}$, we use an united
atom-type force field and model the molecule with two $\exp-6$ sites, each
site standing for a $CH_{3}$ group. We use the following set of parameters:
$\epsilon_{CH_{3}}/k_{B}=129.64$ K, $\sigma_{CH_{3}}=3.679$ Å and
$\alpha_{CH_{3}}=16$. In the case of $CO_{2}$, in addition to the the
repulsion-dispersion interactions, a Coloumbic term is added to account for
the quadrupolar nature of $CO_{2}$. The $CO_{2}$ molecule is modeled with a
distribution of three $\exp-6$ sites and three point charges located on each
of the atoms. We use the following parameters for the $\exp-6$ sites:
$\epsilon_{C}/k_{B}=29.07$ K, $\sigma_{C}=2.753$ Å and $\alpha_{C}=14$,
$\epsilon_{O}/k_{B}=83.20$ K, $\sigma_{O}=3.029$ Å and $\alpha_{O}=14$. We
also have $q_{C}=0.6466$ e and $q_{O}=-03233$ e. Both molecules are considered
to be rigid, with a distance between the two $CH_{3}$ $\exp-6$ sites set to
$1.839$ Å for ethane and a length of the $C-O$ bond fixed to $1.1433$ Å for
$CO_{2}$. The calculation of the interaction energy is performed for distances
up to $13.5$ Å for both systems with, electrostatic interactions calculated
with the Ewald sum method beyond that distance using the parameters given in
previous work Desgranges and Delhommelle (2014b). In line with previous
simulation work on nucleation Ten Wolde and Frenkel (1998), we do not include
any tail corrections for the dispersion-repulsion interactions beyond the
cutoff distance.

###  II.3 Setting up the simulation

####  II.3.1 $Ar-Kr$ mixture

We carry out $\mu_{1}\mu_{2}VT-S$ simulations for the $Ar-Kr$ mixture at
$T=148.15$ K in cubic cells with an edge of $100$ Å (the usual periodic
boundary conditions are applied). We take advantage of the Expanded Wang-
Landau (EWL) method we recently developed Desgranges and Delhommelle (2012a,
b, 2014b, 2016b) to obtain very accurate estimates for $\mu_{1}$ and $\mu_{2}$
both at the vapor-liquid coexistence and for supersaturated vapor phases,
which will be the parent phase for the nucleation events. The values are
presented in Table 1

Table 1: $Ar-Kr$ mixture at $148.15$ K: chemical potentials for $Ar$ $(\mu_{1})$ and $Kr$ $(\mu_{2})$, supersaturations, mole fractions in $Kr$ for the vapor ($y_{Kr}$) and for the liquid ($x_{Kr}$), pressure and entropies for the two coexistence points (coex I and coex II) and for the supersaturated vapors considered in this work (system 1 to 4). | $\mu_{1}$ | $\mu_{2}$ | $\Delta\mu_{1}$ | $\Delta\mu_{2}$ | $x_{Kr}$ | $y_{Kr}$ | $P$ | $P/P_{coex}$ | $S_{l}$ | $S_{v}$
---|---|---|---|---|---|---|---|---|---|---
| $(kJ/mol)$ | $(kJ/mol)$ | $(kJ/mol)$ | $(kJ/mol)$ |  |  | $bar$ |  | $(kJ/mol/K)$ | $(kJ/mol/K)$
coex I | -14.077 | -17.571 | - | - | 0.440 | 0.190 | 28.68 | 1.0 | 0.0848 | 0.1074
system 1 | -13.951 | -17.463 | 0.049 | 0.108 | 0.440 | - | 57.27 | 2.0 | 0.0831 | -
system 2 | -13.930 | -17.439 | 0.070 | 0.132 | 0.440 | - | 63.09 | 2.2 | 0.0828 | -
coex II | -13.720 | -18.165 | - | - | 0.250 | 0.108 | 36.92 | 1.0 | 0.0845 | 0.1021
system 3 | -13.548 | -18.057 | 0.172 | 0.108 | 0.250 | - | 73.83 | 2.0 | 0.0817 | -
system 4 | -13.514 | -18.039 | 0.206 | 0.126 | 0.250 | - | 81.28 | 2.2 | 0.0813 | -

We start by determining the conditions for the vapor-liquid coexistence at the
coex I point. In practice, this is done using the EWL method by finding
numerically the values for $\mu_{1}$ and $\mu_{2}$ which lead to equal
probabilities for the vapor and the liquid phase (see more details in previous
work Desgranges and Delhommelle (2014b)). The EWL method also allows to obtain
all thermodynamic properties for the mixture including the mole fractions,
pressure as well as the entropies for the two coexisting phases. These
entropies, which are also given in Table 1, provide an idea of the range of
entropies that need to be sampled for the system to undergo the vapor $\to$
liquid transition. From the coexistence point, we can increase the value of
the two chemical potentials $\mu_{1}$ and $\mu_{2}$ or, in other words, create
a supersaturated vapor that will serve as a starting point for the nucleation
process. We list in Table 1 the two sets of supersaturations (system 1 and
system 2) generated from the coexistence point coex I. As can be seen from
Table 1, increasing $\Delta\mu$ brings the supersaturated vapor more deeply
into the liquid domain of the phase diagram, resulting in a larger value for
the pressure and a lower value for the entropy of the liquid. We proceed along
the same lines from the second coexistence point (coex II) to define two
supersaturated vapors (system 3 and system 4), this time with a mole fraction
in the liquid ($x_{Kr}$ set to $0.250$). We finally add that there are other
ways of determining the chemical potential at the vapor-liquid coexistence and
for supersaturated vapors Gazenm$\ddot{\mathrm{u}}$ller and Camp (2007); Liu
et al. (2011); Potoff et al. (1999); Nezbeda and Kolafa (1991); Singh and
Errington (2006); Rai et al. (2007); Rane et al. (2013); Escobedo and de Pablo
(1996); Shi and Maginn (2008); Eslami and Müller-Plathe (2007); Vogt et al.
(2001); Widom (1963); Siepmann and Frenkel (1992).

We plot in Fig. 1 the successive umbrella sampling windows carried out during
the $\mu_{1}\mu_{2}VT-S$ simulations on system 1. Fig. 1 shows the histograms
corresponding to the probability according to which a given entropy interval
is visited during the $\mu_{1}\mu_{2}VT-S$ simulations. Each of the umbrella
sampling windows (labeled with an index $i$) is obtained by imposing through
Eq 3 a different value of the target entropy $S_{0,i}$. During the nucleation
of a liquid droplet, the system goes from a supersaturated vapor with a low
density, and thus of high entropy, to a system containing a liquid droplet of
a critical size, i.e. to a much more dense system of lower entropy. To observe
the formation of the liquid droplet, we therefore carry out successive
umbrella sampling windows for decreasing values of the target entropy
$S_{0,i}$ and obtain the histograms, shown in Fig. 1, that cover the entire
nucleation process. The progress of the system towards the formation of a
liquid droplet can be followed by monitoring the relative location of
$S_{max,i}$, the entropy for which the histogram $p_{i}(S)$ reaches its
maximum, and of $S_{0,i}$, the target entropy for the $i^{th}$ umbrella
sampling window. We start with the window located to the right of Fig. 1,
associated with the largest value of $S_{0,1}$ which corresponds to a very
dilute vapor ($S_{0,1}=0.11$ kJ/mol/K). As shown in Fig. 1, for the first
window (starting from the right), we have $S_{max,1}=0.1097$ kJ/kg/K, which is
less than the target value $S_{0,1}=0.11$ kJ/kg/K. This means that the target
value $S_{0,1}$ is greater than the entropy of the metastable supersaturated
vapor for the choice of $(\mu_{1},\mu_{2})$ made for system 1. Gradually
decreasing the target entropy for the next windows allows us to find the value
$S_{0,i}$ which coincides with the maximum for $p_{i}(S)$. This occurs here
for $S_{max,i}=S_{0,i}=0.1091$ kJ/kg/K. At this point, we have the metastable
supersaturated vapor. Then, during the next few umbrella sampling windows, we
observe a change in behavior as the histograms $p_{i}(S)$ now lag behind the
target value for the entropy with $S_{max,i}>S_{0,i}$. This corresponds to the
fact that the system has to overcome the free energy cost in forming the
liquid droplet. Later on, for $S_{0,i}=0.103$ kJ/kg/K, we find again that the
maximum for $p_{i}(S)$ coincides with $S_{0,i}$, indicating that we have
reached the top of the free energy barrier of nucleation and that a liquid
droplet of a critical size has formed. For target values of the entropy
greater than $0.103$ kJ/kg/K, we observe again a change in behavior as the
histograms $p_{i}(S)$ run ahead of the target value for the entropy with
$S_{max,i}<S_{0,i}$, corresponding to the spontaneous growth of the liquid
droplet.

Figure 1: $Ar-Kr$ mixture (system 1). Histograms $p(S)$ collected during the
umbrella sampling simulations for decreasing values for the entropy. Figure 2:
$Ar-Kr$ mixture for a supersaturation of $2$ at $x_{Ar}=0.56$ and at
$T=148.15$ K. (Top) Potential energy of the system during a simulation at
entropy fixed to $S_{0}=0.107$ kJ/mol/K (black) and $S_{0}=0.104$ kJ/mol/K
(red). (Bottom) Number of atoms during the simulation at $S_{0}=0.107$
kJ/mol/K ($N_{Ar}$: light green and $N_{Kr}$: dark green) and at $S_{0}=0.104$
kJ/mol/K ($N_{Ar}$: cyan and $N_{Kr}$: blue).

To provide additional insight in the evolution of the average properties of
the system for different umbrella sampling windows, we show in Fig 2 the
potential energy and the number of atoms for different values of the target
entropy $S_{0,i}=0.107$ kJ/mol/K and $S_{0,i}=0.104$ kJ/mol/K during a
production run. These plots show the impact of decreasing the target value for
the entropy on the system. For example, comparing the results for two
different umbrella sampling windows, $0.107$ kJ/mol/K and $0.104$ kJ/mol/K, we
observe a decrease in the potential energy by about $60$ %. This happens
simultaneously with an increase in the number of atoms of each component in
the mixture by about $15$ % for $Ar$ and by $20$ % for $Kr$. This confirms
that the decrease in the target value for the entropy occurs with a greater
organization and density of the system, that we will analyze in depth in the
’Results’ section. We finally add that, for each umbrella sampling window, we
run an equilibration run of $100\times 10^{6}$ MC steps, followed by a
production run of $200\times 10^{6}$ MC steps. Throughout the
$\mu_{1}\mu_{2}VT-S$ simulations, we also check that the acceptance rate for
the insertion/deletion steps remains high enough to ensure an accurate
sampling of the configurations of the systems. For instance, in the case of
system 1 and for the configurations of highest density (umbrella sampling
window with a critical liquid droplet for $S_{0,i}=0.103$ kJ/mol/K), the
acceptance rates for the insertion/deletion steps are of $45.8$ %.

####  II.3.2 $C_{2}H_{6}-CO_{2}$ mixture

Simulations of droplet nucleation for the $C_{2}H_{6}-CO_{2}$ mixture are
performed at $T=263.15$ K in cubic cells, with an edge of $100$ Å, and with
the usual periodic boundary conditions. We use EWL simulations to determine
the other input parameters for the $\mu_{1}\mu_{2}VT-S$ simulations and list
in Table 2 the sets of chemical potentials $(\mu_{1},\mu_{2})$ used for these
simulations.

Table 2: $C_{2}H_{6}-CO_{2}$ mixture at $263.15$ K: chemical potentials for $C_{2}H_{6}$ $(\mu_{1})$ and $CO_{2}$ $(\mu_{2})$, supersaturations, mole fractions in $CO_{2}$ for the vapor ($y_{CO_{2}}$) and for the liquid ($x_{CO_{2}}$), pressure and entropies for the two coexistence points (Coex III and Coex IV) and for the supersaturated vapors considered in this work (system 5 to 8). | $\mu_{1}$ | $\mu_{2}$ | $\Delta\mu_{1}$ | $\Delta\mu_{2}$ | $x_{CO_{2}}$ | $y_{CO_{2}}$ | $P$ | $P/P_{coex}$ | $S_{l}$ | $S_{v}$
---|---|---|---|---|---|---|---|---|---|---
| $(kJ/mol)$ | $(kJ/mol)$ | $(kJ/mol)$ | $(kJ/mol)$ |  |  | $bar$ |  | $(kJ/mol/K)$ | $(kJ/mol/K)$
coex III | -38.186 | -44.397 | - | - | 0.053 | 0.124 | 27.37 | 1.0 | 0.1397 | 0.1556
system 5 | -38.055 | -44.268 | 0.131 | 0.129 | 0.053 | - | 43.80 | 1.6 | 0.1318 | -
system 6 | -38.036 | -44.250 | 0.150 | 0.147 | 0.053 | - | 46.53 | 1.7 | 0.1316 | -
coex IV | -38.261 | -43.216 | - | - | 0.097 | 0.201 | 29.60 | 1.0 | 0.1409 | 0.1566
system 7 | -38.121 | -43.076 | 0.140 | 0.140 | 0.097 | - | 47.35 | 1.6 | 0.1333 | -
system 8 | -38.099 | -43.061 | 0.162 | 0.155 | 0.097 | - | 50.32 | 1.7 | 0.1331 | -

We identify a first state point, coex III, leading to vapor-liquid coexistence
for the $C_{2}H_{6}-CO_{2}$ mixture. As previously discussed, this is achieved
by finding numerically $\mu_{1}$ and $\mu_{2}$ such that the liquid and the
vapor phases are equally probable. The EWL simulations also yield the mole
fractions as well as the entropies of the two coexisting phases, bracketing
the range of entropies needed to sample the nucleation of the liquid droplet.
From coex III, we increase the chemical potentials of the two mixture
components by $\Delta\mu_{1}$ and $\Delta\mu_{2}$ to obtain thermodynamic
conditions located in the domain of the liquid in the phase diagram and, for
which, we can observe a supersatured vapor. We repeat this step for two
different supersaturations leading to system 5 and 6 in Table 2 (as for $Ar-
Kr$, the choices for $\Delta\mu_{1}$ and $\Delta\mu_{2}$ are made such that
the liquid mole fraction in the second component, here $CO_{2}$, remains
constant). We choose another coexistence point, coex IV, and two corresponding
supersaturated vapors, system 7 and system 8. For each umbrella sampling
window, an equilibration run of $50\times 10^{6}$ MC steps is run, followed by
a production run of $100\times 10^{6}$ MC steps. As with the $Ar-Kr$ mixture,
we check that the acceptance rates for the insertion/deletion steps remains
high enough to ensure an accurate sampling. This is the case here for all
umbrella sampling windows. For instance, in the case of system 5 and for the
window associated with the highest density (umbrella sampling window with a
critical liquid droplet for $S_{0,i}=0.155$ kJ/mol/K), the acceptance rates
for the insertion/deletion steps are of $36.9$ % for $C_{2}H_{6}$ and of
$46.6$ % for $CO_{2}$.

##  III Results and Discussion

###  III.1 $Ar-Kr$ mixture

We start by analyzing the free energy barriers obtained for the $Ar-Kr$
mixture at $T=148.15$ K. The left panel of Fig. 3 shows the free energy
profiles for two different supersaturations, system 1 and system 2, at
$x_{Kr}=0.44$. The supersaturation has a direct effect on the height of the
free energy barrier, with a barrier of $25\pm 2~{}k_{B}T$ for system 1 and of
$15\pm 2~{}k_{B}T$ for a system 2. $\Delta\mu_{1}$ and $\Delta\mu_{2}$ are
greater for system 2, which means that the parent supersaturated vapor of
system 2 is located more deeply into the domain of the liquid in the phase
diagram, and, as such, that the nucleation of a liquid droplet occurs more
easily. This also leads to a lower free energy of nucleation than for system
1. We find that the range of entropies spanned during the nucleation process
is also impacted by the amount of supersaturation. For instance, looking at
system 2 in Fig. 3, we find that, for the higher supersaturation, the passage
from the supersaturated vapor to a system containing a droplet of a critical
size occurs over entropies between $0.1085$ kJ/mol/K, for which the free
energy reaches a minimum for the entropy of the metastable parent phase, and
$0.1055$ kJ/mol/K, for which the free energy reaches a maximum corresponding
to the formation of a droplet of a critical size. On the other hand, for
system 1 (lower supersaturation), we find that liquid nucleation takes place
over a much broader entropy range, with the entropic pathway ranging from
$0.1091$ kJ/mol/K to $0.103$ kJ/mol/K. This can be attributed to the
combination of two effects. For a higher supersaturation, the density for the
parent supersaturated vapor is larger and, therefore, its entropy is lower.
Furthermore, at high supersaturation, the thermodynamic conditions lie further
inside the liquid domain. Thus, the critical size for the liquid droplet
becomes smaller, resulting in a system at the top of the free energy barrier
that has a higher entropy. Both effects account for the narrowing, at high
supersaturation, of the entropy range spanned during the droplet nucleation
along the entropic pathway.

We now turn to the second set of simulations carried out for system 3 and
system 4 at $x_{Kr}=0.25$. The free energy profiles obtained for the two
supersaturations are shown in the right panel of Fig. 3. The behavior observed
for $x_{Kr}=0.25$ exhibits similar qualitative features as for $x_{Kr}=0.44$.
We find that the height of the free energy barrier of nucleation increases as
we go from system 4 to system 3 (as the supersaturation decreases), with the
free energy of nucleation increasing from $7\pm 1~{}k_{B}T$ to $21\pm
2~{}k_{B}T$. The range of entropies spanned during nucleation is also found to
increase as the supersaturation decreases. For system 4, we observe that the
entropy of the parent phase is of $0.1047$ kJ/mol/K while the entropy at the
top of the free energy barrier is of $0.103$ kJ/mol/K. For a lower
supersaturation (system 3), the entropy of the supersaturated vapor is
$0.1055$ kJ/mol/K while it is of $0.100$ kJ/mol/K at the top of the free
energy barrier. This means that a significantly broader range of entropies is
sampled during the nucleation process. There are, however, notable differences
between the two plots, which result from the interplay between pressure and
the chemical composition of the system. Considering the results at fixed
supersaturation, we find that the free energy of nucleation decreases by $16$
% (going from system 1 to system 3), while it decreases by $53$ % between
system 2 and system 4. Similarly, the dependence of the height of the free
energy barrier upon supersaturation is also shown to be impacted, with a
decrease by $40$ % from system 1 to system 2 at $x_{Ar}=0.56$ and a much
larger decrease by $66$ % from system 3 to system 4 at $x_{Ar}=0.75$. Both
findings can be attributed to the fact that liquid mixtures with a lower
fraction of $Kr$ are obtained for higher pressures. This, in turn, implies
that the corresponding parent phases will be supersaturated vapors of larger
densities and, therefore, result in lower free energy barriers of nucleation.

Figure 3: Free energy barriers of nucleation for the $Ar-Kr$ mixture at
$T=148.15$ K. (Left panel) system 1 (black) and system 2 (red) for
$x_{Ar}=0.56$. (Right panel) system 3 (black) and system 4 (red) for
$x_{Ar}=0.75$

Figure 4: $Ar-Kr$ mixture: Snapshots of system 1 during the nucleation process
for a target entropy of $0.1075$ kJ/mol/K (left), $0.105$ kJ/mol/K (middle)
and $0.103$ kJ/mol/K (right). Atoms with a liquid-like environment are shown
with larger spheres, with $Ar$ in red and $Kr$ in blue.

The plot shown in Fig. 3 for the free energy profile of nucleation is a
projection of the multi-dimensional free energy surface, as discussed in prior
work on the nucleation in binary systems Reiss (1950). Here, we do not take
into account any nonisothermal effect during nucleation Wyslouzil and Seinfeld
(1992); Wedekind et al. (2007). To interpret further the results obtained for
the free energy of nucleation as a function of S, we now focus on the
interdependence between the entropy S, the droplet size and its composition.
As we have seen from Fig. 2, the number of atoms steadily increases as the
target entropy $S_{0,i}$ is decreased. To check that the decrease in entropy
undergone by the system leads to an increased organization within the system
and to the formation of a liquid droplet, we show in Fig. 4 snapshots of the
system, obtained during umbrella sampling simulations for decreasing values
for the target entropy. As $S_{0,i}$ decreases, the droplet size increases,
with a larger number of atoms being incorporated to the droplet as the target
entropy decreases from $0.1075$ kJ/mol/K, to $0.105$ kJ/mol/K, and finally to
$0.103$ kJ/mol/K. To assess further this point, we perform a detailed analysis
of the size and composition of the droplet along the entropic pathway. The
atoms belonging to the incipient droplet are identified through a commonly
used geometric criterion Ten Wolde and Frenkel (1998). For this purpose, we
determine the distribution for the number of neighbors within a distance of
$5.4$ Å of a central atom in the vapor and the liquid mixture. We show in Fig.
5 these distributions. Both distributions are sharply peaked around
$N_{nab}=1$ for the vapor and around $N_{nab}=10$ for the liquid. As shown in
Fig. 5, atoms belonging to the vapor always have less than $6$ neighbors
within a distance of $5.4$ Å. This allows us to introduce the following
condition to identify a liquid-like atom, or equivalently an atom belonging to
the developing droplet, as an atom with at least $6$ nearest neighbors within
$5.4$ Å.

Figure 5: $Ar-Kr$ mixture: Distributions for the number of neighbors within a
distance of $5.4$ Å for the liquid (solid line) and for the vapor (dashed
line) at coexistence for $T=148.15$ K. Figure 6: $Ar-Kr$ mixture. (Left panel)
Variation of the total number of atoms in the cluster and of the number of
atoms for $Ar$ and $Kr$ at $x_{Ar}=0.56$ for system 1 (top) and system 2
(bottom). (Right panel) Total number of atoms in the cluster, number of atoms
for $Ar$ and $Kr$ against $S$ at $x_{Ar}=0.75$ for system 3 (top) and system 4
(bottom).

Applying this analysis to the configurations generated during the umbrella
sampling windows, leads to the determination of the evolution of the size of
the droplet during the nucleation process. Furthermore, by keeping track of
the identity (either $Ar$ or $Kr$) of the atoms belonging to the droplet, we
can also shed light on the chemical selectivity during the formation of the
cluster. This means that the interplay between droplet size and composition
during the nucleation process can be directly accessed during the simulations,
and, in turn, shed light on the departure in composition of the critical
droplet from the bulk composition Zeng and Oxtoby (1991b); Braun et al. (2014)
and the possible onset of phase separation in partially miscible mixtures
Talanquer and Oxtoby (1995b); ten Wolde and Frenkel (1998); Napari and
Laaksonen (1999).

We show on the left of Fig. 6, the results obtained for the two
supersaturations of systems 1 and 2. In both cases, throughout the nucleation
process, the overall size of the droplet is shown to gradually increase as the
entropy of the system decreases. The critical size of the droplet is found to
be larger at low supersaturation ($N_{c}=713\pm 60$ for system 1 and
$N_{c}=323\pm 30$ for system 2). It is also reached for a lower value of the
entropy at low supersaturation ($S_{c}=0.103$ kJ/mol/K for system 1) than at
high supersaturation ($S_{c}=0.1055$ kJ/mol/K for system 2), in line with the
results obtained for the free energy barrier in Fig. 3. A closer inspection of
the variation for the number of each type of atoms show that the composition
of the droplet does not remain the same during the entire nucleation process.
For small sizes (high entropy), the droplet is richer in $Ar$ with $55\pm 5$ %
of Ar atoms in the droplet at $S=0.1073$ kJ/mol/K for system 1. There is then
a crossover at $S=0.1067$ kJ/mol/K for which the two types of atoms are
equally present. As the entropy further decreases, the droplet becomes richer
in $Kr$ atoms (with a fraction of $48\pm 4$ % in $Ar$ for a droplet of a
critical size). Looking at the composition of the droplet past the critical
size, we find that the fraction of $Ar$ atoms in the droplet increases again
and reaches $49\pm 4$ % for $S=0.102$ kJ/mol/K. The same mechanism is observed
for system 2, with a crossover point located at $S=0.1068$ kJ/mol/K and a
fraction of $Ar$ atoms of $47\pm 4$ % in the critical droplet. These
fluctuations in the nucleus composition through the nucleation process, as
well as the departure of the composition of the droplet from the composition
of the liquid phase, are consistent with prior work on nucleation, based on a
revision of the classical nucleation theory Wilemski (1987), on classical
density functional theory calculations Zeng and Oxtoby (1991b), on molecular
simulations ten Wolde and Frenkel (1998); Yoo et al. (2001); Braun et al.
(2014), or on approaches based on macroscopic kinetics Alekseechkin (2015).
Our results indicate that the mole fraction in $Ar$ in the nucleus to be less
than for the bulk by $0.08$ for System 1 and of $0.03$ for System 2 with
respect to the bulk composition. These deviations, which are moderate due the
almost ideal nature of the $Ar-Kr$ mixture as noted by Zeng and Oxtoby Zeng
and Oxtoby (1991b), did not give rise to a significant phase separation effect
as reported in simulations of more strongly asymmetric mixtures ten Wolde and
Frenkel (1998); Napari et al. (1999).

Looking now at the evolution of the size of the droplet as a function of
entropy for the other set of conditions ($x_{Ar}=0.75$), we see that the size
of the droplet steadily increases as the entropy of the system decreases. The
droplet reaches a size of $629\pm 90$ atoms for the critical droplet in the
case of system 3, and a size of $273\pm 50$ for system 4. As for the previous
system, the critical droplet has formed once the entropy has reached a
critical value $S_{c}=0.100$ kJ/mol/K for system 3, a value that is notably
lower than its counterpart of $0.103$ for system 4. This reflects the fact
that the range for the entropies spanned during the nucleation event becomes
narrower and narrower as the supersaturation is increased. Unlike for systems
1 and 2, we do not observe any crossover between the mole fractions in $Ar$
and $Kr$ in the droplet, given the very large fraction of $Ar$ in systems 3
and 4. As for systems 1 and 2, we observe, however, fluctuations in the
composition of the droplet with the fraction of $Ar$ atoms in the droplet
varying between $60$ % and $80$ % during nucleation.

###  III.2 $CO_{2}-C_{2}H_{6}$ mixture

Turning to the results obtained for the $CO_{2}-C_{2}H_{6}$ mixture, we plot
in Fig. 7, the free energy barrier of nucleation obtained for systems 5 and 6.
For a liquid fraction of $x_{CO_{2}}=0.053$ (left panel of Fig. 7), we observe
that increasing the supersaturation leads to a decrease in the height of the
free energy barrier of nucleation from $21\pm 2~{}k_{B}T$ (system 5) to $12\pm
1~{}k_{B}T$ (system 6). Similarly, for a mole fraction of $x_{CO_{2}}=0.097$,
a lower supersaturation results in a barrier of $12\pm 1~{}k_{B}T$ (system 7),
while a higher supersaturation yields a free energy of nucleation of $7\pm
1~{}k_{B}T$ (system 8). The range of entropies spanned along the nucleation
pathway is also found to depend strongly on supersaturation, and becomes
broader at lower supersaturations. For instance, for system 5, we obtain an
entropy for the parent (supersaturated vapor) phase of $0.161$ kJ/mol/K and an
entropy of the system of $S_{c}=0.151$ kJ/mol/K, when a droplet of a critical
size has formed. Increasing the supersaturation (system 6) leads to a narrower
entropy range, most notably as a result of a decrease in the entropy
$S_{c}=0.154$ kJ/mol/K for which the critical size of the droplet is reached.
Similar conclusions apply for systems 7 and 8. We find a narrowing of the
entropy range spanned during nucleation at high supersaturation, with the top
of the free energy barrier droplet reached for a higher entropy ($S_{c}=0.157$
kJ/mol/K) than at low supersaturation ($S_{c}=0.155$ kJ/mol/K). The higher
entropy required to form a droplet of a critical size at high supersaturation,
together with the lower free energy barrier of nucleation obtained at high
supersaturation, are most likely due to the smaller size of the critical
droplet at high supersaturation, a point that we discuss in the next section
in our analysis of the relation between size, selectivity and entropy during
the nucleation process.

Figure 7: Free energy barriers of nucleation for the $CO_{2}-C_{2}H_{6}$
mixture at $T=263.15$ K. (Left panel) System 5 (black) and system 6 (red) for
$x_{CO_{2}}=0.053$. (Right panel) System 7 (black) and system 8 (red) for
$x_{CO_{2}}=0.097$

We show in Fig. 8 snapshots of configurations of system 5, obtained throughout
the nucleation process. As for the $Ar-Kr$ system, we start from a metastable
supersaturated vapor and carry out umbrella sampling windows with decreasing
values for the target entropy. During the first few umbrella sampling windows,
the target entropy is high enough, so that the droplet is fairly small (see on
the left of Fig. 8 for $S_{0,i}=0.157$ kJ/kg/K). Then, as $S_{0,i}$ is
decreased further, the droplet starts to become larger and larger, and
eventually reaches its critical size for $S_{0,i}=0.151$ kJ/kg/K. This shows
that decreasing the value for the target entropy not only increases the
density of the system, but also results in an increased level of organization
with the formation of the liquid droplet.

As for the Ar-Kr mixture, Fig. 7 shows a projection of the multi-dimensional
free energy surface in the (entropy, free energy) plane. We pursue our
analysis by characterizing the interdependence between the entropy S, the
droplet size and its composition. We start by determining which molecules in
the system have a liquid-like environment and, as such, belong to the droplet.
We determine, for each molecule, the distributions, shown in Fig. 9, for the
number of neighbors for the vapor and the liquid. Neighboring molecules are
defined as being separated by a distance (between the centers of mass of the 2
neighboring molecules) less than $6.4$ Å. The distributions obtained for the
vapor and liquid reach their maxima for very different numbers of neighbors
(for, on average, a single neighbor in the case of the vapor and for $9$
neighbors in the case of the liquid). This allows us to define a molecule as
having a liquid-like environment. and thus belonging to the droplet, if it has
$6$ or more neighbors within a spherical shell of $6.4$ Å.

Figure 8: $C_{2}H_{6}-CO_{2}$ mixture: Snapshots of system 5 during the
nucleation process for a target entropy of $0.157$ kJ/mol/K (left), $0.154$
kJ/mol/K (middle) and $0.151$ kJ/mol/K (right). Molecules with a liquid-like
environment are shown with larger spheres, with the $CO_{2}$ molecules
identified through their $O$ atoms shown in red. Figure 9: $CO_{2}-C_{2}H_{6}$
mixture: Distributions for the number of neighbors within a distance of $6.4$
Å for the liquid (solid line) and for the vapor (dashed line) at coexistence
for $T=263.15$ K. Figure 10: $C_{2}H_{6}-CO_{2}$ mixture. (Left panel)
Variation of the total number of molecules in the cluster and of the number of
molecules for $C_{2}H_{6}$ and $CO_{2}$ at $x_{CO_{2}}=0.053$ for system 5
(top) and system 6 (bottom). (Right panel) Total number of molecules in the
cluster, number of molecules for $C_{2}H_{6}$ and $CO_{2}$ against $S$ at
$x_{CO_{2}}=0.097$ for system 7 (top) and system 8 (bottom).

We now move on to the analysis of the size of the droplet as a function of the
entropy of the system throughout the nucleation process. Fig. 10 shows that,
for all systems, the total number of molecules within the cluster increases
smoothly as the entropy of the system decreases. Furthermore, we find that the
size of the critical droplet decreases as supersaturation is increased. For a
liquid mole fraction of $x_{CO_{2}}=0.053$, the critical size for the droplet
is of $N_{c}=715\pm 35$ molecules for system 5 and of $N_{c}=422\pm 30$
molecules for system 6. The smaller size of the critical droplet at high
supersaturation (system 6) accounts for the higher value of the entropy for
which the system reaches the top of the free energy barrier. Similarly, when
the liquid mole fraction $x_{CO_{2}}=0.097$, the critical size is of
$N_{c}=494\pm 28$ molecules for a low supersaturation (system 7) and of
$N_{c}=303\pm 21$ molecules for a high supersaturation (system 8). The smaller
critical size and higher critical entropy obtained at the higher
supersaturation are, once again, found to be consistent with the free energy
plot of Fig. 7. Turning to the composition of the droplet, we find that
$C_{2}H_{6}$ remains predominant throughout the nucleation process for all
systems (see Fig. 10). We also find, however, that the composition of the
nucleus depends on its size, as nucleation starts with the formation of a
droplet that has a higher $CO_{2}$ mole fraction than the bulk. It is around
$8$ % for systems 5 and 6 for droplets containing a total of $50-100$
molecules. Similarly, considering the same droplet sizes, it is of about $17$
% for systems 7 and 8. The fraction of $CO_{2}$ then decreases as the size of
the droplet increases. For droplets of a critical size, the fraction of
$CO_{2}$ is of $6$ % for systems 5 and 6, while it is of $13$ % for systems 7
and 8. Despite the small sizes of the critical droplets, which contain only a
few hundred of molecules, the fractions in the critical droplets are
reasonably close to the $CO_{2}$ mole fraction of the liquid, and the
departure from the bulk compositions ($5.3$ % of $CO_{2}$ for systems 5 and 6,
and $9.7$ % of $CO_{2}$ for systems 7 and 8) is small. As for binary mixtures
of atoms, the departures in the droplet composition with respect to that of
the bulk are consistent with prior simulations of droplet nucleation in binary
molecular systems (see e.g. recent simulations of the methane-nonane system
Braun et al. (2014)).

##  IV Conclusion

In this work, we propose a new simulation method to study the nucleation
process in binary mixtures of atomic fluids ($Ar-Kr$) and of molecular fluids
($C_{2}H_{6}-CO_{2}$). The method is based on driving the formation of a
liquid droplet through a series of umbrella sampling simulations where the
bias potential is a function of the entropy $S$ of the system. The resulting
approach is implemented within the grand-canonical ensemble and, since the
entropy serves as the reaction coordinate for the nucleation process, is
called $\mu_{1}\mu_{2}VT-S$. The application of the method to the formation of
liquid droplets in binary mixtures sheds light on the interplay between the
size of the droplet, its composition and the supersaturation at which the
nucleation process occurs. Our findings show that, at low supersaturation, the
range of entropies spanned by the nucleation process becomes broader, as a
result of the combined effect of the larger entropy of the metastable
supersaturated vapor (parent phase) and of the lower entropy associated with
the configurations of the system that contain a liquid droplet of a critical
size. These simulations allow us to characterize the critical droplet in terms
of a critical value reached by the entropy at the top of the free energy
barrier of nucleation. We are also able to obtain the free energy profile
along the entropic pathway underlying the formation of the liquid droplet and
to calculate the free energy of nucleation as a function of the
supersaturation and chemical composition of the system. The analysis of the
composition of the droplet shows that the mole fractions fluctuate throughout
the nucleation process and depart from the composition of the bulk. This
departure is however found to become less and less significant as the size of
the droplet increases and its composition starts to conform more and more to
that predicted by thermodynamics. Finally, while the $\mu_{1}\mu_{2}VT-S$ does
not yield directly the nucleation rate, the method allows to generate and
stabilize configurations of the system close to the top of the free energy
barrier. However, as discussed in previous work Ten Wolde et al. (1999), the
nucleation rate can be obtained by carrying out additional molecular molecular
dynamics simulations, using configurations close to the top of the free energy
barrier as a starting point, and following the Bennett-Chandler scheme Bennett
(1977); Chandler (1978); Carter et al. (1989) to determine the kinetics of the
process. Alternatively, the thereshold method of Yasuoka and Matsumoto can
also be used to determine the nucleation rate Yasuoka and Matsumoto (1998).

Acknowledgements Partial funding for this research was provided by NSF through
CAREER award DMR-1052808.

## References

  * Yasuoka and Matsumoto (1998) K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109, 8451 (1998).
  * Oxtoby (1992) D. W. Oxtoby, J. Phys. Condens. Matter 4, 7627 (1992).
  * Shen and Debenedetti (1999) V. K. Shen and P. G. Debenedetti, J. Chem. Phys. 111, 3581 (1999).
  * Weakliem and Reiss (1993) C. L. Weakliem and H. Reiss, J. Chem. Phys. 99, 5374 (1993).
  * Schenter et al. (1999) G. K. Schenter, S. M. Kathmann, and B. C. Garrett, Phys. Rev. Lett. 82, 3484 (1999).
  * Zeng and Oxtoby (1991a) X. C. Zeng and D. W. Oxtoby, J. Chem. Phys. 94, 4472 (1991a).
  * Yi et al. (2002) P. Yi, D. Poulikakos, J. Walther, and G. Yadigaroglu, Int. J. Heat Mass Tran. 45, 2087 (2002).
  * Kinjo et al. (1999) T. Kinjo, K. Ohguchi, K. Yasuoka, and M. Matsumoto, Comput. Mater. Sci. 14, 138 (1999).
  * Toxvaerd (2001) S. Toxvaerd, J. Chem. Phys. 115, 8913 (2001).
  * Ford (1996) I. Ford, J. Chem. Phys. 105, 8324 (1996).
  * Talanquer and Oxtoby (1995a) V. Talanquer and D. Oxtoby, J. Phys. Chem. 99, 2865 (1995a).
  * Reiss et al. (1990) H. Reiss, A. Tabazadeh, and J. Talbot, J. Chem. Phys. 92, 1266 (1990).
  * Kalikmanov and Van Dongen (1995) V. Kalikmanov and M. Van Dongen, J. Chem. Phys. 103, 4250 (1995).
  * Horsch et al. (2008) M. Horsch, J. Vrabec, and H. Hasse, Phys. Rev. E 78, 011603 (2008).
  * Neimark and Vishnyakov (2005) A. V. Neimark and A. Vishnyakov, J. Chem. Phys. 122, 054707 (2005).
  * Oxtoby and Evans (1988) D. W. Oxtoby and R. Evans, J. Chem. Phys. 89, 7521 (1988).
  * Lutsko (2008) J. F. Lutsko, J. Chem. Phys. 129, 244501 (2008).
  * Wang et al. (2008) Z.-J. Wang, C. Valeriani, and D. Frenkel, J. Phys. Chem. B 113, 3776 (2008).
  * Ten Wolde et al. (1999) P. R. Ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 110, 1591 (1999).
  * Gonzalez et al. (2015) M. A. Gonzalez, J. L. Abascal, C. Valeriani, and F. Bresme, J. Chem. Phys. 142, 154903 (2015).
  * Loeffler et al. (2015) T. D. Loeffler, A. Sepehri, and B. Chen, J. Chem. Theory Comput. 11, 4023 (2015).
  * Sosso et al. (2016) G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chem. Rev. (2016).
  * Xu et al. (2015) W. Xu, Z. Lan, B. Peng, R. Wen, and X. Ma, J. Chem. Phys. 142, 054701 (2015).
  * Keasler and Siepmann (2015) S. J. Keasler and J. I. Siepmann, J. Chem. Phys. 143, 164516 (2015).
  * Wilhelmsen et al. (2015) Ø. Wilhelmsen, T. T. Trinh, S. Kjelstrup, and D. Bedeaux, J. Phys. Chem. C 119, 8160 (2015).
  * van Meel et al. (2015) J. van Meel, Y. Liu, and D. Frenkel, Mol. Phys. 113, 2742 (2015).
  * Hale (1986) B. N. Hale, Phys. Rev. A 33, 4156 (1986).
  * Hale (2005) B. N. Hale, J. Chem. Phys. 122, 204509 (2005).
  * Hale and Thomason (2010) B. N. Hale and M. Thomason, Phys. Rev. Lett. 105, 046101 (2010).
  * Yuhara et al. (2015) D. Yuhara, B. C. Barnes, D. Suh, B. C. Knott, G. T. Beckham, K. Yasuoka, D. T. Wu, and A. K. Sum, Faraday Discuss. 179, 463 (2015).
  * Lauricella et al. (2015) M. Lauricella, S. Meloni, S. Liang, N. J. English, P. G. Kusalik, and G. Ciccotti, J. Chem. Phys. 142, 244503 (2015).
  * Singh and Müller-Plathe (2014) J. K. Singh and F. Müller-Plathe, Appl. Phys. Lett. 104, 021603 (2014).
  * Ni and Dijkstra (2013) R. Ni and M. Dijkstra, Soft Matter 9, 365 (2013).
  * Reinhardt and Doye (2014) A. Reinhardt and J. P. Doye, J. Chem. Phys. 141, 084501 (2014).
  * Ten Wolde and Frenkel (1998) P. R. Ten Wolde and D. Frenkel, J. Chem. Phys. 109, 9901 (1998).
  * Chen et al. (2001) B. Chen, J. I. Siepmann, K. J. Oh, and M. L. Klein, J. Chem. Phys. 115, 10903 (2001).
  * Oh and Zeng (1999) K. Oh and X. C. Zeng, J. Chem. Phys. 110, 4471 (1999).
  * Chen et al. (2002) B. Chen, J. I. Siepmann, K. J. Oh, and M. L. Klein, J. Chem. Phys. 116, 4317 (2002).
  * Zhukhovitskii (1995) D. Zhukhovitskii, J. Chem. Phys. 103, 9401 (1995).
  * Nishi et al. (2015) K. Nishi, S. Inoue, and Y. Matsumura, Chem. Phys. Lett. 634, 194 (2015).
  * Lupi et al. (2016) L. Lupi, B. Peters, and V. Molinero, J. Chem. Phys. 145, 211910 (2016).
  * Santiso and Trout (2015) E. E. Santiso and B. L. Trout, J. Chem. Phys. 143, 174109 (2015).
  * Berryman et al. (2016) J. T. Berryman, M. Anwar, S. Dorosz, and T. Schilling, J. Chem. Phys. 145, 211901 (2016).
  * Zimmermann et al. (2015) N. E. Zimmermann, B. Vorselaars, D. Quigley, and B. Peters, J. Am. Chem. Soc. 137, 13352 (2015).
  * Lam et al. (2015) J. Lam, D. Amans, C. Dujardin, G. Ledoux, and A.-R. Allouche, J. Phys. Chem. A 119, 8944 (2015).
  * Kratzer and Arnold (2015) K. Kratzer and A. Arnold, Soft matter 11, 2174 (2015).
  * Bolhuis and Dellago (2015) P. Bolhuis and C. Dellago, Eur. Phys. J. Special Topics 224, 2409 (2015).
  * Lau et al. (2015) G. V. Lau, P. A. Hunt, E. A. Müller, G. Jackson, and I. J. Ford, J. Chem. Phys. 143, 244709 (2015).
  * Toxvaerd (2016) S. Toxvaerd, J. Chem. Phys. 144, 164502 (2016).
  * Tanaka et al. (2005) K. K. Tanaka, K. Kawamura, H. Tanaka, and K. Nakazawa, J. Chem. Phys. 122, 184514 (2005).
  * Kraska (2006) T. Kraska, J. Chem. Phys. 124, 054507 (2006).
  * Oh and Zeng (2000) K. Oh and X. C. Zeng, J. Chem. Phys. 112, 294 (2000).
  * Senger et al. (1999) B. Senger, P. Schaaf, D. Corti, R. Bowles, D. Pointu, J.-C. Voegel, and H. Reiss, J. Chem. Phys. 110, 6438 (1999).
  * Kulmala and Laaksonen (1990) M. Kulmala and A. Laaksonen, J. Chem. Phys.. 93, 696 (1990).
  * Zeng and Oxtoby (1991b) X. C. Zeng and D. Oxtoby, J. Chem. Phys. 95, 5940 (1991b).
  * Oxtoby and Kashchiev (1994) D. W. Oxtoby and D. Kashchiev, J. Chem. Phys. 100, 7665 (1994).
  * Napari and Laaksonen (1999) I. Napari and A. Laaksonen, J. Chem. Phys. 111, 5485 (1999).
  * Jaecker-Voirol and Mirabel (1988) A. Jaecker-Voirol and P. Mirabel, J. Phys. Chem. 92, 3518 (1988).
  * Talanquer and Oxtoby (1995b) V. Talanquer and D. W. Oxtoby, J. Chem. Phys. 102, 2156 (1995b).
  * ten Wolde and Frenkel (1998) P. R. ten Wolde and D. Frenkel, J. Chem. Phys. 109, 9919 (1998).
  * Laaksonen and Oxtoby (1995) A. Laaksonen and D. W. Oxtoby, J. Chem. Phys. 102, 5803 (1995).
  * Yoo et al. (2001) S. Yoo, K. Oh, and X. C. Zeng, J. Chem. Phys. 115, 8518 (2001).
  * Napari and Laaksonen (2000) I. Napari and A. Laaksonen, Phys. Rev. Lett. 84, 2184 (2000).
  * Braun et al. (2014) S. Braun, V. Kalikmanov, and T. Kraska, J. Chem. Phys. 140, 124305 (2014).
  * Shimizu and Tanaka (2015) R. Shimizu and H. Tanaka, Nature Commun. 6, 7407 (2015).
  * Pinho et al. (2014) B. Pinho, S. Girardon, F. Bazer-Bachi, G. Bergeot, S. Marre, and C. Aymonier, Lab on a Chip 14, 3843 (2014).
  * Gao et al. (2014) X. Gao, D. Fu, B. Xie, Y. Su, and D. Wang, J. Phys. Chem. B 118, 12549 (2014).
  * Alekseechkin (2015) N. V. Alekseechkin, J. Chem. Phys. 143, 054502 (2015).
  * Watson et al. (2011) K. D. Watson, S. T. Nguelo, C. Desgranges, and J. Delhommelle, CrystEngComm 13, 1132 (2011).
  * Desgranges and Delhommelle (2014a) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 136, 8145 (2014a).
  * Desgranges and Delhommelle (2016a) C. Desgranges and J. Delhommelle, J. Chem. Phys. - Part I (accepted for publication) (2016a).
  * McGraw and Laaksonen (1996) R. McGraw and A. Laaksonen, Phys. Rev. Lett. 76, 2754 (1996).
  * Torrie and Valleau (1977) G. M. Torrie and J. P. Valleau, J. Comput, Phys. 23, 187 (1977).
  * Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley, _Computer Simulation of Liquids_ (Clarendon Press, Oxford, 1987).
  * Desgranges and Delhommelle (2009) C. Desgranges and J. Delhommelle, J. Phys. Chem. C 113, 3607 (2009).
  * Desgranges and Delhommelle (2007) C. Desgranges and J. Delhommelle, J. Chem. Phys. 127, 144509 (2007).
  * Vrabec et al. (2001) J. Vrabec, J. Stoll, and H. Hasse, J. Phys. Chem. B 105, 12126 (2001).
  * Potoff et al. (1999) J. J. Potoff, J. R. Errington, and A. Z. Panagiotopoulos, Mol. Phys. 97, 1073 (1999).
  * Errington and Panagiotopoulos (1999a) J. R. Errington and A. Z. Panagiotopoulos, J. Phys. Chem. B 103, 6314 (1999a).
  * Errington and Panagiotopoulos (1999b) J. R. Errington and A. Z. Panagiotopoulos, J. Chem. Phys. 111, 9731 (1999b).
  * Errington (2003) J. R. Errington, J. Chem. Phys. 118, 9915 (2003).
  * Desgranges and Delhommelle (2014b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 140, 104109 (2014b).
  * Desgranges and Delhommelle (2012a) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184107 (2012a).
  * Desgranges and Delhommelle (2012b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184108 (2012b).
  * Desgranges and Delhommelle (2016b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 144, 124510 (2016b).
  * Gazenm$\ddot{\mathrm{u}}$ller and Camp (2007) G. Gazenm$\ddot{\mathrm{u}}$ller and P. J. Camp, J. Chem. Phys. 127, 154504 (2007).
  * Liu et al. (2011) Y. Liu, A. Z. Panagiotopoulos, and P. G. Debenedetti, J. Phys. Chem. B 115, 6629 (2011).
  * Nezbeda and Kolafa (1991) I. Nezbeda and J. Kolafa, Molec. Simul. 5, 391 (1991).
  * Singh and Errington (2006) J. K. Singh and J. R. Errington, J. Phys. Chem. B 110, 1369 (2006).
  * Rai et al. (2007) N. Rai, J. I. Siepmann, N. E. Schultz, and R. B. Ross, J. Phys. Chem. C 111, 15634 (2007).
  * Rane et al. (2013) K. S. Rane, S. Murali, and J. R. Errington, J. Chem. Theory Comput. 9, 2552 (2013).
  * Escobedo and de Pablo (1996) F. Escobedo and J. J. de Pablo, J. Chem. Phys. 105, 4391 (1996).
  * Shi and Maginn (2008) W. Shi and E. J. Maginn, J. Comp. Chem. 29, 2520 (2008).
  * Eslami and Müller-Plathe (2007) H. Eslami and F. Müller-Plathe, J. Comput. Chem. 28, 1763 (2007).
  * Vogt et al. (2001) P. S. Vogt, R. Liapine, B. Kirchner, A. J. Dyson, H. Huber, G. Marcelli, and R. J. Sadus, Phys. Chem. Chem. Phys. 3, 1297 (2001).
  * Widom (1963) B. Widom, J. Chem. Phys. 39, 2808 (1963).
  * Siepmann and Frenkel (1992) J. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992).
  * Reiss (1950) H. Reiss, J. Chem. Phys. 18, 840 (1950).
  * Wyslouzil and Seinfeld (1992) B. Wyslouzil and J. Seinfeld, J. Chem. Phys. 97, 2661 (1992).
  * Wedekind et al. (2007) J. Wedekind, D. Reguera, and R. Strey, J. Chem. Phys. 127, 064501 (2007).
  * Wilemski (1987) G. Wilemski, Journal of Physical Chemistry 91, 2492 (1987).
  * Napari et al. (1999) I. Napari, A. Laaksonen, V. Talanquer, and D. W. Oxtoby, The Journal of chemical physics 110, 5906 (1999).
  * Bennett (1977) C. H. Bennett, _Algorithms for chemical computations_ (ACS symposium Series, vol. 46, American Chemical Society, Washington DC, 1977).
  * Chandler (1978) D. Chandler, J. Chem. Phys. 68, 2959 (1978).
  * Carter et al. (1989) E. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys. Lett. 156, 472 (1989).