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
Molecular simulations

• 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

• Classical MD
• Hardwired potential
• No electronic degrees of freedom
• No chemical reaction
• Accessible length scale 100 Å
• Accessible time scale 10 ns

5
A short intense shock caused the hydrogen to form
a hot plasma and become a conducting metal
Livermores 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)
?(r)
? 0(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 exc
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 Newtons 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
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
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
Merits and Demerits of CPMD

• 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

?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-s
  eries/
• 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 gtgtgtgt 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 Science 2001, 291, 2121
26
Nature of proton transfer in water
Grotthusss 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
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 Science 2001, 291, 2121
28
Order parameter for autoionization
Fluctuations that control routes for proton
No. of hydrogen bond connecting the ions l
l 2 recombination occurs within 100
fsreactant l 0 product l 3
Critical ion separation is 6 Å
At l 2 , sometimes reactant basin Thus l 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 Er(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 Science 2001, 291, 2121
30
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
Proton transport
d ROaH - RObH
For small d equal sharing of excess proton ?
Zundels H5O2
For large d threefold coordinated H3O ?
Eigens H9O4
Free energy
?F(?) -kBT ln ? dROO P(ROO,?)
H5O2 at d 0 0.05Å, Roo 2.46-2.48 Å
?F lt 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
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
H
Exp
C
C
85.0
83.0
O2
C2
C
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
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 e-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-initio MD
  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