Reliable kinetic Monte Carlo simulation based on random set sampling

Abstract

Kinetic Monte Carlo (KMC) method has been widely used in simulating rare events such as chemical reactions or phase transitions. Yet lack of complete knowledge of transitions and the associated rates is one major challenge for accurate KMC predictions. In this paper, a reliable KMC (R-KMC) mechanism is proposed in which sampling is based on random sets instead of random numbers to improve the robustness of KMC results. In R-KMC, rates or propensities are interval estimates instead of precise numbers. A multi-event algorithm based on generalized interval probability is developed. The weak convergence of the multi-event algorithm towards the traditional KMC is demonstrated with a generalized Chapman–Kolmogorov equation.

Appendix: Derivation of the generalized Chapman–Kolmogorov equation

Following the definition of interval derivative by Markov (1979), the derivative of a generalized interval function ${\mathbf{f}}(t) = \left[ {\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{f} (t),\bar{f}(t)} \right]$ is defined as

$$ \begin{aligned} &d{\mathbf{f}}(t)/dt: = \mathop {\lim }\limits_{\Updelta t \to 0} \left( {{\mathbf{f}}(t + \Updelta t) - \text{dual} {\mathbf{f}}(t)} \right)/\Updelta t \hfill \\ &\quad= \left[ {\mathop {\lim }\limits_{\Updelta t \to 0} \frac{{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{f} (t + \Updelta t) - \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{f} (t)}}{\Updelta t},\mathop {\lim }\limits_{\Updelta t \to 0} \frac{{\bar{f}(t + \Updelta t) - \bar{f}(t)}}{\Updelta t}} \right] \hfill \\ \end{aligned} $$

