Search

This paper reports the release of PathSum, a new software suite of state-of-the-art path integral methods for studying the dynamics of single or extended systems coupled to harmonic environments. The package includes two modules, suitable for system-bath problems and extended systems comprising many coupled system-bath units, and is offered in C++ and Fortran implementations. The system-bath module offers the recently developed small matrix path integral (SMatPI) and the well-established iterative quasi-adiabatic propagator path integral (i-QuAPI) method for iteration of the reduced density matrix of the system. In the SMatPI module, the dynamics within the entanglement interval can be computed using QuAPI, the blip sum, time evolving matrix product operators, or the quantum-classical path integral method. These methods have distinct convergence characteristics and their combination allows a user to access a variety of regimes. The extended system module provides the user with two algorithms of the modular path integral method, applicable to quantum spin chains or excitonic molecular aggregates. An overview of the methods and code structure is provided, along with guidance on method selection and representative examples., (© 2023 Author(s). Published under an exclusive license by AIP Publishing.)

Molecular dynamics simulations are useful in obtaining thermodynamic and kinetic properties of bio-molecules, but they are limited by the time scale barrier. That is, we may not obtain properties' efficiently because we need to run microseconds or longer simulations using femtosecond time steps. To overcome this time scale barrier, we can use the weighted ensemble (WE) method, a powerful enhanced sampling method that efficiently samples thermodynamic and kinetic properties. However, the WE method requires an appropriate partitioning of phase space into discrete macrostates, which can be problematic when we have a high-dimensional collective space or when little is known a priori about the molecular system. Hence, we developed a new WE-based method, called the "Concurrent Adaptive Sampling (CAS) algorithm," to tackle these issues. The CAS algorithm is not constrained to use only one or two collective variables, unlike most reaction coordinate-dependent methods. Instead, it can use a large number of collective variables and adaptive macrostates to enhance the sampling in the high-dimensional space. This is especially useful for systems in which we do not know what the right reaction coordinates are, in which case we can use many collective variables to sample conformations and pathways. In addition, a clustering technique based on the committor function is used to accelerate sampling the slowest process in the molecular system. In this paper, we introduce the new method and show results from two-dimensional models and bio-molecules, specifically penta-alanine and a triazine trimer.

Ligand diffusion through a protein interior is a fundamental process governing biological signaling and enzymatic catalysis. A complex topology of channels in proteins leads often to difficulties in modeling ligand escape pathways by classical molecular dynamics simulations. In this paper, two novel memetic methods for searching the exit paths and cavity space exploration are proposed: Memory Enhanced Random Acceleration (MERA) Molecular Dynamics (MD) and Immune Algorithm (IA). In MERA, a pheromone concept is introduced to optimize an expulsion force. In IA, hybrid learning protocols are exploited to predict ligand exit paths. They are tested on three protein channels with increasing complexity: M2 muscarinic G-protein-coupled receptor, enzyme nitrile hydratase, and heme-protein cytochrome P450cam. In these cases, the memetic methods outperform simulated annealing and random acceleration molecular dynamics. The proposed algorithms are general and appropriate in all problems where an accelerated transport of an object through a network of channels is studied.

Stochastic simulation for in silico studies of large biochemical networks requires a great amount of computational time. We recently proposed a new exact simulation algorithm, called the rejection-based stochastic simulation algorithm (RSSA) [Thanh et al., J. Chem. Phys. 141(13), 134116 (2014)], to improve simulation performance by postponing and collapsing as much as possible the propensity updates. In this paper, we analyze the performance of this algorithm in detail, and improve it for simulating large-scale biochemical reaction networks. We also present a new algorithm, called simultaneous RSSA (SRSSA), which generates many independent trajectories simultaneously for the analysis of the biochemical behavior. SRSSA improves simulation performance by utilizing a single data structure across simulations to select reaction firings and forming trajectories. The memory requirement for building and storing the data structure is thus independent of the number of trajectories. The updating of the data structure when needed is performed collectively in a single operation across the simulations. The trajectories generated by SRSSA are exact and independent of each other by exploiting the rejection-based mechanism. We test our new improvement on real biological systems with a wide range of reaction networks to demonstrate its applicability and efficiency.

The method calculating the effective diffusion coefficient in an arbitrary periodic two-dimensional channel, presented in our previous paper [P. Kalinay, J. Chem. Phys. 141, 144101 (2014)], is generalized to 3D channels of cylindrical symmetry, as well as to 2D or 3D channels with particles driven by a constant longitudinal external driving force. The next possible extensions are also indicated. The former calculation was based on calculus in the complex plane, suitable for the stationary diffusion in 2D domains. The method is reformulated here using standard tools of functional analysis, enabling the generalization.

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

In this paper, we propose two methods to carry out the quasi-steady state approximation in stochastic models of enzyme catalytic regulation, based on WKB asymptotics of the chemical master equation or of the corresponding partial differential equation for the generating function. The first of the methods we propose involves the development of multiscale generalisation of a WKB approximation of the solution of the master equation, where the separation of time scales is made explicit which allows us to apply the quasi-steady state approximation in a straightforward manner. To the lowest order, the multi-scale WKB method provides a quasi-steady state, Gaussian approximation of the probability distribution. The second method is based on the Hamilton-Jacobi representation of the stochastic process where, as predicted by large deviation theory, the solution of the partial differential equation for the corresponding characteristic function is given in terms of an effective action functional. The optimal transition paths between two states are then given by those paths that maximise the effective action. Such paths are the solutions of the Hamilton equations for the Hamiltonian associated to the effective action functional. The quasi-steady state approximation is applied to the Hamilton equations thus providing an approximation to the optimal transition paths and the transition time between two states. Using this approximation we predict that, unlike the mean-field quasi-steady approximation result, the rate of enzyme catalysis depends explicitly on the initial number of enzyme molecules. The accuracy and validity of our approximated results as well as that of our predictions regarding the behaviour of the stochastic enzyme catalytic models are verified by direct simulation of the stochastic model using Gillespie stochastic simulation algorithm.

