Forward Flux Sampling Simulations Research Articles

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.

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.

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.

Entropy-reduced Retention Times in Magnetic Memory Elements: A Case of the Meyer-Neldel Compensation Rule.

We compute mean waiting times between thermally activated magnetization reversals in a nanodisk with parameters similar to a free CoFeB layer used in magnetic random access memories. By combining Langer's theory and forward flux sampling simulations, we show that the Arrhenius prefactor can take values up to 10^{21} Hz, orders of magnitude beyond the value of 10^{9} Hz typically assumed, and varies drastically as a function of material parameters. We show that the prefactor behaves like an exponential of the activation energy, which highlights a case of the Meyer-Neldel compensation rule. This suggests that modeling information retention times with a barrier-independent prefactor in such magnetic storage elements is not justified.

An efficient method to reconstruct free energy profiles for diffusive processes in transition interface sampling and forward flux sampling simulations.

We propose a Reweighted Partial Path (RPP) method to compute free energy profiles for diffusive processes in single Transition Interface Sampling (TIS) or Forward Flux Sampling (FFS) simulations. The method employs a partial path reweighting strategy, based on the memory loss assumption for diffusive systems, to derive the equilibrium distribution of states along a chosen order parameter from TIS or FFS trajectories. No additional calculations such as reverse TIS or umbrella sampling are required. The application of the RPP method is demonstrated by calculating the nucleation free energy of early-stage Cu precipitates in a dilute Fe-Cu alloy.

Heterogeneous seeded molecular dynamics as a tool to probe the ice nucleating ability of crystalline surfaces.

Ice nucleation plays a significant role in a large number of natural and technological processes, but it is challenging to investigate experimentally because of the small time scales (ns) and short length scales (nm) involved. On the other hand, conventional molecular simulations struggle to cope with the relatively long time scale required for critical ice nuclei to form. One way to tackle this issue is to take advantage of free energy or path sampling techniques. Unfortunately, these are computationally costly. Seeded molecular dynamics is a much less demanding alternative that has been successfully applied already to study the homogeneous freezing of water. However, in the case of heterogeneous ice nucleation, nature's favourite route to form ice, an array of suitable interfaces between the ice seeds and the substrate of interest has to be built, and this is no trivial task. In this paper, we present a Heterogeneous SEEDing (HSEED) approach which harnesses a random structure search framework to tackle the ice-substrate challenge, thus enabling seeded molecular dynamics simulations of heterogeneous ice nucleation on crystalline surfaces. We validate the HSEED framework by investigating the nucleation of ice on (i) model crystalline surfaces, using the coarse-grained mW model, and (ii) cholesterol crystals, employing the fully atomistic TIP4P/ice water model. We show that the HSEED technique yields results in excellent agreement with both metadynamics and forward flux sampling simulations. Because of its computational efficiency, the HSEED method allows one to rapidly assess the ice nucleation ability of whole libraries of crystalline substrates-a long-awaited computational development in, e.g., atmospheric science.

Arm retraction dynamics of entangled star polymers: A forward flux sampling method study

The study of dynamics and rheology of well-entangled branched polymers remains a challenge for computer simulations due to the exponentially growing terminal relaxation times of these polymers with increasing molecular weights. We present an efficient simulation algorithm for studying the arm retraction dynamics of entangled star polymers by combining the coarse-grained slip-spring (SS) model with the forward flux sampling (FFS) method. This algorithm is first applied to simulate symmetric star polymers in the absence of constraint release (CR). The reaction coordinate for the FFS method is determined by finding good agreement of the simulation results on the terminal relaxation times of mildly entangled stars with those obtained from direct shooting SS model simulations with the relative difference between them less than 5%. The FFS simulations are then carried out for strongly entangled stars with arm lengths up to 16 entanglements that are far beyond the accessibility of brute force simulations in the non-CR condition. Apart from the terminal relaxation times, the same method can also be applied to generate the relaxation spectra of all entanglements along the arms which are desired for the development of quantitative theories of entangled branched polymers. Furthermore, we propose a numerical route to construct the experimentally measurable relaxation correlation functions by effectively linking the data stored at each interface during the FFS runs. The obtained star arm end-to-end vector relaxation functions Φ(t) and the stress relaxation function G(t) are found to be in reasonably good agreement with standard SS simulation results in the terminal regime. Finally, we demonstrate that this simulation method can be conveniently extended to study the arm-retraction problem in entangled star polymer melts with CR by modifying the definition of the reaction coordinate, while the computational efficiency will depend on the particular slip-spring or slip-link model employed.

Primary Fibril Nucleation of Aggregation Prone Tau Fragments PHF6 and PHF6.

We performed replica exchange molecular dynamics and forward flux sampling simulations of hexapeptide VQIINK and VQIVYK systems, also known as, respectively, fragments PHF6* and PHF6 from the tau protein. Being a part of the microtubule binding region, these fragments are known to be aggregation prone, and at least one of them is a prerequisite for fibril formation of the tau protein. Using a coarse-grained force field, we establish the phase behavior of both fragments, and investigate the nucleation kinetics for the conversion into a β-sheet fibril. As the conversion is, in principle, a reversible process, we predict the rate constants for both the fibril formation and melting, and examine the corresponding mechanisms. Our simulations indicate that, while both fragments form disordered aggregates, only PHF6 is able to form β-sheet fibrils. This observation provides a possible explanation for the lack of available steric zipper crystal structures for PHF6*.

Large deviations of Rouse polymer chain: First passage problem.

