Interplay of multiple clusters and initial interface positioning for forward flux sampling simulations of crystal nucleation

Crystal nucleation from the melt is a widespread and important phenomenon. Furthering the understanding of how this phase transition occurs has important implications for, e.g., climate modeling,^1,2 cryopreservation,³ and the oil industry.⁴ Of particular interest is the rate, per unit time and per unit volume, at which nucleation occurs, J. Simulations of nucleation provide an opportunity to study the freezing mechanism and obtain the nucleation rate under physical conditions and at a level of detail often unachievable through experiments.

Under typical conditions of interest, nucleation is a rare event obeying Poisson statistics. Thus, it is often impossible to perform unbiased molecular dynamics (MD) simulations and observe this phase transition directly. Therefore, enhanced sampling techniques are often needed to circumvent the timescale problem. In this work, we focus on forward flux sampling (FFS). FFS is a widely used tool in, e.g., studying transitions in biological systems,^5–7 as well as crystal nucleation.^4,8–14

Despite the simplicity of FFS, there are a number of user-specified parameters that can lead to significantly different implementations. For example, the initial flux, Φ₀ (discussed in Sec. II A), can be calculated by counting every positive crossing of λ₀ or only crossings when the trajectory is coming from within a liquid basin bounded by some value λ_A ≠ λ₀ (illustrated in Fig. 1).^5,8,15 The impact of these different interpretations is often compounded by nomenclature, which makes the exact implementation unclear,^4,9,10,16 exacerbated by the fact that λ_A plays a significant role in calculating transition probabilities. It is therefore possible to define the edge of liquid basin but not consider its role in the calculation of initial flux. One of the aims of the present work is to investigate the robustness and reproducibility of initial flux calculations so as to determine whether or not the variability of the different FFS implementations in the literature is a concern. To this end, we investigate the effects of different placements of the first interface, λ₀, as well as the position of the interface marking the boundary of the liquid basin, λ_A.

Previous work on optimizing interface placement has centered on computational efficiency and error reduction and has neglected the impact of the initial flux (and therefore the locations of λ_A and λ₀) due to the minimal computational cost in computing the initial flux to a low variance compared to the full transition probabilities. In addition, these interface placement schemes rely either on adjusting interface positions as a result of the true transition probability or on the end location of several long trajectories. These can significantly increase the computational cost of FFS, performing several tests of transition probabilities at each interface in order to determine the “optimum” interface placement.^17,18 The placement of λ₀ has been investigated in the work of Velez-Vega et al. as a means to ensure that configurations stored upon crossings of λ₀ are uncorrelated. However, their technique is of limited applicability, especially in the context of crystal nucleation.⁷ To the best of our knowledge, no investigation of the optimum placement of the edge of the liquid basin, λ_A, has been performed. Henceforth, we shall refer to λ_A, λ₀, and, to a lesser extent, λ₁ as the “initial interfaces” in FFS.

For almost all systems of practical interest, the system sizes relevant to real-life applications are several orders of magnitude larger than those that can be modeled using current computational resources. As such, most computational studies of nucleation rely on the use of periodic boundary conditions. These infinite repeats of the same simulation cell lead to spurious effects that would not be present in the corresponding macroscopic systems. These are often dependent on the size of the simulation cell and are known as finite size effects (FSEs). Several investigations of FSEs in the context of simulations of nucleation can be found in the literature. All of these studies have found evidence of finite size effects when a nucleus interacts with one or more periodic replicas of itself.^19–22 More recent studies have even found evidence of finite size effects in simulations including millions of atoms.²³ Of particular relevance to this work are the investigations of FSEs in the nucleation of mW water (on a surface) and Lennard-Jones (LJ) systems (in bulk) performed by Hussain and Haji-Akbari using jumpy FFS.^13,14 For the mW model, they found that, even in systems large enough for the surroundings of each cluster to behave like the bulk liquid phase, the nucleation rate is not constant with respect to volume [the inverse of a proxy for the side length of the simulation cell was found to be ∝log₁₀(J)],¹³ although this could not be corroborated by similar studies on the Lennard-Jones system due to a lack of data.¹⁴ 

In this work, we address some practical aspects of FFS implementations via systematic investigations utilizing the LJ model. We show that the placement of the initial interfaces is especially important. In particular, it is possible to use a simple and relatively cheap analysis of a MD simulation limited to the initial liquid basin as a means to pinpoint specific positions of the initial interfaces yielding consistent results. We propose a simple heuristic for the placement of initial interfaces in the context of crystal nucleation—placing λ_A at the most probable value of the largest cluster size in the metastable liquid, and the first interface at a location where interactions between nuclei are negligible (which reduces the potential to underestimate effective flux through subsequent interfaces due to merging of nuclei).

We also probe the effect of considering clusters other than the largest cluster (often utilized to define the order parameter)—with specific reference to the calculation of initial flux as well as of the subsequent crossing probabilities. It should also be noted that in the conventional implementation of FFS, it is impossible to consider multiple clusters when calculating crossing probabilities. Despite the effects of multiple clusters being accessible for the initial flux, counting only crossings of the single largest cluster appears to be almost universally done in the literature,^4,8–14,24 although the potential impacts of this choice have not, to the best of our knowledge, been investigated. It should be noted that according to work by Cheng and Ceriotti, the thermodynamically stable state of a nucleating liquid will switch from many smaller clusters to a single large cluster as the size of the largest cluster increases.²⁵ However, this switch may occur for cluster sizes much larger than the location of λ₀.

Finally, we provide evidence for the existence of substantial spatial correlations between different clusters. We find that these correlations lead to the merging of smaller clusters into larger ones—an occurrence that has a non-negligible impact at the supercooling condition considered in this work.

The paper is organized as follows: We begin by summarizing the relevant methodology and simulation details in Sec. II. Results are presented in Sec. III. These consist of relevant nuclei populations, an initial flux analysis of the systems with these populations, and a brief consideration of the consistency of the flux through subsequent interfaces. Spatial correlations between clusters are also presented and discussed. The main findings of our work are then summarized in the context of the current literature in Sec. IV.

