Polymorph Selection during the Crystallization of Softly Repulsive

Oct 5, 2007 - Using hybrid Monte Carlo molecular simulations, we study crystallization from the melt of softly repulsive spheres interacting through a...
0 downloads 0 Views 321KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 12257-12262

12257

Polymorph Selection during the Crystallization of Softly Repulsive Spheres: The Inverse Power Law Potential Caroline Desgranges and Jerome Delhommelle* Department of Chemical Engineering, UniVersity of South Carolina, 301 Main Street South, Columbia, South Carolina 29201 ReceiVed: July 7, 2007; In Final Form: August 16, 2007

Using hybrid Monte Carlo molecular simulations, we study crystallization from the melt of softly repulsive spheres interacting through an inverse power law potential. We work at fixed supercooling (i.e., at a temperature 25% below the melting temperature) and consider three systems, defined by different values for the inverse power exponent n: n ) 5, n ) 6.67, and n ) 10. Modifying the value of n allows us to study the onset of crystallization in the domain of stability of the body-centered cubic (bcc) phase (n ) 5 and n ) 6.67) and in the domain of stability of the face-centered cubic (fcc) phase (n ) 10). We show that, for the three systems, polymorph selection does not take place during crystal nucleation since the structure of the critical nuclei obtained for the three systems is not well defined. However, our results demonstrate that polymorph selection takes place during the growth step since growth proceeds either into the stable bcc phase for the two smaller values of n (n ) 5 and n ) 6.67) or into the stable fcc phase for the larger value of n (n ) 10). We also show that we did not achieve complete control of polymorphism for n ) 10. The growth step gives rise to either slowly growing crystallites composed of two blocks of different structures (the stable fcc form and the metastable bcc form) or rapidly growing crystallites of the metastable bcc form.

1. Introduction Polymorphism is the ability of a molecule to crystallize in different structures or polymorphs.1,2 Polymorphs of the same molecule may exhibit very different physical properties such as, for example, the solubility and hence the bioavailability of drug molecules. This phenomenon is thus of great technological importance in the making of drugs or pigments as well as in the food industry.1 The control of polymorphism remains a longstanding and challenging issue, and a complete understanding of this phenomenon still eludes us, even for simple systems of spherical particles. In recent years, molecular simulation studies have shed light on the molecular mechanisms underlying polymorph selection during the early stages of crystallization. Because critical nuclei are small, short-lived, and thus difficult to observe directly in experiments (with the notable exception of colloidal systems, which may be observed using confocal microscopy3), molecular simulation is a natural tool to investigate this phenomenon. For instance, molecular simulations of crystal nucleation in Lennard-Jones fluid have shown that nucleation proceeded first with the formation of small clusters of the metastable body-centered cubic (bcc) form, which then quickly converted into critical nuclei whose structure was predominantly that of the stable face-centered cubic (fcc) phase.4,5 However, polymorph selection does not take always place during the nucleation step. In previous work on the Lennard-Jones system,6 we showed that, at high pressure, predominantly bcc critical nuclei first formed before evolving into either large fcc crystallites or large bcc crystallites at high enough pressure (or, more accurately, if the conditions of crystallization were chosen to lie within the domain of occurrence of the bcc form). Furthermore, we showed that, during the growth of critical nuclei, the cross-nucleation (or hetero* Corresponding author. E-mail: [email protected].

geneous nucleation) of the hexagonal close-packed (hcp) form on the fcc form could occur.7 This behavior is reminiscent of experimental observations of cross-nucleation.8,9 Our simulations on the Lennard-Jones system suggested that cross-nucleation could be controlled by modifying the degree of supercooling at which crystallization was carried out.6 In this work, we focus on understanding how polymorph selection takes place during crystallization from the melt of softly repulsive spheres interacting through an inverse power law (IPL) potential.10-16 This potential provides a continuous path from hard spheres (n f ∞) to the one-component plasma (n f 1). As such, it is an ideal system to understand the effect of the repulsion on the phase diagram and on the crystallization mechanism. The phase diagram and the domains of stability for each of the polymorphs we will come across in this work are well-known for the IPL potential.16 We simulate the entire crystallization process (nucleation and growth) from the melt at fixed supercooling (i.e., at a temperature 25% below the melting temperature) for three different values of the inverse power exponent n: n ) 5, n ) 6.67, and n ) 10. Modifying the value of n allows us to study the onset of crystallization in the domain of stability of the bcc phase (n ) 5 and n ) 6.67) and in the domain of stability of the fcc phase (n ) 10). We add that the free energy of the hcp form is very close to that of the fcc form. This paper is organized as follows: In the next section, we present the simulation methods used in this work. We performed two different types of molecular simulations corresponding to the two mechanistic steps of crystallization, that is, the nucleation and growth steps. To simulate the formation of critical nuclei, we combine a method suited to the study of activated processes, the so-called umbrella sampling technique5,17 together with hybrid Monte Carlo (HMC) simulations.18 To simulate the growth step, we carry out a stochastic molecular

