3.3 Monte Carlo (MC), quantum Monte Carlo (QMC), and kinetic Monte Carlo (KMC) methods
In the limit of steady-state, the master equation in (Eq. 1) can be written in matrix form as WP = P or in the familiar Einstein’s convention of wijPej = Pei , where the summation over the repeated index is implicitly assumed. It is important to recognize is that W is the entire matrix and wij is the ijth element of the matrix. Note that wij is the transition probability of migrating from microstate j to I, consistent with the definition of w in (Eq. 1). Similarly, P is the entire vector of probabilities of each microstate, and the ith element of P is Pi, the probability to access microstate i. Note that here, Pei is the equilibrium distribution. More generally, if we start with a non-equilibrium state P(1), here, (1) is the initial time and the system transitions to later times and is tracked by (2); (3), etc., then: WP(1) = P(2);WP(2) = P(3);WP(3) = P(4); …; WP(n) = P(n+1), and as n becomes large, P(n) = P(n+1) = Pe. In a Monte Carlo simulation, we simulate the system by sampling accessible microstates according to their equilibrium distribution Pe, e.g., as given by the Boltzmann distribution or the appropriate equivalent distribution in different ensembles (for thermodynamic systems at equilibrium). To achieve this task, we need to choose a W or all wij that make up the W, such that Pei= exp(-Ei/kBT)/[Σjexp(-Ej/kBT)] is satisfied. Metropolis et al. [30] recognized that this could be achieved by choosing wij that satisfy equation wijPej = Pei by imposing a stronger criterion, namely, Pemwnm=Penwmn, which leads to the Metropolis Monte Carlo method for sampling microstates of a classical system.
Quantum Monte Carlo (QMC) techniques provide a direct and potentially efficient means for solving the many-body Schrödinger equation of quantum mechanics [31]. The simplest quantum Monte Carlo technique, variational Monte Carlo (VMC), is based on a direct application of Monte Carlo integration to calculate multi-dimensional integrals of expectation values such as the total energy. Monte Carlo methods are statistical, and a key result is that the value of integrals computed using Monte Carlo converges faster than by using conventional methods of numerical quadrature, once the problem involves more than a few dimensions. Therefore, statistical methods provide a practical means of solving the full many-body Schrödinger equation by direct integration, making only limited and well-controlled approximations.
The kinetic Monte Carlo (KMC) method is a Monte Carlo method computer simulation intended to simulate the time evolution of processes that occur with known transition rates among states (such as chemical reactions or diffusion transport). It is essential to understand that 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 [32]. From a mathematical standpoint, solving the master equation for such systems is impossible owing to the combinatorially large number of accessible microstates, even considering that the limited number of accessible states renders the transition probability matrix sparse. The Gillespie algorithm provides an ingenious way out of this issue. The practical idea behind KMC is not to attempt to deal with the entire matrix, but instead to generate stochastic trajectories that propagate the system from state to state (i.e., a Markovian sequence of discrete hops to random states happening at random times). From this, the correct time evolution of the probabilities Pi(t) is then obtained by ensemble averaging over these trajectories. The KMC algorithm does so by selecting elementary processes according to their probabilities to fire, followed by an updating of the time.