Free energy calculations along entropic pathways. I. Homogeneous vapor-liquid nucleation for atomic and molecular systems

In recent years, computer simulations have considerably advanced our knowledge on phase transitions and molecular assembly processes by shedding light on the microscopic mechanisms underlying these phenomena.^1–34 Simulations have allowed to identify key parameters for these activated processes, which allow to measure and analyze the progress of the system towards the formation of a new phase. These parameters, which are often complex functions of, e.g., the coordinates of the atoms or of some features of their collective behavior, are known as order parameters or reaction coordinates (RCs).^35–47 RCs are generally tailored to a type of system or to a type of phenomenon under study. This is the case, for instance, for vapor-liquid nucleation, a process of great interest for many applications such as for atmospheric nucleation. Possible choices of RCs are the radius of the incipient droplet or its size (in terms of numbers of atoms within the droplet),^35,48–53 and the analysis of the underlying pathway is generally made through the determination of the free energy barrier of nucleation and of the critical size of the nucleus beyond which spontaneous growth occurs.

The entropy of the system S is an especially appealing choice for a RC since it is, by definition, directly connected to the onset of order (i.e., S decreases as order increases). It also has the additional advantage of being very general and, as such, can be used for a wide range of systems and phenomena, regardless of the geometry, structure, and phases involved. In the first part of this series, our goal is to establish the use of S as a RC for the homogeneous vapor-liquid nucleation of atomic and molecular systems. For this purpose, we develop a simulation method that allows us to vary the entropy of the system and, as a result, to drive the formation of the liquid droplet. We then calculate the free energy barrier associated with this process. The simulations carried out in this work therefore allow us to identify the entropic pathway underlying the nucleation process, which starts with the metastable supersaturated vapor and ends at the top of the free energy barrier, when a liquid droplet of a critical size has formed. By varying the conditions of nucleation, we are also able to analyze the impact of the supersaturation on the nucleation process and to rationalize the interplay between the free energy of nucleation, the size of the liquid droplet, and the range of entropies spanned during the nucleation process. The validity of the results is assessed through a comparison with the predictions from the classical nucleation theory,^54,55 which provides a relation between the free energy of nucleation and the supersaturation, and with prior simulation work.^35,56 Results are presented for both atomic (Ar) and molecular (N₂ and CO₂) systems to demonstrate the versatility and the transferability of the method to molecular systems.

The paper is organized as follows. In Sec. II, we present the simulation method and discuss how the entropy of the system is used as a RC. We also explain how we calculate the free energy profile for the nucleation process. We then present the molecular force fields used to model the three systems studied in this work, before discussing the technical details. In particular, we give an account of how we obtain and choose the input parameters for the simulations to control the extent of the supersaturation during the nucleation process. Then, we discuss the results obtained here on the three systems considered in this work, Ar, N₂, and CO₂, and focus on quantifying the impact of the supersaturation on the height of the free energy barrier, as well as on the size of the liquid droplet and on the entropy of the system at the top of the free energy barrier. We finally draw the minimum conclusions from this work in Sec. IV.

The $μVT−S$ approach is based on the grand-canonical ensemble, where $μ$ is the chemical potential, V the volume, T the temperature, and S the entropy of the system. Having $μ$ as an input parameter for the simulations leads to a direct comparison with the theoretical predictions from the classical nucleation theory (CNT).^54,55 CNT analyzes the nucleation process, and the free energy barrier that the system has to overcome during this process, as a competition between two contributions, which cancel out when a droplet of a critical size has formed, i.e., when the system has reached the top of the free energy barrier. The first contribution, the surface term, is energetically unfavorable and is the predominant term for small droplets. It results from the cost of creating the interface between the liquid droplet and the vapor and can be evaluated as the product of the area of the incipient spherical droplet by the surface tension. The second contribution, the volume term, is energetically favorable and becomes predominant for larger droplets. The volume term corresponds to the energetic gain in “converting” the parent phase, the metastable supersaturated vapor, into a droplet of the stable phase, the liquid phase. In the CNT approach,^54,55 the volume term, and hence the free energy of nucleation, is proportional to the supersaturation $Δμ=μ−μv$⁠, which is the difference between the chemical potential of the liquid (⁠$μ$⁠) and the chemical potential of the supersaturated vapor $μv$ at the same pressure as the liquid. Carrying out simulations in the grand-canonical ensemble allows us to choose $μ$ and hence the supersaturation $Δμ$ for the nucleation event.

