Car-Parrinello Method and Applications

1. Car-Parrinello Method and Applications Moumita Saharay Jawaharlal Nehru Center for Advanced Scientific Research, Chemistry and Physics of Materials Unit, Bangalore.

2. Outline • Difference between MD and ab initio MD • Why to use ab initio MD ? • Born-Oppenheimer Molecular Dynamics • Car-Parrinello Molecular Dynamics • Applications of CPMD • Disadvantages of CPMD • Other methods • Conclusions

3. DFT, MD, and CPMD • Properties of liquids/fluids depend a lot on configurational entropy • MD with improved empirical potentials • DFT calculation of a frozen liquid configuration • Configurational Entropy part of the free energy will be missing in that case. • Ab initio MD offers a path that mixes the goodness of both MD and of DFT • AIMD is expensive.

4. Classical MD Hardwired potential No electronic degrees of freedom No chemical reaction Accessible length scale ~100 Å Accessible time scale ~ 10 ns Ab initio MD On-the-fly potential Electronic degrees of freedom Formation and breaking of bonds Accessible length scale ~ 20 Å Accessible time scale ~ 10 ps Molecular simulations

5. A short intense shock caused the hydrogen to form a hot plasma and become a conducting metal Livermore’s Nova Laser Sandia National Laboratories Z accelerators The experiments found different compressibilities which could affect the equation of state of hydrogen and its isotope Quantum simulations could give the proper reasons for different results Conditions of the Nova and Z flyer were different : Time scales of the pulse were different

6. Why ab Initio MD ? • Chemical processes • Poorly known inter atomic interactions e.g. at high Pressure and/or Temperature • Properties depending explicitly on electronic states ; IR spectra, Raman scattering, and NMR chemical shift • Bonding properties of complex systems

7. Born-Oppenheimer approximation • Electronic motion and nuclear motion can be separated due to huge difference in mass • Different time scale for electronic and ionic motion • Fast electrons have enough time to readjust and follow the slow ions

8. Born-Oppenheimer MD • Electron quantum adiabatic evolution and classical ionic dynamics Effective Hamiltonian : HoI→ Ionic k.e. and ion-ion interaction 2nd term → Free energy of an inhomogeneous electron gas in the presence of fixed ions at positions (RI) Electronic ground state – electron density ρ(r) – F({RI}) min Born-Oppenheimer Potential Energy Surface

9. Minimization to BO potential surface E{ρ(r)} ρ0(r) ρ(r)

10. Born-Oppenheimer MD Forces on the ions due to electrons in ground state Ionic Potential Energy Ψi (r) one particle electron wave function 1st→ Electronic k.e. ; 2nd → Electrostatic Hartree term 3rd → integral of LDA exchange and correlation energy density εxc 4th→ Electron-Ion pseudopotential interaction ; 5th → Ion-Ion interaction

11. Born-Oppenheimer MD Electronic density ; fi→ occupation number EeI→ Electron-Ion coupling term includes local and nonlocal components Kohn-Sham Hamiltonian operator Time evolution of electronic variables Time dependence of Hks← slow ionic evolution given by Newton’s equations - Uks = minimum of Eks w.r.t. ψi

12. Merits and Demerits of BOMD Advantages Disadvantages • True electronic Adiabatic Evolution on the BO surface • Need to solve the self- consistent electronic-structure problem at each time step • Minimization algorithms require ~ 10 iterations to converge to the BO forces • Poorly converged electronic minimization → damping of the ionic motion Computationally demanding procedure

13. Car-Parrinello MD CP Lagrangian Ψi → classical fields μ → mass like parameter [1 Hartree x 1 atu 2 ] 4th→ orthonormality of the wavefunctions Constraints on the KS orbitals are holonomic No dissipation

14. If the gap between occupied and unoccupied states = 0 Choice of μ Folkmar Bornemann and Christof Schutte demonstrate (Insulators and semiconductors) μ must be small → small integration time step μ ~ 400 au , time step ~ 0.096x 10-15 s If the gap between occupied and unoccupied states = 0 (Metals) Fictitious kinetic energy of the electrons grow without control Use electronic thermostat

15. CP Equations of motion Equations of motion from Lcp : Electronic time evolution Ionic time evolution Constraint equation Boundary conditions

