J . Phys. Chem. 1995, 99, 10446- 10453
10446
Core and Surface Properties of Bcc Clusters of SeF6 and Dynamics of Homogeneous Nucleation to Monoclinic Phase. A Molecular Dynamics Study Lawrence S. Bartell* and Shimin Xu Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109 Received: March 15, 1995; In Final Form: April 20, 1995@ Clusters of SeF6 were observed in molecular dynamics simulations to undergo a solid-state transition when sufficiently undercooled. Analyses were performed on clusters ranging in size from 100 to 500 molecules to examine how faithfully properties of cores of clusters reproduce properties of bulk matter. It was found that configurational energies and coefficients of rotational diffusion of molecules in cores closely approximated those of the bulk. If this had not been true, there would have been no basis for calculating the degree of undercooling at which nucleation was observed. Excess surface energies obtained in the process were roughly half the heat of sublimation per unit area of a molecular layer and perhaps twice the surface free energy. Nucleation from the bcc phase to monoclinic was seen to occur in each of 17 cooling runs made on 150 molecule clusters. Nucleation always began in the interior of the clusters and was mononuclear, always leading to single crystals of the monoclinic phase. From the times of nucleation could be derived nucleation rates at temperatures of 88 and 78 K. Lifetimes of bcc molecules in a given orientation inferred from the runs agreed in order of magnitude with those derived from Raman and NMR spectra and yielded the same activation energy. Frequencies of rotational jumps from the old phase to the new were estimated from the coefficients of diffusion. By applying the classical theory of homogeneous nucleation, it was possible to derive the kinetic parameter o,, for the bcc-monoclinic interface. In nucleation theory it represents the interfacial free energy. This quantity was approximately 0.19 of the heat of transition per unit area from bcc to monoclinic, or about two/thirds that of the corresponding ratio proposed for freezing transitions. W e conclude that clusters serve as convenient and realistic models for studying both surface properties and the dynamics of bulk systems.
TABLE 1: Lennard-JonesParameters for Seven-Site Intermolecular Potential Function for SeFh Molecules
Introduction The concept of homogeneous nucleation as proposed by Turnbull and Fisher’ for the freezing of pure undercooled liquids is now widely accepted for systems free from heterophase contaminants. Moreover, despite the fact that solid-state transitions in the bulk are believed generally to occur by heterogeneous nucleation,’ the homogeneous mechanism has recently been shown to apply to solid-state phase changes in very small, undercooled single Even so, quantitative investigations of homogeneous nucleation have proven to be troublesome to perform. For reasons discussed by Turnbull’.’ and others, experiments are subject to many difficulties, and realistic deductions based on pure theory are hampered by the complexity of the problem. Therefore, it has been worthwhile to turn to computer simulations of phase transitions to explore what insights they might convey. A number of enlightening molecular dynamics (MD) simulations of the freezing of systems of Lennard-Jones spheres and of metal atoms have already been carried o ~ t . One ~ , ~of these,h in particular, has examined the kinetics of production of precritical nuclei in an enormous system of Lennard-Jones spheres and has convincingly supported the idea of homogeneous nucleation as opposed to, for example, spinodal decomposition. It would be desirable, however, to begin to acquire a broader base of information about characteristics of more general systems, including those of polyatomic molecules, and systems undergoing solid-state transitions. The present investigation was stimulated by results of a new experimental m e t h ~ d ~ ,for ’ . ~measuring rates of homogeneous nucleation by studying transformations taking place in undercooled molecular clusters. By this technique both freezing and @Abstractpublished
in
Advance ACS Ahtrruct,,
June 1, 1995
0022-365419512099-10446$09.0010
atom pair
F-F
Se-Se
Se-F
u. A E . kJ/mol
2.92 0.213
3.78 I .38
3.32 0.542
certain solid-state transitions can be followed by electron diffraction while the clusters are in flight in a supersonic beam. In several of the systems studied, projections of experimental results based on the Turnbull-Fisher nucleation theory’.9 suggested that at sufficiently deep undercoolings, nucleation rates would become large enough to make transitions accessible to MD simulations. These projections turned out to be correct. In one example, the case of clusters of SeF6 transforming from bcc to monoclinic crystals, MD simulations provided information that was very helpful in interpreting the experiments. Preliminary results of the simulations were reported in the experimental paper.’ A more refined and comprehensive treatment of the solid-state behavior inferred from the simulations is presented in the following. Procedure Computational Details. Molecules were taken to be rigid octahedra with Se-F internuclear distances of 1.67 A.1° Simulations were carried out with a modestly modified version of the program MDMPOL’ I using a seven-site intermolecular interaction function based on pairwise-additive Lennard-Jones potential functions. No cutoffs of the potential functions were applied in calculations of the clusters. Potential parameters, listed in Table 1, were adapted from Caillat’s Buckingham functions’’ and adjusted to reproduce the experimental bcc lattice constant (5.99 A at 199.3 K”) and energy of sublimation (25 W at 227 K’j) in Monte Carlo simulations (1 28 molecules, with periodic boundary conditions).
0 1995 American Chemical Society
Core and Surface Properties of SeF6 Clusters
J. Phys. Chem., Vol. 99, No. 26, 1995 10447
Our investigation of the kinetics of nucleation concentrates upon the determination of nucleation rates, how they vary with temperature, and what can be inferred from the rates according to nucleation theory. Nucleation rates were derived from constant energy molecular dynamics simulations. In view of the stochastic nature of nucleation, it was necessary to perform enough runs to acquire statistically significant results. For this purpose, 17 independent cooling runs were carried out on approximately spherical bcc clusters of 150 SeF6 molecules to investigate the distribution of times of nucleation. Clusters, whose boundaries are free, were chosen to avoid the strong influence that periodic boundary conditions impose on all but very large systems. It had already been established that molecules in clusters of hexafluorides spontaneously pack into crystalline lattices closely approximating those observed experimentally for bulk systems even when clusters contain as few as four dozen molecule^.'^^'^ Therefore, it was supposed that the interior of a cluster would provide a reasonably realistic model of the bulk. In the following, several additional thermodynamic and kinetic lines of evidence corroborate our contention that the cores of clusters as small as those selected for the present study closely resemble bulk matter. These were acquired from MD computations carried out on quasispherical clusters of 100, 150, 250, 350, and 500 molecules. From these runs could also be deduced the molecular motions in the bcc phase needed in the analysis of the rate data obtained. Simulations were carried out as follows. Starting configurations were perfect bcc structures, initially adjusted to 138 K. Each of the clusters was given a somewhat different initial lattice constant and/or different set of initial molecular velocities to prevent excessive correlation in trajectories in the different runs. At each new temperature clusters were placed in a heat bath, rescaling each of the first 1000 time steps, and then runs were continued at constant energy for another 4000 time steps (for 11 runs) or 9000 time steps (for the other six runs), with all time steps being 1Ofs. During the constant energy stages, thermodynamic averages of such quantities as potential energy and temperature were accumulated, and coordinates, quaternion parameters, etc., of all of the molecules were saved every 20 time steps. In cooling, and also in subsequent heating runs, temperatures were changed in steps of 10 deg. Cluster structures were monitored after each set of 5000 or 10 000 time steps was completed to detect phase changes. These were recognized by viewing the molecular arrangement of the configuration at the last time step with the aid of the proprietary software MACSPIN, by examining Pawley projection^'^^'^ of molecular orientational distributions, by comparing pair-correlation functions of Se atoms, and by studying the angular distribution functions P ( 0 ) of neighboring Se atoms.I6 Further computational details are reported elsewhere.I6 Determination of Nucleation Rate. The procedure to calculate the nucleation rate from observed rates of transformation of an ensemble of clusters was as follows. If the fraction J ( t ) of clusters surviving in the bcc phase at time t is fl(Q =N,/(N,
+ 4)
(1)
with Nl and N2 representing the numbers of clusters in the initial and final phases, then f i ( t ) is expected to evolve with time according to
fIC0 = exp[-JV,(t - to>l
(2)
where to is the time at which N2 is zero, V, is the effective volume of a cluster in which nucleation can occur, and J is the rate of production of critical nuclei per unit volume per unit
-14
-0 E 3
-18
-22
-26 30
110
70
150
230
190
T, K Figure 1. 1. Transitions between the bcc and monoclinic phases of (SeF6)N clusters as indicated by the changes in slopes of the configurational energies per molecule in cooling (solid) and heating (dashed) curves. Filled circles, N = 150; triangles, N = 250; filled squares, N = 350: diamonds, N = 500.
-20.9
1 0
,
,
15
30
,
,
,
1
45
60
75
90
time, ps Figure 2. Evolution of the configurational energy of a 150 molecule cluster Z6.4/10 (Table 3). A critical nucleus appeared 41 ps after the cluster emerged from a bath temperature of 78 K or 5 1 ps after it entered the heat bath.
time once a steady state of embryonic precursors, has been established. Under the conditions of the present calculations our estimate3 of the time lag corresponding to the establishment of the steady state is short compared with the time over which the nucleation takes place and, therefore, is unimportant in analyzing the data. The nucleation rate at a given temperature was calculated from the slope of the curve lnfi(t) where Nl in eq 1 represents the number of clusters surviving in the bcc phase at the beginning of runs at that temperature. Weighting of the data in determinations of the decay rates of the bcc population is outlined in the Appendix. Phase transitions were readily recognized from the plots of configurational energy vs temperature (as illustrated in Figure 1 for the set of runs with clusters of 150-500 molecules), as well as from the other diagnostic techniques mentioned above. Heating curves were reproducible because no specific nucleation events were needed. Surface molecules with their naturally greater disorder were organized sufficiently like those of the warmer phase to serve as sites of the phase change. During warming the change always began at the surface. Cooling curves, on the other hand, display greater or smaller overshoots of the transition temperature depending upon chance nucleation events, and these stochastic events always began in the interior. In order to recognize the time of nucleation in individual clusters, an examination of their temporal behavior is, of course, required. A typical record of the time evolution of configurational energies of a cluster is shown in Figure 2. A complete
10448 J. Phys. Chem., Vol. YY, No. 26, 1995
Bartell and Xu
display of the records is given elsewhere.I6 It can be seen, as expected for a finite system, that the energies fluctuate about a constant mean value until a transition to a lower temperature phase begins to occur. The time at which a critical nucleus appears in each run was somewhat subjectively identified as the time step at which the instantaneous configurational energy began a concerted descent toward that of the more efficiently packed monoclinic phase. In a study of the nucleation in TeF6 reported elsewhere: such an identification was fine-tuned by a more detailed diagnosis. Two questions about the derivation of nucleation rates arise. One involves the rather sparse statistics of the runs. If the number of available nucleation events were very large, any time at which none of the clusters had transformed would be a legitimate value of to (assuming that the nucleation time lag were sufficiently short, a safe assumption, here). When only a few events are available at a given temperature, however, it might be considered that a time somewhere between the entrance into the heat bath at that temperature and the time of emergence 10 ps later might serve as the best time origin, and that the (marginal) information in the time distribution provided by the sample itself in its fit of eq 2 with to a free parameter, would give as valid a fit as any. We chose to reckon nucleation times from time of entry into the heat bath. Unconstrained fits of the MD results via the logarithmic form of eq 2 led to values of to of 0.4 and 3 ps for the 88 and 78 K nucleation events. Arbitrarily. then, we constrained to to be zero. The effect on derived nucleation rates was a small fraction of the statistical error. The other question concerns the appropriate assignment of a value to V,. A plausible answer is suggested by the fact that never in many dozens of MD simulations of freezing and solidstate transitions in cooling runs carried out in this laboratory have transitions been observed to take place in the surface layer of clusters; they have always begun in the interior. This fact accords well with intuition in view of the extra disorder in the surface layer. Therefore, we arbitrarily take V, to be the volume of the cluster interior according to the following criterion. In a cluster of N molecules the fraction of molecules on the surface is very nearly
F = 3(4d3i11"3[1 - 0.5(4~/3N)"~]*
(3)
or almost 213 for a cluster of 150 molecules. However crude this approximation may be, its effect on the quantities to be derived from the MD nucleation rates is modest. Inference of Jump Frequency. Although the classical theory of homogeneous nucleation is not rigorously correct, it gives an excellent qualitative description of nucleation. Therefore, it is worthwhile to interpret derived nucleation rates in terms of its framework. A key quantity in the classical theory is the frequency of molecular jumps from the original phase to the evolving nucleus of the new phase. For some transitions such as freezing the jumps are envisaged to be primarily translational in character. For the nonreconstructive phase transition investigated here, however, the principal change is one of molecular orientation. As shown by Raynerd ef al.'* and by Pawley and c o - w ~ r k e r s ,one-third '~ of the molecules in the bcc phase of hexafluorides must rotate 60" about one of their 3-fold axes in order to generate the new (monoclinic) structure. Rotational frequencies of the appropriate 60" molecular jumps from the bcc orientations in the bcc phase were estimated as follows. Coefficients of rotational diffusion De at each of the temperatures investigated were computed from the internal motions of the core molecules whose coordinates were recorded at intervals in the simulation. These quantities repre-
sented ensemble averages based on the time evolution of directions of a given bond in each molecule. Specifically, De was derived from the slope of the initial mean-square diffusional components of the bond directions from the initial directions, according to
(4) where the librational component had already reached a steady state and where the result had been averaged over as many time origins as the saved data permitted. It is proposed that the diffusion in the bcc structure progresses by jumps from one stable configuration to another of the configurations 90" away from the initial one as suggested by the approximate localization It is of orientations displayed in Pawley projections. lh." reasonable to suppose that the lowest barriers for the jumps involve paths through intermediate orientations in which molecules are rotated 60" about one of their 3-fold axes inasmuch as such orientations direct fluorines of the jumping molecule toward the hollows between fluorines of its neighbors. It was also envisaged that half of such intermediates decay back to the original configuration and half complete the full 90" jump. Because there are four such 3-fold axes, only one of which leads to the specific orientation of the new phase, the rotational jump frequency was taken to be the reciprocal of half the time for molecules to diffuse 90" about one of the four axes, or
v, = 8D$2
(5)
in hertz if Do is expressed in rad%. Inference of Interfacial Free Energy. According to the classical theory of nucleation it is possible, in principle, to determine the free energy of the interface between the bcc and monoclinic phases from the nucleation rate. The theory applied is essentially the theory of homogeneous nucleation developed for freezing by Turnbull and Fisher,'.9 modified as suggested by El Adib et ul.3.20to apply to solid state transitions. That is, the frequencies of translational molecular jumps from the liquid phase to the solid phase in the Tumbull-Fisher theory are replaced by the frequencies of rotational jumps from orientations in the initial crystalline phase to those of the final phase. Our treatment closely parallels that of our less complete study of nucleation in TeF6.4 Before nucleation theory can be applied it is necessary that the degree of undercooling be known. When applying the capillary theory involved, as explained by Reiss, Mirabel, and Whetten,21the undercooling should not be reckoned from the size-dependent transition temperature of the clusters undergoing nucleation, but from the bulk transition temperature. To be consistent in the MD treatment, that temperature should correspond to the one implied by the potential parameters adopted for the simulation. From the mode of calibration of the potential function there is no guarantee that the model's bulk system transforms at the same temperature as observed for the real system. It is not a simple matter to find the transition temperature in simulations employing periodic boundary conditions to model the bulk, partly because the boundary conditions would have to be quite complex because of the change in symmetry of the structure during the transition. Equally troublesome is the fact that the transformation is strongly inhibited when the system is denied a surface (as it is by the imposition of periodic boundary conditions). A possible altemative might seem to be to extrapolate the transition temperature of increasingly large clusters to infinite size according to the classical capillary of size dependence.?' Such a treatment fails, for reasons to be discussed later. For the
Core and Surface Properties of SeF6 Clusters
. .
I
I
CI
u)
c
8
xx
+
- 4
Bo 48
++
-c
%
- 3.5 - 3
XX
P
4.5
I
I
J. Phys. Chem., Vol. 99, No. 26, 1995 10449
x
++
- 2.5
+
+x
+
Bulk-like Behavior of Cluster Cores. In addition to the size independence of coefficients of rotational diffusion is other evidence that might be used to test whether molecules in cluster cores behave as they do in the bulk. Comparisons of rotational diffusion times in the clusters with experimental bulk rotational relaxation times and of core potential energies of clusters with those found in simulations of the bulk (calibrated from enthalpies of sublimation) offer possible tests. According to the picture of rotational diffusion sketched in the foregoing, times for rotational jumps over the angle Oj should be approximately
- 2
(7) 4
5
6
7 8 9 1000IT Figure 3. Upper curve and scale at right, natural logarithms of rotational diffusion coefficients (rad2/ns) of SeF6 molecules in cores of clusters. Filled symbols refer to heating runs and open symbols to cooling. Size-symbol correspondence as in Figure 1. Lower curves from spectroscopic relaxation times which are assumed to be inversely proportional to De,crosses from Raman flg lines, plus signs from NMR spectra. To superpose experiment and simulations via the model of eq 7, angular jumps from one crystalline orientation to another would have to be approximately 140" instead of the expected 90".
purposes of the present paper, then, we shall adopt the experimental value of 128 K22-24 for the bulk transition temperature from which the degree of undercooling is calculated.
Results Molecular Motions. As mentioned in the foregoing, coefficients of rotational diffusion De were determined as implied by eq 4 from the slopes of ([Ae(t)l2) curves in order to be able to compute the temperature dependent frequencies of reorientational jumps of molecules during phase changes. Because the diffusion tends to be greater in surface layers of clusters than in the interior, clusters were divided into shells defined by positions of the minima of the radial distribution functions of the Se atoms. Accordingly, the first six shell boundaries were at 6.9, 8.7, 10.9, 13.25, 15.5, and 17.5 A. The clusters analyzed were composed of 150, 250, 350, and 500 molecules and possessed 5 , 6 , 7 , and 8 shells by this accounting. It was found for the bcc phase that the coefficients De of molecules interior to the outer two shells of the clusters were indistinguishable but substantially lower than those of the outer two shells. Therefore, results from the outer two shells were excluded in calculating De. The outer two "shells" defined as above correspond approximately to the surface layer defined by our working convention, eq 3. In deriving De from ([Ae(t)12) it was necessary to disregard the rapid angular change during the approach to librational equilibrium in the first picosecond and to recognize the curvature due to partial saturation already beginning to show as early as 3 ps. Results were not discernibly dependent upon whether the data were acquired during heating or cooling runs. More importantly, they were pleasingly independent of cluster size, as would have to be true if cluster cores can be considered to be good models of bulk matter. Temperature ranges of stability for the bcc range do, of course, depend upon the cluster size as shown in Figure 1. A plot of In De vs 1/T including points from the various cluster sizes in Figure 3 exhibits satisfactory Arrhenius behavior over the whole bcc range. From this plot via eq 5 can be derived the temperature dependent orientational jump frequency Vj
= 3.30 x 10" exp(-E,lRT)
with an activation energy E, of 2.577 kJImol.
If this time is identified with the conventional Raman relaxation time z2 which is, in turn,2n times the relaxation time @ defined by Gilbert and J h i f f ~ r d ,we ~ ~can compare in Figure 3 the diffusion implied by the Raman and NMR relaxation times plotted by Gilbert and D r i f f ~ r dwith ~ ~ that from the cluster simulations. It is obvious in Figure 1 that configurational energies vary considerably with size and that they are substantially less stabilizing than bulk configurational energies. This was also true in the case of TeF6 clusters where the difference from the bulk was accounted for with modest precision by analyzing potential energies in terms of energies of core and surface molecule^.^ A somewhat more rigorous thermodynamic approach in the present work yields an even better result. If we treat clusters as spherical, each with radius r, and represent their internal energies in terms of interior and surface excess quantities, we can express their configurational energies as
where U, is the potential energy per unit volume of the cluster excluding the surface excess and U , is the surface excess potential energy per unit area at the Gibbsian equimolar surface. We will refer to the component (4/3)n13Uv as Ubulk, the intrinsic volume part excluding effects of the surface. Since a surface excess kinetic energy is not expected, we equate U, to Ea where Eo is the surface excess internal energy per unit area, whence the potential energy per molecule becomes
+
UIN = Ubulk/N ( ~ ~ ~ v , ~ ) " ~ E J N "(9) ~ by virtue of the relation 4n?/3 = Nv,, where V m is the volume per molecule. Therefore, a plot of ULV vs N -'I3 should extrapolate linearly to the intercept Ubulk /N, and this intercept should closely approximate the bulk value if cluster cores do indeed resemble the bulk. This hypothesis is tested in Figure 4 which plots cluster energies per molecule and the least squares line passing through them. In weighting the points it was assumed that standard deviations of points were proportional to (1IN)'" . Table 2 compares the intercepts with bulk values from Monte Carlo runs. An analogous procedure yielded the energies of the monoclinic cores and, therefore, the energy of the transformation from monoclinic to bcc. The result obtained agreed with that estimated from the entropy of the same transition in SF63.25to within the rather large statistical uncertainty. Unfortunately, the simulations of monoclinic clusters covered an inadequate temperature range and gave core energies that were substantially less reliable than those of the bcc phase, and too crude to be of value. Nucleation Rates. In all of the 17 runs the bcc clusters transformed to single crystals of the monoclinic phase: six at 88 K, nine at 78 K, and two at 68 K, as listed in Table 3. All clusters remained monoclinic henceforth during their cooling
Bartell and Xu
10450 J. Phys. Chem., Vol. 99, No. 26, 1995 -10,
I
I
I
I
0'5
z 5 3
7
-1
-2 -2
-301 0
I
I
I
I
I
0.05
0.1
0.15
0.2
0.25
1/
~ 1 / 3
Figure 4. Extrapolation of configurational energies of bcc molecular clusters, per molecule, to infinite size at temperatures (from bottom to top) of 120, 140, 160, and 227 K. Points, from linear fits to bcc regions of Figure I . Lines, weighted least squares fits of points. Intercepts corresponding to the intrinsic (volume) energies are compared in Table 2 with energies from bulk MC runs. Derived surface excess energies are listed in the same table.
TABLE 2: Comparison of Configurational Energies of Bcc SeFs Clusters, Extrapolated to the Bulk as in Figure 4, with Those of Bulk MC Runs, and Surface Excess Energies per Unit Area
T
-UdNh
227 169 140 120
25.0 27.4 27.9 28.6
-U
dN'
25.0 27.5 28.3 29.0
U
