Transcript Introduction into Molecular Dynamics Ralf Schneider, Abha Rai, Amit Rai Sharma

Introduction into Molecular Dynamics
Ralf Schneider,
Abha Rai, Amit Rai Sharma
Outline:
1. Basics
2. Potentials
3. History
4. Numerics
5. Analysis of MD runs
(6. Physics extensions=
7. Numerical extensions
8. Summary
1. Molecular Dynamics


Solve Newton’s equation for a molecular system:


F  ma
Or equivalently solve the classical Hamiltonian
2
equation:
N
Pi
H (p i , ri )  
 V (ri )
i 1 2mi
H
p i  
 fi
ri
H p i
ri 

p i m
1. Molecular Dynamics method
• deterministic method: state of the system at any
future time can be predicted from its current state
• MD cycle for one step:
1) force acting on each atom is assumed to be
constant during the time interval
2) forces on the atoms are computed and
combined with the current positions and velocities
to generate new positions and velocities a short
time ahead
1. Molecular Dynamics method
K. Nordlund
U. Helsinki
1. Molecular Dynamics method
K. Nordlund, U. Helsinki
1. Molecular Dynamics method
K. Nordlund
U. Helsinki
1. Molecular Dynamics method
RuBisCO protein simulations
Important for
converting CO2 to
organic forms of
carbon and in the
photosynthetic
process.
Even though the
pocket is closed, a
CO2 molecule
escapes, which was
a surprise.
Paul Crozier, Sandia
1. Molecular Dynamics method
H. Wang
U. California,
Santa Cruz
ATP Synthase:
within the mitochondria of a cell a rotary engine uses the potential difference
across the bilipid layer to power a chemical transformation of ADP into ATP
1. Motivation: Why atomistic MD simulations?
 MD simulations provide a molecular level
picture of structure and dynamics (biological
systems!)  property/structure relationships
 Experiments often do not provide the molecular
level information available from simulations
 Simulators and experimentalists can have a
synergistic relationship, leading to new insights
into materials properties
1. Motivation: Why atomistic MD simulations?
MD simulations allow prediction of properties for
 Novel materials which have not been synthesized
 Existing materials whose properties are difficult to
measure or poorly understood
 Model validation
1. Timescales
Bond vibrations - 1 fs
Collective vibrations - 1 ps
Conformational transitions - ps or longer
Enzyme catalysis - microsecond/millisecond
Ligand Binding - micro/millisecond
Protein Folding - millisecond/second
Molecular dynamics:
Integration timestep - 1 femtosecond
Set by fastest varying force.
Accessible timescale about 10 nanoseconds.
1. MD dynamics
2. Molecular Dynamics: Potential
We need to know
The motion of the
atoms in a molecule, x(t)
and therefore,
the potential energy, V(x)
2. Molecular Dynamics: Potentials
How do we describe the potential energy V(x) for a
molecule?
Potential Energy includes terms for
Bond stretching
Angle Bending
Torsional rotation
Improper dihedrals
2. Molecular Dynamics: Potentials
Potential energy includes terms for (contd.)
Electrostatic
Interactions
van der Waals
Interactions
2. Molecular Interaction Types –
Non-bonded Energy Terms
• Lennard-Jones Energy.

Coloumb Energy.
2. Molecular Interaction Types –
Bonded Energy Terms
• Bond energy:

Bond Angle Energy:
2. Molecular Interaction Types –
Bonded Energy Terms
• Improper Dihedral
Angle Energy:

Dihedral Angle
Energy:
2. Scaling
 Scaling by model parameters
 size s
 energy e
 mass m
taken from Dr. D. A. Kofke’s lectures on
Molecular Simulation, SUNY Buffalo
http://www.eng.buffalo.edu/~kofke/ce530/index.html
2. L-J: dimensionless form
• Dimensions and units - scaling
 Lennard-Jones potential in
dimensionless form
 1 12  1 6 
u * (r*)  4     
 r * 
 r * 
 Dimensionless properties must also be
parameter independent
 convenient to report properties in
this form, e.g. P*(r*)
 select model values to get actual
values of properties
 Equivalent to selecting unit value for
parameters
taken from Dr. D. A. Kofke’s lectures on Molecular Simulation,
SUNY Buffalo
http://www.eng.buffalo.edu/~kofke/ce530/index.html
3. Historical Perspective on MD
3. First molecular dynamics simulation (1957/59)
Hard disks and spheres
0
uij (r )  

rd
rd
(calculation of collision times)
IBM-704:
solid phase
N=32: 7000 collisions / h
N=500: 500 collisions / h
liquid phase
liquid-vapour-phase
N=32  6.5x105 coll.  4 days
Production run ~20000 steps
N=500  107 coll.
 2.3 years