16. Hellmann-Feynman Theorem If Ψ is an exact eigenfunction of a Hamiltonian H, and E is the corresponding energy eigenvalue : λ is any parameter occurring in H For an approximate wavefunction Ψ For an exact Ψ

17. GI→ constraint force Force on Ions + GI When, ψi is an eigenfunction Force on the ions due to electronic configuration, when electronic wavefunction is an eigen function is zero

18. Constants of motion Vibrational density of states of electronic degrees of freedom Comparison with the highest frequency phonon mode of nuclear subsystem

19. Constants of motion

20. Advantages Fast dynamics compared to BOMD No need to perform the quenching of electronic wave function at each time step Disadvantages Dynamics is different from the adiabatic evolution on BO surface Forces on ions are different from the BO forces Ground state Merits and Demerits of CPMD Ψi ≡ Ψksi → good agreement with the BOMD

21. Velocity Verlet algorithm for CPMD .

22. References • R. Car and M. Parrinello; Phys. Rev. Lett. 55 (22), 2471 (1985) • D. Marx, J. Hutter; http://www.fz-juelich.de/nic-series/ • F. Buda et. al; Phys. Rev. A 44 (10), 6334 (1991) • D.K. Remler, P.A. Madden; Mol. Phys. 70 (6), 921 (1990) • B.M. Deb; Rev. Mod. Phys. 45 (1), 22 (1973) • M. Parrinello; Comp. Chemistry 22, (2000) • M.C. Payne et. al; Rev. Mod. Phys. 64 (4), 1045 (1992)

23. CPMD CPMD code is available at http://www.cpmd.org Michele Parrinello, Jurg Hutter, D. Marx, P. Focher, M. Tuckerman, W. Andreoni, A. Curioni, E. Fois, U. Roethlisberger, P. Giannozzi, T. Deutsch, A. Alavi, D. Sebastiani, A. Laio, J. VandeVondele, A. Seitsonen, S. Billeter and others Code developers : PWscf (Plane Wave Self Consistent field) http://www.pwscf.org http://homepages.nyu.edu/~mt33/PINY_MD/PINY.html PINY-MD

24. Applications

25. Autoionization in Liquid Water pH = - log [H+] pH determination of water by CPMD Intact water molecules dissociate → OH- + H3O+ Rare event ~ 10 hours >>>> fs Transition state separation between the charges ~ 6Å Proposed theory→ Autoionization occurs due to specific solvent structure and hydrogen bond pattern at transition state Diffusion of ions from this transition state Microsecond motion of a system as it crosses transition state can not be resolved experimentally Role of solvent structure in autoionization Diffusion of ions Chandler, Parrinello et. al Science2001, 291, 2121

26. + + Nature of proton transfer in water Grotthuss’s idea :Proton has very high mobility in liquid water which is due to the rearrangement of bonds through a long chain of water molecule; effective motion of proton than the real movement

27. 2 3 4 5 6 Charge separation Dissociation: Fluctuation in solvent electric field ; cleavage of OH bond H3O+ moves by proton transfer within 30 fs 1 Crucial fluctuations carries system to transition state ; breaking of H-bond : 30 fs NO fast ion recombination Conduction of proton through H-bond network 60 fs Chandler, Parrinello et. al Science2001, 291, 2121

28. Order parameter for autoionization Fluctuations that control routes for proton : No. of hydrogen bond connecting the ions : ℓ ℓ = 2 ; recombination occurs within 100 fsreactant ℓ = 0 ; product ℓ ≥ 3 Critical ion separation is 6 Å At ℓ = 2 , sometimes reactant basin ; Thus ℓ is not the only order parameter Potential of proton in H-bonded wire → fluctuation q → configuration description ; q = 1 neutral ; q = 0 charge separated ΔE = E[r(1) – r(0)] → solvent preference for separated ions over neutral molecules

29. Potential of protons in hydrogen bonded wires connecting H3O+ and OH- ions Electric field starts to appear ; metastable state w.r.t. proton motion ; 2kcal/mol higher than neutral state Neutral state, bond destabilizing electric field has not appeared Field fluctuations increase ; stable charge separated state ; 20kcal/mol more stable Chandler, Parrinello et. al Science2001, 291, 2121