FFS is a path sampling technique that can be used for simulating crystallization under conditions where the timescale required for nucleation with brute force MD is intractable. The path between the liquid state, A, and the crystal, B, is described in terms of an order parameter (OP, λ). In FFS, we select a set of isosurfaces of the OP (labeled λ_i), henceforth referred to as interfaces, located along the space of the OP such that moving between successive interfaces leads to a complete path from A to B. Starting from the λ_i−1 interface, when a run crosses λ_i the system configuration is saved. Many statistically independent trajectories are then initialized from these stored configurations. Whether these new trajectories reach the next interface at λ_i+1 or return to a less crystalline state (usually λ_A, defined as the edge of the liquid basin) is then determined. This gives a probability that if the system is in state λ_i, it will progress to state λ_i+1 before returning to state A, P(λ_i+1|λ_i).

The interfaces are typically placed such that the probability of crossing interfaces can be computed with sufficient accuracy. If the probability of crossing the next interface is too low (i.e., the interfaces are too far away from each other), the lack of statistics leads to a high error as a result of insufficient coverage of phase space and sparsity of configurations at interfaces. On the other hand, if almost all trials cross the next interface, little information is gained and there is likely to be a high degree of correlation between configurations at successive interfaces (limiting the accuracy for these new starting configurations). For starting configurations that sufficiently sample the appropriate phase space at the interface, and are uncorrelated, the current consensus in the literature is that the choice of interface positions does not affect the value of the calculated flux although it may affect both the computational efficiency and the calculated uncertainty in the flux.^{16–18,24,26–28}

In this work, we used Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (23 June 2022)³¹ to perform NpT simulations of LJ systems of several different system sizes, ranging from 4000 atoms to 108 000 atoms in cubic periodic boxes. Note that older versions of LAMMPS deal incorrectly with the computation of order parameters via dynamic groups and may therefore not give consistent results.³² Atoms interacted through a shifted LJ forcefield with a cutoff of 3.5 σ and m = σ = ɛ = 1. Thermostatting was applied through a canonical velocity rescaling thermostat³³ set to enforce a temperature T* = 0.86 (here, * shall represent LJ reduced units) with a damping parameter of 0.05t*. This temperature corresponds to an approximate supercooling of 19%. At this supercooling condition, nucleation events can be considered as rare events from a Poisson statistics point of view²² and the re-thermalization of the system should be rapid compared to the timescale of interest.

We also enforced a pressure of p* = 5.68 via a Nosé–Hoover barostat (with a Martyna–Tobias–Klein correction³⁴), where barostat thermalization was applied through a chain of five Nosé–Hoover thermostats. Integration was performed using the velocity Verlet algorithm of time step t* = 0.002.

Crystalline atoms were identified using the sixth-order ten Wolde order parameter (q₆).³⁵ This is a vectorial OP based on the combination of spherical harmonics proposed in the work of Steinhardt et al.³⁶ The maximum distance between neighbors was set to 1.432 σ; a connection was defined as solid if q₆(i) · q₆(j) > 0.5; and an atom with eight or more solid connections was classified as solid. Grouping of crystalline atoms into clusters was performed internally in LAMMPS (see input file). These choices are obviously not unique and may have an impact on the absolute nucleation rates calculated. We do not expect this choice to have a significant effect on the conclusions drawn about initial fluxes, OPs which lead to a difference in the diffusivity of cluster interfaces may, however, change the conclusions about the prevalence of cluster merging. However, as this work is concerned with the self-consistency of FFS calculations, investigations of the impacts of these parameters has not been performed.

It is important to consider the time resolution at which the OP is sampled. Here, we computed the OP at snapshots (or frames) taken every 100 MD time steps. This choice represents a balance between computational expense and time resolution. The consequences of this finite time resolution is discussed later.

To avoid methodological ambiguity originating from the occurrence of clusters of size exactly identical to the value of a given FFS interface, all interfaces were placed at non-integer values of cluster size.

The input files and analysis scripts we have used in this work are available on GitHub at https://github.com/keb721/FFS_Interfaces.

For the initial flux analysis, 24 independent runs were performed for each system size. The systems were all generated by melting a solid system at T* = 2.4 for 10 000 time steps and then cooling via a 1000 time step equilibration at T* = 1.2 (close to but above the melting temperature), followed by a 2000 time step linear quench to the temperature of interest. Data at this temperature were generated by simulating for 1 000 000 time steps. Cluster sizes were output by LAMMPS, with additional post-processing for calculation of fluxes and spatial distribution functions performed by in-house code.

Spatial correlations were estimated by computing a weighted histogram of the minimum distances between clusters above a specified size, denoted g′(r), within the regime where the minimum image convention is applicable (half of the side length of the cubic box). For each cluster of interest, the minimum distance between it and all of the images of all other clusters to consider was computed. The extent of the clusters was considered implicitly, as distances were computed between the location of atoms within the cluster, although atoms are considered point particles. Self-interactions were not included. Weighting was then performed with respect to the volume of the spherical shell of the histogram bin, such that clusters in close proximity were weighted more heavily, and the number of clusters in the frame for which the histogram was being computed—i.e., one less than the number of clusters of interest in the frame, to account for the excluded self-interactions. As the number of clusters of interest for the sizes considered here is not large, the difference between normalizing with respect to the total number of clusters in the frame and the number of considered clusters in the frame could be significant.

Additionally, to test the robustness of the FFS approach, we calculated the flux through a given interface λ′, in two different ways. (1) We calculated the “direct” flux through λ′, i.e., by counting the number of crossings per unit time and volume to yield the flux Φ′. This can also be written as $Φ|λA′$ to make the dependence on liquid basin location explicit. (2) We defined a second interface at λ_0′ < λ′; we computed the flux through λ_0′ by counting the number of crossings per unit time and volume to yield the flux Φ_0′. Then, we calculated the “effective” flux through λ′, by multiplying Φ_0′ by P(λ′|λ_0′) to give Φ_λ′. To account for the effect of λ_A in both Φ_0′ (i.e., $Φ0′|λA$⁠), and P(λ′|λ_0′), the effective flux can also be denoted $Φλ′|λA$⁠. Note that the choice of λ_A will always be applied consistently to both the initial flux and crossing probability components. This so-called “λ′ test” allows us to establish whether the positioning of a given interface has an impact on the effective flux—which has an impact on the overall nucleation rate. Specifically, if the direct and effective fluxes computed via (1) and (2), respectively differ, we have a clear indication of the uncertainty associated with the positioning of λ_A and λ_0′ (and in fact, of any other FFS interface) on the nucleation rate calculated via Eq. (2).