The constrained density functional theory-configuration interaction (CDFT-CI) method has previously been used to calculate ground-state energies and barrier heights, and to describe electronic excited states, in particular conical intersections. However, the method has been limited to evaluating the electronic energy at just a single nuclear configuration, with the gradient of the energy being available only via finite difference. In this paper, we present analytic gradients of the CDFT-CI energy with respect to nuclear coordinates, which gives the potential for accurate geometry optimization and molecular dynamics on both the ground and excited electronic states, a realm which is currently quite challenging for electronic structure theory. We report the performance of CDFT-CI geometry optimization for representative reaction transition states as well as molecules in an excited state. The overall accuracy of CDFT-CI for computing barrier heights is essentially unchanged whether the energies are evaluated at geometries obtained from quadratic configuration-interaction singles and doubles (QCISD) or CDFT-CI, indicating that CDFT-CI produces very good reaction transition states. These results open up tantalizing possibilities for future work on excited states.

Combining simulation and model theories, this paper analyses the impact of pair dynamics on the intermolecular nuclear Overhauser effect (NOE) in liquids. For the first time, we give a distance resolved NOE. When applied to the ionic liquid 1-ethyl-3-methyl-imidazolium tetrafluoroborate the NOE turns out to be of long-range nature. This behaviour translates to the experimentally measured cross- and longitudinal relaxation rates. We were able to calculate the heteronuclear NOE from simulation data, despite the high computational effort. Model theories are computationally less demanding and cover the complete frequency range of the respective spectral density function, they are usually based on a very simple pair distribution function and the solution of the diffusion equation. In order to model the simulated data sufficiently, these simplifications in structure and dynamics have to be generalised considerably.

Network representations are becoming increasingly popular for analyzing kinetic data from techniques like Milestoning, Markov State Models, and Transition Path Theory. Mapping continuous phase space trajectories into a relatively small number of discrete states helps in visualization of the data and in dissecting complex dynamics to concrete mechanisms. However, not only are molecular networks derived from molecular dynamics simulations growing in number, they are also getting increasingly complex, owing partly to the growth in computer power that allows us to generate longer and better converged trajectories. The increased complexity of the networks makes simple interpretation and qualitative insight of the molecular systems more difficult to achieve. In this paper, we focus on various network representations of kinetic data and algorithms to identify important edges and pathways in these networks. The kinetic data can be local and partial (such as the value of rate coefficients between states) or an exact solution to kinetic equations for the entire system (such as the stationary flux between vertices). In particular, we focus on the Milestoning method that provides fluxes as the main output. We proposed Global Maximum Weight Pathways as a useful tool for analyzing molecular mechanism in Milestoning networks. A closely related definition was made in the context of Transition Path Theory. We consider three algorithms to find Global Maximum Weight Pathways: Recursive Dijkstra's, Edge-Elimination, and Edge-List Bisection. The asymptotic efficiency of the algorithms is analyzed and numerical tests on finite networks show that Edge-List Bisection and Recursive Dijkstra's algorithms are most efficient for sparse and dense networks, respectively. Pathways are illustrated for two examples: helix unfolding and membrane permeation. Finally, we illustrate that networks based on local kinetic information can lead to incorrect interpretation of molecular mechanisms.

In this paper, by applying the reduced density matrix (RDM) approach for nonorthogonal orbitals developed in the first paper of this series, efficient algorithms for matrix elements between VB structures and energy gradients in valence bond self-consistent field (VBSCF) method were presented. Both algorithms scale only as nm(4) for integral transformation and d(2)n(β)(2) for VB matrix elements and 3-RDM evaluation, while the computational costs of other procedures are negligible, where n, m, d, and n(β )are the numbers of variable occupied active orbitals, basis functions, determinants, and active β electrons, respectively. Using tensor properties of the energy gradients with respect to the orbital coefficients presented in the first paper of this series, a partial orthogonal auxiliary orbital set was introduced to reduce the computational cost of VBSCF calculation in which orbitals are flexibly defined. Test calculations on the Diels-Alder reaction of butadiene and ethylene have shown that the novel algorithm is very efficient for VBSCF calculations.

Stochastic modeling is essential for an accurate description of the biochemical network dynamics at the level of a single cell. Biochemically reacting systems often evolve on multiple time-scales, thus their stochastic mathematical models manifest stiffness. Stochastic models which, in addition, are stiff and computationally very challenging, therefore the need for developing effective and accurate numerical methods for approximating their solution. An important stochastic model of well-stirred biochemical systems is the chemical Langevin Equation. The chemical Langevin equation is a system of stochastic differential equation with multidimensional non-commutative noise. This model is valid in the regime of large molecular populations, far from the thermodynamic limit. In this paper, we propose a variable time-stepping strategy for the numerical solution of a general chemical Langevin equation, which applies for any level of randomness in the system. Our variable stepsize method allows arbitrary values of the time-step. Numerical results on several models arising in applications show significant improvement in accuracy and efficiency of the proposed adaptive scheme over the existing methods, the strategies based on halving/doubling of the stepsize and the fixed step-size ones.

In biochemical reaction systems dominated by delays, the simulation speed of the stochastic simulation algorithm depends on the size of the wait queue. As a result, it is important to control the size of the wait queue to improve the efficiency of the simulation. An improved accelerated delay stochastic simulation algorithm for biochemical reaction systems with delays, termed the improved delay-leaping algorithm, is proposed in this paper. The update method for the wait queue is effective in reducing the size of the queue as well as shortening the storage and access time, thereby accelerating the simulation speed. Numerical simulation on two examples indicates that this method not only obtains a more significant efficiency compared with the existing methods, but also can be widely applied in biochemical reaction systems with delays.

