Forward Flux Sampling Research Articles

Bridging classical nucleation theory and molecular dynamics simulation for homogeneous ice nucleation.

Water freezing, initiated by ice nucleation, occurs widely in nature, ranging from cellular to global phenomena. Ice nucleation has been experimentally proven to require the formation of a critical ice nucleus, consistent with classical nucleation theory (CNT). However, the accuracy of CNT quantitative predictions of critical cluster sizes and nucleation rates has never been verified experimentally. In this study, we circumvent this difficulty by using molecular dynamics (MD) simulation. The physical properties of water/ice for CNT predictions, including density, chemical potential difference, and diffusion coefficient, are independently obtained using MD simulation, whereas the calculation of interfacial free energy is based on thermodynamic assumptions of CNT, including capillarity approximation among others. The CNT predictions are compared to the MD evaluations of brute-force simulations and forward flux sampling methods. We find that the CNT and MD predicted critical cluster sizes are consistent, and the CNT predicted nucleation rates are higher than the MD predicted values within three orders of magnitude. We also find that the ice crystallized from supercooled water is stacking-disordered ice with a stacking of cubic and hexagonal ices in four representative types of stacking. The prediction discrepancies in nucleation rate mainly arise from the stacking-disordered ice structure, the asphericity of ice cluster, the uncertainty of ice-water interfacial free energy, and the kinetic attachment rate. Our study establishes a relation between CNT and MD to predict homogeneous ice nucleation.

Comparing brute force to transition path sampling for gas hydrate nucleation with a flat interface: comments on time reversal symmetry.

Fluid to solid nucleation is often investigated with the rare event method transition path sampling (TPS). I claim that the inherent irreversibility of solid nucleation, even at stationary conditions, calls into question TPS's applicability for determining solid nucleation mechanisms, especially for pre-critical behavior. Even when applied to a phenomenon which displays time reversal asymmetry like solid nucleation, TPS is a good means of exploring phase space and giving trends in post-critical structure, and its ability to facilitate nucleation rate and free energy calculations remains outstanding. Forward-only splitting and ratcheting methods such as forward flux sampling are more attractive for understanding nucleation mechanisms as they do not require time reversal symmetry, but at low driving forces may suffer from the same limitations as brute force: they may never make it to the first ratchet. Here I briefly summarize the TPS method and gas hydrate nucleation simulation literature, focusing on topics within both to facilitate a comparison of brute force hydrate nucleation to transition path sampling of hydrate nucleation. Perhaps anecdotally, the brute force technique results in more crystalline trajectories despite having higher driving forces than TPS. I maintain this difference is because of the inherent irreversibility of hydrate nucleation, meaning its pre-critical behavior cannot accurately be determined by the melting trajectories that comprise approximately half of the configurations in TPS's path ensemble.

Multivariate alarm systems to recognize rare unpostulated abnormal events

AbstractMost chemical and manufacturing plants have safety/reliability systems in place that are well equipped to handle commonly occurring postulated abnormal events, but often prove to be ineffective in predicting highly unanticipated and randomly occurring unpostulated abnormal events. In this paper, an advisory system is developed to analyze and recognize such events, consisting of novel, multivariate alarm systems and response actions introduced using process modeling and path‐sampling for unpostulated abnormal events. It augments existing safety/reliability systems, suggesting actions when unanticipated abnormal events are approached. Forward‐flux sampling (FFS), developed to discover rare molecular dynamics pathways, is applied. With an approximate process model, for an exothermic continuous‐stirred tank reactor using a perturbed feed concentration, the FFS algorithm is applied to identify rare trajectories between high‐conversion and low‐conversion steady‐states, with key process variables saved at various “crossing points.” Then, committer probabilities, , are computed at each crossing point, yielding a mathematical model that expresses as a function of the key process variables; that is, the reactor temperature (T), cooling‐water flow rate (), and cooling‐water temperature (), selected as choices for the primary alarm variables. For these, alarm thresholds; that is, L—low, LL—low low, and LLL—low low low, are suggested by computing their critical ranges, given the ranges for every alarm threshold.

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

