and barostats CP: Car-Parrinello Atomic units MD and CP textbo Molecular dynamics: Car-Parrinello method Víctor Luaña (†) & Alberto Otero-de-la-Roza (‡) (†) Departamento de Química Física y Analítica, Universidad de Oviedo (‡) University of California at Merced European school on Theoretical Solid State Chemistry ZCAM, Zaragoza, May 13–17, 2013 V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 1 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Content 1 Statistical mechanics 2 MD: Hits 3 MD: Classical mechanics formulations 4 MD: Source of potential energy 5 MD: System and initialization 6 MD: Solving numerically the equations of motion: Verlet (1967) 7 MD: Steps in a MD simulation 8 Ensemble constraints in MD: Thermostats and barostats 9 The Car-Parrinello formulation V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 2 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Statistical mechanics Gibbs ensembles method: microcanonical: Ω(N, V, E) ⇒ S = k ln Ω canonical: Q(N, V, T) = i Ω(N, V, E)e−Ei/kT ⇒ A = −kT ln Q grand canonical: Ξ(µ, V, T) = N Q(N, V, E)eNµ/kT → pV = −kT ln Ξ iso-pT: ∆(N, p, T) = V Q(N, V, E)e−pV/kT ⇒ G = −kT ln ∆ Independent particles: Boltzmann (classical), Maxwell-Boltzmann (velocities), Fermi-Dirac, Bose-Einstein, ... Molecular dynamics: • Integrate the equations of motion for all the particles in the system. • Check the ergodicity of the statistical sampling: A = lim τ→∞ 1 τ τ 0 A(τ)dτ • Formally equivalent to a micro- canonical ensemble, use constraints (thermostats, barostats, ...) for other ensembles. Monte-Carlo methods: Choose randomly a Markov chain of states (the new state depends only on the cur- rent one, the process is memoryless) • Accepting a new state depends on the Metropolis rule: accept always if the energy diminishes (∆E = Ei − Ei−1 < 0) and ac- cept movements such that ∆E > 0 with a probability P(∆E) ∝ e−∆E/kT . • Check the ergodicity of the statis- tical sampling: A = lim N→∞ 1 N N i=1 Ai. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 4 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Molecular Dynamics: classical formulations Any classical formulation can be used: Newton: U = U(RN) ⇒ Fi = −∇Ri U = Mi ¨ Ri , (1) Lagrange: L = T − U ⇒ d dt ∂L ∂ ˙ Ri = ∂L ∂Ri , (2) Hamilton: H = i pi qi − L ⇒ ˙ Ri = ∂H ∂Pi = Pi Mi ; ˙ Pi = ∂H ∂Ri = Fi . (3) pi = ∂L/∂ ˙ qi is a generalized momentum. The Lagrange formulation is, probably, the most common for introducing thermostats and barostats. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 8 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Isaac Newton (1643–1727) Joseph Louis Lagrange (1736–1813) William Rowan Hamilton (1805–1865) V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 9 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD: Source of potential energy I Empirical force ﬁelds (MM: molecular mechanics) U(RN) = bond Ebond + Eangle + Edihedral + nonbonded Eelectrostatics + EvdW (4) = N i Ei(Ri) 1body + 1 2 N i>j Eij(|Ri − Rj|) 2body + 1 6 N i>j>k Eijk(Ri, Rj, Rk) 3body +... Typically developed for particular kinds of compounds: aminoacids, sugars, organometallics, ... Usually developed and validated with an empirical set of data. Some popular: AMBER (As- sisted Model Building and Energy Reﬁnement, http://ambermd.org/), CHARMM (Chem- istry at Harvard macromolecular mechamics, http://www.charmm.org/), MM2, MM3, MM4 (Allinger’s MM), GROMACS (GRÖningen MAchine for Chemical Simulations, http://www. gromacs.org/), UFF (Universal FF), ... Codes: amber (ambermd.org), dl_poly (www.ccp5.ac.uk/DL_POLY_CLASSIC/), gromacs (www.gromacs.org), lammps (lammps.sandia.gov), moldy (www.ccp5.ac.uk/moldy/ moldy.html), namd (www.ks.uiuc.edu/Research/namd/), ... V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 11 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD: System and initialization Periodic boundary conditions: N atoms per cell. Atom prop- erties (i ∈ [1, N]): {R, ˙ R, ˙ F, ...}(i), i ≡ mod(i, N) + 1. (8) Large unit cells are used to simulate ﬁnite systems: Avoid interactions between atoms in different cells. LUC’s can be used to simulate molecules, 1D chains, slabs, surfaces, impurities and defects. A crystal conﬁguration + some randomness can be used to create the initial geometry: fcc, sc, bcc, and hcp are typical start confugurations for gases and liquids. Maxwell-Boltzmann distribution is used for the starting ve- locities: f(v)dv = 4π m 2πkT 3/2 v2 exp − mv2 2kT (9) Simulation can include periods of heating and or compress- ing: Every time the thermodynamic conditions are changed, a period of equilibration must follow, that is discarded for obtaining the average properties. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 14 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Solving numerically the equations of motion: Verlet (1967) I Taylor expansion: r(t ± ∆t) = r(t) + ˙ r(t)∆t + f(t) 2m ∆t2 + 1 3! ∂3r ∂t3 ∆t3 + O(∆t4) (10) Verlet algorithm: [O(∆t4)] r(t + ∆t) = 2r(t) − r(t − ∆t) + f(t) m ∆t2, ˙ r(t) = r(t + ∆t) − r(t − ∆t) 2∆t . (11) Verlet integrator is symplectic: satisﬁes time reversibility and it is robust against time step increase. Velocity Verlet algorithm: r(t + ∆t) = r(t) + ˙ r∆t + f(t) 2m ∆t2, ˙ r(t + ∆t) = ˙ r(t) + f(t) + f(t + ∆t 2m )∆t. (12) Position and velocity Verlet schemes are equivalent. Lyapunov instability: Trajectories through phase space can be very sensitive to the initial conditions, so trajectories that are initially very close will diverge exponentially as time goes. For a correct statistical prediction we only need we only need the trajectory to stay close to the true path long enough compared to the time for the unstability to develop. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 16 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Solving numerically the equations of motion: Verlet (1967) II The stability depends on the time step and the time of simulation. Verlet algorithm conserves energy only on the ∆t → 0 limit. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 17 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Solving numerically the equations of motion: Verlet (1967) III Checking energy conservation is a good test for ∆t. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 18 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Solving numerically the equations of motion: Verlet (1967) IV We are interested in the stability of properties. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 19 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Steps in a MD simulation I -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0 500 1000 1500 2000 2500 3000 3500 4000 175 180 185 190 195 200 205 210 215 220 225 230 Energy Volume Steps fcc Ar 32 (200 K) Equilibration Production E V U K Discard the equlibration period when obtainig the average properties. Temperature: 3 2 Nk T = 1 2 N i=1 Mi ˙ R 2 rms displacement: RMSDα(t) = 1 Nα Nα α=1 (Rα(t) − rα )2 LJ ﬂuid close to the triple point (T = 0.71, ρ = 0.844) Radial correlation function: g(r) = V N2 i>j δ(r − Rij) V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 21 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Steps in a MD simulation II V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 22 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Ensemble constraints in MD: Thermostats and barostats I Thermostats: Velocity rescaling: Every few steps all the velocities are scaled by a common factor (˙ ri → λ˙ ri: T(λ) = λ2T(1)) so the average temperature corresponds to the target value. Easy and practical. Not a true thermostat. Possibly large ﬂuctuations of temperature. Nosé-Hoover (1984): Introduce a heatbath as an integral part of the system, adding an artiﬁcial variable s, associated with a mass Q > 0. The s variable plays the role of a time scaling: dt → sdt, with the rest of dynamic variables being transformed as Ri → Ri, ˙ Ri → s−1 ˙ Ri, ˙ s → s−1˙ s. The extended system is described by a lagrangian L = N i 1 2 mis2 ˙ R2 i + U(R) + 1 2 Q˙ s2 − gkt ln s (15) which includes a kinetic and potential energy associated to s. The time-evolution of s is described by a second-order equation. Heat may ﬂow in and out of the system in an oscilla- tory fashion, leading to nearly periodic temperature ﬂuctuations. Not (always) ergodic. Choose carefully the Q mass, or use a chain of NH thermostats. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 24 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Ensemble constraints in MD: Thermostats and barostats II An ergodic simulation must visit efﬁciently and randomly the phase space: Barostats: V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 25 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo Ensemble constraints in MD: Thermostats and barostats III Andersen (1980): Deﬁne scaled coordinates ρi = Ri/V1/3, declare a enlarged hamiltonian: L = 1 2 Ω2/3 N i=1 mi ˙ ρ2 i + U(Ω1/3ρ) + 1 2 M ˙ Ω2 − p0Ω (16) where Q represents the cell volume, p0 is the target pressure, and Ω is the actual cell volume. Designed for an isotropic deformation of the cell. Parrinello-Rahman (1981): Generalization to allow anisotropic deformations. Let h = (a, b, c) be the row-arranged cell vectors so G = hT h is the cell metric tensor. Atomic positons are RI = hxI, where x = (x, y, z)T is the column vector of cell coordinates, each atom in the main cell having 0 ≤ x, y, z < 1. PR barostat corresponds to the Lagrangian L = 1 2 I Mi ˙ xT I G˙ xI − U[h, {xI}] + 1 2 WTr ˙ hT ˙ h − p0Ω (17) formed by the kinetic and potential energy of the ions, plus the kinetic and potential associated to the deformation and the isotropic expansion of the cell. W is a ﬁctitious mass that controls the barostat ﬂuctuations. It can be used to study reconstructive phase transitions. Wentzcovich (1991): Similar to PR, but it uses the strain instead of h as the dynamical variable. This solves the problem of PR that the trajectory is not uniquely deﬁned. This barostat is included in quantum espresso. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 26 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation I V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 28 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation II "The 2009 Dirac Medal recognizes the joint contributions of Roberto Car and Michele Parrinello in developing the ab initio simulation method in which they combined, elegantly and imaginatively, the quantum mechanical density functional method for the calculation of the electronic properties of matter with molecular dynamics methods for the Newtonian simulation of atomic motions. The Car-Parrinello method has had an enormous impact, joining together the ﬁelds of simulation and of electronic structure theory, and has given rise to a variety of applications well beyond condensed matter physics." V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 29 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation III The idea behind Car-Parrinello dynamics is to put together the ab initio electronic structure and the molecular dynamics methods. Due to the big mass difference between electrons and ions, direct solving the dynamic equations of both kinds of particles is not possible. Rather, CP designed an extended Hamiltonian where the electronic degrees of freedom are introduced through ﬁctitious dynamical variables. At difference with Born-Oppenheimer MD, no electronic minimization is needed at each MP step, both processes taking place simultaneously at CPMD. CPMD starts with an standard electronic minimization that brings the system to the Born-Oppenheimer ground state surface. After that, the ﬁctitious dynamics keeps the system on the electronic ground state. Each MD step moves to a different ionic conﬁguration, the forces being consistently accurate, even if the conﬁguration is far away from equilibrium. The ﬁctitious mass of the electrons must be small enough to keep the trajectory adiabatic, avoiding the transfer of energy from the ionic to the electronic degrees of freedom. This also requires as timestep smaller than the usual values in other MD methods: 1–10 fs. CPMD equations derive from a Lagrangian: L = TI + Te + U + C(ortho) + ... (18) = 1 2 nuc I MI ˙ R2 I + 1 2 orb i µ ˙ ψi| ˙ ψi − E[{RI}, {ψi}] + orb i>j Λij{ ψi|ψj − δij} + other constraints V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 30 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation IV where T kinetic energy (I: ions, e: electrons); U: potential energy (HF, KS, ...); C(ortho): orbital orthogonality constraint. The corresponding equations of state are: MI ¨ RI = −∇IE[{RI}, {ψi}], µ ¨ ψi = − ∂E ∂ψ∗ i + orb j Λijψj. (19) that must be solved in sequence using a variation of Verlet algorithm. Why does CPMD works? Pastore, Simargiassi and Buda [Phys. Rev. A 44 (1991) 6334] show that electrons and ions move on a differente timescale, avoiding reso- nance and making prohibitively slow the trans- fer of energy between both dynamical systems: f(ω) = ∞ 0 dt cos(ωt) i ˙ ψi; t| ˙ ψi; 0 In calculations for insulators and semiconductors using a basis of planewaves, the frequency for the electronic oscillations occurs in the range µω2 e ∈ [Egap, Ecut], where Egap is the band gap and Ecut the planewave cutoff. As the timestep is inversely proportional to the fre- quency, the relation ∆tmax ∝ µ/Ecut governs the largest possible timestep. Diamond phase (Z = 8) Si CPMD(LDA) cal- culation: dt = 0,3 fs, µ = 300 au, τ = 6,3 ps (20 000 steps). Electronic (triangle) and atomic normal modes and Φ-DOS. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 31 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation V Pastore et al. [PRA 44(1991)6334] CPMD calculation of diamond-Si. Oscillations do occur: Econs = orb i 1 2 µ ˙ ψi| ˙ ψi + I 1 2 MI ˙ R2 I + Ψ|H|Ψ , Ephys = I 1 2 MI ˙ R2 I + Ψ|H|Ψ = Econs − Te, Ve = Ψ|H|Ψ , Te = 1 2 i µ ˙ ψi| ˙ ψi . Econs and Ephys are motion constants, Ve and Te show anticorrelated oscillations. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 32 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation VI Ab initio MD techniques (BOMD or CPMD) are required whenever the electronic state is important. For instance, a chemical reaction, a phase transition where chemical bonds are created or broken, or the simulation of spectroscopic phenomena. There are severe limitations on the system size and number of steps of the simulation. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 33 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo The Car-Parrinello formulation VII Classical MD: • hardwired potential, • forces best close to equilibrium, • no electronic degrees of freedom, • not adapted to creation/breaking bonds, • scale: 106 atoms, nm, µs. CPMD: • on the ﬂy potential, • forces correct even far from equilibrium, • electronic degrees of freedom, • adapted to creation/breaking bonds, • no need to solve HΨ = EΨ after each step, • scale: 102 atoms, pm, ps. BOMD: • on the ﬂy potential, • forces correct even far from equilibrium, • electronic degrees of freedom, • adapted to creation/breaking bonds, • solve HΨ = EΨ after each step, • scale: 102 atoms, pm, ps. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 34 / 36
and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD and CP textbooks Statistical thermodynamics: D. A. McQuarrie, Statistical mechanics (Harper Collins, 1976). D. Chandler, Introduction to modern statistical mechanics (Oxford UP, 1987). Molecular dynamics: M. P. Allen and D. J. Tildesley, Computer simulation of liquids (Clarendon Press, 1987). D. Frenkel and B. Smit, Understanding molecular simulation (Academic Press, 2005). Raugei Simone, Introduction to molecular dynamics simulations (people.sissa.it/~raugei/ lecture_notes/md_notes.pdf) G. Bussi, Theory and tips for molecular dynamics (sites.google.com/site/giovannibussi downloads/slides.pdf, 2008). Molecular simulations: A. R. Leach, Molecular modeling (Prentice Hall, 2001). T. Schlick, Molecular modeling and simulation (Springer, 2002). Car-Parrinello and ab initio molecular dynamics: D. Marx and J. Hutter, "Ab initio molecular dynamics: Theory and implementation", in Mod- ern methods and algorithms on quantum chemistry, J. Grotendorst (Ed.), (John von Neu- mann Institute, NIC series vol. 1 & 3, 2000) http://www2.fz-juelich.de/nic-series/ Volume3/marx.pdf. Raugei Simone, Introduction to Car-Parrinello molecular dynamics simulations (http://people sissa.it/~raugei/lecture_notes/cp.pdf) V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 36 / 36