Classical Nucleation Theory (CNT) has recently been used in conjunction with a seeding approach to simulate nucleation phenomena at small-to-moderate supersaturation conditions when large free-energy barriers ensue. In this study, the conventional seeding approach [J. R. Espinosa et al., J. Chem. Phys. 144, 034501 (2016)] is improved by a novel, more robust method to estimate nucleation barriers. Inspired by the interfacial pinning approach [U. R. Pedersen, J. Chem. Phys. 139, 104102 (2013)] used before to determine conditions where two phases coexist, the seed of the incipient phase is pinned to a preselected size to iteratively drive the system toward the conditions where the seed becomes a critical nucleus. The proposed technique is first validated by estimating the critical nucleation conditions for the disorder-to-order transition in hard spheres and then applied to simulate and characterize the highly non-trivial (prolate) morphology of the critical crystal nucleus in hard gyrobifastigia. A generalization of CNT is used to account for nucleus asphericity and predict nucleation free-energy barriers for gyrobifastigia. These predictions of nuclei shape and barriers are validated by independent umbrella sampling calculations.

First order phase transitions from metastable phases to stable phases often take place through a mechanism of nucleation and growth. In particular, crystal nucleation involves the formation of ordered domains of the stable phase (nuclei) that, once sufficiently large, can spontaneously grow to reach the macroscopic size. The condition at which a given nucleus has an equal probability to grow and decay is termed as critical.1,2 To a good approximation, this critical state corresponds to a saddle region in the transition pathway and, hence, where the nucleus attains its highest free energy. The nucleus must overcome this critical or transition state to form the stable phase, and the difference between the free energy of the critical state to the pure metastable state is called the nucleation barrier(ΔG*). Detailed knowledge of the mechanism of how such nuclei form and grow is of great practical importance, for this can open the door to rationally manipulate or bias these processes, e.g., to catalyze the assembly of nanoparticles into ordered structures.3,4

Computer simulations have emerged as a valuable complementary tool to experiments and theory to study the mechanism of nucleation and growth for a variety of phase transitions and systems. Since the formation of an ordered domain in a moderately metastable disordered phase is a rare event, using brute force simulations to sample representative nucleation trajectories is computationally impractical unless the nucleation free-energy barrier is small (say order of a few kBT). For example, for a nucleation barrier of 20 kBT, the chance of sampling the critical condition is nearly one in a billion trial configurations attempting to escape the metastable basin. Hence, a variety of rare-event sampling techniques have been developed and employed to understand nucleation kinetics, such as umbrella sampling,5–7 forward flux sampling,7–11 and metadynamics.12–14 However, application of these techniques is limited in cases where nucleation barriers are steep and/or high where the ratcheting of configurations through milestones along the pathway from the disordered phase becomes inefficient and tedious. This is because, when climbing steep free-energy profiles, the relative frequency with which the high free-energy states are visited is quite low in proportion to the low free-energy states, thus leading to poor sampling and inaccurate estimates of free energy differences. A way to overcome this limitation is the “seeding approach.”15–18 The approach essentially involves insertion of a preformed ordered domain (seed) into the disordered phase and observing its growth/decay behavior. Several simulations are then run for different degrees of supersaturation trying to narrow down the conditions when a transition occurs between growth and decay behaviors,17,18 i.e., the seed grows spontaneously above a critical supersaturation while decays spontaneously below it. Once the critical conditions are found, they are used in combination with classical nucleation theory (CNT) to evaluate the key parameters that characterize the nucleation process such as the nucleation free-energy barrier, interfacial energy, and the nucleation rate. This seeding approach has been successfully demonstrated for nucleation in L-J systems,15 NaCl,18 hard spheres,17 water, lattice models,19 and metals.20 Besides homogenous nucleation, the method has been employed to study heterogeneous nucleation phenomena in nanoparticles,21 including cases with a cubic seed.22 While seeding appears to be a promising method to study nucleation at low-to-moderate supersaturations (and to complement the scope of other methods that can resolve nucleation behavior for moderate-to-high supersaturations without resorting to CNT), its applicability relies on the following preconditions/assumptions:

  • CNT can be used to estimate the free-energy barrier height. Appropriate order parameters to detect the seed and classical (as opposed to non-classical) nucleation pathways are ingredients typically associated with this assumption.

  • The seed growth-to-decay transition behavior is abrupt and detectable with limited sampling.

  • The system can evolve from an initially arbitrary shaped seed (typically spherical) to the inherent seed shape in a timely fashion. The configurations obtained upon insertion and equilibration of the seed are a representative of those that would occur spontaneously in the system.

Assumption (i) underpins the seeding method and the extension to be elaborated in this work. Assumption (ii) can be invalid for systems with a flat or diffusive free energy profile around the top of the barrier (implying a low Zeldovich factor), where the transition from growth-to-decay behavior is gradual, and when even a slightly post-critical nucleus has a significant probability of dissolution. In such a scenario, it would be difficult to precisely pin point the critical condition altogether, and extensive sampling is required to quantify growth vs. decay behavior. We note here that conventional seeding studies have thus far only considered the “first occurrence” of growth to decay to bracket the critical condition.17,18 A more statistically rigorous implementation would entail, e.g., estimating growth/decay probabilities from more extensive sampling of trajectories, such that the critical condition2,16 can be established by the equality of growth and decay probabilities. Assumption (iii) could be violated if the initially assumed seed shape is very different from the inherent seed geometry.23–26 While equilibration of the seed can take place over a shorter time scale than that for seed growth/melting as implemented in conventional seeding,27 such a trend is system and seed-size dependent. Furthermore, it may not be straightforward to establish whether seed equilibration is attained, which may affect the seed growth/decay behavior and lead to unreliable results for both the seed size and critical condition. In fact, even though we can “surgically” insert an ordered seed into a disordered phase and then equilibrate17,18 it, the initial structures of the core and interface may remain far from their “natural” states over very long simulation periods. This is important because the configurations present at the interface can play a crucial role in the overall mechanism of nucleation and growth.3 

In this study, we propose and demonstrate a variant of the seeding approach that improves the handing of the constraints (ii) and (iii) listed above. In essence, we introduce an iterative scheme (similar to interfacial pinning28–30) that converges to the critical nucleation conditions for a seed of a prescribed size without the need to assume the seed shape or to constrain the interface morphology. We call this method nucleus-size pinning (NSP) which we couple to a generalization of CNT to extend its scope to treat seeds of arbitrary geometry. Our main application of this approach is for the nucleation of an aspherical seed for hard gyrobifastigia (GBFs). GBF is a Johnson solid consisting of four square and four triangular faces whose equilibrium phase behavior has been studied earlier, being remarkable for its ability to form a space-filling crystal and inability to form any mesophase.31 GBF can be seen as a test-bed that provides a stringent test to simulation studies of crystal nucleation, given not only the anisotropy of its shape and of its crystal lattice but also the large free-energy barriers that have been associated with their disorder-to-order phase transition.31 

