Transcript Introduction to MPI

Molecular Dynamics
A brief overview
Notes - Websites
• "A Molecular Dynamics Primer", F. Ercolessi
http://www.fisica.uniud.it/~ercolessi/md/
• http://cacs.usc.edu/education/cs596.html
"Scientific computing and visualization", A. Nakano
University of Southern California
• "The Art of Molecular Dynamics Simulation", D. C.
Rapaport, CUP, 1997
• “Computer simulation of liquids”, M. P. Allen, D. J.
Tildesley, OUP, 1990
2
What is molecular dynamics?
• Solving the classical equations of motion
• For a system of N (N>>3) particles
Fi  mi ai
i  1, N
• Which interact through a “given” potential
• And then apply some “tricks” …
• Deterministic technique  Monte Carlo
3
What is molecular dynamics?
•
However, errors in trajectories always accumulate: MD is a statistical mechanics
method  thermodynamic properties
•
To obtain set of configurations according to statistical ensemble
H   E 
-
Microcanonical ensemble (NVE) 
-
Canonical ensemble (NVT)
-
Isobaric-Isothermal Ensemble (NPT)
Grand canonical ensemble (also number of particles can change, mVT)
exp H  / k BT 
•
Ergodic hypothesis
•
Also used for the optimization of structures (simulated annealing)
4
Applications of molecular dynamics
•
•
•
•
•
•
•
•
•
•
Properties of liquids
Plasma physics
Defects in solids
Fracture
Surface properties
Friction
Molecular clusters
Biomolecules
Dynamics of galaxies
Formation of stellar clusters
5
Overview I
• Model
- System Hamiltonian
- Interaction potentials: bonded and non-bonded interactions
- Finite system – infinite system
• Integrator
Symplecticity?
• Statistical ensemble
• Collecting results
6
Overview II
Collecting data
7
First MD simulation (1957/1959)
8
weeks
First MD simulation using
continuous potentials (1960)
IBM 704
9
First MD simulation using
Lennard-Jones potential (1964)
  12   6 
u  r   4      
 r   r  
10
Limitations
• Use of classical forces  quantum effects
- MD only valid if
(at 300 K, t > 6 ps)
• Realism of forces
• Time and size limitations
- Thousands to millions of atoms
- Time step t should be as large as possible while conserving total energy
- In general, t ≈0.01 x fastest behavior of your system
Atoms oscillate about once every 10-12 s in a solid  t ≈10-14 s
- Total time: picoseconds to hundreds of nanoseconds
- Simulation only reliable if simulation time is much longer than relaxation
time of quantities of interest
- Idem for correlation length
11
Initialisation
• Positions
- Random positions
- Regular pattern, e.g. fcc lattice
- Previous run
• Velocities
- Random velocity or from Maxwell distribution
- Previous run
- Linked with temperature
- No drift condition
- Rescale velocities to realize desired temperature
12
Interaction potentials
• Origins are quantum mechanical
• Easiest model:
  12   6 
Lennard-Jones potential u  r   4      
 r   r  
• Truncated (but with energy conservation)

du
r  rc
u (r )  u (rc )  (r  rc )
dr rc
u 'r   
r  rc

0

13
Interaction potentials
• Distinguish between
long range  short range interactions
• Distinguish between
intermolecular  intramolecular forces
14
Interaction potentials
• Stretch energy
• Bending energy
15
Interaction potentials
• Interactions between charge inhomogeneities
Approaches
- Point charges
- Point multipoles
• Screened Coulomb interaction
16
Interaction potentials
• Multi-body interactions
e.g. Tersoff and Brenner potentials
17
Infinite systems
• Periodic boundary conditions
Number of interacting pairs
increases enormously
rc
• Minimum image criterion for
short range potentials
- At most one among all pairs
formed by a particle i in the box
and the set of all periodic
images of another particle j will Central
interact
Simulation
box
Lx , Ly  2rc
18
Infinite systems
• Ewald method for Coulomb interaction
19
Integrators
• How should a good integration scheme look like?
-
High accuracy (reproduces true trajectories well)
Good stability (conservation of energy)
Time reversible
Robust (allow for large time steps)
Conservation of phase space density
(Liouville’s theorem)(symplecticity)
• Simple Euler method is not time reversible and not
symplectic.
20
Integrators
1 2
t a(t )  O(t 3 )
2
1 2
r
(
t


t
)

