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

  • strict warning: Non-static method view::load() should not be called statically in /home/vasp/wwwroot/drupal-6.33/sites/all/modules/views/views.module on line 906.
  • strict warning: Declaration of views_handler_argument::init() should be compatible with views_handler::init(&$view, $options) in /home/vasp/wwwroot/drupal-6.33/sites/all/modules/views/handlers/views_handler_argument.inc on line 744.
  • strict warning: Declaration of views_plugin_row::options_validate() should be compatible with views_plugin::options_validate(&$form, &$form_state) in /home/vasp/wwwroot/drupal-6.33/sites/all/modules/views/plugins/views_plugin_row.inc on line 134.
  • strict warning: Declaration of views_plugin_row::options_submit() should be compatible with views_plugin::options_submit(&$form, &$form_state) in /home/vasp/wwwroot/drupal-6.33/sites/all/modules/views/plugins/views_plugin_row.inc on line 134.

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

KMC方法基本原理

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

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

指数分布与KMC的时间步长

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

因为体系在势能面上无记忆的随机行走,所以任意单位时间内,它找到跃迁途径的概率不变,设为TeX Embedding failed!。因此在TeX Embedding failed!区间内,体系不发生跃迁的概率为

TeX Embedding failed!

类似的,在TeX Embedding failed!区间内,体系不发生跃迁的概率为

TeX Embedding failed!

以此类推,当TeX Embedding failed!时,在TeX Embedding failed!区间内,体系不发生跃迁的概率为

TeX Embedding failed!

因此,当TeX Embedding failed!趋于TeX Embedding failed!时,体系不发生跃迁的概率为

TeX Embedding failed!      (1)

这一行为类似于原子核的衰变方程。从方程(1)我们可以得到单位时间内体系跃迁概率TeX Embedding failed!。从方程(1)的推导过程可以看出体系的跃迁概率是一个随时间积累的物理量,因此TeX Embedding failed!对时间积分到某一时刻TeX Embedding failed!必然等于TeX Embedding failed!,也即TeX Embedding failed!。因此我们立即可以得到 [1]

TeX Embedding failed!       (2)

TeX Embedding failed!是体系处于态TeX Embedding failed!时所有可能的跃迁途径的速率TeX Embedding failed!之和,即

TeX Embedding failed!       (3)

对于每个具体的跃迁途径TeX Embedding failed!,上述讨论均成立。因此,我们可以定义单位时间内体系进行TeX Embedding failed!跃迁的概率TeX Embedding failed!

TeX Embedding failed!       (4)
 

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

TeX Embedding failed!

从而

TeX Embedding failed!       (5)

最后一步是因为TeX Embedding failed!TeX Embedding failed!的分布相同。TeX Embedding failed!也可以通过上述步骤从方程(4)得到。

计算跃迁速率

过渡态理论(TST)

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

TeX Embedding failed!       (6)

方程(6)中TeX Embedding failed!表示在组态A所属态空间里对正则系综的平均。TeX Embedding failed!表示只考虑体系从组态A迁出而不考虑迁入A的情况(后一种情况体系也对通过纵截面的流量有贡献)。根据普遍公式

TeX Embedding failed!

设体系的哈密顿量TeX Embedding failed!TeX Embedding failed!,即可分解为动能和势能,同时设粒子坐标TeX Embedding failed!时体系处于组态A。则方程(6)可写为

TeX Embedding failed!
     TeX Embedding failed!
     TeX Embedding failed!
     TeX Embedding failed!       (7)

上式中无限小量TeX Embedding failed!是为了将TeX Embedding failed!函数全部包含进去。最后一项对于TeX Embedding failed!函数的系综平均可以直接通过Metropolis Monte Carlo方法计算出来:计算粒子落在TeX Embedding failed!范围内的次数相对于Metropolis行走总次数的比例TeX Embedding failed!。方程(7)最后等于

TeX Embedding failed!        (8)

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

