8/7/2019 DMC presentation 08_05_03
1/42
11
Dynamic Monte CarloDynamic Monte CarloMethods:
Theory and Applications
ByByDaniela S. MainardiDaniela S. Mainardi
8/7/2019 DMC presentation 08_05_03
2/42
22
OutlineOutline
1. Theoretical IntroductionDMC MethodProgram Carlos 4.0
2. Applications
Temperature Programmed SimulationsLateral InteractionsVoltammetric Scans
8/7/2019 DMC presentation 08_05_03
3/42
33
THEORYTHEORY
8/7/2019 DMC presentation 08_05_03
4/42
44
Introduction to DMC MethodsIntroduction to DMC Methods
Quantum Mechanical(Ab Initio, DFT)
Classical(Monte Carlo, Molecular Dynamics)
Laws of Quantum Mechanics E H =
Quantum effects determine all the chemistry
1. Do not explicitly include the presence of electrons2. Effective force fields represent the electronic effects3. Physical and thermodynamic properties obtained using
classical statistical mechanics and Newtons Law (MD)
Hardware and software not enough to fully trerelatively large systems!!!
8/7/2019 DMC presentation 08_05_03
5/42
55
Introduction to DMC MethodsIntroduction to DMC Methods
MD simulations provide dynamic phenomena tthe order of nanoseconds.
For phenomena taking place in the range ofseconds,
Coarse-Grained methods
Dynamic Monte Carlo
8/7/2019 DMC presentation 08_05_03
6/42
66
DMC MethodsDMC Methods
Advantages:1.treatment ofnon-equilibrium systems in which
conditions change over time
Importance reaction mechanisms forheterogeneous catalysis 2.It is able tosimultaneously simulate many differe
transitions or reaction types covering events tspan time scales up to theorders of seconds 3.Thereal-time behavior of a system is simulate
8/7/2019 DMC presentation 08_05_03
7/42
77
DMC MethodsDMC Methods
Main Challenge:Reaction rates are needed as input
experiments or ab initio calculations
8/7/2019 DMC presentation 08_05_03
8/42
88
How can we implement theHow can we implement theDMC method?DMC method?
Carlos 4.0Carlos 4.0Copyright Johan J. LukkienCopyright Johan J. LukkienOctober 1995October 1995
If we want to study chemicalreactions on crystal surfaces
8/7/2019 DMC presentation 08_05_03
9/42
99
What is Carlos useful for?What is Carlos useful for?
Systems not at equilibrium like 1. Temperature Programmed systems 2. Voltammetricscans3. Oscillating systems, spatio-temporal patternformation4. Epitaxial growth
Non-homogeneous systems 1. Surface steps/defects
2. Lateral interactionsSystems with slow relaxation
8/7/2019 DMC presentation 08_05_03
10/42
1010
What do we need?What do we need?
1. Specify the chemical System we want tstudy
2. A stochastic model is built:* crystal surface,* the adsorbates,* a law to describe the microscopicreaction steps that change the surfaceconfiguration over time.
The behavior of this system over time isdetermined by the rates of these reactions(specified as probabilities).
8/7/2019 DMC presentation 08_05_03
11/42
1111
Thesurfacesurface is represented by a lattice: eachlattice point corresponds to a surface site
Elements in the model
The crystal surface
fcc (111) surfaceTopBridgeHollow
Atom
Unit cell
8/7/2019 DMC presentation 08_05_03
12/42
8/7/2019 DMC presentation 08_05_03
13/42
1313
Elements in the model
The evolution of the system over time is described by the chemicalmaster equationmaster equationThe law
( ) = cc cccc t cP k t cP k dt t cdP
''' ),(),'(
),(
P(c,t) probability of finding the system inconfigurationc at time t .k cc microscopic rate of the reaction that transfer cinto c
=T k
E r
B
act exp Arrhenius LawArrhenius Law
ratereactionr k
r k
dt
d
A
cc
AAAA
k A
~
*
'
==
8/7/2019 DMC presentation 08_05_03
14/42
1414
We could either give Carlos the macroscopic rates orEact and , and Carlos will compute the microscopicones to solve the Master equation
There are methods which are used internally to solvethis equation: First-Reaction Method (FRM) First-Reaction Method (FRM) Variable Step Size Method (VSSM) Variable Step Size Method (VSSM)
Random Selection Method (RSM) Random Selection Method (RSM)
How Carlos works?
8/7/2019 DMC presentation 08_05_03
15/42
1515
= t
t
t r dt )(exp ''
r r = the time dependent rate of a reaction= the time dependent rate of a reaction = a random number selected uniformly (0,1= a random number selected uniformly (0,1
Solving the Master equation
The way in which this equation issolved gives rise to the differentmethods
R.J. Gelten, R.A. VanSanten, A.P.J. Jansen, in: P.B. Balbuena,J.M. Seminario (Eds.), Molecular dynamics: From classical toquantum methods, Elsevier Science, Amsterdam, 1999, p. 737
8/7/2019 DMC presentation 08_05_03
16/42
1616
FRMFRM
It can be used forBOTHBOTH time-dependent or timeindependent rate constants.
According to this method, when the system igiven configuration , the set of all possiblereactions is determined, and a time of occurrt is generated forEACHEACH reaction
Carloslists all the reactions from the one withthe smallestt to the one with the largest tim
8/7/2019 DMC presentation 08_05_03
17/42
1717
- Then, the reaction with the- Then, the reaction with theSMALLEST SMALLEST t t is selected, the configuration is changeis selected, the configuration is changeaccordingly, and the timeaccordingly, and the timet t is incrementedis incremented
inint t - Finally, the set of possible reactions is- Finally, the set of possible reactions is
generated according to the newgenerated according to the newconfigurationconfiguration , and the procedure is, and the procedure isrepeatedrepeated
FRMFRM
8/7/2019 DMC presentation 08_05_03
18/42
1818
VSSMVSSM
It can be used for time-independent rateconstantsONLY
According to this method, when the system given configuration , the set of all possiblereactions is determined, andONLY ONEONLY ONE time ofoccurrence t is generated forALLALLreactions
Then,ONEONE reaction is selected with a probabiproportional to its rate, the configuration ischanged accordingly, the procedure is repea
8/7/2019 DMC presentation 08_05_03
19/42
1919
It can be used forBOTHBOTH time-dependent or time-independent rate constants.
ASITESITE is selected with a probability 1/N, beingtotal number of sites. Then, a given reaction isselected with probability proportional to its rate
When this reaction is possible on the site, it isexecuted. After each selection of a site, time isincremented by t which is selectedRAMDOMLYRAMDOMLYfroma FIXEDFIXED distribution
the configuration is changed accordingly, the
i
ir Nt exp1
RSM
8/7/2019 DMC presentation 08_05_03
20/42
2020
FRM and VSSM can be implementedvery efficientlyvery efficiently because information generated at some time step cabe re-used at subsequent steps.
The advantage ofVSSM and RSM over FRMis that thecost per generated transition (reaction) isindependentindependentof the sizeof the sizeof the grid used in the simulation, whereain FRM it scales as thelogarithmlogarithmof the number of gridpoints. Due to this, VSSM and RSM are more efficin memory use and execution time.
The advantage ofFRM over VSSM is that it is asuitable method to simulate reactions withtime-time-dependentdependenttransition probabilities.
The advantage ofRSM is that is easy to implement andfor some models (or reaction types) it is a veryfastfastmethod.
Comments
8/7/2019 DMC presentation 08_05_03
21/42
2121
Preparing Input Files...Preparing Input Files...
CARLOS needs two input files to run aCARLOS needs two input files to run asimulation:simulation:1.1. Initial Configuration:Initial Configuration: can be made by CARLOS,
.con file, or can be constructed directly using aeditorname.init name.init
2.2. System FileSystem File:: contains a complete description of simulation model
name.sim name.sim
8/7/2019 DMC presentation 08_05_03
22/42
2222
Geting outputs...Geting outputs...The logfileThe logfile:: namename(information recorded during the simulation)The parameter fileThe parameter file:: name.parname.par
(describing both the simulation model and relevant infoabout the simulation such as timing, memory and methoThe graph fileThe graph file:: name.xxxx.gifname.xxxx.gif(pictures of the lattice in gif format, if g option is not g The configuration fileThe configuration file:: name.conname.con(the final configuration if the -C option is given)
The system fileThe system file:: name.sysname.sys(a copy of the system file if the -S option is present )
8/7/2019 DMC presentation 08_05_03
23/42
2323
APPLICATIONSAPPLICATIONS
8/7/2019 DMC presentation 08_05_03
24/42
2424
OO22 desorption from Pt (111)desorption from Pt (111)
Temperature Programmed
simulations
8/7/2019 DMC presentation 08_05_03
25/42
2525
ONLY O2 desorption reaction isinvestigated:
Prefactor: 1.0 x 1013s-1Activation energy: 0.36 eV
1 O 2 ML
O2 desorption O 2(ads_b) O 2 + b b = empty bridge
O2 O
2O
2O
2O
2O2
O2
Initial configuration: full, i.e, 1 O2 ML is adsorbed on the surfa
T = 100.0 + 2.0 * tI increased T from 100K to 190 K:
8/7/2019 DMC presentation 08_05_03
26/42
2626
Final configuration
O2_des
0
0.2
0.4
0.6
100 120 140 160 180
T (K)
Reactio
ns/unit cell/sec
0
0.2
0.4
0.6
0.8
1
100 120 140 160 180
T (K)
O2 Concentration
Initial configuration
A. Winkler et al.,Surface Science 201(1988), 419 - 443
RESULTS
8/7/2019 DMC presentation 08_05_03
27/42
2727
OO22 adsorption, dissociation,adsorption, dissociation,desorption, and O diffusiondesorption, and O diffusionon Pt (111)on Pt (111)
Temperature Programmed
simulations
8/7/2019 DMC presentation 08_05_03
28/42
2828
TopBridgeHollowPt atom
O2
O2
O2 adsorption on Pt(111) occurs on bridge sites (b). Upon dissociation of O 2 , O atoms occupy hollow sites
(h).
Pt(111) surface
O2 adsorption O 2 + b O 2(ads_b)
O2 dissociation O 2(ads_b) + 2h 2O(ads_h) + b
O diffusion O (ads_h hcp ) + h fcc h hcp + O (ads_h fcc )
O2 desorption O 2(ads_b) O 2 + b
64 x 64 lattice
b = empty bridgeh = empty hollow
Initial configuration: empty Pt surface
8/7/2019 DMC presentation 08_05_03
29/42
2929
Reaction Eact(eV)
Eact/KB (K) (s -1 )
O2 ads 0.09 1044.4 F(T)
O2 diss 0.3 3481.33 1.20E+12O2 des 0.36 4177.6 1.00E+13O diff 0.13 1508.58 2.60E+12
Prefactor calculated using Hertz-Knudsen formula:
MkT P F
2=KTE
0site
a
eFSA += [P] = Pa[k] = J/K[M] = kg (= 0.032 kg/NA)[T] = K[F] = molecules/m2.sec
8/7/2019 DMC presentation 08_05_03
30/42
3030
64 x 64 Pt atoms
In a lattice with 64x64 Pt atoms we hav64x64 unit cells (due to the periodic
boundary conditions)There are 3 bridge, 2 hollow, and 1 topsites per unit cell, therefore,
assuming that all sites have the same ar
N sites per unit cell = 6Total N of sites = 6 x 64 x 64 = 24576
24576 sites have an area of 3.151 x 10-12 cm21 site 1.282 x 10-16 cm2
For my simulations I will consider:
Asite = Area of 1 adsorption site = 1.28 x 10-16 cm2 = 1.28 x 10-20 m2
Assuming a flat surface:
......
.
.
.
.
.
.
rPt = 1.387 = 1.387 x 10-8 cm2DPt = 2.774 x 10-8 cm2
W
W = 64 x DPt = 1.775 x 10-6 cm2 Area = W2 = 3.151 x 10-12 cm2
8/7/2019 DMC presentation 08_05_03
31/42
3131
P = 0.001 PaA = 1.28E-20m2
KTE
0site
a
eFSA +=
T S086.3801091 2.97E-01
92.3057527 2.95E-0196.4545828 2.85E-0198.3604058 2.71E-01
100 272093 2 55E 01
T range (K) O2ads ( s-1 )
100 - 108 4.80 E+03108 - 111 5.72 E+02111 - 120 2.46 E+02
120 - 135 7.88 E+01135 - 148 1.94 E+01148 - 154 7.42 E+00154 - 162 2.47 E+00162 - 170 1.18 E+00
values used in Carlos input
A. Winkler et al.,Surface Science 201(1988), 419 - 443
MkT P F
2=
8/7/2019 DMC presentation 08_05_03
32/42
3232
0.00E+00
1.00E+02
2.00E+02
3.00E+02
120 130 140 150 160 170
T
Prefactor (O2 ads)
0.00E+00
1.00E+03
2.00E+03
3.00E+03
4.00E+03
5.00E+03
6.00E+03
100 105 110 115 120
T
Prefac
tor (O2 ads)
KTE
0site
a
eFSA +
= Prefactor for O2 ads calculated usingPrefactor used in Carlos at each T range
8/7/2019 DMC presentation 08_05_03
33/42
3333
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
100 120 140 160T (K)
Reactions/unit cell/sec
O2_desO2_adsO2_dissO_diff
00.1
0.20.30.4
0.5
0.60.70.8
0.91
100 120 140 160
T (K)
Concentration
O2O
P=10-3 Pa
Results with the previousapproximation
Heating rate = 2K/s
8/7/2019 DMC presentation 08_05_03
34/42
3434
Using Discrete StepsUsing Discrete Steps
(DS)...(DS)...and the following fitting for S0(T) given by Winkler et al
0
0.02
0.040.06
0.08
0.1
0.12
0.14
0.16
100 150 200 250 300
T (K
Rate O2_ads (1/s)
W [O2ads]mathematical fitti
W [O2ads] (1/sec)using S0 Exp.
100 K < T < 111 K
T xS 20 1049.17526.1
=
= T
xS 100
exp1012.3 20
T > 111K 00.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
100 150 200 250 300
Rate O2_ads (1/s)
ReactRate given bCarlos
W [O2ads]mathematical fitting
8/7/2019 DMC presentation 08_05_03
35/42
3535
Using DS I gave Carlos the functional formUsing DS I gave Carlos the functional formof the rate for Oof the rate for O 22 ads as Sads as S 00FAFAsitesite
P = 10-3 Pa
0.3
0.4
0.5
/unit cell/sec O2_des
O2_adsO2_dissO_diff
Heating rate = 2K/s
00.1
0.20.30.40.5
0.6
0.70.8
0.91
100 120 140 160 180 200
T (K)
Concentration
OO2
00.1
0.20.30.40.5
0.60.70.8
0.91
100 120 140 160
T (K)
Concentration
O2O
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.160.18
0.2
100 120 140 160T (K)
Reactions/unit cell/sec O2_des
O2_adsO2_dissO_diff
8/7/2019 DMC presentation 08_05_03
36/42
3636
OO22 electroreduction on Ptelectroreduction on Pt(111)(111)
Voltammetric scans
8/7/2019 DMC presentation 08_05_03
37/42
3737
I O 2(gas) + e - + H + + b HO 2(ads_b)II HO 2(ads_b) + 2h + e
- + H + 2OH(ads_h) + b
III O 2(gas) + b + e- O2
-(ads_b)
IV O 2-(ads_b) + 2h O -(ads_h) + O(ads_h) + b
V O(ads_h hcp ) + h fcc hhcp + O(ads_h fcc )
VI O - (ads_h) + H + OH (ads_h)
VII O (ads_h) + e - + H + OH (ads_h)
VIII OH (ads_h) + t h + OH (ads_t)IX OH (ads_t) + H + + e - H2O(ads_t)
X OH(ads_h) + H + + e - H2O(liq) + h
XI H 2O(ads_t) H2O (liq) + t
Proposed Mechanism
8/7/2019 DMC presentation 08_05_03
38/42
3838
Lateral Interactions
HO2O2- H2O
Up to 3rd neighbors repulsions for O
1st neighbors repulsionsModel 1 Model 2
Model 3
8/7/2019 DMC presentation 08_05_03
39/42
3939
0.8
1
tion O2
OH
Results
Model 1Model 2
0.8
1
tion
8/7/2019 DMC presentation 08_05_03
40/42
4040
( )
+= T K
V eE r
B
a
ii
i 00
exp
e 0 = charge of an electron Eai = activation energy at V=0V = overpotential i0 = pre-exponential factor = transfer coefficient = 0.5 if the reaction consumelectrons = 0.5 if it produces e-
V is taken with respect tthe standard reversiblepotential for O
2
reduction, 1.23 V, on thehydrogen scale.
Voltammetric ScansIn the presence of an overpotentialV (defined as thedifference between the actual and the equilibriumpotential), the rate of a reactionr i is given by:
Current = (e- production rate) (e- consumption rate)The overpotential was decreased linearly from 0 until -1.12 V a rate of 50 mV/sec starting from a clean Pt(111) surface.
8/7/2019 DMC presentation 08_05_03
41/42
4141
0-1
1
1.1
HE)
Tafel plot
Model V0(V) Tafel Slope(mA/cm2)1 -0.36 53 182 -0.39 34 12
3 -0.39 41 15
8/7/2019 DMC presentation 08_05_03
42/42
4242
AcknowledmentsAcknowledments
Dr. Perla BalbuenaDr. Perla BalbuenaMolecular modeling group:Molecular modeling group: Dr. WangDr. Wang Dr. Martinez-limiaDr. Martinez-limia Yingchun ZhangYingchun Zhang Sergio CalvoSergio Calvo Zhihui GuZhihui Gu Francisco Tarazona VasquezFrancisco Tarazona Vasquez Eduardo LamasEduardo Lamas Tyler WattTyler Watt Kyle CorbinKyle Corbin
Charlotte CooperCharlotte Cooper Diego AltomareDiego Altomare
Top Related