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.     

Path ensembles and path sampling in nonequilibrium stochastic systems

Markovian models based on the stochastic master equation are often encountered in single molecule dynamics, reaction networks, and nonequilibrium problems in chemistry, physics, and biology. An efficient and convenient method to simulate these systems is the kinetic Monte Carlo algorithm which generates continuous-time stochastic trajectories. We discuss an alternative simulation method based on sampling of stochastic paths. Utilizing known probabilities of stochastic paths, it is possible to apply Metropolis Monte Carlo in path space to generate a desired ensemble of stochastic paths. The method is a generalization of the path sampling idea to stochastic dynamics, and is especially suited for the analysis of rare paths which are not often produced in the standard kinetic Monte Carlo procedure. Two generic examples are presented to illustrate the methodology.     

Simulating rare events in equilibrium or nonequilibrium stochastic systems

We present three algorithms for calculating rate constants and sampling transition paths for rare events in simulations with stochastic dynamics. The methods do not require a priori knowledge of the phase-space density and are suitable for equilibrium or nonequilibrium systems in stationary state. All the methods use a series of interfaces in phase space, between the initial and final states, to generate transition paths as chains of connected partial paths, in a ratchetlike manner. No assumptions are made about the distribution of paths at the interfaces. The three methods differ in the way that the transition path ensemble is generated. We apply the algorithms to kinetic Monte Carlo simulations of a genetic switch and to Langevin dynamics simulations of intermittently driven polymer translocation through a pore. We find that the three methods are all of comparable efficiency, and that all the methods are much more efficient than brute-force simulation.     

Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap

We describe a replica exchange strategy where trial swap configurations are generated by nonequilibrium switching simulations. By devoting simulation time to the switching simulations, one can systematically increase an effective overlap between replicas, which leads to an increased exchange acceptance rate and less correlated equilibrium samples. In this paper, we derive our method for a general class of stochastic dynamics, and discuss various strategies for enhancing replica overlap through novel dynamical schemes and prudent choices of switching protocols. We then demonstrate our method on a model system of alanine dipeptide in implicit solvent, characterizing decreases in data correlations and gains in sampling efficiency.     

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.     

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.     

Reaction coordinates and transition pathways of rare events via forward flux sampling

A new approach is developed for identifying suitable reaction coordinates to describe the progression of rare events in complex systems. The method is based on the forward flux sampling (FFS) technique and standard least-square estimation (LSE) and it is denoted as FFS-LSE. The FFS algorithm generates trajectories for the transition between stable states as chains of partially connected paths, which can then be used to obtain "on-the-fly" estimates for the committor probability to the final region, p(B). These p(B) data are then used to screen a set of candidate collective properties for an optimal order parameter (i.e., reaction coordinate) that depends on a few relevant variables. LSE is used to find the coefficients of the proposed reaction coordinate model and an analysis of variance is used to determine the significant terms in the model. The method is demonstrated for several test systems, including the folding of a lattice protein. It is shown that a simple approximation to p(B) via a model linear on energy and number of native contacts is sufficient to describe the intrinsic dynamics of the protein system and to ensure an efficient sampling of pathways. In addition, since the p(B) surface found from the FFS-LSE approach leads to the identification of the transition state ensemble, mechanistic details of the dynamics of the system can be readily obtained during a single FFS-type simulation without the need to perform additional committor simulations.     

On equilibrium distributions and detailed balance relations in non-isothermal systems

The significance of the detailed balance principle and equilibrium solutions of the master equation is discussed from a thermodynamic point of view for isolated and isothermal systems. Starting from a master equation for all the time dependent degrees of freedom it is shown that the uniqueness of the equilibrium distribution as a stationary solution is ensured if the detailed rate constants are balanced with the aid of the distribution which maximizes the entropy subject to the thermodynamic constraints. This procedure should precede physical assumptions which simplify the original master equation, e.g. the assumption that rapidly relaxing modes can be described by canonical distribution functions.     

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.     

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.     

Ontology of temperature in nonequilibrium systems

The laws of thermodynamics provide a clear concept of the temperature for an equilibrium system in the continuum limit. Meanwhile, the equipartition theorem allows one to make a connection between the ensemble average of the kinetic energy and the uniform temperature. When a system or its environment is far from equilibrium, however, such an association does not necessarily apply. In small systems, the regression hypothesis may not even apply. Herein, it is shown that in small nonequilibrium systems, the regression hypothesis still holds, though with a generalized definition of the temperature. The latter must now be defined for each such manifestation.     

Phase behavior of polymer blends under equilibrium and nonequilibrium conditions

This contribution embraces two topics related to phase behavior of polymer blends under equilibrium and nonequilibrium. 1. Polymer blends can undergo different phase changes as liquid-liquid phase transition and crystallization. Coupling of demixing and crystallization may occur at the kinetic stage. This is illustrated by blends of poly(ϵ-caprolactone)(PCL) and poly(styrene-co-acrylonitrile)(SAN). 2. Extension of studies to blend systems under flow is necessary for the better understanding of structure formation in polymer blends outside equilibrium. Polymer molecules will be oriented and stretched when subjected to flow. This may result in flow-induced phenomena. Effects of flow on the phase behavior have been studied only for a few blends, as yet. The primary observation was flow-induced miscibility. Apparent shifts of the phase transition temperatures will be discussed qualitatively in terms of a decoupled mode theory.     

Sampling rare events in nonequilibrium and nonstationary systems