This paper revisits the issues surrounding computation of electron transport properties in water vapour as a function of E/n(0) (the ratio of the applied electric field to the water vapour number density) up to 1200 Td. We solve the Boltzmann equation using an improved version of the code of Ness and Robson [Phys. Rev. A 38, 1446 (1988)], facilitating the calculation of transport coefficients to a considerably higher degree of accuracy. This allows a correspondingly more discriminating test of the various electron-water vapour cross section sets proposed by a number of authors, which has become an important issue as such sets are now being applied to study electron driven processes in atmospheric phenomena [P. Thorn, L. Campbell, and M. Brunger, PMC Physics B 2, 1 (2009)] and in modeling charged particle tracks in matter [A. Munoz, F. Blanco, G. Garcia, P. A. Thorn, M. J. Brunger, J. P. Sullivan, and S. J. Buckman, Int. J. Mass Spectrom. 277, 175 (2008)].

In this paper, we present a theoretical approach to construct a fluctuating quantum model of the CO vibration in heme-CO proteins and its interaction with external laser fields. The methodology consists of mixed quantum-classical calculations for a restricted number of snapshots, which are then used to construct a parametrized quantum model. As an example, we calculate the infrared absorption spectrum of carboxy-hemoglobin, based on a simplified protein model, and found the absorption linewidth in good agreement with the experimental results., (© 2011 American Institute of Physics)

In this paper, we present the implementation of efficient approximations to time-dependent density functional theory (TDDFT) within the Tamm-Dancoff approximation (TDA) for hybrid density functionals. For the calculation of the TDDFT/TDA excitation energies and analytical gradients, we combine the resolution of identity (RI-J) algorithm for the computation of the Coulomb terms and the recently introduced "chain of spheres exchange" (COSX) algorithm for the calculation of the exchange terms. It is shown that for extended basis sets, the RIJCOSX approximation leads to speedups of up to 2 orders of magnitude compared to traditional methods, as demonstrated for hydrocarbon chains. The accuracy of the adiabatic transition energies, excited state structures, and vibrational frequencies is assessed on a set of 27 excited states for 25 molecules with the configuration interaction singles and hybrid TDDFT/TDA methods using various basis sets. Compared to the canonical values, the typical error in transition energies is of the order of 0.01 eV. Similar to the ground-state results, excited state equilibrium geometries differ by less than 0.3 pm in the bond distances and 0.5° in the bond angles from the canonical values. The typical error in the calculated excited state normal coordinate displacements is of the order of 0.01, and relative error in the calculated excited state vibrational frequencies is less than 1%. The errors introduced by the RIJCOSX approximation are, thus, insignificant compared to the errors related to the approximate nature of the TDDFT methods and basis set truncation. For TDDFT/TDA energy and gradient calculations on Ag-TB2-helicate (156 atoms, 2732 basis functions), it is demonstrated that the COSX algorithm parallelizes almost perfectly (speedup ~26-29 for 30 processors). The exchange-correlation terms also parallelize well (speedup ~27-29 for 30 processors). The solution of the Z-vector equations shows a speedup of ~24 on 30 processors. The parallelization efficiency for the Coulomb terms can be somewhat smaller (speedup ~15-25 for 30 processors), but their contribution to the total calculation time is small. Thus, the parallel program completes a Becke3-Lee-Yang-Parr energy and gradient calculation on the Ag-TB2-helicate in less than 4 h on 30 processors. We also present the necessary extension of the Lagrangian formalism, which enables the calculation of the TDDFT excited state properties in the frozen-core approximation. The algorithms described in this work are implemented into the ORCA electronic structure system.

Tau leaping methods enable efficient simulation of discrete stochastic chemical systems. Stiff stochastic systems are particularly challenging since implicit methods, which are good for stiffness, result in noninteger states. The occurrence of negative states is also a common problem in tau leaping. In this paper, we introduce the implicit Minkowski-Weyl tau (IMW-τ) methods. Two updating schemes of the IMW-τ methods are presented: implicit Minkowski-Weyl sequential (IMW-S) and implicit Minkowski-Weyl parallel (IMW-P). The main desirable feature of these methods is that they are designed for stiff stochastic systems with molecular copy numbers ranging from small to large and that they produce integer states without rounding. This is accomplished by the use of a split step where the first part is implicit and computes the mean update while the second part is explicit and generates a random update with the mean computed in the first part. We illustrate the IMW-S and IMW-P methods by some numerical examples, and compare them with existing tau methods. For most cases, the IMW-S and IMW-P methods perform favorably.

As originally proposed, the fewest switches surface hopping (FSSH) algorithm does not allow for decoherence between wavefunction amplitudes on different adiabatic surfaces. In this paper, we propose an inexpensive correction to standard FSSH dynamics wherein we explicitly model the decoherence of nuclear wave packets on distinct electronic surfaces. Our augmented fewest switches surface hopping approach is conceptually simple and, thus far, it has allowed us to capture several key features of the exact quantum results. Two points in particular merit attention. First, we obtain the correct branching ratios when a quantum particle passes through more than one region of nonadiabatic coupling. Second, our formalism provides a new and natural approach for rescaling nuclear momenta after a surface hop. Both of these features should become increasingly important as surface hopping schemes are applied to higher-dimensional problems.

The weighted stochastic simulation algorithm (wSSA) was developed by Kuwahara and Mura [J. Chem. Phys. 129, 165101 (2008)] to efficiently estimate the probabilities of rare events in discrete stochastic systems. The wSSA uses importance sampling to enhance the statistical accuracy in the estimation of the probability of the rare event. The original algorithm biases the reaction selection step with a fixed importance sampling parameter. In this paper, we introduce a novel method where the biasing parameter is state-dependent. The new method features improved accuracy, efficiency, and robustness.

