CarParrinello Molecular Dynamics Simulations CPMD Basics Ursula Rothlisberger
Car-Parrinello Molecular Dynamics Simulations (CPMD): Basics Ursula Rothlisberger EPFL Lausanne, Switzerland
Literature Car-Parrinello: • R. Car and M. Parrinello A unified approach for molecular dynamics and density functional Phys. Rev. Lett. 55, 2471 (1985) • 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) • 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
When Quantum Chemistry Starts to Move. . . Traditional QC Methods Classical MD Simulations 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
When Newton meets Schrödinger. . . Sir Isaac Newton (1642 - 1727) Erwin Schrödinger (1887 - 1961)
Newt-dinger The ideal combination for Ab Initio Molecular Dynamics
Atoms, Molecules and Chemical Bonds Atoms + e- N protons & neutrons N electrons Chemical Bonds Chemical Reaction
Basic Principles of Quantum Mechanics
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 d. V around r. The total probability to find the particle anywhere in space integrates to 1.
Classical Mechanics positions and momenta have sharp defined values Continous energy spectrum Epot Quantum Mechanics uncertainty relation energies are quantized n=3 n=2 n=1 n=0 q 0 Newton`s Equations hw q 0 Schroedinger Equation n , E , m , h 0
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
Goal: Computational method that provides us with a microscopic picture of the structural and dynamic properties of complex systems Solution 1: Time-dependent Schrödinger Eq. for a system of N nuclei and n electrons not possible!
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 Product Ansatz for total wavefunction: Electronic Schrödinger Eq. : Electronic Hamiltonoperator:
Solve electronic Schrödinger Eq. for each set of nuclear coordinates potential energy surface (PES) Nuclear Schrödinger. Eq. Nuclear Hamiltonoperator: Nuclear Quantum Dynamics (review: Makri, Ann. Rev. Phys. 50, 167 (1999)
Empirical parameterization → force field based MD Calculate → Car-Parrinello Dynamics Classical Nuclear Dynamics 2) Most atoms are heavy enough so that their motion can be described with classical mechanics • ratio of the de. Broglie wavelength of an electron and a proton: classical approximation is better: m , n , E , T Works surprisingly well in many cases! what cannot be described: • zero point energy effects • (proton) tunneling 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
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
Car - Parrinello Molecular Dynamics (1985) Lagrangian Formulation of Classical Dynamics Euler-Lagrange Equation:
Car - Parrinello Molecular Dynamics (1985) Extended Lagrangian Formulation Roberto Car Michele Parrinello
Equations of Motion Can be integrated simultaneously (e. g. with Verlet, Velocity-Verlet algorithm etc. . ) Verlet algorithm dt ~0. 1 -0. 2 fs
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
What’s the price for it ? • systems sizes: few hundred to few thousands of atoms (CP 2 K) • Time Steps: ~0. 1 fs • Simulation Periods: few tens of ps
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
The Quantum Problem Set of N coupled 1 -particle equations: Basis Set Expansion: Set of algebraic Eqs. Solved iteratively (self-consistent field) Plane-waves: ca. 10’ 000 -100’ 000 FFT Choice of QM method: DFT
DENSITY FUNCTIONAL THEORY
Walter Kohn and John Pople Nobelprize in Chemistry 1998
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.
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:
Theoretical foundations of DFT based on 2 theorems: Hohenberg and Kohn (1964): (Phys. Rev. 136, 864 B) • The ground state energy of a system with N electrons in an external potential Vex is a 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
• Variational principle: The total energy is minimal for the ground state density of the system Kohn and Sham (1965): (Phy. Rev. 1140, 1133 A) 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 -particlepotential due to all the other electrons
(1) (3) (2) (4) (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: (5)
Exchange and Correlation Exchange-Correlation Hole
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
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 %)
cheap and good method for transition metals!: e. g. Cr 2, Mo 2 in good agreement with experiment ( not bound in HF, UHF!) F 2 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, Mg 2, Be 2: factor 2 -4 Re[Å] De (e. V) HF 1. 465 -19. 4 Cr 2 CCSD 1. 560 -2. 9 CCSD(T) 1. 621 0. 5 (Scuseria 1992) DFT 1. 59 1. 5 exp 1. 679 1. 4
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 (B 88, P 86, LYP, PW 91, PBE, B 3 LYP etc. . )
Time-independent electronic Schrödinger Equation: Density-Functional Theory
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 (NMlog. M) • for large systems: orthogonalization ~N 2 • well parallelizable (over number of electronic states and first index of real space grid
Pseudo Potentials Framework ab initio pseudo • Chemical properties determined by valence electrons • perform atomic all electron calculation r > rc smooth fct r < rc rc • invert Schrodinger equation r(a. u. )
ABINIT www. abinit. org CASTEP [ i ] Molecular Simulations Inc. CPMD [ i i ] M. Parrinello, MPI Stuttgart, Germany and IBM Zurich Research Laboratory, Switzerland www. cpmd. org Free software [ i i i ] Fritz-Haber Institute Berlin, Germany fhim@fhi-berlin. mpg. de CP 2 K Fhi 98 md JEEP François Gygi, Lawrence Livermore National Laboratory, USA NWCHEM Pacific Northwest National Laboratory, USA PAW [ i v ] P. E. Blöchl, Clausthal University of Technology, Germany SIESTA [ v ] P. Ordejon, Institut de Ciencia de Materials de Barcelona, Spain VASP [ v i ] J. Hafner, University of Vienna, Austria
CPMD (3. 9) (CP 2 K) 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
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: GROMOS 96 + P 3 M, AMBER (SYBIL, UFF) 1 CPU Interface Region Quantum Region Classical Region A. Laio, J. Vande. Vondele, and U. Rothlisberger, J. Chem. Phys. 116, 6941 (2002); A. Laio, J. Vande. Vondele, and U. Rothlisberger, J. Phys. Chem. B (ASAP article) review in : M. Colombo et al. CHIMIA 56, 11 (2002)
QM/MM Car-Parrinello Simulations QM/MM Lagrangian QM monovalent pseudo potential MM j l k i e. EQM: DFT EMM: Standard biomolecular Force Field + qp - qo included in Vext
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)
QM/MM Car-Parrinello in Combination with Excited State Methods • ROKS HOMO-LUMO single excitations T. Ziegler et al. Theor. Chim. Acta 43, 261 (1977) (sum method) CP-version: I. Frank et al. J. Chem. Phys. 108, 4060 (1998) • 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 m 1 m 2 t 1, 2 E(s) = 2 E(m) - E(t)
Limitations Due to Short Simulation Time • MD as dynamical tool: Real-time simulation of dynamical processes 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!) p. A p. B
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)
Development of Enhanced Sampling Methods Configurational Sampling • multiple time step sampling • classical bias potentials and forces • double thermostatting • parallel tempering Two Dimensional Free Energy Surface with torsional potential bias Sampling of Rare Reactive Events Electronic Bias Potentials • Finite Electronic Temperature • Vibronic Coupling • Charge Restraint T = 500 K EA = 30 kcal/mol Peroxynitrous Acid 48 ps 1 kcal/mol J. Chem. Phys. 113 4863 (2000), J. Chem. Phys. 115 7859 -7864 (2001), J. Phys. Chem. B 106, 203 -208 (2002), J. Am. Chem. Soc. 124, 8163 (2002)
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)= l. RI-RJl:
Umbrella Sampling: Bias Potentials (Torrie&Valleau 1977) (Grubmuller 1995, Voter 1997, Karplus 1997, Wilkins 1998…) ‘Ideal’ Bias: ‘Golden Rules’ • high overlap with original ensemble • close match PES or free energy surface • low dimensionality • computationally inexpensive
Sampling Error in ab initio MD: Methyl Group Rotation in Ethane C 2 H 6 (500 K, 7. 25 ps) Probability Distribution HCCH EA = 2. 8 kcal/mol
Atomic Bias Potentials Methyl Group Rotation in Ethane C 2 H 6 (500 K, 7. 25 ps) Before correction After correction Torsional Bias 0. 0017 au
Bias Potentials from Classical Force Field Peroxynitrous Acid ONOOH Trajectory in biased space (48 ps) Free Energy Surface J. Vande. Vondele, U. R. J. Chem. Phys. 113 4863 (2000)
CAFES: Canonical, Adiabatic Free Energy Sampling • Partitioning into reactive system / environment adiabatic decoupling Slow dynamics of the reactive subsystem different temperatures TR/TE (2 thermostats) Sampling efficiency at TR can be estimated Ea= 20 kcal/mol, TE=300 K, TR=1200 K -> 1013 J. Vande. Vondele, U. R. J. Phys. Chem. B 106, 203 (2002)
Nucleophilic substitution with anchimeric assistance • QMMM SPC/CPMD • CAFES 100 / 2000 K / 300 K • ~22 kcal/mol • shows that the reaction coordinate is not simple
Transition State Path Sampling Given: - initial state A - final state B - one path connecting the two generate the ensemble of ‘reactive paths’ calculate transition rates
Dispersion Interactions in DFT MM QM
Suggested Remedies • add -C 6/r 6 -term (with damping function) (Le. Sar 1984 , Sprik 1996, Scoles 2001, Parrinello 2003, Wang, York 2004… • specially designed (local) functionals (PW 91, PBE, m. PBE, X 3 LYP, …) • density partitioning schemes (Wesolowski 2003…) • nonlocal correlation functionals for special cases (Langreth, Lundvist 2000, 2003…) • perturbation calculation of dispersion forces (Kohn 1998, Szalewicz 2003…)
Optimized Effective Atom Centered Potentials Expansion in linear combination of atom-centered (nonlocal) potentials Analytic pseudopotentials by Goedecker et al.
Optimization Penalty Functional Linear density response calculated via first order perturbation theory with perturbation Hamiltonian For :
s 1 = -0. 00352 s 2 = 3. 280 BLYP OECP MP 2 Is this potential transferable? ? ?
BLYP OECP MP 2
BLYP z = 3. 3 A OECP E=32 me. V/atom exp z = 3. 35 A E=35 me. V/atom
Reference system: Ar 2 1 additional f-channel: s 1 = -0. 00206 s 2 = 2. 902 BLYP OECP MP 2 Klopper et al. J. Chem. Phys. 101, 9747 (1994) OECP MP 2
What about the intramolecular geometry? ? D << 0. 01 A Bond lengths in benzene: What about the electronic properties? ? Dipole moment: benzene-Ar Quadrupole moment: benzene Polarizability: argon axx- ayy benzene azz benzene-Ar BLYP OECP MP 2 0. 047 0. 035 0. 037 -5. 35 -5. 50 -6. 46 12. 30 39. 18 55. 0 12. 31 38. 45 58. 1 11. 15 35. 07 59. 2
Formaldimine. Excited state dynamics after excitation S 0→S 1 The region of conical intersection, CI, is reached only in case of non-thermostatted trajectories. relaxation to product geometry α Φ Φ start on S 1 back to reactant geometry Increasing kinetic energy α ω ω Φ minimum on S 1 α
Landau-Zener SH energy Classical treatment for the derivation of an analytical formula for the transfer rate which is valid for any value of the coupling matrix element spanning the range between adiabatic and nonadiabatic ET. UA UD q* q 0 The asymptotic value for the survival probability of the electron for remaining at the donor The donor survival probability is
Units: atomic units used throughout
Transition Rate Constants Reactive Flux Correlation Function • Can be calculated with trajectories starting at the TS • Is difficult if a RC/TS cannot be defined.
Rate constants in the TPE • C(t) = the fraction of trajectories of length t, starting in A, that arrives in B • Can be calculated with a reversible work calculation. Contrary to direct MD, the computational efficiency does not depend on the height of the barrier.
- Slides: 71