Car-Parrinello Molecular Dynamics Simulations (CPMD): Basics

1. Car-Parrinello Molecular Dynamics Simulations (CPMD): Basics Ursula Rothlisberger EPFL Lausanne, Switzerland

3. Literature Car-Parrinello: • R. Car and M. Parrinello • A unified approach for molecular • dynamics and density functional • Phys.Rev.Lett. 55, 2471 (1985) • P. Carloni, U. Rothlisberger and • M.Parrinello • The role and perspective of ab • initio molecular dynamics in the • study of biological systems • Acc. Chem.Res. 35, 455 (2002) • U Rothlisberger • 15 years of Car-Parrinello • simulations in Physics, Chemistry, • and Biology • Computational Chemistry: Reviews • of Current Trends, J. Leszczynski • (Ed.), World Scientific, Vol. 6, • (2001) p.33 • D. Marx and J. Hutter • Modern Methods and Algorithms • of Quantum • J. Grotendorst (Ed.), NIC • Forschungszentrum Jülich (2000) • p.301 • D. Sebastiani and U. Rothlisberger • Advances in density functional • based modelling techniques: • Recent extensions of the • Car-Parrinello approach • in P. Carloni, F. Alber ‘Medicinal • Quantum Chemistry’, Wiley-VCH, • Weinheim (2003)

4. When Quantum Chemistry Starts to Move... Classical MD Simulations Traditional QC Methods Car-Parrinello MD • improved optimization • finite T effects • thermodynamic & dynamic properties • solids & liquids • parameter-free MD • ab initio force field • no transferability problem • chemical reactions

5. When Newton meets Schrödinger... Sir Isaac Newton Erwin Schrödinger (1642 - 1727) (1887 - 1961)

6. Newt-dinger The ideal combination for Ab Initio Molecular Dynamics

7. Atoms, Molecules and Chemical Bonds Chemical Reaction Chemical Bonds Atoms N protons & neutrons + N electrons e-

8. Basic Principles of Quantum Mechanics

9. Wavefunctions and Probability Distributions Classical Mechanics: The position and velocity of the particle are precisely defined at any instant in time. Quantum Mechanics: The particle is better described via its wave character, with a wave function (r,t). The square of wave function is a measure for the probability P(r) to find the particle in an infinitesimal volume element dV around r. The total probability to find the particle anywhere in space integrates to 1.

