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].
This preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentMajor Issues:
There are some issues that should have been dealt with the two major ones being how to round the
fractional molecules per droplet and the correct way to compute the rates.
A fractional molecule should
have rounded in a probabilisticallycorrect way depending on the fraction.
If there were 0.25 molecules
per droplet, you should have 0 ROOH’s 75% of the time and 1 ROOH 25% of the time.
A more correct
way to deal with this (and droplets that had > 1 molecules) would be to use a probability distribution
(e.g. Poisson distribution) for the number of molecules per droplet.
This is the end of the preview.
Sign up
to
access the rest of the document.
 Fall '04
 ClarkColton
 Standard Deviation, Mean, ROOH

Click to edit the document details