The dynamics of many biological processes of interest, such as the folding of a protein, are slow and complicated enough that a single molecular dynamics simulation trajectory of the entire process is difficult to obtain in any reasonable amount of time. Moreover, one such simulation may not be sufficient to develop an understanding of the mechanism of the process, and multiple simulations may be necessary. One approach to circumvent this computational barrier is the use of Markov state models. These models are useful because they can be constructed using data from a large number of shorter simulations instead of a single long simulation. This paper presents a new Bayesian method for the construction of Markov models from simulation data. A Markov model is specified by (τ,P,T), where τ is the mesoscopic time step, P is a partition of configuration space into mesostates, and T is an N(P)×N(P) transition rate matrix for transitions between the mesostates in one mesoscopic time step, where N(P) is the number of mesostates in P. The method presented here is different from previous Bayesian methods in several ways. (1) The method uses Bayesian analysis to determine the partition as well as the transition probabilities. (2) The method allows the construction of a Markov model for any chosen mesoscopic time-scale τ. (3) It constructs Markov models for which the diagonal elements of T are all equal to or greater than 0.5. Such a model will be called a "consistent mesoscopic Markov model" (CMMM). Such models have important advantages for providing an understanding of the dynamics on a mesoscopic time-scale. The Bayesian method uses simulation data to find a posterior probability distribution for (P,T) for any chosen τ. This distribution can be regarded as the Bayesian probability that the kinetics observed in the atomistic simulation data on the mesoscopic time-scale τ was generated by the CMMM specified by (P,T). An optimization algorithm is used to find the most probable CMMM for the chosen mesoscopic time step. We applied this method of Markov model construction to several toy systems (random walks in one and two dimensions) as well as the dynamics of alanine dipeptide in water. The resulting Markov state models were indeed successful in capturing the dynamics of our test systems on a variety of mesoscopic time-scales.

Gillespie's direct method (DM) [D. Gillespie, J. Chem. Phys. 81, 2340 (1977)] for exact stochastic simulation of chemical reaction systems has been widely adopted. It is easy to implement but requires large computation for relatively large systems. Recently, two more efficient methods, next reaction method (NRM) [M. A. Gibson and J. Bruck, J. Phys. Chem. A 105, 1876 (2000)] and optimized DM (ODM) [Y. Cao et al., J. Chem. Phys. 121, 4059 (2004)], have been developed to improve simulation speed. It has been demonstrated that the ODM is the state-of-the-art most efficient method for exact stochastic simulation of most practical reaction systems. In this paper, we first develop an exact stochastic simulation algorithm named ODMK that is more efficient than the ODM. We then develop an approximate method named K-skip method to further accelerate simulation. Using two chemical reaction systems, we demonstrate that our ODMK and K-skip method can save 20%-30% and 70%-80% simulation time, respectively, comparing to the ODM. We also show that our ODMK and K-skip method provide almost the same simulation accuracy as the ODM.

A new milestoning procedure using Voronoi tessellations is proposed. In the new procedure, the edges of Voronoi cells are used as milestones, and the necessary kinetic information about the transitions between the milestones is calculated by running molecular dynamics (MD) simulations restricted to these cells. Like the traditional milestoning technique, the new procedure offers a reduced description of the original dynamics and permits to efficiently compute the various quantities necessary in this description. However, unlike traditional milestoning, the new procedure does not require to reinitialize trajectories from the milestones, and thereby it avoids the approximation made in traditional milestoning that the distribution for reinitialization is the equilibrium one. In this paper we concentrate on Markovian milestoning, which we show to be valid under suitable assumptions, and we explain how to estimate the rate matrix of transitions between the milestones from data collected from the MD trajectories in the Voronoi cells. The rate matrix can then be used to compute mean first passage times between milestones and reaction rates. The procedure is first illustrated on test-case examples in two dimensions and then applied to study the kinetics of protein insertion into a lipid bilayer by means of a coarse-grained model.

The weighted stochastic simulation algorithm (wSSA) recently introduced by Kuwahara and Mura [J. Chem. Phys. 129, 165101 (2008)] is an innovative variation on the stochastic simulation algorithm (SSA). It enables one to estimate, with much less computational effort than was previously thought possible using a Monte Carlo simulation procedure, the probability that a specified event will occur in a chemically reacting system within a specified time when that probability is very small. This paper presents some procedural extensions to the wSSA that enhance its effectiveness in practical applications. The paper also attempts to clarify some theoretical issues connected with the wSSA, including its connection to first passage time theory and its relation to the SSA.

In this paper we present a moment closure method for stochastically modeled chemical or biochemical reaction networks. We derive a system of differential equations which describes the dynamics of means and all central moments from a chemical master equation. Truncating the system for the central moments at a certain moment term and using Taylor approximation, we obtain explicit representations of means and covariances and even higher central moments in recursive forms. This enables us to deal with the moments in successive differential equations and use conventional numerical methods for their approximations. Furthermore, we estimate the errors in the means and central moments generated by the approximation method. We also find the moments at equilibrium by solving truncated algebraic equations. We show in examples that numerical solutions based on the moment closure method are accurate and efficient by comparing the results to those of stochastic simulation algorithms.

Vanishing of the equilibrium fluctuation expression for the friction coefficient of a massive particle in a finite-volume liquid has been well documented and discussed in literature. This paper investigates the decay of the friction force in the corresponding nonequilibrium situation, when the massive particle moves through a finite volume at a constant velocity. The friction force ultimately vanishes (with the decay form as predicted by the equilibrium integral) because of the finite mass of the rest of the system, which allows it to be dragged by the moving particle. However, it is sufficient to have two infinite masses moving relative to each other in a finite liquid volume for the friction force to be finite at all times.