3. First MD simulation using continuous potentials (1964)
CDC-3600
VACF
RDF
MSD
864 particles
Time / Step ~ 45s
Production run ~20000 steps  10 days !
(standard PC [Pentium 1.2 GHz]: ½ hour)
3. MD – development
(aus T. Schlick, „Molecular Modelling and Simulation“, 2002)
4. Verlet algorithm
d 2 ri 1
1
 [ri (t  h)  2ri (t )  ri (t  h)]  Fi (t )
m
dt 2 h 2
Let
tn  nh,ri  ri (tn ),Fi n  Fi (tn )
0
Starting from ri
then
ri n1  2ri n  ri n1  Fi n h 2 / m (1.1)
1
and ri all subsequent positions are determined
For the kinetic energy we need the velocities
vi n  (ri n1  ri n1) / 2h
(1.2)
Note: the velcoities are one step behind. Therefore:
0
1
r
r
1. Specify positions i and i
n
2. Compute the forces at timestep n: Fi
n  1,2,...
n1
3. Compute the positions at timestep (n+1) as in (1.1): ri
4. Compute velocities at timestep n as in (1.2); then increment n and goto 2.
4. A widely-used algorithm: Leap-frog Verlet
Using accelerations of the current time step,
compute the velocities at half-time step:
v(t+Dt/2) = v(t – Dt/2) + a(t)Dt
v
t-Dt/2
t
t+Dt/2
t+Dt
t+3Dt/2 t+2Dt
4. A widely-used algorithm: Leap-frog Verlet
Using accelerations of the current time step,
compute the velocities at half-time step:
v(t+Dt/2) = v(t – Dt/2) + a(t)Dt
Then determine positions at the next time step:
X(t+Dt) = X(t) + v(t + Dt/2)Dt
v
t-Dt/2
X
t
t+Dt/2
t+Dt
t+3Dt/2 t+2Dt
4. A widely-used algorithm: Leap-frog Verlet
Using accelerations of the current time step,
compute the velocities at half-time step:
v(t+Dt/2) = v(t – Dt/2) + a(t)Dt
Then determine positions at the next time step:
X(t+Dt) = X(t) + v(t + Dt/2)Dt
v
t-Dt/2
v
X
t
t+Dt/2
t+Dt
t+3Dt/2 t+2Dt
4. Verlet algorithm- velocity form
The previous algorithm is not self starting. One needs two sets
of positions: ri 0 and ri1 . If one has the initial velocities vi 0 use
ri1  ri 0  hvi 0  Fi 0 h 2 / 2 m
and proceed as before.
From eq(1.1) and (1.2) we can derive:
ri n 1  ri n  hvi n  Fi n h 2 / 2 m
(1.3) and
vi n 1  vi n  h( Fi n 1  Fi n ) / 2 m
(1.4)
The velocity form of the Verlet algorithm is:
1. Specify the initial positions and velocities ri1 and vi1 .
2. Compute the positions at time step n  1 via (1.3) n  1,2,.. .
3. Compute the velocities at time step n  1 via (1.4)
4. Increment n and go back to 2.
4. Advantages

Positions and velocities are now in step
=> kinetic and potential energies are in step

Numerical stability is enhanced
Eq (1.2) gives velocity as difference of 2 r’s of the
same order of magnitude => round-off errors
 important for long runs


With reasonable h, Verlet’s algorithm
conserves energy
4. Beeman algorithm
 Better
velocities, better energy conservation
 More expensive to calculate
r(t  t )  r(t )  tv(t )  2 t 2a(t )  1 t 2a(t  t )
3
6
v(t  t )  v(t )  1ta(t )  5 ta(t )  1 ta(t  t )
3
6
6
4. Predictor-corrector algorithms



