kinetic Monte Carlo : Self Study Note
Masahiro YamamotoEHCC, Kyoto Univ.
Modiﬁed on March 6, 2009 9 : 28 pm
The method to do the kinetic Monte Carlo simulation is reviewed.

1

通常平衡状態の集団平均をとるために行われる Monte Carlo シミュレーシ
ョンで時間に依存した現象を解析できるのか？

答えは YES です。しかも、分子動力学計算（ps, ns レベル）では対象外の”たまに起きる現象 の kinetics を扱うのに適してます。

1.1

平衡状態のモンテカルロシミュレーション:
復習

In an irreducible Markov chain (An irreducible or ergodic chain is one
where every state can eventually ve reached from another state.), thesteady state probability becomes∑ρi pij=

ρjpij>

0

pij

=1j

∑
j

The unnecessarily strong condition of ’microscopic reversibility’ is usally assumed
ρi pij=ρj pjibecause∑ρi pij

Figure 0: Schematic diagram of MC step and the Metropolis routine∑=iρj pji = ρji∑pij = ρj

1.2(1)N 個の粒子から 1 個をランダム (あるいは 1 から N まで順番) に選び出す。その粒子の番号を i とし, その位置ベクトルを ri = (xi , yi , zi ) とする。粒子 i を ri から ri + ∆ri へ移動する. 変位ベクトル ∆ri = (∆xi , ∆yi , ∆zi )は区問 (0,1) の一様乱数 ξx , ξy , ξz を用意し,

For this condition Metropolis et al. [J. Chem. Phys. 21 (1953) 1087]
proposed{pij=pii=(i ̸= j,(i ̸= j,αij
αij ρj /ρi1−∑andρj /ρi ≥ 1)andρj /ρi < 1)

∆xi

pij

And αij should satisfy
αij

=αji

αij

=1j

If we consider the move of a particle l (Energy is Ei ) at rl to the newposition is the rl + δrl in the R region, we can deﬁne
αij==δ(1 − 2ξx ), ∆yi = δ(1 − 2ξy ), ∆zi = δ(1 − 2ξz )で与える. こうして ri を中心とし直交する 3 辺が x, y, z 座標軸に平行な,一辺が 2δ の立方体内に一様分布する点から点 ri + ∆ri をサンプルする.
(2) 粒子 i の移動後の系のポテンシャルエネルギーを Uk とする.
(i)Ui > Uk であれば, 状態 k を n + 1 回目の状態に採用し (3) へ跳ぶ.
(ii)Ui < Uk であれば, 区間 (0,1) の一様乱数 ξ が ξ < exp(−β[Uk −Ui ]) をみたせば, 状態 k を n + 1 回目の状態として (3) へ跳ぶ. みたさな
ければ n + 1 回目の状態も i のままで (3) へ進む.
(3)n → n + 1 として (1) に戻る.
(i)(ii) をまとめて表せば ξ < min[1, exp{−β(Uk − Ui )}] が成立したとき状態 k を n + 1 回目の状態とし, 成立しなければ、n + 1 回目の状態も iとする, ということになる.
このアルゴリズムでは状態 i から状態 k への遷移確率 Pik は (1) 遷移先を選ぶ確率 Aik (とおく) と、(２) 遷移するか否かを決める確率の積で与えられている。確率 Aik は、体積 (2δ)3 の立方体内から 1 点を選ぶ確率であるが、ri と ri + ∆ri はいずれもこの立方体内の点であるから, 逆に状態 k にあるとき状態 i を遷移先に選ぶ確率とも等しい、つまり i, k について対称 (Aik = Aki ) となっているので式 (7.9) の右辺の分母分子で打ち消されていることを注意しておこう. ・パラメータ δ のとりかたによって系の平衡状態への収束が異なる. 経験的には粒子移動後の状態が採用される割合と棄却される割合が同程度になるようにとるとよいという.
この MC 法ではパラメータとして入力される粒子数 N , 体積 V , 温度 Tが一定であるので N V T 集団の MC 法, 略して N V T -MC 法などと呼ば

j(̸=i)∑

Canonical Metropolis Sampling

i

1/NR

Here NR is the total number that the particle can move in the Rregion.
In the general Metropolis routine,
1) Pick up a paticle (sequential or random sampling)2) move the particle to a point in the R region, the state is j’3) If ρ′j /ρi ≥ 1 accept this move.
4) If ρ′j /ρi < 1, generate random number R1 (0 < R1 < 1)
when ρ′j /ρi ≥ R1 accept this move.
when ρ′j /ρi < R1 reject this move.

1

れる. また N 回のサンプリングをサンプリングの単位に選び、1 モンテカ
ルロ (MC) 時間と呼び習わしている.

1) Determination of all processes j that possibly could takeplace with the actual conﬁguration of the system. Here1≤j≤N

1.32) Calculation of the total rate R = j=1 Γ(j) , where the sum
runs over the possible processes [see step 1)]. Depositionis accounted for in this description by the deposition rateF .

∑N

Kinetic (or dynamical) Monte Carlo

