How accurate are the nonlinear chemical Fokker-Planck and chemical Langevin equations?

The chemical Fokker-Planck equation and the corresponding chemical Langevin equation are commonly used approximations of the chemical master equation. These equations are derived from an uncontrolled, second-order truncation of the Kramers-Moyal expansion of the chemical master equation and hence their accuracy remains to be clarified. We use the system-size expansion to show that chemical Fokker-Planck estimates of the mean concentrations and of the variance of the concentration fluctuations about the mean are accurate to order Ω(-3∕2) for reaction systems which do not obey detailed balance and at least accurate to order Ω(-2) for systems obeying detailed balance, where Ω is the characteristic size of the system. Hence, the chemical Fokker-Planck equation turns out to be more accurate than the linear-noise approximation of the chemical master equation (the linear Fokker-Planck equation) which leads to mean concentration estimates accurate to order Ω(-1∕2) and variance estimates accurate to order Ω(-3∕2). This higher accuracy is particularly conspicuous for chemical systems realized in small volumes such as biochemical reactions inside cells. A formula is also obtained for the approximate size of the relative errors in the concentration and variance predictions of the chemical Fokker-Planck equation, where the relative error is defined as the difference between the predictions of the chemical Fokker-Planck equation and the master equation divided by the prediction of the master equation. For dimerization and enzyme-catalyzed reactions, the errors are typically less than few percent even when the steady-state is characterized by merely few tens of molecules.     

Stochastic simulation of catalytic surface reactions in the fast diffusion limit

The master equation of a lattice gas reaction tracks the probability of visiting all spatial configurations. The large number of unique spatial configurations on a lattice renders master equation simulations infeasible for even small lattices. In this work, a reduced master equation is derived for the probability distribution of the coverages in the infinite diffusion limit. This derivation justifies the widely used assumption that the adlayer is in equilibrium for the current coverages and temperature when all reactants are highly mobile. Given the reduced master equation, two novel and efficient simulation methods of lattice gas reactions in the infinite diffusion limit are derived. The first method involves solving the reduced master equation directly for small lattices, which is intractable in configuration space. The second method involves reducing the master equation further in the large lattice limit to a set of differential equations that tracks only the species coverages. Solution of the reduced master equation and differential equations requires information that can be obtained through short, diffusion-only kinetic Monte Carlo simulation runs at each coverage. These simulations need to be run only once because the data can be stored and used for simulations with any set of kinetic parameters, gas-phase concentrations, and initial conditions. An idealized CO oxidation reaction mechanism with strong lateral interactions is used as an example system for demonstrating the reduced master equation and deterministic simulation techniques.     

Efficient stochastic simulations of complex reaction networks on surfaces

Surfaces serve as highly efficient catalysts for a vast variety of chemical reactions. Typically, such surface reactions involve billions of molecules which diffuse and react over macroscopic areas. Therefore, stochastic fluctuations are negligible and the reaction rates can be evaluated using rate equations, which are based on the mean-field approximation. However, in case that the surface is partitioned into a large number of disconnected microscopic domains, the number of reactants in each domain becomes small and it strongly fluctuates. This is, in fact, the situation in the interstellar medium, where some crucial reactions take place on the surfaces of microscopic dust grains. In this case rate equations fail and the simulation of surface reactions requires stochastic methods such as the master equation. However, in the case of complex reaction networks, the master equation becomes infeasible because the number of equations proliferates exponentially. To solve this problem, we introduce a stochastic method based on moment equations. In this method the number of equations is dramatically reduced to just one equation for each reactive species and one equation for each reaction. Moreover, the equations can be easily constructed using a diagrammatic approach. We demonstrate the method for a set of astrophysically relevant networks of increasing complexity. It is expected to be applicable in many other contexts in which problems that exhibit analogous structure appear, such as surface catalysis in nanoscale systems, aerosol chemistry in stratospheric clouds, and genetic networks in cells.     

Markovian approximation in the relaxation of open quantum systems