1. Predictor. From the positions and their time derivatives up to a
certain order q, all known at time t, one ``predicts'' the same
quantities at time by means of a Taylor expansion. Among these
quantities are, of course, accelerations .
2. Force evaluation. The force is computed taking the gradient of
the potential at the predicted positions. The resulting
acceleration will be in general different from the ``predicted
acceleration''. The difference between the two constitutes an
``error signal''.
3. Corrector. This error signal is used to ``correct'' positions and
their derivatives. All the corrections are proportional to the error
signal, the coefficient of proportionality being a ``magic
number'' determined to maximize the stability of the algorithm.
Fifth-order Gear (requires more calculational effort and memory than Verlet, but
needs only one calculation of the force per time step, wheras Verlet needs two!
4. Evaluate integration methods





Fast, minimal memory, easy to program
Calculation of force is time consuming
Conservation of energy and momentum
Time-reversible
Long time step can be used
4. Choosing the time step

Too small: covering small conformation space

Too large: instability

Suggested time steps



Translation, 10 fs
Flexible molecules and rigid bonds, 2fs
Flexible molecules and bonds, 1fs
4. How do you run a MD simulation?
Get the initial configuration
Assign initial velocities
At thermal equilibrium, the expected value of the kinetic energy of
the system at temperature T is:
Ekin
1 3N
1
2
  mi vi  (3N )k BT
2 i 1
2
This can be obtained by assigning the velocity components vi from
a random Gaussian distribution
with mean 0 and standard deviation (kBT/mi):
2
i
v
k BT

mi
4. Periodic Boundary Conditions





infinite system with small number
of particles
remove surface effects
shaded box represents the
system we are simulating, while
the surrounding boxes are exact
copies in every detail
whenever an atom leaves the
simulation cell, it is replaced by
another with exactly the same
velocity, entering from the
opposite cell face (number of
atoms in the cell is conserved)
rcut is the cutoff radius when
calculating the force between two
atoms
4. Minimum Image


Bulk system modeled via
periodic boundary condition
 not feasible to include
interactions with all
images
 must truncate potential
at half the box length (at
most) to have all
separations treated
consistently
Contributions from distant
separations may be
important
These two are same
distance from central
atom, yet:
Black atom interacts
blue atom does not
Only interactions
considered
4. Potential cut-offs
Bonded interactions: local, therefore O(N),
where N is the number of atoms in the molecule considered)
Non-bonded interactions: involve all pairs of
Atoms, therefore O(N2)
U NB
 R 12  R 6 
qi q j
ij
ij
  e ij    2    
 r   i , j nonbonded 4e er
 rij 
i , j nonbonded
0 ij
 ij  

Reducing the computing time: use of cut-off in UNB
The cutoff distance may be no greater than ½ L (L= box length)
4. Potential truncation
common approach:
cut-off the at a fixed value Rcut
problem: discontinuity in energy and force
possibility of large errors
4. Speed-up
Tamar Schlick, “Molecular Modeling and Simulation”, Springer
4. Cutoff schemes for faster energy computation
U NB
 R 12  R 6 
qi q j
ij
ij
  ij S rij e ij    2     ij S rij 
 r   i, j
4e0erij
 rij 
i, j
ij 



ij : weights (0< ij <1). Can be used to exclude bonded terms,
or to scale some interactions (usually 1-4)
S(r) : cut-off function.
Three types:
1) Truncation:
1 r  b
S (r)  
0 r  b
b
4. Cutoff schemes for faster energy computation
2. Switching
ra
 1

S ( r )  1  y ( r )2 2 y ( r )  3 a  r  b
 0
rb

with
a
r2  a2
y(r)  2
b  a2
b
3. Shifting
  r 2 
S1 ( r )  1    
  b  
2
rb
or
 r
S2 ( r )  1  
 b
2
rb
b
4. Neighbor lists
Verlet requires O(N) operations
Force needs O(N2) operations at each step
Most of these are outside range and hence zero
Time reduced by counting only those within range listed in
a table (needs to be updated)
- Verlet (1967) suggested keeping a list of the near
neighbors for a particular molecule, which is updated
periodically
- Between updates of the list, the program does not check
through all the molecules, just those on the list, to calculate
distances and minimum images, check distances against
cutoff, etc.
4. Simulating at constant T:
the Berendsen scheme
system
Heat bath
Bath supplies or removes heat from the system as
appropriate
Exponentially scale the velocities at each time step by the
factor :
Dt  Tbath 
  1  1 

 
T 
T: “kinetic” temperature
where  determines how strong the bath influences the
system
Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)
4. Simulating at constant P:
Berendsen scheme
system
pressure bath
Couple the system to a pressure bath
Exponentially scale the volume of the simulation box at each
time step by a factor :
  1
Dt
P
P  Pbath 
where
N
2

P   Ekin   xi  Fi 
3v 
i 1

