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