To perform the λ′ test, the direct fluxes were taken from the initial flux analysis, as described above. For calculation of the effective flux via method (2) above, configurations at a range of λ_0′ values between 20.5 and 40.5 were generated. Liquid systems at the temperature of interest were generated by melting a solid system and cooling it to the temperature of interest using the same procedure as above. Relevant configurations were stored every time the size of the largest cluster, λ, reached $(12+λ0′)$⁠. Special care was taken in ensuring that these stored configurations were not correlated with each other. To this end, every time that a configuration at λ_0′ was stored, a simulation of 1000 MD time steps was performed, during which the value of λ was ignored. After these 1000 MD time steps, monitoring of λ was resumed, to check for the condition λ < λ_P. Note that λ_P is system size dependent, as illustrated in Table II. At this stage, we can be certain that the system was sufficiently de-correlated with respect to the previous crossing of λ_0′, and a new configuration was stored when λ reached the appropriate value.

We accumulated 500 configurations at each λ_0′, generated according to the procedure above. Then, we ran 10 000 MD “trials” starting from a random choice of these de-correlated configurations, redrawing particle velocities from a Maxwell–Boltzmann distribution relative to the temperature of interest.

For the initial fluxes obtained at different system sizes to be equivalent, they must exhibit the same population of nuclei in the metastable liquid state. This population is a function of supercooling. Figure 2 represents the populations of nuclei in the simulation (before the onset of the nucleation, monitored by means of the order parameter described in Sec. II B) for all system sizes studied. Figure 2(a) shows, on a log-linear scale, the probability density function (PDF) of all nuclei in a simulation, regardless of the size of any other nuclei in the system. This PDF has been normalized with respect to unit simulation volume (see Table II for the relative volumes of the different simulation cells). The colored lines represent the mean value of the PDF and the shaded regions represent the standard error on the mean. At low cluster sizes, the errors are negligible and the distribution of nuclei is the same for all simulation volumes. As cluster size increases, the errors become more significant as a result of insufficient statistics as large clusters correspond to large fluctuations, which are less accessible to brute force simulations. The longer tails for the larger simulations simply reflect the fact that as there are more nuclei in the system. It is also more likely that one of these nuclei undergoes statistical fluctuations to grow to a larger size. It should be noted that the largest sized cluster in all simulations was under 175 atoms, indicating that no nucleation events have taken place in any of the simulations.

These results demonstrate that although the distribution of total nuclei populations in the liquid basin is volume-independent, the modal size of the primary cluster is not as it is related to the average number of clusters in the simulation volume. This has implications for the definition of the liquid basin in studies that use primary cluster size as an OP. Carefully defining the location of the liquid basin is important in path sampling techniques when determining interface crossing probabilities and, in some cases, the initial flux. Using the same value of λ_A for different system sizes may lead to trajectories committing to the A basin not being treated as such. This can cause both large differences in calculated nucleation rate due to an artificially increased crossing probability and much longer simulation times (see Sec. III C).

As the population of nuclei per unit volume was comparable across different system sizes (see Fig. 2), we also expect the flux through the first interface, Φ₀, to be independent of system size, as it is normalized with respect to simulation volume. In order to compare Φ₀, it is important to ensure collection of sufficient statistics to calculate a stationary flux, the number of crossings through λ₀ should be approximately linear with time. The crossings of λ₀ = 36.5 are presented in Fig. S1, see the supplementary material.

Figure 3 offers an analysis of the flux through the first interface λ₀. Here, the solid line represents the mean, with the error bars (less than the linewidth) representing the standard error on the mean. Figures 3(a) and 3(b) correspond to calculating the flux by counting every positive crossing of λ₀, $Φ0|λ0$ (i.e., any negative crossing of λ₀ is a return to the liquid basin). In contrast, Figs. 3(c) and 3(d) represent the flux, $Φ0|λP$⁠, when the edge of the liquid basin is placed away from λ₀. Again, crossings are counted only if they originate from the liquid basin, but in this implementation this requires the cluster size to fall below λ_P ≠ λ₀ between crossings of λ₀ (see Fig. 1). This is the procedure used in some, but not all, FFS implementations.^5,8,15 The λ_Ps are defined as the peaks in the primary incidences, reported in Table II. Recall that in order to prevent ambiguity regarding the treatment of clusters characterized by size exactly equal to an interface value, both λ_A and λ₀ were chosen to contain a non-integer number of atoms.

When considering the flux of the largest cluster only, the time series of the size of the primary cluster was considered. For the total positive flux [i.e., $Φ0|λ0$⁠, Fig. 3(b)], the number of positive crossings of λ₀ in this time series were counted, before being divided by the length of the simulation and the average volume of the simulation box. A positive crossing is one in which λ₀ is crossed in the direction of increasing OP. As λ₀ defines not only the first interface but also the edge of the liquid basin, this is equivalent to counting all crossings originating therefrom. For Fig. 3(d), the edge of the liquid basin was instead placed at λ_P ≠ λ₀. Under this definition, not all positive crossings of λ₀ originated from the liquid basin. Therefore, more care must be taken to ensure that only the first crossing of λ₀ after leaving the liquid basin is counted. The influence of noise, which could lead to a large increase in the number of positive crossings due to a rapid fluctuation of cluster size, was mitigated by subsampling the time series (as described in Sec. II).

It can be seen in Fig. 3(a) that there is a strong dependence on the initial flux of all clusters with system size (when the edge of the liquid basin is placed at λ₀), despite the comparable nuclei populations. This dependence is most pronounced at small λ₀ and is partly due to the method used for calculating the fluxes of all clusters. At each snapshot, the number of clusters of size $>λ0$ is counted. If, at the next snapshot there are more clusters of size $>λ0$⁠, then the flux is increased by this difference. This is equivalent to counting the total number of net positive crossings between frames. Positive crossings of λ₀ will therefore be missed if negative crossings of λ₀ have also occurred between snapshots.