Although many computational methods for rare event sampling exist, this type of calculation is not usually practical for general nonequilibrium conditions, with macroscopically irreversible dynamics and away from both stationary and metastable states. A novel method for calculating the time-series of the probability of a rare event is presented which is designed for these conditions. The method is validated for the cases of the Glauber-Ising model under time-varying shear flow, the Kawasaki-Ising model after a quench into the region between nucleation dominated and spinodal decomposition dominated phase change dynamics, and the parallel open asymmetric exclusion process. The method requires a subdivision of the phase space of the system: it is benchmarked and found to scale well for increasingly fine subdivisions, meaning that it can be applied without detailed foreknowledge of the physically important reaction pathways.     

Multiple steady states in a CSTR with total condenser: Comparison of equilibrium and nonequilibrium models

Comparison of equilibrium and nonequilibrium models of a CSTR with total condenser focused on the multiple steady states and dynamic behaviour was carried out. The steady-state behaviour of the model system, MTBE synthesis from methanol and isobutene in a reactive distillation column, was studied in terms of the input parameters, i. e. feed flow rate of methanol or butenes, reflux ratio, and mass of catalyst. The dynamic behaviour of the system during the start-up was investigated and perturbations of manipulated variables were found to cause transitions between the parallel steady states. Presented at the 33rd International Conference of the Slovak Society of Chemical Engineering, Tatranské Matliare, 22–26 May 2006.     

Semi-classical generalized Langevin equation for equilibrium and nonequilibrium molecular dynamics simulation

Molecular dynamics (MD) simulation based on Langevin equation has been widely used in the study of structural, thermal properties of matter in different phases. Normally, the atomic dynamics are described by classical equations of motion and the effect of the environment is taken into account through the fluctuating and frictional forces. Generally, the nuclear quantum effects and their coupling to other degrees of freedom are difficult to include in an efficient way. This could be a serious limitation on its application to the study of dynamical properties of materials made from light elements, in the presence of external driving electrical or thermal fields. One example of such system is single molecule dynamics on metal surface, an important system that has received intense study in surface science. In this review, we summarize recent effort in extending the Langevin MD to include nuclear quantum effect and their coupling to flowing electrical current. We discuss its applications in the study of adsorbate dynamics on metal surface, current-induced dynamics in molecular junctions, and quantum thermal transport between different reservoirs.     

Cryofocusing in the combination of gas chromatography with equilibrium headspace sampling

Summary When headspace gas chromatography utilizing capillary columns is used for trace analysis, sample enrichment is often needed. This involves splitless sample injection of fairly large gas volumes and relatively long sampling times. As a result of this, the band of the sample vapor may be too large causing peak distortion and poor resolution. This problem can be easily overcome by the use of cryogenic trapping. While this can be accomplished by cooling the whole column to subambient temperature during sample introduction, a more convenient way is to utilize part of the first coil of the capillary column as a cryogenic trap.The paper discusses the theoretical background and instrumentation of cryogenic trapping and demonstrates the possibilities through a number of examples.Enlarged text of a paper presented at the 37th Pittsburgh Conference on Analytical Chemistry and Applied Spectroscopy, Atlantic City, NJ, March 10–14, 1986.     

Shear viscosity of molten alkali halides from equilibrium and nonequilibrium molecular-dynamics simulations

The shear viscosity of molten NaCl and KCl was calculated through equilibrium (EMD) and nonequilibrium molecular-dynamics (NEMD) simulations in the canonical (N,V,T) ensemble. Two rigid-ion potentials were investigated, namely, the Born-Mayer-Huggins-Tosi-Fumi potential and the Michielsen-Woerlee-Graaf-Ketelaar potential with the parameters proposed by Ladd. The NEMD simulations were performed using the SLLOD equations of motion [D. J. Evans and G. P. Morriss, Phys. Rev. A 30, 1528 (1984)] with a Gaussian isokinetic thermostat and the results are compared with those obtained from Green-Kubo EMD (N,V,T) simulations and experimental shear viscosity data. The NEMD zero strain rate shear viscosity, eta(0), was obtained by fitting a simplified Carreau-type equation and by application of mode-coupling theory, i.e., a eta-gamma(1/2) linear relationship. The values obtained from the first method are found to be significantly lower than those predicted by the second. The agreement between the EMD and NEMD results with experimental data is satisfactory for the two potentials investigated. The ion-ion radial distribution functions obtained with the two rigid-ion potentials for both molten salts are discussed in terms of the differences between the two models.     

Measuring flux distributions for diffusion in the small-numbers limit

For the classical diffusion of independent particles, Fick's law gives a well-known relationship between the average flux and the average concentration gradient. What has not yet been explored experimentally, however, is the dynamical distribution of diffusion rates in the limit of small particle numbers. Here, we measure the distribution of diffusional fluxes using a microfluidics device filled with a colloidal suspension of a small number of microspheres. Our experiments show that (1) the flux distribution is accurately described by a Gaussian function; (2) Fick's law, that the average flux is proportional to the particle gradient, holds even for particle gradients down to a single particle difference; (3) the variance in the flux is proportional to the sum of the particle numbers; and (4) there are backward flows, where particles flow up a concentration gradient, rather than down it. In addition, in recent years, two key theorems about nonequilibrium systems have been introduced: Evans' fluctuation theorem for the distribution of entropies and Jarzynski's work theorem. Here, we introduce a new fluctuation theorem, for the fluxes, and we find that it is confirmed quantitatively by our experiments.