The delay stochastic simulation algorithm (DSSA) by Barrio et al. [Plos Comput. Biol. 2, 117(E) (2006)] was developed to simulate delayed processes in cell biology in the presence of intrinsic noise, that is, when there are small-to-moderate numbers of certain key molecules present in a chemical reaction system. These delayed processes can faithfully represent complex interactions and mechanisms that imply a number of spatiotemporal processes often not explicitly modeled such as transcription and translation, basic in the modeling of cell signaling pathways. However, for systems with widely varying reaction rate constants or large numbers of molecules, the simulation time steps of both the stochastic simulation algorithm (SSA) and the DSSA can become very small causing considerable computational overheads. In order to overcome the limit of small step sizes, various tau-leap strategies have been suggested for improving computational performance of the SSA. In this paper, we present a binomial tau-DSSA method that extends the tau-leap idea to the delay setting and avoids drawing insufficient numbers of reactions, a common shortcoming of existing binomial tau-leap methods that becomes evident when dealing with complex chemical interactions. The resulting inaccuracies are most evident in the delayed case, even when considering reaction products as potential reactants within the same time step in which they are produced. Moreover, we extend the framework to account for multicellular systems with different degrees of intercellular communication. We apply these ideas to two important genetic regulatory models, namely, the hes1 gene, implicated as a molecular clock, and a Her1/Her 7 model for coupled oscillating cells.

The so-called "ultrafast" nuclear magnetic resonance (NMR) methods enable the collection of multidimensional spectra within a single scan. These experiments operate by replacing traditional t(1) time increments, with a series of combined radiofrequency-irradiation/magnetic-field-gradient manipulations that spatially encode the effects of the indirect-domain spin interactions. Barring the presence of sizable displacements, the spatial patterns thus imparted can be read out following a mixing period with the aid of oscillating acquisition gradients, leading to a train of t(2)-modulated echoes carrying in their positions and phases the indirect- and the direct-domain spin interactions. Both the initial spatial encoding as well as the subsequent spatial decoding procedures underlying ultrafast NMR were designed under the assumption that spins remain static within the sample during their execution. Most often this is not the case, and motion-related effects can be expected to affect the outcome of these experiments. The present paper focuses on analyzing the effects of diffusion in ultrafast two-dimensional (2D) NMR. Toward this end both analytical and numerical formalisms are derived, capable of dealing with the nonuniform spin manipulations, macroscopic sample sizes, and microscopic displacements involved in this kind of sequences. After experimentally validating the correctness of these formalisms these were used to analyze the effects of diffusion for a variety of cases, including ultrafast experiments on both rapidly and slowly diffusing molecules. A series of prototypical schemes were considered including discrete and continuous encoding modes, constant- and real-time manipulations, homo- and heteronuclear acquisitions, and single versus multiple quantum modalities. The effects of molecular diffusion were also compared against typical relaxation-driven losses as they happen in these various prototypical situations; from all these situations, general guidelines for choosing the optimal ultrafast 2D NMR scheme for a particular sample and condition could be deduced.

The tau-leap method first developed by Gillespie [D. T. Gillespie, J. Chem. Phys. 115, 1716 (2001)] can significantly speed up stochastic simulation of certain chemically reacting systems with acceptable losses in accuracy. Recently, several improved tau-leap methods, including the binomial, multinomial, and modified tau-leap methods, have been developed. However, in all these tau-leap methods, the mean of the number of times, K(m), that the mth reaction channel fires during a leap is not equal to the true mean. Therefore, all existing tau-leap methods produce biased simulation results, which limit the simulation accuracy and speed. In this paper, we analyze the mean of K(m) based on the chemical master equation. Using this analytical result, we develop unbiased Poisson and binomial tau-leap methods. Moreover, we analyze the variance of K(m), and then develop an unbiased Poisson/Gaussian/binomial tau-leap method to correct the errors in both the mean and variance of K(m). Simulation results demonstrate that our unbiased tau-leap method can significantly improve simulation accuracy without sacrificing speed.

In this paper a Clifford algebra-based method is applied to calculate polymer chain conformations. The approach enables the calculation of the position of an atom in space with the knowledge of the bond length (l), valence angle (theta), and rotation angle (phi) of each of the preceding bonds in the chain. Hence, the set of geometrical parameters {l(i),theta(i),phi(i)} yields all the position coordinates p(i) of the main chain atoms. Moreover, the method allows the calculation of side chain conformations and the computation of rotations of chain segments. With these features it is, in principle, possible to generate conformations of any type of chemical structure. This method is proposed as an alternative for the classical approach by matrix algebra. It is more straightforward and its final symbolic representation considerably simpler than that of matrix algebra. Approaches for realistic modeling by means of incorporation of energetic considerations can be combined with it. This article, however, is entirely focused at showing the suitable mathematical framework on which further developments and applications can be built.

The molecular dynamics expression of heat flux, originally derived by Irving and Kirkwood [J. Chem. Phys. 18, 817 (1950)] for pairwise potentials, is generalized in this paper for systems with many-body potentials. The original formula consists of a kinetic part and a potential part, and the latter term is found in the present study to be expressible as a summation of contributions from all the many-body potentials defined in the system. The energy transfer among a set of sites for which a many-body potential is defined is discussed and evaluated by the rate of increase in the kinetic energy of each site due to the potential, and its accumulation over all the potentials in the system is shown to make up the potential part of the generalized expression. A molecular dynamics simulation for liquid n-octane was performed to demonstrate the applicability of the new expression obtained in this study to measure the heat flux and to elucidate the contributions of inter- and intramolecular potentials to heat conduction.

We present an extension of the molecular mechanics-valence bond (MMVB) hybrid method to study ground and excited states of planar conjugated hydrocarbon cations. Currently, accurate excited state calculations on these systems are limited to expensive ab initio studies of smaller systems: up to 15 active electrons in 16 pi orbitals with complete active space self-consistent field (CASSCF) theory using high symmetry. The new MMVB extension provides a faster, cheaper treatment to investigate larger cation systems with more than 24 active orbitals. Extension requires both new matrix elements and new parameters: In this paper we present both, for the limited planar case. The scheme is tested for the planar radical cations of benzene, naphthalene, anthracene, and phenanthrene. Calculated MMVB relative energies are in good agreement with CASSCF results for equilibrium geometries on the ground and first excited states, and conical intersections.

