Bridging scales: Ab initio molecular dynamics and kinetic Monte Carlo simulations

1
Bridging scales Ab initio molecular dynamics
and kinetic Monte Carlo simulations
Karsten Reuter Fritz-Haber-Institut, Berlin
2
Multiscale modeling
Ab initio atomistic thermodynamics and
statistical mechanics of surface properties and
functions K. Reuter, C. Stampfl and M. Scheffler,
in Handbook of Materials Modeling, Part A.
Methods, (Ed.) Sidney Yip, Springer (Berlin,
2005). http//www.fhi-berlin.mpg.de/th/paper.html
3
I. Straightforward time-evolution Ab initio
molecular dynamics
Understanding Molecular Simulation, D. Frenkel
and B. Smit, Academic Press (2002)
4
MD Numerical integration of Newtons equation
e.g. Verlet algorithm
Solid-state ?t 10-15 s
Alternative algorithms - velocity Verlet -
leapfrog - predictor-corrector -
Runge-Kutta
5
Keeping everything expensive, limited time scale
6
II. First-principles kinetic Monte Carlo
simulations rare events and long time scales
IPAM tutorial by Kristen Fichthorn http//www.ipam
.ucla.edu/publications/matut/matut_5907.ppt
7
First-principles kinetic Monte Carlo simulations
TS
B
A
8
Kinetic Monte Carlo essentially coarse-grained MD
Molecular Dynamics the whole trajectory
Kinetic Monte Carlo coarse-grained hops
ab initio MD up to 50 ps
ab initio kMC up to minutes
9
The crucial ingredients to a kMC simulation
i) Elementary processes Fixed process list vs.
on-the-fly kMC Lattice vs. off-lattice kMC
ii) Process rates PES from density-functional
theory Transition state theory
10
Flowchart of a kinetic Monte-Carlo simulation
START
Get two random numbers r1 , r2 ? 0,1
1
determine all possible processes i for given
configuration of your system and build a
list. Get all rates ? (i)
Calculate R ?i G (i) and find process k
k k-1 ? G
(i) ? r1 R ? ? G (i) i1
i 1
r1 R
k
0
update clock t ? t ln(r2)/R
Execute process number k, i.e. update
configuration
END
11
III. Getting the rates
Chemical kinetics, J.K. Laidler, Harper Row 3rd
ed. (New York, 1987) Methods for finding saddle
points and minimum energy paths, G. Henkelman, G.
Jóhannesson, and H. Jónsson, in Progress on
theoretical chemistry physics, S.D. Schwartz
(Ed.), Kluwer (Amsterdam, 2000)
12
Rare event dynamics
TS
?E
IS
FS
13
Transition state (activated complex) theory I
( Eyring, Evans, Polanyi, 1935 )
Assumptions i) Reaction system passes the
barrier only once (no re- crossings) ii) Energy
distribution of reactant DOF is Boltzmann- like
(many collisions compared to reaction
events yield equilibrium between activated
complex and IS, except with respect to the
reaction coordinate) iii) Passage over barrier
is the motion of only one DOF, the reaction
coordinate, which is independent of all other
motions of the activated complex (no
concerted motions) iv) Passage over barrier is
a classical event (no tunneling)
Derivation see e.g. K.J. Laidler, Chemical
kinetics, Harper Row, New York (1987)
14
Transition state (activated complex) theory II
- If assumptions i)-iv) are
fulfilled, kTST is exact. In general, kTST is an
upper limit to the real rate k fdyn
kTST In principle, one can compute so-called
dynamical corrections. In contrast to liquid
gas phase, fdyn 1 for solid-state processes
(? TST is a rather good approximation)

• problem reduces to locating transition states,
• i.e. saddle points in high dimensional PES

15
Transition state search algorithms I grid and
drag
16
Transition state search algorithms II ridge

• - Initialize with straight line interpolation
• and choose max-energy point Ro
• - Create two replicas slightly displaced from
• Ro on either side of the ridge (side-step)
• - Displace replicas along gradient (downhill-
• step)
• - Find max-energy point Ri along connecting
• line between two replicas
• - Sequentially decrease displacements in
• downhill- and side-steps when approaching
• TS
• ? works nicely on well-defined ridges
• ? difficult to optimize the displacements
• for a given system
• ? then poor performance (many force
• evaluations required)

