A random walk through the dynamics of homogeneous vapor-liquid nucleation

Device optical modeling A method of calculating rates of homogeneous vapor-liquid nucleation based on Langevin dynamics of a few relevant degrees of freedom on a free-energy surface is proposed. The surface is obtained here from simulation and from a semi empirical expression. The mass and friction coefﬁcients are derived from atomistic umbrella-sampling molecular-dynamics simulations. The calculated nucleation rate agrees with atomistic simulations for one particular state point of the Lennard-Jones ﬂuid. The present method is about four orders of magnitude more computationally efﬁcient than the direct atomistic simulation of the transmission coefﬁcient. © 2005 American Institute of Physics

 You may deposit the accepted manuscript immediately after acceptance, using the credit line formatting below  You may deposit the VOR 12 months after publication, with the credit line and a link to the VOR on AIP Publishing's site Format for credit lines  After publication please use: "This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in (citation of published article) and may be found at (URL/link for published article abstract).  Prior to publication please use: "The following article has been submitted to/accepted by [Name of Journal]. After it is published, it will be found at Link."  For Creative Commons licensed material, please use: "Copyright (year) Author(s).
This article is distributed under a Creative Commons Attribution (CC BY) License."

I. INTRODUCTION
Classical nucleation theory ͑CNT͒ provides a simple analytic expression for the rate of homogeneous vapor-liquid nucleation by relying on several assumptions. One such assumption is that the growing liquid droplet is spherical in shape and of uniform liquid density. The free energy of nucleation at a given temperature and supersaturation is a function of a single coordinate, such as the droplet radius or number of molecules in the droplet, and depends on the chemical potential difference between the coexistence and supersaturated vapor phases, and the interfacial tension between the two phases. Further simplifications are made for the nucleation rate, namely, that it depends upon the freeenergy barrier height, a steady flux over the barrier, growth by monomer addition, ideal gas kinetics for the vapor, and complete removal of product.
Although CNT has been found to reproduce experimental estimates of the size of critical nuclei, for many of the same experiments it does not predict the nucleation rate as accurately. 1 This result could be due to the neglect in CNT of additional degrees of freedom that may be relevant to the nucleation process; e.g., fluctuations in the density or shape of the liquid droplet may also be important. A method of calculating the nucleation rate that incorporates other degrees of freedom such as these could therefore be useful.