The inability to identify concurrent growing and shrinking of clusters is also likely to lead to an underestimate of Φ₀ for the flux of primary nuclei, although there are now other important considerations. Comparing Figs. 3(a) and 3(b), it can be seen that (for the same liquid basin, bounded by λ₀) the flux of primary nuclei is lower than the flux of all nuclei for all system sizes at low λ₀. The primary flux is equal to 0 for the 108 000 atom system below λ₀ ≈ 25 despite there being appreciable flux when all clusters are considered. The spread of Φ₀ as a function of system size is also larger at lower λ₀ in Fig. 3(b) than Fig. 3(a), and the stratification persists up to relatively large λ₀. Additionally, it is of interest that the peaks in the flux in Fig. 3(b) do not correspond to the values of the peaks in the primary incidences [Fig. 2(b)], instead occurring at slightly higher cluster sizes. This shows that (for λ_A = λ₀) most flux at small λ₀ occurs as a result of non-primary nuclei for large systems. This is not unexpected, as the probability of having multiple large clusters increases rapidly with the simulation volume. If λ₀ is chosen to be small, and/or the simulation volume is large, the contribution of non-primary clusters is significant.

Figure 3(a) shows only the positive crossings of λ₀, with no restrictions on how far a cluster must shrink back into the liquid basin before a subsequent boundary crossing—i.e., the edge of the liquid basin is taken to be λ₀ itself. In contrast, Fig. 3(c) shows $Φ0|λP$⁠, the initial flux when the edge of the liquid basin is placed away from λ₀, here λ_P (see Table II). Instead of counting the net crossings of λ₀, as in Fig. 3(a), an ordered list of the sizes of all solid clusters is considered. If the size of the Mth largest cluster has fallen below λ_A, then the next crossing of λ₀ by the Mth largest cluster will be counted. Any subsequent crossing of λ₀ by the Mth largest cluster will be ignored until it has again fallen below λ_A. Importantly, the Mth largest cluster is not a fixed identity linked to a certain combination of atoms. Instead, it simply gives the size of the Mth largest cluster in each particular frame. Cluster identities of the largest clusters may therefore vary completely between successive snapshots. In practice, it may be useful to consider only the first N elements of this ordered list, where N is chosen such as to ensure that the last element of the subset never exceeds λ₀ over the course of the simulation, which will be system specific. Note that for the two largest system sizes, the lines do not begin at λ₀ = 20.5 due to λ_P > 20.5.

Separating λ₀ and the liquid basin decreases Φ₀ for all clusters and removes the dependence on simulation volume when all clusters in the simulation volume are considered [Fig. 3(c)]. When the edge of the liquid basin is placed at λ₀ such that all positive crossings are counted [Fig. 3(a)], then fluctuations of a cluster’s size around λ₀ will artificially increase Φ₀. The impact of these fluctuations can be reduced by subsampling, which also reduces crossings as a result of noise in the OP, although this is unlikely to be sufficient for removing variation in the OP as a result of realistic fluctuations in the cluster size. Instead, ensuring that a cluster must vary in size more significantly to return to the liquid basin means that it is much more likely that the next crossing of λ₀ will be as a result of the growth of a different, independent cluster, increasing the robustness of the measure of the flux. This separation of λ_A and λ₀ has the most significant impact on the initial flux for the smallest system sizes. This is likely to be due to a combination of two related effects. First, for a smaller simulation volume there are fewer nuclei and therefore it is more likely that a rapid recrossing of λ₀ occurs as a result of the same cluster changing size. Additionally, for a smaller system the liquid basin is further away (due to the lower number of clusters) and therefore ensuring a cluster has returned there is more likely to remove a greater number of crossings.

Figure 3(d) displays $Φ0|λP$ for only primary clusters. Unlike $Φ0|λP$ for all clusters [Fig. 3(c)], a size-dependence in the initial flux is still observed, despite removing the contributions of small fluctuations in cluster size by ensuring a return to λ_P ≠ λ₀ between subsequent crossings. This demonstrates the importance of non-primary nuclei in accurately calculating the initial flux. Especially in larger simulation volumes and at lower values of λ₀, the difference between primary flux and (size-independent) total flux is pronounced. As the population of primary clusters is size dependent [see Fig. 2(b)], this is not unexpected.

Despite having the same liquid basin (as described by nuclei populations), it can be seen that different methods of calculating the initial flux can lead to different, volume-dependent results as they neglect the, possibly appreciable, fraction of significant non-primary nuclei and because they incorrectly define the edge of the liquid basin. It can be seen that incorrect placement of λ_A leads to unphysical effects in the flux—either artificially increasing it by inclusion of insignificant fluctuations if placed too high, as in Fig. 3(a), or underestimating Φ₀ as a result of real crossings not being counted as the edge of the liquid basin has been placed at a value that makes observing a return unlikely. If the λ_A chosen is not tailored to both the specific volume and composition of the system being simulated, either of these is a possibility and it can be impossible to spot from the value of the initial flux even with rigorous testing.

In the context of initial flux, the error introduced by incorrect placement of interfaces is dependent on both the λ₀ and λ_A values chosen. However, uncertainty in the nucleation rate is usually considered to be dominated by uncertainties in the transition probability, and the uncertainties considered here are small on the scale of variability of the results of FFS calculations. How the of placement of λ_A and λ₀ propagates to changes in transition probabilities through subsequent interfaces is discussed in more detail below.

It should be noted that there are preexisting heuristics for placing the initial interfaces in crystal nucleation, although to the best of our knowledge these are not the conclusions of in-depth studies. A relevant heuristic for crystal nucleation is that used in Ref. 8, which suggests placing λ₀ between the top 1% and top 0.1% of the OP distribution [for our system, this distribution is shown in Fig. 2(b)] and λ_A in the range [μ, μ + σ] where μ and σ are the mean and standard deviation of the OP distribution, respectively. Our results indicate that this placement of λ₀ will ensure that the vast majority of the observed flux will be contributed by primary nuclei only, although the effect of non-primary nuclei on the transition probability through subsequent interfaces will be investigated in Sec. III C.

We also find that placing λ_A at the mode of our observed OP distribution leads to a consistent initial flux (on the chosen scale) when considering all nuclei. However, this is also the case at high λ₀ for λ_A = λ₀ and represents consistency only in a small part of the initial flux, not necessarily the entire FFS calculation. In addition to showing its consistency, we have posited that the rationale for using λ_A as the most probable value of the largest cluster size in the metastable liquid is to best ensure sufficient counting of independent crossings of different clusters. We acknowledge that this is not a perfect method to ensure that all independent crossings are counted, but such a method does not exist. While we have not explicitly studied the case of μ + σ, the possibly significant overestimate of the liquid basin may lead to an assumption of independence between clusters where it does not exist and a premature termination of phase space exploration when computing transition probabilities—similar to what is expected for λ_A = λ₀, although to a lesser extent. In addition, the shape of the OP distribution may be sensitive to changes in cluster definition, and therefore the mean and the mode may not be coincident and the standard deviation may become significant. While changes in cluster definition may also exacerbate other differences between observed clusters, such as arrangement³⁸ and polymorph,³⁹ which may also influence the nucleation rate, this is not always the case.