I.V. Ionova and E.A. Carter, J. Chem. Phys. 98,
6377 (1993)
17
Transition state search algorithms III nudged
elastic band

• - Initialize with several images Ri along a
• straight-line interpolation
• - Minimize
• S(R1, , RN) ?i E(Ri) ?i k/2 (Ri1 - Ri )2
• - Problem
• - elastic band cuts corners
• - images tend to slide down towards
• low-energy IS/FS regions, leaving few
• images for relevant TS region
• - Solution
• - only spring force component parallel
• to path (no corner cutting)
• - only true force component perpendicular
• to path (no down-sliding)
• ? widely applied workhorse

x
x
x
x
x
x
x
x
G. Mills and H. Jónsson, Phys. Rev. Lett. 72,
1124 (1994)
18
Transition state search algorithms IV dimer
- Initialize by putting dimer at an extremum of
a high temperature MD-run - Rotate dimer to
minimize energy (? direction of lowest
frequency normal mode) - Move dimer along
projected gradient perpen- dicular to dimer
axis ? works without any information about
FS
F R

• Generally
• - performance scaling with DOF
• not really known
• not good for rough PES
• high CPU cost
• - no algorithm is fool-proof
• (still lots of room for new ideas)

F ?
G. Henkelman and H. Jónsson, J. Chem. Phys. 111,
7010 (1999)
19
IV. Identifying the processes
Extending the time scale in atomistic simulation
of materials, A.F. Voter, F. Montalenti and T.C.
Germann, Annu. Rev. Mater. Res. 32, 321 (2002)
20
Diffusion at metal surfaces surprises
Hopping mechanism Ag(100) ?E 0.45
eV Au(100) ?E 0.83 eV
B.D. Yu and M. Scheffler, Phys. Rev. B 56, R15569
(1997)
21
Automatized process identification
Accelerated molecular dynamics
Other approaches - metadynamics - dimer
method
22

• If it works, it works
• First-principles kMC simulations for oxidation
  catalysis

23
A materials gap resolved CO oxidation at
Ru(0001) vs. RuO2(110)
K. Reuter et al., Chem. Phys. Lett. 352, 311
(2002)
H. Over and M. Muhler, Prog. Surf. Sci. 72, 3
(2004)
24
kMC events for CO oxidation over RuO2(110)
Adsorption CO - unimolecular, O2
dissociative no barrier rate given by
impingement r So p/(2?mkT) Desorption CO
1st order, O2 2nd order out of DFT
adsorption well ( barrier) prefactor from
detailed balance Diffusion hops to nearest
neighbor sites site and element
specific barrier from DFT (TST) prefactor
1012 s-1 (generic) Reaction site
specific immediate desorption, no
readsorption barrier from DFT (TST) prefactor
from detailed balance
26 elementary processes considered
25
The steady-state of heterogeneous catalysis
T 600 K, pO2 1 atm, pCO 20 atm
Explicit information about fluctuations,
correlations and spatial distribution of
chemicals at the catalyst surface
K. Reuter and M. Scheffler, Phys. Rev. B
(submitted)
26
A ( pCO , pO2 )-map of catalytic activity
600 K
1
10-5
105
10-15
10-10
10-10
1
10-5
105
CObr/COcus
4.0
CObr/COcus
3.0
2.0
1.0
(CO/O)br/ (CO/O)cus
0.0
CObr/-
-1.0
-2.0
-3.0
Obr/Ocus
Obr/ -
Obr/Ocus
log(TOF)
K. Reuter, D. Frenkel and M. Scheffler, Phys.
Rev. Lett. 93, 116105 (2004)
27
and how about experiment?
350 K
J. Wang, C.Y. Fan, K. Jacobi, and G. Ertl, J.
Phys. Chem. B 106, 3422 (2002)
28
Ab initio MD and kMC simulations
Ab initio molecular dynamics - fully dynamics
of the system - straightforward, easy to
implement - times scales up to ns -
acceleration techniques under development
Ab initio kinetic Monte Carlo simulations -
coarse-grained time evolution (rare event
dynamics) - efficient treatment of statistical
interplay of a larger number of elementary
processes - time scales given by process rates,
often seconds or longer - process list
(process identification, lattice models) -
accuracy of rates (DFT-TST), high CPU cost -
low speed-up, if very fast processes present