In this paper, we examine the validity of the Markovian approximation and the slippage scheme used to incorporate short time transient memory effects in the Markovian master equations (Redfield equations). We argue that for a bath described by a spectral function, J(omega), that is dense and smoothly spread out over the range omega(d), a time scale of tau(b) approximately 1/omega(d) exists; for times of t > tau(b), the Markovian approximation is applicable. In addition, if J(omega) decays to zero reasonably fast in both the omega --> 0 and omega --> infinity limits, then the bath relaxation time, tau(b), is determined by the width of the spectral function and is weakly dependent on the temperature of the bath. On the basis of this criterion of tau(b), a scheme to incorporate transient memory effects in the Markovian master equation is suggested. Instead of using slipped initial conditions, we propose a concatenation scheme that uses the second-order perturbation theory for short time dynamics and the Markovian master equation at long times. Application of this concatenation scheme to the spin-boson model shows that it reproduces the reduced dynamics obtained from the non-Markovian master equation for all parameters studied, while the simple slippage scheme breaks down at high temperatures.     

Hopping transport and energy transfer in substitutionally disordered media

We present an alternative derivation of configuration-averaged equations governing transport and trapping in a three-component disordered system. It is assumed that migration of the initial excitation is described by the usual master equation. Using the Zwanzig projection operator technique, we were able, after averaging, to separate the contribution of the transport and transfer processes. We found finally that the averaged master equation including sinks is equivalent to a hierarchy of coupled equations.     

Transition path sampling for discrete master equations with absorbing states

Transition path sampling (TPS) algorithms have been implemented with deterministic dynamics, with thermostatted dynamics, with Brownian dynamics, and with simple spin flip dynamics. Missing from the TPS repertoire is an implementation with kinetic Monte Carlo (kMC), i.e., with the underlying dynamics coming from a discrete master equation. We present a new hybrid kMC-TPS algorithm and prove that it satisfies detailed balance in the transition path ensemble. The new algorithm is illustrated for a simplified Markov State Model of trp-cage folding. The transition path ensemble from kMC-TPS is consistent with that obtained from brute force kMC simulations. The committor probabilities and local fluxes for the simple model are consistent with those obtained from exact methods for simple master equations. The new kMC-TPS method should be useful for analysis of rare transitions in complex master equations where the individual states cannot be enumerated and therefore where exact solutions cannot be obtained.     

Numerical solution methods for large,difficult kinetic master equations

The kinetics of gas-phase reactions, including pressure-dependent weak collision and non-equilibrium effects, can be modelled using a master equation. In this paper, we address the practical computational problem of finding solutions to such kinetic master equations. The mathematical structure of the master equation can be utilised to develop a number of specialised numerical techniques that are capable of solving the master equation in the presence of difficult numerics and for large problems. The former is important for modelling low temperature and pressure systems, and the latter is important for modelling the large networks of isomerising species common in combustion chemistry applications. We focus on numerical methods that exhibit particular practical use because of their robust nature or scalability to many isomers, or both. Recent developments in linear-scaling methods are highlighted.     

Trajectory approach to dissipative quantum phase space dynamics: Application to barrier scattering

The Caldeira-Leggett master equation, expressed in Lindblad form, has been used in the numerical study of the effect of a thermal environment on the dynamics of the scattering of a wave packet from a repulsive Eckart barrier. The dynamics are studied in terms of phase space trajectories associated with the distribution function, W(q,p,t). The equations of motion for the trajectories include quantum terms that introduce nonlocality into the motion, which imply that an ensemble of correlated trajectories needs to be propagated. However, use of the derivative propagation method (DPM) allows each trajectory to be propagated individually. This is achieved by deriving equations of motion for the partial derivatives of W(q,p,t) that appear in the master equation. The effects of dissipation on the trajectories are studied and results are shown for the transmission probability. On short time scales, decoherence is demonstrated by a swelling of trajectories into momentum space. For a nondissipative system, a comparison is made of the DPM with the "exact" transmission probability calculated from a fixed grid calculation.     