Working in the grand-canonical ensemble also has additional advantages. Since the number of atoms/molecules is allowed to fluctuate, this ensemble provides a direct access to thermodynamic quantities that are difficult to calculate in other ensembles, as, e.g., in the NVT or NPT ensemble. This is the case, for instance, for the entropy S of the system. In the grand-canonical ensemble, in which the chemical potential $μ$⁠, the temperature T, and volume V are fixed during the simulations, it is straightforward to evaluate S during the simulations through

where $U¯$ is the internal energy per atom/molecule, obtained by calculating the potential energy due to the interactions between atoms/molecules and by adding the kinetic (ideal gas) contribution to the internal energy of k_BT/2 per degree of freedom.

This means that, in the grand-canonical ensemble, S can be used as the reaction coordinate to measure the onset of order in the system and, in the case of nucleation, its progress towards the formation of a liquid droplet of a critical size. To achieve this, we carry out $μVT−S$ simulations by taking advantage of the well-known umbrella sampling method⁵⁷ to allow the system to follow the RC, in this case S. From a practical standpoint, we add a bias potential energy, a function of the entropy of the system, U_bias(S), to the total potential energy of the system. In this work, we use the following harmonic function of S as the bias potential energy:

where k is the spring constant and S₀ is the target value for the entropy of the system.

$μVT−S$ simulations are carried out within the Monte Carlo (MC) framework, with the usual Metropolis criteria used to accept/reject the different types of MC moves attempted on the atoms/molecules of the system. More specifically, since these simulations are rooted in the grand-canonical ensemble, the conventional Metropolis criteria, as obtained, e.g., by Allen and Tildesley⁵⁸ in the grand-canonical ensemble, are used in $μVT−S$ simulations. For instance, for the MC steps involving the translation of a single atom/molecule, as well as for the rotation of a molecule, we use the following Metropolis criterion:

with $ΔU$ denoting the change in the total potential energy (interaction energy + bias potential) corresponding to the move from the old (o) configuration to the new (n) configuration. Similarly, for the MC moves involving the insertion of a new atom/molecule, we use the following criterion:

where z is the activity $z=exp(βμ)/Λ3$⁠, where Λ is the de Broglie wavelength and N is the current number of atoms/molecules in the system. In the case of molecular fluids, one also needs to take into account the terms related to the other degrees of freedom. For linear molecules like N₂ and CO₂, we include for the rotation a factor of $8πIkBT/2h2$⁠, where I is the moment of inertia of the linear molecule, k_B is the Boltzmann constant, and h is the Planck’s constant.

For all systems (Ar, CO₂, and N₂), we carry out $μVT−S$ simulations in cubic cells with an edge of 100 Å  and apply the usual 3D periodic boundary conditions. To simulate the entire nucleation process, we carry out a series of umbrella sampling simulations with decreasing values for the target entropy S₀. This allows us to sample the underlying entropic pathway, which goes from the high entropy parent phase (supersaturated vapor) to the low entropy phase (liquid). During the course of the simulation, we collect histograms for the number of times each entropy interval is visited. These histograms are then used to build the free energy profile of nucleation using the techniques developed to analyze the results from umbrella sampling simulations (more details can be found in previous work^57–61). For each value of the target entropy S₀, simulations are first run for $100×106$ MC steps to allow the system to relax. Then, a second run, the production run, is performed for $200×106$ MC steps during which the data are collected and the calculation of average properties is carried out. The different types of MC steps are attempted with the following rates. In the case of Ar, 75 of the MC steps consist in the translation of a single atom, while the 25% remaining moves are equally split between the MC steps corresponding to the insertion and deletion of Ar atoms. In terms of computational efficiency, working in the grand-canonical ensemble allows the number of atoms to vary, which means that the simulations for the higher S values will be performed on systems with few atoms. This is advantageous when compared to methods based on the NPT ensemble, which simulate systems with a constant number of atoms. Furthermore, calculating the reaction coordinate S through Eq. (1) after every time step is very fast. The CPU time for a $μVT−S$ run is therefore very close to a conventional grand-canonical run. For instance, a production run carried out on systems containing an Ar droplet of a critical size (i.e., the largest system sizes during nucleation) takes 30 CPU hours (CPU time given for a single Intel Xeon processor clocked at 3.2 GHz). For CO₂ and N₂, we take into account the MC steps for the rotation of a single molecule (37.5%), in addition to the translation (37.5 %) and the insertion/deletion (25 %). Throughout the simulations, we check that we have sufficiently large acceptance rates for the insertion/deletion steps to ensure an efficient sampling. For all systems and all target entropies S₀, we obtain acceptance rates greater than 40 % for the insertion/deletion steps.