In a real FFS run, the initial flux Φ₀ is simply a component that is then combined with the probabilities of crossing subsequent interfaces to give the nucleation rate. The location of λ₀ is of no consequence to the crossing probabilities starting at subsequent interfaces (i.e., λ₁ and above). In contrast, the choice of λ_A is relevant to all interfaces, as it defines the point at which a trajectory is considered to have returned to the liquid basin.

The current consensus is that the placement of interfaces in FFS is important only in terms of the efficiency of and uncertainty in the rate calculation—it does not affect the value of the calculated rate, assuming that interfaces are sufficiently well-spaced that stored configurations are uncorrelated and that they are sufficient in number to sample the relevant phase space.^{16–18,24,26–28} If this is true, then it should be possible to retrieve the direct flux through λ₁ by multiplying the flux through λ₀ < λ₁ by the probability that (uncorrelated) configurations with λ = λ₀ will progress to λ₁ before returning to λ_A for any consistent choice of λ_A, λ₀, and λ₁. The methodology used in this section is outlined in detail in Sec. II B 2. For consistency and clarity, we shall use the nomenclature outlined there, although it is worth noting that in practice the interface λ′ corresponds to both the λ₀ (when computing the direct flux) and λ₁ interfaces (when computing the effective flux).

The direct and effective fluxes of primary nuclei through two different values of λ′, for two choices of λ_A, are given in Fig. 4. Similar to the initial fluxes, these results should be volume-independent. In the initial flux, consideration of non-primary nuclei was required to eliminate significant stratification with system size. Conversely, for the effective flux, including the contribution of non-primary nuclei in Φ_λ′ reduces agreement between system sizes (see Fig. S2 in the supplementary material, noting the larger shaded region of the effective flux, indicating that there is more spread between different volumes when considering all nuclei). Unlike in the initial flux stage, non-primary nuclei cannot be explicitly considered when calculating transition probabilities. Although they cannot be directly accounted for, significant non-primary nuclei are likely to be present in stored configurations with λ = λ_0′ and may therefore influence the flux. Any growth of an initially non-primary cluster, X, above the size of the initially primary cluster, Y, will be attributed to cluster Y. This will spuriously increase the transition probability as a result of the growth of clusters smaller than those present at the λ_0′ interface, although these are ostensibly the only relevant clusters in the system. As previously discussed, the prevalence of significant non-primary nuclei is more likely in larger simulation volumes. The greatest increase in transition probabilities—implicit in effective flux calculations—will be under the conditions where the prevalence of non-primary nuclei leads to the greatest decrease in initial flux, as shown in Fig. 3 (i.e., large simulation volumes and low λ₀). It is encouraging that the impact of non-primary nuclei in both of these computations appears to offset each other, leading to similar effective fluxes for all simulation volumes. It should also be noted that at sufficiently large primary clusters, the size of any non-primary cluster becomes negligible,²⁵ and therefore these need not be considered.

The impact of atomic velocities on the effective fluxes was also investigated (Fig. S3 in the supplementary material). There was little difference observed when the atomic velocities used when initializing trajectories at λ_0′ were those stored during positive crossings of λ_0′ during the flux stage or reinitialized randomly, indicating that the OP dynamics are overdamped with respect to particle momenta on the timescale of interface crossing.

Figure 4(a) corresponds to the flux through λ′ = 40.5 using λ_A = λ_0′, whereas Fig. 4(b) corresponds to the flux through λ′ = 60.5, again with λ_A = λ_0′. Note that in Fig. 4(b), the variation of effective flux with $λ0′$ is too small to be seen on this scale and is lower than the variation in Fig. 4(a) due to fewer crossings of 60.5 compared to 40.5. For both values of λ′, the direct flux is significantly above the effective flux for all values of λ_0′. This choice of λ_A leads to an artificially reduced transition probability—small deviations of the OP below λ_0′ lead to immediate termination of trajectories, underestimating the true likelihood of reaching λ′. This underestimation is likely to be most pronounced at interfaces in the vicinity of A, where the amount of natural nuclear variation possible before a cluster has been assigned as committed to the liquid basin is greatly reduced. In addition, the low probability of crossing λ′ may lead to a lack of stored configurations at the next interface, which has severe implications for the accuracy of the results of further stages.

Conversely, if λ_A were chosen to be significantly below the modal value of the primary cluster, the transition probability would be overestimated (shown in the supplementary material in Fig. S4). In order for a crossing to λ′ to be determined to have failed, all other clusters in the simulation must also be below λ_A when the primary cluster shrinks back into the melt. For a choice where this is unlikely (due to concurrent growing and shrinking of clusters), there is a large probability that another cluster will grow to exceed λ′ before all clusters fall below the chosen λ_A value. As this is impossible to account for, there will be a significant and unrepresentative increase in the calculated transition probability as a result of this. In addition, there is likely to be an increase in simulation time for calculating transition probabilities—either due to waiting for the statistically improbable situation of all clusters being below λ_A or due to waiting for a new cluster to grow to cross λ_i+1, possibly from below λ_i. The overcounting of crossings as a result of requiring over-commitment to the liquid basin is likely to be an insidious problem that will significantly affect all interfaces below the critical nuclear size.

Figure 4(d) shows the direct and effective fluxes through λ′ = 60.5 using λ_A = λ_P. In comparison to Fig. 4(b), there is now a good agreement between direct and effective flux observed through λ′ = 60.5. However, in Fig. 4(c), where λ_A = λ_P but λ′ = 40.5 does not display this consistency. There is a volume-independent (as indicated by a lack of significant broadening in the shaded region representing the spread of results at different volumes) and systematic decrease of effective flux with increasing λ_0′. It should be noted that this is also visible toward the right-hand side of Fig. 4(d). There are four possible reasons for this.

