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

1
Introduction into Molecular DynamicsRalf
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
2
1. Molecular Dynamics

• Solve Newtons equation for a molecular system

• Or equivalently solve the classical Hamiltonian
  equation

3
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

4
1. Molecular Dynamics method

• K. Nordlund
• U. Helsinki

5
1. Molecular Dynamics method

• K. Nordlund, U. Helsinki

6
1. Molecular Dynamics method

• K. Nordlund
• U. Helsinki

7
1. Molecular Dynamics method

• RuBisCO protein simulations Paul Crozier,
  Sandia

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.
8
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

9
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

10
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

11
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.
12
1. MD dynamics
13
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)

14
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

15
2. Molecular Dynamics Potentials

• Potential energy includes terms for (contd.)
• Electrostatic
• Interactions
• van der Waals
• Interactions

16
2. Molecular Interaction Types Non-bonded
Energy Terms

• Lennard-Jones Energy.

• Coloumb Energy.

17
2. Molecular Interaction Types Bonded Energy
Terms

• Bond energy

• Bond Angle Energy

18
2. Molecular Interaction Types Bonded Energy
Terms

• Improper Dihedral Angle Energy

• Dihedral Angle Energy

19
2. Scaling

• Scaling by model parameters
• size s
• energy e
• mass m

taken from Dr. D. A. Kofkes lectures on
Molecular Simulation, SUNY Buffalo http//www.eng
.buffalo.edu/kofke/ce530/index.html
20
2. L-J dimensionless form

• Dimensions and units - scaling

• Lennard-Jones potential in dimensionless form
• 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. Kofkes lectures on
Molecular Simulation, SUNY Buffalo http//www.eng.
buffalo.edu/kofke/ce530/index.html
21
3. Historical Perspective on MD
22
3. First molecular dynamics simulation (1957/59)
Hard disks and spheres (calculation of
collision times)
IBM-704
solid phase
liquid phase
liquid-vapour-phase
N32 ? 6.5x105 coll. ? 4 days N500 ? 107
coll. ? 2.3 years
N32 7000 collisions / h N500 500
collisions / h
Production run 20000 steps
23
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)
24
3. MD development
(aus T. Schlick, Molecular Modelling and
Simulation, 2002)
25
4. Verlet algorithm
then
Let
Starting from
and
all subsequent positions are determined
For the kinetic energy we need the velocities

• Note the velcoities are one step behind.
  Therefore
• Specify positions and
• Compute the forces at timestep n
• 3. Compute the positions at timestep (n1) as in
  (1.1)
• 4. Compute velocities at timestep n as in (1.2)
  then increment n and goto 2.

26
4. A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt

v
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
27
4. A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt
• Then determine positions at the next time step
• X(tDt) X(t) v(t Dt/2)Dt

v
X
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
28
4. A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt
• Then determine positions at the next time step
• X(tDt) X(t) v(t Dt/2)Dt

v
v
X
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
29
4. Verlet algorithm- velocity form
30
4. Advantages

• Positions and velocities are now in step
• gt kinetic and potential energies are in step
• Numerical stability is enhanced
• Eq (1.2) gives velocity as difference of 2 rs of
  the same order of magnitude gt round-off errors
• important for long runs
• With reasonable h, Verlets algorithm conserves
  energy

31
4. Beeman algorithm

• Better velocities, better energy conservation
• More expensive to calculate

32
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!
33
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

34
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

35
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
• This can be obtained by assigning the velocity
  components vi from a random Gaussian distribution
• with mean 0 and standard deviation (kBT/mi)

36
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

37
4. Minimum Image
These two are same distance from central atom,
yet Black atom interactsblue atom does not

• 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

Only interactions considered
38
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)
Reducing the computing time use of cut-off in
UNB The cutoff distance may be no greater than ½
L (L box length)
39
4. Potential truncation
common approach cut-off the at a fixed value
Rcut problem discontinuity in energy and force
possibility of large errors
40
4. Speed-up
Tamar Schlick, Molecular Modeling and
Simulation, Springer
41
4. Cutoff schemes for faster energy computation
wij weights (0lt wij lt1). Can be used to exclude
bonded terms, or to scale some interactions
(usually 1-4)
S(r) cut-off function. Three types 1)
Truncation
b
42
4. Cutoff schemes for faster energy computation
2. Switching
a
b
with
3. Shifting
or
b
43
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.

44
4. Simulating at constant T the Berendsen
scheme

• Bath supplies or removes heat from the system as
  appropriate
• Exponentially scale the velocities at each time
  step by the factor ?
• where ? determines how strong the bath influences
  the system