30. H9O4+ H5O2+ + + + Nature of the hydrated excess proton in liquid water Two proposed theories : 1. Formation of H9O4+ (by Eigen) 2. Formation of H5O2+ (by Zundel) Ab initio calculations show that transport of H+ and OH- are significantly different Charge migration happens in a few picoseconds Hydrogen bonds in solvation shells of the ions break and reform and the local environment reorders Tuckermann, Parrinello et. al J. Chem. Phys.1997, 275, 817

31. Proton transport Proton diffusion does not occur via hydrodynamic Stokes diffusion of a rigid complex Continual interconversion between the covalent and hydrogen bonds Tuckermann, Parrinello et. al Nature 1997, 275, 817

32. Ob H Oa + Proton transport δ = ROaH - RObH For small δ ; equal sharing of excess proton → Zundel’s H5O2+ For large δ ; threefold coordinated H3O+ → Eigen’s H9O4+ Free energy : ΔF(ν) = -kBT ln [ ∫ dROO P(ROO,ν) ] H5O2+ : at δ = 0 ± 0.05Å, Roo ~ 2.46-2.48 Å ΔF < 0.15 kcal/mol, thermal energy = 0.59 kcal/mol Numerous unclassified situations exists in between these two limiting structures Tuckermann, Parrinello et. al Nature 1997, 275, 817

33. Breaking bonds by mechanical stress Reactions induced by mechanical stress in PEG 1. Formation of ions corresponds to heterolytic bond cleavage 2. Motion of electrons during the reaction Polymer is expanded with AFM tip Unconstrained reactions can not be observed by classical MD Quantum chemical approaches are more powerful in describing the general chemical reactivity of complex systems Frank et. al J. Am. Chem. Soc. 2002, 124, 3402

34. H O H Breaking bonds by mechanical stress Method ΔE (C-O) kcal/mol ΔE (C-C) kcal/mol Radicaloid bond breaking BLYP Solvent 79.1 83.9 H H H H O1 Exp H C C 85.0 83.0 O2 C C2 H H H H H Small piece of PEG in water After equilibration, distance between O1 and O2 was increased continuously by 0.0001 au/time Reaction started at 250 K ; C2O1 ~ 3.2 Å

35. H - O H H H H H H H H H H O O O O O O O O O O O O O H H H H H H H H H + H H H O O O H H H O H H H H O - O H O O O H O O O H H H O H O O + O Snapshots of the reaction mechanisms 250 K 320 K Frank et. al J. Am. Chem. Soc. 2002, 124, 3402

36. Hydrogen bond driven chemical reaction Beckmann rearrangement of Cyclohexanone Oxime into ε-Caprolactam in SCW SCW accelerates and make selective synthetic organic reactions System description : • CPMD simulation , BLYP exchange correlation • MT norm conserving pseudo potential • Plane wave cut-off 70 Ry, Nose-Hoover thermostat • T = 673K, 300K • 64 H2O + 1 solute, 18 ps analysis + 11 ps equil. Disrupted hydrogen bond network of SCW alters the solvation of O and N Parrinello et. al J. Am. Chem. Soc.2004, 126, 6280

37. Proton attack on the Cyclohexanone Oxime Parrinello et. al J. Am. Chem. Soc.2004, 126, 6280

38. Problems • Computationally costly • Can not simulate slow chemical processes that take place beyond time scales of 10 ps • Inaccuracy in the assumption of exchange and correlation potential Limitation in the number of atoms and time scale of simulation Inaccurate van der Waals forces, height of the transition energy barrier • BOMD not applicable for photochemistry; transition between different electronic energy levels

39. Other methods • QM/MM – quantum mechanics / molecular mechanics e.g. catalytic part in enzyme AIMD Classical MD • Path-sampling approach combined with ab-initioMD for slow chemical processes • Metadynamics, for slow processes

40. Conclusions • CPMD : nuclear and electronic degrees of freedom • Interaction potential is evaluated on-the-fly • Bond formation and breaking is accessible in CPMD : direct access to the chemistry of materials • Transferability over different phases of matter • CPMD is computationally expensive

41. Acknowledgement Prof. S. Balasubramanian Dr. M. Krishnan, Bhargava, Sheeba, Saswati

42. THANK YOU