320 likes | 502 Views
Bridging scales: Ab initio atomistic thermodynamics. Karsten Reuter Fritz-Haber-Institut, Berlin. gas phase. surface. bulk. General idea. Motivation: extend length scale consider finite temperature effects. Approach: separate system into sub-systems
E N D
Bridging scales:Ab initio atomistic thermodynamics Karsten Reuter Fritz-Haber-Institut, Berlin
gas phase surface bulk General idea • Motivation: • extend length scale • consider finite temperature effects • Approach: • separate system into sub-systems • (exploit idea of reservoirs!) • calculate properties of sub-systems • separately (cheaper…) • connect by implying equilibrium • between sub-systems Drawback: - no temporal information („system properties after infinite time“) - equilibrium assumption Ab initio atomistic thermodynamics and statistical mechanics of surface properties and functions K. Reuter, C. Stampfl and M. Scheffler, in: Handbook of Materials Modeling Vol. 1, (Ed.) S. Yip, Springer (Berlin, 2005). http://www.fhi-berlin.mpg.de/th/paper.html
Connecting thermodynamics, • statistical mechanics and density-functional theory Statistical Mechanics, D.A. McQuarrie, Harper Collins Publ. (1976) Introduction to Modern Statistical Mechanics, D. Chandler, Oxford Univ. Press (1987) M. Scheffler in Physics of Solid Surfaces 1987, J. Koukal (Ed.), Elsevier (1988)
Thermodynamics in a nutshell Internal energy (U) Etot(S,V) Enthalpy H(S,p) = Etot + pV (Helmholtz) free energy F(T,V) = Etot - TS Gibbs free energy G(T,p) = Etot - TS + pV Potential functions • Equilibrium state of system minimizes corresponding potential function • In its set of variables the total derivative of each potential function is simple • (derive from 1st law of ThD: dEtot = dQ + dW, dW = -pdV, dQ = TdS) • dE = TdS – pdV • dH = TdS + Vdp • dF = -SdT – pdV • dG = -SdT + Vdp • These expressions open the gate to a whole set of general relations like: S = - (F/T)V , p = - (F/V)T Etot = - T 2 (/T)V (F/T) Gibbs-Helmholtz eq. (T/V)S = - (p/S)V etc. Maxwell relations - Chemical potential µ = (G/ n)T,p is the cost to remove a particle from the system. Homogeneous system: µ = G/N (= g) i.e. Gibbs free energy per particle
Link to statistical mechanics A many-particle system will flow through its huge phase space, fluctuating through all microscopic states consistent with the constraints imposed on the system. For an isolated system with fixed energy E and fixed size V,N (microcanonic ensemble) these microscopic states are all equally likely at thermodynamic equilibrium (i.e. equilibrium is the most random situation). • Partition function Z = Z(T,V) = i exp(-Ei / kBT) Boltzmann-weighted sum • over all possible system states • F = - kBT ln( Z ) • If groups of degrees of freedom are decoupled from each other (i.e. if the energetic • states of one group do not depend on the state within the other group), then • Ztotal = (i exp(-EiA / kBT) ) (i exp(-EiB / kBT) ) = ZAZB • Ftotal = FA + FBe.g. electronic nuclear (Born-Oppenheimer) • rotational vibrational • N indistinguishable, independent particles: Ztotal = 1/N! (Zone particle)N
Computation of free energies: ideal gas I X Z= 1/N! (Znucl Zel Ztrans Zrot Zvib)N µ(T,p) = G / N = (F + pV) / N = ( - kBT ln( Z ) + pV ) / N i) Electr. free energyZel = i exp(-Eiel / kBT) Typical excitation energies eV >> kBT, only (possibly degenerate) ground state Fel Etot – kBT ln( Ispin ) contributes significantly Required input: Internal energyEtot Ground state spin degeneracy Ispin ii) Transl. free energyZtrans = k exp(-ħk2 / 2mkBT) Particle in a box of length L = V1/3 (L) Ztrans V ( 2 mkBT /ħ2 )3/2 Required input: Particle mass m
iv) Vibrational free energyZvib = i=1n exp(-(n + ½)ħi / kBT) Harmonic oscillator µvib(T) = i=1 ½ ħi + kBT ln(1 - exp(-ħi/kBT) ) Required input: M fundamental vibr. modes i M M Calculate dynamic matrix Dij = (mimj)-½ (2Etot/rirj)req Solve eigenvalue problem det(D – 1 i2) Computation of free energies: ideal gas II iii) Rotational free energyZrot = J (2J+1)exp(-J(J+1)Bo / kBT) Rigid rotator (Diatomic molecule) Zrot - kBT ln(kBT/Bo ) = 2 (homonucl.), = 1 (heteronucl.) Bo ~md2 (d = bond length) Required input: Rotational constant Bo (exp: tabulated microwave data)
vibr. + electr. + rotat. + rotat. Δµ(T, p = 1 atm) (eV) O2 CO + transl. + transl. Temperature (K) Temperature (K) Computation of free energies: ideal gas III O2 CO m (amu) 32 28 stretch (meV) 196 269 Bo (meV) 0.18 0.24 2 1 Ispin 3 1 µ = µ(T,p) = Etot + Δμ(T,p) Alternatively: Δ(T, p) = Δ(T, po) + kT ln(p/po) and Δ(T, po = 1 atm) tabulated in thermochem. tables (e.g. JANAF)
1/M 0 0 for p < 100 atm • Trouble maker… DFT phonon band structure Computation of free energies: solids G(T,p) = Etot + Ftrans + Frot + Fvib + Fconf + pV Ftrans Translational free energy Frot Rotational free energy pVV = V(T,p) from equation of state, varies little Fconf Configurational free energy EtotInternal energy Fvib Vibrational free energy Etot, Fvib use differences use simple models to approx. Fvib (Debye, Einstein) Solids (low T): G(T,p) ~ Etot + Fconf
II. Starting simple: Equilibrium concentration of point defects Solid State Physics, N.W. Ashcroft and N.D. Mermin, Holt-Saunders (1976)
Config. entropy: Fconf = kBT ln Z(n) with Z = = N (N-1) … (N-n-1) N! 1·2· …· n (N-n)!n! n/N = exp(-ED/kBT) Forget pV, use Stirling: ln N! N(lnN-1) Isolated point defects and bulk dissolution On entropic grounds there will always be a finite concentration of defects at finite temperature, even though the creation of a defect costs energy (ED > 0). How large is it? N sites, n defects(n <<N) Internal energy: Etot = nED Minimize free energy: (G/n)T,p = /nT,p (Etot – Fconf + pV) = 0
III. Slightly more involved: Effect of a surrounding gas phase on the surface structure and composition E. Kaxiras et al., Phys. Rev. B 35, 9625 (1987) X.-G. Wang et al., Phys. Rev. Lett. 81, 1038 (1998) K. Reuter and M. Scheffler, Phys. Rev. B 65, 035406 (2002)
Surface thermodynamics solid – gas solid – liquid solid – solid (“interface”) … A surface can never be alone: there are always “two sides” to it !!! Phase I / phase II alone (bulk): GI = NII GII = NIIII Total system (with surface): GI+II = GI+ GII + Gsurf (T,p) Phase II Phase I A = 1/A(GI+II - iNii) Surface tension (free energy per area)
i) Use reservoirs: ii) M = gM bulk i) O from ideal gas ii) Forget about Fvib and Fconf for the moment: (T,p) ( Esurf. – NMEM )/A – NOO(T,p) /A bulk (slab) Example: Surface in contact with oxygen gas phase surf.= 1/A[Gsurf.(NO, NM) – NOO - NMM] O2 gas surface bulk
( Esurf. – NMEM )/A – NOO/A (slab) bulk - clean (meV/Å2) Oxide formation on Pd(100) p(2x2) O/Pd(100) (√5 x √5)R27° PdO(101)/Pd(100) M. Todorova et al., Surf. Sci. 541, 101 (2003); K. Reuter and M. Scheffler, Appl. Phys. A 78, 793 (2004)
30 vib = Fvib/A = = 1/A d Fvib(T,)[ surf.( )- NRu bulk ( )] 20 (meV) Only the vibrational changes at the surface contribute to the surface free energy 10 Pd bulk 0 wPd(bulk) ~ 25 meV Vibrational contributions to the surface free energy Fvib(T,V) = d Fvib(T,) () • Use simple models for order of magnitude estimate e.g. Einstein model: ( ) = ( -)
20 10 + 50% 0 vib (meV/Å2) - 50% -10 -20 0 200 400 600 800 1000 Temperature (K) Surface induced variations of substrate modes < 10 meV/Å2 for T < 600 K - in this case!!!
H2O 40 30 vib. (meV/Å2) 20 OH 10 0 0 200 400 600 800 1000 Temperature (K) Surface functional groups Q. Sun, K. Reuter and M. Scheffler, Phys. Rev. B 67, 205424 (2003)
(meV/Å2) Fconf = kBT ln (N+n)!/(N!n!) No lateral interactions: Langmuir adsorption-isotherm <(Ocus)> = 1.0 T = 300 K T = 600 K 0.5 < (O) > 0.0 -1.5 -1.0 -0.5 O (eV) 1 1 + exp((Ebind - O)/kBT) Configurational entropy smears out phase transitions Configurational entropy and phase transitions clean surface O(1x1)
IV. Exploration of configuration space: Monte Carlo simulations and lattice gas Hamiltonians Understanding Molecular Simulation, D. Frenkel and B. Smit, Academic Press (2002) A Guide to Monte Carlo Simulations in Statistical Physics, D.P. Landau and K. Binder, Cambridge Univ. Press (2000)
Configuration space and configurational free energy Fconf = - kBT ln( Zconf ) • Canonic ensemble (constant temperature): • Partition function Z = Z(T,V) = i exp(-Ei / kBT) Boltzmann-weighted sum • over all possible system states • In general, the configuration space is spanned by all possible (continuous) positions • rN of the N atoms in the sample: • Z = ∫ drNexp(- E(r1,r2,…,rN) / kBT) • The average value of any observable A at temperature T in this ensemble is then • <A> = 1/Z ∫ drN A(r1,r2,…,rN) exp( -E(r1,r2,…,rN) / kBT)
Alternative: - random sampling (Monte Carlo) Example for integration by simple sampling Evaluating high-dimensional integrals: Monte Carlo techniques <A> = 1/Z ∫ drN A(r1,r2,…,rN) exp( -E(r1,r2,…,rN)/ kBT) • Problem: - numerical quadrature (on a grid) rapidly unfeasible • scales with: (no. of grid points)N • e.g.: 10 atoms in 3D, 5 grid points: 530 ~ 1021 evaluations
Measuring the depth of the Nile Thank you, Daan!!! • - Many evaluations where • integrand vanishes • Need extremely fine grid • (very inefficient) • Do random walk, and reject all • moves that bring you out of the water • - Provides average depth of Nile, but • NOT the total volume!! Finding a needle in a haystack: Importance sampling <A> = 1/Z ∫ drN A(r1,r2,…,rN) exp( -E(r1,r2,…,rN)/ kBT)
? ? Specifying “getting out of the water”:The Metropolis algorithm Etrial < Epresent: accept Etrial > Epresent: accept with probability exp[- (Etrial-Epresent) / kBT] • Some remarks: - With this definition, Metropolis fulfills „detailed • balance“ and thus samples a canonic ensemble • - If temperature T is steadily decreased during • simulation, upward moves become less likely and • one ends up with an efficient ground state search • („simulated annealing“)
Modern importance sampling Monte Carlo techniques allow to - efficiently evaluate the high-dimensional integrals needed for evaluation of canonic averages - properly explore the configuration space, and thus configurational entropy is intrinsically accounted for in MC simulations Major limitations: - still need easily 105 – 106 total energy evaluations - this is presently an unsolved issue. First steps in the direction of true „ab initio Monte Carlo“ are only achieved using lattice models In short:
A very simple lattice system: O / Ru(0001) • Consider only adsorption into hcp sites (for simplicity) • Simple hexagonal lattice, one adsorption site per unit cell • Questions: which ordered phases exist ? • order-disorder transition at which temperature ? Configuration space comprises: ordered structures (arbitrary periodicity) disordered structures How can we then sample the configuration space? BUT: only periodic structures accessible to direct DFT, and supercell size quite limited
Lattice gas Hamiltonians / Cluster expansions Expand total energy of arbitrary configuration in terms of lateral interactions Elatt = åiEo + 1/2 åi,j Vpair(dij)sisj + 1/3 åi,j,k Vtrio(dij,djk,dki)sisjsk + … • Algebraic sum (very fast to evaluate) • Ising, Heisenberg models • Conceptually easily generalized to • multiple adsorbate species • more complex lattices • (different site types etc.) …but how can we get the lateral interactions from DFT?
e.g. O / Ru(0001) LGH parametrization through DFT Since isolated clusters not compatible with supercell approach, exploit instead the interaction with supercell images in a systematic way: 3 • Compute many ordered structures • Write total energy as LGH expansion, • e.g. • E(3x3) = 2Eo + 2V1pair + 2V3pair • Set up system of linear equations • „Invert“ to get lateral interactions 3
clean p(2x2) p(2x1) clean clean p(2x2) p(2x1) p(2x1) disordered lattice-gas clean p(2x1) T (K) p(2x1) p(2x2) (meV/Å2) O (eV)
DFT parametrized lattice gas Hamiltonians enable - efficient sampling of configurational space - parameter-free prediction of phase diagrams - first treatment of disordered structures Major limitations: - systematics / convergence of LGH expansion - restricted to systems that can be mapped onto a lattice - expansion rapidly very cumbersome for complex lattices, multiple adsorbates, at defects/steps/etc. In short:
Ab initio atomistic thermodynamics Use DFT in the computation of free energies Suitably exploit equilibria and concept of reservoirs • allows any general thermodynamic reasoning • concentration of point defects at finite T • surface structure and composition in • realistic environments major limitations Vxc vs. kBT sampling of configurational space „only“ equilibrium Lecture 2 tomorrow: kinetics, time scales