10.1021/jp0753023 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/05/2007

12258 J. Phys. Chem. B, Vol. 111, No. 42, 2007 dynamics simulation.19 We then discuss the results obtained for the different values of the inverse power exponent. As pointed out by Prestipino et al.,16 the phase diagram of the IPL potential is similar to that of the Yukawa potential, often used to model charge-stabilized colloidal suspensions. We therefore compare the results obtained in this work for the crystallization of the IPL potential to those obtained in previous work for the Yukawa potential.20,21 Finally, we draw our conclusions in the last section. 2. Simulation Methods We use an IPL potential to model the interactions between particles.10-15 The potential is defined as

(σr )

Un(r) ) 

n

(1)

where we take the exponent n to be equal to 5, 6.67, or 10. Once n is fixed, the thermodynamic properties can be expressed in terms of one single quantity γ ) F*T *-3/n, where F* ) Fσ3 is the reduced density, and T * ) kBT/ is the reduced temperature (throughout this work, we will use a reduced system of units in which  is the unit of energy, σ is the unit of length, and the mass of a particle is the unit of mass). The value taken by the inverse power exponent n dictates whether the fcc polymorph or the bcc polymorph is stable.16 For large values of n, such as, for example, n ) 10, fcc is the stable phase and bcc is metastable. As n decreases, the potential becomes increasingly softer and longer ranged. The bcc phase is mechanically stable when n reaches a value of approximately 8. For smaller values of n, the looser packing of the bcc solid (with respect to the fcc solid) makes the entropy of this phase higher, and hence, at finite temperature, the bcc phase becomes more stable. In recent work,16 the bcc phase was found to be thermodynamically stable for values of the inverse power exponent less than 7.1. In this work, all simulations of crystal nucleation and growth are carried out in the NPT ensemble. We study the crystallization process in IPL liquids at a temperature 25% below the melting temperature (T * ) 1) for the pressures calculated by Prestipino et al.16 We thus perform NPT simulations at T * ) 0.75 and P* ) 203.334 for n ) 5, at T * ) 0.75 and P* ) 71.151 for n ) 6, and at T * ) 0.75 and P* ) 29.747 for n ) 10. We also use a cutoff radius rc beyond which interactions between particles are neglected, and apply the standard long-range corrections.22 We fix rc ) 5σ for n ) 5, rc ) 3.5σ for n ) 6.67, and rc ) 2.5σ for n ) 10. For a degree of supercooling of 25%, nucleation is an activated process.4,5,23-27 The system has to overcome a free energy barrier of nucleation to form a critical nucleus. In line with previous work,6,7,20,21,28,29 we choose to perform HMC simulations16 together with the umbrella sampling technique5,17 in order to study this activated, or rare, event. The umbrella sampling technique is a non-Boltzmann sampling method in which a so-called umbrella bias potential is applied. The bias potential is a harmonic function of the global parameter Q6, introduced by Steinhardt et al.30 Q6 gives a measure of the amount of crystalline order in the system studied. It is equal to 0 for a liquid and takes similar values for all the structures we will encounter in this work.30 Therefore, applying a reaction coordinate like Q6 in the bias potential will not favor the formation of a specific polymorph over another.6,7,20 Q6 can be defined as follows:5

