Efficient formulation of the stochastic simulation algorithm for chemically reacting systems

In this paper we examine the different formulations of Gillespie's stochastic simulation algorithm (SSA) [D. Gillespie, J. Phys. Chem. 81, 2340 (1977)] with respect to computational efficiency, and propose an optimization to improve the efficiency of the direct method. Based on careful timing studies and an analysis of the time-consuming operations, we conclude that for most practical problems the optimized direct method is the most efficient formulation of SSA. This is in contrast to the widely held belief that Gibson and Bruck's next reaction method [M. Gibson and J. Bruck, J. Phys. Chem. A 104, 1876 (2000)] is the best way to implement the SSA for large systems. Our analysis explains the source of the discrepancy.     

A fast algorithm for simulating flow-induced nucleation in polymers

We present a fast computer simulation algorithm for high dimensional barrier crossing simulations. The algorithm is described with reference to the Graham and Olmsted (GO) model of flow-induced nucleation in polymers [R. S. Graham and P. D. Olmsted, Phys. Rev Lett. 103, 115702 (2009)]. Inspired by Chandler's barrier crossing algorithm [D. Chandler, J. Chem. Phys 68, 2959 (1978)], our algorithm simulates only the region around the top of the nucleation barrier, where the system deviates most strongly from equilibrium. When applied to the kinetic Monte Carlo (kMC) routine used in the GO model, our algorithm has two advantages: it requires very little additional coding; and it is simple enough to be applied to any barrier crossing problem that can be written in terms of a kMC simulation. Our fast nucleation algorithm is shown to vastly decrease the computer time required to perform the kMC simulations of high barrier crossing.     

Nested stochastic simulation algorithm for chemical kinetic systems with disparate rates

An efficient simulation algorithm for chemical kinetic systems with disparate rates is proposed. This new algorithm is quite general, and it amounts to a simple and seamless modification of the classical stochastic simulation algorithm (SSA), also known as the Gillespie [J. Comput. Phys. 22, 403 (1976); J. Phys. Chem. 81, 2340 (1977)] algorithm. The basic idea is to use an outer SSA to simulate the slow processes with rates computed from an inner SSA which simulates the fast reactions. Averaging theorems for Markov processes can be used to identify the fast and slow variables in the system as well as the effective dynamics over the slow time scale, even though the algorithm itself does not rely on such information. This nested SSA can be easily generalized to systems with more than two separated time scales. Convergence and efficiency of the algorithm are discussed using the established error estimates and illustrated through examples.     

Monte Carlo simulations of critical cluster sizes and nucleation rates of water

We have calculated the critical cluster sizes and homogeneous nucleation rates of water at temperatures and vapor densities corresponding to experiments by Wolk and Strey [J. Phys. Chem B 105, 11683 (2001)]. The calculations have been done with an expanded version of a Monte Carlo method originally developed by Vehkamaki and Ford [J. Chem. Phys. 112, 4193 (2000)]. Their method calculates the statistical growth and decay probabilities of molecular clusters. We have derived a connection between these probabilities and kinetic condensation and evaporation rates, and introduce a new way for the calculation of the work of formation of clusters. Three different interaction potential models of water have been used in the simulations. These include the unpolarizable SPC/E [J. Phys. Chem. 91, 6269 (1987)] and TIP4P [J. Chem. Phys. 79, 926 (1983)] models and a polarizable model by Guillot and Guissani [J. Chem. Phys. 114, 6720 (2001)]. We show that TIP4P produces critical cluster sizes and a temperature and vapor density dependence for the nucleation rate that agree well with the experimental data, although the magnitude of nucleation rate is constantly overestimated by a factor of 2 x 10(4). Guissani and Guillot's model is somewhat less successful, but both the TIP4P and Guillot and Guissani models are able to reproduce a much better experimental temperature dependency of the nucleation rate than the classical nucleation theory. Using SPC/E results in dramatically too small critical clusters and high nucleation rates. The water models give different average binding energies for clusters. We show that stronger binding between cluster molecules suppresses the decay probability of a cluster, while the growth probability is not affected. This explains the differences in results from different water models.     

Binomial distribution based tau-leap accelerated stochastic simulation