where dual is defined as in Eq. (2). Note that all boldface symbols in this paper are generalized intervals. We define the derivative of generalized interval probability ${\mathbf{p}}(x,t|y,t')$ with respect to time t as follows. With state variables $x,y \in R^{n}$,

$$ \frac{\partial }{\partial t}{\mathbf{p}}(x,t|y,t'): = \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ \begin{gathered} {\mathbf{p}}(x,t + \Updelta t|y,t') \hfill \\ - \text{dual} {\mathbf{p}}(x,t|y,t') \hfill \\ \end{gathered} \right\} $$

(15)

Because of the logic coherent constraint in generalized interval probability, we have

$$ \int {dz{\mathbf{p}}(z,t + \Updelta t|x,t)} = 1 $$

(16)

where state variable $z \in R^{n}$. With the Markovian property, this leads to

$$ {\mathbf{p}}(x,t + \Updelta t|y,t') = \int {dz{\mathbf{p}}(x,t + \Updelta t|z,t){\mathbf{p}}(z,t|y,t')} $$

(17)

Replace the first term in Eq. (15) with Eq. (17) and multiply the second term by Eq. (16). Eq. (15) now becomes

$$ \begin{aligned} &\frac{\partial }{\partial t}{\mathbf{p}}(x,t|y,t') \hfill \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {\int {dz\left[ \begin{gathered} {\mathbf{p}}(x,t + \Updelta t|z,t){\mathbf{p}}(z,t|y,t') \hfill \\ - \text{dual} {\mathbf{p}}(z,t + \Updelta t|x,t){\mathbf{p}}(x,t|y,t') \hfill \\ \end{gathered} \right]} } \right\} \hfill \\ \end{aligned} $$

(18)

Consider two subdomains $\left\| {x - z} \right\| \le \varepsilon$ and $\left\| {x - z} \right\| > \varepsilon$ separately where $\varepsilon$ is small enough. Eq. (18) can be regarded as

$$ \frac{\partial }{\partial t}{\mathbf{p}}(x,t|y,t') = {\mathbf{I}}_{{\left\| {x - z} \right\| \le \varepsilon }} + {\mathbf{I}}_{{\left\| {x - z} \right\| > \varepsilon }} $$

where ${\mathbf{I}}_{{\left\| {x - z} \right\| \le \varepsilon }}$ is the right-side integral in Eq. (18) within the small neighborhood that captures continuous diffusion processes, whereas ${\mathbf{I}}_{{\left\| {x - z} \right\| > \varepsilon }}$ is the one outside the neighborhood that represents jump processes.

Within the small neighborhood $\left\| {x - z} \right\| \le \varepsilon$, let $x - z = h$ and define ${\mathbf{f}}(x,h): = {\mathbf{p}}(x + h,t + \Updelta t|x,t){\mathbf{p}}(x,t|y,t')$. Then

$$ {\mathbf{p}}(z,t + \Updelta t|x,t){\mathbf{p}}(x,t|y,t') = {\mathbf{f}}(x, - h) $$

and

$$ \begin{aligned} &{\mathbf{p}}(x,t + \Updelta t|z,t){\mathbf{p}}(z,t|y,t') \hfill \\ &\quad= {\mathbf{p}}(x,t + \Updelta t|x - h,t){\mathbf{p}}(x - h,t|y,t') = {\mathbf{f}}(x - h,h). \hfill \\ \end{aligned} $$

Apply the Taylor alike expansion

$$ \begin{gathered} {\mathbf{f}}(x - h,h) = {\mathbf{f}}(x,h) - \text{dual} \sum\limits_{i = 1}^{n} {\partial {\mathbf{f}}/\partial x_{i} h_{i} } \hfill \\ + \frac{1}{2}\sum\limits_{i} {\sum\limits_{j}^{{}} {\partial^{2} {\mathbf{f}}/\partial x_{i} \partial x_{j} h_{i} h_{j} } } + O(\varepsilon^{3} ) \hfill \\ \end{gathered} $$

where

$$ \partial {\mathbf{f}}(x,h)/\partial x_{i} : = \mathop {\lim }\limits_{{\Updelta x_{i} \to 0}} \left\{ {\frac{{{\mathbf{f}}(x + \Updelta x_{i} ,h) - \text{dual} {\mathbf{f}}(x,h)}}{{\Updelta x_{i} }}} \right\} $$

and similarly.

$$ \begin{aligned} \frac{{\partial^{2} {\mathbf{f}}(x,h)}}{{\partial x_{i} \partial x_{j} }}: = \hfill \\ \mathop {\lim }\limits_{\begin{subarray}{l} \Updelta x_{i} \to 0 \\ \Updelta x_{j} \to 0 \end{subarray} } \left\{ {\frac{\begin{aligned} {\mathbf{f}}(x + \Updelta x_{i} + \Updelta x_{j} ,h) - \text{dual} {\mathbf{f}}(x + \Updelta x_{i} ,h) \hfill \\ - \text{dual} {\mathbf{f}}(x + \Updelta x_{j} ,h) + {\mathbf{f}}(x,h) \hfill \\ \end{aligned} }{{\Updelta x_{i} \Updelta x_{j} }}} \right\} \hfill \\ \end{aligned} $$

Then

$$ \begin{aligned} &{\mathbf{I}}_{{\left\| {x - z} \right\| \le \varepsilon }} \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {\int\limits_{{\left\| {x - z} \right\| \le \varepsilon }} {dz\left[ \begin{gathered} {\mathbf{p}}(x,t + \Updelta t|z,t){\mathbf{p}}(z,t|y,t') \hfill \\ - \text{dual} {\mathbf{p}}(z,t + \Updelta t|x,t){\mathbf{p}}(x,t|y,t') \hfill \\ \end{gathered} \right]} } \right\} \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {\int\limits_{\left\| h \right\| \le \varepsilon } {dh\left[ {{\mathbf{f}}(x - h,h) - \text{dual} {\mathbf{f}}(x, - h)} \right]} } \right\} \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {\int\limits_{\left\| h \right\| \le \varepsilon } {dh\left[ \begin{gathered} {\mathbf{f}}(x,h) - \text{dual} \sum\limits_{i = 1}^{n} {\frac{{\partial {\mathbf{f}}}}{{\partial x_{i} }}h_{i} } \hfill \\ + \frac{1}{2}\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {\frac{{\partial^{2} {\mathbf{f}}}}{{\partial x_{i} \partial x_{j} }}h_{i} h_{j} } } \hfill \\ - \text{dual} {\mathbf{f}}(x, - h) + O(\varepsilon^{3} ) \hfill \\ \end{gathered} \right]} } \right\} \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left[ {\int\limits_{\left\| h \right\| \le \varepsilon } {dh{\mathbf{f}}(x,h)} - \text{dual} \int\limits_{\left\| h \right\| \le \varepsilon } {dh{\mathbf{f}}(x, - h)} } \right] \\ &\quad- \text{dual} \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\sum\limits_{i = 1}^{n} {\frac{{\partial {\mathbf{f}}}}{{\partial x_{i} }}h_{i} } } \\ &\quad+ \frac{1}{2}\mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {\frac{{\partial^{2} {\mathbf{f}}}}{{\partial x_{i} \partial x_{j} }}h_{i} h_{j} } } } \\ &\quad+ \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}O(\varepsilon^{4} ) \\ &= - \text{dual} \sum\limits_{i = 1}^{n} {\left[ {\mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\frac{{\partial {\mathbf{f}}}}{{\partial x_{i} }}h_{i} } } \right]} \\ &\quad+ \frac{1}{2}\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {\left[ {\mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\frac{{\partial^{2} {\mathbf{f}}}}{{\partial x_{i} \partial x_{j} }}h_{i} h_{j} } } \right]} } \\ &\quad+ \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}O(\varepsilon^{4} ) \\ \end{aligned} $$