TeX Embedding failed!       (9)

其中TeX Embedding failed!是纵截面方程,TeX Embedding failed!代表3维情况中粒子流动方向与截面TeX Embedding failed!法向不平行对于计数的影响。

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

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

TeX Embedding failed!       (10)

其中TeX Embedding failed!为在跃迁TeX Embedding failed!中体系在鞍点和态TeX Embedding failed!处的自由能之差

TeX Embedding failed!

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

TeX Embedding failed!       (11)

hTST认为体系在稳态附近的振动可以用谐振子表示,因此其配分函数是经典谐振子体系的配分函数。分别写出体系在态TeX Embedding failed!和鞍点处的配分函数TeX Embedding failed!TeX Embedding failed!

TeX Embedding failed!
TeX Embedding failed!

根据Boltzmann公式,

TeX Embedding failed!       (12)

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

TeX Embedding failed!       (13)

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

KMC几种不同的实现算法
 

点阵映射
 

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

无拒绝方式
 

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

I. 直接法

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

  1. 根据方程(4)计算体系处于态TeX Embedding failed!时的总跃迁速率TeX Embedding failed!
  2. 选择随机数TeX Embedding failed!
  3. 寻找途径TeX Embedding failed!,满足TeX Embedding failed!
  4. 体系移动到态TeX Embedding failed!,同时模拟时间前进TeX Embedding failed!
  5. 重复上述过程。

需要指出的是,虽然一般步骤4中的TeX Embedding failed!根据方程(5)生成,但是如果将其换为TeX Embedding failed!并不会影响模拟结果。在文献[5]和[6]中均采用这种方式。

II. 第一反应法

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

  1. 设共有TeX Embedding failed!条反应途径,生成TeX Embedding failed!个随机数TeX Embedding failed!
  2. 根据公式TeX Embedding failed!,给出每条路径的预计发生时间;
  3. 找出TeX Embedding failed!的最小值TeX Embedding failed!
  4. 体系移动到态TeX Embedding failed!,同时模拟时间前进TeX Embedding failed!
  5. 重复上述过程。

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

III. 次级反应法

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

试探-接受/拒绝方式
 

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

IV. 选择直接法

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

  1. 设共有TeX Embedding failed!条反应途径,选择反应速率最大值TeX Embedding failed!,设为TeX Embedding failed!。生成在TeX Embedding failed!均匀分布的随机数TeX Embedding failed!
  2. TeX Embedding failed!
  3. 如果TeX Embedding failed! < TeX Embedding failed!,则体系跃迁至新态TeX Embedding failed!,否则保持在态TeX Embedding failed!
  4. 模拟时间前进TeX Embedding failed!
  5. 重复上述过程。

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

V. 恒定步长法

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

  1. 给定恒定时间步长TeX Embedding failed!
  2. 将所有途径TeX Embedding failed!(共有TeX Embedding failed!个)设为长度恒为TeX Embedding failed!的线段,生成在TeX Embedding failed!均匀分布的随机数TeX Embedding failed!,选择途径TeX Embedding failed!
  3. 生成在TeX Embedding failed!均匀分布的随机数TeX Embedding failed!,如果TeX Embedding failed! < TeX Embedding failed!,则体系跃迁至新态TeX Embedding failed!,否则保持在态TeX Embedding failed!
  4. 模拟时间前进TeX Embedding failed!
  5. 重复上述过程。