Forward flux sampling (FFS) is a path sampling technique widely used in computer simulations of crystal nucleation from the melt. In such studies, the order parameter underpinning the progress of the FFS algorithm is often the size of the largest crystalline nucleus. In this work, we investigate the effects of two computational aspects of FFS simulations, using the prototypical Lennard-Jones liquid as our computational test bed. First, we quantify the impact of the positioning of the liquid basin and first interface in the space of the order parameter. In particular, we demonstrate that these choices are key to ensuring the consistency of the FFS results. Second, we focus on the frequently encountered scenario where the population of crystalline nuclei is such that there are multiple clusters of size comparable to the largest one. We demonstrate the contribution of clusters other than the largest cluster to the initial flux; however, we also show that they can be safely ignored for the purposes of converging a full FFS calculation. We also investigate the impact of different clusters merging, a process that appears to be facilitated by substantial spatial correlations-at least at the supercooling considered here. Importantly, all of our results have been obtained as a function of system size, thus contributing to the ongoing discussion on the impact of finite size effects on simulations of crystal nucleation. Overall, this work either provides or justifies several practical guidelines for performing FFS simulations that can also be applied to more complex and/or computationally expensive models.

On the challenge of sampling multiple nucleation pathways: A case study of heterogeneous ice nucleation on FCC (211) surface.

Heterogeneous ice nucleation is ubiquitous but its microscopic mechanisms can be extraordinarily complex even on a simple surface. Such complexity poses a challenge in modeling nucleation using advanced sampling methods. Here, we investigate heterogeneous ice nucleation on an FCC (211) surface by a forward flux sampling (FFS) method to understand how the complexity in nucleation pathways affects the accuracy of FFS. We first show the commonly adopted, size-based order parameter fails to describe heterogeneous ice nucleation on the FCC (211) surface. Inclusion of geometric anisotropy of ice nucleus as an additional descriptor is found to significantly improve the quality of the size-based order parameter for the current system. Subsequent application of this new order parameter in FFS identifies two competing ice nucleation pathways in the system: a primary-prism-planed (PPP) path and a secondary-prism-planed (SPP) path, both leading to the formation of hexagonal ice but with different crystalline orientations. Although the PPP pathway dominates ice nucleation on the FCC (211) surface, the occurrence of the less efficient SPP pathway, particularly its strong presence at the early stage of FFS, is found to yield a significant statistical uncertainty in the calculated FFS rate constant. We develop a two-path model that enables gaining a general, quantitative understanding of the impact of initial finite sampling on the reliability of FFS calculations in the presence of multiple nucleation pathways. Our study also suggests a few general strategies for improving the accuracy of FFS when exploring unknown but complex systems.

Nonadiabatic Forward Flux Sampling for Excited-State Rare Events.

We present a rare event sampling scheme applicable to coupled electronic excited states. In particular, we extend the forward flux sampling (FFS) method for rare event sampling to a nonadiabatic version (NAFFS) that uses the trajectory surface hopping (TSH) method for nonadiabatic dynamics. NAFFS is applied to two dynamically relevant excited-state models that feature an avoided crossing and a conical intersection with tunable parameters. We investigate how nonadiabatic couplings, temperature, and reaction barriers affect transition rate constants in regimes that cannot be otherwise obtained with plain, traditional TSH. The comparison with reference brute-force TSH simulations for limiting cases of rareness shows that NAFFS can be several orders of magnitude cheaper than conventional TSH and thus represents a conceptually novel tool to extend excited-state dynamics to time scales that are able to capture rare nonadiabatic events.

Classical nucleation theory and molecular dynamics simulation: Cavitation

Based on the results of molecular dynamics simulations of spontaneous cavitation in a stretched Lennard-Jones liquid at temperatures close to the triple point temperature, the possibilities and the limits of applicability of approximations of the classical nucleation theory were analyzed. In systems containing from 5×103 to 3×105 particles, the nucleation rate, radii and work of formation of critical bubbles, and the diffusion coefficient of nuclei were calculated. Methods such as direct forward flux sampling, seeding, and well-tempered metadynamics were used. From the obtained data, the surface free energy of critical bubbles was determined, and their sizes and shapes were estimated. The asymptotic laws of the behavior of nucleation parameters when approaching the spinodal of a stretched liquid are discussed.

Kinetic control of competing nuclei in a dimer lattice-gas model.

