Optimizing the sampling and staging for simulations of rare events via forward flux sampling schemes

In this work, we present an adaptive algorithm to optimize the phase space sampling for simulations of rare events in complex systems via forward flux sampling (FFS) schemes. In FFS, interfaces are used to partition the phase space along an order parameter lambda connecting the initial and final regions of interest. Since the kinetic "bottleneck" regions along the order parameter are not usually known beforehand, an adaptive procedure is used that first finds these regions by estimating the rate constants associated with reaching subsequent interfaces; thereafter, the FFS simulation is reset to concentrate the sampling on those bottlenecks. The approach can optimize for either the number and position of the interfaces (i.e., optimized lambda phase staging) or the number M of fired trial runs per interface (i.e., the {M(i)} set) to minimize the statistical error in the rate constant estimation per simulation period. For example, the optimization of the lambda staging leads to a net constant flux of partial trajectories between interfaces and hence a constant flux of connected paths throughout the region between the two end states. The method is demonstrated for several test systems, including the folding of a lattice protein. It is shown that the proposed approach leads to an optimized lambda staging and {M(i)} set which increase the computational efficiency of the sampling algorithm.     

Non-stationary forward flux sampling

We present a method, Non-Stationary Forward Flux Sampling, that allows efficient simulation of rare events in both stationary and non-stationary stochastic systems. The method uses stochastic branching and pruning to achieve uniform sampling of trajectories in phase space and time, leading to accurate estimates for time-dependent switching propensities and time-dependent phase space probability densities. It is suitable for equilibrium or non-equilibrium systems, in or out of stationary state, including non-Markovian or externally driven systems. We demonstrate the validity of the technique by applying it to a one-dimensional barrier crossing problem that can be solved exactly, and show its usefulness by applying it to the time-dependent switching of a genetic toggle switch.     

Adaptive importance sampling Monte Carlo simulation of rare transition events

We develop a general theoretical framework for the recently proposed importance sampling method for enhancing the efficiency of rare-event simulations [W. Cai, M. H. Kalos, M. de Koning, and V. V. Bulatov, Phys. Rev. E 66, 046703 (2002)], and discuss practical aspects of its application. We define the success/fail ensemble of all possible successful and failed transition paths of any duration and demonstrate that in this formulation the rare-event problem can be interpreted as a "hit-or-miss" Monte Carlo quadrature calculation of a path integral. The fact that the integrand contributes significantly only for a very tiny fraction of all possible paths then naturally leads to a "standard" importance sampling approach to Monte Carlo (MC) quadrature and the existence of an optimal importance function. In addition to showing that the approach is general and expected to be applicable beyond the realm of Markovian path simulations, for which the method was originally proposed, the formulation reveals a conceptual analogy with the variational MC (VMC) method. The search for the optimal importance function in the former is analogous to finding the ground-state wave function in the latter. In two model problems we discuss practical aspects of finding a suitable approximation for the optimal importance function. For this purpose we follow the strategy that is typically adopted in VMC calculations: the selection of a trial functional form for the optimal importance function, followed by the optimization of its adjustable parameters. The latter is accomplished by means of an adaptive optimization procedure based on a combination of steepest-descent and genetic algorithms.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

Folding kinetics of a lattice protein via a forward flux sampling approach

We implement a forward flux sampling approach [R. J. Allen et al., J. Chem. Phys. 124, 194111 (2006)] for calculating transition rate constants and for sampling paths of protein folding events. The algorithm generates trajectories for the transition between the unfolded and folded states as chains of partially connected paths, which can be used to obtain the transition-state ensemble and the properties that characterize these intermediates. We apply this approach to Monte Carlo simulations of a model lattice protein in open space and in confined spaces of varying dimensions. We study the effect of confinement on both protein thermodynamic stability and folding kinetics; the former by mapping free-energy landscapes and the latter by the determination of rate constants and mechanistic details of the folding pathway. Our results show that, for the range of temperatures where the native state is stable, confinement of a protein destabilizes the unfolded state by reducing its entropy, resulting in increased thermodynamic stability of the folded state. Relative to the folding in open space, we find that the kinetics can be accelerated at temperatures above the temperature at which the unconfined protein folds fastest and that the rate constant increases with the number of constrained dimensions. By examining the statistical properties of the transition-state ensemble, we detect signs of a classical nucleation folding mechanism for a core of native contacts formed at an early stage of the process. This nucleus acts as folding foci and is composed of those residues that have higher probability to form native contacts in the transition-state intermediates, which can vary depending on the confinement conditions of the system.     