We model argon with a Lennard-Jones (LJ) potential

where $ϵ= 117.05K$ is the depth of the potential well, $σ=3.4$ Å  is the atom diameter, and r is the distance between two interacting atoms. The interactions between Ar atoms are calculated for all distances up to 13.6 Å  and neglected beyond that cutoff distance. We use this large cutoff distance (⁠$4σ$⁠) and do not apply tail corrections, as in previous work on the vapor-liquid nucleation.³⁵ We add that an alternative approach, relying on the calculation of long-range corrections as a function of the local density, has been developed in recent years for inhomogeneous systems.⁶² 

For N₂ and CO₂, we consider both molecules to be rigid and model the interactions with a distribution of point charges to account for the molecular quadrupole and with a distribution of LJ sites to model the dispersion-repulsion interactions. For N₂, we consider two LJ sites per N₂ molecule (one on each atom) and three point charges (one negative charge on each atom and one positive charge at the center of the N–N bond) to model the quadrupole of the N₂ molecule. We use the following parameters:⁶³ for the LJ sites, $ϵNN/kB=36$ K and $σNN=3.30$ Å , and for the point charges, $qcenter=0.966e$ and q_N = −q_center/2. The N–N bond length is set to 0.549 Å  and uses a spherical cutoff set to $4σNN=13.2$ Å  for the LJ part of the potential and to 15 Å  for the quadrupole-quadrupole interactions.

For CO₂, we use the TraPPE force field⁶⁴ which is based on a distribution of three LJ sites (located on the atoms of the molecule) and of three atomic charges, with the following parameters: $ϵCC/kB=27K$⁠, $ϵOO/kB=79K$⁠, $σCC=2.80$ Å  and $σOO=3.05$ Å , and $qc=0.7e$ and $qO=−0.35e$⁠. The parameters for the LJ interactions between unlike atoms are given by the Lorentz-Berthelot mixing rules. The bond length is set to 1.16 Å  as in the original paper.⁶⁴ We use a spherical cutoff set to $4σOO=12.2$ Å  for the LJ part of the potential and to 14.6 Å  for the quadrupole-quadrupole interactions.

$μVT−S$ simulations for Ar are carried out at T = 128.76 K. We present in Table I the conditions for the three supersaturations studied in this work. They are obtained through the Expanded Wang-Landau (EWL) simulation method we recently developed^65–68 (other methods can also be used to determine this information^64,69–80). EWL simulations provide the grand-canonical partition function and the probability distribution for the number of Ar atoms in the system for any value of $μ$⁠. From there, the value of the chemical potential at the vapor-liquid coexistence can be directly obtained by finding the value of $μ$ leading to equal probabilities for the vapor and liquid phases (more details are provided in previous work⁶⁵). The EWL simulations also give access to all other thermodynamic properties, including the pressure P and entropies at coexistence, S_l and S_v, also given in Table I. This provides an estimate for the range of entropies that needs to be sampled during the $μVT−S$ simulations to cover the entire vapor → liquid transition. We also give in Table I the supersaturation $Δμ$ that we apply to obtain the supersaturated vapors, which serve as parent phases for the nucleation processes. Increasing the value of the chemical potential by $Δμ$ (⁠$Δμ>0$⁠) pushes the system further into the domain of the phase diagram where the liquid phase is stable. We add that $μv$ that appears in the definition of the supersaturation $Δμ$ is very close to $μcoex$ for the low pressures studied here (of the order of 10⁻³ kJ/kg). Here, we consider 3 different supersaturations and label each set of conditions as system 1, 2, and 3.