II. BEYOND CLASSICAL NUCLEATION THEORY
According to classical nucleation theory, the Gibbs free energy ⌬G CNT ͑R͒ of a spherical droplet of liquid density l and radius R in a supersaturated vapor at a pressure p is This is measured relative to the uniform vapor phase, and ⌬ = l ͑p͒ − v ͑p͒ is the difference between the chemical potential of the liquid l and that of the vapor v at the pres-sure p, and ␥ † is the liquid-vapor surface tension for a planar interface at liquid-vapor coexistence. The second term in Eq. ͑1͒, which is positive, dominates for small R, while the first term, which is negative, dominates for large R. Therefore, ⌬G CNT ͑R͒ starts off positive, has a maximum, and eventually becomes negative for very large values of R. The CNT expression for the nucleation free energy, Eq. ͑1͒, can be generalized to include changes in the droplet density ͑or some other degree of freedom͒ as well as the droplet radius R.
In order to do this, we consider the constrained thermodynamic potential 2 F͑N,V,A͉,p,␥,T͒ ϵ F͑N,V,A,T͒ + pV − N + ␥A, ͑2͒ where = ͑ v , T͒ and p = p͑ v , T͒ are, respectively, the chemical potential and pressure of a reservoir of density v at a temperature of T, and ␥ ϵ ␥͑ , v ͒ is the surface tension between the reservoir and a subsystem of volume V containing N particles ͑␥ is assumed to be independent of geom-etry͒. This thermodynamic potential is equal to difference in Helmholtz free energy between a total system of volume V tot and particle number N tot consisting of a subsystem of volume V ͑and surface area A͒ containing N particles and a reservoir of density v at a temperature T and a total system of uniform density v at temperature T.
We can show that and, therefore, that But ␥A is identically zero if both the subsystem and reservoir have density v , and so is true for all A. By setting N = V = 0, it is easy to show that this constant is zero. In other words, if N / V = v , the total system consists entirely of vapor and therefore this constrained thermodynamic potential, which is the total Helmholtz free-energy relative to a total system of uniform density v at temperature T, must be zero.
To make it clear that we are considering changes in free energy in going from a uniform system of density v to a nonuniform system consisting of N particles in a volume V surrounded by a reservoir of density v , we define the quantity If we take = l , the liquid density, ⌬F͑N = l V , V , A ͉ , p , ␥ , T͒ is the free-energy change to nucleate a droplet of volume V ͑and surface area A͒ in a vapor of density v . We now make a couple assumptions. Firstly, we assume that the free energy of the subsystem is equal to the freeenergy density of a uniform system of density at temperature T multiplied by the volume of the system, Secondly, we assume that the surface tension is given by the interpolation formula 3,4 where S 0 ϵ ‡ / v † and S ϵ v / v † , and the superscript " ‡" denotes the spinodal while the superscript " †" denotes liquidvapor coexistence. This equation gives ␥͑ v , v ͒ =0, ␥͑ l † , v † ͒ = ␥ † , and ␥͑ , ‡ ͒ = 0. We also specialize to a spherical geometry; i.e., V =4R 3 / 3 and A =4R 2 , which leads to With f͑ , T͒ from an equation of state for the fluid and ␥ † from simulation or experiment, ⌬F͑ , R͒ can be calculated as a function of and R. For the Lennard-Jones ͑LJ͒ fluid, for example, accurate equations of state are available 5 and so is simulation data for the liquid-vapor surface tension. The pathway to nucleation in classical nucleation theory ͑CNT͒ is the slice through this surface at = l , although the CNT free-energy values use ␥ † instead of ␥͑ l , v ͒.
The free energy ⌬F͑ , R͒ is plotted for the LJ fluid under two different sets of thermodynamic conditions in Fig. 1. Allowing fluctuations in the spherical droplet's density as well as its radius R, it can be seen from Fig. 1 that close to the critical point-the same occurs at high supersaturationchanges in at constant R cost little energy near the barrier to nucleation. This result indicates the potential importance of in determining the nucleation rate.
Accurate empirical free-energy surfaces ⌬F͑q͒, in terms of a few relevant degrees of freedom, q = ͑q 1 , q 2 , ... ,q N ͒, such as those in Fig. 1, may provide a more accurate means than CNT for calculating nucleation rates while avoiding costly atomic simulations, provided a way of carrying out dynamics on the surface is available.
One shortcoming with this method is that it does not provide a way of determining a priori which degrees of freedom are relevant to the reaction. However, comparisons with the experimental or simulation results may indicate whether a particular choice of coordinates is suitable.

III. LANGEVIN DYNAMICS
One way of carrying out dynamics on the free-energy surface ⌬F͑q͒ is to use the Langevin equation, 6,7 where m i is the mass corresponding to the coordinate q i and v i = q i . There are two contributions to the force on the coordinates q. The first is a force due to variations in the freeenergy surface ⌬F, The second is a random force which has a mean of zero and a delta-function correlation in time, ͗⌫ i ͑t͒͘ = 0, ͗⌫ i ͑t͒⌫ j ͑tЈ͒͘ = 2␦ ij ␦͑t − tЈ͒.

