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

动力学蒙特卡洛方法(KMC)及相关讨论星期二，2010-05-1101:05satchel1979动态模拟在目前的计算科学中占据着非常重要的位置。随着计算能力和第一原理算法的开展，复杂的动态参数(扩散势垒、缺陷相互作用能等)均可利用第一原理计算得出。因此，局部复杂的体系动态变化，如外表形貌演化或辐射损伤中缺陷集团的聚合-分解演变等，已可以较为精确的予以研究。KMC动力学蒙特卡洛方法(kineticMonteCarlo)原理简单，适应性强，因此在很多情况下都是研究人员的首选。此外，KMC在复杂体系或复杂过程中的算法开展也非常活泼。本文试图介绍KMC方法的基础理论和若干进展。KMC方法根本原理在原子模拟领域内，分子动力学(moleculardynamics,MD)具有突出的优势。它可以非常精确的描述体系演化的轨迹。一般情况下MD的时间步长在飞秒(10-%)量级，因此足以追踪原子振动的具体变化。但是这一优势同时限制了MD在大时间尺度模拟上的应用。现有的计算条件足以支持MD到10ns,运用特殊的算法可以到达10的尺度。即便如此，很多动态过程，如外表生长或材料老化等，时间跨度均在S以上，大大超出了MD的应用范围。有什么方法可以克服这种局限呢？当体系处于稳定状态时，我们可以将其描述为处于NV维势能函数面的一个局域极小值(阱底)处。有限温度下，虽然体系内的原子不停的进行热运动，但是绝大局部时间内原子都是在势能阱底附近振动。偶然情况下体系会越过不同势阱间的势垒从而完成一次“演化"，这类小概率事件才是决定体系演化的重点。因此，如果我们将关注点从“原子”升格到“体系”，同时将"原子运动轨迹”粗化为"体系组态跃迁”，那么模拟的时间跨度就将从原子振动的尺度提高到组态跃迁的尺度。这是因为这种处理方法抵弃了与体系穿越势垒无关的微小振动，而只着眼于体系的组态变化。因此，虽然不能描绘原子的运动轨迹，但是作为体系演化，其“组态轨迹"仍然是正确的。此外，因为组态变化的时间间隔很长，体系完成的连续两次演化是独立的，无记忆的，所以这个过程是一种典型的马尔可夫过程(MarkoVPrOCess),即体系从组态;到组态力,一j这一过程只与其跃迁速率和j有关。如果精确地知道自j,我们便可以构造一个随机过程，使得体系按照正确的轨迹演化。这里''正确''的意思是某条给定演化轨迹出现的几率与MD模拟结果完全一致(假设我们进行了大量的MD模拟，每次模拟中每个原子的初始动量随机给定)。这种通过构造随机过程研究体系演化的方法即为动力学蒙特卡洛方法(kineticMonteCarlo,KMC)1o指数分布与KMC的时间步长在KMC模拟中，构造呈指数分布的随机数是一个相当重要的步骤。这一节中我们对此进行讨论。因为体系在势能面上无记忆的随机行走，所以任意单位时间内，它找到跃迁途径的概率不变，设为人因此在f,f+Af)区间内,体系不发生跃迁的概率为FStW(Af)=1KotAf+。(Af)之类似的，在核，f+2Uf)区间内，体系不发生跃迁的概率为t-(2f)=(1一fctot+O(A£产)=1_2totf+(f)2以此类推，当T=KAt时，在f+)区间内，体系不发生跃迁的概率为Kg(T)=(-fctot+(c-2)因此，当T趋于3C时，体系不发生跃迁的概率为FStW(T)=Iim-8(I-&Ot+£+0(A',)=exp(-fctot)这一行为类似于原子核的衰变方程。从方程我们可以得到单位时间内体系跃迁概率双打。从方程的推导过程可以看出体系的跃迁概率是一个随时间积累的物理量，因此P(D对时间积分到某一时刻步必然等于1-Eitay(),也即p(£)=a(-尸"伍)/”因此我们立即可以得到田p(0=fctotexp(-tz>tO(2)Kot是体系处于态i时所有可能的跃迁途径的速率和，之和，即对于每个具体的跃迁途径自j,上述讨论均成立。因此，我们可以定义单位时间内体系进行ij跃迁的概率Pij(。为Pij(t)=k,jexp(-kijt)(4)单位时间内体系的跃迁概率呈指数分布这一事实说明KMC的时间步长田也应是指数分布。因此我们需要产生一个指数分布的随机数序列。这一点可以非常容易的通过一个(0,1平均分布的随机数序列转化得到：F=1Ktay(ZtoMf)从而=-ln(1-r)=-ln(r)最后一步是因为1-r和T的分布相同。粒也可以通过上述步骤从方程得到。计算跃迁速率过渡态理论(TST)Lij决定了KMC模拟的精度甚至准确性。为避开通过原子轨迹来确定3j的做法(这样又回到了MD的情况)，一般情况下采用过渡态理论(transitionstatetheory,TST)进行计算2o在TST中，体系的跃迁速率决定于体系在鞍点处的行为，而平衡态(势阱)处的状态对其影响可以忽略不计。如果大量的相同的体系组成正则系综，则在平衡状态下体系在单位时间内越过某个垂直于1T/跃迁途径的纵截面的流量即为自上简单起见，假设有大量相同的一维双组态(势阱)体系，平衡状态下鞍点所在的假想面(对应于流量最小的纵截面)为工=。，则TST给出该体系从组态A迁出到B的速率为5,6a-b=j(rf(j-g)i)A方程(6)中)A表示在组态A所属态空间里对正则系综的平均。!表示只考虑体系从组态A迁出而不考虑迁入A的情况(后一种情况体系也对通过纵截面的流量有奉献)。根据普遍公式/)_0e靖TdXdP设体系的哈密顿量“为P22m+1'(工)，即可分解为动能和势能，同时设粒子坐标ZWq时体系处于组态A。则方程可写为_'dpJ+五附工一。)团expp32m+l'(N)LTAT-XXrXa<xpP22m+V(x)fcT_1/XX同厂了二川八日“、心-q)La)%focace-tj,23m)fcdpI1e_yq)/*1Jretr=bW>("±-q)A/x1/2=l(¾f)M(Z-G)A(7)上式中无限小量£是为了将A函数全部包含进去。最后一项对于3函数的系综平均可以直接通过MetropolisMonteCarlo方法计算出来：计算粒子落在卜一3,9+3范围内的次数相对于MetroPoIiS行走总次数的比例/瓦方程最后等于fcAB=()v2()(8)将上述讨论扩展到3维情况非常直接，这里只给出结果，详细讨论请参阅文献5:心-B="2f)同/(RHIWI)A其中/(R)是纵截面方程，el代表3维情况中粒子流动方向与截面/法向不平行对于计数的影响。简谐近似下的过渡态理论(hTST)虽然上一节已经给出了TST计算跃迁速率的方法，但是在具体工作中，曲j更多地是利用简谐近似下的过渡态理论(harmonicTST,hTST)通过解析表达式给出。根据TST,跃迁速率依，为3自j=exp-AFjj/(LBT)(io)其中F,j为在跃迁ij中体系在鞍点和态i处的自由能之差AHj=Effd-TS乎TELTsJ=AE1J-T5b将上式代入方程(10),可以得到%=军exp(-A%/强T)exp(S*b/逅)(11)hTST认为体系在稳态附近的振动可以用谐振子表示，因此其配分函数是经典谐振子体系的配分函数。分别写出体系在态i和鞍点处的配分函数Zo和ZSad:/,3-l3-1ZSad=11根据Boltzmann公式，S=kInZ(12)并将配分函数代入，则方程(11)得3JV(13)n*=舟一exp11*=1方程(13)在通常的文献上经常可以见到。声子谱可以通过HeSSian矩阵对角化或者密度泛函微扰法(DFPT)求出，而AE”就是iJ的势垒，可以通过NEB或者drag方法求出。因此，方程(13)保证了可以通过原子模拟(MD或者DFT方法)解析地求出岛人事实上这个方程有两点需要注意。首先虽然方程(10)中出现了普朗克常数/?.，但是在最终结果中/l被抵消了。这是因为TST本质上是一个经典理论，所以充分考虑了统计效应后/l不会出现1»其次，方程(13)说明对于每一个跃迁过程，鞍点处的声子谱应该单独计算。这样会大大增加计算量，因此在绝大局部计算中均设前置因子为常数，不随跃迁过程而变化。具体数值取决于体系，对于金属而言，一般取IO12Hz。KMC几种不同的实现算法点阵映射到目前为止，进行KMC模拟的所有理论基础均已具备。但是前面所进行的讨论并没有联系到具体的模型。KMC在固体物理中的应用往往利用点阵映射将原子与格点联系起来。从而将跃迁(事件)具象化为原子一格点关系的变化。比方空位(团)/吸附原子(岛)迁移等等。虽然与实际情况并不完全一致，但这样做在很多情况下可以简化建模的工作量，而且是非常合理的近似。很多情况下体系中的原子虽然对理想格点均有一定的偏离，但是并不太大(001o),因此这种原子一点阵映射是有效的。这种做法的另一个好处是可以对跃迁进行局域化处理。每条跃迁途径只与其近邻的体系环境有关，这样可以极大的减少跃迁途径的数目，从而简化计算lo需要指出的是，这种映射对于KMC模拟并不是必须的。比方化学分子反应炉或者生物分子的生长等等，这些情况下根本不存在点阵。无拒绝方式KMC的实现方法有很多种，这些算法大致可以分为拒绝(rejection)和无拒绝(rejection-free)两种范畴。每种范畴之下还有不同的实现方式。本文只选择几种最为常用的方法加以介绍。I.直接法直接法(directmethod)是最常用的一种KMC算法，其效率非常高。每一步只需要产生两个在(°，门之间平均分布的随机数”和2。其中Tl被用来选定跃迁途径，2确定模拟的前进时间。设体系处于态£,将每条跃迁途径,想象成长度与跃迁速率Qj成正比的线段。将这些线段首尾相连。如果riKot落在线段Jk中，这个线段所代表的跃迁途径就被选中，体系移动到态Jk,同时体系时间根据方程(5)前进。总结其算法如下：1. 根据方程计算体系处于态:时的总跃迁速率ha=23M;2. 选择随机数-l.ktjTiktQtlt;£%3. 寻找途径Jk,满足5=1J=I；4. 体系移动到态Jk,同时模拟时间前进"=ln(r2)15. 重复上述过程。需要指出的是，虽然一般步骤4中的科根据方程(5)生成，但是如果将其换为擀=自并不会影响模拟结果。在文献5和6中均采用这种方式。II.第一反应法第一反应法(firstreactionmethod,FRM)在思路上比直接法更为自然。前面说过，对于处于稳态i的体系而言，它可以有不同的跃迁途径J可以选择。每条途径均可以根据方程给出一个指数分布的“发生时间"Bf”,也即从当前算起iTj第一次发生的时间。然后从3f,j中选出最小值(最先发生的"第一反应")，体系跃迁到相应的组态Jmin,模拟时间相应地前进左5”;“。总结其算法如下:1. 设共有A/条反应途径，生成A/个随机数n,r2,rM;2. 根据公式优'=一亡由代力，给出每条路径的预计发生时间；3. 找出解，的最小值能兀；”；4. 体系移动到态Jmin,同时模拟时间前进Mijm;“；5. 重复上述过程。可以看出，这种算法的效率比直接法低下，因为每一步KMC模拟需要生成A/个随机数。通常情况下KMC模拟需要IO7步来到达较好的统计性质，如果每一步都需要生成A/个随机数，则利用这种方法需要一个高质量的伪随机数发生器，这一点在A/比较大时尤为重要。III .次级反应法次级反应法(nextreactionmethod,NRM)是FRM方法的一种衍生方法，其核心思想是假设体系的一次跃迁并不会导致处于新态的体系对于其他跃迁途径的舍弃(比方充满可以发生A1种化学反应的分子，第一种反应发生并不会造成别的反应物的变化)，这样体系还可以选择储八中的次小值能，从而跃迁到态,2nd,模拟时间前进也小小一能工如果这次跃迁还可以满足上述假设条件，再重复上述过程。理想情况下，平均每一步KMC模拟只需要生成1个随机数。这无疑会大大提高效率以及时间跨度。但是实际上NRM的假设条件很难在体系每次跃迁之后都得到满足，在固体物理的模拟中尤其如此，因此其应用范围集中于研究复杂化学环境下的反应过程。试探-接受/拒绝方式这一大类算法虽然在效率上不如直接法，但是它们所采用的试探-接受/拒绝在形式上更接近MetropolisMC方法，而且可以很方便的引入恒定步长，即兴固定。因此有必要进行详细的介绍。IV .选择直接法选择直接法在决定体系是否跃迁方面和MetropolisMC方法形式上非常相像，均是通过产生随机数和预定的阈值比较决定事件是否被采纳。具体算法如下：1. 设共有A/条反应途径，选择反应速率最大值ArmaX,设为比生成在0.均匀分布的随机数2. 设j=Int(r)+L3. 如果j-r<3jh,则体系跃迁至新态工否则保持在态K4. 模拟时间前进摄=ln(0s5. 重复上述过程。这种方法的长处在于每一步只需要生成一个随机数。但是缺点也很明显，对于反应速率相差太大，尤其是只有一个低势垒途径(与其他途径相比自，过大)的体系来讲，这种方法的效率会非常低下。某些情况下，这种低效率问题可以通过如下方法改良：将全部途径按照Aij的大小分为几个亚组，每个亚组选定一个上限京h但是这一步骤在整个KMC模拟过程中可能需要重复很屡次，因此并不能完全解决问题。事实上低势垒在KMC中是个普遍的问题。这一点在后面还要简要提及。V.恒定步长法与上述四种方法不同，恒定步长法(COnStanttimestepmethod,CTSM)中体系的前进时间是个给定的参数citedawnkaski。在理想情况下，CTSM与直接法效率相同，每一步只需产生两个随机数。具体算法如下：1. 给定恒定时间步长战；2. 将所有途径，(共有M个)设为长度恒为1/M的线段，生成在l°1均匀分布的随机数ri,选择途径j=Int(rM)+l;3. 生成在°，1均匀分布的随机数万，如果2自#,则体系跃迁至新态,，否则保持在态£；4. 模拟时间前进用；5. 重复上述过程。实际模拟中，务需要满足小于(见''第一反应法")，以及(2)对于自j最大的途径，接受率大致在0.5。其中第一个条件保证了所有的迁移途径发生概率都小于1,第二个条件则保证体系演化的效率不会过于低下。CTSM是非常行之有效的一类KMC算法，但是选择济时需要特别的注意以保证效率。自决定于具体体系以及模拟温度。这在一定程度上增加了CTSM的实现及使用难度。低势垒问题前面已经指出，低势垒的途径需要特别注意。如果体系在演化过程中一直存在着势垒较其他途径低很多的一个或几个途径，会对模拟过程产生不利的影响。这个问题被称之为低势垒问题。低势垒途径对于KMC模拟最直接的影响就是大大缩短了模拟过程所涵盖的时间跨度。这一点可以从方程(5)中看出。更为深刻的影响在于，这些由低势垒的途径联系起来的组态会组成一个近似于封闭的族。体系会频繁的访问这些态，而其他的对于体系演化更为重要的高势垒途径被选择的概率非常低，这显然会降低KMC的模拟效率。例如，吸附原子在高指数金属外表扩散，其沿台阶的迁移所对应的势垒要远低于与台阶别离的移动。这样，KMC模拟的绝大局部时间内吸附原子都在台阶处来回往复，而不会选择离开台阶在平台上扩散。这显然不是我们希望看到的情形。一种解决方法是人为地将这些低势垒加高以降低体系访问这些组态的几率，但是无法预测这种干扰是否会造成体系对于真实情况的严重偏离。另一种选择是利用NRM或者CTSM进行模拟，但是其效果如何尚待检测。如果考察体系的势能面，这类低势垒的途径一般处在一个超势阱''之中。体系在这个超势阱中可以很快的到达热平衡，所需时间要短于从其中逸出的时间。如果可以明确的知道超势阱所包含的组态以及从超势阱逸出的所有途径，我们就可以按照Boltzmann分布合理的选择其中一条途径，使得体系向前演化。但是如何确定哪些组态包含在超势阱之中以及体系是否已在其中到达热平衡本身就是两个难题。对于第一点，MaSon提出可以利用ZobriSt密钥法标定访问过于频繁的组态7；NoVOtny则提出通过建立及对角化一个描述体系在这些组态间演化的传递矩阵来解决第二点8»对这个问题的详细讨论已超出了本文的讨论范围，请参阅文献7以及8。实体动力学蒙特卡洛方法OKMC上述的KMC都假设任何时候原子均处于其理想点阵格子上。但是很多情况下这种点阵映射是无效的，比方间隙原子或者位错。这类结构缺陷的运动在材料的辐射损伤和老化过程中扮演着非常重要的角色。而且与单个原子或者空位的运动相比，这类缺陷的运动时间跨度更长，也更为复杂，比方间隙原子团和空穴的湮没，间隙原子团的解构/融合，或者位错的攀移/交滑移等等。传统的KMC算法很难有效的处理这类问题，一方面是因为时间跨度太大，另一方面这类缺陷各自均可视为独立的实体(ObjeCt),其运动更近似于系统激发，因此单个或几个原子运动的积累效果很多情况下并不能有效地反应这些实体的整体运动。实体动力学蒙特卡洛方法(ObjectkineticMonteCarlo,OKMC)就是为了处理这类问题而被提出的。OKMC在算法上与普通的KMC完全一样。需要注意的地方是在OKMC中并不存在原子点阵。所有的实体在一个真空的箱子中按照其物理实质离散化运动，比方位错环的最小移动距离是其BUrgerS矢量大小，方向则为BUrgerS矢量方向；空位的移动距离为第一近邻或第二近邻的原子间距，等等。模拟过程中我们需要追踪该实体的形心，从而决定其位置、移动距离等等。此外，OKMC中对于跃迁速率A确实定也和普通的KMC有所区别。本文前面已经指出，可以表达为自j=exp(AEij/LBT)的形式。普通的KMC假定为常数，不同途径的岛j由AE,j决定。但是在OKMC的模拟中,AE,j的直接确定非常困难，因此一般的策略是对于特定的事件(包括实体自身的运动以及不同实体间的反应等)，跃迁势垒AE”保持恒定，而将前置因子视为实体规模(所包含的原子/空位数目)的函数，通过MD模拟得出，一般而言可以表示为形如M)(m)=M)(q-)m-a的表达式，其中q和S是拟合参量，nr是实体规模。最后需要注意的是在OKMC的模型中，实体有空间范围，因此需要一个额外的参数凡ntity来表征其空间半径(假设为球形分布，否则凡ntity的数目多于一个)。在模拟不同实体间的反应时，需要特别考虑其形心的间距，如果小于''反应距离"，即反应一定进行，否则认为两个实体互相独立。DOnIain利用OKMC研究了Fe-CU合金的辐射损伤9,在模拟中考虑了间隙原子(空位)的聚合、间隙原子(空位)团的发射、间隙原子-空位湮没、空位团对杂质的捕获、外表对于空位(团)的捕获、甚至辐射轰击引起的间隙原子(空位)萌生、增殖等等事件。从中可以看出，对于OKMC,一个棘手的问题是需要预先想到所有的事件。此外，0KMC所需要的所有参量根本上不可能通过原子模拟直接获得，人为的设定参数不可防止。这些参数会在多大程度上决定OKMC的准确程度无法预先得知。需要根据现有的实验数据进行修改、调试。这些困难都限制了OKMC的普及。但是如前所述，这种方法可以有效地进行大尺度的时间(天)和空间模拟Mn以上)，而且对于缺陷的描述更为直接和符合直观，因此在材料研究中同样占有重要的地位。KMC的若干进展等时蛙跳算法(rTeapKMC)引入这类算法前，我们先简要介绍两个常用的离散分布：泊松分布(POiSSe)nDistribution,PD)以及二项式分布(BinOmialDistribution,BD)。泊松随机数尸(向j，)定义为给定事件发生率自j以及观测时间T下事件发生的数目。如果用“代表给定的发生数H,则尸(Lij，。恰好等于n的概率是一个泊松分布：P(n;M,)=PrPkij,)=n=。"，普"(14)也即如果产生一个泊松随机数序列0(k”,),则这个序列符合泊松分布PD。需要指出，P(Lj,)是无界的，范围是任意非负整数。与其类似，二项式随机数6(N.p)定义为重复N次独立的成功率均为P的伯努利实验的成功数。如果给定成功数n,则6(NN恰好等于m的概率是一个二项式分布：BgN,p)=Pr6(N,p)=n=7j11pn(ln=0,1.，N(为了和本文中的标号一致，我们将跃迁i-j的成功率P表示为¾“v,将方程(15)重新写为8(n;NMm=Ny(I-%N)Nf(与PD不同，BD中的门是有界的，为O到N之间的任意整数。可以看出，如果将这两种随机数理解为给定跃迁路径(发生率为ej)在一定的时间步长(T)内发生的次数，则可以立即运用于粒子数空间N内的KMC中，其时间范围可以得到很大提高。这就是等时蛙跳算法T-IeaPKMC10,11。rTeapKMC方法最早由GiIleSPie提出，通过pd方程(14),在给定时间步长r下决定每个跃迁途径，发生的次数nj,然后将体系移到这些跃迁累计发生后产生的新态。因为每一步模拟体系不止发生一次跃迁，所以模拟的速度可以大大加快。我们以多种反应物在化学反应炉中的演化为例加以详细说明。设在炉内共有L种分子85l,在时亥加各自的个数为X"。，则在粒子数空间N中X(f)XMf)T构成一个矢量x(a,或称为一个组态K总共有a/种反应路径Aj(j=1对于给定的Rj,反应速率后：x(f)是占据态(a的函数。此外，我们单独定义一个矢量Vj=X1-X0,其中i由xO通过反应Aj而得，即x02xl因此4的元素Wj代表反应R7所引起的多种分子的数目变化。由此建立算法如下：VI.PD-F-IeapKMC101.给定恒定时间步长丁;2.对于每条反应途径Aj按照方程(14)生成泊松随机数序列,按照模拟步数从序列中找出每种反应发生的次数否;MX(f十T)=X+£与叼3. 按照5=1更新体系;4. 模拟时间前进r；5. 重复上述过程。GilleSPie仔细考虑了丁的选择条件，称为蛙跳条件(IeaPcondition)：=minJ1.43凶4ujX)j(X)2(X)(17)其中m(X)=-J,m=l.，M1=1Mj(X)=fjm(X)km(X),j=1,M711=lM响X)=(X)fcn(X),j=1,M如前所述，尸伙5丁)没有上限，因此即使r满足方程(17),在模拟过程中也可能会出现某种分子总数为负数的情况，这显然不符合实际，也是PD-TTeapKMC的一个弱点。Tian和BUrrage提出可以用二项式分布BD取代PD,因为B(N,p)有上限，所以可以有效的解决这个问题。此外，他们对于某种分子参与多种反应的情况也进行了考虑，从而提高了TTeapKMC的稳定性和普适性。其算法如下：VII.BD-T-IeapKMC11O<ICJ(X)T<1给定恒定时间步长丁，满足U二Nj-1；对于每条反应途径RJ按照方程(16)生成二项式随机数序列®M'L'X"",按照模拟步数从序列中找出每种反应发生的次数,叮；如果有某种分子Sr同时参与了Hj和R,n,则首先生成叩J小J(X)+%(X)然后通过"m=njm-Rj确定RnI的发生次数;MXH十T)=X+£与否按照J=I更新体系;模拟时间前进r；重复上述过程。步骤1、2中出现的NJ一是参与反应Rj的各类分子的个数的最小值，即Nj=min.me此外Gillespie,Tian和Burrage还考虑用预测17?时刻体系状态的方法来进一步提高精度。具体请参阅文献10,11。如果TTeaP算法和OKMC结合起来可以进一步加大模拟的时间尺度，但是目前还没有这方面工作的介绍。基于即时动态分析的KMC方法(On-the-flyKMC)到目前为止，所有的KMC都是在模拟之前建立好所有可能的跃迁途径。但是实际上所有"是很难到达的目标。因为很多途径远离一般的直觉，而且在演化过程中体系有可能寻找到新的途径。因此，跃迁途径应该随着体系的演化而不断更新，是动态的过程。Henkehnan和J6nsson将途径搜索和KMC结合起来，提出了即时动态的KMC方法on-the-flyKMC12：在每一个稳态(势阱)处，选定一个激活原子(一般是近邻不饱和的原子)，在以其为中心的局部区域内引入呈高斯分布的随机位移，即参加扰动，然后利用dimer方法13寻找所有可能的跃迁途径。建立起即时的途径库之后再通过普通KMC算法进行模拟。显然，这种方法的计算量非常大，需要一个有效的标识方法来识别所有已经遇到过的途径以防止重复计算。TrUShin提出可以利用包括至激活原子第三壳层的所有格点(顺时针排列)的占据与否(分别标记为1和0)来构建二进制数，从而根据始态和终态的标号来唯一地标识某条途径14,例如，激活原子标为"1",其第一壳层的原子标记为"2","3",.,''AY',依此类推，然后将原子的标号”作为二进制的数位2工，这样，每一个稳态都有唯一的一个二进制数与之对应。虽然仍不完善，但是这种方法具有非常清晰的逻辑结构，具有良好的扩展性。OaQg2N)和O(I)KMC方法一般情况下KMC的大局部时间花费在选择途径上。如果采用普通的方法，即循环叠加自，直至°3=1”从而选择Jk,这种情况下计算用时与途径数目M呈线性增长，即O(N)算法。按照二叉树安排不同数目的柘之和可以改良到5bg2N)15；将所有自j作为树叶(缺乏2整数次嘉的叶子由0填补)，每两片叶子之和作为父节点，依次类推直至树根卜ct。一株二叉树构建完毕后，生成一个随机数ri,由树根开始寻找S=rRot,若S不大于左子节点向Ofb沿左分支向下寻找；否则设6=S-Koft,沿右分支向下寻找，直至树叶上体系按途径亍演化。Slepoy和ThOnIPSon等进一步提出分流-拒绝(COnIPoSitiOn-rejection,CR)方法以实现搜索用时与途径总数无关的°(D算法16：(1)先找出*min和Lmax,按照(Amin,2Ain,4hnin,及max)将A/条途径分为G个组,g=1°S2(,(2)然后生成随机数”，按照上述二叉树寻找6=rRst所落入的组别Gj,再生成两个随机数2和3,设'=Int(r2NGJ,其中NG,为该组中包含的途径数，f=rax,如果fh,则选择途径£,否则重复步骤(3),直至有一条途径被选中为止。可以看出，CR算法虽然搜索速度很快，但是每一步KMC需要产生至少4个随机数(以用于确定前进时间)，因此需要高质量的随机数发生器。不过对于跃迁途径复杂的体系演化而言，CR的°(D效率无疑是很有吸引力的。1 A.F.Voter,itRadiationEffectsinSolids(Springer2006)p.1-24.2 H.Eyring,J.Chem.Phys.3,107(1935).3 P.Kratzer,MultiscaleSimulationMethodinMolecularScience(NICSerices,Vol.42,Forschungszentrum,Julich2009)p.51-76.4 E.J.Dawnkaski,D.SrivastavaandB.J.Gamson,J.Chem.Phys.102,9401(1995).5 A.F.VoterandJ.D.Doll,J.Chem.Phys.80,5832(1984).6 A.F.Voter,Phys.Rev.B34,6819(1986).7 D.R.Mason,T.S.HudsonandA.P.Sutton,Comp.Phys.Comm.165,37(2005).8 M.A.Novotny,Phys.Rev.Lett.74,1(1994);Erratum75,1424(1995).9 C.Domain,C.S.BecquartandL.Malerba,J.Nucl.Mater.335,121(2004).10 D.T.Gillespie,J.Chem.Phys.115,1716(2001).11 T.TianandK.Burrage,J.Chem.Phys.121,10356(2004).12 G.HenkelmanandH.J,onsson,J.Chem.Phys.115,9657(2001).13 G.HenkelmanandH.J,onsson,J.Chem.Phys.Ill,7010(1999).14 0.Trushin,A.Karim,A.KaraandT.S.Rahman,Phys.Rev.B72,115401(2005).15 M.A.GibsonandJ.Bruck,J.Phys.Chem.A104,1876(2000).16 A.Slepoy,A.P.ThompsonandS.J.Plimpton,J.Chem.Phys.128,205101(2008).资料来源量子化学网