Q6 )

[ | 4π 13

6

∑

m)-6

|]

Desgranges and Delhommelle

∑i ∑j Y6m(rˆ ij)R(rij) ∑i ∑j R(rij)

2 1/2

(2)

where rij is the vector joining two neighboring atoms (two atoms are considered as neighbors if they are separated by a distance rij smaller than a cutoff distance rq, corresponding to the first minimum of the pair correlation function in the liquid phase), rij is its norm, rˆ ij is the corresponding unit vector, Y6m(rˆ ij) is a spherical harmonics, and R(rij) is a weight function, which goes smoothly to 0 to ensure that Q6 is a continuously differentiable function (R(rij) ) (rij - rq)2 for rij < rq and R(rij) ) 0 elsewhere). We take rq ) 1.05σ for n ) 5, rq ) 1.18σ for n ) 6.67, and rq ) 1.35σ for n ) 10. We finally add that the success of this approach relies on an appropriate choice for the reaction coordinate. While previous work obtained on the nucleation of spherical particles indicate that Q6 is a suitable choice, using a single reaction coordinate to describe the nucleation event is only an approximation. This may result, for example, in an overestimate of the free energy barrier of nucleation. On the other hand, transition path sampling31 or transition interface sampling32,33 methods allow the system to sample its own, often complex and multidimensional reaction coordinate and provide unbiased insight into the nucleation process. In future work, we will check the accuracy of our results, obtained with the umbrella sampling technique, against those obtained with transition path sampling methods. There are two types of steps in the HMC simulations. An HMC step is either a molecular dynamics trajectory of 20 time steps of 5 × 10-3 reduced units (67% of the HMC steps) or a Monte Carlo volume change (33% of the HMC steps: the amplitude of the volume change is adjusted in the course of the simulation so that 50% of these moves are accepted). Velocities are drawn from a Gaussian distribution at the beginning of each molecular dynamics trajectory. The Newtonian equations of motion are then integrated using a velocityVerlet integrator. Once we have formed a critical nucleus, we study its spontaneous evolution in the absence of the bias potential. This allows us to (i) check that the critical nucleus we have formed is genuine and (ii) analyze polymorph selection during the growth step. For this purpose, we embed the system of 3000 atoms containing the critical nucleus in a supercooled liquid matrix of 22 000 atoms and obtain a system containing a total of 25 000 atoms. We then equilibrate this system according to the following two successive steps. First, the central region of 3000 particles is frozen (i.e., kept fixed), and we let the rest of the system relax during an HMC simulation. Second, we carry out an HMC simulation of the whole system while applying the umbrella sampling bias potential to the central subsystem of 3000 particles. Every 1000 HMC steps, we store an equilibrated configuration of the system (we repeat this procedure 10 times). Each of these 10 configurations will then be used as a starting point of an unbiased trajectory (i.e., with the bias potential switched off) in the NPT ensemble. Each unbiased trajectory is generated using the stochastic molecular dynamics simulation method proposed by Attard19 for the isothermal-isobaric ensemble. This method consists of deterministic time steps (using Newton’s equations of motion), alternated with stochastic steps based upon the Boltzmann distribution (it is very similar to the popular method proposed by Andersen34). The introduction of stochastic steps could, in principle, alter the growth and dissolution mechanisms of the

Crystallization of Softly Repulsive Spheres

J. Phys. Chem. B, Vol. 111, No. 42, 2007 12259

