1 / 28

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

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,

jillian-lara
Download Presentation

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

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  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 Newton’s equation e.g. Verlet algorithm Solid-state: Δt ~ 10-15 s Alternative algorithms: - velocity Verlet - leapfrog - predictor-corrector - Runge-Kutta …

  5. 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) 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. rB→A rA→B ΔEA ΔEB kinetic Monte Carlo equilibrium Monte Carlo N B A EB B EA A Molecular Dynamics t < N > 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. Ocus x CObr x 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 numbersr1 ,r2 [0,1[ 1 determine all possible processes i for given configuration of your system and build a list. Get all rates  (i) CalculateR = iG (i) and find process “k”: k k-1 G (i) r1 R   G (i) i=1 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. Brute force approach to rate constants: • 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! Rare event dynamics TS E IS FS

  13. X X eq. IS FS kT h  kTST = [ () e S/k ] e-E/kT ISFS 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. - Attempt frequency/preexponential factor koTST = e S/k ~ 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) kT h 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 = fdynkTST 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. i) Grid method: - Compute PES on a (regular) grid - Scales like: (no. of points) dim - e.g. 59 ~ 2 million grid points  often unfeasable ii) Drag method: - Choose appropriate reaction coordinate q - Constrain q and relax all other DOF - Move from IS to FS  highly dependent on good reaction coordinate  hysteresis! FS x x x x x x x x x x IS x TS x x Transition state search algorithms I: grid and drag

  16. FS x Ro x x x x R1 x x x x R2 x x x x x x x IS TS x 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. FS x x IS TS x F  F || spring true Transition state search algorithms III: nudged elastic band • - Initialize with several images {Ri} along a • straight-line interpolation • - Minimize • S(R1, …, RN) = iE(Ri) + ik/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) •  widely applied workhorse •  has problems, if energy varies largely along path, but very little perpendicular to it (kinky PES) x x x x x x x x G. Mills and H. Jónsson, Phys. Rev. Lett. 72, 1124 (1994)

  18. FS x x IS TS x 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. Exchange mechanism Ag(100) E = 0.73 eV Au(100) E = 0.65 eV 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. TAD Hyperdynamics 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. cus site bridge site O Ru 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. pO (atm) pO (atm) 2 2 105 105 1 1 pCO(atm) pCO(atm) 10-5 10-5 A ( pCO ,pO2 )-map of catalytic activity 600 K 10-15 10-10 10-5 1 10+5 10-10 10-5 1 10+5 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. pO (atm) 2 T = 350 K pCO = 10-10 atm T = 350 K pO2 = 10-10 atm 1 CObr/COcus pCO(atm) 10-10 CObr/- -4.0 - / - Obr/Ocus -5.0 10-20 -6.0 10-30 10-20 10-10 1 -7.0 -8.0 -9.0 -10.0 -11.0 log(TOF) …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

More Related