High Performance Parallel Programming

1
High Performance Parallel Programming

• Dirk van der Knijff
• Advanced Research Computing
• Information Division

2
High Performance Parallel Programming

• Lecture 12 Particle-Particle Methods

3
Introduction

• Particle-Particle Methods are used when the
  number of particles is low enough to consider
  each particle and its interactions with the
  other particles
• Generally dynamic systems - i.e. we either watch
  them evolve or reach equilibrium
• Forces can be long or short range
• Numerical accuracy can be very important to
  prevent error accumulation
• Also non-numerical problems

4
Non-numeric problems

• e.g. Economic simulation of stock exchange
• Follows similar rules to numeric examples but
  uses different rules of interaction
• Often mix of Forces to be calculated
• Forces may be very complex
• Rules should be examined carefully to determine
  best ways of implementing
• Many systems oscillatory and/or chaotic

5
Molecular dynamics

• Many physical systems can be represented as
  collections of particles
• Physics used depends on system being
  studiedthere are different rules for different
  length scales
• 10-15-10-9m - Quantum Mechanics
• Particle Physics, Chemistry
• 10-10-1012m - Newtonian Mechanics
• Biochemistry, Materials Science, Engineering,
  Astrophysics
• 109-1015m - Relativistic Mechanics
• Astrophysics

6
Examples

• Galactic modelling
• need large numbers of stars to model galaxies
• gravity is a long range force - all stars
  involved
• very long time scales (varying for universe..)
• Crack formation
• complex short range forces
• sudden state change so need short timesteps
• need to simulate large samples
• Weather
• some models use particle methods within cells

7
General Picture

• Models are specified as a finite set of particles
  interacting via fields.
• The field is found by the pairwise addition of
  the potential energies, e.g. in an electric
  field
• The Force on the particle is given by the field
  equations

8
Simulation

• Define the starting conditions, positions and
  velocities
• Choose a technique for integrating the equations
  of motion
• Choose a functional form for the forces and
  potential energies. Sum forces over all
  interacting pairs, using neighbourlist or similar
  techniques if possible
• Allow for boundary conditions and incorporate
  long range forces if applicable
• Allow the system to reach equilibrium then
  measure properties of the system as it involves
  over time

9
Starting conditions

• Choice of initial conditions depends on knowledge
  of system
• Each case is different
• may require fudge factors to account for unknowns
• a good starting guess can save equibliration time
• but many physical systems are chaotic..
• Some overlap between choice of equations and
  choice of starting conditions

10
Integrating the equations of motion

• This is an O(N) operation. For Newtonian dynamics
  there are a number of systems
• Euler (direct) method - very unstable, poor
  conservation of energy

11
Direct method example

• Consider a rock thrown straight up

12
Direct method example
13
Verlet Method

• Higher order method in which velocities cancel
• From Taylor expansion about rt
• Note for mathematicians dt is normalized to 1

14
Verlet example

• Note Perfect agreement is due to constant
  acceleration.
• This is not the case in most situations...

15
Leap-frog algorithm

• Verlet algorithm can give rise to numerical
  inaccuracy
• Velocities must be calculated seperately if
  required
• Leapfrog uses velocities dt/2 out of phase
• Velocity if required is easily found

16
Comparison of Verlet schemes
t-1
t
t1
t-1
t
t1
t-1
t
t1
Original Verlet Algorithm
t-1
t
t1
t-1
t
t1
t-1
t
t1
t-1
t
t1
Leapfrog Algorithm
17
Taylors Theorem

• Mean Value Theorem
• Taylors theorem

18
Gear Predictor-Corrector Algorithm 1

• The predictor rt1, vt1, at1,.. from values at
  time t
• Does not contain information from the Newtonian
  eqn. of motion Evaluate the forces using
• Calculate the corrected acceleration

19
Gear Predictor-Corrector 2

• The corrector - correct the predicted values
• Coefficients depend on the order of the
  differentialequations and of the Taylor series

20
Comparison of integration schemes

• Original Verlet algorithm
• 9N words of storage
• Good energy conservation
• Numerical imprecision - add small term to diff of
  large terms
• Leapfrog algorithm
• 6N words of storage (if accelerations not saved)
• Gear-predictor algorithm
• 15N words of storage (for 4-value algorithm)
• can be easily modified to deal with modified
  equations of motion (e.g. damping effects)

21
Timesteps

• Note Ive used a normalized value of 1 in the
  examples
• The choice of real dt depends on the curvature of
  the potential energy surface V. Important
  properties are
• Energy conservation, Etotal Ekinetic
  Epotential
• Comvergence - agreement with trajectories of the
  same system with smaller timesteps
• The minimum timestep is around 1/10 of the
  oscillation frequency of the potential
• tmin 0.1x2p/we - around 1 femtosecond (10-15s)
  for molecular systems

22
Force Evaluation

• This is the most crucial part of the simulation,
  and for most cases will be the most time
  consuming. It is an O(N) operation if no
  interactions can be excluded.