The purpose of this paper is to investigate several analytical methods of solving first passage (FP) problem for the Rouse model, a simplest model of a polymer chain. We show that this problem has to be treated as a multi-dimensional Kramers' problem, which presents rich and unexpected behavior. We first perform direct and forward-flux sampling (FFS) simulations and measure the mean first-passage time τ(z) for the free end to reach a certain distance z away from the origin. The results show that the mean FP time is getting faster if the Rouse chain is represented by more beads. Two scaling regimes of τ(z) are observed, with transition between them varying as a function of chain length. We use these simulation results to test two theoretical approaches. One is a well known asymptotic theory valid in the limit of zero temperature. We show that this limit corresponds to fully extended chain when each chain segment is stretched, which is not particularly realistic. A new theory based on the well known Freidlin-Wentzell theory is proposed, where dynamics is projected onto the minimal action path. The new theory predicts both scaling regimes correctly, but fails to get the correct numerical prefactor in the first regime. Combining our theory with the FFS simulations leads us to a simple analytical expression valid for all extensions and chain lengths. One of the applications of polymer FP problem occurs in the context of branched polymer rheology. In this paper, we consider the arm-retraction mechanism in the tube model, which maps exactly on the model we have solved. The results are compared to the Milner-McLeish theory without constraint release, which is found to overestimate FP time by a factor of 10 or more.

Simulation of protein association: Kinetic pathways towards crystal contacts

We conducted molecular dynamics simulations combined with distance-based umbrella sampling and forward flux sampling to investigate the early stages of protein crystallization. Formation of contacts with long-range interactions and/or an exposed position on the protein surface was kinetically preferred over more stable hydrophobic contacts with a shorter attractive range, while the thermodynamic stability of the protein crystal was provided by hydrophobic interactions. Contacts with a large interaction area showed complex dissociation pathways that were not detected by distance-based umbrella sampling. Instead, forward flux sampling simulations of contact dissociation identified long-range attractive interactions.

Automatic, optimized interface placement in forward flux sampling simulations

Forward flux sampling (FFS) provides a convenient and efficient way to simulate rare events in equilibrium or non-equilibrium systems. FFS ratchets the system from an initial state to a final state via a series of interfaces in phase space. The efficiency of FFS depends sensitively on the positions of the interfaces. We present two alternative methods for placing interfaces automatically and adaptively in their optimal locations, on-the-fly as an FFS simulation progresses, without prior knowledge or user intervention. These methods allow the FFS simulation to advance efficiently through bottlenecks in phase space by placing more interfaces where the probability of advancement is lower. The methods are demonstrated both for a single-particle test problem and for the crystallization of Yukawa particles. By removing the need for manual interface placement, our methods both facilitate the setting up of FFS simulations and improve their performance, especially for rare events which involve complex trajectories through phase space, with many bottlenecks.

Optimizing transition interface sampling simulations

We demonstrate that a recently proposed adaptive optimization algorithm for forward flux sampling simulations [E. E. Borrero and F. A. Escobedo, J. Chem. Phys. 129, 024115 (2008)] can be easily applied within the framework of transition interface sampling. This optimization algorithm systematically identifies the kinetic bottlenecks along the order parameter used to partition phase space via interfaces and improves the statistical accuracy of the reaction rate constant estimate. In different versions of the algorithm, the number or the placement of the interfaces (or both) are varied in order to allocate the numerical effort in a balanced way. The algorithm is demonstrated for a simple two-dimensional model and for the dipole flip transition of icelike structures inside carbon nanotubes. For these test systems, the optimization yielded an efficiency increase by a factor of 2-15.

Kinetics and Reaction Coordinates of the Reassembly of Protein Fragments Via Forward Flux Sampling

We studied the mechanism of the reassembly and folding process of two fragments of a split lattice protein by using forward flux sampling (FFS). Our results confirmed previous thermodynamics and kinetics analyses that suggested that the disruption of the critical core (of an unsplit protein that folds by a nucleation mechanism) plays a key role in the reassembly mechanism of the split system. For several split systems derived from a parent 48-mer model, we estimated the reaction coordinates in terms of collective variables by using the FFS least-square estimation method and found that the reassembly transition is best described by a combination of the total number of native contacts, the number of interchain native contacts, and the total conformational energy of the split system. We also analyzed the transition path ensemble obtained from FFS simulations using the estimated reaction coordinates as order parameters to identify the microscopic features that differentiate the reassembly of the different split systems studied. We found that in the fastest folding split system, a balanced distribution of the original-core amino acids (of the unsplit system) between protein fragments propitiates interchain interactions at early stages of the folding process. Only this system exhibits a different reassembly mechanism from that of the unsplit protein, involving the formation of a different folding nucleus. In the slowest folding system, the concentration of the folding nucleus in one fragment causes its early prefolding, whereas the second fragment tends to remain as a detached random coil. We also show that the reassembly rate can be either increased or decreased by tuning interchain cooperativeness via the introduction of a single point mutation that either strengthens or weakens one of the native interchain contacts (prevalent in the transition state ensemble).

DNA looping provides stability and robustness to the bacteriophage λ switch

The bistable gene regulatory switch controlling the transition from lysogeny to lysis in bacteriophage lambda presents a unique challenge to quantitative modeling. Despite extensive characterization of this regulatory network, the origin of the extreme stability of the lysogenic state remains unclear. We have constructed a stochastic model for this switch. Using Forward Flux Sampling simulations, we show that this model predicts an extremely low rate of spontaneous prophage induction in a recA mutant, in agreement with experimental observations. In our model, the DNA loop formed by octamerization of CI bound to the O(L) and O(R) operator regions is crucial for stability, allowing the lysogenic state to remain stable even when a large fraction of the total CI is depleted by nonspecific binding to genomic DNA. DNA looping also ensures that the switch is robust to mutations in the order of the O(R) binding sites. Our results suggest that DNA looping can provide a mechanism to maintain a stable lysogenic state in the face of a range of challenges including noisy gene expression, nonspecific DNA binding, and operator site mutations.