实际模拟中,TeX Embedding failed!需要满足(1)小于TeX Embedding failed!(见"第一反应法"),以及(2)对于TeX Embedding failed!最大的途径,接受率大致在0.5。其中第一个条件保证了所有的迁移途径发生概率都小于1,第二个条件则保证体系演化的效率不会过于低下。CTSM是非常行之有效的一类KMC算法,但是选择TeX Embedding failed!时需要特别的注意以保证效率。TeX Embedding failed!决定于具体体系以及模拟温度。这在一定程度上增加了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中对于跃迁速率TeX Embedding failed!的确定也和普通的KMC有所区别。本文前面已经指出,TeX Embedding failed!可以表达为TeX Embedding failed!的形式。普通的KMC假定TeX Embedding failed!为常数,不同途径的TeX Embedding failed!TeX Embedding failed!决定。但是在OKMC的模拟中,TeX Embedding failed!的直接确定非常困难,因此一般的策略是对于特定的事件(包括实体自身的运动以及不同实体间的反应等),跃迁势垒TeX Embedding failed!保持恒定,而将前置因子TeX Embedding failed!视为实体规模(所包含的原子/空位数目)的函数,通过MD模拟得出,一般而言可以表示为形如TeX Embedding failed!的表达式,其中TeX Embedding failed!TeX Embedding failed!是拟合参量,TeX Embedding failed!是实体规模。最后需要注意的是在OKMC的模型中,实体有空间范围,因此需要一个额外的参数TeX Embedding failed!来表征其空间半径(假设为球形分布,否则TeX Embedding failed!的数目多于一个)。在模拟不同实体间的反应时,需要特别考虑其形心的间距,如果小于"反应距离",即TeX Embedding failed!,反应一定进行,否则认为两个实体互相独立。

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

KMC的若干进展

等时蛙跳算法(TeX Embedding failed!-leap KMC)
 

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

泊松随机数TeX Embedding failed!定义为给定事件发生率TeX Embedding failed!以及观测时间TeX Embedding failed!下事件发生的数目。如果用TeX Embedding failed!代表给定的发生数目,则TeX Embedding failed!恰好等于TeX Embedding failed!的概率是一个泊松分布:

TeX Embedding failed!       (14)

也即如果产生一个泊松随机数序列TeX Embedding failed!,则这个序列符合泊松分布PD。需要指出,TeX Embedding failed!是无界的,范围是任意非负整数。

与其类似,二项式随机数TeX Embedding failed!定义为重复TeX Embedding failed!次独立的成功率均为TeX Embedding failed!的伯努利实验的成功数。如果给定成功数TeX Embedding failed!,则TeX Embedding failed!恰好等于TeX Embedding failed!的概率是一个二项式分布:

TeX Embedding failed!       (15)

为了和本文中的标号一致,我们将跃迁TeX Embedding failed!的成功率TeX Embedding failed!表示为TeX Embedding failed!,将方程(15)重新写为

TeX Embedding failed!       (16)

与PD不同,BD中的TeX Embedding failed!是有界的,为0到TeX Embedding failed!之间的任意整数。

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

设在炉内共有TeX Embedding failed!种分子TeX Embedding failed!,在时刻TeX Embedding failed!各自的个数为TeX Embedding failed!,则在粒子数空间TeX Embedding failed!TeX Embedding failed!构成一个矢量TeX Embedding failed!,或称为一个组态TeX Embedding failed!。总共有TeX Embedding failed!种反应路径TeX Embedding failed!。对于给定的TeX Embedding failed!,反应速率TeX Embedding failed!是占据态TeX Embedding failed!的函数。此外,我们单独定义一个矢量
TeX Embedding failed!,其中TeX Embedding failed!TeX Embedding failed!通过反应TeX Embedding failed!而得,即TeX Embedding failed!。因此TeX Embedding failed!的元素TeX Embedding failed!代表反应TeX Embedding failed!所引起的TeX Embedding failed!种分
子的数目变化。由此建立算法如下:

VI. PD-TeX Embedding failed!-leap KMC [10]

  1. 给定恒定时间步长TeX Embedding failed!
  2. 对于每条反应途径TeX Embedding failed!按照方程(14)生成泊松随机数序列TeX Embedding failed!,按照模拟步数从序列中找出每种反应发生的次数TeX Embedding failed!
  3. 按照TeX Embedding failed!更新体系;
  4. 模拟时间前进TeX Embedding failed!
  5. 重复上述过程。