Master equation approach to vibrational energy transfer between two diatomic molecules

In this paper a semiclassical non-markovian master equation is derived. We begin by using the well-known tetradic form of the Liouville equation for a reduced density operator. By projecting the diagonal matrix elements of the operator, we obtain an infinite-order master equation. This equation is then applied in the lowest-order approximation to collinear collisions between the diatomic molecules: H₂H₂, N₂N₂ and Cl₂Cl₂. With an assumed form of the interaction potential for such a problem we have also derived an analytical expression for the V—V transition probabilities. They are then calculated over a wide range of velocities of the colliding molecules and compared with exact semiclassical ones. An excellent agreement of the results is found for small velocities (i.e. υ ≈ 10⁴ cm/s). For larger values of υ (≈ 10⁵ cm/s) the results obtained from the master equation approach agree with the exact ones only in the low-velocity range for light molecules and low oscillatory states.     

A comparative study of the master equation approach and the discrete galerkin method for the simulation of free radical polymerization

The present article deals with the mathematical treatment of free radical polymerization reactions. As a typical example the synthesis of poly(methyl methacrylate) under realistic experimental conditions is investigated. Since the mathematical treatment of the kinetic rate equations raises severe numerical problems, alternative approaches are required. In this paper two of these methods, i.e. the discrete Galerkin method and the master equation approach, are compared. The discrete Galerkin method circumvents difficulties encountered by the direct integration of the kinetic rate equations but requires much a priori knowledge of the chemical reaction system. Within the framework of the master equation approach the polymerization reaction is regarded as a stochastic process. For the simulation of this stochastic process a modified algorithm is presented. The example of the polymerization of methyl methacrylate shows that the master equation approach is an efficient tool in the simulation of free radical polymerization reactions.     

A numerical solution of the pair equation of a model two‐electron diatomic system

It has been well‐documented that about 90% of the total correlation energy of atomic systems can be obtained by solving so‐called pair equations. For atoms, this approach requires solving partial differential equations (PDE) in two variables. In case of a diatomic molecule, we face devising a method for treating PDEs in five variables. This article shows how a well‐established finite difference method used to solve Hartree–Fock equations for diatomic molecules can be extended to solve numerically a model two‐electron Schrödinger equation for such systems. We show that using less than 100 grid points in each variable, it is possible to obtain the total energy of the helium atom and hydrogen molecule with a chemical accuracy and the S energy of the helium atom and hydride ion as accurately as the best results available. © 2015 Wiley Periodicals, Inc.     

Stochastic theory of direct simulation Monte Carlo method

A treatment of direct simulation Monte Carlo method as a Markov process with a master equation is given and the corresponding master equation is derived. A hierarchy of equations for the reduced probability distributions is derived from the master equation. An equation similar to the Boltzmann equation for single particle probability distribution is derived using assumption of molecular chaos. It is shown that starting from an uncorrelated state, the system remains uncorrelated always in the limit N→∞, where N is the number of particles. Simple applications of the formalism to direct simulation money games are given as examples to the formalism. The formalism is applied to the direct simulation of homogenous gases. It is shown that appropriately normalized single particle probability distribution satisfies the Boltzmann equation for simple gases and Wang Chang–Uhlenbeck equation for a mixture of molecular gases. As a consequence of this development we derive Birds no time counter algorithm. We extend the analysis to the inhomogeneous gases and define a new direct simulation algorithm for this case. We show that single particle probability distribution satisfies the Boltzmann equation in our algorithm in the limit N→∞, V _k →0, Δt→0 where V _k is the volume of kth cell. We also show that our algorithm and Bird’s algorithm approach each other in the limit N _k →∞ where N _k is the number of particles in the volume V _k .     

Coarse master equation from Bayesian analysis of replica molecular dynamics simulations