Recently, Gillespie introduced the tau-leap approximate, accelerated stochastic Monte Carlo method for well-mixed reacting systems [J. Chem. Phys. 115, 1716 (2001)]. In each time increment of that method, one executes a number of reaction events, selected randomly from a Poisson distribution, to enable simulation of long times. Here we introduce a binomial distribution tau-leap algorithm (abbreviated as BD-tau method). This method combines the bounded nature of the binomial distribution variable with the limiting reactant and constrained firing concepts to avoid negative populations encountered in the original tau-leap method of Gillespie for large time increments, and thus conserve mass. Simulations using prototype reaction networks show that the BD-tau method is more accurate than the original method for comparable coarse-graining in time.     

Structures and correlation functions of multicomponent and polydisperse hard-sphere mixtures from a density functional theory

The structures of nonuniform binary hard-sphere mixtures and the correlation functions of uniform ternary hard-sphere mixtures were studied using a modified fundamental-measure theory based on the weight functions of Rosenfeld [Rosenfeld, Phys. Rev. Lett. 63, 980 (1989)] and Boublik-Mansoori-Carnahan-Starling-Leland equation of state [Boublik, J. Chem. Phys. 53, 471 (1970); Mansoori et al., J. Chem. Phys. 54, 1523 (1971)]. The theoretical predictions agreed very well with the molecular simulations for the overall density profiles, the local compositions, and the radial distribution functions of uniform as well as inhomogeneous hard-sphere mixtures. The density functional theory was further extended to represent the structure of a polydisperse hard-sphere fluid near a hard wall. Excellent agreement was also achieved between theory and Monte Carlo simulations. The density functional theory predicted oscillatory size segregations near a hard wall for a polydisperse hard-sphere fluid of a uniform size distribution.     

Kinetic lattice grand canonical Monte Carlo simulation for ion current calculations in a model ion channel system

An algorithm in which kinetic lattice grand canonical Monte Carlo simulations are combined with mean field theory (KLGCMC/MF) is presented to calculate ion currents in a model ion channel system. In this simulation, the relevant region of the system is treated by KLGCMC simulations, while the rest of the system is described by modified Poisson-Boltzmann mean field theory. Calculation of reaction field due to induced charges on the channel/water and membrane/water boundaries is carried out using a basis-set expansion method [Im and Roux, J. Chem. Phys. 115, 4850 (2001)]. Calculation of ion currents, electrostatic potentials, and ion concentrations, as obtained from the KLGCMC/MF simulations, shows good agreement with Poisson-Nernst-Planck (PNP) theory predictions when the channel and membrane have the same dielectric constant as water. If the channel and membrane have a lower dielectric constant than water, however, there is a considerable difference between the KLGCMC/MF and PNP predictions. This difference is attributed to the reaction field, which is missing in PNP theory. It is demonstrated that the reaction field as well as fixed charges in the channel play key roles in selective ion transport. Limitations and further development of the current KLGCMC/MF approach are also discussed.     

Automated estimation of rare event probabilities in biochemical systems

In biochemical systems, the occurrence of a rare event can be accompanied by catastrophic consequences. Precise characterization of these events using Monte Carlo simulation methods is often intractable, as the number of realizations needed to witness even a single rare event can be very large. The weighted stochastic simulation algorithm (wSSA) [J. Chem. Phys. 129, 165101 (2008)] and its subsequent extension [J. Chem. Phys. 130, 174103 (2009)] alleviate this difficulty with importance sampling, which effectively biases the system toward the desired rare event. However, extensive computation coupled with substantial insight into a given system is required, as there is currently no automatic approach for choosing wSSA parameters. We present a novel modification of the wSSA--the doubly weighted SSA (dwSSA)--that makes possible a fully automated parameter selection method. Our approach uses the information-theoretic concept of cross entropy to identify parameter values yielding minimum variance rare event probability estimates. We apply the method to four examples: a pure birth process, a birth-death process, an enzymatic futile cycle, and a yeast polarization model. Our results demonstrate that the proposed method (1) enables probability estimation for a class of rare events that cannot be interrogated with the wSSA, and (2) for all examples tested, reduces the number of runs needed to achieve comparable accuracy by multiple orders of magnitude. For a particular rare event in the yeast polarization model, our method transforms a projected simulation time of 600 years to three hours. Furthermore, by incorporating information-theoretic principles, our approach provides a framework for the development of more sophisticated influencing schemes that should further improve estimation accuracy.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. I. Theory and phase-space analysis