͑13͒
The friction coefficients ij can be related to the strength of the random force via and to the diffusion coefficients D ij via In principle, m i and ij can vary with q. It should also be noted that m i and ij are masses and friction coefficients associated with generalized coordinates and do not have a simple interpretation in terms of atomic masses or the friction acting on individual particles.
The parameters m i and ij can, however, be estimated from molecular dynamics ͑MD͒ in which the time dependence of the coordinates q are measured. In order to obtain accurate statistics for m i and ij as functions of q, it is also necessary to confine the system to a narrow range of values of q, particularly regions of configuration space that the system rarely explores, such as near free-energy barriers. This can be achieved, for example, with umbrella-sampling molecular dynamics.
From the computational point of view, the Langevin dynamics simulations for a few relevant degrees of freedom q are trivial compared to carrying out nucleation rate calculations directly in fully atomistic simulations. The relative efficiency in obtaining the parameters required for the Langevin equation from an atomistic simulation compared to simulating the rate directly is less clear. We note however that the quantities we report below were obtained from a MD trajectory of 50 000 time steps, and we estimate the relative statistical error in them and in the nucleation rate itself as less than 10%. This can be compared with the direct atomistic simulations of the rate constant using the reactive flux method by ten Wolde et al. 17 who used 5000 independent trajectories of 4000 steps, or 2 ϫ 10 7 MD time steps in total, which gave a nucleation rate with a relative statistical error of 67%. On the basis of these numbers, to get the same statistical accuracy one would require 2 ϫ 10 4 times as many MD steps for the direct simulation of the transmission coefficient than for the diffusion constant and inertial mass. In other words, the present Langevin method is more than four orders of magnitude more efficient than the conventional direct simulation of the nucleation rate. Further it may be possible to estimate the Langevin parameters by analytic approximations or other means as functions of the coordinates q as well as of the thermodynamic conditions. If this can be done, the Langevin dynamics method offers even greater savings. As a first step towards implementing the method and understanding how m i and ij vary with q and the thermodynamic conditions, in this work these quantities have been obtained from atomistic simulations.

A. General equations
For a system with total potential energy where U 0 ͓r N ͔ is the potential energy of the reference ͑unbi-ased͒ system as a function of atomic coordinates r N and V͓q͑r N ͔͒ is a biasing ͑umbrella͒ potential, which is a function of the coordinates q, the force acting on a particle i in the system is

B. Continuous cluster definition
In the nucleation problem, we are interested in coordinates like q = ͑ , R͒, or equivalently q = ͑N , R͒, where , R, and N are the droplet density, radius, and particle number, respectively.
In order to apply Eq. ͑17͒, N and R must be differentiable functions of the atomic coordinates r N . Therefore, continuous analogs of the more conventional "discrete" cluster definitions 8,9 are required. We have taken for the cluster particle number and for the cluster radius, where the center of mass is the weight of an atom in the cluster is and the continuous analog of the Heaviside function is The weight w i is a measure of the "contribution" of particle i to the cluster particle number ͑w i Ϸ 1 for a particle near the center of the cluster and 0 Ͻ w i Ͻ 1 for a particle near the surface͒.
The upper limit on the sums, N d2 , is the number of particles in the largest discrete cluster in the system, defined according to Stillinger's criterion 8 with cluster cutoff distance q c2 . According to this criterion, two particles belong to the same cluster if separated by less than the cutoff distance ͑see Fig. 2͒. The quantities q c , d, n c , and M are parameters to be chosen. To obtain a physically meaningful value for N, n c is chosen such that w i Ϸ 1 for particles in the center of the cluster. q c and q c2 are chosen such that the discontinuity in the force is small when a particle joins the discrete cluster.
This definition of the continuous cluster particle number N is the continuous analog of the "liquidlike" cluster particle number used by ten Wolde and Frenkel. 9 Their liquidlike cluster criterion is almost the same as Stillinger's criterion, except that only liquidlike particles are counted when determining whether two particles are part of the same cluster. A particle is liquidlike if it has more than a certain number of nearest neighbors ͑in their work 9 this number was four͒, with particles defined as being nearest neighbors if separated by less than the cluster cutoff distance. According to this liquidlike cluster criterion, particles near the center of the cluster contribute a value of 1 to the particle number while particles near the surface may not contribute at all. Similarly, according to our continuous cluster definition, particles near the center contribute w i Ϸ 1 to the particle number while particles near the surface contribute 0 Ͻ w i Ͻ 1, with w i approaching zero as the particle moves away from the cluster.

C. Force on the cluster
For the quadratic biasing potential the contribution from the biasing potential to the force on particle i is where and Here, r ij ϵ r i − r j ,r ic ϵ r i − r cm ͑and similarly for r jc and r kc ͒, r ij ϵ͉r ij ͉, and fЈ͑r ij ͒ = df͑r ij ͒ / dr ij . Note that the force only acts on the N d2 particles in largest discrete cluster in the system.