Overcoming barriers in trajectory space: mechanism and kinetics of rare events via Wang-Landau enhanced transition path sampling

Within the framework of transition path sampling (TPS), activation energies can be computed as path ensemble averages without a priori information about the reaction mechanism [C. Dellago and P. G. Bolhuis, Mol. Simul. 30, 795 (2004)]. Activation energies computed for different conditions can then be used to determine by numerical integration the rate constant for a system of interest from the rate constant known for a reference system. However, in systems with complex potential energy surfaces, multiple reaction pathways may exist making ergodic sampling of trajectory space difficult. Here, we present a combination of TPS with the Wang-Landau (WL) flat-histogram algorithm for an efficient sampling of the transition path ensemble. This method, denoted by WL-TPS, has the advantage that from one single simulation, activation energies at different temperatures can be determined even for systems with multiple reaction mechanisms. The proposed methodology for rate constant calculations does not require the knowledge of the reaction coordinate and is generally applicable to Arrhenius and non-Arrhenius processes. We illustrate the applicability of this technique by studying a two-dimensional toy system consisting of a triatomic molecule immersed in a fluid of repulsive soft disks. We also provide an expression for the calculation of activation volumes from path averages such that the pressure dependence of the rate constant can be obtained by numerical integration.     

Computing stationary distributions in equilibrium and nonequilibrium systems with forward flux sampling

We present a method for computing stationary distributions for activated processes in equilibrium and nonequilibrium systems using forward flux sampling. In this method, the stationary distributions are obtained directly from the rate constant calculations for the forward and backward reactions; there is no need to perform separate calculations for the stationary distribution and the rate constant. We apply the method to the nonequilibrium rare event problem proposed by Maier and Stein, to nucleation in a 2-dimensional Ising system, and to the flipping of a genetic switch.     

Efficient transition path sampling for systems with multiple reaction pathways

A new algorithm is developed for sampling transition paths and computing reaction rates. To illustrate the use of this method, we study a two-dimensional system that has two reaction pathways: one pathway is straight with a relatively high barrier and the other is roundabout with a lower barrier. The transition rate and the ratio between the numbers of the straight and roundabout transition paths are computed for a wide range of temperatures. Our study shows that the harmonic approximation for fluctuations about the steepest-descent paths is not valid even at relatively low temperatures and, furthermore, that factors related to entropy have to be determined by the global geometry of the potential-energy surface (rather than just the local curvatures alone) for complex reaction systems. It is reasonable to expect that this algorithm is also applicable to higher dimensional systems.     

Forward flux sampling-type schemes for simulating rare events: efficiency analysis

We analyze the efficiency of several simulation methods which we have recently proposed for calculating rate constants for rare events in stochastic dynamical systems in or out of equilibrium. We derive analytical expressions for the computational cost of using these methods and for the statistical error in the final estimate of the rate constant for a given computational cost. These expressions can be used to determine which method to use for a given problem, to optimize the choice of parameters, and to evaluate the significance of the results obtained. We apply the expressions to the two-dimensional nonequilibrium rare event problem proposed by Maier and Stein [Phys. Rev. E 48, 931 (1993)]. For this problem, our analysis gives accurate quantitative predictions for the computational efficiency of the three methods.     

An efficient self-optimized sampling method for rare events in nonequilibrium systems

