动力学蒙特卡罗模拟方法简介

踪原子振动的具体变化。
但这也限制了其在大时间尺度模拟上的应用（现有计算条件可支 持时间步长达到10ns，运用特殊算法可达到10μ s，但很多动态过程 的时间跨度在秒数量级以上）
• 体系处于稳定状态时，可将其描述为处于3N维势能函数面的一个局域 最小值（势阱底）处。
• 有限温度下，虽然体系内的原子不停进行热运动，但绝大部分时间内
j 1
2.2 第一反应法 对处于稳态i的体系，其每条跃迁途径j均可给出一个指数分布的“发生时间”δtij，即 从当前算起i→j第一发生的时间。然后从{δtij}中选出最小值，体系跃迁到相应的
组态jmin，模拟时间相应地前进δtijmin。
算法： (1)设共有M条反应途径，生成M个随机数r1,r2,„,rM;（效率） (2)计算出每条路径的预计发生时间； (3)找出{δtij}的最小值δtijmin；
其中<...>A表示在组态A所属态空间里对正则系综的平均。1/2表示只考虑体 系从组态A迁出而不考虑迁入A的情况。根据普遍公式：
H ˆ exp k T dxdp B ˆ H exp k T dxdp B
2、简谐近似下的过渡态理论
2 Pstay lim1 ktot K exp ktot K

k
pt ktot exp ktott
其中ktot是体系处于组态i时所用可能的跃迁途径的速率kij之和。
因此，对于单位时间内体系进行某一个具体的跃迁途径kij的概率则可定义为：
则系综，则在平衡状态下体系在单位时间内越过某个垂直于i→j跃迁途径的
纵截面的流量即为kij。
假设有大量相同的一维双组态（势阱）体系，平衡状态下鞍点所在的假想
面（对应流量最小的纵截面）为x=q，则过渡态理论给出该体系从组态A迁
出到组态B的速率为：
k A B
1 x q x 2
A
算法： (1)计算体系处于组态i时的各条路径跃迁速率kij，以及总跃迁速率ktot； (2)选择随机数r1； (3)寻找途径jk，满足 kij r1ktot kij ；
jk 1 j 1 jk
1 (4)体系移动到态jk，同时模拟时间前进 t Inr2； ktot
(5)重复上述步骤。
2、无拒绝方法
直接法、第一反应法、次级反应法等。
2.1 直接法
效率高，最常用
每一步需要产生两个在(0,1]上平均分布的随机数r1和r2，分别用于选定跃
迁途径和确定模拟的前进时间。设体系处于态i，将每条跃迁途径j想象
成长度与跃迁速率kij成正比的线段。将这些线段首尾相连。如果r1ktot落
在线段jk中，这个线段所代表的跃迁途径jk就被选中，体系移动到态jk， 同时体系时间根据时间步长方程前进。
目
录
1 2 3
4
KMC 的基本原理
指数分布与 KMC 的时间步长
跃迁速率的计算 KMC 的实现算法
1、过渡态理论
跃迁速率决定了KMC模拟的精度甚至准确性。为避开通过原子轨迹来确定 kij的做法，一般采用过渡态理论进行计算。 过渡态理论中，体系的跃迁速率取决于体系在鞍点处的行为，而平衡态 （势阱）处的状态对其影响很小，可以忽略。如果大量相同的体系组成正

1
sad ij
结合玻尔兹曼公式 S k B InZ
3N
可得：
可通过原子模拟（MD算法或
kij
i 1 3 N 1 i 1
0 ij
sad ij
Eij exp k T B
DFT方法）解析求出kij。 前置因子设为常数。（金属： 约1012Hz）
1 Inr； ktot
其是只有一个低势垒途径
的体系来讲，效率很低。
实际模拟中，δt需满足：
(1)小于δtmin，以保证所有的
迁移途径发生的概率都小 于 1； (2)对于kij最大的途径，接 受率大致为50%，以保证 体系演化的效率不会过低。
3.2 恒定步长法（CTSM）
前进时间是给定的参数 理想情况下，其效率与选择路径法相同，每一步只需要产生两个随机数 算法： (1)给定恒定时间步长δt；
Pstay t 1 ktot t t
2
在区间[t,t+2△t)上，体系不发生跃迁的概率：
Pstay 2t 1 ktot t t

2 2

1 2ktot t t
2
以此类推，当τ=K△t时，在区间[t,t+τ)上，体系不发生跃迁的概率：
(2)将所有途径j(共有M个)设为长度恒为1/M的线段，生成在区间[0,1]上均匀分布的随
机数r1，选择途径j=INT(r1M)+1； (3)生成区间[0,1]上均匀分布的随机数r2，如果r2＜kijδt，则体系跃迁至新态j，否则保 持在态i； (4)模拟时间前进δt； (5)重复上述过程。
Thank You
算法：
ˆ 。生成在[0,M) (1)设共有M条反应途径，选择反应速率最大值kmax，设为 k
区间内均匀分布的随机数r；
(2)设j=INT(r)+1;
• 每一步只需要生成一个随 机数； • 对反应速率相差太大，尤
ˆ ,则体系跃迁至新态j，否则保持在组态i； (3)如果j-r＜kij/ k
(4)模拟时间前进 t (5)重复上述过程。