Differing perspectives on the accuracy of three-electron reduced-density-matrix (3-RDM) reconstruction in nonminimal basis sets exist in the literature. This paper demonstrates the accuracy of cumulant-based reconstructions, developed by Valdemoro (V) [F. Colmenero et al., Phys. Rev. A 47, 971 (1993)], Nakatsuji and Yasuda (NY) [Phys. Rev. Lett. 76, 1039 (1996)], Mazziotti (M) [Phys. Rev. A 60, 3618 (1999)], and Valdemoro-Tel-Perez-Romero (VTP) [Many-electron Densities and Density Matrices, edited by J. Cioslowski (Kluwer, Boston, 2000)]. Computationally, we extend previous investigations to study a variety of molecules, including LiH, HF, NH(3), H(2)O, and N(2), in Slater-type, double-zeta, and polarized double-zeta basis sets at both equilibrium and nonequilibrium geometries. The reconstructed 3-RDMs, compared with 3-RDMs from full configuration interaction, demonstrate in nonminimal basis sets the accuracy of the first-order expansion (V) as well as the important role of the second-order corrections (NY, M, and VTP). Calculations at nonequilibrium geometries further show that cumulant functionals can reconstruct the 3-RDM from a multireferenced 2-RDM with reasonable accuracy, which is relevant to recent multireferenced formulations of the anti-Hermitian contracted Schrodinger equation (ACSE) and canonical diagonalization. Theoretically, we perform a detailed perturbative analysis of the M functional to identify its second-order components. With these second-order components we connect the M, NY, and VTP reconstructions for the first time by deriving both the NY and VTP functionals from the M functional. Finally, these 3-RDM reconstructions are employed within the ACSE [D. Mazziotti, Phys. Rev. Lett. 97, 143002 (2006)] to compute ground-state energies which are compared with the energies from the contracted Schrodinger equation and several wave function methods.

This paper extends our previous work on obtaining reaction coordinates from aimless shooting and likelihood maximization. We introduce a simplified version of aimless shooting and a half-trajectory likelihood score based on the committor probability. Additionally, we analyze and compare the absolute log-likelihood score for perfect and approximate reaction coordinates. We also compare the aimless shooting and likelihood maximization approach to the earlier genetic neural network (GNN) approach of Ma and Dinner [J. Phys. Chem. B 109, 6769 (2005)]. For a fixed number of total trajectories in the GNN approach, the accuracy of the transition state ensemble decreases as the number of trajectories per committor probability estimate increases. This quantitatively demonstrates the benefit of individual committor probability realizations over committor probability estimates. Furthermore, when the least squares score of the GNN approach is applied to individual committor probability realizations, the likelihood score still provides a better approximation to the true transition state surface. Finally, the polymorph transition in terephthalic acid demonstrates that the new half-trajectory likelihood scheme estimates the transition state location more accurately than likelihood schemes based on the probability of being on a transition path.

The thermodynamics and kinetics of a many-body system can be described in terms of a potential energy landscape in multidimensional configuration space. The partition function of such a landscape can be written in terms of a density of states, which can be computed using a variety of Monte Carlo techniques. In this paper, a new self-consistent Monte Carlo method for computing density of states is described that uses importance sampling and a multiplicative update factor to achieve rapid convergence. The technique is then applied to compute the equilibrium quench probability of the various inherent structures (minima) in the landscape. The quench probability depends on both the potential energy of the inherent structure and the volume of its corresponding basin in configuration space. Finally, the methodology is extended to the isothermal-isobaric ensemble in order to compute inherent structure quench probabilities in an enthalpy landscape.

This paper reports the first computational estimation of the comb polymers' third virial coefficients. The number of the chains in the comb polymers range from 5 to 11. An algorithm that counts the contributing terms of the third virial coefficients in an accelerated manner is presented along with its efficiency dependence on the polymers' size. In addition, the second virial coefficients are estimated for the comb polymers and compared to previously reported results.

Finding coarse-grained, low-dimensional descriptions is an important task in the analysis of complex, stochastic models of gene regulatory networks. This task involves (a) identifying observables that best describe the state of these complex systems and (b) characterizing the dynamics of the observables. In a previous paper [R. Erban et al., J. Chem. Phys. 124, 084106 (2006)] the authors assumed that good observables were known a priori, and presented an equation-free approach to approximate coarse-grained quantities (i.e., effective drift and diffusion coefficients) that characterize the long-time behavior of the observables. Here we use diffusion maps [R. Coifman et al., Proc. Natl. Acad. Sci. U.S.A. 102, 7426 (2005)] to extract appropriate observables ("reduction coordinates") in an automated fashion; these involve the leading eigenvectors of a weighted Laplacian on a graph constructed from network simulation data. We present lifting and restriction procedures for translating between physical variables and these data-based observables. These procedures allow us to perform equation-free, coarse-grained computations characterizing the long-term dynamics through the design and processing of short bursts of stochastic simulation initialized at appropriate values of the data-based observables.

Molecular ground-state energies and two-electron reduced density matrices (2-RDMs) have recently been computed without the many-electron wave function by constraining the 2-RDM to satisfy a complete set of three-positivity conditions for N representability [D. A. Mazziotti, Phys. Rev. A 74, 032501 (2006)]. Energies at both equilibrium and nonequilibrium geometries are obtained within 0.3% of the correlation energy. In this paper the authors extend this work to examine the accuracy of molecular properties, including multipole moments and components of the ground-state energy, relative to full configuration interaction (FCI). Comparisons are also made with 2-RDM methods with two-positivity conditions and two-positivity plus the generalized T1T2 conditions as well as several approximate wave function methods. Using the 2-RDM method with three-positivity conditions, the authors obtain dipole, quadrupole, and octupole moments for BeH2, BH, H2O, CO, and NH3 at equilibrium geometries that are within 0.04% of their FCI values. In addition, for the potential energy surface of N2, the 2-RDM method with three-positivity yields not only accurate total ground-state energies but also accurate expectation values of the kinetic energy operator, the electron-nuclei potential, and electron-electron repulsion.