The $μVT−S$ simulations are initiated from a very low density vapor, obtained, e.g., by placing an Ar atom in the simulation cell. We then carry out a series of umbrella sampling simulations for decreasing values of the target entropy S_0,i, where i is the index labeling the ith umbrella sampling simulation. A reasonable estimate for the starting value for the target entropy is provided by S_v, the value of the vapor at coexistence. We show in Fig. 1 the histograms p_i(S) for the number of times a given entropy interval is visited during the $μVT−S$ simulations for system 2. For each umbrella sampling simulation, comparing the relative values of S_max,i, the entropy for which the histogram p_i(S) reaches its maximum, to the target entropy S_0,i, imposed during the ith umbrella sampling simulation, provides a great deal of insight into the nucleation process. For instance, at the beginning of the nucleation process, the system tries to overcome the free energy barrier of nucleation. This means that increasing the size of the liquid droplet (achieved here by decreasing the entropy of the system) is associated with the large cost, in terms of free energy, due to the predominant surface term. As a result, the histogram p_i(S) for the entropy of the system lags behind the target value for the entropy and that $Smax,i>S0,i$⁠. On the other hand, when the system is past the top of the free energy barrier, the droplet has a larger probability to keep growing and we have the reverse situation with $Smax,i<S0,i$⁠. The same applies to very low density vapors, when the target value for the entropy is greater than the entropy of the metastable supersaturated vapor. Therefore, the entropies S_max,i and S_0,i only coincide when the system reaches a point corresponding to an extremum of the free energy profile. This extremum can be either a maximum (i.e., the top of the free energy barrier, when a liquid droplet of a critical size has formed) or a minimum (i.e., the metastable supersaturated vapor or the stable liquid).

Looking at the first umbrella sampling window, starting from the right of Fig. 1, we find that S_max,1 (2.494 (kJ/kg)/K) is lower than S_0,1 (2.5 (kJ/kg)/K). This shows that, at this stage, the entropy of the system is greater than the entropy of the metastable supersaturated vapor, or, in other words, that the density of the system has not yet reached the density of the supersaturated vapor. Decreasing the target entropy allows us to identify the first value of the entropy for which S_max,i = S_0,i, where we have the metastable supersaturated vapor for an entropy of 2.49 (kJ/kg)/K. Then, during the next umbrella sampling simulations for decreasing S_0,i, the system climbs up the free energy barrier of nucleation and we have $Smax,i>S0,i$⁠, until we observe again the coincidence of S_max,i with S_0,i for an entropy of 2.3 (kJ/kg)/K. This occurs when the system reaches the top of the free energy barrier and a liquid droplet of a critical size has formed. These results show that the $μVT−S$ method allows us to explore the entropic pathway underlying the nucleation process and identify the extrema of the free energy profile. We also show in Fig. 2 the impact of changing the target entropy on the properties of the system. Comparing the results for target entropies of 2.43 (kJ/kg)/K and 2.38 (kJ/kg)/K shows a decrease in the potential energy of the system by close to 30% and an increase in the number of Ar atoms in the system by about 10%. Both results can be accounted for by the increase in the density of the system, triggered by the decrease in the target entropy, which leads to an increased number of attractive interactions between Ar atoms. This shows that the decrease in the entropy achieved through the $μVT−S$ simulations results in a transition towards the liquid phase. We leave the detailed discussion of the energetics and characteristics of the liquid droplet to Sec. III.

The formation of a liquid droplet of N₂ from a supersaturated vapor is studied at $T=100K$⁠. As for Ar, we use the EWL simulation method to obtain the properties at coexistence (coex) and at three supersaturations (systems 4, 5, and 6), which are increasingly deeper inside the liquid domain of the phase diagram of N₂. We list in Table II the data for these points, from which we simulate the nucleation process.

In line with the two previous systems, we start from the conditions for vapor-liquid coexistence (coex) provided by the EWL simulations. We then define three supersaturations, for increasing $Δμ$⁠, leading to systems 7, 8, and 9. The thermodynamic data at coexistence and for the three supersaturations are given in Table III.