First, FFS assumes that the time taken to move between successive interfaces is negligible compared to the time taken for the initial crossings, although this will not always be true. For the 4000 atom system with λ_A = λ_P, λ_0′ = 38.5, and λ′ = 40.5 (where the discrepancy between direct and effective flux is largest), the average time taken between subsequent crossings of λ_0′ is 18.2 ± 0.3 t*, with the average transition time between λ_0′ and λ′ being 0.774 ± 0.015 t*. This shows that the time taken to move between interfaces is of the same order of magnitude as the uncertainty in the mean time between crossings of λ_0′. It can hence be neglected. Therefore, effects from slow dynamics in the crossing stage are unlikely to be the cause of the systematic decrease in effective flux with λ_0′.

Second, this discrepancy may indicate that the delay between snapshots is too long and important crossings are therefore being missed. As the separation between λ_0′ and λ′ decreases, fewer additions to the nucleus are required for the cluster to cross λ′. Assuming that the average rate of monomer addition is proportional to cluster surface area, then not only do clusters with a small λ′ − λ_0′ need fewer monomer additions to exceed λ′, but (assuming constant λ′) these additions will happen on a faster timescale. If there is a long interval between snapshots, it is probable that these short-time fluctuations will be missed, while longer timescale fluctuations are still observed. Increasing the separation between λ_0′ and λ′ does not eliminate the problem of missing crossings due to sampling intervals, although it may become less meaningful as a result of several factors: the larger chance to explore the phase space; the larger expected crossing time; and the decreased probability of being able to cross λ′ at all. However, decreasing the sampling interval not only results in capturing legitimate interface crossings as a result of growing clusters; fluctuations in the cluster size due to an imperfect OP are also likely to be captured and counted as they are indistinguishable from real transitions. Again, these fluctuations are less important for a trajectory truly returning to the A basin than for one crossing λ′ (or λ_0′)—small variations are more meaningful in the rarer region of OP space away from the metastable basin and therefore OP fluctuations are more liable to unavoidable misinterpretation. In addition to any impact on crossing probabilities, changing the sampling interval will necessarily also influence the observed flux. Although this change will be most pronounced for λ_A = λ_0′, due to the inclusion of more observed rapid and statistically meaningless fluctuations of nucleus size (whether occurring as a result of noise on the OP or not), it will have an effect in all other cases as well, e.g., by capturing additional forays back into the liquid basin. Decreasing the sampling interval was not able to resolve the systematic decrease in effective flux with increasing λ_0′ (see supplementary material, Fig. S6), although it does exhibit the expected increase in crossing probability.

Third, in a realistic FFS run with a finite sampling frequency, the configurations stored at interfaces are likely to be those that have crossed an interface,^{7,8,16,17,24,28} which do not necessarily lie exactly on that interface. This means that the ensemble of stored configurations at an interface is likely to include clusters significantly larger than the interface size. Numerical work by Haji-Akbari has indicated that neglecting these configurations is likely to lead to an underestimation of the flux.⁴⁰ For our work, including configurations that have crossed λ_0′ but have a cluster size above $λ0′+12$ at λ_0′ = 30.5 (shown in Fig. S3 in the supplementary material) did not resolve the discrepancy between direct and effective flux, despite the large proportion of stored configurations with sizes much larger than 31, deemed to be configurations “at the boundary.”

Finally, for crystal nucleation in LJ there is evidence that cluster growth occurs not only as a result of addition of single particles to a cluster but also through merging of multiple clusters (two or more) into a larger cluster.^41–43 This merging of clusters would violate the assumption of FFS that the only pathway between A and B is one that passes through every intermediate interface sequentially.⁴⁴ If this sequential crossing were not the case, then the effective flux would underestimate the true flux, as it would not be able to take into account the possible pathway of crossing λ′ by merging of multiple clusters of a size less than λ_0′. Such an event would be included in the calculation of direct flux through λ′ but omitted from the effective flux as at no time is there a cluster with λ = λ_0′ along this pathway. This would explain the steady decrease in effective flux in Fig. 4(c) with increasing λ_0′—the smaller the difference between λ_0′ and λ′, the more likely it is that a merging event would allow for a crossing of λ′ without requiring a crossing of λ_0′. This would also explain the small deviation at high λ_0′ in Fig. 4(d). The possibility of merging clusters in our system is explored in more detail in Sec. III D.

The potential influence of merging clusters on the effective flux can be minimized in a number of ways. First, λ_0′ can be placed at a sufficiently rare value of the OP that there is only ever one significant cluster at that size. This also means that the contribution of non-primary nuclei in the initial flux can be neglected, although crossings of λ_0′ may then be so rare that there is a large uncertainty in the flux, and that it is difficult to store configurations that sample a sufficiently broad region of phase space. Alternatively, λ_0′ and λ′ can be spaced such that the likelihood of not crossing λ_0′ on the way to λ′ is negligible. Although this is unlikely to cause sampling issues at λ_0′, the low probability of reaching λ′ may lead to insufficient sampling at subsequent interfaces. Additionally, instead of using conventional FFS to determine the nucleation rate, the jumpy forward flux sampling (jFFS) algorithm proposed by Haji-Akbari can be used in order to take into account an OP that may vary significantly enough that not all interfaces may be crossed. However, to minimize computational cost, the suggested implementation of jFFS involves adaptively placing interfaces, such that they are sufficiently far apart for the jumpiness of the OP to be neglected.⁴⁰ Again, care must be taken to ensure adequate sampling at higher interfaces.

We have shown that consistency in the choices of λ_A, λ₀, and λ₁ does not automatically lead to agreement between direct and effective flux. Although we do not explicitly consider the case of considering a different liquid basin when computing crossing probabilities and the initial fluxes, the differences in initial fluxes for λ₀ = 20.5 for $Φ0|λ0$ and $Φ0|λP$ are of the order of 2 (see Fig. 3), but this difference is significantly less pronounced for λ₀ = 60.5. It is therefore obvious that, as well as the methodological inconsistency in the definition of the liquid basin, the results obtained when using $Φ0|λ0$ and λ_A = λ_P for the crossing stage are also inconsistent. We also expect this conclusion to hold for other inconsistent choices.