We use Bayesian inference to derive the rate coefficients of a coarse master equation from molecular dynamics simulations. Results from multiple short simulation trajectories are used to estimate propagators. A likelihood function constructed as a product of the propagators provides a posterior distribution of the free coefficients in the rate matrix determining the Markovian master equation. Extensions to non-Markovian dynamics are discussed, using the trajectory "paths" as observations. The Markovian approach is illustrated for the filling and emptying transitions of short carbon nanotubes dissolved in water. We show that accurate thermodynamic and kinetic properties, such as free energy surfaces and kinetic rate coefficients, can be computed from coarse master equations obtained through Bayesian inference.     

A kinetic approach to the theory of heterogeneous nucleation on soluble particles during the deliquescence stage

Deliquescence is the dissolution of a solid nucleus in a liquid film formed on the nucleus due to vapor condensation. Previously, the kinetics of deliquescence was examined in the framework of the capillarity approximation which involves the thermodynamic interfacial tensions for a thin film and the approximation of uniform density therein. In the present paper we propose a kinetic approach to the theory of deliquescence which avoids the use of the above macroscopic quantities for thin films. The rates of emission of molecules from the liquid film into the vapor and from the solid core into the liquid film are determined through a first passage time analysis whereas the respective rates of absorption are calculated through the gas kinetic theory. The first passage time is obtained by solving the single-molecule master equation for the probability distribution of a "surface" molecule moving in a potential field created by the cluster. Furthermore, the time evolution of the liquid film around the solid core is described by means of two mass balance equations which involve the rates of absorption and emission of molecules by the film at its two interfaces. When the deliquescence of an ensemble of solid particles occurs by means of large fluctuations, the time evolution of the distribution of composite droplets (liquid film+solid core) with respect to the independent variables of state is governed by a Fokker-Planck kinetic equation. When both the vapor and the solid soluble particles are single component, this equation has the form of the kinetic equation of binary nucleation. A steady-state solution for this equation is obtained by the method of separation of variables. The theory is illustrated with numerical calculation regarding the deliquescence of spherical particles in a water vapor with intermolecular interactions of the Lennard-Jones kind. The new approach allows one to qualitatively explain an important feature of experimental data on deliquescence, namely the occurrence of nonsharp deliquescence, a feature that the previous deliquescence theory based on classical thermodynamics could not account for.     

Relative Boltzmann entropy, evolution equations for fluctuations of thermodynamic intensive variables, and a statistical mechanical representation of the zeroth law of thermodynamics

Generalized thermodynamics or extended irreversible thermodynamics presumes the existence of thermodynamic intensive variables (e.g., temperature, pressure, chemical potentials, generalized potentials) even if the system is removed from equilibrium. It is necessary to properly understand the nature of such intensive variables and, in particular, of their fluctuations, that is, their deviations from those defined in the extended irreversible thermodynamic sense. The meaning of temperature is examined by means of a kinetic theory of macroscopic irreversible processes to assess the validity of the generalized (or extended) thermodynamic method applied to nonequilibrium phenomena. The Boltzmann equation is used for the purpose. Since the relative Boltzmann entropy has been known to be intimately related to the evolution of the aforementioned fluctuations in the intensive thermodynamic variables, we derive the evolution equations for such fluctuations of intensive variables to lay the foundation for investigating the physical implications and evolution of the relative Boltzmann entropy, so that the range of validity of the thermodynamic theory of irreversible processes can be elucidated. Within the framework of this work, we examine a special case of the evolution equations for the aforementioned fluctuations of intensive variables, which also facilitate investigation of the molecular theory meaning of the zeroth law of thermodynamics. We derive an evolution equation describing the relaxation of temperature fluctuations from its local value and present a formula for the temperature relaxation time.     

On the origins of approximations for stochastic chemical kinetics

This paper considers the derivation of approximations for stochastic chemical kinetics governed by the discrete master equation. Here, the concepts of (1) partitioning on the basis of fast and slow reactions as opposed to fast and slow species and (2) conditional probability densities are used to derive approximate, partitioned master equations, which are Markovian in nature, from the original master equation. Under different conditions dictated by relaxation time arguments, such approximations give rise to both the equilibrium and hybrid (deterministic or Langevin equations coupled with discrete stochastic simulation) approximations previously reported. In addition, the derivation points out several weaknesses in previous justifications of both the hybrid and equilibrium systems and demonstrates the connection between the original and approximate master equations. Two simple examples illustrate situations in which these two approximate methods are applicable and demonstrate the two methods' efficiencies.     