理論的な基礎付けは, K. A. Fichtorn and W. H. Weinberg, ”Theoretical foundations of synamical Monte Carlo sumulations”, J.
Chem. Phys. 95 (1991) 1090-1096. にある。
すべての過程はまれにしか起こらない（同時に起こることのない）
Poisson 過程とする。２つのイベントがおこる待ち時間 τ の分布はP (τ ) = R exp(−Rτ ) である。ここで R は total rate である。〈τ 〉 =1/R である。0 から 1 までの一様乱数 r2 を発生させると時間ステップ∆t = − ln r2 /R, r2 = exp(−R∆t) となる。これは，ある時間 t とt + ∆t の間にイベントが起こらない確率が exp(−R∆t) であることに由来する。[数学的な裏付けは、D. A. McQuarrie ”Mathematicalmethods for Scientists and Engineers” Univ. Sci. Books (2003)
pp.1063-1067 を参照のこと。わかりやすく説明してある。]
より具体的に示すために、 Fritz-Haber-Institut der MaxPlanck-Gesellschaft の P. Ruggerone, A. Kley, and M. Scheﬄer”Bridging the length and time scales: from ab initio electronic
structure calculations to macroscopic proportions”の一部を抜粋する。

Figure 0: The diﬀerent atomistic processes for adatoms on a surface:(a) deposition, (b) diﬀusion at ﬂat regions, (c) nucleation of an island,(d) diﬀusion towards and capture by a step edge, (e) detachmentfrom an island, (f ) diﬀusion parallel to a step edge, (g) diﬀusion
down from an upper to a lower terrace, (h) nucleation of an island
on top of an already existing island, and (i) diﬀusion of a dimer (or abigger island). For the processes (a), (c), (g) and (h) also the reverse
direction is possible, but typically less likely.

=(j)

(j)

Γ0 exp(−Ed /kB T )4) Find the integer number l for whichl−1∑

Γ(j) ≤ ρ1 R <

j=1l
∑

Γ(j).

j=15) Let process l to take place. Here 1 ≤ l ≤ N . (l = 1 means0 ≤ ρ1 R < Γ(1) )6) Update the simulation time t := t + ∆t with ∆t =− ln(ρ2 )/R .
7) Go back to step 1).
KMC simulations have been used to study crystal growth ofsemiconductors and metals. However, most of these studies havebeen based on restrictive approximations. For example, the input parameters, such as activation barriers, have been treated aseﬀective parameters determined rather indirectly, e.g. from theﬁtting of experimental quantities. Thus, the connection between
these parameters and the microscopic nature of the processes
may be somewhat uncertain. Often even the surface structure
was treated incorrectly, i.e. the simulation was done on a simple quadratic lattice while the system of interest had an fcc or
bcc structure. Despite these approximations such studies haveprovided signiﬁcant qualitative quantitative insight into growthphenomena. The next better approach is to use semi-empiricalcalculations such as the embedded atom method or eﬀective
medium theory to evaluate the input parameters for KMC simulations of growth. The best, but also most elaborate procedureto obtain these input parameters is to exploit the accuracy ofDFT, and this is our approach. .....

以下の論文を紹介したが、kineticMC の計算の詳細の記述はないので省略します。脱離過程をバトラーボルマーで、チオール間の最近接（引力）相互作用のみが考慮されてます。
論文紹介（山本）2003 Dec. 22√√( 3 × 3)-R30˚アルカンチオールの Au(111) 面のテラスおよびステップからの脱離のモンテカルロシミュレーション[A Monte√ √Carlo simulation for the stripping of the 3× 3R30˚alkanethiol
lattice from Au(111) terraces and steps ]J. Chem. Phys. 117, 2293-2298, 2002
H. Martin, C. Vericat, G. Andreasen, M. E. Vela, and R. C.
Salvarezza
Instituto de Investigaciones Fisicoquicas Teicas y Aplicadas(INIFTA) (UNLP-CONICET)
Sucursal 4, Casilla de Correo
16, 1900 La Plata, Argentina√ Carlo simulations of a model for the stripping of√ Monte3 × 3R30˚alkanethiol lattices from terraces and steps of a(111) metal face in aqueous solutions are presented. In themodel the stripping probability of an adsorbed alkanethiolate
molecule depends on the applied potential, on intermolecularforces that stabilize the alkanethiol layer, and on the presenceof substrate defects. Stabilizing intermolecular forces are also
responsible for alkanethiolate aggregate formation during thestripping process. Snapshots and voltammograms derived fromthe model reproduce experimental STM images and electrochemical data for alkanethiol stripping from the Au(111) surface.

A severe problem in describing crystal growth is the timescale that it is too long for MD simulations. MD runs can provide important insight into elementary microscopic mechanisms
but typically they can cover at most picoseconds, possibly somenanoseconds. Since growth patterns usually develop on a timescale of seconds, the inadequacy of MD is evident. Moreover,the growth structures involve large numbers of particles (from∼ 102 to ∼ 104 ) far beyond the reach of MD (simulations with∼ 102 atoms are hardly feasible, and only for very short times).
Instead, the method of choice for studying the spatial and temporal development of growth is KMC. The key idea behind KMC
is to describe stochastic processes (such as deposition, diﬀusion,desorption, etc.) on the microscopic scales by rates usually given
by:
Γ(j)

3) Choose two random numbers ρ1 and ρ2 in the range (0, 1].

.

(j)Here, Γ0 is an eﬀective attempt frequency, T the temperature,(j)kB the Boltzmann constant, and Ed is the energy barrier thatneeds to be overcome for the event j to take place. In this way
one avoids the explicit calculation of unsuccessful attempts, andKMC simulations therefore can describe phenomena with timescales of seconds. Furthermore, large systems can be treated
without great diﬃculties. The basic procedure of a KMC simulation can be sketched as follows:2