Gillespie仔细考虑了TeX Embedding failed!的选择条件,称为蛙跳条件(leap condition):

TeX Embedding failed!       (17)

其中

TeX Embedding failed!
TeX Embedding failed!
TeX Embedding failed!
 

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

VII. BD-TeX Embedding failed!-leap KMC [11]
 

  1. 给定恒定时间步长TeX Embedding failed!,满足TeX Embedding failed!
  2. 对于每条反应途径TeX Embedding failed!按照方程(16)生成二项式随机数序列TeX Embedding failed!,按照模拟步数从序列中找出每种反应发生的次数TeX Embedding failed!;如果有某种分子TeX Embedding failed!同时参与了TeX Embedding failed!TeX Embedding failed!,则首先生成
    TeX Embedding failed!

    然后通过TeX Embedding failed!确定TeX Embedding failed!的发生次数;

  3. 按照TeX Embedding failed!更新体系;
  4. 模拟时间前进TeX Embedding failed!
  5. 重复上述过程。
  6. 步骤1、2中出现的TeX Embedding failed!是参与反应TeX Embedding failed!的各类分子的个数的最小值,即

    TeX Embedding failed!

    此外Gillespie,Tian和Burrage还考虑用预测TeX Embedding failed!时刻体系状态的方法来进一步提高精度。具体请参阅文献[10,11]。如果TeX Embedding failed!-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",TeX Embedding failed!,"TeX Embedding failed!",依此类推,然后将原子的标号"TeX Embedding failed!"作为二进制的数位TeX Embedding failed!,这样,每一个稳态都有唯一的一个二进制数与之对应。虽然仍不完善,但是这种方法具有非常清晰的逻辑结构,具有良好的扩展性。

    TeX Embedding failed!TeX Embedding failed!KMC方法

    一般情况下KMC的大部分时间花费在选择途径上。如果采用普通的方法,即循环叠加TeX Embedding failed!直至TeX Embedding failed!从而选择TeX Embedding failed!,这种情况下计算用时与途径数目TeX Embedding failed!呈线性增长,即TeX Embedding failed!算法。按照二叉树安排不同数目的TeX Embedding failed!之和可以改进到TeX Embedding failed! [15]:
    将所有TeX Embedding failed!作为树叶(不足2整数次幂的叶子由0填补),每两片叶子之和作为父节点,依次类推直至树根TeX Embedding failed!。一株二叉树构建完毕后,生成一个随机数TeX Embedding failed!,由树根开始寻找TeX Embedding failed!,若TeX Embedding failed!不大于左子节点TeX Embedding failed!,沿左分支向下寻找;否则设TeX Embedding failed!,沿右分支向下寻找,直至树叶TeX Embedding failed!,体系按途径TeX Embedding failed!演化。

    Slepoy和Thompson等进一步提出分流-拒绝(composition-rejection, CR)方法以实现搜索用时与途径总数无关的TeX Embedding failed!算法 [16]:(1)先找出TeX Embedding failed!TeX Embedding failed!,按照(TeX Embedding failed!)将TeX Embedding failed!条途径分为TeX Embedding failed!个组,TeX Embedding failed!,(2)然后生成随机数TeX Embedding failed!,按照上述二叉树寻找TeX Embedding failed!所落入的组别TeX Embedding failed!,(3)再生成两个随机数TeX Embedding failed!TeX Embedding failed!,设TeX Embedding failed!,其中TeX Embedding failed!为该组中包含的途径数,TeX Embedding failed!,如果TeX Embedding failed!,则选择途径TeX Embedding failed!,否则重复步骤(3),直至有一条途径被选中为止。可以看出,CR算法虽然搜索速度很快,但是每一步KMC需要产生至少4个随机数(TeX Embedding failed!用于确定前进时间),因此需要高质量的随机数发生器。不过对于跃迁途径复杂的体系演化而言,CR的TeX Embedding failed!效率无疑是很有吸引力的。

     [1] A.F. Voter, 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).