Current constant pressure molecular-dynamics (MD) algorithms are not consistent with the recent reformulation of the isothermal-isobaric (NpT) ensemble. The NpT ensemble partition function requires the use of a "shell" molecule to identify uniquely the volume of the system, thereby avoiding the redundant counting of configurations [e.g., G. J. M. Koper and H. Reiss, J. Phys. Chem. 100, 422 (1996); D. S. Corti, Phys. Rev. E, 64, 016128 (2001)]. So far, only the NpT Monte Carlo method has been updated to allow the system volume to be defined by a shell particle [D. S. Corti, Mol. Phys. 100, 1887 (2002)]. A shell particle has yet to be incorporated into MD simulations. The proper modification of the NpT MD algorithm is therefore the subject of this paper. Unlike Andersen's method [H. C. Andersen, J. Chem. Phys. 72, 2384 (1980)] where a piston of unknown mass serves to control the response time of volume fluctuations, the newly proposed equations of motion impose a constant external pressure via the introduction of a shell particle of known mass. Hence, the system itself sets the time scales for pressure and volume fluctuations. The new algorithm is subject to a number of fundamentally rigorous tests to ensure that the equations of motion sample phase space correctly. We also show that the Hoover NpT algorithm [W. G. Hoover, Phys. Rev. A. 31, 1695 (1985); 34, 2499 (1986)] does sample phase correctly, but only when periodic boundary conditions are employed.     

Computation of the free energy of solids

We present a procedure to compute the absolute free energy of solid phases by Monte Carlo simulation. The method is based on the so-called "Einstein-crystal" method of Frenkel and Ladd [J. Chem. Phys. 81, 3188 (1984)]. The new technique is more general and simplifies the calculation for systems with hard core interactions. In addition, the reference Einstein crystal is built up to fulfill translational invariance, which seems to reduce the system size dependence of the results.     

Low-temperature nucleation in a kinetic Ising model under different stochastic dynamics with local energy barriers

Using both analytical and simulational methods, we study low-temperature nucleation rates in kinetic Ising lattice-gas models that evolve under two different Arrhenius dynamics that interpose between the Ising states a transition state representing a local energy barrier. The two dynamics are the transition-state approximation [T. Ala-Nissila, J. Kjoll, and S. C. Ying, Phys. Rev. B 46, 846 (1992)] and the one-step dynamic [H. C. Kang and W. H. Weinberg, J. Chem. Phys. 90, 2824 (1989)]. Even though they both obey detailed balance and are here applied to a situation that does not conserve the order parameter, we find significant differences between the nucleation rates observed with the two dynamics, and between them and the standard Glauber dynamic [R. J. Glauber, J. Math. Phys. 4, 294 (1963)], which does not contain transition states. Our results show that great care must be exercised when devising kinetic Monte Carlo transition rates for specific physical or chemical systems.     

The free energy of the metastable supersaturated vapor via restricted ensemble simulations

Pressure, excess chemical potential, and excess free energy, with respect to ideal gas data at different densities of the supersaturated Lennard-Jones particle vapor at the reduced temperature 0.7 are obtained by the restricted canonical ensemble Monte Carlo simulation method [D. S. Corti and P. Debenedetti, Chem. Eng. Sci. 49, 2717 (1994)]. The excess free energy values depend upon the constraints imposed on the system with local minima exhibited for densities below the spinodal density and monotonic variation for densities larger than the spinodal density. The results are compared with a molecular dynamics simulation study [A. Linharton et al., J. Chem. Phys. 122, 144506 (2005)] on the same system. The current study verifies the conclusion drawn by the simulation work that clustering of Lennard-Jones atoms exists even in the vicinity of spinodal. Our method gives an alternative to molecular dynamic simulations for the determination of equilibrium properties of a metastable fluid, especially close to the spinodal, and does not require a very large system to carry out the simulation.     

Kinetic Monte Carlo study of nucleation processes on patterned surfaces

