5636
J. Phys. Chem. A 1999, 103, 5636-5644
Freezing of Small SeF6 Clusters: Simulations, Nucleation Statistics When Events Are Few, and Effects of Laplace Pressure Yaroslav Chushak, Prakriteswar Santikary, and Lawrence S. Bartell* Department of Chemistry, UniVersity of Michigan, Ann Arbor, Michigan 48109 ReceiVed: March 15, 1999; In Final Form: May 31, 1999
Simulations were carried out on 138-molecule clusters freezing isothermally at 130, 120, and 80 K. At 120 K, the nucleation rate was the same as in our prior simulations performed adiabatically but the final product was different. During the nanosecond period of the runs, clusters transforming adiabatically had frozen to bcc crystals while warming from 120 K to about 130 K. On the other hand, isothermal clusters at 120 and 130 K changed to monoclinic clusters after passing through the bcc phase. Clusters cooled to 80 K froze to a variety of structures. The number of molecules whose Voronoi polyhedra qualified them as being in bcc embryos grew in size erratically, and in most runs it was difficult to use the Voronoi information by itself to identify a well-characterized nucleation time. Therefore, a more discriminating criterion for the onset of nucleation was devised. The 138-molecule clusters proved to be too small to yield definitive profiles of the several order parameters characterizing the change from the liquid phase to the critical nucleus. Even though the sizes of the nuclei were not established accurately, it was clear that critical nuclei were considerably larger than the five-molecule size forecast by the classical theory of homogeneous nucleation. At the deep supercooling of the simulations, precritical and critical nuclei were extremely ramified and haphazard in molecular orientation, but the chaotically organized nuclei at 120 and 130 K quickly annealed and grew to single crystals in most clusters. Clues were found suggesting that surface molecules may participate in the formation of critical nuclei, contrary to our long-standing belief. From nucleation rates were derived the kinetic parameters σsl, the solid-liquid interfacial free energy of the classical nucleation theory, and δ, the interface thickness of Gra´na´sy’s diffuse interface theory (DIT). In addition, the effect of pressure on the DIT, a new treatment of errors, and an improved weighted least-squares analysis of nucleation data were developed.
Introduction Homogeneous nucleation of a new phase in a system of condensed matter is a process of fundamental importance in the natural sciences and technology. Although such nucleation has been the subject of scientific investigation for the half century following Turnbull’s classic studies,1-6 what happens on the molecular scale has been uncertain. Theoretical ideas have been incorporated into formulations accounting qualitatively for observations. Moreover, recent density functional treatments7-11 appear to correct some of the inadequacies of the original classical (capillary) model. Nevertheless, because nucleation is a complex process, it is not clear that all important aspects are adequately taken into account in current theoretical approaches, particularly at deep supercooling. Giving rise to such doubts are molecular dynamics (MD) simulations which are able to examine how molecules behave spontaneously in nucleation events, uninfluenced by possible theoretical biases. Until comparatively recently, only the most powerful supercomputers could address the problem realistically. With the advent of desktop workstations that are faster than yesterday’s supercomputers, it is now feasible to examine details of nucleation that have previously been obscure. The systems studied in most previous simulations of nucleation in condensed phases have been freezing atomic liquids,12-17 although a few studies of polyatomic systems undergoing freezing18-20 or solidstate transitions21 in small clusters at deep supercooling have been reported. The present investigation follows up several of these.
Hexafluorides of the group VI elements have proven to be attractive systems to examine because of their simplicity and high symmetry, because their phase behavior is interesting,22,23 and because their interaction potentials are fairly well established.24 In the preceding paper of this series,20 we reported constant energy simulations of the freezing of small supercooled liquid clusters of SeF6. In the present paper, we compare the behavior of the same clusters when simulations are carried out at constant temperature instead of constant energy. The production of undissipated heat of crystallization in our prior study did have consequences. We also compare the freezing at several temperatures and comment on effects related to the size of clusters. Computational Details Simulations. Molecular dynamics simulations were carried out on clusters of SeF6 as described in the previous paper20 of this series with a few exceptions. The principal change in conditions in the current study was to keep the system in a heat bath at the temperature of interest instead of performing a constant energy simulation. This had the effect of dissipating the heat of crystallization in a plausible way, making the process more like that which would occur in a much larger droplet that has more mass to serve as a heat sink. The algorithm for the heat bath seeks to maintain a constant temperature by adjusting molecular speeds proportionally upward or downward, as is needed to account for changes in potential energy. Therefore, it reduces speeds of the hotter molecules more rapidly than those
10.1021/jp9908911 CCC: $18.00 © 1999 American Chemical Society Published on Web 07/02/1999
Freezing of Small SeF6 Clusters
J. Phys. Chem. A, Vol. 103, No. 29, 1999 5637
of the cooler molecules. This mimics the outward flow of thermal energy from a heat source such as a crystallizing nucleus. In most cases, the clusters were made up of 138 molecules, but several clusters composed of 611 or 1722 molecules were examined to investigate size effects. A more complete report on the larger systems will be forthcoming. The interaction potential adopted was a seven-site model function described elsewhere.24 It included partial charges on atoms implied by the proprietary program Biograph/Polygraph,25 and these charges were found to be necessary to reproduce crystallographic data. Simulations were performed on clusters instead of bulk systems to avoid the imposition of periodic boundary conditions which have been found to interfere with phase transitions unless systems are very large.13,14 As will be shown, however, a certain price is paid for this choice of system. Runs were begun with the same quasispherical clusters generated in the previous paper20 by heating a cluster of 150 molecules to 220 K, then cooling (at 10 K per ns) to 140 K, a procedure shown to yield highly supercooled liquid clusters with no crystalline seeds large enough to initiate nucleation. In the process of preparation of the systems for the present runs, 12 molecules were lost by evaporation. For each of the runs at 130, 120, and 80 K, 15 independent liquid configurations were prepared from the initial 140 K melt by running an additional 5000 time steps from the previous configuration and then cooling abruptly. For technical reasons, one of the runs at 120 K was unsatisfactory. For the sake of example, a few runs were also carried out by quenching to 20 K, a temperature corresponding to a kinetic energy in the classical simulation lower than the zero-point energy of the quantum system. Durations of runs at 120 and 130 K were 1 ns (100 000 time steps). Most colder runs were shorter. Molecular coordinates were saved every 0.5 ps. Diagnoses of Phase Change. The principal crystalline phases encountered in the runs were bcc and monoclinic (space group C2/m). Diagnoses of the phase change from liquid to crystal were based on Voronoi polyhedra18,26-28 and Pawley projections,29 as described previously. Translational differences between solid and liquid phases are identified by the Voronoi polyhedra, whereas Pawley projections distinguish monoclinic from bcc. One difference from prior papers was that a reference set of Voronoi polyhedra was generated for monoclinic clusters to complement the reference set existing for bcc clusters. Because the translational differences between the bcc and monoclinic structures are minor, either set can be used with little penalty. Although Pawley projections clearly distinguish liquid from crystalline clusters and are very effective in differentiating bcc from monoclinic once the phase changes are well under way, they are not nearly as sensitive in identifying crystalline embryos as Voronoi polyhedra. In the present analyses, an additional diagnosis of transformation was also made, based on the number of “bulklike solid” molecules, molecules which pass the Voronoi test and which, in addition, are surrounded closely by at least 12 molecules also passing the test for bcc or monoclinic environments. Reasons for adopting this criterion will become evident. Application of Nucleation Theory. It is assumed that the nucleation rate J can be represented by the conventional expression1
J ) A exp(-∆G*/kBT)
(1)
As discussed in detail in the previous paper of this series, we apply two variants of the prefactor A and two variants of the
free energy cost of producing a critical nucleus, ∆G*. For the former, we invoke both the classical (molecular diffusion)1,30,31 and the Grant-Gunton (thermal diffusion)32 formulations, and for the latter, those based on the classical capillary nucleation theory (CNT)1,30 and Gra´na´sy’s diffuse interface theory (DIT).33-36 For the CNT, the free energy barrier to nucleation is given by
∆G* ) 16πσsl3/[3(∆Gv + w′)2]
(2)
if it is assumed that the critical nucleus is spherical, where σsl is a kinetic parameter supposed to represent the interfacial free energy per unit area for the solid-liquid boundary, ∆Gv is the free energy of freezing per unit volume of the solid, and w′ arises from the change in free energy accompanying a change in the surface area of the freezing droplet, radius r0, during nucleation,37 or
w′ ) PL(Fl - Fs)/Fl
(3)
where PL is the Laplace pressure 2σl/r0 inside the cluster and the F’s represent densities. A value of σsl can be derived from the nucleation rate via eqs 1 and 2. Our procedure to determine Gra´na´sy’s interface thickness parameter δ by applying his DIT was described in the previous paper of this series.20 That treatment was based on Gra´na´sy’s formulation for bulk systems neglecting the effect of any external pressure imposed on the liquid phase. Our revised treatment includes the correction for the substantial Laplace pressure exerted on the liquid phase as outlined in the Appendix. Inference of Nucleation Rate from Simulations. It is assumed that the fraction N/N0 of unfrozen clusters in which a critical nucleus has not yet formed is given by the first-order rate law
N(t)/N0 ) e-JVc(t-t0)
(4)
where t0 is the time lag to achieve a steady state of precritical embryos31 and Vc is the volume of the cluster considered to be effective in nucleation. In all prior computations on SeF6 clusters in this laboratory, this volume has been taken to be the total volume minus the volume occupied by surface molecules, because nucleation has always been initiated in the interior of the clusters, never on the surface itself. Until the present paper, the fraction of molecules assigned to the interior (core) of a cluster of N molecules had been based on the formula
Fcore ) 1 - 3(4π/3N)1/3[1 - 0.5(4π/3N)1/3]2
(5)
but in the present paper, an alternative convention was adopted of counting as core the number of molecules possessing Voronoi polyhedra. The resultant count is very nearly the same as that implied by eq 5. As will be explained subsequently, in the present work with very small clusters, we found suggestive clues that, even though a critical nucleus almost certainly begins life in the interior of a cluster, it is likely to incorporate some of the surface molecules before it has grown to the critical size. Therefore, our criterion of considering only core molecules probably underestimates Vc and, hence, overestimates J. To apply eq 4, we plot ln(Nn/N0) vs the time tn at which the nth nucleation event in the set of N0 clusters has taken place. From the slope of the plot is derived the quantity VcJ, and from the intercept, the time lag t0. In prior papers, we took Nn to be the number of clusters remaining liquid after the nth nucleation event, a nonoptimal choice of no importance when N0 is large
5638 J. Phys. Chem. A, Vol. 103, No. 29, 1999
Chushak et al.
TABLE 1: Volumes Per Moleculea in the Liquid, bcc, and Monoclinic Phases of 138-Molecule Clusters of SeF6 T, K
a
phase
130
120
80
liquid bcc monoclinic
107.5 102 98
106 101 97.2
101.5 96.5
In angstroms3.
but nontrivial when N0 is small. A more appropriate value for Nn is
Nn ) N0 - n + ∆
(6)
a number greater than that adopted in our prior papers by the quantity ∆ whose value we take to be unity (but see the Appendix). Among other advantages, the revised procedure permits all of the N0 events, not just N0 - 1 of them, to be applied in the derivation of J. A justification of this counting and a suitable weighted least squares procedure to derive t0 and J are given in the Appendix. The volumes per molecule for the cores in the liquid, bcc, and monoclinic phases corresponding to the present potential function and cluster size are listed in Table 1. Because of the limited sampling available, they are less precise than volumes determined for bulk phases, but this source of uncertainty in the determinations of nucleation rates is minor in comparison with the statistical uncertainty. Results To illustrate the time evolution of the number of molecules in bcc or monoclinic aggregates in the various clusters studied, Figures 1-3 display Voronoi plots showing the increase in numbers of solid molecules with time at 130, 120, and 80 K. For comparison, the same figures also illustrate the evolution of bulklike solid molecules (defined in the foregoing section). It is evident that the onset of growth of nuclei is more distinct when the second indicator is used. To understand the figures, it should be recalled that only the interior molecules have defined Voronoi polyhedra, and there are only about 50 molecules inside the surface layer. Pawley plots at 130 and 120 K reveal that the liquid freezes first to bcc (with a few exceptions mentioned below), then transforms to monoclinic as one might expect from Ostwald’s step rule.38 This is illustrated in Figure 4 for a representative cluster at 130 K. An alternative monitoring of the transitions is given by the evolution of configurational energy of the cluster with time in Figure 5. When the liquid clusters were quenched to 80 K, however, they froze more chaotically as the Voronoi plots in Figure 3 suggest by the greater jitter in number of solidlike molecules and the leveling off of some of the plots at values considerably less than the number of core molecules. Numbers of bulklike bcc or monoclinic solid molecules plotted in the same figure give additional evidence of irregular behavior. Corroborating this erratic behavior in Figure 6 are the Pawley projections of the final products produced. Some of the clusters transformed to ordered monoclinic single crystals, but others ended up as rhombohedral, orthorhombic, apparently glassy bcc, or mixtures of structures. During the nanosecond runs, all 15 runs at 130 K froze, as did all 14 of the valid runs at 120 K and 12 of those at 80 K where the criteria based on bcc Voronoi polyhedra failed to apply. Those at the higher temperatures froze to single crystals, although one of the clusters at 120 K became a rather disordered
Figure 1. Bold lines represent the time evolution of the number of molecules identified by Voronoi polyhedra as being in bcc or monoclinic nuclei at 130 K (left-hand ordinate). Note that only core molecules, of which there are roughly 50, possess Voronoi polyhedra. For comparison, the dashed lines (right-hand ordinate) show the corresponding evolution of “bulklike solid” molecules, namely those molecules surrounded by at least 12 close bcc or monoclinic neighbors.
monoclinic crystal and another froze to orthorhombic instead of to bcc/monoclinic. The last behavior is also seen in our experimental studies of hexafluoride clusters if the cooling is fast and deep.23 How rapidly the clusters froze is illustrated by the curves of ln(Nn/N0) vs nucleation time, tn, plotted in Figure 7. Following the convention adopted in nucleation experiments, freezing was counted whenever a cluster froze, whatever its final structure, and its nucleation time was estimated from the onset of “bulklike solid” molecules. If, in addition, we retain our convention of excluding surface molecules from the cluster volume susceptible of initiating a nucleus, we obtain the time lags and rates listed in Table 2. Rates based on the total cluster volume are also listed. As will be discussed presently, the validity of the convention excluding the surface layer from the nucleating volume is called into question by results for larger clusters. Clusters quenched to 20 K behaved quite differently from those cooled more moderately. They became glassy, retaining their random structures according to Pawley projections and MACSPIN images. Values of the interfacial free energy and interface thickness derived from the nucleation rates are listed in Table 3. Physical quantities adopted in their determination are those listed in ref 20 with the exception of the liquid and bcc volumes which are listed in Table 1. No attempt was made to include in Table 3 any results from the runs at 80 K because there was no consistency in the phase formed upon freezing. Nucleation rates calculated at the three temperatures are plotted in Figure 8 where they are compared with rates
Freezing of Small SeF6 Clusters
Figure 2. Same as in Figure 1 except at 120 K.
calculated with the two different prefactors (classical and GrantGunton) and two different nucleation barriers (CNT and Gra´na´sy’s DIT), with certain common simplifications. Critical nuclei were taken to be spherical which is certainly not the case for the present MD runs, but the interfacial free energy parameter in the CNT and the δ parameter in the DIT were adjusted to make the calculated rates agree with the MD nucleation rate at 120 K. Fortuitously, this rate was virtually identical with that found in the prior constant energy simulation, the difference being considerably smaller than the statistical uncertainty in either set of runs. If another shape factor had been introduced into the formalisms, it would not have affected the trends illustrated. The differences portrayed are real differences in the natures of the different theories. Following Gra´na´sy’s suggestion39 about the essential difference between the CNT and the DIT, we constrained σsl to be constant in the former and δ to be constant in the latter. Error bars plotted for the MD rates reflect only the statistical uncertainties discussed in the Appendix. For the runs at 120 and 130 K, this is not unreasonable. For that at 80 K, the purely statistical uncertainties do not take into account the fact that the nucleation was to a variety of different solid structures and the diagnostic tool applied was valid at most for the bcc and monoclinic structures. Therefore, no confidence can be attached to the exact value of the 80 K result. In view of the uncertainties, only the DIT in combination with the Grant-Gunton prefactor can be ruled out by the simulations. Discussion Constant Energy vs Constant Temperature. It was found that the nucleation rates for runs in which clusters froze while
J. Phys. Chem. A, Vol. 103, No. 29, 1999 5639
Figure 3. Same as in Figure 1 except at 80 K.
Figure 4. Pawley projections showing the transformation of a typical cluster at 130 K from liquid to bcc to the monoclinic phase. Pawley’s dots plot the projections of bond directions of all the molecules upon a hemisphere far over the cluster. The actual disposition of the clumps of dots depends, of course, on the angle at which the plane of the hemisphere’s horizon cuts through the cluster.
kept at constant temperature in a heat bath were not discernibly different from the rate for the runs carried out adiabatically.20 What was different was the phase at the end of the runs. Clusters
5640 J. Phys. Chem. A, Vol. 103, No. 29, 1999
Chushak et al.
Figure 5. Typical time evolution of configurational energy of a cluster (no. 2) at 130 K. The energy falls as the Voronoi count of bcc molecules increases, followed by a sharp drop as the “bulklike solid” nucleus grows and transforms to monoclinic.
Figure 7. Plots of ln(Nn/N0) vs time of nucleation, tn, at 130, 120, and 80 K.
TABLE 2: Estimated Times of Nucleation in Clusters at Three Different Temperatures and Derived Nucleation Rates and Time Lagsa Nucleation Time (ps) temperature (K)
Figure 6. Pawley projections showing the varied final products of the clusters frozen at 80 K.
in which the heat of crystallization was removed rapidly froze to the bcc phase, then transformed to the more stable monoclinic habit. On the other hand, when clusters initially at 120 K froze at constant energy, they warmed to about 130 K as the heat of transition evolved and they remained in the bcc phase for the nanosecond duration of the runs. Had they transformed to monoclinic, the additional evolution of heat would have warmed them further, of course, probably by over 10° according to a prior simulation.40 This local heating inhibits the transition to monoclinic. On Inferences about Critical Nuclei. In Figures 1-3, it can be seen that in some clusters nuclei appearing to be large enough
run
130
120
80
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
570 700 280 750 940 195 450 995 425 205 500 525 620 400 500
363 350 310 690 280 560 530 310 220 420 115 740 80 580
340 400 100 175
basis core volume total volume time lag
180 90 410 240 210 70 430 300
Nucleation Rate (m-3 s-1) 5.6(1.8) × 1035 6.9(2.5) × 1035 10.8(4.2) × 1035 2.13(0.7) × 1035 2.6(1.0) × 1035 4.1(1.6) × 1035 257(50) 193(46) 98(33)
aUncertainties are standard deviations based solely on the counting statistics of eqs 13-16.
to be critical dissolve back to the liquid instead of growing immediately. This is to be expected because the fortuitous structural fluctuations that produced them in the first place can destroy them. Even though the data provide hints about the size of critical nuclei, it will be shown in later papers on much larger clusters that the present small clusters are simply too small to allow a proper characterization of sizes and properties of critical nuclei. At the outset of this study, when less computing power was available and only the classical estimate of the size of
Freezing of Small SeF6 Clusters
J. Phys. Chem. A, Vol. 103, No. 29, 1999 5641
TABLE 3: Parameters Derived from MD Nucleation Rates via the CNT and DIT with the Classical and Grant-Gunton Prefactorsa T 120 130 Tm
nucl. volume core total core total
CNT σsl(Acl) σsl(AGG) 13.24 13.86 13.36 13.92
17.01 17.42 16.47 16.87
δ(Acl) 1.79 1.88 1.68 1.76 (1.81)
DIT σsl(Acl) δ(AGG) 2.30 2.36 2.08 2.13 20.97b
ξ 1.17 1.30 2.86
Interfacial free energues σsl in interfacial thicknesses δ and correlation lengths ξ in angstroms. The parenthesized interface thickness was imposed. Standard devieations for the σsl and δ values corresponding to statistical uncertainties in nucleation rates are 0.2 mJ/m2 and 0.03 Å, respectively. b Calculated for a very large drop, not the present small cluster, because, at the small supercooling where the DIT and CNT are equivalent and therebye allow σsl to be computed via the DIT, the critical nucleus would exceed the size of our cluster. a
mJ/m2,
Figure 8. Temperature dependence of nucleation rate (m-3 s-1) for the freezing of 138-molecule SeF6 clusters, taking the nucleation volumes to be those of the cores of the clusters. Points give MD simulations with 2σ error bars including only the statistical uncertainties. Curves show calculated nucleation rates. Bold curves represent classical nucleation theory, and light curves represent Gra´na´sy’s diffuse interface theory. Solid curves are based upon the classical molecular diffusion prefactor and dashed curves upon the thermal diffusion prefactor of ref 32. The MD point at 80 K cannot be reliably compared with the other points and curves because the clusters froze to a variety of structures and, moreover, there was no fully satisfactory guide to identify the times of nucleation. Error bars indicate 3 times the 1σ errors in log J. This is because 3σ for J is as large as or larger than J itself, which would make log(J - 3σ) pathological.
critical nuclei was known (∼five molecules, far too small a value),20 it was supposed that a cluster considerably larger than five molecules would be adequate. For that matter, even though the Gra´na´sy diffuse interface theory does to some extent include the interface thickness absent from the CNT, it turrns out that the DIT predicts only about nine molecules per critical nucleus at 120 K. This also appears to be far too small. It is reasonable to ask why Voronoi polyhedra themselves are poorer indicators of the onset of freezing than our alternative “bulklike” criterion. What happens is that as the bcc embryos materialize, the Voronoi numbers climb erratically and do not always show an abrupt change. These embryos are initially composed entirely of surface molecules in filaments and sheets with all molecules in them in direct contact with the liquid. Such embryos fluctuate in shape and coalesce and break apart, and they contain many more molecules than critical nuclei were envisaged to possess. Moreover, in much larger clusters, the precritical embryos may contain many more molecules than the entire cores of the present 138 molecule clusters. However suggestive this behavior is of the approach to a spinodal, we
do not believe our results support such an interpretation. We encounter stochastic onsets of nucleation following first-order kinetics, even if the nuclei are much larger than implied by classical nucleation theory. Moreover, the free energy barriers implied by the nucleation rates are a quite a bit larger than kBT. Since the nucleating molecules identified by Voronoi polyhedra do not display a clear signal of the onset of nucleation, we choose an alternative diagnosis for the onset of nucleation. We define “bulklike bcc” molecules as those molecules passing the Voronoi test which are surrounded by at least 12 (Voronoi) bcc neighbors within the distance of the first minimum of the pair correlation function. These bulklike pockets do appear abruptly, usually followed by rapid growth. We associate this appearance with the time of nucleation. It is interesting to note that as these nuclei grow they ingest nearby embryos with entirely different orientations of the arrays of molecules. Yet, after the disparate parts have been in contact with each other, they rapidly rearrange and grow to a single crystal more often than not. Further evidence of the annealing, once the freezing is complete, is shown by the steady decrease in volume and configurational energy until a plateau is reached. If our criterion for nucleation time does not seem to be rigorously associated with the actual formation of a genuine critical nucleus, it is fair to point out that the criterion surely accords more closely with the real process of nucleation than those employed in all experimental measurements of nucleation rates in normal liquid drops. Such experimental determinations are based on the times at which whole droplets are frozen, it being assumed that the time it takes for a drop to freeze once a nucleus appears is short compared with the span of time over which the stochastic nucleation events are spread out. Possible Dependence of Nucleation Rate on Cluster Size. In all simulations carried out in this laboratory on clusters whose melts wet the solid phase, it has been seen when heating the solid that the surface melts well before the core does. Conversely, when the melt is cooled, freezing is always initiated in the interior of the cluster, not on the (more disordered) surface. Therefore, when we calculated J, the time rate of appearance of critical nuclei per unit volume, we adopted the convention that the volume considered to participate in the nucleation was the total volume of the cluster minus the volume of the surface molecules. For clusters as small as the present 138-molecule clusters, the correction is far from trivial since the correction leaves only 1/3 of the total volume. Evidence that our convention is flawed came to light when we began to obtain nucleation rates for much larger clusters (to be reported in detail in subsequent papers). When we applied our convention to clusters of 611 and 1722 molecules, the apparent rate was found to decrease markedly, falling almost 4-fold when clusters increased in size from 138 to 1722 molecules. On the other hand, if we assumed that the effective volume was the total volume of clusters, the calculated nucleation rates all fell, those of the smallest clusters decreasing the most, and the rates became roughly the same, within the statistically expected scatter. This suggests that surface molecules may well participate in bringing embryos initially formed in the interior to a size large enough to qualify as critical. Such surface molecules, if they do indeed contribute, are unfortunately not recognized by our Voronoi or “bulklike solid” counting schemes. It is not obvious that the physical space available for nucleation should have any effect upon the nucleation rate, just as long as the number of molecules available for nucleation greatly exceeds the size of critical nuclei. If the number of molecules were only marginally greater than the number
5642 J. Phys. Chem. A, Vol. 103, No. 29, 1999 required for a critical nucleus, it seems intuitively plausible that this would result in a lower, not a higher, nucleation rate. Therefore, the question to be answered is whether other sizedependent effects over and above the space available may help to account for our observations. The most obvious effect is that imposed by the Laplace pressure, a pressure increasing as the reciprocal of the cluster radius and reaching nearly 400 bar for the 138-molecule cluster. Since the solid is denser than the liquid, this pressure would lower the free energy barrier to nucleation, in conformity with eq 2 and the equations in the Appendix via the w′ term. Uncertainties in nucleation theory prevent an accurate calculation of the effect on the rate, among other reasons because the computation of the barrier from the nucleation rate depends on the prefactor adopted in its calculation. Large uncertainties in prefactors in current use are conspicuous in Figure 8. For the classical prefactor, the rate for the 138-molecule cluster at 120 K is enhanced by the Laplace pressure over that for the 1722-molecule cluster by a factor of 1.3 (CNT) or 1.7 (DIT). For the Grant-Gunton prefactor, the factor is 1.8 (CNT) or 2.8 (DIT). The truth is probably in between, closer to the result for the classical prefactor than to that for the Grant-Gunton prefactor, and perhaps closer to the DIT result than to the CNT. In any event, the Laplace pressure by itself does not appear to be large enough to account for the full differences in the nucleation rates based in our convention for the nucleation volume. We conclude that surface molecules can participate to some extent in nucleation, even if they are unlikely to be the sites at which embryos begin to form. Kinetic Parameters Derived From Nucleation Rates. Chief among these parameters are the interfacial free energy parameter σsl (derived via the CNT) and the interface thickness parameter δ (derived via the DIT). The former is of special interest to surface scientists because of the great difficulty of deriving interfacial free energies from thermodynamic measurements. The latter attracts attention because of the great difficulty in determining interfacial thicknesses by other techniques. Also determined were the time lags for nucleation at the three temperatures. Such time lags depend significantly upon the initial state of the systems. If there were no existing embryos at the time of entry into the heat bath at the various temperatures, common wisdom has time lags getting longer, the deeper the supercooling, opposite to the direction in Table 3. The existing embryos at 140 K just before immersion in the heat bath might be linked to the reversal of the lags, but the differences are not significant at the 3σ level. Unresolved problems in the theoretical and practical treatments of nucleation are evident in the values of the kinetic parameters listed in Table 3. The values tabulated depend significantly upon the assumptions made in their determination. It is fair to mention that the derivation of the Grant-Gunton prefactor was not designed to be applicable to freezing at such deep supercoolings, and the prefactor appears to be much too large, at least at the temperatures of the present simulations. An excessive prefactor requires a large nucleation barrier to bring the calculated nucleation rate down to the observed rate and, hence, leads to excessive values of σsl. On the other hand, there is also evidence that the classical prefactor is too small, particularly at deep supercoolings.41 Therefore, the truth is probably somewhere in between, probably closer to that associated with the classical prefactor. In assessing the values of σsl in Table 3, one might for sake of comparison refer to the value implied by Turnbull’s empirical relation4 (see the previous paper in this series for its application to SeF6 clusters).20 This value, depending on assumptions, ranges
Chushak et al. from 15.3 to 17.4 mJ/m2 and, hence, is closer to those in Table 3 derived via the Grant-Gunton prefactor than to the classical ones. Because Turnbull’s relation was originally based on kinetic data treated by applying the classical prefactor, not the larger Grant-Gunton prefactor, it would be inconsistent to use the values in Table 3 to lend support to the Grant-Gunton treatment. It might be argued that the lower values given at deep supercooling by the classical prefactor are consistent with the projected decrease in σsl with decreasing temperature owing to the postulated negative interfacial entropy.42,43 The DIT thickness parameter δ is of a plausible magnitude, being of the same magnitude as the interface correlation length ξ of the Grant-Gunton treatment32 estimated as described in the previous paper of this series.20 Since it is defined differently from the correlation length and other thickness parameters such as the Tolman length,44 it is not appropriate to try to identify it closely with other measures of interfacial thicknesses, particularly since one of its main virtues (of simplicity) is that it is claimed to be independent of temperature (in all cases studied so far, excluding water). The other measures of thickness surely depend on temperature. The disagreement between our values at 120 and 130 K can be attributed to an imprecision in the derived nucleation rates. It is, of course, the assumed constancy of δ that enables the DIT to provide an extrapolated value of σsl at the melting temperature where the DIT and CNT are equivalent, for ∆G* can be extrapolated to Tm if δ is constant, enabling σsl to be derived from this extrapolated barrier. Comparison With Other Simulations. A few words should be said about the substantial differences between our investigation and the excellent studies of critical nuclei by Ba´ez and Clancy16 and by ten Wolde, et al.17 Quite apart from the fact that these authors studied monatomic systems (Lennard-Jones spheres) while ours were of polyatomic molecules, the nuclei were formed by different procedures. Those of Ba´ez and Clancy were preformed with quasispherical crystalline particles implanted into the liquid. This technique imposes a bias on the character of the nuclei. Those of ten Wolde, et al. were generated and melted by adjusting a biasing potential. By contrast, nuclei generated in our investigation formed spontaneously, any way the accidental structural fluctuations happened to determine the result, and the outcome was far more haphazard and rather more difficult to model theoretically. At least part of the difference and perhaps a major part of the difference between nuclei was the much deeper supercooling in our work, a supercooling required for our systems to nucleate at a rate sufficiently rapidly for freezing to take place spontaneously on the time scale accessible to current workstations. That nuclei become increasingly ramified as the degree of supercooling increases was shown by Yang, et al.15 who worked with systems of atoms interacting with purely repulsive potentials. The degrees of supercooling in the various MD studies were about 20% for ten Wolde et al., 24-34% for Ba´ez and Clancy, and 45% and higher for our runs. It is worth mentioning that, although supercoolings in experiments tend to be low, in some cases they have exceeded 28% for metal drops,45 and in our experiments with submicroscopic drops of nonmetals (clusters of ∼104 molecules), they are typically over 30% and, in the case of benzene, a supercooling of 37% was achieved without freezing the liquid.46,47 The degree to which droplets can be supercooled before they freeze depends, of course, on the drop size, the cooling rate, and physical properties of the substance being cooled as discussed elsewhere.48 Although ten Wolde et al. stated that at a supercooling of more than 40% “one should expect the free energy barrier to
Freezing of Small SeF6 Clusters
J. Phys. Chem. A, Vol. 103, No. 29, 1999 5643
vanish for essentially all possible crystalline phases”, others, including Skripov49 who has carried out some of the most careful studies of nucleation in freezing in the last quarter of a century, assert that there is no spinodal in freezing. Certainly, our simulations at supercoolings well over 40% have shown no evidence of the vanishing of the barrier to nucleation. Concluding Remarks Improved procedures for analyzing nucleation data from simulations were devised. The most significant findings included information about the character of embryos of the new phase at deep supercooling and hints about the fraction of the volume of small clusters that is effective in nucleation. Precritical nuclei, identified as bcc aggregates by their Voronoi polyhedra, appeared that were enormously larger than the critical nuclei forecast by the classical and diffuse interface nucleation theories. These nuclei, however, consisted of fluctuating filaments and sheets so thin that all molecules were interfacial. At some point in the simulations, a bulklike ordered aggregate of molecules happened to materialize in an embryo, often followed by a rapid buildup of more solidlike molecules, signaling the onset of nucleation. Provocative evidence was also found in studies of clusters over a range of sizes that surface molecules, despite their tendency to greater disorder and contrary to our previous belief, can participate in the formation of critical nuclei. It has become clear that the present 138-molecule clusters are too small to afford definitive information about the true size and characteristics of critical nuclei. Current analyses of the results of MD simulations involving considerably larger clusters promise to provide much more detailed information about critical nuclei. Results of these analyses will be reported in the next papers in this series. Appendix Correction of Diffuse Interface Theory for Effect of Laplace Pressure. Gra´na´sy formulated his diffuse interface theory (DIT) for bulk phases. When this theory is applied to small clusters in which the Laplace pressure may be hundreds of atmospheres, it is desirable to include a correction for w′, the work per unit volume defined by eq 3. In the DIT, the free energy cost of forming a spherical critical nucleus in bulk matter at ambient pressure is given by33,34
∆G* ) -4πδ3∆Gvψ/3
(7)
where, letting η ) ∆Gfus/∆Hfus and q ) (1 - η)1/2, the quantity ψ turns out to be
ψ ) 2(1 + q)η-3 - (3 + 2q)η-2 + η-1
(8)
Equation 7 was derived from the expression for ∆G(r), the free energy of producing a nucleus, by finding the maximum of ∆G(r) with respect to the size of the nucleus. If to this expression is added the interfacial free energy change in the spherical liquid drop containing the nucleus from its initial radius of R0 to its current radius of R, or
4π(R2 - R02)σl
(9)
where σl is the surface tension of the liquid, and the maximum of the resultant equation is calculated, the barrier to nucleation is found to be
∆G* ) -4πδ3∆GvΨ/3
(10)
where, replacing Gra´na´sy’s low-pressure variables expressed by lower case Greek and Roman letters with the corresponding new variables denoted by capital letters, and defining ζ as
ζ ) w′/∆Gv
(11)
Ψ ) [2(1 + Q)H-2 - (3 + 2Q)H-1 + 1]/η
(12)
we obtain
where Η ) η(1 + ζ) and Q ) (1 - H)1/2. The expression for w′ is unchanged from that of eq 3. It can be seen that, if the solid is denser than the liquid, the free energy barrier is decreased by w′ in accord with le Chatelier’s principle. If w′ vanishes, Gra´na´sy’s original nucleation barrier is recovered. Statistical Considerations when Nucleation Events Are Few in Number. The procedure adopted to derive the nucleation rate differs from that described in our prior papers, including our revision of the convention for assigning Nn, the number of surviving liquid clusters associated with the nucleation time, tn. Two-thirds of a century ago, Peierls50 analyzed the statistics of radioactive decay when the number of disintegrating nuclei was small, focusing on the number of events expected vs the number actually found as a function of time. We have chosen the alternative approach, placing the burden of uncertainty upon the times of the nucleation events and not on Nn. The latter quantity is definite, whereas the time, being decided by chance, is uncertain. Therefore, in applying a standard weighted linear least-squares procedure, we take the time tn to be the uncertain “y” variable and ln(Nn/N0) to be the exactly known “x” variable. Although the basis is different from that of Peierls, the errors implied turn out to be quite similar. To check our selection of the most appropriate value for Nn as well as to find a suitable weight function and uncertainty to be expected for a small set of events, we resorted to a numerical analysis. A procedure was devised to construct synthetic sets of N0 events and stochastically generated times of nucleation. To do this, a set of equally probable nucleation time bins was constructed, with time ranging from zero to infinity. For a set of NT time bins, the mean time for the nth bin is
tn ) t0 - ln[1 - (n - 0.5)/NT]/JVc
(13)
Results for small N0 turned out to be insensitive to whether times were selected from 100 or 5000 time bins. For our final tests on synthetic runs, we chose 5000, from which each individual time was taken to be the center of a randomly picked time bin identified with the aid of a pseudorandom number generator. After picking a set of N0 times, the times were sorted and thereby paired with the corresponding variable ln(Nn/N0) and subjected to a weighted least-squares analysis. Various weight functions were tested with various values of N0 ranging from 5 to 75, and in each case, 30 000 or more runs of independently constructed sets of N0 events were analyzed. It was verified that the choice of ∆ to be used in eq 6 in order to recover the value of VcJ ) 1/τ built into the model runs was not the old choice with ∆ ) 0. If one focuses on the mean value of VcJ ) 1/τ, the optimum value of ∆ in eq 6 is unity, although a value of 0.62 is better if it is the mean value of τ, itself, that one desires to recover. Of the many weight functions tested, the function arctan(τ/t) gave the smallest variance of the least-squares values of 1/τ from the average 〈τ-1〉, although a number of functions weighting large t values less heavily than small gave very nearly the same variance. Moreover, the uncertainty in the τ value for an arbitrary set of N0 events was found to be
5644 J. Phys. Chem. A, Vol. 103, No. 29, 1999
Chushak et al.
στ/τ ) στ-1/τ-1 ≈ 1.10/xN0 - 3
(14)
although a denominator of xN0-2 might have been expected in view of the fact that two, not three, parameters, i.e., τ and t0, were derived in the least-squares fitting. The result, eq 14, is not far from the uncertainty
στ/τ ≈ 1.24/xN0
(15)
derived by Peierls for radioactive decay where only a single parameter, τ, needed to be determined. Peierls also pointed out that there are slightly more accurate procedures for determining the radioactive decay constant than plotting ln(Nn/N0) vs tn, but in his analysis no time lag variable t0 had to be determined simultaneously. The expectation value of the intercept, t0, associated with the least-squares fitting was found to be overestimated by the small amount
t0(least-squares) ≈ t0(true) + 1.8τ/N01.4
(16)
However, the implied correction is appreciably smaller than the uncertainty in t0 given by
σt0/τ ≈ 0.526/xN0 - 2
(17)
The almost trivial overshoot in eq 16 could be reduced significantly by applying a different weight factor, but at the cost of increasing the uncertainly in the nucleation rate. References and Notes (1) Turnbull, D.; Fisher, J. C. J. Chem Phys. 1949, 17, 71. (2) Turnbull, D. J. Met. 1950, 188, 1144. (3) Turnbull, D.; Cech, R. E. J. Appl. Phys. 1950, 21, 804. (4) Turnbull, D J. Appl. Phys. 1950, 21, 1022. (5) Turnbull, D. J. Chem Phys. 1952, 20, 411. (6) Cargilll, G. S., Spaepen, F., Tu, K.-N., Eds. Phase Transitions in Condensed SystemssExperiments and Theory (Turnbull Festschrift); Materials Research Society: Pittsburgh, PA, 1987. (7) Klupach, T. Phys. Status Solidi 1982, 109, 535. (8) Oxtoby, D. W. Nature 1990, 347, 725. (9) Harrowell, P.; Oxtoby, D. W. J. Chem. Phys. 1984, 80, 1639. Oxtoby, D. W.; Harrowell, P. J. Chem. Phys. 1994, 86, 3834. (10) Shen, Y. C.; Oxtoby, D. W. J. Chem. Phys. 1996, 104, 4233. (11) Iwamatsu, M.; Horii, K. J. Phys. Soc. Jpn. 1996, 65, 2311; 1997, 66, 1210.
(12) Mandell, M. J.; McTague, J. P.; Rahman, A. J. Chem. Phys. 1976, 64, 3699; 1977, 66, 3070. (13) Honeycutt, J. D.; Andersen, H. C. Chem. Phys. Lett. 1984, 108, 535. (14) Swope, W. C.; Andersen, H. C. Phys. ReV. B 1990, 41, 7042. (15) Yang, J.-X.; Gould, H.; Klein, W.; Mountain, R. D. J. Chem. Phys. 1990, 93, 711. (16) Ba´ez, L. A.; Clancy, P. J. Chem. Phys. 1995, 102, 8138. (17) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. J. Chem. Phys. 1996, 104, 9932. (18) Xu, S. Ph.D. Thesis, University of Michigan, 1993. Kinney, K. E.; Xu, S.; Bartell, L. S. J. Phys. Chem. 1996, 100, 6935. (19) Svischcev, I. M.; Kusalik, P. G. Phys. ReV. Lett. 1995, 75, 3289. (20) Santikary, P.; Kinney, K. E.; Bartell, L. S. J. Phys. Chem. 1998, 102, 10324. (21) Bartell, L. S.; Chen, J. J. Phys. Chem. 1995, 99, 12444. (22) Mu¨ller, U. Acta Crystallogr. 1979, A35, 188. (23) Bartell, L. S.; Hovick, J. W.; Dibble, T. S.; Lennon, P. J. J. Phys. Chem. 1993, 97, 230. (24) Kinney, K. E.; Xu, S.; Bartell, L. S. J. Phys. Chem. 1996, 100, 6935. (25) Biograph Molecular Simulations, Inc., Waltham, Mass. (26) Tanemura, M.; Hiwatari, Y.; Ogawa, T.; Ogita, N.; Ueda, A. Prog. Theor. Phys. 1977, 58, 1079. (27) Hsu, C. S.; Rahman, A. J. Chem. Phys. 1979, 71, 4974. (28) Cape, J. N.; Finney, J. L.; Woodstock, L. V. J. Chem. Phys. 1981, 75, 2366. (29) Fuchs, A. H.; Pawley, G. S. J. Phys. Fr. 1988, 49, 41. (30) Buckle, E. R. Proc. R. Soc. A (London) 1961, 261, 189. (31) Kelton, K. F.; Greer, A. L.; Thompson, C. V. J. Chem. Phys. 1983, 79, 6261. (32) Grant, M.; Gunton, J. D. Phys. ReV. B 1985, 32, 7299. (33) Gra´na´sy, L. Europhys. Lett. 1993, 24, 121. (34) Gra´na´sy, L. J. Non-Cryst. Sol. 1994, 162, 301. (35) Gra´na´sy, L. Mater. Sci. Eng. 1994, A178, 121. (36) Gra´na´sy, L. J. Phys. Chem. 1995, 99, 14183. (37) Fukuta, N. In Lecture Notes in Physics, Atmospheric Aerosols and Nucleation; Wagner, P. E., Vali, G., Eds.; Springer-Verlag: Berlin, 1989; p 504. (38) Ostwald, W. Lehrbuch der Algemeinen Chem; Engelmann: Leipzig, 1893; Vol. II-2, p 445. (39) Gra´na´sy, L. Private communication. (40) Bartell, L. S.; Xu, S. J. Phys. Chem. 1995, 99, 10446. (41) Bartell, L. S. J. Phys. Chem. 1995, 99, 1080. (42) Turnbull, D. In Physics of Non-Crystalline Solids; Prins, J. A., Ed.; North-Holland: Amsterdam, 1964; p 43. (43) Spaepen, F. Acta Metall. 1975, 23, 729. (44) Tolman, R. J. Chem. Phys. 1949, 17, 118. (45) Vinet, B.; Cortella, I.; Favier, J. J.; Desre´, P. Appl. Phys. Lett. 1991, 58, 97. (46) Valente, E. J.; Bartell, L. S. J. Chem. Phys. 1984, 80, 1451. (47) Bartell, L. S.; Sharkey, L. R.; Shi, X. J. Am. Chem. Soc. 1988, 110, 7006. (48) Bartell, L. S. J. Phys. Chem. 1992, 96, 108; 1996, 100, 8197. (49) Skripov, VC. P. First International Workshop Nucleation and NonLinear Problems in First-Order Phase Transitions, St. Petersburg, Russia, July 1, 1998. (50) Peierls, R. Proc. Royal Soc. (London) 1935, A149, 1467.