nuclei. We tested this method against the purely deterministic molecular dynamics used in previous work7 and obtained similar results for the structure of the crystallites. During each deterministic time step, we integrate the equation of motion with a velocity-Verlet integrator and a time step of 5 × 10-3 reduced units. We generate 10 unbiased evolutions of a critical nucleus. We observe the dissolution of the nucleus in the surrounding liquid for five of them and the growth of the nucleus in the remaining five free evolutions. The 5:5 ratio so obtained demonstrates that the crystal nuclei we have formed are critical nuclei. Throughout nucleation and growth, we need to identify the structural identity of each atom. In this work, we consider three different structures: fcc, hcp, or bcc. In order to carry out this analysis, we first need to compute the local order parameters q6, q4, and wˆ 4 for each atom30 (in this calculation, we use a shorter cutoff distance to define pairs of neighboring particles: the cutoff distance is equal to 0.9σ for n ) 5, 1.1σ for n ) 6.67, and 1.25σ for n ) 10). We now briefly present how we proceed (this analysis is detailed in previous work7,20,21). Once we have computed the local order parameters, we need to identify which atoms are solid-like and which atoms are liquidlike. This is done according to the method proposed by ten Wolde et al.,5 on the basis of the amount of correlation between the local parameters of two neighboring molecules. When an atom is identified as solid-like, we then assign a structural identity to this atom, according to the value of q4 and wˆ 4 for this atom. An atom is identified as fcc-like if q4(i) > 0.16 and w4(i) < -0.05, bcc-like if q4(i) < 0.08 or if q4(i) < 0.11 and w4(i) < -0.025, and hcp-like if q4(i) > 0.1 and w4(i) > 0.04. We use the same rules for the three systems studied. We will compare the results obtained on the IPL potential to those reported in previous work on the Yukawa potential.20,21 The Yukawa (or screened-Coulomb) potential is often used to model the interactions between colloidal particles in chargestabilized colloidal suspensions.35 If we take a ) n-1/3 (where n ) N/V is the macroion number density) as the unit of length, then the Yukawa potential can be written as

U(r) )

Φ e-λr/a a r/a

(3)

where λ is the dimensionless screening parameter. Experiments have shown that adding more salt to the suspension destabilizes the bcc structure and stabilizes the fcc structure. We can mimic this effect by tuning the value of the dimensionless screening parameter λ of the Yukawa potential. The phase diagram and the relative stability of the polymorphs of the Yukawa system have been determined in previous work.36,37 For a small value of the screening λ, the stable phase is the bcc phase, and the fcc phase is metastable. For a large value of λ, we have the inverse situation: the fcc phase is stable, whereas the bcc phase is metastable. We add that the IPL potential (when n f 1) and the Yukawa potential (when λ f 0) tend toward the same limit: the one-component plasma model. 3. Results and Discussion We qualitatively obtain similar results for the two lower values of n (n ) 5 and n ) 6.67). Therefore, in this section, we will only present the results obtained for n ) 6.67 together with those obtained for n ) 10. We first comment on the results obtained for the nucleation step. We recall that we combined the umbrella sampling technique with HMC simulations to drive the formation of critical nuclei for both systems. Using this

Figure 1. Outside view of the critical nucleus for (a) n ) 6 and (b) n ) 10. Only the particles identified as solid-like are shown. Gray: fcc; yellow: hcp; red: bcc.