We start by discussing the results obtained for the free energy profiles of nucleation for argon at T = 128.76 K. We show in Fig. 3 the free energy barriers obtained from the $μVT−S$ simulations for systems 1, 2, and 3 against the entropy of the system, which serves as the reaction coordinate for the nucleation process. Two main trends appear in this plot. First, the height of the free energy barrier increases as the supersaturation $Δμ$ decreases. When $Δμ=2.84$ kJ/kg (system 3), we obtain a free energy barrier of nucleation of 26$±$2 $kBT$⁠, while for $Δμ=2.47$ kJ/kg (system 2), we have a free energy barrier of 41$±$3 $kBT$⁠. Moreover, as the supersaturation further decreases to $Δμ=2.04$ kJ/kg (system 1), the free energy barrier becomes $62±5kBT$⁠. This behavior is consistent with the predictions from the classical nucleation theory and with the findings from previous simulation work. Changing the supersaturation $Δμ$ has a direct impact on the volume term, which is the gain, in free energy, resulting from the formation of the thermodynamically stable liquid phase from the metastable vapor phase. This volume term is proportional to $Δμ$ and it therefore follows that for a small $Δμ$ (system 1), the top of the free energy barrier is reached later, for a larger liquid droplet, leading to a higher free energy barrier. The values obtained in this work for the free energy barriers of nucleation are also in good agreement with those predicted by the classical nucleation theory^54,55 and with those found for the Lennard-Jones system.³⁵ For instance, CNT leads to the following free energy of nucleation:^54,55 $W= 16πγ3/(3ρl2Δμ2)$⁠, where $γ$ is the surface tension for a flat interface and $ρl$ is the density of the (bulk) liquid. Applying this to, e.g., system 1 with $γ=3.535×10−10kJ/cm2$ (determined from the EWL simulations⁶⁵) and $ρl=1.0882g/cm3$ leads to a free energy barrier of $70kBT$⁠, in reasonably good agreement with the barrier found here of $62±5kBT$⁠. Our results are also consistent with ten Wolde and Frenkel³⁵ who reported free energy barriers of nucleation $5–15kBT$ lower than the CNT predictions for the Lennard-Jones system. The second main trend is the following. The entropic pathway spanned during nucleation covers a range of entropy that depends on the supersaturation. More specifically, for a large supersaturation, the entropy range becomes narrower. This is due to the fact that for a low supersaturation (i.e., not too deep into the liquid domain of the phase diagram), the metastable supersaturated vapor is less dense. This means that the starting point for the nucleation process is associated with a larger entropy. Furthermore, at low supersaturation, the size of the critical nucleus is larger. This implies that, when the supersaturation is low, the entropy of the system at the top of the free energy barrier is smaller. Both factors account for the extended range of entropy spanned during the nucleation process at low supersaturation. We finally add that the range of entropy spanned during nucleation is, for all supersaturations, well within the interval defined by the entropies of the vapor (2.701 (kJ/kg)/K) and liquid (1.856 (kJ/kg)/K) at coexistence. This is the expected behavior, since the metastable supersaturated vapors are more dense, with a smaller entropy, than the vapor at coexistence, and since systems containing a critical liquid droplet are considerably less dense, with a larger entropy, than a uniform liquid phase like the liquid at coexistence.

The plots in Fig. 3 also allow us to characterize the state of the system at the top of the free energy barrier of nucleation in terms of the critical entropy associated with the formation of a critical nucleus. This critical entropy S_c is shown to depend on the supersaturation, with the value for S_c being obtained for the lowest supersaturation $Sc=2.21±0.02$ kJ/kg (system 1) and the value $Sc=2.35±0.01$ kJ/kg for the highest supersaturation (system 3). To analyze further the characteristics of the liquid droplet, we examine the relation between the entropy and the size of the liquid droplet. Starting from the approach developed by ten Wolde and Frenkel,³⁵ we first determine the distribution for the number of neighbors in the vapor and in the liquid. Here, neighboring particles are defined as particles that are less than $1.6σ$ apart. This distance of $1.6σ$ corresponds to the first minimum of the radial pair distribution function of the liquid. Fig. 4(a) shows the distribution so obtained for the vapor and for the liquid phases. The two distributions show little overlap, with the vast majority of the liquid atoms having at least 6 neighbors, while atoms in the vapor phase very rarely have more than 5 neighbors. This allows us to define liquid-like atoms as having at least 6 neighbors and vapor-like atoms having fewer than 6 neighbors and to follow the variations of the number of liquid-like atoms during nucleation. Applying this analysis during the $μVT−S$ simulations leads us to draw a correspondence between the value of the entropy and the number of liquid-like atoms of the system as nucleation proceeds. We present in Fig. 4(b) the result for this correspondence for system 2. The number of liquid-like atoms increases steadily as the nucleation process advances, implying that the size of the liquid droplet steadily increases as the entropy of the system decreases. This can best be seen by examining snapshots of the configurations of the system during the nucleation process. Fig. 5 shows the formation and development of the liquid droplet along the entropic pathway. As can be seen on the snapshots, the droplet size steadily increases with the entropy and reaches a critical size of $1715±50$ atoms for system 1. The snapshots for the other 2 supersaturations reveal a similar behavior, and critical sizes of $980±25$ and $649±15$ atoms are found for the two other supersaturation systems 2 and 3, respectively. We obtain the expected correlation between the height of the free energy barrier of nucleation and the size of the critical nucleus, as the size of the critical droplet decreases with the free energy of nucleation as we go from system 1 to system 2 and finally to system 3.