Nucleation is a key step in the synthesis of a new material from a solution. The well-established lattice-gas models can be used to gain insight into the basic physics of nucleation pathways involving a single nucleus type. In many situations, a solution is supersaturated with respect to more than one precipitating phase. This can generate a population of both stable and metastable nuclei on similar timescales and, hence, complex nucleation pathways involving a competition between the two. In this study, we introduce a lattice-gas model based on two types of interacting dimers representing the particles in a solution. Each type of dimer nucleates to a specific space-filling structure. Our model is tuned such that stable and metastable phases nucleate on a similar timescale. Either structure may nucleate first, with a probability sensitive to the relative rate at which a solute is replenished from their respective reservoirs. We calculate these nucleation rates via forward flux sampling and demonstrate how the resulting data can be used to infer the nucleation outcome and pathway. Possibilities include direct nucleation of the stable phase, domination of long-lived metastable crystallites, and pathways in which the stable phase nucleates only after multiple post-critical nuclei of the metastable phase have appeared.

Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory.

We calculate bubble nucleation rates in a Lennard-Jones fluid through explicit molecular dynamics simulations. Our approach-based on a recent free energy method (dubbed reweighted Jarzynski sampling), transition state theory, and a simple recrossing correction-allows us to probe a fairly wide range of rates in several superheated and cavitation regimes in a consistent manner. Rate predictions from this approach bridge disparate independent literature studies on the same model system. As such, we find that rate predictions based on classical nucleation theory, direct brute force molecular dynamics simulations, and seeding are consistent with our approach and one another. Published rates derived from forward flux sampling simulations are, however, found to be outliers. This study serves two purposes: First, we validate the reliability of common modeling techniques and extrapolation approaches on a paradigmatic problem in materials science and chemical physics. Second, we further test our highly generic recipe for rate calculations, and establish its applicability to nucleation processes.

Ideal Conductor Model: An Analytical Finite-Size Correction for Nonequilibrium Molecular Dynamics Simulations of Ion Transport through Nanoporous Membranes.

Modulating ion transport through nanoporous membranes is critical to many important chemical and biological separation processes. The corresponding transport timescales, however, are often too long to capture accurately using conventional molecular dynamics (MD). Recently, path sampling techniques, such as forward-flux sampling (FFS), have emerged as attractive alternatives for efficiently and accurately estimating arbitrarily long ionic passage times. Here, we use non-equilibrium MD and FFS to explore how the kinetics and mechanisms of pressure-driven chloride transport through a nanoporous graphitic membrane are affected by its lateral dimensions. We not only find ionic passage times and free energy barriers to decrease dramatically upon increasing the membrane surface area but also observe an abrupt and discontinuous change in the locus of the transition state. These strong finite size effects arise due to the cumulative effect of the periodic images of the leading ion entering the pore on the distribution of the induced excess charge at the membrane surface in the feed. By assuming that the feed is an ideal conductor, we analytically derive a finite size correction term that can be computed from the information obtained from a single simulation and successfully use it to obtain corrected free energy profiles with no dependence on the system size. We then estimate ionic passage times in the thermodynamic limit by assuming an Eyring-type dependence of rates on barriers with a size-independent prefactor. This approach constitutes a universal framework for removing finite size artifacts in molecular simulations of ion transport through nanoporous membranes and biological channel proteins.

Practical guide to replica exchange transition interface sampling and forward flux sampling.

Path sampling approaches have become invaluable tools to explore the mechanisms and dynamics of the so-called rare events that are characterized by transitions between metastable states separated by sizable free energy barriers. Their practical application, in particular to ever more complex molecular systems, is, however, not entirely trivial. Focusing on replica exchange transition interface sampling (RETIS) and forward flux sampling (FFS), we discuss a range of analysis tools that can be used to assess the quality and convergence of such simulations, which is crucial to obtain reliable results. The basic ideas of a step-wise evaluation are exemplified for the study of nucleation in several systems with different complexities, providing a general guide for the critical assessment of RETIS and FFS simulations.

Molecular dynamics simulation-guided toehold mediated strand displacement probe for single-nucleotide variants detection.

