Transcript maws3 5886

Bridging scales:
Ab initio molecular dynamics and
kinetic Monte Carlo simulations
Karsten Reuter
Fritz-Haber-Institut, Berlin
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
I. Straightforward time-evolution:
Ab initio molecular dynamics
Understanding Molecular Simulation,
D. Frenkel and B. Smit, Academic Press (2002)
MD: Numerical integration of Newton’s equation

..



mi ri   ri V (r1 , r2 , , rN )
e.g. Verlet algorithm



1  
r (t  t )  2r (t )  r (t  t )  [( )V (r (t ))]t 2  error (t 4 )
m
Solid-state:
Δt ~ 10-15 s
Alternative algorithms:
- velocity Verlet
- leapfrog
- predictor-corrector
- Runge-Kutta
…
Keeping everything: expensive, limited time scale
Example:
O2 dissociation at Al(111)
Total time of trajectory: 0.5 ps
Time step: 2.5 fs (200 steps)
CPU cost: 45 days on
1 Compaq
ES45 processor
J. Behler et al., Phys. Rev. Lett. 94, 036104 (2005)
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
First-principles kinetic Monte Carlo simulations
rB→A
ΔEA rA→B ΔEB
equilibrium
Monte Carlo
TS
B
EB
EA
A
<N>
kinetic
Monte Carlo
N
A
B
B
Molecular
Dynamics
A
t
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
The crucial ingredients to a kMC simulation

dP(r , t )
 


 
   (r  r ' ) P(r , t )    (r '  r ) P(r ' , t )


dt
r'
r'

P(r , t )
i) Elementary processes
 
 (r  r 'O
) cus
CObr
x
x
ii) Process rates
PES from density-functional theory
Transition state theory
Fixed process list vs. „on-the-fly“ kMC
Lattice vs. off-lattice kMC
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  (i)
and find process “k”:
k
 
i=1
(i)
k
r1 R
k-1
 r1 R    (i)
i =1
0
update clock
t  t – ln(r2)/R
END
Execute process number “k”,
i.e. update configuration
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)
Rare event dynamics
TS
Brute force approach to rate constants:
i)
ii)
E
IS
FS
Have accurate potential energy surface (forces)
Run MD trajectory so long, that it establishes
equilibrium, crossing the barrier many, many
times back and forth:
 k=
no. of crossings IS  FS per unit time
fraction of time system has spent in IS
Yet:
- Relevant time step in MD run is fs (vibrations)
- Typical barrier E for surface reactions ~ 1 eV  10-2 - 102 reactions per second (TOF!)
- Requires to run trajectory over about 1015 – 1020 time steps
 unfeasable…
…and essentially 99,9999% of the time, the system will
just vibrate around IS basin (short time dynamics)
 require approximate theories to obtain process rates!
Transition state (activated complex) theory I
( Eyring, Evans, Polanyi, ~1935 )
Assumptions:
i)
Reaction system passes the barrier only once (no recrossings)
ii)
Energy distribution of reactant DOF is Boltzmannlike (many collisions compared to reaction events
yield equilibrium between activated complex and IS,
except with respect to the reaction coordinate)
iii)
iv)
X
X
eq.
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)
IS
FS
Passage over barrier is a classical event (no tunneling)
Derivation: see e.g.
K.J. Laidler, Chemical kinetics,
Harper & Row, New York (1987)

kTST
ISFS
=
[(
kT
h
)e
S/k
] e-E/kT
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)
- Attempt frequency/preexponential factor koTST =
kT S/k
e
h
~ 1013 sec-1
~1
- In harmonic TST, koTST is given by the harmonic normal modes at the IS and TS
(as explained in talk by Christian Ratsch on Tuesday)
 problem reduces to locating transition states,
i.e. saddle points in high dimensional PES
Transition state search algorithms I: grid and drag
i) Grid method:
- Compute PES on a (regular) grid
- Scales like: (no. of points) dim
- e.g. 59 ~ 2 million grid points
FS
x
x
x
 often unfeasable
xx
ii) Drag method:
x
- Choose appropriate reaction coordinate q
- Constrain q and relax all other DOF
- Move from IS to FS
 highly dependent on good
reaction coordinate
 hysteresis!
x
IS
x
x
x
xx
x
TS
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 (downhillstep)
- Find max-energy point Ri along connecting
line between two replicas
FS
x
Ro
- 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)
IS
x
x x
x x R
1
x x
x
x R
x
x
2
x
x
x
x x
TS
I.V. Ionova and E.A. Carter,
J. Chem. Phys. 98, 6377 (1993)
Transition state search algorithms III: nudged elastic band
- Initialize with several images {Ri} along a
straight-line interpolation
- Minimize
S(R1, …, RN) =
FS
x
i E(Ri) + i k/2 (Ri+1 - 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)
x
x
x
x
x
IS
x
 widely applied workhorse
 has problems, if energy varies largely along path,
but very little perpendicular to it (kinky PES)
x F ||
spring
x

Ftrue
xx
TS
G. Mills and H. Jónsson,
Phys. Rev. Lett. 72, 1124 (1994)
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 perpendicular to dimer axis
FS
x
 works without any information about FS
FR
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
IS
x
x
TS
G. Henkelman and H. Jónsson,
J. Chem. Phys. 111, 7010 (1999)
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)
Diffusion at metal surfaces: surprises…
Hopping mechanism
Ag(100) E = 0.45 eV
Au(100) E = 0.83 eV
Exchange mechanism
Ag(100) E = 0.73 eV
Au(100) E = 0.65 eV
B.D. Yu and M. Scheffler, Phys. Rev. B 56, R15569 (1997)
Automatized process identification
Accelerated molecular dynamics:
TAD
Hyperdynamics
Other approaches:
- metadynamics
- dimer method
…
V. If it works, it works:
First-principles kMC simulations for oxidation catalysis
A materials gap resolved: CO oxidation at Ru(0001) vs. RuO2(110)
K. Reuter et al., Chem. Phys. Lett. 352, 311 (2002)
cus site
H. Over and M. Muhler,
Prog. Surf. Sci. 72, 3 (2004)
bridge site
O
Ru
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
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)
A ( pCO , pO2 )-map of catalytic activity
600 K
pO2 (atm)
10-15 10-10 10-5
10+5
br/CO
br/CO
cus cus
CO
CO
CObr/-
Obr/
-
1
10+5
105
1
10-5
10-10 10-5
pCO (atm)
pCO (atm)
105
1
pO2 (atm)
br/O
cus
cus
O/O
Obr
1
10-5
log(TOF)
4.0
3.0
2.0
1.0
0.0
-1.0
-2.0
-3.0
K. Reuter, D. Frenkel and M. Scheffler, Phys. Rev. Lett. 93, 116105 (2004)
…and how about experiment?
pO2 (atm)
350 K
10-30 10-20 10-10
pCO (atm)
1
1
TT==350
350KK
-10 atm
ppCO
= 10
10-10
atm
O2 =
CObr/COcus
10-10
CObr/-/-
Obr/Ocus
10-20
log(TOF)
-4.0
-5.0
-6.0
-7.0
-8.0
-9.0
-10.0
-11.0
J. Wang, C.Y. Fan, K. Jacobi, and G. Ertl,
J. Phys. Chem. B 106, 3422 (2002)
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