We now turn to the assessment of the versatility of the method and test it on a molecular fluid. We apply the $μVT−S$ approach to study the nucleation of a liquid droplet of N₂ at 100 K. We follow the same procedure as for Ar and perform a series of umbrella sampling simulations spanning the entropic pathway underlying the formation of the droplet. We present in Fig. 6 the free energy profiles of nucleation of N₂ at T = 100 K for three supersaturations (systems 4, 5, and 6). The results show features that are qualitatively similar to those found for Ar. In terms of the entropy range, the entropies sampled during the process lie within the wider interval defined by the entropy of the liquid and of the vapor at coexistence (from S_liq = 3.374 (kJ/kg)/K to S_vap = 4.912 (kJ/kg)/K). For instance, for the lowest supersaturation (system 4), the nucleation starts from a supersaturated vapor with S = 4.6 (kJ/kg)/K (lower than the entropy of the coexisting vapor that has a lower density) and ends with configurations containing a critical nucleus at S = 4.08 (kJ/kg)/K (higher than the entropy of the liquid at coexistence, which has a much higher density). The same reasoning applies to systems 5 and 6, for which the nucleation process is associated with an even narrower entropy range due to the higher supersaturations involved for these systems. In particular, the entropy S_c for which a liquid droplet of a critical size has formed is shown to increase with supersaturation, starting from $Sc=4.08±0.02$ (kJ/kg)/K (system 4) and increasing to $Sc=4.16±0.02$ (kJ/kg)/K (system 5) and finally to $Sc=4.20±0.02$ (kJ/kg)/K (system 6) for the highest supersaturation. The free energy barriers of nucleation for N₂ are also shown to decrease gradually as the supersaturation is increased. The free energy barrier of nucleation reaches $73±5kBT$ for system 4 (lowest supersaturation), $60± 4kBT$ for system 5, and $52± 4kBT$ for system 6 (highest supersaturation).

We further refine our analysis by determining the size of the droplet during the nucleation process. For this purpose, we adopt a similar procedure to that followed for argon. We compute the distributions for the number of neighbors for N₂ molecules in the supersaturated vapor and in the liquid, using the first minimum of the pair distribution function (based on the distance between the centers-of-mass of 2 N₂ molecules) at 5.7 Å  as a cutoff distance used to define neighbors. The results plotted in Fig. 7(a) show that these distributions allow us to identify vapor-like molecules (with less than 6 neighbors) and liquid-like molecules (with 6 or more neighbors). The variation of the number of liquid-like molecules in the system as a function of the entropy during the nucleation process is shown in Fig. 7(b). As with Ar, the number of liquid-like molecules exhibits a steady increase, corresponding to the increase in the size of the droplet as nucleation proceeds. An examination of the snapshots (see Fig. 8) obtained during the $μVT−S$ simulations sheds light on the formation of the liquid droplet as the entropy of the system decreases. This allows us to identify the size of the critical droplet, which contains $1006±30$ N₂ molecules for system 4, $735±25$ molecules for system 5, and $638±20$ molecules for system 6. We observe the expected correlation between the size of the critical droplet and the height of the free energy barrier of nucleation, with the lowest supersaturation of system 4 leading to the largest free energy of nucleation and to the biggest critical droplet, which are reached for the lowest critical entropy. This set of results show that S can be used as the reaction coordinate for the nucleation process in a molecular fluid and that no adjustment or redefinition of the RC is necessary.

We now apply the $μVT−S$ method to a second molecular fluid and study the nucleation of a liquid droplet in CO₂ at $260K$⁠. CO₂ is similar to N₂ in many aspects, having a similar shape and exhibiting similar types of intermolecular interactions (with a quadrupole moment for CO₂ more than 6 times greater than for N₂). Carrying out the $μVT−S$ simulations for the formation of the liquid droplet at 260 K allows us to draw a comparison with the results on N₂. This is because the conditions of temperature for N₂ and CO₂ can be considered as similar, if one considers reduced temperatures, with respect to the critical temperature, for each compounds (0.85 for CO₂ and 0.83 for N₂). As for Ar and N₂, we perform a series of $μVT−S$ simulations to sample the entropic pathway leading towards the formation of a liquid droplet of CO₂ of a critical size for 3 supersaturations. We show in Fig. 9 the free energy profiles of nucleation so obtained. Similarly to the nucleation of a liquid droplet in N₂, the profiles show that decreasing the entropy of the system during the $μVT−S$ simulations allows the system to overcome the free energy barrier and that S can be reliably used to drive the nucleation process in the molecular fluid of CO₂.

