Derivation of the residence time for kinetic Monte Carlo simulations

0.11 MB·English

by  Clinton DeW. Van Siclen

Derivation of the residence time for kinetic Monte Carlo simulations
Clinton DeW. Van Siclen∗
Idaho National Laboratory, Idaho Falls, Idaho 83415, USA
(Dated: November16, 2007)
The kinetic Monte Carlo method is a standard approach for simulating physical systems whose
dynamics are stochastic or that evolve in a probabilistic manner. Here we show how to calculate
thesystem time for such simulations.
8 PACSnumbers: 02.50.Ey,02.70.Uu,05.10.-a
0
0
2 I. INTRODUCTION physical states (for example, the adatom on a finite sur-
face may be foundat sites in proportionto the residence
n
a In a kinetic Monte Carlo (kMC) simulation, a model times at the corresponding states in the model system).
J More generally the kMC method is used to model the
system is moved from state to state according to prob-
3 evolution of a non-equilibrium system (that is, where a
abilistic rules. A nice example is adatom diffusion on a
1 stateonceleftcanneverbereturnedto),sothiscasewill
surface: in this case the “system” is the adatom plus
motivate the derivation of the residence time presented
surface, and a “state” has the adatom at a particu-
] below. However, the mathematics apply to equilibrium
h lar surface site. The probabilistic rules for moving the
p adatomfromonesitetoanother(i.e.,thejumprates)are systems as well, as will be made apparent at the end of
- the derivation.
p obtained from molecular dynamics (MD) calculations,
m which give the energy barriers over which the moves oc-
cur. Subsequently,aseparatekMCsimulationmovesthe
o II. DERIVATION OF THE RESIDENCE TIME
adatom from site to site over the surface, and thereby
c
. enables calculation of the adatom diffusion coefficient,
cs D = x2 /(4t), from distance x traveled by the adatom Consider a collection {j} of elements, where each el-
ement j has an expected lifetime τ . The kMC method
i over the time interval t. j
s (cid:10) (cid:11) follows from the assumption that the values {τ } are
y This example is typical in that the kMC simulation j
h givesa macro-ormesoscopicquantity ofinterest(in this time-independent (i.e., that the transition rates {τj−1}
p from the current state to accessible states are time-
case the adatom diffusion coefficient), using information
[ independent). Thenthe probabilityp (t)dt thatthe par-
obtained from the microscopic or atomic scale. In fact, j
ticular element j will fail during the infinitesimal time
2 moleculardynamicscalculationsaloneareimpracticalfor
interval[t,t+dt]isgivenbytheprobabilitythatelement
v simulating adatom diffusion over a surface: the adatom
4 at room temperature attempts a jump every 10−13 sec- j will not fail prior to time t, multiplied by the proba-
6 bility τ−1dt that it will fail during the subsequent time
onds,whileitissuccessfulonceevery10−5 seconds(soin j
4 interval dt; that is,
MDsimulationsitisthenumerousunsuccessfulattempts
2
. at many sites, rather than the rare successful one, that t
2 dt
1 enables the energy landscape to be mapped and jump pj(t)dt= 1− pj(t′)dt′ . (1)
7 rates determined).  Z τj
0
0 The great virtue of the kMC method is that it en-  
: ables simulation of physical processes that involve very This equation simplifies to
v
Xi dtiimspearinattehteimpehyssciaclaels.syItstiesmthbene ceossrernectitalyl trheaptrothdeucaecdtubayl t 1
p (t)= 1− p (t′)dt′ . (2)
r the residence times in the states of the model system, j  j τ
a or equivalently, by the time between events in the model Z0 j
system.  
TakingthederivativeofeachsideofEq. (2)withrespect
While equivalent statements, the different semantics
totandintegratingproducestheexponentialprobability
suggests the ability of the kMC method to obtain prop-
density function (PDF)
erties of an equilibrium system as well as to model the
dynamic evolution of a non-equilibrium system. In the 1 t
p (t)= exp − . (3)
case of equilibrium systems, each state may be visited j τ τ
j (cid:20) j(cid:21)
many times, so that the residence times in the states
are proportional to the equilibrium populations of the Specifically, pj(t) is the distribution of failure times for
element j; the average value of the failure time is ht i=
j
∞
tp (t)dt=τ .
0 j j
Toperformasimulation,weneedtorandomlyselecta
∗Electronicaddress: [email protected] tRimetfromthisdistribution. Theformulaforconverting
2
a random number x taken from the uniform probability A different expression for △t may be found by noting
distributionP(x)=1(suchxvaluesareproducedbythe thattheprobabilityp(△t)dtthatthenextelementtofail
standardrandomnumber generators)to the correspond- willdosoduringtheinfinitesimaltimeinterval[△t,△t+
ing t value is derived as follows. The probabilities p(t)dt dt] is
and P(x)dx must be equal, so p(t)dt=dx. Thus
N N
t
x(t)= p(t′)dt′ =1−exp − t . (4) p(△t)dt= pj(△t)dt [1−pk(△t)] (9)
Z (cid:20) τj(cid:21) Xj=1 k(6=Yj)=1 
0
where the contentof the curly brackets is the probabil-
Inverting this expression then gives the desired relation
ity that element j will fail during the infinitesimal time
between x and t,
interval [△t,△t+dt], multiplied by the probability that
t=τj[−ln(1−x)] (5) no other element will fail during △t. Then
with x randomly chosen from the interval [0,1).
N N
Inthiswaywegettheset{t }ofelementfailuretimes. 1 △t △t
j p(△t)= exp − exp −
If the elements {j} are completely independent of one τ τ τ 
another, they fail according to the time ordering of the Xj=1 j (cid:20) j (cid:21)k(6=Yj)=1 (cid:20) k (cid:21)
setIf{,thj}ow. ever, the elements are not independent, the re- = N τ1 exp −△t Nk=1τk−1 
maining values {τj} are altered with each failure. Con- Xj=1(cid:26) j h P i(cid:27)
sider that the smallest member of the set {tj} is T1, so N
the first element to fail does so at time T . Which is = τ−1 exp −△t N τ−1 . (10)
1 j j=1 j
the next element to fail, and when does that occur? If j=1 ! h i
P P
element j is still functioning at time T , the PDF
1 As above,a value △t may be randomly chosenfrom this
1 (t−T ) distribution by use of the relation
1
p (t−T )= exp − (6)
j 1 τ′ τ′
j " j # −ln(1−x)
△t= (11)
where τ′ is the new lifetime ofelement j for time t>T . τ−1
j 1 j j
Equation (6) gives the distribution of times t > T at
1 P
which element j fails, or equivalently, the distribution of where x ∈ [0,1). The average value is given by Eq. (8),
time intervals △t ≡ t−T at the end of which element as expected.
1
j fails. As above, a value △t may be randomly chosen Theevolutionofthesystemofelementsisthusaccom-
from this distribution by use of the relation plished by selecting an element to fail, and incrementing
time accordingly. With each failure the system enters a
△t=τ′[−ln(1−x)] (7)
j new state; obviously it cannot return to any old states.
where x ∈ [0,1). Thus the next element to fail is that For this reason it has been convenient to use the set of
producing the smallest member of the set {△tj}. “lifetimes” {τj} corresponding to the surviving elements
Analternativeapproachtosimulatingthe systemevo- j,ratherthanthe setoftransitionrates{ki→j} fortran-
lution is to randomly choose the next element to fail sitions from the current state i to the accessible states j
(the system will enter state j if element j is the next to
according to the set of probabilities τ−1/ τ−1 ,
k j j fail). The connectionbetweenthe setsis madeby noting
where element k is one of the set {j},nτk−1/Pj (cid:0)τj−1(cid:1)ois that jτj−1 = jki→j, where the system is currently in
theprobabilitythatelementk isthenexttofail,andthe state i in either case.
sumisoverallN currentlysurvivingelementPsj ((cid:0)reme(cid:1)m- ByPmaking tPhis replacement in the denominators of
ber that the values {τj} may change after every failure). Eqs. (8) and (11), those equations for △t may be used
This failure occurs at the end of the time increment △t, forequilibrium systems(wherethesetofstatesandtran-
thatmaybetakentobetheaverage valueofthesmallest sition rates between states don’t change) as well. More
member of the set {△tj} (were the set to be calculated generally, τj in the equations above may be replaced by
innumerable times), k−1 when it is understood that the system is currently
i→j
in state i.
τ−1 1 N τ−1
△t= △t k = △t k For computational efficiency, it is most typical for a
* k jτj−1+ N k=1 k jτj−1 simulationtoselectthetransitionfromcurrentstateiac-
X
= 1 kP[−ln(1−xk)] = h−ln(P1−xk)i cording to the set of probabilities ki→j′/ jki→j and
N τ−1 τ−1 calculatethe transitiontime intervnal△t usPingeithoer Eq.
P j j j j
(8) or (11), rather than to calculate the set {△t } from
1 j
= .P P (8) the relation△t =k−1 [−ln(1−x)] and choose the des-
jτj−1 tinationstatejj′ fromi→thje smallestofthosevalues. When
P
3
thedistribution ofsimulationcompletiontimesisdesired where the lower limit x(T) = 1−exp[−(T/τ′)γ′] is ob-
(forexample,thetimeatwhichallelementshavefailed), tained from Eq. (15), and γ′ and τ′ are the parameter
it is necessary to perform a large number of nominally values for t>T. Then the counterpart to Eq. (4) is
identicalsimulations,using Eq. (11) ratherthan Eq. (8)
for every transition time interval △t. (Note that this T γ′ t γ′
x(t)=x(T)+exp − −exp −
distributionwill convergeto a Gaussiandistribution, ac- τ′ τ′
cording to the central limit theorem.) " (cid:18) (cid:19) # " (cid:18) (cid:19) #
More generally, evolution of a dynamic system is not t γ′
=1−exp − (18)
dsteasnctrsib.edIbnytahissetc{aksei→tjh}eoref tiismen-oindseetpeonfdepnrtobraatbeilcitoines- " (cid:18)τ′(cid:19) #
ki→j′/ jki→j from which to select the destination so the failure time for element j is
sbnteatcealjc′u.PlaHtoewdefvroeorm, aasPetD{Ft,j}givofintgratnhseitsieoqnuetinmceesocfaenvesntitlsl tj =τj′[−ln(1−x)]1/γj′ (19)
from an initial time. For example, consider that the dis-
where x∈[x(T),1).
tribution of failure times for element j is given by the
two-parameter (γ and τ) Weibull distribution
tγ−1 t γ III. DISCUSSION
p(t)=γ exp − (12)
τγ τ
(cid:20) (cid:18) (cid:19) (cid:21) In physics and materials research,the kMC method is
ratherthanbyEq. (3). Forγ >1thisPDFdecaysfaster mostoftenusedforsimulationofnon-deterministic,ther-
than exponentially and so might be appropriate for an mallyactivatedprocesses,asthesegivetime-independent
element that “ages”, or accumulates damage over time. rate constants (time-independent probabilities) for tran-
(Indeed, the probability π(t)dt that the element, having sitions between states. A similar, but more complex,
survivedto time t, will immediately fail is γ(tγ−1/τγ)dt, example to the adatom diffusion calculation mentioned
meaning that the probability of immediate failure in- in the Introduction is diffusion of vacancies and intersti-
creases with time.) The first and second moments of tial atoms in a grainboundary [1], where the transitions
this distribution are, respectively, between states may occur by very complex mechanisms
(for example, by concerted motion of atoms). However
1 1 the rates k are of the simple form
hti=τ Γ (13)
γ γ
(cid:18) (cid:19) E
m
k =νexp − (20)
i→j
and k T
(cid:20) B (cid:21)
2 2 where E is the energy barrier between the initial state
t2 =τ2 Γ (14) m
γ γ i and the final state j (which the system may overcome
(cid:18) (cid:19)
(cid:10) (cid:11) with thermal energy supplied by local temperature fluc-
where Γ is the Gamma function. The counterpartto Eq. tuations),k istheBoltzmannconstant,T istheaverage
B
(4) is temperature of the system, and ν is the frequency with
which the system attempts to make the transition. The
t
t γ latter quantity gives the time scale for the microscopic
x(t)= p(t′)dt′ =1−exp − , (15) process, that is needed for the kMC simulation.
τ
Z (cid:20) (cid:18) (cid:19) (cid:21) Besides its natural application to atomic-scale pro-
0
cesses, the kMC method has been used to simulate the
so the failure time for element j is time-dependent damage and failure of material under
stress. For example,Curtinet al. [2]considereda spring
t =τ [−ln(1−x)]1/γj (16)
j j network model where the rate r of failure at site j is
where x∈[0,1). r (t)=Aσ (t)η (21)
j j
If the parameter values {γ ,τ } change with each ele-
j j
mentfailure(reflecting,say,morerapidaging),thePDFs whereσ (t)isthelocalstressattimet,andtheexponent
j
{pj(t)} given by Eq. (12) change accordingly. To calcu- η > 1 ensures a nonlinear relationship between damage
late a new set {tj} of transition times following an ele- rateandstressasisthecaseforcreep. Thetimebetween
ment failure at time T, we note that x and t must be failures was calculatedby Eq. (8), after which all the lo-
related by calstresses{σ }andfailurerates{r }wererecalculated.
j j
In contrast, Harlow et al. [3] considered a polycrystal
x t
understress,andrandomlychose,fromaprobabilitydis-
dx′ = p(t′)dt′ (17) tributionthat included the localstress asa parameter,a
Z Z time-to-failure t for each grain boundary facet j. Then
x(T) T j
4
time was advanced by the amount of the smallest mem- tem has equilibrated in the basin (which is to say, there
ber of the set {t }, the load supported by that failed is no correlationbetween the entry and exit states).
j
facet was redistributed to adjacent facets, and a new set Toconclude,thekineticMonteCarlomethodisapow-
{tj} was obtained from new distributions. In a similar erful technique for simulating the dynamics or evolution
wayAndersenandSornette[4]consideredafiberrupture of non-deterministic systems. Thus it should have ap-
modelwherefailureofeachfiberaffectedthefailure-time plications beyond the physical sciences, for example to
probability distributions for the remaining fibers. biology, ecology, risk assessment, and finance.
Effective use of the kMC method obviously requires
thatallimportantstatesofthesystem,andallimportant
transitions between them, be identified. How to ensure
this, and other technical and computational(i.e., practi- Acknowledgments
cal)issuesareaddressedinareviewbyVoter[5]. Another
problemispresentedby“basins”ofstates,whicharesub-
This work was supported in part by the INL Labora-
sets of mutually accessible states from which escape is a
tory Directed Research and Development Program un-
very rare event. This of course defeats the purpose of
der DOE Idaho Operations Office Contract DE-AC07-
the kMC method. However Van Siclen [6] has recently
05ID14517.
shown how to calculate the residence time for a system
trapped in a basin, under the assumption that the sys-
[1] M. A. Sorensen, Y. Mishin, and A. F. Voter, “Diffusion colation andhierarchicalfiberbundles,”Phys.Rev.E72,
mechanisms in Cu grain boundaries,” Phys. Rev. B 62, 056124 (2005).
3658-3673 (2000). [5] A. F. Voter, “Introduction to the Kinetic Monte Carlo
[2] W.A.Curtin,M.Pamel,andH.Scher,“Time-dependent Method,” in Radiation Effects in Solids, edited by K. E.
damage evolution and failure in materials. II. Simula- Sickafus, E. A. Kotomin, and B. P. Uberuaga (Springer,
tions,” Phys. Rev.B 55, 12051-12061 (1997). NATO Publishing Unit, Dordrecht, The Netherlands,
[3] D.G.Harlow,H.-M.Lu,J.A.Hittinger,T.J.Delph,and 2005), in press.
R.P.Wei,“Athree-dimensionalmodelfortheprobabilis- [6] C.DeW.VanSiclen,“Stochasticmethodforaccommoda-
tic intergranular failure of polycrystalline arrays,” Mod- tion of equilibrating basins in kinetic Monte Carlo sim-
elling Simul. Mater. Sci. Eng. 4, 261-279 (1996). ulations,” J. Phys.: Condens. Matter 19, 072201 (2007).
[4] J.V.AndersenandD.Sornette,“Predictingfailure using
conditioning on damage history: Demonstration on per-