Heat bath
T kinetic temperature
Berendsen et al. Molecular dynamics with coupling
to an external bath. J. Chem. Phys. 813684 (1984)
45
4. Simulating at constant P Berendsen scheme

• Couple the system to a pressure bath
• Exponentially scale the volume of the simulation
  box at each time step by a factor ?
• where ? isothermal compressibility
• ?P coupling constant

pressure bath
where
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. 813684 (1984)
46
4. MD scheme
47
5. Analysis of MD
Configurations Averages Fluctuations Time
Correlations
48
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
49
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

50
5. Statistical Mechanics Extracting properties
from simulations

• static properties such as structure, energy, and
  pressure are obtained from pair (radial)
  distribution functions

DME-water and DMP-water solutions
51
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) -gt 0 as r -gt 0 i.e. atoms are inpenetrable
• g(r) tends to 1 as r goes to infinity

52
5. Pair correlation function
53
5. Outputs
54
5. Equilibration of energy
55
5. Time variation of energies

• kinetic energies
• potential energies

56
5. Time variation of pressure

• Equilibration of pressure with time

57
5. Statistical Mechanics Extracting properties
from simulations

• dynamic and transport properties are obtained
  from time correlation functions

ó
ô
dt
T

lt
A

(t)
A

(0)gt
ô
õÿ
0
ltA(t) - A(0)2gt

lim


/2t





t
Polybutadiene, 353 K Mw 1600
58
5. Outputs
- Hydrogen in perfect crystal graphite
150K
900K
59
5. Outputs
two diffusion channels
no diffusion across graphene layers (150K 900K)
Lévy flights?
60
5. Outputs
Non-Arrhenius temperature dependence
61
7. Bottlenecks in Molecular Dynamics

• Long-range electrostatic interactions O(N2) fast
  electrostatics algorithms (each method still
  needs fine-tuning 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 mat-vec) 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

62
7. Ewald Sum Method
63
7. Ewald Sum Method
64
7. Ewald Sum Method

• additional corrections
• arises from a gaussian acting on its own site
• (self-energy correction)
• or from a surface in vacuum

65
7. Particle Mesh Ewald

• 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

66
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)
67
7. Multigrid
68
7. Multiple time step dynamics

• Reversible reference system propagation algorithm
  (r-RESPA)
• 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

69
7. Multiple Time Step Algorithm
Reversible Reference System Propagator Algorithm
(r-RESPA)
Tuckerman, Berne, Martyna, 1992
The Liouville Operator
r(tDt), v(tDt)
r(t), v(t)
The Liouville Propagator and the state of system
is given by
Trotter expansion of the Liouville Propagator
Reference System Propagator dt Correction
Propagator Dt n dt
70
7. RESPA for Biosystems
5-stage r-RESPA decomposition for biological
systems
Zhou Berne, 1997
The Liouville Operator decomposition for
biosystems
Reference Correction
time step
0.50 fs 1.0 fs 2.0 fs 4.0 fs 8.0 fs
71
7. FMM/RESPA Performance
72
7. P3ME/RESPA Performance
73
7. P3ME/RESPA Performance
74
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

75
Multi-scales
sputtered and backscattered species and fluxes
Plasma-wall interaction
Molecular dynamics
Kinetic Monte Carlo
Kinetic model
Fluid model
Binary collision approximation
impinging particle and energy fluxes
76
From atoms to W7-X
77
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, Garching)
78
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)
79
5. Outputs
- Hydrogen in perfect crystal graphite
150K
900K
80
5. Outputs
two diffusion channels
no diffusion across graphene layers (150K 900K)
Lévy flights?
81
5. Outputs
Non-Arrhenius temperature dependence
82
8. Multi-scale approach
Macroscales KMC and Monte Carlo Diffusion (MCD)
Mesoscales Kinetic Monte Carlo (KMC)
Microscales Molecular Dynamics (MD)
83
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

84
Multi-scale results
standard graphites
highly saturated graphite
85
Effect of voids
B 20 voids
C 20 voids
A 10 voids
Larger voids
Higher diffusion
Longer jumps
86

TRIM, TRIDYN much faster than MD (simplified
physics binary collision approximation)

• very good match of physical sputtering
• dynamical changes of surface composition

87
2D-TRIDYN
88
2D-TRIDYN
89
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)

90
Max-Planck-Institut für Plasmaphysik, EURATOM
Association
Extension of TRIDYN
91
Max-Planck-Institut für Plasmaphysik, EURATOM
Association
Extension of TRIDYN
92
Max-Planck-Institut für Plasmaphysik, EURATOM
Association
Extension of TRIDYN