since $\int_{\left\| h \right\| \le \varepsilon } {dh{\mathbf{f}}(x,h)} = \int_{\left\| h \right\| \le \varepsilon } {dh{\mathbf{f}}(x, - h)}$. Furthermore,

$$ \begin{aligned} &\frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\frac{{\partial {\mathbf{f}}}}{{\partial x_{i} }}h_{i} } \\ \quad&= \frac{\partial }{{\partial x_{i} }}\left[ {\frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} {\mathbf{f}}dh} } \right] \\ \quad&= \frac{\partial }{{\partial x_{i} }}\left[ {{\mathbf{p}}(x,t|y,t')\frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} } \right] \\ \end{aligned} $$

Similarly,

$$ \begin{gathered} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {dh\frac{{\partial^{2} {\mathbf{f}}}}{{\partial x_{i} \partial x_{j} }}h_{i} h_{j} } = \hfill \\ \frac{{\partial^{2} }}{{\partial x_{i} \partial x_{j} }}\left[ {{\mathbf{p}}(x,t|y,t')\frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} h_{j} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} } \right] \hfill \\ \end{gathered} $$

In addition, $\mathop {\lim }\nolimits_{\Updelta t \to 0} O\left( {\varepsilon^{4} } \right)/\Updelta t = \mathop {\lim }\nolimits_{\Updelta t \to 0} O(\varepsilon^{3} )/1 = O(\varepsilon^{3} )$

We define

$$ {\mathbf{A}}_{i} (x,t): = \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} $$

(19)

$$ {\mathbf{B}}_{ij} (x,t): = \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} h_{j} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} $$

(20)

If $X(t) = (x_{1} (t), \ldots ,x_{n} (t))$ represents the state of the system at time $t$, then

$$ \begin{gathered} \int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} \hfill \\ = E\left[ {x_{i} (t + \Updelta t) - x_{i} (t)|X(t)} \right] \hfill \\ \end{gathered} $$

is the expected state value change in the $i$th direction. Therefore, Eq. (19) is interpreted as the state value change rate along one direction, known as the drift vector. Similarly,

$$ \begin{aligned} &\int\limits_{\left\| h \right\| \le \varepsilon } {h_{i} h_{j} dh{\mathbf{p}}(x + h,t + \Updelta t|x,t)} \hfill \\ &\quad= E\left[ {(x_{i} (t + \Updelta t) - x_{i} (t))(x_{j} (t + \Updelta t) - x_{j} (t))|X(t)} \right]. \hfill \\ \end{aligned} $$

Equation (20) is the combined area change rate in two directions, known as the diffusion matrix.

We now have

$$ \begin{aligned} &{\mathbf{I}}_{{\left\| {x - z} \right\| \le \varepsilon }} \\ &= - \text{dual} \sum\limits_{i = 1}^{n} {\frac{\partial }{{\partial x_{i} }}\left[ {{\mathbf{A}}_{i} (x,t){\mathbf{p}}(x,t|y,t')} \right]} \\ &\quad+ \frac{1}{2}\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {\frac{\partial }{{\partial x_{i} \partial x_{j} }}\left[ {{\mathbf{B}}_{ij} (x,t){\mathbf{p}}(x,t|y,t')} \right]} } + O(\varepsilon^{3} ) \\ \end{aligned} $$

for the diffusion process.

For the jump process

$$ \begin{aligned} &{\mathbf{I}}_{{\left\| {x - z} \right\| > \varepsilon }} \\ &= \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {\int\limits_{{\left\| {x - z} \right\| > \varepsilon }} {dz\left[ \begin{gathered} {\mathbf{p}}(x,t + \Updelta t|z,t){\mathbf{p}}(z,t|y,t') \hfill \\ - \text{dual} {\mathbf{p}}(z,t + \Updelta t|x,t){\mathbf{p}}(x,t|y,t') \hfill \\ \end{gathered} \right]} } \right\} \\ &= \int\limits_{{\left\| {x - z} \right\| > \varepsilon }} {dz\left[ {\mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {{\mathbf{p}}(x,t + \Updelta t|z,t)} \right\}{\mathbf{p}}(z,t|y,t')} \right]} \\ &\quad- \text{dual} \int\limits_{{\left\| {x - z} \right\| > \varepsilon }} {dz\left[ {\mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {{\mathbf{p}}(z,t + \Updelta t|x,t)} \right\}{\mathbf{p}}(x,t|y,t')} \right]} \\ \end{aligned} $$