V. APPLICATION TO THE LJ FLUID
The umbrella-sampling MD method of the previous section were applied to the LJ fluid with intermolecular potential truncated and shifted at 2.5. All quantities hereafter are in reduced LJ units unless otherwise stated.
The MD simulations of 864 particles at a constant temperature of T = 0.741 and a constant pressure of P = 0.01202 were carried out using the extended system methods of Nosé, 10 Hoover, 11 and Andersen, 12 and the reversible integration algorithm of Martyna et al. 13 The MD time step used was 0.01, where = ͱ / ͑m LJ 2 ͒.
It should be noted that Eq. ͑7͒ for the internal pressure P int in Ref. 13 must be modified to include the contribution of the biasing potential ⌬P int . For V͑N , R͒ given by Eq. ͑23͒, this contribution is where V is the volume of the cubic simulation box and d is the dimensionality.

A. Equilibrium free energies
The umbrella-sampling MD method was tested by comparing the change in the Gibbs free energy for nucleation as a function of cluster particle number N with results obtained by ten Wolde and Frenkel 9 from Monte Carlo ͑MC͒ simulations under the same thermodynamic conditions. For these calculations, the biasing potential acted only on N and not on R; i.e., c R = 0 in Eq. ͑23͒ or As explained in Sec. IV B, ten Wolde and Frenkel used a different, liquidlike, definition of the cluster particle number than the continuous definition used here. The liquidlike cluster number will be denoted hereafter as N l .
The parameters in Eqs. ͑21͒ and ͑22͒ were chosen so that N and N l in our simulations matched as closely as possible over the entire range on N ͓n c = 9.0, q c = 1.5 ͑same as the cutoff used by ten Wolde and Frenkel͒, q c2 = 2.0, and d = 0.02͔. Figure 3 shows a typical configuration near the freeenergy barrier to nucleation from our simulations.
Following an equilibration period of 1 000 000 MD cycles, simulations of 3 000 000 MD cycles in length were carried out in each of 17 approximately evenly spaced windows for N spanning the range from 0 to approximately 350. where q c Ͻ q c2 . Here, the discrete cluster particle numbers are N d = 5 and N d2 = 6, respectively.
The umbrella-sampling MD algorithm does not work for small N, where the identity of the largest cluster can fluctuate between different clusters. So, for small N, a 10 000 000cycle simulation without the biasing potential was carried out. The length of these simulations was chosen to match the statistical accuracy of the MC simulations of ten Wolde and Frenkel. 9 They carried out simulations of 250 000 MC cycles in length in approximately the same number of umbrellasampling windows as we used. However, the "diffusion constant" in the N l coordinate in the MC simulations was estimated to be 7-8 times as large as that in MD simulations with the same time step as ours. 9 Therefore, it was necessary to carry out MD simulations approximately eight times as long to obtain the same accuracy. Figure 4 shows contour plots of two-dimensional histograms n͑N , N l ͒ from two different simulations, one with the quadratic biasing potential centered at small N and the other at large N. It can be seen from Fig. 4 that N and N l match very closely at small N. At large N, N is generally larger than N l and for a given value of N, N l can take quite a large range of values.
The joint probability distribution function p 2 ͑N , N l ͒ for the reference system without the biasing potential was calculated from the distribution function with the biasing potential, p 2 ͑N , N l ͒, using where subscript "1" notes an average for the biased system. The full joint probability distribution function p 2 ͑N , N l ͒ was obtained for the entire range of N and N l by using the multiple-histogram method 14 generalized to two dimensions. The singlet distribution functions were then obtained by using The Gibbs free energy was calculated using ⌬G͑q͒ = − k B T ln͓p 1 ͑q͔͒, q = N or N l . ͑31͒ Almost identical results were obtained when the singlet distribution function for the reference system, p 1 ͑q͒, where q = N or N l was calculated directly from the singlet distribution with the biasing potential, p 1 ͑q͒, and joined together using the multiple-histogram method in one dimension. Since the joint distribution function p 2 ͑N , N d ͒ was not calculated in the simulations while the singlet distribution function p 1 ͑N d ͒ was, this latter method was used to obtain the free energy for N d . The Gibbs free energy as a function of N, N l , and N d is shown in Fig. 5. The free energy as a function of N does not extend to zero because we only calculated N for the largest cluster in the system. It was instead assumed that ⌬G͑N͒ coincided with ⌬G͑N l ͒ for N = N l = 15.  Figure 5 shows that there is a discrepancy between our results for ⌬G͑N l ͒ and those of ten Wolde and Frenkel. 9 However, this discrepancy is likely due to the extreme sensitivity of ⌬G͑N l ͒ to supersaturation under the particular set of thermodynamic conditions studied. 9 ͑Fitting their data for the free-energy barrier ␤⌬G͑N l * ͒ as a function of supersaturation S to ␤⌬G͑N l * ͒ = −397.59+ 2341.71/ S − 4635.03/ S 2 + 3240.05/ S 3 , it can be seen that changing S by just 2% from S = 1.53 ͑where P = 0.01202͒ to S = 1.50 changes ␤⌬G͑N l * ͒ from 57.5 to 63.5. Another possible contribution to the discrepancy between our results of ⌬G͑N l ͒ and those of ten Wolde and Frenkel 9 could be the accumulation of small errors arising from the joining of nonoverlapping umbrellasampling windows by eye in Ref. 9. This could be particularly problematic at small N l due to the steepness of the free-energy curve and to the statistical error in the data.
The crucial point to note from Fig. 5 is that ⌬G͑N l ͒ Ϸ ⌬G͑N͒, for all N and N l . On the other hand ⌬G͑N d ͒ is significantly different. The free-energy barrier as a function of N d is also lower than that for N and N l . There is, in fact, no reason why ⌬G͑N d * ͒, ⌬G͑N l * ͒, and ⌬G͑N * ͒ should be the same, since there is no one-to-one correspondence between these coordinates.

