Transition path sampling

Last updated

Transition path sampling (TPS) is a rare-event sampling method used in computer simulations of rare events: physical or chemical transitions of a system from one stable state to another that occur too rarely to be observed on a computer timescale. Examples include protein folding, chemical reactions and nucleation. Standard simulation tools such as molecular dynamics can generate the dynamical trajectories of all the atoms in the system. However, because of the gap in accessible time-scales between simulation and reality, even present supercomputers might require years of simulations to show an event that occurs once per millisecond without some kind of acceleration.

Contents

Transition path ensemble

TPS focuses on the most interesting part of the simulation, the transition. For example, an initially unfolded protein will vibrate for a long time in an open-string configuration before undergoing a transition and fold on itself. The aim of the method is to reproduce precisely those folding moments.

Consider in general a system with two stable states A and B. The system will spend a long time in those states and occasionally jump from one to the other. There are many ways in which the transition can take place. Once a probability is assigned to each of the many pathways, one can construct a Monte Carlo random walk in the path space of the transition trajectories, and thus generate the ensemble of all transition paths. All the relevant information can then be extracted from the ensemble, such as the reaction mechanism, the transition states, and the rate constants.

Given an initial path, TPS provides some algorithms to perturb that path and create a new one. As in all Monte Carlo walks, the new path will then be accepted or rejected in order to have the correct path probability. The procedure is iterated and the ensemble is gradually sampled.

A powerful and efficient algorithm is the so-called shooting move. [1] Consider the case of a classical many-body system described by coordinates r and momenta p. Molecular dynamics generates a path as a set of (rt, pt) at discrete times t in [0,T] where T is the length of the path. For a transition from A to B, (r0, p0) is in A, and (rT, pT) is in B. One of the path times is chosen at random, the momenta p are modified slightly into p + δp, where δp is a random perturbation consistent with system constraints, e.g. conservation of energy and linear and angular momentum. A new trajectory is then simulated from this point, both backward and forward in time until one of the states is reached. Being in a transition region, this will not take long. If the new path still connects A to B it is accepted, otherwise it is rejected and the procedure starts again.

Rate constant computation

In the Bennett–Chandler procedure, [2] [3] the rate constant kAB for the transition from A to B is derived from the correlation function

,

where hX is the characteristic function of state X, and hX(t) is either 1 if the system at time t is in state X or 0 if not. The time-derivative C'(t) starts at time 0 at the transition state theory (TST) value kABTST and reaches a plateau kABkABTST for times of the order of the transition time. Hence once the function is known up to these times, the rate constant is also available.

In the TPS framework C(t) can be rewritten as an average in the path ensemble

,

where the subscript AB denotes an average in the ensemble of paths that start in A and visit B at least once. Time t' is an arbitrary time in the plateau region of C(t). The factor C(t') at this specific time can be computed with a combination of path sampling and umbrella sampling.

Transition interface sampling

The TPS rate constant calculation can be improved in a variation of the method called Transition interface sampling (TIS). [4] In this method the transition region is divided in subregions using interfaces. The first interface defines state A and the last state B. The interfaces are not physical interfaces but hypersurfaces in the phase space.

The rate constant can be viewed as a flux through these interfaces. The rate kAB is the flux of trajectories starting before the first interface and going through the last interface. Being a rare event, the flux is very small and practically impossible to compute with a direct simulation. However, using the other interfaces between the states, one can rewrite the flux in terms of transition probabilities between interfaces

,

where PA(i + 1|i) is the probability for trajectories, coming from state A and crossing interface i, to reach interface i + 1. Here interface 0 defines state A and interface n defines state B. The factor Φ1,0 is the flux through the interface closest to A. By making this interface close enough, the quantity can be computed with a standard simulation, as the crossing event through this interface is not a rare event any more.

Remarkably, in the formula above there is no Markov assumption of independent transition probabilities. The quantities PA(i + 1|i) carry a subscript A to indicate that the probabilities are all dependent on the history of the path, all the way from when it left A. These probabilities can be computed with a path sampling simulation using the TPS shooting move. A path crossing interface i is perturbed and a new path is shot. If the path still starts from A and crosses interface i, is accepted. The probability PA(i + 1|i) follows from the ratio of the number of paths that reach interface i + 1 to the total number of paths in the ensemble.

Theoretical considerations show that TIS computations are at least twice as fast as TPS, and computer experiments have shown that the TIS rate constant can converge up to 10 times faster. A reason for this is due to TIS using paths of adjustable length and on average shorter than TPS. Also, TPS relies on the correlation function C(t), computed by summation of positive and negative terms due to recrossings. TIS instead computes the rate as an effective positive flux, the quantity kAB is directly computed as an average of only positive terms contributing to the interface transition probabilities.