When attempting to minimize the errors in the crossing probabilities in FFS, the major issue considered is simply one of uncertainty due to a lack of crossing statistics. This can be minimized by increasing the number of trial runs performed. However, here we demonstrate the existence of significant systematic errors due to poor consideration of interface placement, which cannot be eliminated by generating more statistics. In the worst case, these systematic errors will be present for every subsequent interface, meaning that a relatively small error on a transition probability can be compounded into several orders of magnitude for the overall nucleation rate. Assuming a relatively modest difference in observed transition probability of 10%, similar to that which we have found for different choices of λ_A for the first four interfaces (see the supplementary material, Fig. S5), could lead to an uncertainty in the nucleation rate of 45%, which cannot be rectified by the inclusion of more trials at interfaces. This can be very significant in some cases, e.g., when comparing the dominance of two process with closely competing rates, as in Ref. 39.

Again, we have shown that a correct definition of the liquid basin is critical to ensuring consistency. We suggest that merging of clusters is a potentially important pathway to be considered when performing FFS of nucleation, the effect of which can be minimized in several ways. We also give evidence that the contributions of non-primary nuclei to transition probability are implicitly considered.

In Sec. III C, we discuss the possibility that the systematic difference between effective and direct flux for small separations between λ_0′ and λ′ may arise as a result of clusters merging. In order for clusters to merge within a timescale that would impact the flux, they must be located near each other to some degree. As such, we investigated potential spatial correlations between the crystalline clusters in the system. To do so, we have computed a weighted histogram of all of the minimum distances between clusters, here denoted as g′(r), for all clusters above a given size. Note that the actual extent of the clusters has been taken into account when determining these distances. These were binned into spherical shells of thickness 0.1 σ, with the furthest extent of the final bin being given in Table III to ensure that configurations that violate the minimum image convention are never included. Using this procedure, the incidence of single solid particles in the melt gives a g′(r) similar to the radial distribution function, g(r), expected for a liquid (see the supplementary material, Fig. S7). This shows that single solid particles are reasonably common and will occur with approximately the same probability at any location in the simulation box—as expected.

These weighted histograms are given in Fig. 5 for several different cluster sizes. As the size of clusters increases, the distribution retains some of the character of the g′(r) of single solid particles, still showing well-defined peaks, which become sharper and more pronounced due to the decrease in number of clusters. Unlike single solid particles, nuclei containing 20 or more atoms are expected to be sufficiently rare that there should be no obvious correlation between their observed locations. As can be seen in Fig. 5(a), we observe distinct peaks within the short/medium-range order region of the g′(r) across all system sizes. It should be noted that these g′(r) are intrinsically biased toward smaller distances, as only the minimum distances between clusters are considered. However, if the peaks in the g′(r) were occurring solely as a result of this, we would observe a flattening of the peaks as the simulation volume increases due to rise in the maximum possible distance between nuclei. This is not the case, as the position of the peaks does not change as a function of the system size. In fact, this result indicates that there is some degree of clustering of nuclei around other nuclei, which can also be clearly seen when visualizing the trajectories [see Fig. 6, obtained via the visual molecular dynamics (VMD)⁴⁵ software]. Although there would inevitably be incidences of clustering for truly randomly distributed clusters, the combination of peaks in the g′(r) and the apparent frequency of the observed clustering in visualized trajectories indicates that these clusters are not completely independent. In addition, the observed shape of the peaks maps reasonably closely to a hard sphere model with two preferred, related cluster separations (see the supplementary material, Fig. S8).

When the size of the clusters increases, the peaks become less pronounced. Considering Fig. 5(b), it is clear that although the g′(r) of the smaller systems is much noisier, the peaks are still well defined for the larger simulations. As the cluster sizes increase still further to 60 and 80 atoms [Figs. 5(c) and 5(d) respectively], the noise continues to increase and the g′(r) goes to 0 for the smaller volumes, as there are no frames containing more than one cluster of the relevant size. The cluster sizes considered in Figs. 5(c) and 5(d) are very large in the context of a unbiased simulation under these conditions. Although having multiple large clusters in the same frame becomes more likely as simulation volume increases, the probability of this happening is still very small. This makes it unclear to what extent the apparent loss of clustering is simply due to a lack of statistics collected and how much is due to other, unknown effects.

Although there is evidence of clustering (see Figs. 5 and 6), the exact origin of this is unknown. A possible explanation is that clustering arises due to the breakup of clusters. If an originally dumbbell-like solid cluster were to lose solid atoms from the bridging section of the nucleus, the result would be two nearby clusters of appreciable size. This is difficult to directly observe without detailed tracking of cluster identities.

Alternatively, clustering may suggest structuring of the LJ melt around an established nucleus, potentially due to the diffuse interface of crystalline clusters, in a way that increases the propensity for other clusters to form in the vicinity. Hussain and Haji-Akbari found evidence for nuclei structuring the liquid around critical nuclei, showing distinct peaks and troughs in the liquid density before reaching a plateau.¹⁴ Although the results are not comparable due to different simulation conditions and the fact that the nuclei considered here are far from critical, this indicates that the presence of nuclei does structure the surrounding liquid, perhaps making it easier for other nuclei to form.

In addition, the increased presence of significant solid clusters in the liquid around other solid clusters may indicate a propensity for nuclei to grow by amalgamation, since growth or movement of the cluster in any direction may lead to an encounter with another nucleus, with which it would then merge. This is shown in Fig. 6. It should be noted that this occurs at too fast a timescale to be observed when sampling at every 100 time steps (here additional snapshots have been generated with single time step separation), although the effects of merging occurring at faster timescales are still important to results obtained at higher sampling separations.

We have demonstrated the importance of careful interface placement when studying crystal nucleation using FFS. Our data clearly show that the correct placement of λ_A is critical to ensuring a consistent and volume-independent flux. We show that the placement of λ_A and λ₀ in nucleation studies is and should be simulation size dependent. We have given evidence for cluster merging with significant influences on effective flux, although this can be minimized by judicious placement of λ₁ and λ₀. We also illustrated the importance of non-primary nuclei but indicated that their contribution is likely to be implicitly considered when performing a complete FFS run. The importance of these findings also extends to other enhanced sampling techniques using interface placement, e.g., transition interface sampling.⁴⁶ 