Single nucleotide variant (SNV) has become an emerging biomarker for various diseases such as cancers and infectious diseases. Toehold-mediated strand displacement (TMSD), the core reaction of DNA nanotechnology, has been widely leveraged to identify SNVs. However, inappropriate choice of mismatch location results in poor discrimination ability. Here, we comprehensively investigate the effect of mismatch location on TMSD kinetics by molecular dynamic simulation tool oxDNA through umbrella sampling and forward flux sampling disclosing that mismatches at the border of the toehold and branch migration domain yield the lowest TMSD reaction rate. Nine disease-related SNVs (SARS-CoV-2-D614G, EGFR-L858R, EGFR-T790M, KRAS-G12R, etc.) were tested experimentally showing a good agreement with simulation. The best choice of mismatch location enables high discrimination factor with a median of 124 for SNV and wild type. Coupling with a probe-sink system, a low variant allele frequency of 0.1% was detected with 3 S/N. We successfully used the probes to detect SNVs with high confidence in the PCR clones of constructed plasmids. This work provides mechanistic insights into TMSD process at the single-nucleotide level and can be a guidance for the design of TMSD system with fine-tuning kinetics for various applications in biosensors and nanotechnology.

Homogeneous ice nucleation rate at negative pressures: The role of the density anomaly

Recent experiments suggest a role for pressure fluctuations in nucleation. Homogeneous ice nucleation rates for the ML-mW and mW water models are evaluated at pressures ranging from atmospheric to -1000 atm, using forward flux sampling and constant cooling simulations. Results indicate that the density difference Δνls between water and ice exhibited by these models plays a central role in controlling the change in nucleation rate with pressure. A linear function is found to be a reasonable approximation for lines of constant nucleation rate, which can be useful in making experimental predictions to advance the study of ice nucleation mechanisms.

Single-stage direct Langevin dynamic simulations of transitions over arbitrarily high energy barriers: Concept of energy-dependent temperature

In this paper we present an algorithm which allows single-stage direct Langevin dynamics simulations of transitions over arbitrary high energy barriers employing the concept of the energy-dependent temperature (EDT). In our algorithm, simulation time required for the computation of the corresponding switching rate does not increase with energy barrier. This is achieved by using in simulations an effective temperature which depends on the system energy: around the energy minima this temperature is high and tends towards the room temperature when the energy approaches the saddle point value. Switching times computed via our EDT algorithm show an excellent agreement with results obtained with the established forward flux sampling (FFS) method. As the simulation time required by our method does not increase with the energy barrier, we achieve a very large speedup when compared even to the highly optimized FFS version. In addition, our method does not suffer from stability problems occurring in multi-stage algorithms (like FFS and 'energy bounce' methods) due to the multiplication of a large number of transition probabilities between the interfaces.

Concept of the Energy-Dependent Temperature for Direct Langevin Dynamics Simulations of Rare Events in Systems with Arbitrary High Energy Barriers.

We suggest an algorithm which allows single-stage direct Langevin dynamics simulations of transitions over arbitrarily high energy barriers. For this purpose, we propose a concept of the energy-dependent temperature (EDT): near the energy minima this temperature is high, but it tends toward room temperature for energies approaching the barrier value. In the resulting algorithm simulation time required for the computation of the escape rate over the barrier does not increase with barrier height. Switching times computed via our EDT algorithm agree very well with those obtained with the forward flux sampling (FFS). As the simulation time required by our method does not increase with the energy barrier, we achieve a very large speed-up compared even to our strongly optimized version of FFS (and all other multistage algorithms). In addition, our approach is free from the instability occurring in all multistage "climbing" methods where a product of a large number of transition probabilities between the interfaces must be computed.

The interplay of diffusion and heterogeneity in nucleation of the networked Ising model

Nucleation underlies the dynamics of most first-order phase transitions in natural and man-made systems. However, most of the systems of interest are out of equilibrium. Little is known on the effect of nonequilibrium factors on the dynamics of nucleation. Here, we use the forward flux sampling method to investigate the effect of nonequilibrium diffusion on nucleation in small-world Ising networks, wherein spins can be exchanged between nearest-neighboring nodes. We introduce a parameter α to quantify the difference of nucleation rate with and without diffusion. We find that α shows a nonmonotonic dependence on the rewiring probability p of small-world networks. In particular, for different diffusion probability D, a crossover happens at p ≃ 0.17, below which the nucleation rate decreases as D increases, suggesting that the diffusion is against nucleation; while above which the nucleation rate increases with D, indicating that the diffusion is in favor of nucleation. By identifying the distinct features of nucleating clusters along the pathways for different randomness of networks, we reveal the underlying mechanism of such a nontrivial dependence.

