10.34 Fall 2006
HW Set #11 Solution
Kinetic Monte Carlo
Notes:
The code
p1_main
will not generate all of the output for the problem, but enough for parts 1 and 2.
The remaining part can be done with
p1_main_part3
. However, once these are run (or if you
download the *.mat files posted), you can simply use the
p1_postprocessing
file to generate all
the plots and output. If you change the name of the saved *.mat files, you will need to also change the
“load” commands.
Algorithm (Gillespie):
The way to perform kinetic MC is relatively simple in theory.
The idea is that you calculate a
characteristic time for a reaction to occur based on the rates of all possible reaction.
This characteristic
time is defined as:
−
1
τ
=
⎜
⎛
∑
rate
i
⎟
⎞
units of time
(
per reaction event
)
⎝
reactions i
⎠
It is of the utmost importance to have the correct units on tau, otherwise the results will be meaningless.
Tau turns out to be an extensive property and will ultimately depend on the system size, as well as the
other parameters that determine the rate of a reaction (T, C
i
’s, parameters).
You can rationalize this
because if you have a larger volume (more molecules), the time between reaction events should be less.
This is why you never want volume terms in you rates when computing tau.
The
units
on the rates
should be in terms of
reactions events per time
for the entire system
.
Once tau is calculated, you estimate the time until the next reaction event by using a random number:
t
rxn
=
t
current
−
⋅
ln
(
rand
[
0,1
])
Now, you must determine which reaction will occur when
t
rxn
is reached.
At this time, you have
determined that some reaction will occur, and the probability of which reaction it is will be determined
by the relative rates.
A reaction with a fast rate will have a higher probability of occurring than one with
a low rate. This is typically done by using another random number between 0 and 1, and giving each
reaction a chunk of the interval [0,1] based on the value of
rate
i
⋅
. Taking all reactions will span the
entire interval, and reactions with larger rate will happen more often.
After this is done, it is simply a matter of keeping track of all of the appropriate variables and optimizing
your code to run in an efficient manner, which is not necessarily a trivial task.
Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to
Chemical Engineering, Fall 2006. MIT OpenCourseWare (http://ocw.mit.edu),
Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].