Rare events such as nucleation processes are of ubiquitous importance in real systems.The most popular method for nonequilibrium systems,forward flux sampling(FFS),samples rare events by using interfaces to partition the whole transition process into sequence of steps along an order parameter connecting the initial and final states.FFS usually suffers from two main difficulties:low computational efficiency due to bad interface locations and even being not applicable when trapping into unknown intermediate metastable states.In the present work,we propose an approach to overcome these difficulties,by self-adaptively locating the interfaces on the fly in an optimized manner.Contrary to the conventional FFS which set the interfaces with equal distance of the order parameter,our approach determines the interfaces with equal transition probability which is shown to satisfy the optimization condition.This is done by firstly running long local trajectories starting from the current interface i to get the conditional probability distribution Pc(>i|i),and then determining i+1by equaling Pc(i+1|i)to a give value p0.With these optimized interfaces,FFS can be run in a much more efficient way.In addition,our approach can conveniently find the intermediate metastable states by monitoring some special long trajectories that neither end at the initial state nor reach the next interface,the number of which will increase sharply from zero if such metastable states are encountered.We apply our approach to a two-state model system and a two-dimensional lattice gas Ising model.Our approach is shown to be much more efficient than the conventional FFS method without losing accuracy,and it can also well reproduce the two-step nucleation scenario of the Ising model with easy identification of the intermediate metastable state.     

Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: a comparison of simulation techniques

Over the last number of years several simulation methods have been introduced to study rare events such as nucleation. In this paper we examine the crystal nucleation rate of hard spheres using three such numerical techniques: molecular dynamics, forward flux sampling, and a Bennett-Chandler-type theory where the nucleation barrier is determined using umbrella sampling simulations. The resulting nucleation rates are compared with the experimental rates of Harland and van Megen [Phys. Rev. E 55, 3054 (1997)], Sinn et al. [Prog. Colloid Polym. Sci. 118, 266 (2001)], Sch?tzel and Ackerson [Phys. Rev. E 48, 3766 (1993)], and the predicted rates for monodisperse and 5% polydisperse hard spheres of Auer and Frenkel [Nature 409, 1020 (2001)]. When the rates are examined in units of the long-time diffusion coefficient, we find agreement between all the theoretically predicted nucleation rates, however, the experimental results display a markedly different behavior for low supersaturation. Additionally, we examined the precritical nuclei arising in the molecular dynamics, forward flux sampling, and umbrella sampling simulations. The structure of the nuclei appears independent of the simulation method, and in all cases, the nuclei contains on average significantly more face-centered-cubic ordered particles than hexagonal-close-packed ordered particles.     

Exploring reaction pathways with transition path and umbrella sampling: application to methyl maltoside

The transition path sampling (TPS) method is a powerful approach to study chemical reactions or transitional properties on complex potential energy landscapes. One of the main advantages of the method over potential of mean force methods is that reaction rates can be directly accessed without knowledge of the exact reaction coordinate. We have investigated the complementary nature of these two differing approaches, comparing transition path sampling with the weighted histogram analysis method to study a conformational change in a small model system. In this case study, the transition paths for a transition between two rotational conformers of a model disaccharide molecule, methyl beta-D-maltoside, were compared with a free energy surface constrained by the two commonly used glycosidic (phi,psi) torsional angles. The TPS method revealed a reaction channel that was not apparent from the potential of mean force method, and the suitability of phi and psi as reaction coordinates to describe the isomerization in vacuo was confirmed by examination of the transition path ensemble. Using both transition state theory and transition path sampling methods, the transition rate was estimated. We have estimated a characteristic time between transitions of approximately 160 ns for this rare isomerization event between the two conformations of the carbohydrate. We conclude that transition path sampling can extract subtle information about the dynamics not apparent from the potential of mean force method. However, in calculating the reaction rate, the transition path sampling method required 27.5 times the computational effort than was needed by the potential of mean force method.     

