Molecular dynamics: Car-Parinello method

Transcript

1. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
2. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
3. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
4. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD: History 1872–1875 Boltzmann equation (S = k ln Ω) 1902 Gibbs book on statistical mechanics 1952 MANIAC operational at Los Alamos 1953 [Metropolis, Rosenbluth (2), Teller (2)] Metropolis Monte Carlo method 1954 [Fermi, Pasta, Ulam] experiment on ergodicity 1956 [Livermore: Alder, Wainwright] hard spheres dynamics 1959 [Brookhaven: Vineyard] radiation damage in Cu 1964 [Argonne: Rahman] liquid Argon 1967 [Verlet] correlation functions 1970 [Livermore: Alder, Wainwright] long time tails 1977 [Ryckaert, Ciccoti and Berendsen] SHAKE: include constraints avoiding matrix inversion 1980 [Anderson, Rahman, Parrinello] constant pressure MD 1983 [Nosé, Hoover] constant temperature MD 1985 [Sissa: Car, Parrinello] CPMD 1997 [Berendsen] LINCS: non-iterative efficient way to apply bond constraints V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 6 / 36
5. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
6. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
7. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD: Source of potential energy I Empirical force fields (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 Refinement, 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
8. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo MD: Source of potential energy II Born-Oppenheimer dynamics (BOMD) ˆ HΨ(x|R) = E(R)Ψ(x|R) HF,... =⇒ DFT:KS SCF: ˆ Fiψi = iψi ⇒ FI = ∇IE(R) (5) LBO(Ψ, R, ˙ R) = N i=1 1 2 Mi ˙ Ri − EBO(Ψ, R) (6) Car-Parrinello dynamics (CPMD) LCP({ψ}, { ˙ ψ}, R, ˙ R) = N I=1 1 2 MI ˙ RI + occ i µ ˙ ψi| ˙ ψi − EKS BO ({ψi}, R) (7) Codes: cpmd (www.cpmd.org), cp2k (www.cp2k.org), siesta (www.icmab.es/dmmis/ leem/siesta/), quantum espresso (www.pwscf.org), ... Rare event methods (metadynamics, ...): Designed to study phenomena that would require imposibly large simulation times to occur using conventional MD techniques. Example: protein folding. Codes: plumed (www.plumed-code.org), ... Visualization codes: jmol (jmol.org), rasmol (www.umass.edu/microbio/rasmol/), vmd (www.ks.uiuc.edu/Research/vmd/), ... V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 12 / 36
9. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 finite 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 configuration + 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
10. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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: satisfies 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
11. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
12. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
13. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
14. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 fluid 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
15. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
16. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 fluctuations of temperature. Nosé-Hoover (1984): Introduce a heatbath as an integral part of the system, adding an artificial 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 flow in and out of the system in an oscilla- tory fashion, leading to nearly periodic temperature fluctuations. 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
17. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo Ensemble constraints in MD: Thermostats and barostats II An ergodic simulation must visit efficiently and randomly the phase space: Barostats: V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 25 / 36
18. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo Ensemble constraints in MD: Thermostats and barostats III Andersen (1980): Define 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 fictitious mass that controls the barostat fluctuations. 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 defined. 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
19. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
20. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 fields 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
21. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 fictitious 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 fictitious dynamics keeps the system on the electronic ground state. Each MD step moves to a different ionic configuration, the forces being consistently accurate, even if the configuration is far away from equilibrium. The fictitious 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
22. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
23. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
24. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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
25. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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 fly 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 fly 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
26. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo and barostats CP: Car-Parrinello Atomic units MD and CP textbo Atomic units Length (a0 ): 1 bohr = 0,529 177 208 59 = 5,291 772 085 9 · 10−11 m. Energy (Eh ): 1 hartree = 27,210 7 eV = 2 625,500 kJ/mol = 219 474,63 cm−1 = 6,579 684 PHz; 1 rydberg = 0,5 hartree. Time ( /Eh ): 1 aut = 0,024 188 843 265 05 fs. Velocity (a0Eh/ ): 1 auv = 2,187 691 263 3 · 106 m/s. Temperature (Eh/k): 1 auT = 3,157 746 4 · 105 K. Pressure (Eh/a3 0 ): 1 auP = 29,421 912 TPa. V. Luaña & A. Otero-de-la-Roza () Molecular dynamics: Car-Parrinello method ZCAM, Zaragoza 2013 35 / 36
27. Formulations Potential energy Initialization Verlet algorithm MD: Steps MD: Thermo 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