where  : isothermal compressibility
P : coupling constant
u : volume
xi: position of particle i
Fi : force on particle i
Berendsen et al. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:3684 (1984)
4. MD scheme
5. Analysis of MD
Configurations
Averages
Fluctuations
Time Correlations
5. Time averages and ensemble averages
macroscopic numbers of atoms or molecules (of the order
of 1023, Avogadro's number is 6.02214199 × 1023 ):
impossible to handle for MD
statistical mechanics (Boltzmann, Gibbs): a single system
evolving in time is replaced by a large number of
replications of the same system that are considered
simultaneously
time average is replaced by an ensemble average:
A    dp N dr N A( p N , r N ) r ( p N , r N )
5. Ergodic hypothesis



Classical statistical mechanics integrates over all
of phase space {r,p}.
The ergodic hypothesis assumes that for sufficiently
long time the phase trajectory of a closed system
passes arbitrarily close to every point in phase
space.
Thus the two averages are equal:  O O
5. Statistical Mechanics
Extracting properties from simulations
 static
4.0
3.5
pair distribution functions
properties such
as structure, energy,
and pressure are
obtained from pair
(radial) distribution
functions
Oether-Hw
3.0
DME
DMP
2.5
2.0
1.5
Oether-Ow
1.0
0.5
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
r (Å)
DME-water and DMP-water solutions
5. Pair correlation function




g(r)dr is the probability of finding a particle in
volume d3r around r given one at r =0
For isotropic system g(r) depends on r only
g(r) -> 0 as r -> 0 i.e. atoms are inpenetrable
g(r) tends to 1 as r goes to infinity
rg(r) 
 (r
j0

 (r0  rj)
rg(r)d d r  N  1
1 is the reference particle
5. Pair correlation function
5. Outputs
Simplest quantity is the Kinetic Energy:
1 t0  t 1
K  lim t  mvi 2 ( t' )dt'
t t 0
i 2
According to the equipartition theorem we have
d
1
2
where d is the dimensionality
< mvi  Nk B T
2
2
This defines our temperature. For the Potential Energy V:
Nr
1 t0  t
d
r
d
)
r
(
g
)
r
(
V
V T  lim t  V ( ri ( t' )  rj ( t' ))dt' 

t t 0
2
i j
where g ( r ), the pair correlation function is defined by:
rg( r ) 
  ( r  ( r0  rj )
j0
and r  N / V is the density
5. Equilibration of energy
5. Time variation of energies

kinetic
energies

potential
energies
5. Time variation of pressure

Equilibration of pressure with time
5. Statistical Mechanics
Extracting properties from simulations
 dynamic
0
2>
<[A(t)
A(0)]
= lim
/2t
t
2
> (Angstrom )
cm
 
T = dt <A(t) A(0)>




2

600
<R
and
transport properties
are obtained from
time correlation
functions
Rouse time
500
400
300
200
slope = 6 X D
100
0
0
5
10
15
20
t (ns)
Polybutadiene, 353 K Mw = 1600
5. Outputs
- Hydrogen in perfect crystal graphite
150K
900K
5. Outputs
two diffusion channels
no diffusion across graphene layers (150K – 900K)
Lévy flights?
5. Outputs
Non-Arrhenius temperature dependence
7. Bottlenecks in Molecular Dynamics

Long-range electrostatic interactions O(N2): fast
electrostatics algorithms (each method still needs finetuning for each system!)
Ewald summation
O(N 3/2 )
Ewald, 1921
Fast Multipole
Method
O(N)
Greengard, 1987
Particle Mesh Ewald
O(N log N)
Darden, 1993
Multi-grid summation
(dense mat-vec as a
sum of sparse matvec)
O(N)
Brandt et al., 1990
Skeel et al., 2002
Izaguirre et al., 2003

Intrinsic different timescales, very small time step needed:
multiple-time step methods
7. Ewald Sum Method
7. Ewald Sum Method
7. Ewald Sum Method
additional corrections:
• arises from a gaussian acting on its own site
(self-energy correction)
• or from a surface in vacuum
7. Particle Mesh Ewald


(1)
(2)
(3)
(4)
(5)
Similar to Ewald method except that it uses FFT
P3ME method has a very similar spirit with PME
Assigning charges onto grids
Use Fast Fourier Transform to speed up the k-space
evaluation
Differentiation to determine forces on the grids
Interpolating the forces on the grid back to particles
Calculating the real-space potential as normal Ewald
7. Fast Multipole Method




Represent charge
distributions in a
hierarchically structured
multipole expansion
Translate distant multipoles
into local electric field
Particles interact with local
fields to count for the
interactions from distant
charges
Short-range interactions are
evaluated pairwise directly
CPU: O(N)
7. Multigrid
7. Multiple time step dynamics

Reversible reference system propagation algorithm (rRESPA)


Forces within a system classified into a number of groups
according to how rapidly the force changes
Each group has its own time step, while maintaining accuracy
and numerical stability
7. Multiple Time Step Algorithm
Reversible Reference System Propagator Algorithm (r-RESPA)
Tuckerman, Berne, Martyna, 1992
The Liouville Operator:
{r(t+Dt), v(t+Dt)}
{r(t), v(t)}
 
  3N  
 



iL    xi
 pi
 Fi ( x) 
    xi

x

p

x
pi 
i 1 
i 1 
i
i
i
3N
The Liouville Propagator and the state of system is given by:
(0)  x(0), p(0)
(Dt )  eiLDt (0)
Trotter expansion of the Liouville Propagator:
L  L1  L2
e
iLDt
e
i ( L1  L2 ) Dt
e
iL2 Dt / 2
Reference System Propagator: t
Correction Propagator:
Dt= n t
e  e
iL1t n
iL2 Dt / 2
7. RESPA for Biosystems
5-stage r-RESPA decomposition for biological systems
Zhou & Berne, 1997
The Liouville Operator decomposition for biosystems:
L  L1  L2  L3  L4  L5


 F1 ( x)
x
p

iL

F
(
x
)
Correction:
i
i
p
F1 ( x)  Fbond ( x)
Reference:
iL1  x
F2 ( x)  Fangle ( x)  Ftors ( x)  FHbond ( x)
time step
0.50 fs
1.0 fs
short
short
F3 ( x)  Fvdw
( x)  Felec
( x)
2.0 fs
med
med
F4 ( x)  Fvdw
( x)  Felec
( x)
4.0 fs
long
long
F5 ( x)  Fvdw
( x)  Felec
( x)
8.0 fs
7. FMM/RESPA Performance
7. P3ME/RESPA Performance
7. P3ME/RESPA Performance
8. MD as a tool

MD is a powerful tool with different levels of
sophistication in terms of physics and numerics

EVERYTHING DEPENDS ON THE
POTENTIAL!

Time-step limitations require combinations with
other methods, e.g. Kinetic Monte Carlo
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Multi-scales
sputtered and backscattered species and fluxes
Plasma-wall interaction
Molecular Binary collision Kinetic
dynamics approximation Monte Carlo
Kinetic
model
Fluid
model
impinging particle and energy fluxes
Max-Planck-Institut für Plasmaphysik, EURATOM Association
From atoms to W7-X
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Diffusion in graphite
Carbon deposition in divertor regions of JET and
ASDEX UPGRADE
Major topics: tritium codeposition
chemical erosion
JET
Paul Coad
(JET)
ASDEX
UPGRADE
Achim von Keudell (IPP,
Garching)
V. Rohde (IPP,
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Diffusion in graphite
Real structure of the material
needs to be included
Internal Structure of Graphite
Granule sizes ~ microns
Void sizes ~ 0.1 microns
Crystallite sizes ~ 50-100 Ångstroms
Micro-void sizes ~ 5-10 Ångstroms
Multi-scale problem in space (1cm to
Ångstroms) and time (pico-seconds to
seconds)
5. Outputs
- Hydrogen in perfect crystal graphite
150K
900K
5. Outputs
two diffusion channels
no diffusion across graphene layers (150K – 900K)
Lévy flights?
5. Outputs
Non-Arrhenius temperature dependence
8. Multi-scale approach
Macroscales
KMC and Monte Carlo
Diffusion (MCD)
Mesoscales
Kinetic Monte Carlo (KMC)
Microscales
Molecular Dynamics (MD)
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Kinetic Monte Carlo
0 = jump attempt frequency (s-1)
Em = migration energy (eV)
T = trapped species temperature (K)
Assumptions:
- Poisson process (assigns real time to the jumps)
- jumps are not correlated
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Multi-scale results
standard
graphites
D(cm2 / s)
highly saturated
graphite
1000 / T (K
1
)
Large variation in observed diffusion coefficients
Diffusion in voids dominates
Strong dependence on void sizes and not void fraction
Saturated H: 0~105s-1 and step sizes ~1Å (QM?)
Diffusion coefficients without knowledge of structure are meaningless
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Effect of voids
A: 10 % voids
Larger voids
B: 20 % voids
Longer jumps
C: 20 % voids
Higher diffusion
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Downgrading
TRIM, TRIDYN: much faster than MD (simplified physics:
binary collision approximation)
- very good match of physical sputtering
- dynamical changes of surface composition
Max-Planck-Institut für Plasmaphysik, EURATOM Association
2D-TRIDYN
Max-Planck-Institut für Plasmaphysik, EURATOM Association
2D-TRIDYN
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Extension of TRIDYN
Bombardment of tungsten with carbon (6 keV)
• steepening of surface structures
(Ivan Biyzukov, IPP Garching)
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Extension of TRIDYN
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Extension of TRIDYN
Max-Planck-Institut für Plasmaphysik, EURATOM Association
Extension of TRIDYN