The rest of this paper is organized as follows. In Sec. II A, we present the generalization of CNT for arbitrary seed geometries. In Sec. II B, we describe the details of the NSP method, followed by a description of simulation details and order parameters in Secs. II C and II D, respectively. In Sec. III A, we first validate the applicability of NSP by replicating established nucleation barrier results for hard spheres,6,7,17 and then in Sec. III B, we apply this method to GBF crystallization. Finally, in Sec. IV, we present our conclusions and an outlook for the method.

For a system with an ordered (solid) seed of volume Vs and surface area As, the extensive free energy of the system ϕT in reference to a disordered (liquid) bulk phase is given as

ϕT=ϕbulk+ϕi,
(1)

where ϕi is the contribution of the interface between the seed and the liquid, and ϕbulk is the difference between contributions of bulk solid and liquid phases for a volume Vs. ϕi is proportional to As, while the bulk contributions are proportional to Vs or to the number of particles in the seed Ns, thus Eq. (1) becomes

ϕT=ΔμNs+γAs,
(2)

where Δμ is the chemical potential difference between liquid and solid, and γ is the interfacial free energy to maintain the solid-liquid interface. For systems where the liquid phase is metastable, Δμ is positive (implying that ϕ decreases with Vs because the solid is thermodynamically more stable). For a given seed of fixed geometry, Ns is proportional to the volume Vs, which in turn is related to the length scale of the seed (L) as

NsVSL3
(3)
LNs13.
(4)

Similarly, for a given geometry, the surface area As is proportional to L2 so that, according to Eq. (4),

AsL2Ns23.
(5)

Equation (5) can be used to rewrite Eq. (2) with a proportionality constant (kA) as

ϕT=ΔμNs+kANs23.
(6)

In order to find the critical nucleus size Ns*, the expression for the derivative of ϕT with respect to Ns is set to zero,

dϕTdNs=Δμ+2kANs*133=0,
(7)
Ns*=2kA3Δμ3.
(8)

Which when substituted in Eq. (6) gives

ϕT*=Δμ2kA3Δμ3+kA2kA3Δμ2=4kA327Δμ2=Ns*Δμ2.

This critical free energy barrier will also be denoted as ΔG* so that this equation is

ϕT*ΔG*=Ns*Δμ/2.
(9)

This expression is the same as the one associated with CNT and has been reported in earlier studies.17,18,32,33 The key feature of Eq. (9) is its independence of seed geometry; as it only assumes the validity of the square-cube law relating the surface area and the volume. Thus, Eq. (9) can be used to calculate the nucleation barrier (ΔG*) once a critical seed (of size Ns*) of appropriate geometry is obtained, which need not be spherical.

In this section, we shall derive the iterative expression used by nucleus-size pinning for estimating the critical supersaturation conditions. For concreteness, this derivation is targeted to estimate the critical pressure (at constant temperature), but it can be readily generalized for any thermodynamic field other than pressure (e.g., temperature) as described in our earlier publications.28 For a pure component system at a fixed temperature T, the fundamental thermodynamic equation can be written as

dμ=vdp.
(10)

In an isothermal-isobaric (NPT) ensemble with N total particles at pressure p and temperature T, the probability p(Ns) of a state with Ns particles in the solid phase seed can be written as

pNsexp(ϕT(Ns)).
(11)

The total free energy can again be split up as

ϕTNs=ϕbulk,solidNsϕbulk,liquidNs+ϕiNs.
(12)

For critical nucleation conditions, the slope ST=dϕTdNs should be zero, thus

ST=dϕbulk,soliddNsdϕbulk,liquiddNs+dϕidNs=0.
(13)

But since at constant pressure and temperature ϕbulk,j=μjNs, Eq. (13) becomes

ST=μsolidμliquid+dϕidNs.
(14)

Let p* denote the critical pressure corresponding to a preset, target seed size Ns,t. Assume that an initial isothermal isobaric simulation is conducted at a pressure pA corresponding to a point A on the free energy curve where the seed size is Ns,A. The goal is to devise a scheme that, given simulation data at pA, it returns a guess pressure pB that is closer to p*, a process that can be iterated till convergence to p*, in a manner analogous to interfacial pinning.28–30 For this purpose, we approximate selected terms in Eq. (14) by truncated Taylor expansions around point A to obtain an extrapolation to a point of interest B as follows:

μj,B=μj,A+μj,ApΔpBA,
(15)
dϕi,BdNs=CA,
(16)

which are the first and zeroth order approximations, respectively, with μj,A and CA being the values of properties at point A, and ΔpBA=pBpA. Note that for Eq. (16), we assume that for small changes in pressure the interfacial contribution to the slope is not a strong function of pressure17 and hence can be approximated as a constant (whose value is CA). This assumption may be relaxed with the use of a suitable modification to CNT that describes the variation of interfacial free energy with pressure; a more detailed analysis of this assumption, including a CNT-based estimation of γ, is given in the  Appendix. Substituting Eq. (10) into Eq. (15), we get

μj,B=μj,A+v¯AΔpBA,
(17)

where v¯A is the ensemble average of the specific volume at point A. Furthermore,

Δμsl,B=μsolid,Bμliquid,B=Δμsl,A+Δv¯sl,AΔpBA,
(18)

where Δμsl,A and Δv¯sl,A are the differences in chemical potential and specific volume at point A, respectively. Substituting Eqs. (16) and (18) into Eq. (14), we obtain an expression for ST,B,

ST,B=Δμsl,A+CA+Δv¯sl,AΔpBA=ST,A+Δv¯sl,AΔpBA.
(19)

For a given point A, if ST,A and Δv¯sl,A are known, the needed change in pressure ΔpBA to reach the desired value of ST,B = 0 (critical condition) can be obtained as

ΔpBA=ST,AΔv¯sl,A
(20)
pB=pAST,AΔv¯sl,A.
(21)

This estimate of pB can be further improved by iterating the entire process, eventually converging toward p*, the critical pressure. Such a convergence is expected for systems having a monotonically varying slope for free energy (with respect to the nucleus size), as was the case for the systems studied in this work and is likely the case for many nucleation and growth phenomena.