method, we are able to determine the free energy barrier of nucleation as well as the structure of the critical nuclei. Snapshots of critical nuclei are shown in Figure 1 for n ) 6.67 and n ) 10. The nucleation mechanism for the two systems share several common features. First, nucleation starts in both cases with the formation of small bcc clusters, in agreement with the theoretical predictions of Alexander and McTague.38 As the crystal nuclei grow and get closer to the critical size, fcc-like and hcp-like particles appear (see Figure 1). This nucleation mechanism is similar to that reported for simple (Lennard-Jones) fluids in previous density functional theory (DFT)39-41 and computational studies.4,5,7 Both types of studies concluded that small nuclei have a predominantly bcc structure and evolve into stable fcc crystals as they grow larger. DFT studies39-41 also showed that a substantial amount of bcc order always remained at the interface between the fcc crystal and the liquid. Second, in both cases, the structural identity of the critical nuclei is not well defined. Our analysis reveals that we have roughly the same numbers of particles identified as fcclike, bcc-like, or hcp-like. Third, the height of the free energy barrier of nucleation is estimated to be the same for the two systems (17 ( 2kBT). These results show that the nucleation mechanism is similar for both systems. Clearly, polymorph selection has not yet taken place at the end of the nucleation process. We now turn to the results obtained for the growth step. The two systems exhibit very different behaviors. Let us start with the results obtained for n ) 6.67. We plot in Figures 2 and 3

12260 J. Phys. Chem. B, Vol. 111, No. 42, 2007

Figure 2. n ) 6. Evolution of the number of fcc, bcc, and hcp particles during the growth of a critical nucleus (t ) 0 denotes the time at which the bias potential is switched off and the system starts to evolve freely).

Figure 3. n ) 6. Crystallite obtained at the end of the growth trajectory shown in Figure 2: (a) outside view and (b) cross-section. Same legend as that in Figure 1.

the results obtained for one of the growth trajectories (all growth trajectories for n ) 6.67 led to the same qualitative result). Figure 2 shows the evolution with time of the numbers for each type of particles (fcc-like, bcc-like, and hcp-like) within the cluster. Throughout the growth step, the number of bcc-like particles increases at a much faster rate than the number of hcplike or fcc-like particles. Snapshots of the crystallite obtained at the end of the growth trajectory (t ) 100) are also presented

Desgranges and Delhommelle in Figure 3. The crystallite is mostly composed of bcc-like particles, with hcp-like and fcc-like particles scattered on the surface (a few of them appear as defects in the core of the crystallite). These results demonstrate that, for n ) 6.67 (as well as for n ) 5), polymorph selection takes place during the growth step since the structure of the crystallites becomes predominantly that of the bcc form as they grow. Why does the structure of all the crystallites evolve toward bcc rather than fcc? Indeed, while the bcc polymorph is the stable phase for conditions close to coexistence, the fcc form happens to be the stable phase at the pressure and temperature of crystallization. According to Prestipino et al., fcc becomes the stable phase for a ratio P*/T * > 76.945 for n ) 6.67, and the simulations presented here for n ) 6.67 are carried out at a ratio P*/T * ) 95. The fact that the bcc phase has a lower chemical potential than the liquid phase under the conditions of crystallization is more significant and more relevant to the understanding of the growth process. In other words, the conditions of crystallization lie within the occurrence domain of the bcc form. As seen during the nucleation process, the formation of bcc-like particles is kinetically favored. Since the bcc form is more stable than the surrounding liquid, the crystallite keeps growing in the bcc form. We move on to the results obtained for n ) 10. Figure 4 shows the two different behaviors observed during the growth step. Growth trajectories lead either to slowly growing crystallites whose structure is predominantly that of the stable fcc phase (Figure 4a-d) or to rapidly growing crystallites whose structure is predominantly that of the metastable bcc phase (Figure 4e,f). We first examine Figure 4a-d. Panels a and c of Figure 4 show that, after a time of approximately 40 reduced units, polymorph selection starts to take place as the number of fcc-like particles begins to increase at a faster rate than those for hcp-like and bcc-like particles. The two snapshots presented in Figure 4b,d show that the crystallites can be described as being composed of an fcc block, with planes of hcp particles trapped within this core, and a bcc block (on the left of Figure 4b and on the right of Figure 4d). We attribute the presence of hcp planes to successive heterogeneous nucleations (or cross-nucleations) of the hcp polymorph on the stable fcc polymorph and vice-versa. Cross-nucleation in this system is favored by two factors. First, the two forms have almost the same free energy. Second, the two forms share common geometrical features: cross-nucleation takes place on the structurally compatible (111) planes of the fcc form. The overall growth rate of the crystallite is, however, limited by the fact that two blocks of different structures (fcc and bcc) coexist. We now examine Figure 4e,f. In this case, polymorph selection takes place earlier (around t ) 25), and growth proceeds in the metastable (and kinetically favored) bcc phase. The number of bcc-like particles increases at a fast rate, and the growth mechanism is similar to that observed for lower values of n. The fcc-like and hcp-like particles remain scattered on the surface, and we do not observe the formation of an fcc block as in Figure 4a-d. As a result, the overall growth rate is much faster. While we observed polymorph selection to take place during the growth step (most often with the formation of the stable fcc form), we did not achieve complete control of polymorphism for n ) 10 since the metastable bcc form also formed for kinetic reasons. We finally compare the results obtained on the IPL potential to those reported in previous work on the Yukawa potential.20,21 As discussed above, there are several features common to the phase diagrams for the two models. Our results show that the crystallization mechanisms followed by the two models are very similar. For a small value of n for the IPL potential (or a small

Crystallization of Softly Repulsive Spheres

J. Phys. Chem. B, Vol. 111, No. 42, 2007 12261

Figure 4. n ) 10. (a,c,e) Evolution of the number of fcc, bcc, and hcp particles during the growth of the crystal nucleus. (b,d,f) Cross-sections of the crystallite obtained at the end of the growth trajectory shown in panels a, c, and e, respectively. Same legend as that in Figure 1.