10. Epot 0 n=3 n=2 n=1 h w n=0 q 0 Classical Mechanics Quantum Mechanics positions and momenta uncertainty have sharp defined relation values Continous energy spectrum energies are quantized q Newton`s Equations Schroedinger Equation n, E, m, h0

11. Classical Mechanics: Particle Motion ro,vo r(t),v(t) Position and velocity of a particle can be calculated exactly at any time t. Continuous energy

12. Solution 1: Time-dependent Schrödinger Eq. for a system of N nuclei and n electrons  not possible! Goal: Computational method that provides us with a microscopic picture of the structural and dynamic properties of complex systems

13. Product Ansatz for total wavefunction: Electronic Schrödinger Eq.: Electronic Hamiltonoperator: Approximations: 1) Born-Oppenheimer Approximation (1927): mel <<< mp  electronic and nuclear motion are separable Exceptions: Jahn-Teller instabilities, strong electron-phonon coupling, molecules in high intensity laser fields  nonadiabatic dynamics

14. Nuclear SchrödingerEq. Nuclear Hamiltonoperator: Nuclear Quantum Dynamics (review: Makri, Ann. Rev. Phys. 50, 167 (1999) Solve electronic Schrödinger Eq. for each set of nuclear coordinates potential energy surface (PES)

15. Empirical parameterization → force field based MD Calculate → Car-Parrinello Dynamics of an electron and a • ratio of the deBroglie wavelength proton:  classical approximation is better: m, n, E, T  Works surprisingly well in many cases! • zero point energy effects • (proton) tunneling  what cannot be described: Classical Nuclear Dynamics 2) Most atoms are heavy enough so that their motion can be described with classical mechanics  quantum corrections to classical results (Wigner&Kirkwood)  classical MD extended to quantum effects on equilibrium properties and to some extend also to quantum dynamics  path integral MD and centroid dynamics

16. First-Principles Molecular Dynamics • How do we do that? • 1) straight-forward: • solve electronic structure problem for a set of ionic coordinates • evaluate forces • move atoms Born-Oppenheimer Dynamics

17. Car - Parrinello Molecular Dynamics (1985) Lagrangian Formulation of Classical Dynamics Euler-Lagrange Equation:

18. Car - Parrinello Molecular Dynamics (1985) Extended Lagrangian Formulation Roberto Car Michele Parrinello

19. Can be integrated simultaneously (e.g. with Verlet, Velocity-Verlet algorithm etc..) Verlet algorithm dt ~0.1-0.2 fs Equations of Motion

20. Does this fictitious classical dynamics described via the extended Lagrangian have anything to do with the real physical dynamics??? • if •  total energy of the system becomes the real physical total energy  can be checked via energy conservation After initial wfct optimization, system is propagated adiabatically and moves within finite thickness Ke over the potential energy surface

21. What’s the price for it ? • systems sizes: • few hundred to few thousands of atoms (CP2K) • Time Steps: ~0.1 fs • Simulation Periods: few tens of ps

22. The Quantum Problem Stationary Solutions: Time-independent Schrödinger Eq. Variable Separation: Electronic Schrödinger Eq.: Electronic Hamiltonoperator: Product Ansatz for the wavefunction: Effective 1-particle model

23. The Quantum Problem Set of N coupled 1-particle equations: Basis Set Expansion:  Set of algebraic Eqs. Solved iteratively (self-consistent field) ca. 10’000-100’000 Plane-waves:  FFT Choice of QM method: DFT

24. DENSITY FUNCTIONAL THEORY

25. Walter Kohn and John Pople Nobelprize in Chemistry 1998

26. Literature on DFT: • Original Papers: • P.Hohenberg, W.Kohn, Phys.Rev.B 1964, 136, 864-871. • W.Kohn, L.J.Sham, Phys.Rev.A 1965, 140, 1133-1138. • Textbooks: • W.Kohn, P.Vashista, in Theory of the Inhomogeneous Electron Gas, N.H.March and S.Lundqvist (Eds), Plenum, New York 1983 • R.G.Parr, W.Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, New York 1989. • R.M.Dreizler, E.K.U.Gross, Density-Functional Theory, Springer, Berlin 1990. • W.Kohn, Rev.Mod.Phys. 1999, 71.

27. Density Functional Theory (DFT) Like Hatree-Fock: effective 1-particle Hamiltonian Let’s define a new central variable: Electron density Total electron density integrates to the number of electrons:

28. Theoretical foundations of DFT based on 2 theorems: • Hohenberg and Kohn (1964): • (Phys.Rev. 136, 864B) • The ground state energy of a system with N electrons in an external potential Vex isa unique functional of the electron density •  Vex determines the exact •  vice versa: Vex is determined within an additive constant by •  gs expectation value of any observable (i.e. the H) is a unique functional of the gs density

29. Variational principle: The total energy is minimal for the ground state density of the system • Kohn and Sham (1965): • (Phy. Rev. 1140, 1133A) • The many-electron problem can be mapped exactly onto: • an auxiliary noninteracting reference system with the same density (i.e. the exact gs density) • where each electrons moves in an effective 1-particle-potential due to all the other electrons

30. (1) (2) (4) (5) (3) (1) Kinetic energy of the non interacting system (2) External potential due to ionic cores (3) Hartree-term ~ classical Coulomb energy (4) exchange-correlation energy functional (5) Core -core interaction Kohn-Sham eqs:

31. Exchange and Correlation Exchange-Correlation Hole

32. Universal exchange-correlation functional, exact form not known!  local density approximation can be determined exactly: Exchange: (P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26, 376 (1930), E.P. Wigner, Trans. Fraraday Soc. 34, 678 (1987)) Correlation: (D.M. Ceperly, B.J. Alder, Phys. Rev. Lett. 45, 566 (1980), G.Ortiz, P. Ballone, Phys. Rev. B 50, 1391 (1994)) exact (numerical) results from Quantum Monte Carlo simulations

33. Parametrized analytic forms that interpolate between different density regimes are available (e.g. J.P. Perdew, A. Zunger, Phys. Rev. B. 23, 5084 (1981)) - in principle very crude approximation! - Exc of a non uniform system locally ~ uniform electron gas results - should ‘work’ only for systems with slowly varying density but: atoms and molecules are inhomogeneous systems! - works remarkably well in practice: Performance of LDA/LSDA  in general good structural properties:  bond lenghts up to 1-2%  bond angles ~ 1-2 degrees  torsional angles ~ a few degrees  vibrational frequencies ~ 10% ( phonon modes up to few %)

34.  cheap and good method for transition metals!: e.g. Cr2, Mo2 in good agreement with experiment ( not bound in HF, UHF!)  F2 re within 3% (not bound in HF)  atomization, dissociation energies over estimated (mainly due to errors for atoms), typically by 10-20%  hydrogen-bonding overestimated  van der Waals-complexes: strongly overestimated binding (e.g. noble gas dimers, Mg2, Be2: factor 2-4 Re[Å] De (eV) HF 1.465 -19.4 CCSD 1.560 -2.9 CCSD(T) 1.621 0.5 DFT 1.59 1.5 exp 1.679 1.4 Cr2 (Scuseria 1992)

35. Generalized Gradient Approximation (GGA) • correction function chosen to fulfill formal conditions for the properties of the ex-corr hole • Determination of parameters: • fully non empirical • fit to exact Ex-Corr energies for atoms • fit to experimental data (empirical) •  man different forms (B88, P86, LYP, • PW91, PBE, B3LYP etc..)

36. Time-independentelectronic Density-Functional Theory Schrödinger Equation:

37. Practical Implementation • periodic boundary conditions • plane wave basis set up to a given kinetic energy cutoff Ecut •  use of FFT techniques • convenient evaluation of different terms in real space • (Eex-corr, Eext) or in reciprocal space (Ekin, Ehartree) • typical real space grid: ~1003, ~10000-100000 pws • most of the time: FFT most time consuming step (NMlogM) • for large systems: orthogonalization ~N2 • well parallelizable (over number of electronic states and • first index of real space grid

38. Pseudo Potentials Framework • Chemical properties determined • by valence electrons • perform atomic all electron • calculation ab initio pseudo r > rc smooth fct r < rc rc • invert Schrodinger equation r(a.u.)

40. CPMD (3.9) (CP2K) www.cpmd.org • Features (see also online manual): • plane wave basis, pseudopotentials, pbc and isolated systems • LDA, LSD, GGAs (single point hybrid fct calcs possible) • geometry optimization • MD (NVE, NVT, NPT, Parrinello-Rahman) • path integral MD • different types of constraints and restraints • Property calculations: population analysis, multipole moments, • atomic charges, Wannier fcts, Fukui fcts etc.. Runs on essentially all platforms.. Most Recent Features: • QM/MM interface • Response function calculations: • NMR Chemical shifts, electronic spectra, vibrational spectra • Time Dependent DFT MD in excited states • History dependent Metadynamics

41. Mixed Quantum-Classical QM/MM- Car-Parrinello Simulations • Fully Hamiltonian • QM/MM Car-Parrinello • hybrid code • QM-Part: CPMD 3.8 • pbc, PWs, pseudo potentials • (n-1) CPUs • MM-Part: GROMOS96 + P3M, • AMBER (SYBIL, UFF) • 1 CPU Interface Region Quantum Region Classical Region A. Laio, J. VandeVondele, and U. Rothlisberger, J. Chem. Phys. 116, 6941 (2002); A. Laio, J. VandeVondele, and U. Rothlisberger, J. Phys. Chem. B (ASAP article) review in : M. Colombo et al. CHIMIA 56, 11 (2002)

42. - qp + qo included in Vext QM/MM Car-Parrinello Simulations monovalent pseudo potential QM/MM Lagrangian QM MM j l k i e- EQM: DFT EMM: Standard biomolecular Force Field

43. QM/MM Car-Parrinello in Combination with Response Properties • Variational Perturbation Theory: A. Putrino, D. Sebastiani, M. Parrinello, 113, 7103 (2000) • IR and Raman Spectra • Fukui Functions R. Vuilleumier, M. Sprik J.Chem.Phys. 115, 3454 (2001) • Chemical Shifts D. Sebastiani, M. Parrinello, J. Phys. Chem. A 105, 1951 (2001) • TDDFT: Spectra and Dynamics J. Hutter J.Chem.Phys. 118, 3928 (2003)

44. T. Ziegler et al. Theor. Chim. Acta 43, 261 (1977) (sum method) m1 m2 t1,2 QM/MM Car-Parrinello in Combination with Excited State Methods • ROKS HOMO-LUMO single excitations CP-version: I. Frank et al. J. Chem. Phys. 108, 4060 (1998) E(s) = 2E(m) - E(t) • LR-TDDFT-MD (Tamm-Dancoff Approximation) J. Hutter J. Chem.Phys. 118, 3928 (2003) L. Bernasconi et al. J. Chem.Phys. 119, 12417 (2003) • P-TDDFT-MD I. Tavernelli (to be published) Landau-Zener Surface Hopping Ehrenfest Dynamics

45. pB pA • MD as dynamical tool: Real-time simulation of • dynamical processes Limitations Due to Short Simulation Time  many processes lie outside time range • MD as sampling tool: • only small portion of phase space is sampled • relevant parts might be missed, • especially if there exist large • barriers between different • important regions • (e.g. different conformers) • ensemble average have large statistical errors • (e.g. relative free energies!)

46. Techniques from Classical MD: • Sampling at enhanced temperature • Rescaling of atomic mass(es) • Constraints • (Ryckaert, Ciccotti, Berendsen 1977) • (Sprik & Ciccotti 1998) • Umbrella Sampling • (Torrie&Valleau 1977) • Quasi-Harmonic Analysis • (Karplus, Jushick 1981) • Reaction Path Method • (Elber & Karplus 1987) • ‘Hypersurface Deformation’ • (Scheraga 1988, Wales 1990) • Multiple Time Step MD • (Tuckerman, Berne 1991) • (Tuckerman, Parrinello 1994) • Subspace Integration Method • (Rabitz 1993) • Local Elevation • (van Gunsteren 1994) • Conformational Flooding • (Grubmuller 1995) • Essential Dynamics • (Amadei&Berendsen 1996) • Path Optimization • (Olender & Elber 1996) • Multidimensional Adaptive • Umbrella Sampling • (Bartels, Karplus 1997) • Hyperdynamics • (Voter 1997) • (Steiner, Genilloud, Wilkins 1998) • (Gong & Wilkins 1999) • Transition Path Sampling • (Dellago, Bolhuis, Csajka, • Chandler 1998) • Adiabatic Bias MD • (Marchi, Ballone 1999) • Metadynamics • (Laio, Iannuzzi, Parrinello PNAS 99, • 12562 (2002), PRL 90, 23802 (2003)

47. Configurational Sampling Sampling of Rare Reactive Events Development of Enhanced Sampling Methods • multiple time step sampling • classical bias potentials and forces • double thermostatting • parallel tempering • Electronic Bias Potentials • Finite Electronic Temperature • Vibronic Coupling • Charge Restraint Two Dimensional Free Energy Surface with torsional potential bias T = 500K EA = 30 kcal/mol Peroxynitrous Acid 48ps  1kcal/mol J. Chem. Phys. 113 4863 (2000), J. Chem. Phys. 1157859-7864 (2001), J. Phys. Chem. B 106, 203-208 (2002), J. Am. Chem. Soc. 124, 8163 (2002)

48. Constraints Lagrangian: Equations of motion: • freeze out fast motions increase integration time step • ( linear speed up) • constrain slowest motion guide system ‘manually’ over barrier • (condition: slowest part of reaction coordinate  is known, all other • degrees of freedom have time to equilibrate along the path) • ( free energy differences via thermodynamic integration)  integral replaced by a discrete set of points (R)= ’ for a simple distance constraint (R)= lRI-RJl:

49. Umbrella Sampling: Bias Potentials(Torrie&Valleau 1977) (Grubmuller 1995, Voter 1997, Karplus 1997, Wilkins 1998…) • high overlap with original • ensemble • close match PES or free • energy surface • low dimensionality • computationally inexpensive ‘Ideal’ Bias: ‘Golden Rules’

50. Sampling Error in ab initio MD: Methyl Group Rotation in Ethane C2H6 (500K, 7.25 ps) Probability Distribution  HCCH EA = 2.8 kcal/mol 