首页

动力学蒙特卡洛方法(KMC)及相关讨论

星期二, 2010-05-11 01:05 — satchel1979

动态模拟在目前的计算科学中占据着非常重要的位置。随着计算能力和第一原理算法的发展，复杂的动态参数（扩散势垒、缺陷相互作用能等）均可利用第一原理计算得出。因此，部分复杂的体系动态变化，如表面形貌演化或辐射损伤中缺陷集团的聚合-分解演变等，已可以较为精确的予以研究。KMC——动力学蒙特卡洛方法（kinetic Monte Carlo）原理简单，适应性强，因此在很多情况下都是研究人员的首选。此外，KMC在复杂体系或复杂过程中的算法发展也非常活跃。本文试图介绍KMC方法的基础理论和若干进展。

KMC方法基本原理

在原子模拟领域内，分子动力学(molecular dynamics, MD)具有突出的优势。它可以非常精确的描述体系演化的轨迹。一般情况下MD的时间步长在飞秒($ 10^{-15} $s)量级，因此足以追踪原子振动的具体变化。但是这一优势同时限制了MD在大时间尺度模拟上的应用。现有的计算条件足以支持MD到10 ns，运用特殊的算法可以达到10 $ \mu $s的尺度。即便如此，很多动态过程，如表面生长或材料老化等，时间跨度均在s以上，大大超出了MD的应用范围。有什么方法可以克服这种局限呢？

当体系处于稳定状态时，我们可以将其描述为处于$ 3N $维势能函数面的一个局域极小值（阱底）处。有限温度下，虽然体系内的原子不停的进行热运动，但是绝大部分时间内原子都是在势能阱底附近振动。偶然情况下体系会越过不同势阱间的势垒从而完成一次“演化”，这类小概率事件才是决定体系演化的重点。因此，如果我们将关注点从“原子”升格到“体系”，同时将“原子运动轨迹”粗化为“体系组态跃迁”，那么模拟的时间跨度就将从原子振动的尺度提高到组态跃迁的尺度。这是因为这种处理方法摈弃了与体系穿越势垒无关的微小振动，而只着眼于体系的组态变化。因此，虽然不能描绘原子的运动轨迹，但是作为体系演化，其“组态轨迹”仍然是正确的。此外，因为组态变化的时间间隔很长，体系完成的连续两次演化是独立的，无记忆的，所以这个过程是一种典型的马尔可夫过程(Markov process)，即体系从组态$ i $到组态$ j $，$ i\rightarrow j $这一过程只与其跃迁速率$ k_{ij} $有关。如果精确地知道$ k_{ij} $，我们便可以构造一个随机过程，使得体系按照正确的轨迹演化。这里``正确''的意思是某条给定演化轨迹出现的几率与MD模拟结果完全一致(假设我们进行了大量的MD模拟，每次模拟中每个原子的初始动量随机给定)。这种通过构造随机过程研究体系演化的方法即为动力学蒙特卡洛方法(kinetic Monte Carlo, KMC) [1]。

指数分布与KMC的时间步长

在KMC模拟中，构造呈指数分布的随机数是一个相当重要的步骤。这一节中我们对此进行讨论。

因为体系在势能面上无记忆的随机行走，所以任意单位时间内，它找到跃迁途径的概率不变，设为$ k_{\rm tot} $。因此在$ [t,t+\Delta t) $区间内，体系不发生跃迁的概率为

$ F_{\rm stay}(\Delta t)=1-k_{\rm tot}\Delta t+\mathcal{O}(\Delta t)^2 $

类似的，在$ [t,t+2\Delta t) $区间内，体系不发生跃迁的概率为

$ F_{\rm stay}(2\Delta t)=\bigg(1-k_{\rm tot}\Delta t+\mathcal{O}(\Delta t)^2\bigg)^2=1-2k_{\rm tot}\Delta t+\mathcal{O}(\Delta t)^2 $

以此类推，当$ \tau = K\Delta t $时，在$ [t,t+\tau) $区间内，体系不发生跃迁的概率为

$ F_{\rm stay}(\tau)=\bigg(1-k_{\rm tot}\frac{\tau}{K}\Delta t+\mathcal{O}(K^{-2})\bigg)^K $

因此，当$ \tau $趋于$ \infty $时，体系不发生跃迁的概率为

$ F_{\rm stay}(\tau) = \lim_{\tau\rightarrow\infty}\bigg(1-k_{\rm tot}\frac{\tau}{K}\Delta t+\mathcal{O}(K^{-2})\bigg)^K=\exp(-k_{\rm tot}\tau) $      (1)