value of the screening parameter λ for the Yukawa potential), the stable polymorph is the bcc form, and crystal growth leads to the formation of large bcc crystallites. For a large value of n for the IPL potential (or a large value of λ for the Yukawa potential), the stable polymorph is the fcc form, and crystal growth leads to the formation of predominantly fcc crystallites partly composed of a block of the metastable bcc form. 4. Conclusions In this work, we used HMC molecular simulations to study crystallization from the melt of softly repulsive spheres interacting through an IPL potential. We worked at fixed supercooling (i.e., at a temperature 25% below the melting temperature) and considered three systems, defined by different values for the inverse power exponent n: n ) 5, n ) 6.67, and n ) 10. This allowed us to study the crystallization of IPL fluids in the domain of stability of the bcc phase (n ) 5 and n ) 6.67) and in the domain of stability of the fcc phase (n ) 10). We showed that, for the three systems, polymorph selection did not take place during crystal nucleation since the structure of the critical nuclei obtained for the three systems was not well defined. For the two smaller values of n (n ) 5 and n ) 6.67), our results demonstrated that polymorph selection took place during the growth step since growth proceeded into the stable bcc phase.

Similarly, for the larger value of n (n )10), polymorph selection took place during the growth step since growth proceeded into the stable fcc phase. However, we did not achieve complete control of polymorphism. The growth step gave rise either to slowly growing crystallites composed of two blocks of different structures (the stable fcc form and the metastable bcc form) or to rapidly growing crystallites of the metastable bcc form. The crystallization mechanisms observed for the IPL potential were found to be in excellent agreement with those previously observed for the Yukawa potential. Acknowledgment. Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund, for partial support of this research. References and Notes (1) Bernstein, J. Polymorphism in Molecular Crystals; Oxford University Press: Oxford, 2002. (2) Bernstein, J.; Davey, R. J.; Henck, J. O. Angew. Chem., Int. Ed. 1999, 38, 3441. (3) Gasser, U.; Weeks, E. R.; Schofield, A.; Pusey, P. N.; Weitz, D. A. Science 2001, 292, 258. (4) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. Phys. ReV. Lett. 1995, 75, 2714. (5) ten Wolde, P. R.; Ruiz-Montero, M. J.; Frenkel, D. J. Chem. Phys. 1996, 104, 9932.