Understanding rare safety and reliability events using forward-flux sampling

This paper presents improvements upon methods that explore process trajectories leading to rare safety and reliability events. It applies the forward-flux sampling (FFS) method from the family of sampling algorithms developed to discover rare molecular dynamics pathways. For a relatively simple, dynamic exothermic CSTR model with noisy feed concentration, it shows how to apply the FFS algorithm to simulate and analyze rare trajectories between high- and low-conversion steady states. First, it compares results with a less efficient brute-force (BF) method, and then with a transition-path sampling (TPS) method. The effects of varying key process parameters; i.e., the residence time, τ, noise variance, ση2, and the controller gain, KC, which impact the rareness of an event, are investigated. Rates of rare-path transition between high- and low-conversion steady states, forward and backward, are shown to exhibit equilibrium ratios independent of ση2, with the forward rates decreasing with τ and KC, and the backward rates increasing with KC, whereas both increase with ση2. These unanticipated relationships should lead to improved alarm notifications being developed in this research work.

How to quantify and avoid finite size effects in computational studies of crystal nucleation: The case of heterogeneous ice nucleation.

Computational studies of crystal nucleation can be impacted by finite size effects, primarily due to unphysical interactions between crystalline nuclei and their periodic images. It is, however, not always feasible to systematically investigate the sensitivity of nucleation kinetics and mechanism to system size due to large computational costs of nucleation studies. Here, we use jumpy forward flux sampling to accurately compute the rates of heterogeneous ice nucleation in the vicinity of square-shaped model structureless ice nucleating particles (INPs) of different sizes and identify three distinct regimes for the dependence of rate on the INP dimension, L. For small INPs, the rate is a strong function of L due to the artificial spanning of critical nuclei across the periodic boundary. Intermediate-sized INPs, however, give rise to the emergence of non-spanning "proximal" nuclei that are close enough to their periodic images to fully structure the intermediary liquid. While such proximity can facilitate nucleation, its effect is offset by the higher density of the intermediary liquid, leading to artificially small nucleation rates overall. The critical nuclei formed at large INPs are neither spanning nor proximal. Yet, the rate is a weak function of L, with its logarithm scaling linearly with 1/L. The key heuristic emerging from these observations is that finite size effects will be minimal if critical nuclei are neither spanning nor proximal and if the intermediary liquid has a region that is structurally indistinguishable from the supercooled liquid under the same conditions.

Nucleation rate in the two dimensional Ising model in the presence of random impurities.

Nucleation phenomena are ubiquitous in nature and the presence of impurities in every real and experimental system is unavoidable. Yet numerical studies of nucleation are nearly always conducted for entirely pure systems. We have studied the behaviour of the droplet free energy in two dimensional Ising model in the presence of randomly positioned static and dynamic impurities. We have shown that both the free energy barrier height and critical nucleus size monotonically decreases with increasing the impurity density for the static case. We have compared the nucleation rates obtained from the Classical Nucleation Theory and the Forward Flux Sampling method for different densities of the static impurities. The results show good agreement. In the case of dynamic impurities, we observe preferential occupancy of the impurities at the boundary positions of the nucleus when the temperature is low. This further boosts enhancement of the nucleation rate due to lowering of the effective interfacial free energy.

Impact of Surface Roughness on Crystal Nucleation

We examine the effect of rough surfaces on crystal nucleation by means of kinetic Monte Carlo simulations. Our work makes use of three-dimensional kMC models, explicit representation of transport in solution and rough surfaces modeled as randomly varying height fluctuations (roughness) with exponentially decaying correlation length (topology). We use Forward-Flux Sampling to determine the nucleation rate for crystallization for surfaces of different roughness and topology and show that the effect on crystallization is a complex interplay between the two. For surfaces with low roughness, small clusters form on the surface but as clusters become larger they are increasingly likely to be found in the bulk solution while rougher surfaces eventually favor heterogeneous nucleation on the surface. In both cases, the rough surface raises the local supersaturation in the solution thus leading to another mechanism of enhanced nucleation rate.