Time Dependent Processes

TPS/TIS as normally implemented can be acceptable for non-equilibrium calculations provided that the interfacial fluxes are time-independent (stationary). To treat non-stationary systems in which there is time dependence in the dynamics, due either to variation of an external parameter or to evolution of the system itself, then other rare-event methods may be needed, such as stochastic-process rare-event sampling. [5]

Cited references

  1. Dellago, Christoph; Bolhuis, Peter G.; Chandler, David (1998). "Efficient transition path sampling: Application to Lennard-Jones cluster rearrangements". The Journal of Chemical Physics. 108 (22): 9236. Bibcode:1998JChPh.108.9236D. doi:10.1063/1.476378.
  2. Chandler, David (1978). "Statistical mechanics of isomerization dynamics in liquids and the transition state approximation". The Journal of Chemical Physics. 68 (6): 2959–2970. Bibcode:1978JChPh..68.2959C. doi:10.1063/1.436049.
  3. Bennett, C. H. (1977). Christofferson, R. (ed.). Algorithms for Chemical Computations, ACS Symposium Series No. 46. Washington, D.C.: American Chemical Society. ISBN   978-0-8412-0371-6.
  4. Van Erp, Titus S.; Moroni, Daniele; Bolhuis, Peter G. (2003). "A novel path sampling method for the calculation of rate constants". The Journal of Chemical Physics. 118 (17): 7762. arXiv: cond-mat/0210614 . Bibcode:2003JChPh.118.7762V. doi:10.1063/1.1562614. S2CID   94328349.
  5. Berryman, Joshua T.; Schilling, Tanja (2010). "Sampling rare events in nonequilibrium and nonstationary systems". The Journal of Chemical Physics. 133 (24): 244101. arXiv: 1001.2456 . Bibcode:2010JChPh.133x4101B. doi:10.1063/1.3525099. PMID   21197970. S2CID   34154184.

More references

For a review of TPS:

For a review of TIS

Related Research Articles

<span class="mw-page-title-main">BQP</span> Computational complexity class of problems

In computational complexity theory, bounded-error quantum polynomial time (BQP) is the class of decision problems solvable by a quantum computer in polynomial time, with an error probability of at most 1/3 for all instances. It is the quantum analogue to the complexity class BPP.

The thermal conductivity of a material is a measure of its ability to conduct heat. It is commonly denoted by , , or and is measured in W·m−1·K−1.

In physics, a partition function describes the statistical properties of a system in thermodynamic equilibrium. Partition functions are functions of the thermodynamic state variables, such as the temperature and volume. Most of the aggregate thermodynamic variables of the system, such as the total energy, free energy, entropy, and pressure, can be expressed in terms of the partition function or its derivatives. The partition function is dimensionless.

<span class="mw-page-title-main">Rabi cycle</span> Quantum mechanical phenomenon

In physics, the Rabi cycle is the cyclic behaviour of a two-level quantum system in the presence of an oscillatory driving field. A great variety of physical processes belonging to the areas of quantum computing, condensed matter, atomic and molecular physics, and nuclear and particle physics can be conveniently studied in terms of two-level quantum mechanical systems, and exhibit Rabi flopping when coupled to an optical driving field. The effect is important in quantum optics, magnetic resonance and quantum computing, and is named after Isidor Isaac Rabi.

The classical XY model is a lattice model of statistical mechanics. In general, the XY model can be seen as a specialization of Stanley's n-vector model for n = 2.

The Green–Kubo relations give the exact mathematical expression for transport coefficients in terms of integrals of time correlation functions:

<span class="mw-page-title-main">Wigner quasiprobability distribution</span> Wigner distribution function in physics as opposed to in signal processing

The Wigner quasiprobability distribution is a quasiprobability distribution. It was introduced by Eugene Wigner in 1932 to study quantum corrections to classical statistical mechanics. The goal was to link the wavefunction that appears in Schrödinger's equation to a probability distribution in phase space.

The Eyring equation is an equation used in chemical kinetics to describe changes in the rate of a chemical reaction against temperature. It was developed almost simultaneously in 1935 by Henry Eyring, Meredith Gwynne Evans and Michael Polanyi. The equation follows from the transition state theory, also known as activated-complex theory. If one assumes a constant enthalpy of activation and constant entropy of activation, the Eyring equation is similar to the empirical Arrhenius equation, despite the Arrhenius equation being empirical and the Eyring equation based on statistical mechanical justification.