12262 J. Phys. Chem. B, Vol. 111, No. 42, 2007 (6) Desgranges, C.; Delhommelle, J. Phys. ReV. Lett. 2007, 98, 235502. (7) Desgranges, C.; Delhommelle, J. J. Am. Chem. Soc. 2006, 128, 10368. (8) Yu, L. J. Am. Chem. Soc. 2003, 125, 6380. (9) Chen, S.; Xi, H.; Yu, L. J. Am. Chem. Soc. 2005, 127, 17439. (10) Hoover, W. G.; Ross, M.; Johnson, K. W.; Henderson, D.; Barker, J. A.; Brown, B. C. J. Chem. Phys. 1970, 52, 4931. (11) Hoover, W. G.; Gray, S. C.; Johnson, K. W. J. Chem. Phys. 1971, 55, 1128. (12) Laird, B. B.; Haymet, A. D. J. Mol. Phys. 1992, 75, 71. (13) Agrawal, R.; Kofke, D. A. Phys. ReV. Lett. 1995, 74, 122. (14) Agrawal, R.; Kofke, D. A. Mol. Phys. 1995, 85, 23. (15) Dubin, D. H. E.; Dewitt, H. Phys. ReV. B 1994, 49, 3043. (16) Prestipino, S.; Saija, F.; Giaquinta, P. V. J. Chem. Phys. 2005, 123, 144110. (17) Torrie, G. M.; Valleau, J. P. Chem. Phys. Lett. 1974, 28, 578. (18) Mehlig, B.; Heerman, D. W.; Forrest, B. M. Phys. ReV. B 1992, 45, 679. (19) Attard, P. J. Chem. Phys. 2002, 116, 9616. (20) Desgranges, C.; Delhommelle, J. J. Am. Chem. Soc. 2006, 128, 15104. (21) Desgranges, C.; Delhommelle, J. J. Chem. Phys. 2007, 126, 054501. (22) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: Oxford, 1987. (23) Trudu, F.; Donadio, D.; Parrinello, M. Phys. ReV. Lett. 2006, 97, 105701. (24) Leyssale, J. M.; Delhommelle, J.; Millot, C. Chem. Phys. Lett. 2003, 375, 612.

Desgranges and Delhommelle (25) Leyssale, J. M.; Delhommelle, J.; Millot, C. J. Am. Chem. Soc. 2004, 126, 122286. (26) Leyssale, J. M.; Delhommelle, J.; Millot, C. J. Chem. Phys. 2005, 122, 104510. (27) Leyssale, J. M.; Delhommelle, J.; Millot, C. J. Chem. Phys. 2005, 122, 184518. (28) Desgranges, C.; Delhommelle, J. J. Am. Chem. Soc. 2007, 129, 7012. (29) Desgranges, C.; Delhommelle, J. J. Phys. Chem. B 2007, 111, 1465. (30) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Phys. ReV. B 1983, 28, 784. (31) Bolhuis, P. G.; Dellago, C.; Chandler, D. Faraday Discuss. 1998, 110, 421. (32) Allen, R. J.; Warren, B.; ten Wolde Phys. ReV. Lett. 2005, 94, 018104. (33) Moroni, D.; ten Wolde, P. R.; Bolhuis, P. G. Phys. ReV. Lett. 2005, 94, 235703. (34) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (35) Robbins, M. O.; Kremer, K.; Grest, G. S. J. Chem. Phys. 1988, 88, 3286. (36) Hoy, R. S.; Robbins, M. O. Phys. ReV. E 2004, 69, 056103. (37) Hamaguchi, S.; Farouki, R. T.; Dubin, D. H. E. Phys. ReV. E 1997, 56, 4671. (38) Alexander, S.; McTague, J. P. Phys. ReV. Lett. 1978, 41, 702. (39) Oxtoby, D. W. Acc. Chem. Res. 1998, 31, 91. (40) Shen, Y. C.; Oxtoby, D. W. J. Chem. Phys. 1996, 105, 6517. (41) Shen, Y. C.; Oxtoby, D. W. Phys. ReV. Lett. 1996, 77, 3585.