The initial implementation of FFS did not include a λ_A interface.⁴⁴ This was instead included in a later paper offering refinements, although it was stated that λ₀ = λ_A was always a valid choice, and often a convenient one.²⁸ Some implementations of FFS have indicated that there is an advantage to placing λ_A away from λ₀ but usually only as a way of removing correlations between stored configurations at λ₀.^27,47 In contrast to this earlier work, we have shown that the choice of λ_A is crucial to ensure a consistent effective flux across λ₁, which is likely to extend to every subsequent interface. Unfortunately, the nomenclature often used in the literature to describe how the initial flux was calculated is ambiguous^4,9—sometimes with reference to counting crossings leaving the A basin^10,16—which makes it unclear if the crossings of λ₀ being counted are simply all positive crossings or positive crossings that have occurred after a return to the A basin (i.e., if what is being computed is $Φ0|λ0$ or $Φ0|λA$ for λ_A ≠ λ₀). This is often compounded by the fact that the value of λ_A used, if any, is omitted. This then extends the methodological confusion to the treatment of crossing of subsequent interfaces, as the condition for determining an unsuccessful transition attempt is undefined.

Previous investigations of interface placement have centered on computational efficiency and reducing the variance in the final nucleation rate and have neglected the location of λ_A and λ₀ due to the minimal computational cost in computing Φ₀ to within a small variance.^17,18 In the work of Velez-Vega et al., the optimum placement of λ₀ has been explored by finding the minimum of the average simulation time needed to generate uncorrelated crossings. However, the need for an additional OP that can be cheaply computed and is distinct from λ significantly limits the applicability of their technique, especially in the context of crystal nucleation where OP choices are generally variants of the Steinhardt parameters.^8–14 In addition, they do not attempt to refine the placement of λ_A, even using λ₀ < λ_A, which is likely to lead to increased correlations between configurations stored at λ₀.⁷ In this work, we have proposed a simple alternative heuristic for the placement of initial interfaces in the context of crystal nucleation—placing the edge of the liquid basin at the peak of the incidences of primary nuclei, and the first interface at a location where the effects of non-primary nuclei are negligible (which reduces the potential underestimate of effective flux through subsequent interfaces due to merging of nuclei). We note that the value of the edge of the liquid basin can continue to have a significant effect on transition probabilities at interfaces beyond λ₀, as it marks the failure criterion for a cluster to reach the next interface (see the supplementary material, Fig. S5). We find that the flux through λ₃ can differ by approximately an order of magnitude depending on where λ_A is placed, which is significant on the scale of uncertainty in a FFS calculation.

In systems where nuclei do not merge, it is only important that λ₀ is not placed in a region where primary flux is negligible, and in these cases placing λ₀ between the top 1% and the top 0.1% of primary incidences (as in Ref. 8 and suggested in Ref. 27) will be more than sufficient. However, when nuclei can merge, it is not only nuclei of a comparable size to the primary nucleus that are important—any significantly sized nucleus can be involved in a merging event that can bypass an interface and thus influence the rate. Therefore, for only primary nuclei to be relevant, λ₀ should be placed after it becomes more thermodynamically favorable to have a single nucleus rather than several smaller nuclei.²⁵ Depending on the properties of the system of interest, this system-specific condition is likely to be difficult to confirm and occur at a larger value of λ₀ than is practicable. In these scenarios, the effects of merging can be mitigated in a number of other ways.

While this work was performed under conditions where multiple nuclei are common, we believe that the results are applicable to systems where this is not the case—e.g., at low supercooling, low interfacial free energy, or for heterogeneous nucleation. A λ_A that is placed too high will count rapid, unstatistical recrossings of λ₀ regardless of why these occur. Similarly, it ensures that the counted flux is indeed the flux of an attempted nucleation event. It is interesting to note that the conclusions presented here regarding interface placement are similar to those presented in recent work of Zhao and Li in heterogeneous mW ice nucleation occurring through two pathways of substantially different rate. In their case, positioning λ₀ at higher values was to increase sampling of the reaction pathway that dominated at high λ, which is somewhat analogous to avoiding growth through the pathway of merging which is impossible at sufficiently large λ.³⁸ 

We also investigated the effects of non-primary nuclei on direct and effective flux. Although a contribution from non-primary nuclei is not unexpected, to the best of our knowledge it has not been explored or discussed in the literature. This is somewhat concerning as the size of primary crystalline nucleus is a commonly used OP in studies of nucleation. In this work, we have demonstrated that there is an appreciable FSE in the initial flux that is caused by non-primary clusters and that non-primary clusters can affect crossing probabilities as they are present in stored configurations. However, we have shown that these effects combine destructively and explicit consideration of non-primary clusters is therefore not necessary.

Finally, we have presented strong evidence for clustering of nuclei around one another. We have found indications that these clusters merge and that this can have an appreciable impact on the observed flux through an interface.

It is important to note that this detailed investigation of FFS implementations is possible due to the low computational cost associated with the LJ potential, especially compared to more realistic models whose results are of greater interest. We have been able to systematically investigate the effects of initial interface placement and, from this, have given a procedure for determining the placement we have found to be necessary without requiring costly additional computation, which is of practical use when investigating more complex models.

We also acknowledge that the choice of OP made may also influence the nucleation rate, although the potential effects of this are significantly less important than a lack of self-consistency within the FFS calculation utilizing the same OP.

The supplementary material for this paper includes graphs of the number of boundary crossings as a function of time and the radial distribution function of single solid particles; several graphs showing different implementations of the λ′ test—the influence of λA placement, the effect of using all nuclei for the fluxes, and the impact of keeping velocities and of storing configurations not “on the interface”; and change in direct and effective fluxes as a result of sampling interval. In addition, it also shows the influence of different values of λA on interfaces up to λ3.

The authors would like to acknowledge the use of the computational facilities provided by the University of Warwick Scientific Computing Research Technology Platform. We would also like to acknowledge the EPSRC Centre for Doctoral Training in Modeling of Heterogeneous Systems (EPSRC Grant No. EP/S022848/1). D.Q. is funded by EPSRC Program Grant No. EP/R018820/1. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission.

The authors have no conflicts to disclose.

Katarina E. Blow: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Gareth A. Tribello: Methodology (supporting); Writing – review & editing (equal). Gabriele C. Sosso: Conceptualization (equal); Formal analysis (supporting); Methodology (equal); Supervision (equal); Writing – review & editing (equal). David Quigley: Conceptualization (equal); Formal analysis (supporting); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data generated by the simulations reported in this manuscript are openly available at http://wrap.warwick.ac.uk/175784.