At any given pressure, the value of Δv¯sl can be obtained using the equations of state for either phases. The slope ST of the free energy profile can be obtained using an approach similar to interfacial pinning.28–30 However, instead of pinning the position of the interface, we pin the nucleus size by biasing the system with a harmonic potential of the form

ΔUH=κ2NsNs,t2.
(22)

We note that while this biasing potential is similar to the one traditionally used in umbrella sampling simulations, the purpose here is not to map the local free energy profile but rather to provide the balancing potential to pin the seed size and iteratively converge to the underlying free energy maxima. The ideal probability distribution associated with such a potential is a Gaussian distribution centered at Ns = Ns,t, noting that both Ns,t and the pinning force constant κ are user defined values. When applied to a given point A on the total free energy curve, the resulting probability distribution is given by

pNsexp((ϕT+ΔUH)).
(23)

In the neighborhood of Ns,t, ϕT can be assumed to be a linear function with slope ST. With this assumption, it can be shown that the p(Ns) is also a quasi-Gaussian distribution with a shifted mean value of Ns such that

ST=κ(NsNs,t).
(24)

Thus, by calculating the average value of Ns in a biased simulation [under the external potential of Eq. (22)], we can determine the value ST needed for the iterative scheme of Eq. (21) to obtain p*. Convergence proceeds as the pressure correction becomes negligible and ST approaches zero.

1. Model

For both hard spheres and gyrobifastigia, we used the hard pair potential given by

βUij=0ifnooverlapifoverlap.
(25)

For spheres, overlaps are detected by evaluating if the center-to-center distance rij is less than the diameter (σ). For gyrobifastigia, overlaps are checked by employing the separating axis theorem.34 

2. Metropolis Monte Carlo

Metropolis Monte Carlo (MC) simulations were performed in an isothermal-isobaric (NPT) ensemble where the number of particles, the pressure, and the temperature of the system were kept constant. As is customary,6,7 for spheres, the pressure was scaled as p=βpaσ3, where β=1kBT, pa is the unscaled pressure, kB is the Boltzmann constant, and σ is the diameter of the sphere. For GBF, the pressure was scaled as p=βpaac3, where ac is the radius of the circumscribing sphere. Simulations were conducted using periodic boundary conditions and cycles consisting of N translational moves, 2 isotropic volume moves, and for GBF, additional N rotation moves. A volume move attempt involves the rescaling of the simulation box size while maintaining constant the reduced coordinates of the particles’ centers of mass (and particle orientations for GBF).31 For spheres, the simulation box was cubic, while for GBF, the box was cuboidal in shape to accommodate equal number of particle layers of the anisotropic solid-phase lattice vector. For hard spheres, the degree of supersaturation was taken and interpolated from the literature.17 For GBF, the equations of state for the various phases were obtained and selected degrees of supersaturation (DSS) (= βΔμ) were evaluated upon integration of suitable branches of those equations of state, by using the same formulas discussed in earlier work.3,31

3. Nucleus-size pinning

The biased NPT MC simulations were performed in a manner similar to those described in prior studies.28 Following the conventional seeding16,17 methods, an ordered spherical seed of a chosen size is inserted into a hole carved in a disordered phase by eliminating any overlapping particles. To avoid finite size effects, a system sufficiently larger than the seed size was simulated to ensure that the seed does not interact with its periodic images. For this purpose, the size and shape of the simulation box were altered as necessary. Both the ordered seed and the disordered configurations were at equilibrated volume fractions corresponding to the pressure of interest as per the equations of state. To save computational time, the order parameter calculations were performed once every 100 MC cycles, upon which the modified Metropolis criteria with the biased potential were used to accept or reject the trajectory of the last 100 cycles (if rejected, the old trajectory was restored). Once key system properties such as the cluster size, the specific volume, and, in the case of GBF, the aspect ratio of the nucleus, have fully converged (usually in about 3 × 106 MC cycles), the run is stopped and statistics for the cluster size are evaluated (note that this equilibration process eliminates potential artifacts from the original seed preparation). Using Eq. (21), a new estimate of the sought-after pressure is obtained, and a new iteration of the process is started.

For the method to work efficiently, the strength of harmonic potential (κ) needs to be chosen carefully. If too high, the system would be too constrained to be able to explore diverse configurations, and statistics thus obtained would not be representative. If too low, the assumptions regarding the local behavior of various system properties as implied by the Taylor series expansions in Eqs. (15) and (16) may not be correct. An insufficiently strong potential can lead to multimodal p(Ns) distribution, with peaks on either side of the barrier. This is possible when the unbiased free energy profile exhibits a high variation in slope near the target nucleus size, allowing for solutions to Eq. (13) on either side of the peak. This issue can be remedied by increasing κ and a suitable value can be found using the following criterion:28 

κ>1σNs2,
(26)

where σNs2 is the variance in the Ns values for a short unbiased simulation around the target conditions. For the purpose of the estimation of barrier heights, since the method uses CNT, it is more indicated for lower DSS where larger nuclei ensue. Thus, it is a complementary technique to conventional methods such as umbrella sampling and forward flux sampling which tend to be inefficient at low DSS.

4. Umbrella sampling (US)

To provide an independent estimation of the homogeneous solid-nucleation free-energies for GBFs for comparing estimates from the proposed NSP method, we also implemented umbrella sampling (US) simulations.3 Using as a reaction coordinate the size of the greatest cluster of ordered particles (Ns), the transition path from the disordered to the ordered state is divided into overlapping windows. Each window is simulated separately with reflective walls and Ns is recorded every 2 MC cycles. Reflective walls are implemented such that any trajectory moving outside the window at the end of 2 MC cycles is returned to the configuration prior to the 2 MC cycles, which is counted again. Statistics obtained from each window are used to obtain relative free energies for the Ns states within a window. Finally, individual sections are stitched together by matching values at the boundaries of the windows, keeping the value for the most frequent entry near Ns = 0 as the reference. All US calculations were performed with a system of 1728 particles.

1. Hard spheres

As in previous crystal nucleation simulations of hard spheres,6,7 we use the size of the largest translationally ordered cluster as the reaction coordinate. We use the q6 Steinhardt35 order parameter which is defined as follows: For every particle i, the local bond order parameter, ql,m(i), is defined by

ql,mi=1Nbij=1lYl,m(θi,j,ϕi,j),
(27)