r
(
t
)

v
(
t
)

t

t a(t )  O(t 3 )
+
2
• Verlet algorithm
- Positions
r (t  t )  r (t )  v(t )t 
r (t  t )  2r (t )  r (t  t )  t 2 a(t )  O(t 4 )
- Velocities
v(t )  r(t ) 
1
(r (t  t )  r (t  t ))  O(t 2 )
2t
- Properties
•
•
•
•
Time reversible
Symplectic
Does not suffer from energy drift
But no info on velocity untill the next step is made
21
Integrators
• Comparison between Euler and Verlet
Test system consists of 7 Lennard-Jones atoms (Ar)
Time step
is 10 fs
Time step
is 1 fs
22
Integrators
• Velocity Verlet
ri  t  t   ri  t   vi  t  t  1/ 2  a (t )t 2
vi  t  t / 2   vi  t   1/ 2  a (t )t

ai  t  t    1/ m  V r N  t  t 

vi  t  t   vi  t  t / 2   1/ 2  a (t  t )t



vi (t  t / 2)  vi (t )  (1 / 2)ai (t )t
Velocity calculated explicitly

Possible to control the temperature


T0 
Stable
v
(
t


t
/
2
)

v
(
t
)

(
1
/
2
)
a
i
i (t ) t
Most commonly used algorithm i
T (t )
- Properties
•
•
•
•
23
Neighbour lists
• Complexity of force calculations ~O(N2)
But there are only
often n << N
interactions,
24
Neighbour lists
• Verlet lists
- Idea: introduce a list, where particles are included which are
located within interaction sphere
- Also introduce a “reservoir”, where particles outside Rc are
stored, so that unknown particles cannot become neighbors
in next steps
- Do an update of the list
every n steps
• Either statically with fixed n
• Or dynamically with an update
criterion
- There exists an optimal Rskin
25
Neighbour lists
• Linked lists method ~ O(N)
Interacts with atoms
in 26 neighbour cells
26
Measuring
•
1
A 
NT
NT
 A t 
t 1
• Kinetic energy + Potential energy = Total energy
2
1
K  t    mi vi  t   V  t    u ri  t   rj  t 
2 i
i j i


• Temperature kBT / 2 per degree of freedom
• The caloric curve E(T)
• Mean square displacement

 2
!Periodic boundary conditions
MSD  r (t )  r (0)
 2
1 
• Diffusion coefficient D  lim t 
r (t )  r (0)
6t
• Correlation functions
27
MD (and MC) as optimization tool
• Simulated annealing
Cooling schemes
- Start at high T, decrease T in small steps (cooling schedule)
- Easy to understand & implement
- Drawback: might be easily trapped in local minima
28
Parallel strategies
• Atom decomposition
-
Atoms are distributed among processors
All coordinates are exchanged before computing forces
OK for long range interactions
Easy to implement
• Force decomposition
- Each processor calculates the interactions for certain atom pairs
• Spatial decomposition
- Subdivides space and assigns each processor a particular subregion
- Atoms are allowed to move from one processor to its neighbours
- More complex to implement (similar to linked lists)
29
Analyzing the sequential code
• md.c and lmd.c
"Scientific computing and visualization", A. Nakano
University of Southern California
• Code description, details at
http://cacs.usc.edu/education/cs596.html
- “Basic molecular dynamics algorithms”
- “Linked-list cell MD algorithm”
30
Guidelines for final report
• Test energy conservation as function of t
• Compare speed of md.c and lmd.c
• Add possibility to save configuration, and to start from
a previous run. Allow temperature rescaling and
equilibration.
• Study the behavior of the caloric curve E(T) by means
of constant energy runs at a fixed density, starting
from a crystalline arrangement (r=0.6-0.8, Tmax=1.5).
• Insert calculations of the total linear and angular
momenta. Check their conservation.
• Insert calculation of the mean square displacement.
• Remove periodic boundary conditions and study a free
cluster.
31