9th November 2010, Bratislava
LCC Laboratory of Computational
Chemistry
Petr Kulhánek
kulhanek@chemi.muni.cz
National Centre for Biomolecular Research, Masaryk University
Faculty of Science, Kotlářská 2
CZ-61137 Brno
Metadynamics
&
CPMD
-1-
Contents
 Free Energy Calculations
 Sampling Problem
 Potential of Mean Force Methods
 Metadynamics
 Constrained Dynamics
 Adaptive Biasing Force
 Multiple Walkers Approach
 Molecular Dynamics and Reactions
 Empirical Valence Bond
 CPMD versus BOMD
9th November 2010, Bratislava -2-
Free Energy Calculations
9th November 2010, Bratislava -3-
Free Energy
9th November 2010, Bratislava -4-
KRTAŕ
ln
RT
A
B
e
h
Tk
k1
r
A
A
Free energy is related to
equilibrium and rate constants.
Knowledge of free energy allows to quantify:
• chemical reactivity (e.g. enzymatic activities)
• thermodynamics (e.g. binding affinities)
Free Energy Calculations
9th November 2010, Bratislava -5-
1
2
21
lnRTA
0
)(ln)( ARTA
Density of state (probability)
It can be calculated from molecular
dynamics or Monte Carlo
simulations, but …
reaction coordinate
collective variables
Sampling Problem
9th November 2010, Bratislava -6-
d1
10 ns long simulation is able to discover
free energy landscape with depth only
about 3 kcal/mol.
Free Energy Calculations
9th November 2010, Bratislava -7Available
methods:
constrained dynamics
system is biased by constraining reaction coordinate
adaptive biasing force
system is biased by force which is opposite to potential of mean force
umbrella sampling
system is biased by restraining reaction coordinate
metadynamics
system is biased by Gaussian hills, which fill free energy landscape
A system has to be biased achieving efficient sampling in the region of interest. We
need to know how to obtain the unbiased free energy from such biased simulation.
Free Energy Calculations
9th November 2010, Bratislava -8Alchemical
Transformation
one system is slowly changed to another one (changes are very often unrealistic,
atoms are created and/or annihilated)
what: mostly changes in binding free energies:
how: thermodynamic integration (TI), free energy perturbation (FEP)
 Potential of Mean Force
system is changed along reaction coordinate
what: free energy of conformation changes, reaction free energies
how: constrained dynamics, adaptive biasing force, umbrella sampling,
metadynamics, steered dynamics
 End-points Methods
free energy of every state is calculated independently
what: mostly binding free energies
how: MM/XXSA; XX=PB, GB, LRA
Metadynamics
9th November 2010, Bratislava -9Implemented
in CPMD
Metadynamics, theory
9th November 2010, Bratislava -10-
i
i
i
x
xV
=
t
tx
m 2
2
Free energy landscape is filled by Gaussian hills.
ix,V+xV
x
=
t
tx
m h
i
i
i 2
2
Equations of motion MTD history potential
Equations of MTD motion (direct approach)
2
σ
ss
H=is,V
t
i
=t
th
2
exp
2
1
History-dependent term converges to FES
si,V=sA h
i
lim
Metadynamics, example
9th November 2010, Bratislava -11DIS
(distance)
hill height 0.01 kcal/mol, width 0.5 x 0.5 Å
MTD frequency 500 fs
2 ns long simulation
300 K, vacuum, GAFF force field, time step 0.5 fs
Metadynamics, example
9th November 2010, Bratislava -12DIS
(distance)
hill height 0.01 kcal/mol, width 0.5 x 0.5 Å
MTD frequency 500 fs
2 ns long simulation
300 K, vacuum, GAFF force field, time step 0.5 fs
Constrained Dynamics
9th November 2010, Bratislava -13Implemented
in CPMD
Constrained Dynamics, theory
9th November 2010, Bratislava -14-
i
i
i
x
xV
=
t
tx
m 2
2
Reaction coordinate is fixed (constrained) at the value of interest.
xσtλ+xV
x
=
t
tx
m k
k
k
i
i
i 2
2
00
=ξxξ=xσ
Equations of motion Constraint condition
method of Lagrange multipliers
Equations of constrained motion
= Lagrange multipliers
holonomic constraint
k
λ
Constrained Dynamics, theory
9th November 2010, Bratislava -15d
ξ
ξ
d ξ
ξdA
=ΔG
ξ
uc2
1
Derivative of unbiased free energy is also given by (concise formulation):
0
2/1
0
ln
ξξ
ucccuc
Z
d ξ
d
RTλ=
d ξ
ξdA
+
d ξ
ξdA
=
d ξ
ξdA
Final free energy is obtained by numerical integration:
second derivatives are not required
Constrained Dynamics, example
9th November 2010, Bratislava -16Small
numbers are calculated by averaging big numbers.
70 points
d177 points, 0.1 Å
Method B: 5 ps shift, 5 ps equilibration, 20 ps production
300 K, vacuum, GAFF force field, time step 1 fs / 0.5 fs
DIS (distance)
Constrained Dynamics, example
9th November 2010, Bratislava -17d177
points, 0.1 Å
Method B: 5 ps shift, 5 ps equilibration, 20 ps production
300 K, vacuum, GAFF force field, time step 1 fs / 0.5 fs
DIS (distance)
Adaptive Biasing Force
9th November 2010, Bratislava -18Implemented
in CPMD
ABF, theory
9th November 2010, Bratislava -19-
i
i
i
x
xV
=
t
tx
m 2
2
Movement along reaction coordinate is the subject of diffusion process.
Equations of motion Free energy and force along RC
Equations of ABF motion
ξF+
x
xV
=
t
tx
m ABF
i
i
i 2
2
force along reaction coordinate is subtracted
from the system
dx
dξ
dξ
ξdA
=xF ABF
ABF, theory
9th November 2010, Bratislava -20-
ξ
ξ
dt
d ξ
Zdt
d
=
ξ
A 1
Free energy is given by:
it contains the second derivatives of reaction
coordinate if treated analytically
equation is solved numerically
Procedure:
• range of reaction coordinate is divided into bins
• a value of reaction coordinate determine a bin
• a contribution to the derivative of free energy is accumulated
into a bin
• ABF force calculated from accumulated free energy derivative
is applied to the system
• accumulated free energy derivative very rapidly converges
ABF, example
9th November 2010, Bratislava -21-
DD
(difference of
distances)
range from -6.5 to 6.5 Å, 130 bins
2 ns long simulation
300 K, vacuum, GAFF force field, time step 1 fs
Show movie?
Multiple Walkers Approach
9th November 2010, Bratislava -22-
Server
Client 1 Client 2 Client 3
server collects information about free
energy surface (FES) and redistributes it
among clients
Applicable to:
Metadynamics
Adaptive biasing force
Advantages:
 Faster convergence
 Easy to implement
 Parallel scaling is almost linear
FES data are exchanged every Nupdate MD steps
Nupdate ~ 500-1000 steps
Multiple Walkers Approach
9th November 2010, Bratislava -23Two
walkers
(from reactants and products)
One walker
(from reactants)
Nucleophilic substitution reaction (test case)
PMFLib
A Toolkit for Free Energy Calculations
developed by
Petr Kulhánek
9th November 2010, Bratislava -24-
PMFLib, functionality
9th November 2010, Bratislava -25Implemented
methods:
 Constrained dynamics (BM)
 Adaptive biasing force (ABF)
 Metadynamics (MTD)
 Umbrella sampling
 Multiple walker MTD
 Multiple walker ABF
 Replica Exchange Dynamics
 String Method
Reaction coordinates:
 DIS, DS
 DD, DDS
 ODISM, ODSM
 ANG, ANGM, CANG, CANGM
 DIH, DIHM
 AC, GC
 RGYR
 RMSDT, RMSDL
 EPOT, EGAP and variants
Tools:
 MTD energy
 BM integration
 ABF integration
 Multiple walker MTD – server and administration client
 Multiple walker ABF – server and administration client
PMFLib, design
9th November 2010, Bratislava -26-
HiPoLySciMaFicPRMFile
tools
cpmffpmf
PMFLib
drivers
sander pmemdxdynbp
PMFLib
• Fortran 90 (162 files / ~25000 lines)
• C/C++ (178 files / ~24000 lines)
cpmd libatoms
Applications
MM
QM/MM
EVB
MM
QM/MM
semiempirical
QM
DFT
QM/MM
DFT
general
LOTF
theory precision
NetLib
Molecular Dynamics
and
Reactions
9th November 2010, Bratislava -27-
Description of Chemical Reactions
9th November 2010, Bratislava -28EVB
X QM/MM X QM
theory precision
complexity, sampling problems
too complex for
biomolecular systems
problem with boundaries
between QM and MM parts
requires
parametrizations
QM
MM QM
products
reactants
EVB
CPMD
Car-Parrinello
Molecular Dynamics
9th November 2010, Bratislava -29-
Method
9th November 2010, Bratislava -30Equations
of motion:
fictitious mass of wavefunction
(ca 300-700 a.u. , typical value is about 600 a.u.
ions
wavefunction
constraints due to
orthonormality of wavefunction
CPMD versus BOMD
9th November 2010, Bratislava -31-
CPMD
• no SCF procedure
• motion of ions in time
• motions of wavefunction in time
• time step ~ 0.1 fs (5 a.u.)
• DFT only (in CPMD)
• hybrid functional possible but very slow
• dispersion correction available
• planewaves wavefunction (periodicity!}
• wavefunction quality is determined by cutoff (single value)
• pseudopotentials required (core electrons)
BOMD
• SCF procedure
• motion of ions in time only
• wavefunction follows ions by SCF (BO)
• time step max ~ 1 fs
• gradients require very tight convergence of
wavefunction optimization
Practicals …
9th November 2010, Bratislava -32•
read manual (it was very improved in the last version) !
• read two chapters from NIC books about CPMD (about 150 pages), freely downloadable
• be veeeery patient
Typical setup:
• time step 5 a.u.
• WF mass 600 a.u.
• pseudopotentials: Troulier-Martins normconserving
• WF cutoff: 70 Rydbergs
• charged/isolated systems (also in QM/MM calculations)
• add 2-3A to box dimensions, molecule has to be centered!!!
• Poission solver: Tuckerman
• heating: Berendsen thermostat
• production: Nose-Hoover thermostat
• ultrasoft Vanderbilt pseudopotential are problematic for Mg
Practicals in smaller group.
Acknowledgements
 Igor Tvaroška, Stanislav Kozmon
 Jaroslav Koča, Zora Střelcová, Jakub Štěpán
 Monika Fuxreiter
 POSTBIOMIN (financial support)
9th November 2010, Bratislava -33-