Reaction pathways leading to arylene ethynylene macrocycles via alkyne metathesis

Mechanistic studies on the direct formation of arylene ethynylene macrocycles via alkyne metathesis catalyzed by a molybdenum complex are reported. Gel permeation chromatography (GPC) and matrix-assisted laser desorption ionization (MALDI) mass spectrometry on the products from metathesis of monomer 1 show the initial formation of linear oligomers and large macrocycles (n > 6), followed by their transformation into the thermodynamically most stable product distribution-mainly the cyclic hexamer. Variable temperature and scrambling experiments reveal the reversibility of macrocycle formation. Nearly identical product distributions are observed from the cross metathesis of hexacycle 2 with diphenylacetylene and from the metathesis of bis(phenylethynyl) substituted monomer 4, demonstrating that macrocycle formation is thermodynamically rather than kinetically controlled. The metathesis byproduct, 3-hexyne, is shown to inhibit the catalyst. It is suggested that the relative metathesis rates of dialkylalkynes versus diarylalkynes trap the catalyst in a nonproductive manifold, rendering it unavailable for the productive metathesis of aryl alkylalkyne substrates. This finding indicates that dialkyl-substituted alkyne byproducts should be avoided (or efficiently removed) if the metatheses of aryl substrates, especially those with electron-withdrawing groups, are to proceed to high conversion.     

Selective sampling of transition paths

In this short paper, we introduce an approximate method for the quick estimate of rate constants based on a simple sampling method of reactive transition paths over high energy barriers. It makes use of the previously introduced accelerated molecular dynamics (MD) simulation method to generate initial points for trajectory shooting. The accelerated MD simulations, although with the loss of real dynamics, lead to a quick calculation of thermodynamic properties and at the same time produce an ensemble of configurations with an enhanced sampling over the phase space that is more "reactive." The forward/backward trajectory shooting as that used in the transition path sampling method is then initiated from the configurations obtained from accelerated MD simulations to generate transition paths on the original unbiased potential. This method selectively enhances sampling of successful trajectories and at the same time accelerates significantly the calculation of rate constants.     

Reaction pathways of propene pyrolysis

The gas‐phase reaction pathways in preparing pyrolytic carbon with propene pyrolysis have been investigated in detail with a total number of 110 transition states and 50 intermediates. The structure of the species was determined with density functional theory at B3PW91/6‐311G(d,p) level. The transition states and their linked intermediates were confirmed with frequency and the intrinsic reaction coordinates analyses. The elementary reactions were explored in the pathways of both direct and the radical attacking decompositions. The energy barriers and the reaction energies were determined with accurate model chemistry method at G3(MP2) level after an examination of the nondynamic electronic correlations. The heat capacities and entropies were obtained with statistical thermodynamics. The Gibbs free energies at 298.15 K for all the reaction steps were reported. Those at any temperature can be developed with classical thermodynamics by using the fitted (as a function of temperature) heat capacities. It was found that the most favorable paths are mainly in the radical attacking chain reactions. The chain was proposed with 26 reaction steps including two steps of the initialization of the chain to produce H and CH₃ radicals. For a typical temperature (1200 K) adopted in the experiments, the highest energy barriers were found in the production of C₃ to be 203.4 and 193.7 kJ/mol. The highest energy barriers for the production of C₂ and C were found 174.1 and 181.4 kJ/mol, respectively. These results are comparable with the most recent experimental observation of the apparent activation energy 201.9 ± 0.6 or 137 ± 25 kJ/mol. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

On the detection of rare,and moderately rare,nuclear events

Some of the more important developments in science and practical demands in commerce have been linked to attempts to detect rare events and rare contaminants, ranging from the early “counting” of solar neutrinos to the occurrence of dodder seeds in clover. For moderately rare events (≈5 to 50 counts) we consider limitations of the Poisson-normal approximation, together with the apparent problem of excessive false positives when a common expression is (mis-)used for detection decisions. For very rare events, rigorous approaches published more than half a century ago are applicable to such current problems as trace actinide contamination and nuclear treaty monitoring.     