A Note on the Initial Condition of the Differential Equation Which Defines Proper Quantum Topological Subspaces

We analyze Baders variational procedure (BVP) to define Proper Quantum Topological Subspaces (PQTS) in multiatomic systems (MAS). In particular, we explicitly treat the problem of the initial condition of the partial differential equation resulting from BVP and defining the boundary of PQTS. We show that in general for MAS where nuclei correspond to spherical regions of constant (uniform) distribution the proper initial condition associated with the equation defining single atoms automatically emerges from the fundamental physical and mathematical hypothesis on which BVP is based. Our results justify the uniqueness of Baders partition of atoms in molecules on the basis of an a priori mathematical argument implicitly contained in the theory rather than on an a posteriori chemical one as done so far.     

Wigner function approach to the quantum Brownian motion of a particle in a potential

Recent progress in our understanding of quantum effects on the Brownian motion in an external potential is reviewed. This problem is ubiquitous in physics and chemistry, particularly in the context of decay of metastable states, for example, the reversal of the magnetization of a single domain ferromagnetic particle, kinetics of a superconducting tunnelling junction, etc. Emphasis is laid on the establishment of master equations describing the diffusion process in phase space analogous to the classical Fokker-Planck equation. In particular, it is shown how Wigner's [E. P. Wigner, Phys. Rev., 1932, 40, 749] method of obtaining quantum corrections to the classical equilibrium Maxwell-Boltzmann distribution may be extended to the dissipative non-equilibrium dynamics governing the quantum Brownian motion in an external potential V(x), yielding a master equation for the Wigner distribution function W(x,p,t) in phase space (x,p). The explicit form of the master equation so obtained contains quantum correction terms up to o(h(4)) and in the classical limit, h --> 0, reduces to the classical Klein-Kramers equation. For a quantum oscillator, the method yields an evolution equation coinciding in all respects with that of Agarwal [G. S. Agarwal, Phys. Rev. A, 1971, 4, 739]. In the high dissipation limit, the master equation reduces to a semi-classical Smoluchowski equation describing non-inertial quantum diffusion in configuration space. The Wigner function formulation of quantum Brownian motion is further illustrated by finding quantum corrections to the Kramers escape rate, which, in appropriate limits, reduce to those yielded via quantum generalizations of reaction rate theory.     

Stochastic theory of large-scale enzyme-reaction networks: finite copy number corrections to rate equation models

Chemical reactions inside cells occur in compartment volumes in the range of atto- to femtoliters. Physiological concentrations realized in such small volumes imply low copy numbers of interacting molecules with the consequence of considerable fluctuations in the concentrations. In contrast, rate equation models are based on the implicit assumption of infinitely large numbers of interacting molecules, or equivalently, that reactions occur in infinite volumes at constant macroscopic concentrations. In this article we compute the finite-volume corrections (or equivalently the finite copy number corrections) to the solutions of the rate equations for chemical reaction networks composed of arbitrarily large numbers of enzyme-catalyzed reactions which are confined inside a small subcellular compartment. This is achieved by applying a mesoscopic version of the quasisteady-state assumption to the exact Fokker-Planck equation associated with the Poisson representation of the chemical master equation. The procedure yields impressively simple and compact expressions for the finite-volume corrections. We prove that the predictions of the rate equations will always underestimate the actual steady-state substrate concentrations for an enzyme-reaction network confined in a small volume. In particular we show that the finite-volume corrections increase with decreasing subcellular volume, decreasing Michaelis-Menten constants, and increasing enzyme saturation. The magnitude of the corrections depends sensitively on the topology of the network. The predictions of the theory are shown to be in excellent agreement with stochastic simulations for two types of networks typically associated with protein methylation and metabolism.