If you've ever tried to simulate a deposition process for film growth, you'll quickly find that a lot of different events take place. Moreover, these events occur at different timescales. In a deposition process, for example, the particles/species are arriving at the surface at a much faster timescale than, say, surface diffusion. To put it in context, say we have a 1 mm$^2$ area, and our particle flux is $10^{21} \frac{\#}{\text{m}^2 \cdot \text{s}}$, while the diffusion coefficient is $10^{-4} \frac{\text{cm}^2}{\text{s}}$. If we assume some desired number of particles to deposit to create a film on this area, then the timescales to see events are something like microseconds for deposition and thousands of seconds for diffusion.
These timescales are vastly different and hence will need to be handled differently. In the deposition process, we could treat it directly using atomistic simulations, but it turns out that for most practical fluxes, the timescales are just out of reach. To solve this problem, we can turn to a technique called kinetic Monte Carlo (KMC). For diffusion, it's clear that we need to use KMC or solve mass transport equations.
Kinetic Monte Carlo (KMC)
Most people are familiar with Monte Carlo methods. These are stochastic approaches to solving problems where, instead of computing the solutions directly, we sample from probable events. The probability is determined by the probability mass or density function. Monte Carlo methods can be used to solve all kinds of problems, for example performing integration. Say we have a function we want to integrate but don't know a deterministic algorithm or analytical solution. We can use Monte Carlo to randomly sample the area under the curve and estimate a definite integral. The ultimate accuracy depends on how much sampling you decide to perform.
Monte Carlo methods let us perform random sampling to estimate the behavior of a system or process; if we think back to our integral example, we can think of the system or process being described by a high-dimensional probability density or mass function that we are trying to integrate to get the expected behavior (i.e., value).
So what does it mean when we add the prefix "kinetic" to Monte Carlo? It simply designates that the process has a time component to the evolution. Think of a chemical reaction or diffusive dynamics. In KMC, the events that are described are derived from rates, $\nu$, which are usually in inverse time units. The probability of an event is determined by the rate, and this in turn informs the elapsed time. To expand, let's look at the math.
Say we have some event, $\Omega_i$, that corresponds to a diffusive hop. There is some energy required to overcome for the diffusive hop to occur, given as $E_{a}^i$. We can describe the rate as $\nu_i = \nu_0^i \exp\left(-\frac{E_{a}^i}{k_B T}\right)$, where $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\nu_0^i$ is the prefactor (attempt frequency) for event $i$. So now assume we have $N$ events that describe some diffusive hop to a different site/location. We can describe the probability of any one event occurring by:
\begin{equation} P(\Omega_i; E_a^i, \nu_0^i) = \frac{\nu_0^i \exp\left(-\dfrac{E_{a}^i}{k_B T}\right)}{\sum\limits_{j=1}^N \nu_0^j \exp\left(-\dfrac{E_{a}^j}{k_B T}\right)} \label{eq:probability} \end{equation}
I intentionally include that the probability is parameterized based on selecting the activation energy and the prefactor. However, we could easily consider this as a conditional probability if we have some knowledge of the range of values for these terms.
Thus, you have a way to sample from the event catalog now. So how do you define time based on the sampling? If an event $\Omega_i$ is selected, then we can calculate the time increment associated with this event. In kinetic Monte Carlo, we generate a second random number to calculate the time increment. Specifically, the time increment $\Delta t$ is given by:
\begin{equation} \Delta t = -\dfrac{1}{R_N} \ln r \end{equation}
where $r$ is a uniform random number between 0 and 1, and $R_N$ is the cumulative rate for all events/processes, i.e., $R_N = \sum\limits_{j=1}^N \nu_j$. This time increment reflects the stochastic time elapsed before the occurrence of event $\Omega_i$.
It's a pretty simple formalism but very powerful to simulate kinetic processes at the right timescales. The challenge is always two-fold:
-
Treating a global event table with all process timescales and rates can lead to the need for very large sampling trials. In other words, if we treat the deposition and diffusive processes on the same event table, diffusive events will have a lower chance of being sampled. The easiest way to address this is simply by decoupling them and having two KMC times simultaneously evolve.
-
A priori knowledge of what should happen, i.e., this atom deposits at this site and reaction $\alpha$ occurs if it has species A but does reaction $\beta$ if it has species B. For well-understood systems that have simple species, this is manageable, but it should be clear that the combinatorics of events for multi-species deposition and diffusion becomes prohibitive, and more importantly, we have many model parameters to decide on.
Lattice KMC
One way to simplify the KMC approach and to simulate crystalline materials is to use a lattice formalism. This means we can describe lattice points where events can occur and thus it somewhat restricts what can happen. In particular, diffusion hops are constrained to occur between neighboring lattice sites. This works well for crystalline materials with a single basis, but one may need to fall back on mesh methods to treat more complicated crystals.
🧠 Thinking out Loud
I haven't really thought too much about meshing approaches, but the idea would be to define nested lattices where each basis point is defined. One could then define the events that are allowed within a lattice and between lattices. Take a Cesium Chloride crystal; we would have two interpenetrating simple cubic lattices.
Some of the challenges with lattice KMC are that you can't treat amorphous materials and grain boundaries. And it gets complicated for multiple species in my view.
You're also still plagued with the fact that you need to enumerate all the events and assign the rates. Maybe not too difficult for a monoatomic species, where you can find a good amount of information for the surface deposition/diffusion and bulk diffusion.
On-the-Fly Off-Lattice KMC
First, imagine we remove the lattice constraint and let points (representing atoms or molecules) occupy any point in space. Then, say we have a way to determine for a species the value of $E_a^i$ for some motion it could have in space. Think of a surface atom hopping on top of another surface atom. If you could find and enumerate all the possible motions such an atom could have, you would come up with a list of activation energies and therefore rates. Thus, you would have cataloged an event table you could then sample from to perform KMC. This is on-the-fly KMC.
Basically, you need to have a model/description for the energy barrier associated with different dynamical processes. We do this through the potential energy surface (PES) model, which is the functional surface describing, for a given configuration of atoms in space, what is the potential energy and the associated gradients (i.e., forces).
Once you have a way to describe the PES, then we need a methodology to sample different pathways (e.g., reaction coordinates) to determine what pathways through space an atom(s) could take and the associated energetic barrier. Then, through transition state theory, we can determine the rate of each pathway; that is, we get an Arrhenius expression for the rate of a reaction coordinate:
\begin{equation} \nu_i = \nu_0 \exp\left(-\dfrac{E_{a}^i}{k_B T}\right) \end{equation}
Challenges with Transition State Theory
The two major hurdles with transition state theory are:
-
You need to have an accurate PES description for the chemical system you're trying to describe. This usually requires force-field or interatomic potential development (see my posts discussing such topics).
-
The pathway along each reaction coordinate either needs an intelligent way to search around an initial point or requires an initial and final point to be connected. The former are called single-ended methods, and techniques like the dimer method or growing string method are used. These essentially find the saddle points of the PES, and thus the final state is just the "downhill" path from the saddle point. The other method is dual-ended, where the initial and final coordinates are known, but possible paths are not. We then can find the possible barriers and usually only tabulate the one that is the lowest.
Most of the effort is always centered around proper PES description and setting up the TST calculations, which can be truly challenging. This is why we probably don't see many on-the-fly KMC codes in practice. The only one that I came across was EON from the Henkelman group at UT Austin [1]. It uses a client-server approach where the server runs the TST calculations and then hands the event table back to the KMC client.
Putting It All Together
In practice, I think the outline for an on-the-fly off-lattice KMC code that simulates deposition and diffusion would be something like:
As this blog was more a discussion on KMC, just so I can jot down my thoughts, I'll leave it at this for now. If you want to read more on lattice KMC, you can read the chapter in the book by LeSar on it; it's a pretty good starting point [2].
References
[1] S.T. Chill, M. Welborn, R. Terrell, L. Zhang, J.-C. Berthet, A. Pedersen, H. Jónsson, G. Henkelman, EON: software for long time simulations of atomic scale systems, Modelling Simul. Mater. Sci. Eng. 22 (2014) 055002. DOI.
[2] R. LeSar, Introduction to Computational Materials Science: Fundamentals to Applications, Cambridge University Press, Cambridge; New York, 2013. DOI
