Nucleus-size pinning for determination of nucleation free-energy barriers and nucleus geometry

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 k_BT). For example, for a nucleation barrier of 20 k_BT, 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,¹⁵ NaCl,¹⁸ hard spheres,¹⁷ water, lattice models,¹⁹ and metals.²⁰ Besides homogenous nucleation, the method has been employed to study heterogeneous nucleation phenomena in nanoparticles,²¹ including cases with a cubic seed.²² 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 condition^2,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,²⁷ 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 equilibrate^17,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.³ 

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 pinning^28–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.³¹ 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.³¹ 

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 V_s and surface area A_s, the extensive free energy of the system $ϕT$ in reference to a disordered (liquid) bulk phase is given as

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 V_s. $ϕi$ is proportional to A_s, while the bulk contributions are proportional to V_s or to the number of particles in the seed N_s, thus Eq. (1) becomes

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 V_s because the solid is thermodynamically more stable). For a given seed of fixed geometry, N_s is proportional to the volume V_s, which in turn is related to the length scale of the seed (L) as

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

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

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

Which when substituted in Eq. (6) gives

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

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.²⁸ For a pure component system at a fixed temperature T, the fundamental thermodynamic equation can be written as

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

The total free energy can again be split up as

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

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

Let p^* denote the critical pressure corresponding to a preset, target seed size N_s,t. Assume that an initial isothermal isobaric simulation is conducted at a pressure p_A corresponding to a point A on the free energy curve where the seed size is N_s,A. The goal is to devise a scheme that, given simulation data at p_A, it returns a guess pressure p_B 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:

which are the first and zeroth order approximations, respectively, with $μj,A$ and C_A being the values of properties at point A, and $ΔpBA=pB−pA$⁠. 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 pressure¹⁷ and hence can be approximated as a constant (whose value is C_A). 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

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

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 S_T,B,

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

This estimate of p_B 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 S_T 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

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 N_s = N_s,t, noting that both N_s,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

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

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

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

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.³⁴ 

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$⁠, p_a is the unscaled pressure, k_B is the Boltzmann constant, and σ is the diameter of the sphere. For GBF, the pressure was scaled as $p=βpaac3$⁠, where a_c 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).³¹ 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.¹⁷ 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

The biased NPT MC simulations were performed in a manner similar to those described in prior studies.²⁸ Following the conventional seeding^16,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 × 10⁶ 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(N_s) 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:²⁸ 

where $σNs2$ is the variance in the N_s 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.

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.³ Using as a reaction coordinate the size of the greatest cluster of ordered particles (N_s), the transition path from the disordered to the ordered state is divided into overlapping windows. Each window is simulated separately with reflective walls and N_s 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 N_s 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 N_s = 0 as the reference. All US calculations were performed with a system of 1728 particles.

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 q₆ Steinhardt³⁵ order parameter which is defined as follows: For every particle i, the local bond order parameter, q_l,m(i), is defined by

where N_b(i) is the number of neighbors of particle i, Y_l,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 r_c = 1.4σ of particle i. The translational-order correlation between particle i and its neighbor j, d₆(i,j), is given by

where the asterisk (*) denotes the complex conjugate. Two particles i and j are defined as translationally connected if d₆(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 r_c belong to the same cluster. The tunable parameters in such order parameters need to be optimized using criteria described in prior studies.^3,7

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.³¹ 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 d₆ as follows. For a given pair of neighboring particles i and j (i.e., within the cutoff distance r_c = 1.4a, where a is the side length of GBF), we define the orientational order correlation d_z(i,j) as

where z(i) and z(j) are unit vectors pointing along the long-axis of particles i and j. d_z(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 d_z(i,j) > d_z,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 US⁶ and seeding.¹⁷ 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.

For N_t = 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.¹⁷ 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.

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.⁷ 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 volume³⁶ and maximizing packing entropy,³⁷ 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.⁴ Such poisoning can even preclude nucleation as is the case in the isotropic-nematic transition of colloidal rods;³⁸ 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.

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.

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 N_s ≈ 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 ∼10k_BT, 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 r_c and d_z,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.⁷ 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.

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).

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.

In this study, we introduced the nucleus-size pinning (NSP) method which, as an extension of the seeding technique,¹⁷ 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 normality³⁹ 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.⁴³ 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 C_A 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*<1$¹⁷ 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 (S_i) as

Thus, for a simulation performed at a given pressure, Eq. (A1) can be used to evaluate S_i, 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 N_t = 1300.

It can be seen in Fig. 7 that the major contribution to the variation of the total slope (S_T) 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 S_i = S_i(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 S_i. Assuming for concreteness that the nucleus is spherical, then $kA=γ36πNtρs213$ and $Si=2kA3N13$⁠, from which it follows that

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.⁶ 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

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.

Finally, once we have an estimate for S_i, and hence k_A, 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.