简谐近似下的过渡态理论认为体系在稳态附近的振动可以用谐振子表示， 故可视为经典谐振子体系。则体系在态i和鞍点处的配分函数Z0和Zsad为：
k BT 0 Z h
3N 3N

i 1
1
0 ij
Z
sad
k BT h
3 N 1 3 N 1 i 1
全一致）
• 这种通过随机过程研究体系演化的方法即为KMC方法。
目
录
1 2 3
4
KMC 的基本原理
指数分布与 KMC 的时间步长
跃迁速率的计算 KMC 的实现算法
体系在势能面上无记忆地随机行走，因此其在任意单位时间内找到跃迁途径 的概率是恒定的，设为ktot。则： 在区间[t,t+△t)上，体系不发生跃迁的概率：
目
录
1 2 3
4
KMC 的基本原理
指数分布与 KMC 的时间步长
跃迁速率的计算 KMC 的实现算法
目
录
1 2 3
4
KMC 的基本原理
指数分布与 KMC 的时间步长
跃迁速率的计算 KMC 的实现算法
分子动力学在原子模拟领域具有突出的优势。其可以精确描述体 系演化的轨迹。分子动力学的时间步长通常在飞秒数量级，这足以追
2 Pstay 1 ktot K K

k
故：当τ区域∞时，体系不发生跃迁的概率为：
由此即可得到单位时间内体系跃迁的概率p(t)。由之前的推导过程可知，体 系的跃迁概率是一个随时间积累的物理量，因此p(t)对时间积分到某一时刻 t'必然等于1-Pstay(t')，即 pt 1 Pstay t t ，于是有：
都是在势阱底附近振动。 • 偶然情况下，体系会越过不同势阱间的势垒而完成一次“演化”（决定 体系演化的重点） • 关注点：原子 体系
粗粒化
• 原子运动轨迹
体系组态跃迁
模拟的时间跨度
• 组态变化的时间间隔很长，完成的连续两次演化是独立的、无记 忆的，因此其为一种Markov过程，即体系从组态i到组态j这一过 程只与其跃迁速率kij有关。 • 精确知道kij，便可构造一个随机过程，使得体系按照正确的轨迹 演化（“正确”是指某条给定演化轨迹出现的概率与MD模拟结果完
较选择路径法更自然，但效率更低 通常KMC模拟需要107步来达到较好的统
(4)体系移动到态jmin，同时模拟时间前进δtijmin； 数，则利用这种方法需要一个高质量的伪 (5)重复上述过程。
随机数发生器，M较大时尤为重要。
计性质，如果每一步都需要生成M个随机
2.3 次级反应法 假设体系的一次跃迁并不会导致处于新态的体系对于其他跃迁途径的取舍 （比如充满可以发生M种化学反应的分子，第一种反应发生并不会造成别 的反应物的变化），这样体系还可以选择{δtij}中的次小值δtij2nd，从而跃迁 到态j2nd，模拟时间前进δtij2nd-δtij2nd。如果此次跃迁还可以满足上述假设，再 重复此过程。
理想情况下，平均每一步KMC模拟只需要生成一个随机数，因此能大大提
高效率并加大时间跨度。
应用范围集中于研究复杂化学环境下的反应过程。
3、试探-接受/拒绝方法
效率低于无拒绝方法，但其形式更接近蒙托卡 罗方法，且可方便地引入恒定步长，即δt固定。
3.1 选择路径法
选择路径法在决定体系是否跃迁方面与蒙托卡罗方法形式上非常相像， 均是通过产生随机数和预定的阈值比较决定事件是否被采纳。
pij t kij exp kijt
即单位时间内体系的跃迁概率呈指数分布，这说明KMC的时间步长δt也呈
指数分布，因此需要产生一个按指数分布的随机数序列：
通过一个在区间(0,1]上平均分布的随机数序列r转化得
r 1 Pstay ktott
由于1-r和r的分布相同，从而有 1 1 t In1 r Inr ktot ktot
根据过渡态理论，跃迁速率为：
k BT kij exp Fij k BT h
其中


vib ij ，表示跃迁
Fij E
sad ij
TS
vib ij
Ei TSi Eij TS
i→j中体系在鞍点和态i处的自由能之差。 于是：
k BT vib kij exp Eij k BT exp Sij kB h
பைடு நூலகம் 目
录
1 2 3
4
KMC 的基本原理
指数分布与 KMC 的时间步长
跃迁速率的计算 KMC 的实现算法
1、点阵映射
点阵映射
KMC在固体物理领域的应用中，常利用点阵映射将原子与格点联系起来，
从而将跃迁（事件）具象化为原子 格点关系的变化。
与实际情况不完全一致，但很多情况下都可以简化建模工作量，且是非 常合理的近似。很多情况下体系中的原子虽然对理想格点有一定偏离， 但并不太大（约0.01a0），因此这种映射是有效的。 可以对跃迁进行局域化处理。每条跃迁途径只与其近邻的体系环境有关， 这样可以极大地减少跃迁途径的数目，从而简化计算。 非必需