Spatially distributed problems are often approximately modeled in terms of partial differential equations (PDEs) for appropriate coarse-grained quantities (e.g., concentrations). The derivation of accurate such PDEs starting from finer scale, atomistic models, and using suitable averaging is often a challenging task; approximate PDEs are typically obtained through mathematical closure procedures (e.g., mean field approximations). In this paper, we show how such approximate macroscopic PDEs can be exploited in constructing preconditioners to accelerate stochastic computations for spatially distributed particle-based process models. We illustrate how such preconditioning can improve the convergence of equation-free coarse-grained methods based on coarse timesteppers. Our model problem is a stochastic reaction-diffusion model capable of exhibiting Turing instabilities.

The equilibrium free energy difference between two long-lived molecular species or "conformational states" of a protein (or any other molecule) can in principle be estimated by measuring the work needed to shuttle the system between them, independent of the irreversibility of the process. This is the meaning of the Jarzynski equality (JE), which we test in this paper by performing simulations that unfold a protein by pulling two atoms apart. Pulling is performed fast relative to the relaxation time of the molecule and is thus far from equilibrium. Choosing a simple protein model for which we can independently compute its equilibrium properties, we show that the free energy can be exactly and effectively estimated from nonequilibrium simulations. To do so, one must carefully and correctly determine the ensemble of states that are pulled, which is more important the farther from equilibrium one performs simulations; this highlights a potential problem in using the JE to extract the free energy from forced unfolding experiments. The results presented here also demonstrate that the free energy difference between the native and denatured states of a protein measured in solution is not always equal to the free energy profile that can be estimated from forced unfolding simulations (or experiments) using the JE.

Determining the functional conformation of a protein from its amino acid sequence remains a central problem in computational biology. In this paper, we establish the mathematical optimal model of protein folding problem (PFP) on two-dimensional space based on the minimal energy principle. A novel hybrid of elastic net algorithm and local search method (ENLS) is applied successfully to simulations of protein folding on two-dimensional hydrophobic-polar (HP) lattice model. Eight HP benchmark instances with up to 64 amino acids are tested to verify the effectiveness of proposed approach and model. In several cases, the ENLS method finds new lower energy states. The numerical results show that it is drastically superior to other methods in finding the ground state of a protein.

Special structures in a chemical reaction network can give rise to bistability, oscillations, and chaos. It has been shown recently [A. Sensse and M. Eiswirth, J. Chem. Phys. 122, 044516 (2005)] that the introduction of an additional species in a supplementary feedback loop to a minimal autocatalytic oscillator gives rise to chaotic dynamics in a certain range of parameters, independent of the particular realization of the additional loop. This provides a possibility to decide if chaos may occur just by analyzing the network structure of an existing model. Here, we apply this concept to analyze the complex dynamics in several essential subsystems of the peroxidase-oxidase reaction system. The aim of the present paper is to determine the nature of the occurring chaos and its location in the parameter space by numerical bifurcation analysis and simulations.

The evaluation of interactions between nearby particles constitutes the majority of the computational workload involved in classical molecular dynamics (MD) simulations. In this paper, we introduce a new method for the parallelization of range-limited particle interactions that proves particularly suitable to MD applications. Because it applies not only to pairwise interactions but also to interactions involving three or more particles, the method can be used for evaluation of both nonbonded and bonded forces in a MD simulation. It requires less interprocessor data transfer than traditional spatial decomposition methods at all but the lowest levels of parallelism. It gains an additional practical advantage in certain commonly used interprocessor communication networks by distributing the communication burden more evenly across network links and by decreasing the associated latency. When used to parallelize MD, it further reduces communication requirements by allowing the computations associated with short-range nonbonded interactions, long-range electrostatics, bonded interactions, and particle migration to use much of the same communicated data. We also introduce certain variants of this method that can significantly improve the balance of computational load across processors.

When a minimum on the potential energy surface is surrounded by multiple saddle points with similar energy barriers, the transition pathways with greater prefactors are more important than those that have similar energy barriers but smaller prefactors. In this paper, we present a theoretical formulation for the prefactors, computing the probabilities for transition paths from a minimum to its surrounding saddle points. We apply this formulation to a system of 2 degrees of freedom and a system of 14 degrees of freedom. The first is Brownian motion in a two-dimensional potential whose global anharmonicities play a dominant role in determining the transition rates. The second is a Lennard-Jones (LJ) cluster of seven particles in two dimensions. Low lying transition states of the LJ cluster, which can be reached directly from a minimum without passing through another minimum, are identified without any presumption of their characteristics nor of the product states they lead to. The probabilities are computed for paths going from an equilibrium ensemble of states near a given minimum to the surrounding transition states. These probabilities are directly related to the prefactors in the rate formula. This determination of the rate prefactors includes all anharmonicities, near or far from transition states, which are pertinent in the very sophisticated energy landscape of LJ clusters and in many other complex systems.

In this paper, we present simulations of the decay of quantum coherence between vibrational states of I(2) in its ground (X) electronic state embedded in a cryogenic Kr matrix. We employ a numerical method based on the semiclassical limit of the quantum Liouville equation, which allows the simulation of the evolution and decay of quantum vibrational coherence using classical trajectories and ensemble averaging. The vibrational level-dependent interaction of the I(2)(X) oscillator with the rare-gas environment is modeled using a recently developed method for constructing state-dependent many-body potentials for quantum vibrations in a many-body classical environment [J. M. Riga, E. Fredj, and C. C. Martens, J. Chem. Phys. 122, 174107 (2005)]. The vibrational dephasing rates gamma(0n) for coherences prepared between the ground vibrational state mid R:0 and excited vibrational state mid R:n are calculated as a function of n and lattice temperature T. Excellent agreement with recent experiments performed by Karavitis et al. [Phys. Chem. Chem. Phys. 7, 791 (2005)] is obtained.