Using redundant internal coordinates to optimize equilibrium geometries and transition states

A redundant internal coordinate system for optimizing molecular geometries is constructed from all bonds, all valence angles between bonded atoms, and all dihedral angles between bonded atoms. Redundancies are removed by using the generalized inverse of the G matrix; constraints can be added by using an appropriate projector. For minimizations, redundant internal coordinates provide substantial improvements in optimization efficiency over Cartesian and nonredundant internal coordinates, especially for flexible and polycyclic systems. Transition structure searches are also improved when redundant coordinates are used and when the initial steps are guided by the quadratic synchronous transit approach. © 1996 by John Wiley & Sons, Inc.     

Reaction pathways of the Simmons-Smith reaction

The cyclopropanation reaction of an alkene with a metal carbenoid has been studied by means of the B3LYP hybrid density functional method. The cyclopropanation of ethylene with a lithium carbenoid or a zinc carbenoid [Simmons-Smith (SS) reagent] goes through two competing pathways, methylene transfer and carbometalation. Both processes are fast for the lithium carbenoid, while, for the zinc carbenoid, only the former is fast enough to be experimentally feasible. The reaction of an SS reagent (ClZnCH(2)Cl) with ethylene and an allyl alcohol in the presence of ZnCl(2) was also studied. The allyl alcohol reaction was modeled with an SS reagent/alkoxide complex (ClCH(2)ZnOCH(2)CH=CH(2)) formed from the SS reagent and allyl alcohol. Two modes of acceleration were found. The first involves the well-accepted mechanism of 1,2-chlorine migration, and the second involves a five-centered bond alternation. The latter was found to be more facile than the former and to operate equally well both with ethylene and with aggregates of SS reagent/alkoxide complexes. Calculations on the SS reaction with 2-cyclohexen-1-ol offer a reasonable model for the hydroxy-directed diastereoselective SS reaction, which has been used for a long time in organic synthesis.     

Reaction pathways and mechanisms of photodegradation of pesticides

The photodegradation of pesticides is reviewed, with particular reference to the studies that describe the mechanisms of the processes involved, the nature of reactive intermediates and final products. Potential use of photochemical processes in advanced oxidation methods for water treatment is also discussed. Processes considered include direct photolysis leading to homolysis or heterolysis of the pesticide, photosensitized photodegradation by singlet oxygen and a variety of metal complexes, photolysis in heterogeneous media and degradation by reaction with intermediates generated by photolytic or radiolytic means.     

Sampling of rare events using hidden restraints

A method to enhance sampling of rare events is presented. It makes use of distance or dihedral-angle restraints to overcome an energy barrier separating two metastable states or to stabilize a transition state between the two metastable states. In order not to perturb these metastable end states themselves, a prefactor is introduced into the restraining energy function, which smoothly increases the weight of this function from zero to one at the transition state or on top of the separating energy barrier and then decreases the weight again to zero at the final state. The method is combined with multi-configurational thermodynamic integration and applied to two biomolecular systems, which were difficult to treat using standard thermodynamic integration. As first example the free energy difference of a cyclic alpha-aminoxy-hexapeptide-ion complex upon changing the ion from Cl- to Na+ was calculated. A large conformational rearrangement of the peptide was necessary to accommodate this change. Stabilizing the transition state by (hidden) restraints facilitates that. As a second example, the free energy difference between the 4C1 and the 1C4 conformation of beta-D-glucopyranoside was calculated. In unrestrained simulations the change from the 4C1 into the 1C4 conformation was never observed because of the high energy barrier separating the two states. Using (hidden) restraints, the transition from the 4C1 into the 1C4 state and back could be enforced without perturbing the end states. As comparison, for the same transitions the potential of mean force as obtained by using dihedral-angle constraints is provided.