The kinetic Monte Carlo (KMC) method is a Monte Carlo method computer simulation intended to simulate the time evolution of some processes occurring in nature. Typically these are processes that occur with known transition rates among states. These rates are inputs to the KMC algorithm; the method itself cannot predict them.
The KMC method is essentially the same as the dynamic Monte Carlo method and the Gillespie algorithm.
One possible classification of KMC algorithms is as rejection-KMC (rKMC) and rejection-free-KMC (rfKMC).
A rfKMC algorithm, often only called KMC, for simulating the time evolution of a system, where some processes can occur with known rates r, can be written for instance as follows: [1] : 13–14
(Note: because the average value of is equal to unity, the same average time scale can be obtained by instead using in step 9. In this case, however, the delay associated with transition i will not be drawn from the Poisson distribution described by the rate , but will instead be the mean of that distribution.[ citation needed ])
This algorithm is known in different sources variously as the residence-time algorithm or the n-fold way or the Bortz-Kalos-Lebowitz (BKL) algorithm. It is important to note that the timestep involved is a function of the probability that all events i, did not occur. [1] : 13–14
Rejection KMC has typically the advantage of an easier data handling, and faster computations for each attempted step, since the time consuming action of getting all is not needed. On the other hand, the time evolved at each step is smaller than for rfKMC. The relative weight of pros and cons varies with the case at hand, and with available resources.
An rKMC associated with the same transition rates as above can be written as follows:
(Note: can change from one MC step to another.) This algorithm is usually called a standard algorithm.
Theoretical [2] and numerical [1] [3] comparisons between the algorithms were provided.
If the rates are time dependent, step 9 in the rfKMC must be modified by: [4]
The reaction (step 6) has to be chosen after this by
Another very similar algorithm is called the First Reaction Method (FRM). It consists of choosing the first-occurring reaction, meaning to choose the smallest time , and the corresponding reaction number i, from the formula
where the are N random numbers.
The key property of the KMC algorithm (and of the FRM one) is that if the rates are correct, if the processes associated with the rates are of the Poisson process type, and if different processes are independent (i.e. not correlated) then the KMC algorithm gives the correct time scale for the evolution of the simulated system. There was some debate about the correctness of the time scale for rKMC algorithms, but this was also rigorously shown to be correct. [2]
If furthermore the transitions follow detailed balance, the KMC algorithm can be used to simulate thermodynamic equilibrium. However, KMC is widely used to simulate non-equilibrium processes, [5] in which case detailed balance need not be obeyed.
The rfKMC algorithm is efficient in the sense that every iteration is guaranteed to produce a transition. However, in the form presented above it requires operations for each transition, which is not too efficient. In many cases this can be much improved on by binning the same kinds of transitions into bins, and/or forming a tree data structure of the events. A constant-time scaling algorithm of this type has recently been developed and tested. [6]
The major disadvantage with rfKMC is that all possible rates and reactions have to be known in advance. The method itself can do nothing about predicting them. The rates and reactions must be obtained from other methods, such as diffusion (or other) experiments, molecular dynamics or density-functional theory simulations.
KMC has been used in simulations of the following physical systems:
To give an idea what the "objects" and "events" may be in practice, here is one concrete simple example, corresponding to example 2 above.
Consider a system where individual atoms are deposited on a surface one at a time (typical of physical vapor deposition), but also may migrate on the surface with some known jump rate . In this case the "objects" of the KMC algorithm are simply the individual atoms.
If two atoms come right next to each other, they become immobile. Then the flux of incoming atoms determines a rate rdeposit, and the system can be simulated with KMC considering all deposited mobile atoms which have not (yet) met a counterpart and become immobile. This way there are the following events possible at each KMC step:
After an event has been selected and carried out with the KMC algorithm, one then needs to check whether the new or just jumped atom has become immediately adjacent to some other atom. If this has happened, the atom(s) which are now adjacent needs to be moved away from the list of mobile atoms, and correspondingly their jump events removed from the list of possible events.
Naturally in applying KMC to problems in physics and chemistry, one has to first consider whether the real system follows the assumptions underlying KMC well enough. Real processes do not necessarily have well-defined rates, the transition processes may be correlated, in case of atom or particle jumps the jumps may not occur in random directions, and so on. When simulating widely disparate time scales one also needs to consider whether new processes may be present at longer time scales. If any of these issues are valid, the time scale and system evolution predicted by KMC may be skewed or even completely wrong.
The first publication which described the basic features of the KMC method (namely using a cumulative function to select an event and a time scale calculation of the form 1/R) was by Young and Elcock in 1966. [8] The residence-time algorithm was also published at about the same time. [10]
Apparently independent of the work of Young and Elcock, Bortz, Kalos and Lebowitz [1] developed a KMC algorithm for simulating the Ising model, which they called the n-fold way. The basics of their algorithm is the same as that of Young, [8] but they do provide much greater detail on the method.
The following year Dan Gillespie published what is now known as the Gillespie algorithm to describe chemical reactions. [11] The algorithm is similar and the time advancement scheme essentially the same as in KMC.
There is as of the writing of this (June 2006) no definitive treatise of the theory of KMC, but Fichthorn and Weinberg have discussed the theory for thermodynamic equilibrium KMC simulations in detail. [12] A good introduction is given also by Art Voter, [13] and by A.P.J. Jansen, [14] , and a recent review is (Chatterjee 2007) [15] or (Chotia 2008). [16] The justification of KMC as a coarse-graining of the Langevin dynamics using the quasi-stationary distribution approach has been developed by T. Lelièvre and collaborators. [17] [18]
In March 2006 the, probably, first commercial software using Kinetic Monte Carlo to simulate the diffusion and activation/deactivation of dopants in Silicon and Silicon like materials is released by Synopsys, reported by Martin-Bragado et al. [19]
The KMC method can be subdivided by how the objects are moving or reactions occurring. At least the following subdivisions are used:
In statistics and statistical physics, the Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. New samples are added to the sequence in two steps: first a new sample is proposed based on the previous sample, then the proposed sample is either added to the sequence or rejected depending on the value of the probability distribution at that point. The resulting sequence can be used to approximate the distribution or to compute an integral.
In finance, the binomial options pricing model (BOPM) provides a generalizable numerical method for the valuation of options. Essentially, the model uses a "discrete-time" model of the varying price over time of the underlying financial instrument, addressing cases where the closed-form Black–Scholes formula is wanting.
In plasma physics, the particle-in-cell (PIC) method refers to a technique used to solve a certain class of partial differential equations. In this method, individual particles in a Lagrangian frame are tracked in continuous phase space, whereas moments of the distribution such as densities and currents are computed simultaneously on Eulerian (stationary) mesh points.
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.
In physical organic chemistry, a kinetic isotope effect (KIE) is the change in the reaction rate of a chemical reaction when one of the atoms in the reactants is replaced by one of its isotopes. Formally, it is the ratio of rate constants for the reactions involving the light (kL) and the heavy (kH) isotopically substituted reactants (isotopologues):
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.
A stochastic simulation is a simulation of a system that has variables that can change stochastically (randomly) with individual probabilities.
Umbrella sampling is a technique in computational physics and chemistry, used to improve sampling of a system where ergodicity is hindered by the form of the system's energy landscape. It was first suggested by Torrie and Valleau in 1977. It is a particular physical application of the more general importance sampling in statistics.
In thermodynamics, the excess chemical potential is defined as the difference between the chemical potential of a given species and that of an ideal gas under the same conditions . The chemical potential of a particle species is therefore given by an ideal part and an excess part.
The Wang and Landau algorithm, proposed by Fugao Wang and David P. Landau, is a Monte Carlo method designed to estimate the density of states of a system. The method performs a non-Markovian random walk to build the density of states by quickly visiting all the available energy spectrum. The Wang and Landau algorithm is an important method to obtain the density of states required to perform a multicanonical simulation.
The BFM is a lattice model for simulating the conformation and dynamics of polymer systems. There are two versions of the BFM used: The earlier version was first introduced by I. Carmesin and Kurt Kremer in 1988, and the later version by J. Scott Shaffer in 1994. Conversion between models is possible.
The Swendsen–Wang algorithm is the first non-local or cluster algorithm for Monte Carlo simulation for large systems near criticality. It has been introduced by Robert Swendsen and Jian-Sheng Wang in 1987 at Carnegie Mellon.
The Nosé–Hoover thermostat is a deterministic algorithm for constant-temperature molecular dynamics simulations. It was originally developed by Nosé and was improved further by Hoover. Although the heat bath of Nosé–Hoover thermostat consists of only one imaginary particle, simulation systems achieve realistic constant-temperature condition. Therefore, the Nosé–Hoover thermostat has been commonly used as one of the most accurate and efficient methods for constant-temperature molecular dynamics simulations.
The Monte Carlo method for electron transport is a semiclassical Monte Carlo (MC) approach of modeling semiconductor transport. Assuming the carrier motion consists of free flights interrupted by scattering mechanisms, a computer is utilized to simulate the trajectories of particles as they move across the device under the influence of an electric field using classical mechanics. The scattering events and the duration of particle flight is determined through the use of random numbers.
Local elevation is a technique used in computational chemistry or physics, mainly in the field of molecular simulation. It was developed in 1994 by Huber, Torda and van Gunsteren to enhance the searching of conformational space in molecular dynamics simulations and is available in the GROMOS software for molecular dynamics simulation. The method was, together with the conformational flooding method, the first to introduce memory dependence into molecular simulations. Many recent methods build on the principles of the local elevation technique, including the Engkvist-Karlström, adaptive biasing force, Wang–Landau, metadynamics, adaptively biased molecular dynamics, adaptive reaction coordinate forces, and local elevation umbrella sampling methods. The basic principle of the method is to add a memory-dependent potential energy term in the simulation so as to prevent the simulation to revisit already sampled configurations, which leads to the increased probability of discovering new configurations. The method can be seen as a continuous variant of the Tabu search method.
The Hamiltonian Monte Carlo algorithm is a Markov chain Monte Carlo method for obtaining a sequence of random samples which converge to being distributed according to a target probability distribution for which direct sampling is difficult. This sequence can be used to estimate integrals with respect to the target distribution.
In statistics and physics, multicanonical ensemble is a Markov chain Monte Carlo sampling technique that uses the Metropolis–Hastings algorithm to compute integrals where the integrand has a rough landscape with multiple local minima. It samples states according to the inverse of the density of states, which has to be known a priori or be computed using other techniques like the Wang and Landau algorithm. Multicanonical sampling is an important technique for spin systems like the Ising model or spin glasses.
Classical nucleation theory (CNT) is the most common theoretical model used to quantitatively study the kinetics of nucleation.
In mathematics, the walk-on-spheres method (WoS) is a numerical probabilistic algorithm, or Monte-Carlo method, used mainly in order to approximate the solutions of some specific boundary value problem for partial differential equations (PDEs). The WoS method was first introduced by Mervin E. Muller in 1956 to solve Laplace's equation, and was since then generalized to other problems.
In mathematics and physics, surface growth refers to models used in the dynamical study of the growth of a surface, usually by means of a stochastic differential equation of a field.