Triplet-triplet (TT) energy transfer requires two molecular fragments to exchange electrons that carry different spin and energy. In this paper, we analyze and report values of the electronic coupling strengths for TT energy transfer. Two different methods were proposed and tested: (1) Directly calculating the off-diagonal Hamiltonian matrix element. This direct coupling scheme was generalized from the one used for electron transfer coupling, where two spin-localized unrestricted Hartree-Fock wave functions are used as the zero-order reactant and product states, and the off-diagonal Hamiltonian matrix elements are calculated directly. (2) From energy gaps derived from configuration-interaction-singles (CIS) scheme. Both methods yielded very similar results for the systems tested. For TT coupling between a pair of face-to-face ethylene molecules, the exponential attenuation factor is 2.59 A(-1)(CIS6-311+G(**)), which is about twice as large as typical values for electron transfer. With a series of fully stacked polyene pairs, we found that the TT coupling magnitudes and attenuation rates are very similar irrespective of their molecular size. If the polyenes were partially stacked, TT couplings were much reduced, and they decay more rapidly with distance than those of full-stacked systems. Our results showed that the TT coupling arises mainly from the region of close contact between the donor and acceptor frontier orbitals, and the exponential decay of the coupling with separation depends on the details of the molecular contacts. With our calculated results, nanosecond or picosecond time scales for TT energy-transfer rates are possible.

In this paper we present the results of a stochastic model of reversible biochemical reaction networks that are being driven through an open boundary, such that the system is interacting with its surrounding environment with explicit material exchange. The stochastic model is based on the master equation approach and is intimately related to the grand canonical ensemble of statistical mechanics. We show that it is possible to analytically calculate the joint probability function of the random variables describing the number of molecules in each state of the system for general linear networks. Definitions of reaction chemical potentials and conductances follow from inherent properties of this model, providing a description of energy dissipation in the system. We are also able to suggest novel methods for experimentally determining reaction fluxes and biochemical affinities at nonequilibrium steady state as well as the overall network connectivity.

The finite-temperature string method proposed by E, et al. [W. E, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 (2002)] is a very effective way of identifying transition mechanisms and transition rates between metastable states in systems with complex energy landscapes. In this paper, we discuss the theoretical background and algorithmic details of the finite-temperature string method, as well as the application to the study of isomerization reaction of the alanine dipeptide, both in vacuum and in explicit solvent. We demonstrate that the method allows us to identify directly the isocommittor surfaces, which are approximated by hyperplanes, in the region of configuration space where the most probable transition trajectories are concentrated. These results are verified subsequently by computing directly the committor distribution on the hyperplanes that define the transition state region.

We describe the dynamic behavior of a 1-stearoyl-2-oleoyl-phosphatidylethanolamine (SOPE) bilayer from a 20 ns molecular dynamics simulation. The dynamics of individual molecules are characterized in terms of (2)H spin-lattice relaxation rates, nuclear overhauser enhancement spectroscopy (NOESY) cross-relaxation rates, and lateral diffusion coefficients. Additionally, we describe the dynamics of hydrogen bonding through an analysis of hydrogen bond lifetimes and the time evolution of clusters of hydrogen bonded lipids. The simulated trajectory is shown to be consistent with experimental measures of internal, intermolecular, and diffusive motion. Consistent with our analysis of SOPE structure in the companion paper, we see hydrogen bonding dominating the dynamics of the interface region. Comparison of (2)H T(1) relaxation rates for chain methylene segments in phosphatidylcholine and phosphatidylethanolamine bilayers indicates that slower motion resulting from hydrogen bonding extends at least three carbons into the hydrophobic core. NOESY cross-relaxation rates compare well with experimental values, indicating the observed hydrogen bonding dynamics are realistic. Calculated lateral diffusion rates (4 +/ -1 x 10(-8) cm(2)s) are comparable, though somewhat lower than, those determined by pulsed field gradient NMR methods.

In this paper, we present a new analytical approach for describing the spin dynamics of synchronous and asynchronous time-dependent modulations in solid-state nuclear magnetic resonance experiments. The approach, based on multimode Floquet theory, employs the multipole operator basis of Sanctuary for spin description and illustrates the time evolution in the Floquet-Liouville space using the effective Hamiltonians obtained from the contact (or van Vleck) transformation procedure. Since the Hamiltonian and the density operator are expressed in terms of irreducible tensor operators, extensions to higher spin magnitudes (I>12) and multiple spins are quite straightforward and permit analytical treatments for many problems. We outline the general underlying principles involved in this approach with a brief mention of its potential application in other branches of spectroscopy.

The dynamical solution of a well-mixed, nonlinear stochastic chemical kinetic system, described by the Master equation, may be exactly computed using the stochastic simulation algorithm. However, because the computational cost scales with the number of reaction occurrences, systems with one or more "fast" reactions become costly to simulate. This paper describes a hybrid stochastic method that partitions the system into subsets of fast and slow reactions, approximates the fast reactions as a continuous Markov process, using a chemical Langevin equation, and accurately describes the slow dynamics using the integral form of the "Next Reaction" variant of the stochastic simulation algorithm. The key innovation of this method is its mechanism of efficiently monitoring the occurrences of slow, discrete events while simultaneously simulating the dynamics of a continuous, stochastic or deterministic process. In addition, by introducing an approximation in which multiple slow reactions may occur within a time step of the numerical integration of the chemical Langevin equation, the hybrid stochastic method performs much faster with only a marginal decrease in accuracy. Multiple examples, including a biological pulse generator and a large-scale system benchmark, are simulated using the exact and proposed hybrid methods as well as, for comparison, a previous hybrid stochastic method. Probability distributions of the solutions are compared and the weak errors of the first two moments are computed. In general, these hybrid methods may be applied to the simulation of the dynamics of a system described by stochastic differential, ordinary differential, and Master equations.