The range of entropies sampled along the entropic pathway is shown to depend on the extent of the supersaturation. As for Ar and N₂, the entropy for the starting point, the metastable supersaturated vapor, decreases at high supersaturations as the supersaturated vapor becomes more dense. The entropy for which the top of the free energy barrier is reached also increases as the supersaturation becomes high, as a result of the smaller size of the critical droplet. Both factors account for the reduced range of entropies sampled for system 9 (high supersaturation) when compared to system 7 (low supersaturation). This leads to the following values for the entropies S_c corresponding to systems containing a liquid droplet of a critical size, with $Sc= 3.55±0.02$ (kJ/kg)/K (system 7) and increasing to $Sc=3.59± 0.02$ (kJ/kg)/K (system 8) and finally to $Sc=3.61± 0.02$ (kJ/kg)/K (system 9) for the highest supersaturation. This increase in S_c at high supersaturations is correlated with a decrease in the height of the free energy barrier of nucleation which reaches $26± 3kBT$ for system 7 (lowest supersaturation), $21±4kBT$ for system 8, and $14±2kBT$ for system 9 (highest supersaturation). Since the reduced temperatures and the supersaturations for the nucleations of N₂ and CO₂ are similar, we attribute the much smaller free energy

barriers of nucleation, when compared to those found for N₂, to the stronger intermolecular interactions in CO₂.

We now analyze the nucleation process in terms of the number of liquid-like molecules in the system. Closely following the process used to define this order parameter for the other systems, we determine the distributions for the number of neighbors in the supersaturated vapor and in the liquid using a cutoff distance of 5.85 Å  to characterize neighboring molecules. The distributions plotted in Fig. 10(a) allow us to distinguish between vapor-like molecules (6 neighbors or less) and liquid-like molecules (7 neighbors or more). Using this criterion, we determine the number of liquid-like molecules during the nucleation process as the entropy changes along the nucleation pathway. The results are plotted in Fig. 10(b) for system 7 and show a smooth increase in the number of liquid-like molecules as a function of the entropy during the formation of the liquid droplet. The formation of the liquid droplet can be best captured by looking at the snapshots of Fig. 11, which show the system for decreasing values of the entropy during the nucleation process. These snapshots show the increase in the size of the nucleus and reveal that the critical droplet contains $395±33$ molecules for system 7. The critical sizes obtained for the other systems are the following: $290± 28$ molecules for system 8 and $242±20$ molecules for system 9. As for Ar and N₂, the critical sizes are correlated with the free energy barriers of nucleation, with the highest barrier of system 7 corresponding to the biggest critical droplet and the lowest supersaturation.

In this work, we develop the $μVT−S$ simulation method to calculate the free energy of the nucleation of atomic and molecular fluids along entropic pathways. We achieve this by working in the grand-canonical ensemble, which allows for the direct evaluation of the entropy during the simulations and for the study of the nucleation process at a given supersaturation $Δμ$⁠, thereby providing a direct connection with the predictions from the classical nucleation theory. The simulation protocol for $μVT−S$ consists of performing a series of umbrella sampling simulations that use S as the reaction coordinate for the nucleation process. The results presented here show that the method provides a picture of the nucleation process that is consistent with the classical nucleation theory and with that found in previous work, including, e.g., the dependence of the height of the free energy barrier and of the size of the critical nucleus on supersaturation. They also establish that the $μVT−S$ approach can be applied to atomic as well as molecular systems without any modification. Furthermore, our simulations reveal the correlation between the entropy of the system and the droplet size and show how supersaturation impacts the entropic pathway visited during the nucleation process, with a narrowing of the range of the entropy visited at high supersaturation. Our findings also allow us to identify the critical value for the entropy associated with the formation of a nucleus of a critical size. Further work presented in the next two papers of the series will focus on the development of the method for multi-component systems and for heterogeneous nucleation processes.

Partial funding for this research was provided by NSF through CAREER Award No. DMR-1052808.