这一行为类似于原子核的衰变方程。从方程(1)我们可以得到单位时间内体系跃迁概率$ p(t) $。从方程(1)的推导过程可以看出体系的跃迁概率是一个随时间积累的物理量，因此$ p(t) $对时间积分到某一时刻$ t' $必然等于$ 1-F_{\rm stay}(t') $，也即$ p(t)=\partial\big(1-F_{\rm stay}(t')\big)/\partial t $。因此我们立即可以得到 [1]

$ p(t)=k_{\rm tot}\exp(-k_{\rm tot}t) $       (2)

$ k_{\rm tot} $是体系处于态$ i $时所有可能的跃迁途径的速率$ k_{ij} $之和，即

$ k_{\rm tot} = \sum_{j}k_{ij} $       (3)

对于每个具体的跃迁途径$ k_{ij} $，上述讨论均成立。因此，我们可以定义单位时间内体系进行$ i\rightarrow j $跃迁的概率$ p_{ij}(t) $为

$ p_{ij}(t)=k_{ij}\exp(-k_{ij}t) $       (4)
 

单位时间内体系的跃迁概率呈指数分布这一事实说明KMC的时间步长$ \delta t $也应是指数分布。因此我们需要产生一个指数分布的随机数序列。这一点可以非常容易的通过一个(0,1]平均分布的随机数序列$ {r} $转化得到：

$ r=1-F_{\rm stay}(k_{\rm tot}\delta t) $

从而

$ \delta t = -\frac{1}{k_{\rm tot}}\ln(1-r) = -\frac{1}{k_{\rm tot}}\ln(r) $       (5)

最后一步是因为$ 1-r $和$ r $的分布相同。$ \delta t $也可以通过上述步骤从方程(4)得到。

计算跃迁速率

过渡态理论(TST)

$ k_{ij} $决定了KMC模拟的精度甚至准确性。为避开通过原子轨迹来确定$ k_{ij} $的做法(这样又回到了MD的情况)，一般情况下采用过渡态理论(transition state theory, TST)进行计算 [2]。在TST中，体系的跃迁速率决定于体系在鞍点处的行为，而平衡态(势阱)处的状态对其影响可以忽略不计。如果大量的相同的体系组成正则系综，则在平衡状态下体系在单位时间内越过某个垂直于$ i\leftrightarrow j $跃迁途径的纵截面的流量即为$ k_{ij} $。简单起见，假设有大量相同的一维双组态(势阱)体系，平衡状态下鞍点所在的假想面(对应于流量最小的纵截面)为$ x=q $，则TST给出该体系从组态A迁出到B的速率为 [5,6]

$ k_{{\rm A}\rightarrow{\rm B}}=\frac{1}{2}\langle\delta(x-q)|\dot{x}|\rangle_{\rm A} $       (6)

方程(6)中$ \langle\cdots\rangle_{\rm A} $表示在组态A所属态空间里对正则系综的平均。$ \frac{1}{2} $表示只考虑体系从组态A迁出而不考虑迁入A的情况(后一种情况体系也对通过纵截面的流量有贡献)。根据普遍公式

$ \langle\hat{O}\rangle = \frac{\int\hat{O}e^{-\frac{H}{k_BT}}dxdp}{\int e^{-\frac{H}{k_BT}}dxdp} $

设体系的哈密顿量$ H $为$ p^2/2m+V(x) $，即可分解为动能和势能，同时设粒子坐标$ x\leq q $时体系处于组态A。则方程(6)可写为

$ k_{{\rm A}\rightarrow{\rm B}}=\frac{\frac{1}{2}\int^{\infty}_{-\infty}dp\int^{q+\epsilon}_{-\infty}dx[\delta(x-q)|\dot{x}|\exp\{[p^2/2m+V(x)]/k_BT\}]}{\int^{\infty}_{-\infty}dp\int^{q+\epsilon}_{-\infty}dx\exp\{[p^2/2m+V(x)]/k_BT\}} $
     $ =\frac{1}{2}\Bigg(\frac{\int^{\infty}_{-\infty}|\dot{x}|e^{-(p^2/2m)/k_BT}dp}{\int^{\infty}_{-\infty}e^{-(p^2/2m)/k_BT}dp}\Bigg)\frac{\int^{q+\epsilon}_{-\infty}\delta(x-q)e^{-V(x)/k_BT}dx}{\int^{q+\epsilon}_{-\infty}e^{-V(x)/k_BT}dx} $
     $ =\frac{1}{2}\langle|\dot{x}|\rangle\langle\delta(x-q)\rangle_{\rm A} $
     $ =\frac{1}{2}\Big(\frac{2k_BT}{\pi m}\Big)^{1/2}\langle\delta(x-q)\rangle_{\rm A} $       (7)

上式中无限小量$ \epsilon $是为了将$ \delta $函数全部包含进去。最后一项对于$ \delta $函数的系综平均可以直接通过Metropolis Monte Carlo方法计算出来：计算粒子落在$ [q-\omega,q+\omega] $范围内的次数相对于Metropolis行走总次数的比例$ f_B $。方程(7)最后等于

$ k_{{\rm A}\rightarrow{\rm B}}=\frac{1}{2}\Big(\frac{2k_BT}{\pi m}\Big)^{1/2}\Big(\frac{f_B}{\omega}\Big) $        (8)

将上述讨论扩展到3维情况非常直接，这里只给出结果，详细讨论请参阅文献 [5]:

$ k_{{\rm A}\rightarrow{\rm B}}=\frac{1}{2}\Big(\frac{2k_BT}{\pi m}\Big)^{1/2}\langle\delta[f(\mathbf{R})]|\nabla f|\rangle_{\rm A} $       (9)

其中$ f(\mathbf{R}) $是纵截面方程，$ |\nabla f| $代表3维情况中粒子流动方向与截面$ f $法向不平行对于计数的影响。

简谐近似下的过渡态理论(hTST)

虽然上一节已经给出了TST计算跃迁速率的方法，但是在具体工作中，$ k_{ij} $更多地是利用简谐近似下的过渡态理论(harmonic TST, hTST)通过解析表达式给出。根据TST，跃迁速率$ k_{ij} $为 [3]

$ k_{ij}=\frac{k_BT}{h}\exp[-\Delta F_{ij}/(k_BT)] $       (10)

其中$ \Delta F_{ij} $为在跃迁$ i\rightarrow j $中体系在鞍点和态$ i $处的自由能之差

$ \Delta F_{ij}=E^{\rm sad}_{ij}-TS^{\rm vib}_{ij}-[E_i-TS_i]=\Delta E_{ij}-T\Delta S^{\rm vib}_{ij} $

将上式代入方程(10)，可以得到

$ k_{ij}=\frac{k_BT}{h}\exp(-\Delta E_{ij}/k_BT)\exp(\Delta S^{\rm vib}_{ij}/k_B) $       (11)

hTST认为体系在稳态附近的振动可以用谐振子表示，因此其配分函数是经典谐振子体系的配分函数。分别写出体系在态$ i $和鞍点处的配分函数$ Z^0 $和$ Z^{\rm sad} $：

$ Z^0 = \bigg(\frac{k_BT}{h}\bigg)^{3N}\thickspace\overset{3N}{\underset{i=1}{\prod}}\frac{1}{\nu_{ij}^0} $
$ Z^{\rm sad} = \bigg(\frac{k_BT}{h}\bigg)^{3N-1}\thickspace\overset{3N-1}{\underset{i=1}{\prod}}\frac{1}{\nu_{ij}^{\rm sad}} $

根据Boltzmann公式，

$ S = k_B\ln Z $       (12)

并将配分函数代入，则方程(11)得

$ k_{ij}=\frac{\overset{3N}{\underset{i=1}{\prod}}\nu_{ij}^0}{\overset{3N-1}{\underset{i=1}{\prod}}\nu_{ij}^{\rm sad}}\exp\Big(-\frac{\Delta E_{ij}}{k_BT}\Big) $       (13)

方程(13)在通常的文献上经常可以见到。声子谱可以通过Hessian矩阵对角化或者密度泛函微扰法(DFPT)求出，而$ \Delta E_{ij} $就是$ i\rightarrow j $的势垒，可以通过NEB或者drag方法求出。因此，方程(13)保证了可以通过原子模拟(MD或者DFT方法)解析地求出$ k_{ij} $。事实上这个方程有两点需要注意。首先虽然方程(10)中出现了普朗克常数$ h $，但是在最终结果中$ h $被抵消了。这是因为TST本质上是一个经典理论，所以充分考虑了统计效应后$ h $不会出现 [1]。其次，方程(13)表明对于每一个跃迁过程，鞍点处的声子谱应该单独计算。这样会大大增加计算量，因此在绝大部分计算中均设前置因子为常数，不随跃迁过程而变化。具体数值取决于体系，对于金属而言，一般取$ \sim 10^{12} $ Hz。

KMC几种不同的实现算法
 

点阵映射
 

到目前为止，进行KMC模拟的所有理论基础均已具备。但是前面所进行的讨论并没有联系到具体的模型。KMC在固体物理中的应用往往利用点阵映射将原子与格点联系起来。从而将跃迁(事件)具象化为原子$ \leftrightarrow $格点关系的变化。比如空位(团)/吸附原子(岛)迁移等等。虽然与实际情况并不完全一致，但这样做在很多情况下可以简化建模的工作量，而且是非常合理的近似。很多情况下体系中的原子虽然对理想格点均有一定的偏离，但是并不太大($ \sim 0.01a_0 $)，因此这种原子$ \leftrightarrow $点阵映射是有效的。这种做法的另一个好处是可以对跃迁进行局域化处理。每条跃迁途径只与其近邻的体系环境有关，这样可以极大的减少跃迁途径的数目，从而简化计算 [1]。需要指出的是，这种映射对于KMC模拟并不是必须的。比如化学分子反应炉或者生物分子的生长等等，这些情况下根本不存在点阵。

无拒绝方式
 

KMC的实现方法有很多种，这些算法大致可以分为拒绝(rejection)和无拒绝(rejection-free)两种范畴。每种范畴之下还有不同的实现方式。本文只选择几种最为常用的方法加以介绍。

I. 直接法

直接法(direct method)是最常用的一种KMC算法，其效率非常高。每一步只需要产生两个在$ (0,1] $之间平均分布的随机数$ r_1 $和$ r_2 $。其中$ r_1 $被用来选定跃迁途径，$ r_2 $确定模拟的前进时间。设体系处于态$ i $，将每条跃迁途径$ j $想象成长度与跃迁速率$ k_{ij} $成正比的线段。将这些线段首尾相连。如果$ r_1k_{\rm tot} $落在线段$ j_k $中，这个线段所代表的跃迁途径$ j_k $就被选中，体系移动到态$ j_k $，同时体系时间根据方程(5)前进。总结其算法如下：

  1. 根据方程(4)计算体系处于态$ i $时的总跃迁速率$ k_{\rm tot} = \sum_{j}k_{ij} $；
  2. 选择随机数$ r_1 $；
  3. 寻找途径$ j_k $，满足$ \overset{j_k-1}{\underset{j=1}{\sum}}k_{ij}\leq r_1k_{\rm tot}<\overset{j_k}{\underset{j=1}{\sum}}k_{ij} $；
  4. 体系移动到态$ j_k $，同时模拟时间前进$ \delta t = -\frac{1}{k_{\rm tot}}\ln(r_2) $；
  5. 重复上述过程。

需要指出的是，虽然一般步骤4中的$ \delta t $根据方程(5)生成，但是如果将其换为$ \delta t = \frac{1}{k_{\rm tot}} $并不会影响模拟结果。在文献[5]和[6]中均采用这种方式。

II. 第一反应法

第一反应法(first reaction method, FRM)在思路上比直接法更为自然。前面说过，对于处于稳态$ i $的体系而言，它可以有不同的跃迁途径$ j $可以选择。每条途径均可以根据方程(4)给出一个指数分布的"发生时间"$ \delta t_{ij} $，也即从当前算起$ i\rightarrow j $第一次发生的时间。然后从$ \{\delta t_{ij}\} $中选出最小值(最先发生的"第一反应")，体系跃迁到相应的组态$ j_{\rm min} $，模拟时间相应地前进$ \delta t_{ij_{\rm min}} $。总结其算法如下：

  1. 设共有$ M $条反应途径，生成$ M $个随机数$ r_1, r_2, \cdots r_M $；
  2. 根据公式$ \delta t_{ij} = -\frac{1}{k_{ij}}\ln(r_j) $，给出每条路径的预计发生时间；
  3. 找出$ \{\delta t_{ij}\} $的最小值$ \delta t_{ij_{\rm min}} $；
  4. 体系移动到态$ j_{\rm min} $，同时模拟时间前进$ \delta t_{ij_{\rm min}} $；
  5. 重复上述过程。

可以看出，这种算法的效率比直接法低下，因为每一步KMC模拟需要生成$ M $个随机数。通常情况下KMC模拟需要$ 10^7 $步来达到较好的统计性质，如果每一步都需要生成$ M $个随机数，则利用这种方法需要一个高质量的伪随机数发生器，这一点在$ M $比较大时尤为重要。

III. 次级反应法

次级反应法(next reaction method, NRM)是FRM方法的一种衍生方法，其核心思想是假设体系的一次跃迁并不会导致处于新态的体系对于其他跃迁途径的舍弃(比如充满可以发生$ M $种化学反应的分子，第一种反应发生并不会造成别的反应物的变化)，这样体系还可以选择$ \{\delta t_{ij}\} $中的次小值$ \delta t_{ij_{\rm 2nd}} $，从而跃迁到态$ j_{\rm 2nd} $，模拟时间前进$ \delta t_{ij_{\rm 2nd}}-\delta t_{ij_{\rm min}} $。如果这次跃迁还可以满足上述假设条件，再重复上述过程。理想情况下，平均每一步KMC模拟只需要生成1个随机数。这无疑会大大提高效率以及时间跨度。但是实际上NRM的假设条件很难在体系每次跃迁之后都得到满足，在固体物理的模拟中尤其如此，因此其应用范围集中于研究复杂化学环境下的反应过程。

试探-接受/拒绝方式
 

这一大类算法虽然在效率上不如直接法，但是它们所采用的试探-接受/拒绝在形式上更接近Metropolis MC方法，而且可以很方便的引入恒定步长，即$ \delta t $固定。因此有必要进行详细的介绍。

IV. 选择直接法

选择直接法在决定体系是否跃迁方面和Metropolis MC方法形式上非常相像，均是通过产生随机数和预定的阈值比较决定事件是否被采纳。具体算法如下：

  1. 设共有$ M $条反应途径，选择反应速率最大值$ k^{\rm max} $，设为$ \hat{k} $。生成在$ [0,M) $均匀分布的随机数$ r $；
  2. 设$ j={\rm Int}(r)+1 $；
  3. 如果$ j-r $ < $ k_{ij}/\hat{k} $，则体系跃迁至新态$ j $，否则保持在态$ i $；
  4. 模拟时间前进$ \delta t = -\frac{1}{k_{\rm tot}}\ln(r) $；
  5. 重复上述过程。

这种方法的长处在于每一步只需要生成一个随机数。但是缺点也很明显，对于反应速率相差太大，尤其是只有一个低势垒途径(与其他途径相比$ k_{ij} $过大)的体系来讲，这种方法的效率会非常低下。某些情况下，这种低效率问题可以通过如下方法改进：将全部途径按照$ k_{ij} $的大小分为几个亚组，每个亚组选定一个上限$ \hat{k}^{n} $。但是这一步骤在整个KMC模拟过程中可能需要重复很多次，因此并不能完全解决问题。事实上低势垒在KMC中是个普遍的问题。这一点在后面还要简要提及。

V. 恒定步长法

与上述四种方法不同，恒定步长法(constant time step method, CTSM)中体系的前进时间是个给定的参数\cite{dawnkaski}。在理想情况下，CTSM与直接法效率相同，每一步只需产生两个随机数。具体算法如下：

  1. 给定恒定时间步长$ \delta t $；
  2. 将所有途径$ j $(共有$ M $个)设为长度恒为$ 1/M $的线段，生成在$ [0,1] $均匀分布的随机数$ r_1 $，选择途径$ j={\rm Int}(r_1M)+1 $；
  3. 生成在$ [0,1] $均匀分布的随机数$ r_2 $，如果$ r_2 $ < $ k_{ij}\delta t $，则体系跃迁至新态$ j $，否则保持在态$ i $；
  4. 模拟时间前进$ \delta t $；
  5. 重复上述过程。

实际模拟中，$ \delta t $需要满足(1)小于$ \delta t_{ij_{\rm min}} $(见"第一反应法")，以及(2)对于$ k_{ij} $最大的途径，接受率大致在0.5。其中第一个条件保证了所有的迁移途径发生概率都小于1，第二个条件则保证体系演化的效率不会过于低下。CTSM是非常行之有效的一类KMC算法，但是选择$ \delta t $时需要特别的注意以保证效率。$ \delta t $决定于具体体系以及模拟温度。这在一定程度上增加了CTSM的实现及使用难度。

低势垒问题

前面已经指出，低势垒的途径需要特别注意。如果体系在演化过程中一直存在着势垒较其他途径低很多的一个或几个途径，会对模拟过程产生不利的影响。这个问题被称之为低势垒问题。低势垒途径对于KMC模拟最直接的影响就是大大缩短了模拟过程所涵盖的时间跨度。这一点可以从方程(5)中看出。更为深刻的影响在于，这些由低势垒的途径联系起来的组态会组成一个近似于封闭的族。体系会频繁的访问这些态，而其他的对于体系演化更为重要的高势垒途径被选择的概率非常低，这显然会降低KMC的模拟效率。例如，吸附原子在高指数金属表面扩散，其沿台阶的迁移所对应的势垒要远低于与台阶分离的移动。这样，KMC模拟的绝大部分时间内吸附原子都在台阶处来回往复，而不会选择离开台阶在平台上扩散。这显然不是我们希望看到的情形。一种解决办法是人为地将这些低势垒加高以降低体系访问这些组态的几率，但是无法预测这种干扰是否会造成体系对于真实情况的严重偏离。另一种选择是利用NRM或者CTSM进行模拟，但是其效果如何尚待检测。

如果考察体系的势能面，这类低势垒的途径一般处在一个"超势阱"之中。体系在这个超势阱中可以很快的达到热平衡，所需时间要短于从其中逸出的时间。如果可以明确的知道超势阱所包含的组态以及从超势阱逸出的所有途径，我们就可以按照Boltzmann分布合理的选择其中一条途径，使得体系向前演化。但是如何确定哪些组态包含在超势阱之中以及体系是否已在其中达到热平衡本身就是两个难题。对于第一点，Mason提出可以利用Zobrist密钥法标定访问过于频繁的组态 [7]；Novotny则提出通过建立及对角化一个描述体系在这些组态间演化的传递矩阵来解决第二点 [8]。对这个问题的详细讨论已超出了本文的讨论范围，请参阅文献[7]以及[8]。

实体动力学蒙特卡洛方法 OKMC

上述的KMC都假设任何时候原子均处于其理想点阵格子上。但是很多情况下这种点阵映射是无效的，比如间隙原子或者位错。这类结构缺陷的运动在材料的辐射损伤和老化过程中扮演着非常重要的角色。而且与单个原子或者空位的运动相比，这类缺陷的运动时间跨度更长，也更为复杂，比如间隙原子团和空穴的湮没，间隙原子团的解构/融合，或者位错的攀移/交滑移等等。传统的KMC算法很难有效的处理这类问题，一方面是因为时间跨度太大，另一方面这类缺陷各自均可视为独立的实体(object)，其运动更近似于系统激发，因此单个或几个原子运动的积累效果很多情况下并不能有效地反应这些实体的整体运动。实体动力学蒙特卡洛方法(Object kinetic Monte Carlo, OKMC)就是为了处理这类问题而被提出的。OKMC在算法上与普通的KMC完全一样。需要注意的地方是在OKMC中并不存在原子点阵。所有的实体在一个真空的箱子中按照其物理实质离散化运动，比如位错环的最小移动距离是其Burgers矢量大小，方向则为Burgers矢量方向；空位的移动距离为第一近邻或第二近邻的原子间距，等等。模拟过程中我们需要追踪该实体的形心，从而决定其位置、移动距离等等。此外，OKMC中对于跃迁速率$ k $的确定也和普通的KMC有所区别。本文前面已经指出，$ k_{ij} $可以表达为$ k_{ij}=\nu_0\exp\big(-\Delta E_{ij}/k_BT\big) $的形式。普通的KMC假定$ \nu_0 $为常数，不同途径的$ k_{ij} $由$ \Delta E_{ij} $决定。但是在OKMC的模拟中，$ \Delta E_{ij} $的直接确定非常困难，因此一般的策略是对于特定的事件(包括实体自身的运动以及不同实体间的反应等)，跃迁势垒$ \Delta E_{ij} $保持恒定，而将前置因子$ \nu_0 $视为实体规模(所包含的原子/空位数目)的函数，通过MD模拟得出，一般而言可以表示为形如$ \nu_0(m)=\nu_0(q^{-1})^{m-s} $的表达式，其中$ q $和$ s $是拟合参量，$ m $是实体规模。最后需要注意的是在OKMC的模型中，实体有空间范围，因此需要一个额外的参数$ R_{\rm entity} $来表征其空间半径(假设为球形分布，否则$ R_{\rm entity} $的数目多于一个)。在模拟不同实体间的反应时，需要特别考虑其形心的间距，如果小于"反应距离"，即$ R^1+R^2 $，反应一定进行，否则认为两个实体互相独立。

Domain利用OKMC研究了Fe-Cu合金的辐射损伤[9]，在模拟中考虑了间隙原子(空位)的聚合、间隙原子(空位)团的发射、间隙原子-空位湮没、空位团对杂质的捕获、表面对于空位(团)的捕获、甚至辐射轰击引起的间隙原子(空位)萌生、增殖等等事件。从中可以看出，对于OKMC，一个棘手的问题是需要预先想到所有的事件。此外，OKMC所需要的所有参量基本上不可能通过原子模拟直接获得，人为的设定参数不可避免。这些参数会在多大程度上决定OKMC的准确程度无法预先得知。需要根据现有的实验数据进行修改、调试。这些困难都限制了OKMC的普及。但是如前所述，这种方法可以有效地进行大尺度的时间(天)和空间模拟($ \mu $m以上)，而且对于缺陷的描述更为直接和符合直观，因此在材料研究中同样占有重要的地位。

KMC的若干进展

等时蛙跳算法($ \tau $-leap KMC)
 

引入这类算法前，我们先简要介绍两个常用的离散分布：泊松分布(Poisson Distribution, PD)以及二项式分布(Binomial Distribution, BD)。

泊松随机数$ \mathcal{P}(k_{ij},\tau) $定义为给定事件发生率$ k_{ij} $以及观测时间$ \tau $下事件发生的数目。如果用$ n $代表给定的发生数目，则$ \mathcal{P}(k_{ij},t) $恰好等于$ n $的概率是一个泊松分布：

$ P(n;k_{ij},\tau)=Pr\{\mathcal{P}(k_{ij},\tau)=n\}=\frac{e^{-k_{ij}\tau}(k_{ij}\tau)^k}{k!} $       (14)

也即如果产生一个泊松随机数序列$ \{\mathcal{P}(k_{ij},\tau)\} $，则这个序列符合泊松分布PD。需要指出，$ \mathcal{P}(k_{ij},\tau) $是无界的，范围是任意非负整数。

与其类似，二项式随机数$ \mathcal{B}(N,p) $定义为重复$ N $次独立的成功率均为$ p $的伯努利实验的成功数。如果给定成功数$ n $，则$ \mathcal{B}(N,p) $恰好等于$ n $的概率是一个二项式分布：

$ B(n;N,p)=Pr\{\mathcal{B}(N,p)=n\}=\frac{N!}{n!(N-n)!}p^n(1-p)^{N-n}, n=0,1,\ldots,N $       (15)

为了和本文中的标号一致，我们将跃迁$ i\rightarrow j $的成功率$ p $表示为$ k_{ij}\tau/N $，将方程(15)重新写为

$ B(n;N,k_{ij},\tau)=\frac{N!}{n!(N-n)!}\bigg(k_{ij}\tau/N\bigg)^n(1-k_{ij}\tau/N)^{N-n} $       (16)

与PD不同，BD中的$ n $是有界的，为0到$ N $之间的任意整数。

可以看出，如果将这两种随机数理解为给定跃迁路径(发生率为$ k_{ij} $)在一定的时间步长($ \tau $)内发生的次数，则可以立即运用于粒子数空间$ \mathbb{N} $内的KMC中，其时间范围可以得到很大提高。这就是等时蛙跳算法$ \tau $-leap KMC [10,11]。$ \tau $-leap KMC方法最早由Gillespie提出，通过PD[方程(14)]，在给定时间步长$ \tau $下决定每个跃迁途径$ j $发生的次数$ n_j $，然后将体系移到这些跃迁累计发生后产生的新态。因为每一步模拟体系不止发生一次跃迁，所以模拟的速度可以大大加快。我们以多种反应物在化学反应炉中的演化为例加以详细说明。

设在炉内共有$ L $种分子$ \{S_1,\ldots,S_L\} $，在时刻$ t $各自的个数为$ X_l(t) $，则在粒子数空间$ \mathbb{N} $中$ [X_1(t),\ldots,X_L(t)]^{\rm T} $构成一个矢量$ \mathbf{X}(t) $，或称为一个组态$ i $。总共有$ M $种反应路径$ R_j(j=1,\ldots,M) $。对于给定的$ R_j $，反应速率$ k_j[\mathbf{X}(t)] $是占据态$ \mathbf{X}(t) $的函数。此外，我们单独定义一个矢量
$ \bm\nu_j=\mathbf{X}^1-\mathbf{X}^0 $，其中$ \mathbf{X}^1 $由$ \mathbf{X}^0 $通过反应$ R_j $而得，即$ \mathbf{X}^0\xrightarrow{R_j}\mathbf{X}^1 $。因此$ \bm\nu_j $的元素$ \nu_{lj} $代表反应$ R_j $所引起的$ S_l $种分
子的数目变化。由此建立算法如下：

VI. PD-$ \tau $-leap KMC [10]

  1. 给定恒定时间步长$ \tau $；
  2. 对于每条反应途径$ R_j $按照方程(14)生成泊松随机数序列$ \bigg\{\mathcal{P}^j\big\{k_j[\mathbf{X}(t)],\tau\big\}\bigg\} $，按照模拟步数从序列中找出每种反应发生的次数$ n_j $；
  3. 按照$ \mathbf{X}(t+\tau)=\mathbf{X}(t)+\overset{M}{\underset{j=1}{\sum}}\bm\nu_j n_j $更新体系；
  4. 模拟时间前进$ \tau $；
  5. 重复上述过程。

Gillespie仔细考虑了$ \tau $的选择条件，称为蛙跳条件(leap condition)：

$ \tau=\underset{j\in[1,M]}{\min}\bigg\{\frac{\epsilon k_{\rm tot}(\mathbf{X})}{\mu_j(\mathbf{X})},\frac{\epsilon^2k^2_{\rm tot}(\mathbf{X})}{\sigma^2_j(\mathbf{X})}\bigg\} $       (17)

其中

$ f_{jm}(\mathbf{X})=\overset{L}{\underset{l=1}{\sum}}\frac{\partial k_j(\mathbf{X})}{\partial X_l}\nu_{lk}, j,m=1,\ldots,M $
$ \mu_j(\mathbf{X})=\overset{M}{\underset{m=1}{\sum}}f_{jm}(\mathbf{X})k_m(\mathbf{X}), j=1,\ldots,M $
$ \sigma^2_j(\mathbf{X})=\overset{M}{\underset{m=1}{\sum}}f^2_{jm}(\mathbf{X})k_m(\mathbf{X}), j=1,\ldots,M $
 

如前所述，$ \mathcal{P}(k_{ij},\tau) $没有上限，因此即使$ \tau $满足方程(17)，在模拟过程中也可能会出现某种分子总数为负数的情况，这显然不符合实际，也是PD-$ \tau $-leap KMC的一个弱点。Tian和Burrage提出可以用二项式分布BD取代PD，因为$ \mathcal{B}(N,p) $有上限，所以可以有效的解决这个问题。此外，他们对于某种分子参与多种反应的情况也进行了考虑，从而提高了$ \tau $-leap KMC的稳定性和普适性。其算法如下：

VII. BD-$ \tau $-leap KMC [11]
 

  1. 给定恒定时间步长$ \tau $，满足$ 0\leq\frac{k_j(\mathbf{X})\tau}{N_j}\leq1 $；
  2. 对于每条反应途径$ R_j $按照方程(16)生成二项式随机数序列$ \bigg\{\mathcal{B}^j\big\{N_j;k_j[\mathbf{X}(t)],\tau\big\}\bigg\} $，按照模拟步数从序列中找出每种反应发生的次数$ n_j $；如果有某种分子$ S_l $同时参与了$ R_j $和$ R_m $，则首先生成
    $$\mathcal{B}\bigg[n_{jm};\frac{k_j(\mathbf{X})}{k_j(\mathbf{X})+k_m(\mathbf{X})}\bigg],$$

    然后通过$ n_m = n_{jm}-n_j $确定$ R_m $的发生次数；

  3. 按照$ \mathbf{X}(t+\tau)=\mathbf{X}(t)+\overset{M}{\underset{j=1}{\sum}}\bm\nu_j n_j $更新体系；
  4. 模拟时间前进$ \tau $；
  5. 重复上述过程。

    6. 步骤1、2中出现的$ N_j $是参与反应$ R_j $的各类分子的个数的最小值，即

    $$N_j=\underset{m\in R_j}{\min}\big\{\{X_m\}\big\}.$$

    此外Gillespie，Tian和Burrage还考虑用预测$ \tau/2 $时刻体系状态的方法来进一步提高精度。具体请参阅文献[10,11]。如果$ \tau $-leap算法和OKMC结合起来可以进一步加大模拟的时间尺度，但是目前还没有这方面工作的介绍。

    基于即时动态分析的KMC方法(on-the-fly KMC)

    到目前为止，所有的KMC都是在模拟之前建立好所有可能的跃迁途径。但是实际上"所有"是很难达到的目标。因为很多途径远离一般的直觉，而且在演化过程中体系有可能寻找到新的途径。因此，跃迁途径应该随着体系的演化而不断更新，是动态的过程。Henkelman和Jönsson将途径搜索和KMC结合起来，提出了即时动态的KMC方法 on-the-fly KMC [12]：在每一个稳态(势阱)处，选定一个激活原子(一般是近邻不饱和的原子)，在以其为中心的局部区域内引入呈高斯分布的随机位移，即加入扰动，然后利用dimer方法 [13]寻找所有可能的跃迁途径。建立起即时的途径库之后再通过普通KMC算法进行模拟。显然，这种方法的计算量非常大，需要一个有效的标识方法来识别所有已经遇到过的途径以避免重复计算。Trushin提出可以利用包括至激活原子第三壳层的所有格点(顺时针排列)的占据与否(分别标记为1和0)来构建二进制数，从而根据始态和终态的标号来唯一地标识某条途径 [14]，例如，激活原子标为"1"，其第一壳层的原子标记为"2","3",$ \ldots $,"$ N_{1} $"，依此类推，然后将原子的标号"$ i $"作为二进制的数位$ 2^i $，这样，每一个稳态都有唯一的一个二进制数与之对应。虽然仍不完善，但是这种方法具有非常清晰的逻辑结构，具有良好的扩展性。

    $ \mathcal{O}(\log_2N) $和$ \mathcal{O}(1) $KMC方法

    一般情况下KMC的大部分时间花费在选择途径上。如果采用普通的方法，即循环叠加$ k_{ij} $直至$ r_1k_{\rm tot}&lt;\overset{j_k}{\underset{j=1}{\sum}}k_{ij} $从而选择$ j_k $，这种情况下计算用时与途径数目$ M $呈线性增长，即$ \mathcal{O}(N) $算法。按照二叉树安排不同数目的$ k_{ij} $之和可以改进到$ \mathcal{O}(\log_2N) $ [15]：
    将所有$ k_{ij} $作为树叶(不足2整数次幂的叶子由0填补)，每两片叶子之和作为父节点，依次类推直至树根$ k_{\rm tot} $。一株二叉树构建完毕后，生成一个随机数$ r_1 $，由树根开始寻找$ s=r_1k_{\rm tot} $，若$ s $不大于左子节点$ k_{\rm left} $，沿左分支向下寻找；否则设$ s=s-k_{\rm left} $，沿右分支向下寻找，直至树叶$ j $，体系按途径$ j $演化。

    Slepoy和Thompson等进一步提出分流-拒绝(composition-rejection, CR)方法以实现搜索用时与途径总数无关的$ \mathcal{O}(1) $算法 [16]：(1)先找出$ k_{\rm min} $和$ k_{\rm max} $，按照($ k_{\rm min},2k_{\rm min},4k_{\rm min},\ldots,k_{\rm max} $)将$ M $条途径分为$ G $个组，$ G=\log_2(\frac{k_{\rm max}}{k_{\rm min}}) $，(2)然后生成随机数$ r_1 $，按照上述二叉树寻找$ s=r_1k_{\rm tot} $所落入的组别$ G_j $，(3)再生成两个随机数$ r_2 $和$ r_3 $，设$ l={\rm Int}(r_2/N_{G_j}) $，其中$ N_{G_j} $为该组中包含的途径数，$ t=r_3k_{\rm max} $，如果$ t\leq k_l $，则选择途径$ l $，否则重复步骤(3)，直至有一条途径被选中为止。可以看出，CR算法虽然搜索速度很快，但是每一步KMC需要产生至少4个随机数($ r_4 $用于确定前进时间)，因此需要高质量的随机数发生器。不过对于跃迁途径复杂的体系演化而言，CR的$ \mathcal{O}(1) $效率无疑是很有吸引力的。

     [1] A.F. Voter, {\it Radiation Effects in Solids} (Springer 2006) p. 1-24.
     [2] H. Eyring, J. Chem. Phys. 3, 107 (1935).
     [3] P. Kratzer, Multiscale Simulation Method in Molecular Science (NIC Serices, Vol. 42, Forschungszentrum, Jülich 2009) p. 51-76.
     [4] E.J. Dawnkaski, D. Srivastava and B.J. Gamson, J. Chem. Phys. 102, 9401 (1995).
     [5] A.F. Voter and J.D. Doll, J. Chem. Phys. 80, 5832 (1984).
     [6] A.F. Voter, Phys. Rev. B 34, 6819 (1986).
     [7] D.R. Mason, T.S. Hudson and A.P. Sutton, Comp. Phys. Comm. 165, 37 (2005).
     [8] M.A. Novotny, Phys. Rev. Lett. 74, 1 (1994); Erratum 75, 1424 (1995).
     [9] C. Domain, C.S. Becquart and L. Malerba, J. Nucl. Mater. 335, 121 (2004).
     [10] D.T. Gillespie, J. Chem. Phys. 115, 1716 (2001).
     [11] T. Tian and K. Burrage, J. Chem. Phys. 121, 10356 (2004).
     [12] G. Henkelman and H. J\'{o}nsson, J. Chem. Phys. 115, 9657 (2001).
     [13] G. Henkelman and H. J\'{o}nsson, J. Chem. Phys. 111, 7010 (1999).
     [14] O. Trushin, A. Karim, A. Kara and T.S. Rahman, Phys. Rev. B 72, 115401 (2005).
     [15] M.A. Gibson and J. Bruck, J. Phys. Chem. A 104, 1876 (2000).
     [16] A. Slepoy, A.P. Thompson and S.J. Plimpton, J. Chem. Phys. 128, 205101 (2008).