The properties of template-directed nucleation are studied in the transition region where full nucleation control is lost and additional nucleation beyond the prepatterned structure is observed. To get deeper insight into the microscopic mechanisms, Monte Carlo simulations were performed. In this context, the previously used continuous algorithm [F. Kalischewski, J. Zhu, and A. Heuer, Phys. Rev. B 77, 155401, (2008)] was replaced by a discrete one to reduce simulation time and to allow more detailed calculations. The applied method is based on the assumption that the molecules on the surface occupy the sites of a simple fcc lattice. It is shown that a careful mapping of the continuous Monte Carlo technique onto the discrete algorithm leads to a good reproduction of the former results by means of the latter method. Furthermore, the new method facilitates the calculation of the spatial distribution of nuclei on the surface. This provides a detailed comparison with experimental data.     

Direct correlation functions of binary mixtures of hard Gaussian overlap molecules

We study the direct correlation function (DCF) of a classical fluid mixture of nonspherical molecules. The components of the mixture are two types of hard ellipsoidal molecules with different elongations, interacting through the hard Gaussian overlap (HGO) model. Two different approaches are used to calculate the DCFs of this fluid, and the results are compared. Here, the Pynn approximation [J. Chem. Phys. 60, 4579 (1974)] is extended to calculate the DCF of the binary mixtures of HGO molecules, then we use a formalism based on the weighted density functional theory introduced by Chamoux and Perera [J. Chem. Phys. 104, 1493 (1996)]. These results are fairly in agreement with each other. The pressure of this system is also calculated using the Fourier zero components of the DCF. The results are in agreement with the Monte Carlo molecular simulation.     

Using gamma distributions to predict self-diffusivities and density of states of fluids confined in carbon nanotubes

The density of states of rare gas atoms confined in carbon nanotubes is analyzed using a recently proposed model based on gamma distributions [Krishnan and Ayappa, J. Chem. Phys., 124 144503 (2006)]. The inputs into the model are the 2nd and 4th frequency moments that are obtained from molecular dynamics simulations. The predicted density of states, velocity autocorrelation functions and self-diffusivities are compared with those obtained from molecular dynamics simulations, for different nanotube loadings and temperatures. All results are reported for argon confined in a (16,16) carbon nanotube. The model predictions are extremely accurate at intermediate reduced densities of rhosigma(3) = 0.3, 0.4, where the majority of the self-diffusivity predictions lie within 10% of the simulation results. Since the frequency moments can be also obtained from Monte Carlo simulations, the study suggests an alternate route to the system dynamics of strongly confined fluids.     

Electronic conductivity in polyaromatic hydrocarbon glasses: A theoretical perspective

Based upon Monte Carlo simulations of amorphous molecular glasses, we have computed the electronic structure of five prototypical polyaromatic hydrocarbons using an extended Su-Schrieffer-Heeger model [J. R. Schrieffer, W. P. Su, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979)]. In the presence of excess charges, the resulting potential energy surfaces have been analyzed using Marcus' [Annu. Rev. Phys. Chem. 15, 155 (1964)] theory of charge transfer to yield reaction coefficients and--via the application of linear response theory--local conductivities. Applying Kirchhoff's rules, the emerging random resistor network problem leads to global conductivities of the order of 10(-1)-1 Scm, which correlate with the structural characteristics of the underlying geometry.     

Grand canonical Monte Carlo simulations of water in protein environments

The grand canonical simulation algorithm is considered as a general methodology to sample the configuration of water molecules confined within protein environments. First, the probability distribution of the number of water molecules and their configuration in a region of interest for biochemical simulations, such as the active site of a protein, is derived by considering a finite subvolume in open equilibrium with a large system serving as a bulk reservoir. It is shown that the influence of the bulk reservoir can be represented as a many-body potential of mean force acting on the atoms located inside the subvolume. The grand canonical Monte Carlo (GCMC) algorithm, augmented by a number of technical advances to increase the acceptance of insertion attempts, is implemented, and tested for simple systems. In particular, the method is illustrated in the case of a pure water box with periodic boundary conditions. In addition, finite spherical systems of pure water and containing a dialanine peptide, are simulated with GCMC while the influence of the surrounding infinite bulk is incorporated using the generalized solvent boundary potential [W. Im, S. Berneche, and B. Roux, J. Chem. Phys. 114, 2924 (2001)]. As a last illustration of water confined in the interior of a protein, the hydration of the central cavity of the KcsA potassium channel is simulated.