动力学蒙特卡洛方法(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 .次