VI. NUCLEATION RATE
As a simple first test of the Langevin dynamics method for calculating nucleation rates, the calculations were restricted to those of the Langevin equation for one degree of freedom ͑N͒ with constant coefficients; i.e., mN ͑t͒ = − m␥Ṅ ͑t͒ + f͓N͑t͔͒ + ␣⌫͑t͒, ͑32͒ where ␣ = k B Tm␥ and ␥ = / m. The coefficients m and ␥ were estimated from the autocorrelation function for the velocity of N, Ṅ = dN / dt, It should be noted that this equation does not hold for t =0, where d͗Ṅ ͑t͒Ṅ ͑0͒͘ / dt =0. For the constant NPT MD simulations carried out in this work using the algorithm in Ref. 13, Ṅ can be calculated analytically by using where v ij = dr i / dt − dr j / dt and p ⑀ and W are the momentum and mass, respectively, of the barostat in the MD algorithm ͑see p. 1119 of Ref. 13 for notation͒. An advantage of using a continuous coordinate such as N rather than a discrete coordinate like N l or N d is that the velocity can be calculated analytically. This avoids any ambiguity or error that may arise from taking a numerical difference to calculate the velocity of the discrete coordinate.
The simulations in which m and ␥ were calculated were carried out for various values of c N and N 0 in the biasing potential ͑see Table I͒. The simulations were 50 000 MD cycles in length. They were effectively carried out at con-stant volume by setting the inertia on the barostat at a very large number. This was done because the fluctuations in the simulation box volume were found to produce a lot of noise in ͗Ṅ ͑t͒Ṅ ͑0͒͘. The simulation box volume was set at its average value in constant NPT simulations using the same biasing potential. Since these simulations were effectively at constant volume, then the second term in Eq. ͑34͒ can be ignored.
The autocorrelation functions ͗Ṅ ͑t͒Ṅ ͑0͒͘ obtained from these simulations are shown in Fig. 6. k B T / m was taken to be the value of ͗Ṅ ͑t͒Ṅ ͑0͒͘ at t = 0, while ␥ was obtained by fitting ͗Ṅ ͑t͒Ṅ ͑0͒͘ for small but nonzero t ͑t = 0.02-0.06͒ to an exponential function of the form of Eq. ͑33͒. The values of k B T / m and ␥ obtained from the simulations are shown in Table I. The data for k B T / m from Table I are plotted as a function of ͗N͘ in Fig. 7.
It can be seen from Table I that the friction coefficient ␥ acting on the N coordinate is roughly independent of N and approximately equal to 50 −1 . On the other hand, from Fig.  7, it can be seen that k B T / m varies linearly with N ͑or at least  with ͗N͘ in our simulations͒. Therefore, the "mass" of the N coordinate is inversely proportional to the actual mass of the cluster.
The reason for this dependence follows directly from the definition Eq. ͑34͒. Ignoring the barostat term, it is straightforward to show that where f 1 ϵ͐drfЈ͑r͒, f 2 ϵ͐drfЈ͑r͒ 2 , and the pair and first Legendre transform of the triplet distribution function appear. This result uses the facts ͗v i ͘ = 0 and ͗v i · v j ͘ =3␦ ij k B T / m LJ , and assumes that fЈ͑r͒ is sharply peaked at q c , that the droplet is of liquid density l , that the correlation functions are bulk quantities, and that the contribution from the surface of the droplet may be ignored. The particular consequence of this last assumption is that the mass of the coordinate N is inversely proportional to N, in agreement with Fig. 7. One concludes from these results that the dominant contribution to Ṅ at small times comes from changing the number of bonds of atoms in the interior of the droplet rather than from atoms joining or leaving the surface of the droplet. This in turn suggests that the behavior of the liquidlike cluster definition and its continuum analog used here may not be as appealing as the original Stillinger definition. In order to calculate the rate coefficient for nucleation k, the reactive flux method 15,16 was used, in conjunction with Langevin dynamics. In this method, the rate coefficient is given by a time-dependent function, which should be approximately constant and equal to the forward rate of reaction for time t intermediate between the molecular relaxation time and the reaction time. Here, q is an order parameter that discriminates between the reactant and product states, ͑in our nucleation calculations, q ϵ N͒, q * is the transition state separating the reactant and product states, ͑q Ͻ q * denotes reactants and q Ͼ q * denotes products͒, and the subscript "eq" denotes an equilibrium average. The transition state theory ͑TST͒ estimate of the rate coefficient k TST , which assumes that trajectories that cross the transition state q * from reactants to products never recross the transition state surface, is the t → 0 + limit of k͑t͒: The transmission coefficient is defined as the ratio of k to k TST , The rate k was taken as the plateau value of k͑t͒. The transmission coefficient should be governed largely by diffusive behavior near the free-energy barrier. So from Langevin dynamics simulations on different free-energy surfaces, ⌬G͑N͒ that are the same near the barrier but differ in the stable states should be approximately the same.
Since we had accurate free-energy surfaces for ⌬G͑q͒, where q = N or N l from both our atomistic simulations and those of ten Wolde and Frenkel, 9 these surfaces were used in the Langevin dynamics calculations rather than a free-energy surface derived from the equation of state in Eq. ͑10͒. The barrier region of our simulation curve for ⌬G͑N͒ was fitted to a couple of different polynomial functions that differed in the stable states ͓curves ͑4͒ and ͑5͒ in Fig. 8͔. Due to the discrepancy between our results for ⌬G͑N l ͒ and those of ten Wolde and Frenkel, 9 we did the same for their simulation curve for ⌬G͑N l ͒ ͓curves ͑1͒-͑3͒ in Fig. 8͔. In order to use Eq. ͑32͒, both ␥ and m had to be fixed at constant values. ␥ was fixed at its approximately constant value of 50.0. Since should depend mostly on what happens at the barrier, m was fixed at its value at the barrier in the free-energy curves from the atomistic simulations. Therefore, we took m = 7.0ϫ 10 −4 2 in the simulations on curves ͑1͒-͑3͒ ͑since N * Ϸ 280͒ and m = 6.0ϫ 10 −4 2 in the simulations on curves ͑4͒ and ͑5͒ ͑since N * Ϸ 330͒.
The calculated values of k, k TST , and on the freeenergy surfaces ͑1͒-͑6͒ in Fig. 8 are shown in Table II. As expected, is approximately independent of the differences in the stable states. The similarity between the value for surfaces ͑1͒-͑3͒ and ͑4͒ and ͑5͒ is because, despite the discrepancy between our results for ⌬G͑N l ͒ and those of ten Wolde and Frenkel, 9 the curvature of our surface and theirs near the barrier is approximately the same.
The sensitivity of k, k TST , and to the values of m and ␥ chosen can be estimated from Kramers' rate theory in the high-friction limit, which predicts that k ϰ 1/m␥, k TST ϰ 1/ ͱ m, and ϰ 1/ ͱ m␥.
The total nucleation rate, ͑number of droplets produced per unit volume per unit time͒, as shown by ten Wolde et al., 17 can be calculated by using The required input for this equation are and ͉͗Ṅ * ͉͘ from the Langevin simulations, the vapor density v , and the freeenergy barrier with respect to the vapor phase, ⌬G͑N * ͒, from the full empirical free-energy surface. For almost the same thermodynamic conditions of the LJ fluid studied in this work ͑T = 0.741 and P = 0.012 rather than P = 0.01202 as used in this work͒, ten Wolde et al. 17 carried out fully atomistic MD simulations to estimate the nucleation rate. For the liquidlike cluster particle number coordinate N l , they found that = 0.003± 0.002 and ͉͗Ṅ l * ͉͘ = 76.2 −1 . The errors in their value of are very large, because its value is very small and it takes a lot of costly simulation time to improve the statistics. ͑In fact, using several different generalizations of the reactive flux method, they calculated values of that included = 0.03± 0.03, 0.004± 0.004, 0.011± 0.009, and 0.004± 0.01.͒ In comparison, we obtained for the continuous cluster particle number coordinate N the values of Ϸ 0.014 and ͉͗Ṅ * ͉͘ = 26.0 −1 . If it is assumed that ⌬G͑N l * ͒Ϸ⌬G͑N * ͒, which we found from our simulations, and the values of ␤⌬G͑N l * ͒Ϸ59.4 and v = 0.0188 are taken from ten Wolde et al., 17 it is found that the total nucleation rate calculated by our method, 5.4ϫ 10 −29 −3 −1 , agrees with the value of ͑3.5± 2.3͒ ϫ 10 −29 −3 −1 obtained by ten Wolde et al. 17 from their atomistic MD simulations.