where Nb(i) is the number of neighbors of particle i, Yl,m(θ,ϕ) are the spherical harmonics, θi,j and ϕi,j are polar and azimuthal angles between particle i and its neighbor j, respectively, l is the symmetry index, and the value of m ranges from −l to l. In this work, we use l = 6. The neighbors of particle i are those particles which are within the cutoff distance rc = 1.4σ of particle i. The translational-order correlation between particle i and its neighbor j, d6(i,j), is given by

d6i,j=m=66q6,miq6,m*jk=66q6,ki212l=66q6,lj212,
(28)

where the asterisk (*) denotes the complex conjugate. Two particles i and j are defined as translationally connected if d6(i,j) > 0.7. A particle with at least 7 translational connections is classified as translationally ordered or solid-like. Solid clusters are identified by the condition that any two solid-like particles within the rc belong to the same cluster. The tunable parameters in such order parameters need to be optimized using criteria described in prior studies.3,7

2. Gyrobifastigia

We found that the Steinhardt order parameters were not effective to capture the local translational symmetry present in the ABCD cubic lattice formed by GBFs.31 Hence, we devised an orientational order parameter that proved to be suitable to capture the inherent crystal symmetry. It is based on the observation that in the crystal lattice all constituent particles have their long-axes aligned. We can then define an orientational order correlation or a “dot-product” similar to d6 as follows. For a given pair of neighboring particles i and j (i.e., within the cutoff distance rc = 1.4a, where a is the side length of GBF), we define the orientational order correlation dz(i,j) as

dz(i,j)=zi.zj,
(29)

where z(i) and z(j) are unit vectors pointing along the long-axis of particles i and j. dz(i,j) is essentially the magnitude of the cosine of the angle included between the two axis vectors and hence varies from zero to 1 for orthogonal to aligned configurations, respectively. Note that the correlation is unaffected by whether the vectors are flipped 180° and the operation is commutative. Particles i and j are considered orientationally correlated if dz(i,j) > dz,min = 0.95, implying a tolerance of ∼18° from a perfect alignment. If a given particle has at least 4 aligned neighbors, it is considered solid-like. If two solid-like particles are orientationally correlated, then they belong to the same cluster. The size of the largest orientationally ordered cluster is used as the reaction coordinate for the nucleation process.

To validate the NSP method for a well-studied system, we apply it to hard spheres (HS). The free energy barriers ΔG* for disorder to order transition for HS have been determined via a variety of methods such as US6 and seeding.17 Thus, our first objective was to test the barriers and critical nuclei sizes found by NSP by comparing with those reported in the literature. The results for various target nuclei sizes are tabulated in Table I and shown in Fig. 1(a). It can be seen that the NSP-estimated values provide a good interpolation to the reported estimates, hence showing a suitable consistency with prior data. It typically takes 5-6 iterations to obtain convergence within one percent relative error in the pressure value.

TABLE I.

Nucleus-size pinning for hard spheres. N is the total number of particles in the system, Nt is the target seed size, and pi is the initial guess of the critical pressure (whose converged value is p*) and κ is the strength of the biasing potential as described in Eq. (22).

NNtpiκp*
55 183 400 14 0.005 14.42 
55 183 600 14 0.0005 14.066 
55 143 800 14 0.0005 13.783 
55 116 1300 12.5 0.0005 13.37 
55 116 1300 14 0.0005 13.42 
NNtpiκp*
55 183 400 14 0.005 14.42 
55 183 600 14 0.0005 14.066 
55 143 800 14 0.0005 13.783 
55 116 1300 12.5 0.0005 13.37 
55 116 1300 14 0.0005 13.42 
FIG. 1.

Nucleus-size pinning simulation results for hard spheres. (a) Critical nucleation pressures (p*) for various critical nucleus sizes Ns* as determined by NSP (blue) and by the seeding method of Espinosa et al.17 (red). (b) Convergence toward the critical pressure for a target critical nucleus size of 1300 for different initial pressures of 12.5 (blue) and 14 (red). The values asymptotically converge toward a common critical pressure value, with the last step correction on each run being less than Δp*0.02.

FIG. 1.

Nucleus-size pinning simulation results for hard spheres. (a) Critical nucleation pressures (p*) for various critical nucleus sizes Ns* as determined by NSP (blue) and by the seeding method of Espinosa et al.17 (red). (b) Convergence toward the critical pressure for a target critical nucleus size of 1300 for different initial pressures of 12.5 (blue) and 14 (red). The values asymptotically converge toward a common critical pressure value, with the last step correction on each run being less than Δp*0.02.

Close modal

For Nt = 1300, the convergence of the algorithm was tested with different initial guesses for pressure (one higher and one lower than the reported value), as shown in Fig. 1(b). The values obtained from either direction are asymptotically convergent. In practice, one may stop iterations as successive values are within the desired tolerance, and yet the estimates from the different starting points may differ by a non-negligible amount. This is possible if the profile is relatively flat (i.e., the corresponding Zeldovich factor is low). Convergence to a narrower range of p* could be improved by using a weaker biasing potential that can better resolve weaker slopes near the top of the barrier.

The associated free energy barriers can be calculated using Eq. (9), which shall be in agreement with reported results for seeding.17 These estimates will become increasingly accurate for lower DSS as the critical nucleus becomes larger and the system approaches the macroscopic phase properties as assumed in the CNT approach.

1. Umbrella sampling

We performed US calculations to generate reference results for the crystal critical nucleus size (Ns*) and free-energy barrier (βΔG*) of GBF to test the predictions of NSP (for N*) and its combination with the CNT (for βΔG*). These simulations generate the complete free energy profile along the reaction coordinate, from which one can readily obtain Ns* and βΔG*.

The results from US are tabulated in Table II. Figure 2(a) shows the estimated free energy profiles for GBFs at various pressures. Figure 2(b) shows a comparison of nucleation barriers for GBFs and hard spheres.7 It is observed that for the same DSS value, GBF has a much higher βΔG* than hard spheres. These large ΔG* for GBF may be due to the tendency of neighboring particles to pack in configurations that, while minimizing excluded volume36 and maximizing packing entropy,37 are inconsistent with the final ABCD crystal structure of GBFs. For example, two particles may tend to match their square facets while maintaining their long axes at a 60° angle [Fig. 2(d)], excluding a volume (to other particles) larger than that in the ABCD lattice where the particle long-axes align [Fig. 2(c)]. Such misaligned configurations can lead to “dead ends” in the nucleation pathway and require particles to first disassemble to further progress toward the crystalline global order.4 Such poisoning can even preclude nucleation as is the case in the isotropic-nematic transition of colloidal rods;38 we speculate that the coupling of orientational and translational order in GBFs helps the system overcome similar local configurational traps. Particle-alignment defects near the crystal interface likely increase the interfacial tension; in fact, our estimates shown in Fig. 8 in the  Appendix indicate that γ of GBFs is nearly twice larger than that of hard spheres. In addition, since for a given size, an aspherical seed has an associated larger surface area, a larger effective interfacial contribution to the free energy and a higher barrier height would also be expected. Since estimating larger βΔG* via US or another standard simulation method becomes more challenging, its calculation for GBFs for DSS = βΔμ < 0.5 would be an especially suitable application for the NSP method.