We define

$$ {\mathbf{W}}(y|x,t): = \mathop {\lim }\limits_{\Updelta t \to 0} \frac{1}{\Updelta t}\left\{ {{\mathbf{p}}(y,t + \Updelta t|x,t)} \right\} $$

(21)

as the rate of transition from $x$ to $y$. Then

$$ \begin{aligned} {\mathbf{I}}_{{\left\| {x - z} \right\| > \varepsilon }} &= \int\limits_{{\left\| {x - z} \right\| > \varepsilon }} {dz\left[ {{\mathbf{W}}(x|z,t){\mathbf{p}}(z,t|y,t')} \right]} \\ &\quad- \text{dual} \int\limits_{{\left\| {x - z} \right\| > \varepsilon }} {dz\left[ {{\mathbf{W}}(z|x,t){\mathbf{p}}(x,t|y,t')} \right]} \\ \end{aligned} $$

As $\varepsilon \to 0$, the subdomain $\left\| {x - z} \right\| > \varepsilon$ becomes the complete domain.

Therefore, with the consideration of both $\left\| {x - z} \right\| \le \varepsilon$ and $\left\| {x - z} \right\| > \varepsilon$, the generalized differential Chapman–Kolmogorov equation is

$$ \begin{aligned} &\frac{\partial }{\partial t}{\mathbf{p}}(x,t|y,t') = \hfill \\ &\quad- \text{dual} \sum\limits_{i = 1}^{n} {\frac{\partial }{{\partial x_{i} }}\left[ {{\mathbf{A}}_{i} (x,t){\mathbf{p}}(x,t|y,t')} \right]} \hfill \\ &\quad+ \frac{1}{2}\sum\limits_{i = 1}^{n} {\sum\limits_{j = 1}^{n} {\frac{{\partial^{2} }}{{\partial x_{i} \partial x_{j} }}\left[ {{\mathbf{B}}_{ij} (x,t){\mathbf{p}}(x,t|y,t')} \right]} } \hfill \\ &\quad+ \int_{{}} {dz{\mathbf{W}}(x|z,t){\mathbf{p}}(z,t|y,t')} \hfill \\ &\quad- \text{dual} \int_{{}} {dz{\mathbf{W}}(z|x,t){\mathbf{p}}(x,t|y,t')} \hfill \\ \end{aligned} $$

(22)