VII. DISCUSSION
This simple first test of the Langevin dynamics method for calculating the nucleation rate demonstrates its potential utility. In this work the free-energy surface and mass and friction coefficients used in the Langevin calculations were obtained from atomistic computer simulations. Compared to the direct atomistic simulation of the rate coefficient using the reactive flux method, it appears that more than four orders of magnitude less computer time is required to obtain the Langevin coefficients for comparable statistical accuracy. In addition, there are potentially other means for estimating the free-energy surface and mass and friction coefficients that may even avoid the atomistic simulations.
Our first test of the method for the LJ fluid demonstrated that, for at least one set of thermodynamic conditions, the nucleation rate is accurately estimated by assuming the diffusive behavior of a single coordinate, the number of particles in the cluster. In simulations of nucleation dynamics of the Ising model under similar thermodynamic conditions to the LJ fluid studied in this work ͑away from the critical point and at low supersaturation͒, Pan and Chandler 18 found that the number of particles in the cluster was "good" reaction coordinate and obtained a reasonable estimate of with a simple equation based on diffusive behavior of N near barrier. ͑This analysis was analogous to ours in terms of Langevin dynamics͒.
However, as discussed in Sec. II, other degrees of freedom may be important under other conditions ͑e.g., near the critical point and at high supersaturation͒. These degrees of freedom can readily be incorporated into our Langevin dynamics calculations, along with the coordinate dependence of mass and friction coefficients. 7