In probability theory, the Gillespie algorithm generates a statistically correct trajectory of a stochastic equation system for which the reaction rates are known. It was created by Joseph L. Doob and others, presented by Dan Gillespie in 1976, and popularized in 1977 in a paper where he uses it to simulate chemical or biochemical systems of reactions efficiently and accurately using limited computational power. As computers have become faster, the algorithm has been used to simulate increasingly complex systems. The algorithm is particularly useful for simulating reactions within cells, where the number of reagents is low and keeping track of every single reaction is computationally feasible. Mathematically, it is a variant of a dynamic Monte Carlo method and similar to the kinetic Monte Carlo methods. It is used heavily in computational systems biology.

Photon polarization is the quantum mechanical description of the classical polarized sinusoidal plane electromagnetic wave. An individual photon can be described as having right or left circular polarization, or a superposition of the two. Equivalently, a photon can be described as having horizontal or vertical linear polarization, or a superposition of the two.

A stochastic simulation is a simulation of a system that has variables that can change stochastically (randomly) with individual probabilities.

Resonance fluorescence is the process in which a two-level atom system interacts with the quantum electromagnetic field if the field is driven at a frequency near to the natural frequency of the atom.

Thermodynamic integration is a method used to compare the difference in free energy between two given states whose potential energies and have different dependences on the spatial coordinates. Because the free energy of a system is not simply a function of the phase space coordinates of the system, but is instead a function of the Boltzmann-weighted integral over phase space, the free energy difference between two states cannot be calculated directly from the potential energy of just two coordinate sets. In thermodynamic integration, the free energy difference is calculated by defining a thermodynamic path between the states and integrating over ensemble-averaged enthalpy changes along the path. Such paths can either be real chemical processes or alchemical processes. An example alchemical process is the Kirkwood's coupling parameter method.

The Bennett acceptance ratio method (BAR) is an algorithm for estimating the difference in free energy between two systems . It was suggested by Charles H. Bennett in 1976.

Rare event sampling is an umbrella term for a group of computer simulation methods intended to selectively sample 'special' regions of the dynamic space of systems which are unlikely to visit those special regions through brute-force simulation. A familiar example of a rare event in this context would be nucleation of a raindrop from over-saturated water vapour: although raindrops form every day, relative to the length and time scales defined by the motion of water molecules in the vapour phase, the formation of a liquid droplet is extremely rare.

Stochastic-process rare event sampling (SPRES) is a rare-event sampling method in computer simulation, designed specifically for non-equilibrium calculations, including those for which the rare-event rates are time-dependent. To treat systems in which there is time dependence in the dynamics, due either to variation of an external parameter or to evolution of the system itself, the scheme for branching paths must be devised so as to achieve sampling which is distributed evenly in time and which takes account of changing fluxes through different regions of the phase space.

The min-entropy, in information theory, is the smallest of the Rényi family of entropies, corresponding to the most conservative way of measuring the unpredictability of a set of outcomes, as the negative logarithm of the probability of the most likely outcome. The various Rényi entropies are all equal for a uniform distribution, but measure the unpredictability of a nonuniform distribution in different ways. The min-entropy is never greater than the ordinary or Shannon entropy and that in turn is never greater than the Hartley or max-entropy, defined as the logarithm of the number of outcomes with nonzero probability.

Surface hopping is a mixed quantum-classical technique that incorporates quantum mechanical effects into molecular dynamics simulations. Traditional molecular dynamics assume the Born-Oppenheimer approximation, where the lighter electrons adjust instantaneously to the motion of the nuclei. Though the Born-Oppenheimer approximation is applicable to a wide range of problems, there are several applications, such as photoexcited dynamics, electron transfer, and surface chemistry where this approximation falls apart. Surface hopping partially incorporates the non-adiabatic effects by including excited adiabatic surfaces in the calculations, and allowing for 'hops' between these surfaces, subject to certain criteria.

<span class="mw-page-title-main">Newton-X</span> Molecular dynamics simulation software

Newton-X is a general program for molecular dynamics simulations beyond the Born-Oppenheimer approximation. It has been primarily used for simulations of ultrafast processes in photoexcited molecules. It has also been used for simulation of band envelops of absorption and emission spectra.

Feynman's algorithm is an algorithm that is used to simulate the operations of a quantum computer on a classical computer. It is based on the Path integral formulation of quantum mechanics, which was formulated by Richard Feynman.