• A range of forms can be used for the potential
  with different degrees of accuracy and different
  evaluation times
• Inverse square laws
• Lennard-Jones potentials
• Harmonic opscillators
• Morse potentials
• Quantum mechanical laws

23
Inverse Square Laws

• Common in both astrophysical systems
  (gravitation) and molecular systems
  (electrostatic)
• Has several problems
• requires a square root calculation for each force
  calculation so is time consuming
• converges very slowly at large r, so distance
  dependent cutoffs are unreliable

24
Lennard-Jones potentials

• The most widely used intermolecular potential. It
  is quick to evaluate, can be altered to model a
  range of properties, and has some physical
  justification..
• The force on particle i due to j is thus

25
Computational Aspects

• Tricks of the trade
• Avoid square roots
• Forces 1/rn where n is odd Fxrx/rn1 (only r2
  terms required)
• Lennard-Jones cutoffs
• Simple distance cutoff will lead to integration
  problems
• Add a constant to make sure potential is zero at
  cutoff, or..
• ..add a constant linear function
• Interpolate unavoidable functions
• SQRT(r) in gravitational simulations (polynomial
  approximations)
• Exp(r) in morse potentials (table look-up)

26
Nearest Neighbour techniques

• Most of the computational effort in a simulation
  is spent in force evaluation. Full force
  evaluation requires 1/2N(N-1) seperate
  calculations.
• any technique which reduces this overhead will be
  useful, even if it complicates the code
  considerably
• different techniques are applicable to different
  systems, depending on the mobility of the
  particles
• Methods include
• Neighbourlists
• Multicell algorithms
• Multiple timestep methods

27
Verlet Neighbourlists

• If the force has a cutoff distance, a set of
  neighbours can be isolated for each particle
• All particles within the cutoff distance must be
  included
• Particles which might move within the cutoff
  distance after n steps must be included
• rneighbourlistn.vmax.dt rcutoff
• n is typically in the range 10-20
• Problems - large memory requirements, and for
  parallel implementations parallelising force
  calculations requires distribution of large
  irregular data structure every n steps

28
Neighbourlists
P1 P3, P4, P5 P2 P5 ...
P5
P3
P4
P1
P2
29
Multicell algorithms

• For short range forces, a spatial decomposition
  is possible
• Each particle is assigned to a box, as small as
  possible but larger than rcutoff
• Force evaluation then simply involves pairs of
  atoms in the same or adjacent boxes
• The box contents are tracked using a linked list
  - the contents of each box are treated as if
  threaded on a single chain, and an array of
  length N contains the connectivity information
• Used for solids without electrostatic potential
  terms

30
Multiple Timestep methods

• Related to Verlet neighbourlists, this recognises
  that long range interactions are less sensitive
  to motions than short range ones.
• Successive shells of coordinates are stored, the
  inner one updated each timestep, the outer ones
  updated progressively less often
• Useful for fluids with long range forces such as
  solutions of water with solvated ions
• Requires storage of extra sets of coordinates,
  and is very sensitive to timestep.

31
Multiple timestep
10dt
dt
32
Boundary conditions

• Most molecular dynamics simulations are used to
  sample the properties of some infinite system...
• Exceptions may include engineering simulations of
  sets of objects or astrophysical simulations, but
  in most cases some sort of boundary condition
  must be applied
• The simplest case is to have the particles in a
  box with completely elastic walls
• velocities are reversed when the trajectory leads
  out of the box, and the particle is reflected
  back into the simulation area
• this may be hard to use with a high order
  integration scheme a simple steep potential can
  be provided at the edges of the system
• choice depends on system being modelled

33
Periodic boundaries

• The most common solution is to have particles
  which leave one face of the simulation box
  re-enter through the opposite side.
• The simulation then becomes a simulation of an
  infinite crystal in which each cell exhibits the
  same dynamics
• Overcomes/avoids surface effects
• For short range forces, any Fij pair can have j
  taken from the simulation box or its nearest
  neighbours
• Long range forces may contain artefacts due to
  the periodicity

34
Periodic boundaries
1
1
t1
1
t
F2,3
F3,2
2
2
3
3
35
Long range forces

• Forces such as gravity and electrostatic
  potential do not converge quickly to bulk average
  values
• Simulations based on these forces may not scale
  well without 0(N2) calculations
• If the system damps out the forces a distance
  dependent cutoff may be used
• Faster treatments of large systems require
  particle mesh methods (later)

36
Observations

• The aim of particle-particle simulation is
  usually to find the statistical average
  properties of a system for which the full phase
  space cannot be integrated analytically
• The hope is that a long enough simulation will
  sample enough of the phase space (position,
  momentum) that the properties will be close to
  the thermodynamic averages
• Systems such as proteins, have many constraints
  which make the system non-ergodic. Phase space
  will contain closed regions, barriers and
  bottlenecks
• Care must be taken in these cases that starting
  conditions are as realistic as possible

37
High Performance Parallel Programming

• Tomorrow - Particle-Particle Methods continued