《KM知识管理动力学蒙特卡洛方法KMC及相关讨论》由会员分享，可在线阅读，更多相关《KM知识管理动力学蒙特卡洛方法KMC及相关讨论（14页珍藏版）》请在金锄头文库上搜索。
1、（KM知识管理)动力学蒙特卡洛方法(KMC)及相关讨论动力学蒙特卡洛方法(KMC)及相关讨论星期二,2010-05-1101:05satchel1979动态模拟在目前的计算科学中占据着非常重要的位置。随着计算能力和第一原理算法的发展，复杂的动态参数（扩散势垒、缺陷相互作用能等）均可利用第一原理计算得出。因此，部分复杂的体系动态变化，如表面形貌演化或辐射损伤中缺陷集团的聚合-分解演变等，已可以较为精确的予以研究。KMC动力学蒙特卡洛方法（kineticMonteCarlo）原理简单，适应性强，因此在很多情况下都是研究人员的首选。此外，KMC在复杂体系或复杂过程中的算法发展也非常活跃。本文试图介绍
2、KMC方法的基础理论和若干进展。KMC方法基本原理在原子模拟领域内，分子动力学(moleculardynamics,MD)具有突出的优势。它可以非常精确的描述体系演化的轨迹。一般情况下MD的时间步长在飞秒(s)量级，因此足以追踪原子振动的具体变化。但是这一优势同时限制了MD在大时间尺度模拟上的应用。现有的计算条件足以支持MD到10ns，运用特殊的算法可以达到10s的尺度。即便如此，很多动态过程，如表面生长或材料老化等，时间跨度均在s以上，大大超出了MD的应用范围。有什么方法可以克服这种局限呢？当体系处于稳定状态时，我们可以将其描述为处于维势能函数面的一个局域极小值（阱底）处。有限温度下，虽然体
3、系内的原子不停的进行热运动，但是绝大部分时间内原子都是在势能阱底附近振动。偶然情况下体系会越过不同势阱间的势垒从而完成一次“演化”，这类小概率事件才是决定体系演化的重点。因此，如果我们将关注点从“原子”升格到“体系”，同时将“原子运动轨迹”粗化为“体系组态跃迁”，那么模拟的时间跨度就将从原子振动的尺度提高到组态跃迁的尺度。这是因为这种处理方法摈弃了与体系穿越势垒无关的微小振动，而只着眼于体系的组态变化。因此，虽然不能描绘原子的运动轨迹，但是作为体系演化，其“组态轨迹”仍然是正确的。此外，因为组态变化的时间间隔很长，体系完成的连续两次演化是独立的，无记忆的，所以这个过程是一种典型的马尔可夫过程(
4、Markovprocess)，即体系从组态到组态，这一过程只与其跃迁速率有关。如果精确地知道，我们便可以构造一个随机过程，使得体系按照正确的轨迹演化。这里正确的意思是某条给定演化轨迹出现的几率与MD模拟结果完全一致(假设我们进行了大量的MD模拟，每次模拟中每个原子的初始动量随机给定)。这种通过构造随机过程研究体系演化的方法即为动力学蒙特卡洛方法(kineticMonteCarlo,KMC)1。指数分布与KMC的时间步长在KMC模拟中，构造呈指数分布的随机数是一个相当重要的步骤。这一节中我们对此进行讨论。因为体系在势能面上无记忆的随机行走，所以任意单位时间内，它找到跃迁途径的概率不变，设为。因此
5、在区间内，体系不发生跃迁的概率为类似的，在区间内，体系不发生跃迁的概率为以此类推，当时，在区间内，体系不发生跃迁的概率为因此，当趋于时，体系不发生跃迁的概率为(1)这一行为类似于原子核的衰变方程。从方程(1)我们可以得到单位时间内体系跃迁概率。从方程(1)的推导过程可以看出体系的跃迁概率是一个随时间积累的物理量，因此对时间积分到某一时刻必然等于，也即。因此我们立即可以得到1(2)是体系处于态时所有可能的跃迁途径的速率之和，即(3)对于每个具体的跃迁途径，上述讨论均成立。因此，我们可以定义单位时间内体系进行跃迁的概率为(4)单位时间内体系的跃迁概率呈指数分布这一事实说明KMC的时间步长也应是指数
6、分布。因此我们需要产生一个指数分布的随机数序列。这一点可以非常容易的通过一个(0,1平均分布的随机数序列转化得到：从而(5)最后一步是因为和的分布相同。也可以通过上述步骤从方程(4)得到。计算跃迁速率过渡态理论(TST)决定了KMC模拟的精度甚至准确性。为避开通过原子轨迹来确定的做法(这样又回到了MD的情况)，一般情况下采用过渡态理论(transitionstatetheory,TST)进行计算2。在TST中，体系的跃迁速率决定于体系在鞍点处的行为，而平衡态(势阱)处的状态对其影响可以忽略不计。如果大量的相同的体系组成正则系综，则在平衡状态下体系在单位时间内越过某个垂直于跃迁途径的纵截面的流量
7、即为。简单起见，假设有大量相同的一维双组态(势阱)体系，平衡状态下鞍点所在的假想面(对应于流量最小的纵截面)为，则TST给出该体系从组态A迁出到B的速率为5,6(6)方程(6)中表示在组态A所属态空间里对正则系综的平均。表示只考虑体系从组态A迁出而不考虑迁入A的情况(后一种情况体系也对通过纵截面的流量有贡献)。根据普遍公式设体系的哈密顿量为，即可分解为动能和势能，同时设粒子坐标时体系处于组态A。则方程(6)可写为(7)上式中无限小量是为了将函数全部包含进去。最后一项对于函数的系综平均可以直接通过MetropolisMonteCarlo方法计算出来：计算粒子落在范围内的次数相对于Metropol
8、is行走总次数的比例。方程(7)最后等于(8)将上述讨论扩展到3维情况非常直接，这里只给出结果，详细讨论请参阅文献5:(9)其中是纵截面方程，代表3维情况中粒子流动方向与截面法向不平行对于计数的影响。简谐近似下的过渡态理论(hTST)虽然上一节已经给出了TST计算跃迁速率的方法，但是在具体工作中，更多地是利用简谐近似下的过渡态理论(harmonicTST,hTST)通过解析表达式给出。根据TST，跃迁速率为3(10)其中为在跃迁中体系在鞍点和态处的自由能之差将上式代入方程(10)，可以得到(11)hTST认为体系在稳态附近的振动可以用谐振子表示，因此其配分函数是经典谐振子体系的配分函数。分别写
9、出体系在态和鞍点处的配分函数和：根据Boltzmann公式，(12)并将配分函数代入，则方程(11)得(13)方程(13)在通常的文献上经常可以见到。声子谱可以通过Hessian矩阵对角化或者密度泛函微扰法(DFPT)求出，而就是的势垒，可以通过NEB或者drag方法求出。因此，方程(13)保证了可以通过原子模拟(MD或者DFT方法)解析地求出。事实上这个方程有两点需要注意。首先虽然方程(10)中出现了普朗克常数，但是在最终结果中被抵消了。这是因为TST本质上是一个经典理论，所以充分考虑了统计效应后不会出现1。其次，方程(13)表明对于每一个跃迁过程，鞍点处的声子谱应该单独计算。这样会大大增加
10、计算量，因此在绝大部分计算中均设前置因子为常数，不随跃迁过程而变化。具体数值取决于体系，对于金属而言，一般取Hz。KMC几种不同的实现算法点阵映射到目前为止，进行KMC模拟的所有理论基础均已具备。但是前面所进行的讨论并没有联系到具体的模型。KMC在固体物理中的应用往往利用点阵映射将原子与格点联系起来。从而将跃迁(事件)具象化为原子格点关系的变化。比如空位(团)/吸附原子(岛)迁移等等。虽然与实际情况并不完全一致，但这样做在很多情况下可以简化建模的工作量，而且是非常合理的近似。很多情况下体系中的原子虽然对理想格点均有一定的偏离，但是并不太大()，因此这种原子点阵映射是有效的。这种做法的另一个好处
11、是可以对跃迁进行局域化处理。每条跃迁途径只与其近邻的体系环境有关，这样可以极大的减少跃迁途径的数目，从而简化计算1。需要指出的是，这种映射对于KMC模拟并不是必须的。比如化学分子反应炉或者生物分子的生长等等，这些情况下根本不存在点阵。无拒绝方式KMC的实现方法有很多种，这些算法大致可以分为拒绝(rejection)和无拒绝(rejection-free)两种范畴。每种范畴之下还有不同的实现方式。本文只选择几种最为常用的方法加以介绍。I.直接法直接法(directmethod)是最常用的一种KMC算法，其效率非常高。每一步只需要产生两个在之间平均分布的随机数和。其中被用来选定跃迁途径，确定模拟的
12、前进时间。设体系处于态，将每条跃迁途径想象成长度与跃迁速率成正比的线段。将这些线段首尾相连。如果落在线段中，这个线段所代表的跃迁途径就被选中，体系移动到态，同时体系时间根据方程(5)前进。总结其算法如下：1. 根据方程(4)计算体系处于态时的总跃迁速率；2. 选择随机数；3. 寻找途径，满足；4. 体系移动到态，同时模拟时间前进；5. 重复上述过程。需要指出的是，虽然一般步骤4中的根据方程(5)生成，但是如果将其换为并不会影响模拟结果。在文献5和6中均采用这种方式。II.第一反应法第一反应法(firstreactionmethod,FRM)在思路上比直接法更为自然。前面说过，对于处于稳态的体系
13、而言，它可以有不同的跃迁途径可以选择。每条途径均可以根据方程(4)给出一个指数分布的发生时间，也即从当前算起第一次发生的时间。然后从中选出最小值(最先发生的第一反应)，体系跃迁到相应的组态，模拟时间相应地前进。总结其算法如下：1. 设共有条反应途径，生成个随机数；2. 根据公式，给出每条路径的预计发生时间；3. 找出的最小值；4. 体系移动到态，同时模拟时间前进；5. 重复上述过程。可以看出，这种算法的效率比直接法低下，因为每一步KMC模拟需要生成个随机数。通常情况下KMC模拟需要步来达到较好的统计性质，如果每一步都需要生成个随机数，则利用这种方法需要一个高质量的伪随机数发生器，这一点在比较大
14、时尤为重要。III.次级反应法次级反应法(nextreactionmethod,NRM)是FRM方法的一种衍生方法，其核心思想是假设体系的一次跃迁并不会导致处于新态的体系对于其他跃迁途径的舍弃(比如充满可以发生种化学反应的分子，第一种反应发生并不会造成别的反应物的变化)，这样体系还可以选择中的次小值，从而跃迁到态，模拟时间前进。如果这次跃迁还可以满足上述假设条件，再重复上述过程。理想情况下，平均每一步KMC模拟只需要生成1个随机数。这无疑会大大提高效率以及时间跨度。但是实际上NRM的假设条件很难在体系每次跃迁之后都得到满足，在固体物理的模拟中尤其如此，因此其应用范围集中于研究复杂化学环境下的反
15、应过程。试探-接受/拒绝方式这一大类算法虽然在效率上不如直接法，但是它们所采用的试探-接受/拒绝在形式上更接近MetropolisMC方法，而且可以很方便的引入恒定步长，即固定。因此有必要进行详细的介绍。IV.选择直接法选择直接法在决定体系是否跃迁方面和MetropolisMC方法形式上非常相像，均是通过产生随机数和预定的阈值比较决定事件是否被采纳。具体算法如下：1. 设共有条反应途径，选择反应速率最大值，设为。生成在均匀分布的随机数；2. 设；3. 如果，则体系跃迁至新态，否则保持在态；4. 模拟时间前进；5. 重复上述过程。这种方法的长处在于每一步只需要生成一个随机数。但是缺点也很明显，对于反应速率相差太大，尤其是只有一个低势垒途径(与其他途径相比过大)的体系来讲，这种方法的效率会非常低下。某些情况下，这种低效率问题可以通过如下方法改进：将全部途径按照的大小分为几个亚组，每个亚组选定一个上限。但是这一步骤在整个KMC模拟过程中可能需要重复很多次，因此并不能完全解决问题。事实上低势垒在KMC中是个普遍的问题。这一点在后面还要简要提及。V.恒定步长法与上述四种方法不同，恒定步长法(constanttimestepmethod,CTSM)中体系的前进时间是个给定的参数citedawnkaski。在理想情况下，CTSM与直接法效率相同，每一步只需产生两个随机数。具体算法如
