Transcript Protein Kinase Inhibitors: Insight into Drug Design from

Car-Parrinello Method and
Applications
Moumita Saharay
Jawaharlal Nehru Center for Advanced Scientific Research,
Chemistry and Physics of Materials Unit,
Bangalore.
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
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.
Molecular simulations

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
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
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
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
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
E{ρ(r)}
Minimization to BO potential surface
ρ 0(r)
ρ(r)
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
Born-Oppenheimer MD
; fi → occupation number
Electronic density
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
Merits and Demerits of BOMD
Advantages
Disadvantages
True electronic Adiabatic
Evolution on the BO surface
Need to solve the selfconsistent 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
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
Choice of μ
Folkmar Bornemann and Christof Schutte demonstrate
If the gap between occupied and unoccupied states = 0
(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
CP Equations of motion
Equations of motion from Lcp :
Electronic time evolution
Ionic time evolution
Constraint equation
Boundary conditions
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 Ψ
FI  
E
R I
Force on Ions
+ GI
GI →
constraint force
When, ψi is an eigenfunction
Force on the ions due to electronic configuration, when
electronic wavefunction is an eigen function is zero
Constants of motion
Vibrational density of states of electronic
degrees of freedom
Comparison with the highest frequency
phonon mode of nuclear subsystem
Constants of motion
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
Ground state

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
Velocity Verlet algorithm for CPMD
.
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)
CPMD
CPMD code is available at http://www.cpmd.org
Code developers :
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
PWscf (Plane Wave Self Consistent field) http://www.pwscf.org
PINY-MD http://homepages.nyu.edu/~mt33/PINY_MD/PINY.html
Applications
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 Science 2001, 291, 2121
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
+
+
Charge separation
Dissociation:
Fluctuation in solvent
electric field ;
cleavage of OH bond
1
H3O+ moves by
proton transfer
within 30 fs
2
Crucial
NO fast ion
fluctuations carries recombination
system to transition
state ; breaking of
H-bond : 30 fs
Conduction of
proton through
H-bond network
60 fs
4
3
5
Chandler, Parrinello et. al Science 2001, 291, 2121
6
Order parameter for autoionization
Fluctuations that control routes for proton :
No. of hydrogen bond connecting the ions : ℓ
ℓ = 2 ; recombination occurs within 100 fs
reactant ℓ = 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
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
Nature of the hydrated excess proton in
liquid water
Two proposed theories : 1. Formation of H9O4+ (by Eigen)
2. Formation of H5O2+ (by Zundel)
H9O4+
H5O2+
+
Charge migration happens in a few picoseconds
Ab initio calculations show that
transport of H+ and OH- are
significantly different
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
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
Proton transport
δ = ROaH - RObH
For small δ ; equal sharing of excess proton → Zundel’s H5O2+
For large δ ; threefold coordinated H3O+ → Eigen’s H9O4+
H
+
Ob
Oa
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
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
Breaking bonds by mechanical stress
Method
H
O
H
O1
C
H
H
79.1
H H
H
Exp
C
C
O2
C2
H
Radicaloid bond breaking
ΔE (C-C) kcal/mol
83.9
H
H
BLYP
Solvent
ΔE (C-O) kcal/mol
H
85.0
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 Å
83.0
Snapshots of the reaction mechanisms
250 K
H
H
-O
O
H
O
H
H
O
O
O
H
O
O
O
+
320 K
H
H
H
-O
O
H
H
O
O
H
O
H
O
O
O
+
O
Frank et. al J. Am. Chem. Soc. 2002, 124, 3402
H
O
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
Proton attack on the Cyclohexanone Oxime
Parrinello et. al J. Am. Chem. Soc. 2004, 126, 6280
Problems

Computationally costly
Limitation in the number of atoms and time scale of simulation


Can not simulate slow chemical processes that
take place beyond time scales of 10 ps
Inaccuracy in the assumption of exchange and
correlation potential
Inaccurate van der Waals forces, height of the transition energy
barrier

BOMD not applicable for photochemistry;
transition between different electronic energy levels
Other methods

QM/MM – quantum mechanics / molecular
mechanics
Classical MD


AIMD
e.g. catalytic part in enzyme
Path-sampling approach combined with ab-initio
MD for slow chemical processes
Metadynamics, for slow processes
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
Acknowledgement
Prof. S. Balasubramanian
Dr. M. Krishnan, Bhargava, Sheeba, Saswati
THANK YOU