TABLE II.

Results for barrier height (βΔG*) and critical nucleus size (Ns*) from umbrella sampling (US) calculations for gyrobifastigia for various pressures (p), along with their corresponding degrees of supersaturation (βΔμ).

Reduced pressure (p=βpavp)DSS (βΔμ)βΔG*Ns*
12.5 0.684 73.0 252 
13 0.873 46.1 153 
14 1.244 19.0 48 
15 1.608 8.6 27 
Reduced pressure (p=βpavp)DSS (βΔμ)βΔG*Ns*
12.5 0.684 73.0 252 
13 0.873 46.1 153 
14 1.244 19.0 48 
15 1.608 8.6 27 
FIG. 2.

Nucleation free-energy barriers for gyrobifastigia using umbrella sampling. (a) Representative free energy (βΔG) vs. nucleus size (Ns) (for different pressures (p*). (b) Nucleation barriers βΔG* for gyrobifastigia (diamonds) and for hard spheres7 (circles) at various DSSs. Configurations of gyrobifastigia particle pairs in (c) the ABCD crystal and (d) a non-crystalline high local packing entropy state.4,37

FIG. 2.

Nucleation free-energy barriers for gyrobifastigia using umbrella sampling. (a) Representative free energy (βΔG) vs. nucleus size (Ns) (for different pressures (p*). (b) Nucleation barriers βΔG* for gyrobifastigia (diamonds) and for hard spheres7 (circles) at various DSSs. Configurations of gyrobifastigia particle pairs in (c) the ABCD crystal and (d) a non-crystalline high local packing entropy state.4,37

Close modal

Interestingly, the morphology of the nuclei formed for GBF is significantly aspherical as shown in Fig. 3(a). This behavior likely stems from the difference in interfacial energies of the nucleus along directions parallel to the particles’ long axes and perpendicular to it. A feature emerging out of the geometry of the GBF lattice is that the facets of the anisotropic nucleus are slanted at an angle of ≈ 73.9° as revealed in Fig. 3(a), inset.

FIG. 3.

Morphology of nuclei for gyrobifastigia. (a) Near-critical nucleus (Ns*150) obtained by umbrella sampling at p* = 13, evidencing its aspherical, elongated shape along the particle long-axis. Inset: The nucleus is approximately faceted at a specific angle = tan−1(23). In the inset, a is the side length of a gyrobifastigium. (b) A quasi-octahedron shaped (post-critical) nucleus corresponding to Ns ≈ 8000. Particles are colored according to the deviation from the director of the nucleus—from blue to red as the alignment decreases.

FIG. 3.

Morphology of nuclei for gyrobifastigia. (a) Near-critical nucleus (Ns*150) obtained by umbrella sampling at p* = 13, evidencing its aspherical, elongated shape along the particle long-axis. Inset: The nucleus is approximately faceted at a specific angle = tan−1(23). In the inset, a is the side length of a gyrobifastigium. (b) A quasi-octahedron shaped (post-critical) nucleus corresponding to Ns ≈ 8000. Particles are colored according to the deviation from the director of the nucleus—from blue to red as the alignment decreases.

Close modal

We also investigated how the shape anisotropy of the crystal nucleus evolved as it grew beyond its critical size. While the nucleus aspect ratio could increase indefinitely with size to become filamentous in shape, it was conjectured that it maintained a quasi-octahedral shape, based on the characteristic slant of the nucleus facets. To test this conjecture, we simulated the growth of a post critical nucleus up to a size of Ns ≈ 8000 (over 50 times its critical size) at which point local features are smoothed out and macroscopic behavior could be considered emergent. We observed that the shape anisotropy persists even for such large nuclei and the aspect ratio remains about the same [Fig. 3(b)]. The presence of a nearly invariant nucleus geometry validates the use of the CNT equation (9) to estimate ΔG*. We thus estimate ΔG* from the known critical nucleus sizes (from US) using Eq. (9); the results and their error bars are shown in Fig. 4. We find that CNT tends to overestimate ΔG* by a relatively constant deviation of ∼10kBT, but the relative deviations decrease as DSS is reduced. These trends are expected, as the nuclei tend to have a larger size for lower DSS and hence approach the macroscopic limiting behavior that CNT assumes. Note that one could use a ΔG* value obtained from US, ideally corresponding to large Ns* (and small Δμ), to fine-tune the order parameter definition (i.e., parameters rc and dz,min in Sec. II D 2) until the Ns* calculated in the US simulation agrees with the Ns* calculated via the CNT equation (9). While in this way the US and CNT curves in Fig. 4 would be brought in closer agreement, this would not necessarily provide an independent test of the quality of the order parameter or of the general applicability of CNT.7 Regardless of such refinements, ΔG* estimates from the critical nucleus size [via Eq. (9)] are expected to have large uncertainties for smaller nucleus sizes (when errors in estimates of interfacial particles can contribute significantly to the nucleus size) and become increasingly accurate as lower DSS values (and hence larger cluster sizes) are probed.

FIG. 4.

Comparison of CNT-estimated nucleation barriers (βΔG*) with those calculated using umbrella sampling (US) as a function of degree of supersaturation (DSS). The CNT values are obtained using Eq. (9) by inputting Ns* values from umbrella sampling.

FIG. 4.

Comparison of CNT-estimated nucleation barriers (βΔG*) with those calculated using umbrella sampling (US) as a function of degree of supersaturation (DSS). The CNT values are obtained using Eq. (9) by inputting Ns* values from umbrella sampling.

Close modal

2. Nucleus-size pinning for gyrobifastigia

Having established the presence of anisotropic solid nuclei and the applicability of the CNT to GBF in Sec. III B 1, we now demonstrate how the NSP method can be implemented as an alternative and complement to US to study the shape and size of the critical nuclei and to estimate nucleation barriers. The simulations were performed for a range of nucleus sizes that overlaps and extends beyond the upper bound explored before using US. The results are shown in Table III. As shown in Fig. 5(a), the relationship between critical nucleation pressure and Ns* as determined by NSP agrees well with that obtained from our US calculations performed earlier. Similarly, Fig. 5(b) shows that NSP with CNT produces free-energy barrier heights which slightly overestimate those found by US, as was earlier observed (in Fig. 4).

TABLE III.

Results from nucleus-size pinning method for GBF.

piNNtκp*DSS = βΔμβΔG*
14 3373 100 0.01 13.34 0.977 48.875 
11.5 1666 100 0.05 13.35 0.9813 49.065 
14 1664 50 0.05 13.9004 1.1844 29.61 
14 3373 150 0.01 12.9083 0.8155 61.163 
14 3344 300 0.01 12.41 0.6275 94.125 
14 3334 500 0.005 12.1205 0.5149 128.725 
piNNtκp*DSS = βΔμβΔG*
14 3373 100 0.01 13.34 0.977 48.875 
11.5 1666 100 0.05 13.35 0.9813 49.065 
14 1664 50 0.05 13.9004 1.1844 29.61 
14 3373 150 0.01 12.9083 0.8155 61.163 
14 3344 300 0.01 12.41 0.6275 94.125 
14 3334 500 0.005 12.1205 0.5149 128.725 
FIG. 5.

Comparison of NSP (blue) and umbrella sampling (US, red) for GBF. (a) Comparison of critical conditions: p* and Ns* are the critical pressure and nucleus size, respectively. (b) Comparison of free energy barrier height values (βΔG*) as a function of degree of supersaturation (βΔμ).

FIG. 5.

Comparison of NSP (blue) and umbrella sampling (US, red) for GBF. (a) Comparison of critical conditions: p* and Ns* are the critical pressure and nucleus size, respectively. (b) Comparison of free energy barrier height values (βΔG*) as a function of degree of supersaturation (βΔμ).

Close modal

The converged values obtained from iterations started from different initial pressure guesses were quite close in the cases examined, presumably due to a sharper peak of the free energy profile as compared to the calculations for hard spheres discussed in Sec. III A. Since there is no constraint on the shape of the nucleus but only on its size, its shape should eventually tend to the “natural” one. Indeed, Fig. 6 shows that the shape of the GBF nuclei as determined by NSP has geometries consistent with those found earlier via US, becoming increasingly quasi-octahedron-like as the size increases and having an aspect ratio of ∼2 found from the principal moments of the nucleus inertial tensor, as a ratio of the longest axis to the mean of the shorter ones. This aspect ratio would slightly increase for much larger nucleus sizes [it is still 2.1 for the one in Fig. 3(b)] if the aberrations at the top and bottom vertices of the quasi-octahedron were to smooth out.

FIG. 6.

Various critical nuclei of GBF as obtained using the NSP method for Ns ≈ (a) 50, (b) 100, (c) 150, (d) 300, and (e) 500, with aspect ratios close to 2.0 (1.9, 2.0, 1.9, 2.0, and 2.1, respectively). Particles are colored as in Fig. 3.

FIG. 6.

Various critical nuclei of GBF as obtained using the NSP method for Ns ≈ (a) 50, (b) 100, (c) 150, (d) 300, and (e) 500, with aspect ratios close to 2.0 (1.9, 2.0, 1.9, 2.0, and 2.1, respectively). Particles are colored as in Fig. 3.

Close modal

In this study, we introduced the nucleus-size pinning (NSP) method which, as an extension of the seeding technique,17 is designed to iteratively converge to the critical state in a nucleation process. We first validated the method by reproducing the known results for hard spheres. We then demonstrated the robustness of the method by identifying non-trivial anisotropic growth of nuclei during disorder-to-order transition in gyrobifastigia. By combining NSP with a formulation of classical nucleation theory that allows the treatment of an aspherical nucleus, we were able to determine nucleation barriers for gyrobifastigia in a range of small DSS where conventional methods would be computationally very expensive. The method allows the “natural” seed geometry to be attained and the critical state to be readily identified, as demonstrated here for gyrobifastigium. The proposed method may even prove useful in studying the transition state in situations where a global order parameter rather than a local order parameter (as in classical nucleation) is used to describe the phase transition process. The proposed NSP approach can be easily applied to other systems for the determination of the critical states provided suitable order parameters are employed. For example, applications should be readily implementable to nucleation processes where the degree of metastability (or DSS) is defined by departures of temperature (as opposed to pressure) from the coexistence conditions. Moreover, NSP can help reveal important microstructural features of transition states, which might be non-trivial to hypothesize a priori. Altogether, NSP provides a means to effectively and unambiguously target the simulation of the critical transition state in a nucleation process, one that complements existing rare event sampling techniques. The iterative process can be further automated such that the equilibration (at each iteration step) is checked via statistical tests for normality39 and by pre-specifying a tolerance for global convergence.

While the NSP method provides a more robust way for establishing the conditions for cluster criticality than the traditional seeding approach, it retains certain limitations. First, it requires a valid model (such as CNT) to relate the critical nucleus size to the nucleation free-energy barrier. Only in the presence of an applicable model can we determine other properties of interest such as interfacial energy. Nevertheless, even in the absence of such a model, the method is still able to describe the critical nuclei at different DSS, whose characterization may be of interest in various contexts, e.g., to identify differences in structural order between the core and interfacial regions of potential configurational bottlenecks associated with transition states.40–42 Second, the convergence of the algorithm is dependent upon the “flatness” of the target free energy profile. If there is a broad region (in comparison to the spread of the biasing harmonic potential) with nearly constant free energy around the barrier top, then the method can end up identifying any of those states as the critical one (within a user-specified tolerance). Note, however, that such flat barrier tops would pose challenges to many simulation methods and entail an intrinsic ambiguity to the definition of the critical nucleus. Third, and in common with other methods (like umbrella sampling), NSP requires the use of order parameters that are a good approximation to the reaction coordinate that best describes the transition process.43 In fact, the choice of nucleation order parameter (e.g., the criteria to decide which particles are solid-like) affects not only the effective height of the free-energy barrier but also its flatness (as pertaining to the second limitation listed before). Finally, variants of the NSP method can be implemented to improve the convergence of the iterative scheme, e.g., to more explicitly account for the dependence of interfacial free energy on supersaturation conditions (as discussed in the  Appendix).

Funding support from NSF Award No. CBET-1402117 is gratefully acknowledged. The authors are grateful to Prajwal B. Prakash and Unmukt Gupta for valuable discussions and sharing the base simulation code.

In Sec. II B, we assumed that the interfacial-energy contribution to the slope CA in Eq. (16) can be neglected for small changes in pressure Δp*<1, an assumption which is consequential to the convergence properties of the NSP iterations to the critical conditions but not to the final results. We note that this assumption is consistent with earlier studies that found the variation of interfacial energy (γ) to be small (i.e., less than 5%) for changes Δp*<117 and that in NSP iterations the usual step corrections in pressure are much smaller than Δp*<0.5. Nevertheless, in this section, we re-examine this assumption by first rearranging Eq. (14) to obtain the interfacial contribution to the slope (Si) as

Si=dϕidNs=ST(μsolidμliquid)=STΔμsl.
(A1)

Thus, for a simulation performed at a given pressure, Eq. (A1) can be used to evaluate Si, the interfacial contribution to the slope (note that this involves a contribution from interfacial tension and from nucleus geometry, either of which may vary with pressure). One such evaluation is illustrated in Table IV and Fig. 7 for a series of iterations in the case of hard spheres for a target cluster size of Nt = 1300.

TABLE IV.

Calculation of interfacial contribution to the slope from data of total slope and chemical potential difference for various iterations of NSP for hard spheres with Nt = 1300.

IterationpiΔμslSTSi
14.00 −0.235 −0.054 0.180 
13.72 −0.207 −0.029 0.178 
13.58 −0.193 −0.017 0.176 
13.49 −0.184 −0.005 0.179 
13.46 −0.181 −0.003 0.178 
13.45 −0.180 −0.003 0.178 
13.44 −0.179 −0.003 0.176 
13.42 −0.177 0.001 0.179 
IterationpiΔμslSTSi
14.00 −0.235 −0.054 0.180 
13.72 −0.207 −0.029 0.178 
13.58 −0.193 −0.017 0.176 
13.49 −0.184 −0.005 0.179 
13.46 −0.181 −0.003 0.178 
13.45 −0.180 −0.003 0.178 
13.44 −0.179 −0.003 0.176 
13.42 −0.177 0.001 0.179 
FIG. 7.

Key contributions to the variation of the slope ST of the free energy as a function of pressure. The interfacial contribution Si is nearly constant.

FIG. 7.

Key contributions to the variation of the slope ST of the free energy as a function of pressure. The interfacial contribution Si is nearly constant.

Close modal

It can be seen in Fig. 7 that the major contribution to the variation of the total slope (ST) is through the variation of Δμsl=DSS, while the interfacial contribution remains nearly constant ∼0.178, consistent with the assumption of Eq. (16).

Alternative implementations of the NSP could explicitly take into account the variation of γ with pressure. For instance, one can use the known values of Si=dϕidNs at two or more nearby pressures (from the iterative steps) to obtain an extrapolation model Si = Si(p*), and use such a model to improve the estimate of p* for the next iteration.

As a related extension, we note that the NSP method and CNT can also be combined to extract the interfacial tension γ from Si. Assuming for concreteness that the nucleus is spherical, then kA=γ36πNtρs213 and Si=2kA3N13, from which it follows that

γ=Si3Ntρs232π13.
(A2)

Equation (A2) can be used to estimate γ at successive steps in the NSP iteration process, as shown in Fig. 8(a) for hard spheres. The results are in line with what has been reported in prior studies using other techniques.6,17,44,45 The intermediate calculations, however, may be less reliable than that for the converged critical condition as a CNT fit may be more suitable near the top of the free-energy barrier than away from it.6 Similar calculations for γ could be performed for aspherical nuclei, perusing explicit expressions for volume and surface area of a nucleus of known geometry. For example, for GBF, if we approximate the nuclei as a prolate spheroid with an aspect ratio ∼2 (see Fig. 6), the interfacial tension can be written as

γ=Si4π33+13Ntρs28π13.
(A3)

In the reduced units employed, γ values for GBFs [Fig. 8(b)] appear to be larger than those for hard-spheres [Fig. 8(a)] for comparable DSS, a difference that correlates with the difference in free energy barrier heights shown earlier in Fig. 2(b). The intermediate-step γ estimates using this approach are also shown in Fig. 8(b). Note, however, that isolating an average, crystal-facet independent γ for the GBF nucleus should be seen rather as an approximate calculation, especially for small nuclei (large p*) where their shape may depart vastly from a spheroid, possibly causing the non-monotonous trend for p* > 13. Furthermore, the error bars were found to be rather large, i.e., ∼±0.05 for the initial pressure p* = 14, which should be representative of those of nearby points where equilibration was less converged. Within error bars then, the data available does not allow us to ascertain whether for p* > 13 the γ trend is non-monotonous or simply approaches a plateau.

FIG. 8.

Interfacial tension (γ) vs. pressures (p*) from various simulation methods. (a) Hard spheres, comparing NSP + CNT results for (converged) final pressure values and intermediate step values (excluding the first iteration), seeding method,17 US+CNT,6 cleaving walls (CW) approach,45 and capillary fluctuation (CF) method.44 (b) GBFs for NSP+CNT results.

FIG. 8.

Interfacial tension (γ) vs. pressures (p*) from various simulation methods. (a) Hard spheres, comparing NSP + CNT results for (converged) final pressure values and intermediate step values (excluding the first iteration), seeding method,17 US+CNT,6 cleaving walls (CW) approach,45 and capillary fluctuation (CF) method.44 (b) GBFs for NSP+CNT results.

Close modal

Finally, once we have an estimate for Si, and hence kA, we can estimate Ns* for the corresponding pressure using Eq. (8). Note that we can apply such a CNT-based approach to the estimate Ns* in a single simulation, irrespective of whether convergence to the critical condition has been reached. We tested this idea in Fig. 9 for hard spheres and GBF, by comparing the NSP+CNT based estimates of Ns* from intermediate steps and from the converged values and with data from other methods. We observe a good agreement in the results among the various techniques, suggesting the validity of CNT for the systems under study.

FIG. 9.

Comparison of estimates of critical conditions in the pressure (p*) vs. nucleus size (Ns*) plane, obtained from seeding,17 US = umbrella sampling, the NSP described in the main text, and the NSP+CNT approach described in the  Appendix. (a) Hard spheres. (b) Hard gyrobifastigia.

FIG. 9.

Comparison of estimates of critical conditions in the pressure (p*) vs. nucleus size (Ns*) plane, obtained from seeding,17 US = umbrella sampling, the NSP described in the main text, and the NSP+CNT approach described in the  Appendix. (a) Hard spheres. (b) Hard gyrobifastigia.

Close modal
1.
G. M.
White
,
J. Chem. Phys.
50
,
4672
(
1969
).
2.
D.
Kashchiev
,
J. Chem. Phys.
127
,
064505
(
2007
).
3.
V.
Thapar
and
F. A.
Escobedo
,
Phys. Rev. Lett.
112
,
048301
(
2014
).
4.
A.
Sharma
,
V.
Thapar
, and
F.
Escobedo
,
Soft Matter
14
,
1996
(
2018
).
5.
S.
Auer
and
D.
Frenkel
,
Nature
409
,
1020
(
2001
).
6.
S.
Auer
and
D.
Frenkel
,
J. Chem. Phys.
120
,
3015
(
2004
).
7.
L.
Filion
,
M.
Hermes
,
R.
Ni
, and
M.
Dijkstra
,
J. Chem. Phys.
133
,
244115
(
2010
).
8.
R. J.
Allen
,
C.
Valeriani
, and
P. R.
ten Wolde
,
J. Phys.: Condens. Matter
21
,
463102
(
2009
).
9.
V.
Thapar
and
F. A.
Escobedo
,
J. Chem. Phys.
143
,
244113
(
2015
).
10.
R. J.
Allen
,
D.
Frenkel
, and
P. R.
ten Wolde
,
J. Chem. Phys.
124
,
194111
(
2006
).
11.
R. J.
Allen
,
P. B.
Warren
, and
P. R.
ten Wolde
,
Phys. Rev. Lett.
94
,
018104
(
2005
).
12.
F.
Giberti
,
M.
Salvalaglio
, and
M.
Parrinello
,
IUCrJ
2
,
256
(
2015
).
13.
A.
Laio
and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
99
,
12562
(
2002
).
14.
F.
Trudu
,
D.
Donadio
, and
M.
Parrinello
,
Phys. Rev. Lett.
97
,
105701
(
2006
).
15.
X. M.
Bai
and
M.
Li
,
J. Chem. Phys.
123
,
151102
(
2005
).
16.
J. H.
ter Horst
and
D.
Kashchiev
,
J. Chem. Phys.
119
,
2241
(
2003
).
17.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
,
J. Chem. Phys.
144
,
034501
(
2016
).
18.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
,
J. Chem. Phys.
142
,
194709
(
2015
).
19.
Y.
Lifanov
,
B.
Vorselaars
, and
D.
Quigley
,
J. Chem. Phys.
145
,
211912
(
2016
).
20.
Y.
Sun
,
H.
Song
,
F.
Zhang
,
L.
Yang
,
Z.
Ye
,
M. I.
Mendelev
,
C.-Z.
Wang
, and
K.-M.
Ho
,
Phys. Rev. Lett.
120
,
085703
(
2018
).
21.
D.
Suh
and
K.
Yasuoka
,
J. Phys. Chem. B
115
,
10631
(
2011
).
22.
D.
Suh
and
K.
Yasuoka
,
J. Phys. Chem. B
116
,
14637
(
2012
).
23.
U.
Gasser
,
E. R.
Weeks
,
A.
Schofield
,
P. N.
Pusey
, and
D. A.
Weitz
,
Science
292
,
258
(
2001
).
24.
U.
Gasser
,
J. Phys.: Condens. Matter
21
,
203101
(
2009
).
25.
B.
O’Malley
and
I.
Snook
,
Phys. Rev. Lett.
90
,
085702
(
2003
).
26.
P.
Koß
,
A.
Statt
,
P.
Virnau
, and
K.
Binder
,
Phys. Rev. E
96
,
042609
(
2017
).
27.
A.
Zaragoza
,
M. M.
Conde
,
J. R.
Espinosa
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
,
J. Chem. Phys.
143
,
134504
(
2015
).
28.
V.
Thapar
and
F. A.
Escobedo
,
J. Chem. Phys.
141
,
124117
(
2014
).
29.
U. R.
Pedersen
,
J. Chem. Phys.
139
,
104102
(
2013
).
30.
U. R.
Pedersen
,
F.
Hummel
,
G.
Kresse
,
G.
Kahl
, and
C.
Dellago
,
Phys. Rev. B
88
,
094101
(
2013
).
31.
U.
Agarwal
and
F. A.
Escobedo
,
Nat. Mater.
10
,
230
(
2011
).
32.
N. E. R.
Zimmermann
,
B.
Vorselaars
,
D.
Quigley
, and
B.
Peters
,
J. Am. Chem. Soc.
137
,
13352
(
2015
).
33.
D.
Zahn
,
ChemPhysChem
16
,
2069
(
2015
).
34.
S.
Gottschalk
,
M. C.
Lin
, and
D.
Manocha
, in
Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques - SIGGRAPH ‘ 96
(ACM, New York,
1996
), p.
171
.
35.
P.
Steinhardt
,
D.
Nelson
, and
M.
Ronchetti
,
Phys. Rev. B
28
,
784
(
1983
).
36.
L.
Onsager
,
Ann. N. Y. Acad. Sci.
51
,
627
(
1949
).
37.
G.
van Anders
,
D.
Klotsa
,
N. K.
Ahmed
,
M.
Engel
, and
S. C.
Glotzer
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
E4812
(
2014
).
38.
T.
Schilling
and
D.
Frenkel
,
Phys. Rev. Lett.
92
,
085505
(
2004
).
39.
S. S.
Shapiro
and
M. B.
Wilk
,
Biometrika
52
,
591
(
1965
).
40.
P. R.
ten Wolde
,
M. J.
Ruiz-Montero
, and
D.
Frenkel
,
J. Chem. Phys.
104
,
9932
(
1996
).
41.
P. R.
ten Wolde
,
M. J.
Ruiz-Montero
, and
D.
Frenkel
,
Phys. Rev. Lett.
75
,
2714
(
1995
).
42.
D.
Moroni
,
P. R.
ten Wolde
, and
P. G.
Bolhuis
,
Phys. Rev. Lett.
94
,
235703
(
2005
).
43.
B.
Peters
and
B. L.
Trout
,
J. Chem. Phys.
125
,
054108
(
2006
).
44.
R. L.
Davidchack
,
J. R.
Morris
, and
B. B.
Laird
,
J. Chem. Phys.
125
,
094710
(
2006
).
45.
R. L.
Davidchack
,
J. Chem. Phys.
133
,
234701
(
2010
).