Autumn School, Sofia
Bulgaria 2012

DL_POLY/FIELD/MESO TUTORIAL
Software Solutions in Atomistic Modelling
ILIAN TODOROVCHIN YONG, MICHAEL SEATON
CSED, STFC DARESBURY LABORATORY, DARESBURYWARRINGTON WA4 4AD, CHESHIRE, ENGLAND, UK

Daresbury Laboratory
Autumn School, Sofia
Bulgaria 2012

Alice’s Wonderland (1865)
Lewis Carroll (Charles Lutwidge Dodgson)

Multiple Scales of Materials Modelling
Autumn School, Sofia
Bulgaria 2012

Autumn School, Sofia
Bulgaria 2012

Part 1
DL_POLY Project BackgroundDL_POLY Trivia
Autumn School, Sofia
Bulgaria 2012

• General purpose parallel (classical) MD simulation software• It was conceived to meet needs of CCP5 - The Computer
Simulation of Condensed Phases (academic collaboration
community)• Written in modularised Fortran90 (NagWare & FORCHECK
compliant) with MPI2 (MPI1+MPI-I/O) & fully self-contained• 1994 – 2012: DL_POLY_2 (RD) by W. Smith & T.R. Forester(funded for 6 years by EPSRC at DL). Now moved to a BSD
open source licence as DL_POLY_Classic.
• 2003 – 2012: DL_POLY_3 (DD) by I.T. Todorov & W. Smith(funded for 4 years by NERC at Cambridge). Up-licensed to
DL_POLY_4 in 2010 - free of charge to academic
researchers and at cost to industry (provided as source).
• ~ 14,000 licences taken out since 1994 (~1,600 annually)• ~ 1,200 registered FORUM members since 2005

Current Versions
Autumn School, Sofia
Bulgaria 2012

Written in modularised free formatted F90 (+MPI) with rigorous
code syntax (FORCHECK and NAGWare verified) and no external
library dependencies• DL_POLY_4 (version 3.2)
– Domain Decomposition parallelisation, based on domain
decomposition (no dynamic load balancing), limits: up to≈2.1×109 atoms with inherent parallelisation– Parallel I/O (amber netCDF) and radiation damage features– Free format (flexible) reading with some fail-safe featuresand basic reporting (but not fully fool-proofed)• DL_POLY Classic (version 1.9)
– Replicated Data parallelisation, limits up to ≈30,000 atoms
with good parallelisation up to 100 (system dependent)
processors (running on any processor count)– Hyper-dynamics, Temperature Accelerated Dynamics,
Solvation Dynamics, Path Integral MD– Free format reading (somewhat rigid)

DL_POLY on the Web
Autumn School, Sofia
Bulgaria 2012

WWW:
http://www.ccp5.ac.uk/DL_POLY/
FTP:
ftp://ftp.dl.ac.uk/ccp5/DL_POLY/
DEV:

http://ccpforge.cse.rl.ac.uk/gf/project/dl-poly/
http://ccpforge.cse.rl.ac.uk/gf/project/dl_poly_classic/
FORUM:
http://www.cse.stfc.ac.uk/disco/forums/ubbthreads.php/

Further Information
Autumn School, Sofia
Bulgaria 2012

W. Smith and T.R. Forester,J. Molec. Graphics (1996), 14, 136
W. Smith, C.W. Yong, P.M. Rodger,
Molecular Simulation (2002), 28, 385
I.T. Todorov, W. Smith, K. Trachenko, M.T. Dove,
J. Mater. Chem. (2006), 16, 1611-1618
W. Smith (Guest Editor),
Molecular Simulation (2006), 32, 933
I.J. Bush, I.T. Todorov and W. Smith,Comp. Phys. Commun. (2006), 175, 323-329

DL_POLY_DD Development Statistics
Autumn School, Sofia
Bulgaria 2012

Published lines of code [x 1,000]

140
120
100
80
60
40
20

0
2002

2004

2006

2008

Year

2010

2012

2014

DL_POLY Licence Statistics
Autumn School, Sofia
Bulgaria 2012

Licences

Downloads

1600

1200

800

400

0
1992

1996

2000

2004

Year

2008

2012

DL_POLY Licence Statistics
Autumn School, Sofia
Bulgaria 2012

DL_POLY Licences
2011 by Sub-Areas

DL_POLY Licences
2010 by Sub-Areas
Latin America
4%Europe-REST
5%

Africa2%

UK
11%

Australia and
New Zealandand Philipnes
1%

Europe-REST
4%

Latin America
8%

Africa
3%

Asia
38%

Asia
39%

UK
10%

North America
22%

North America
15%

EU-REST
17%

Australia and
New Zealandand Philipnes
2%

EU-REST
19%

DL_POLY Licence Statistics
Autumn School, Sofia
Bulgaria 2012

DL_POLY Licences
2010 by Science Domain
Software5%

DL_POLY Licences
2011 by Science Domain

Bio-Chemsitry5%

Bio-Chemsitry
4%

Physics
28%
Physics
33%

Chemistry
33%

Chemistry
38%

Other
3%

Other
1%

Mechanics
10%

Engineering
13%

Mechanics
3%

Materials
7%
Engineering17%Autumn School, Sofia
Bulgaria 2012

Part 2The Molecular Dynamics Method

Molecular Dynamics: Definitions
Autumn School, Sofia
Bulgaria 2012

• Theoretical tool for modelling the detailed microscopic behaviourof many different types of systems, including; gases, liquids, solids,
polymers, surfaces and clusters.
• In an MD simulation, the classical equations of motion governing
the microscopic time evolution of a many body system are solved
numerically, subject to the boundary conditions appropriate for thegeometry or symmetry of the system.
• Can be used to monitor the microscopic mechanisms of energy
and mass transfer in chemical processes, and dynamical properties
such as absorption spectra, rate constants and transport properties
can be calculated.
• Can be employed as a means of sampling from a statisticalmechanical ensemble and determining equilibrium properties.
These properties include average thermodynamic quantities(pressure, volume, temperature, etc.), structure, and free energies
along reaction paths.

Molecular Dynamics for Beginners
Autumn School, Sofia
Bulgaria 2012

MD simulations are used for:• Microscopic insight: we can follow the
motion of a single molecule (glass ofwater)• Investigation of phase change (NaCl)• Understanding of complex systems likepolymers (plastics – hydrophilic andhydrophobic behaviour)Example: Simulation of Argon
Autumn School, Sofia
Bulgaria 2012

Pair Potential:r cut  12   6 
V (r )  4       r   r  
Lagrangian:
N 1
1N
 
L(ri , vi )   mi vi 2    V (rij )2 i
i j iLennard -Jones Potential
Autumn School, Sofia
Bulgaria 2012  12   6 
V (r )  4       r   r  

V(r)




r
r cu tPair-wise radial distanceEquations of Motion
Autumn School, Sofia
Bulgaria 2012d L   L   
dt   vi   ri

mi ai  FiN 
Fi   f ij
j i



f ij   iV (rij )Lagrange Equation –time evolution

Force Evaluation –particle interactions

Boundary Conditions
Autumn School, Sofia
Bulgaria 2012

2D cubic periodic• None – biopolymer
simulations• Stochastic boundaries– biopolymers• Hard wall boundaries– pores, capillaries• Periodic boundaries –most MD simulations

Periodic Boundary Conditions
Autumn School, Sofia
Bulgaria 2012

Triclinic

Truncated octahedron

Hexagonal prism

Rhombic dodecahedron

System Properties: Static (1)Autumn School, Sofia
Bulgaria 2012

• Kinetic Energy:

1N 2
K . E.   mi vi
2 i• Temperature:2T
K . E.
3 Nk BN

Uc    V (rij )• Configuration Energy:• Pressure:i j i

1 N  PV  Nk BT   ri  f i
3 i

3 2 2  3Nk B  (U c ) Nk BT 1 
NVE
2
 2Cv 2

• Specific heat:System Properties: Static (2)Autumn School, Sofia
Bulgaria 2012

Structural Properties– Pair correlation (Radial Distribution Function):n(r )Vg(r ) 
 2
2
4  r r NN

   (r  rij )i j i– Structure factor:sin(kr)S (k )  1  4  
g (r )  1 r 2 dr
0kr
– Note: S(k) available from x-ray diffraction

Radial Distribution Function (RDF)
Autumn School, Sofia
Bulgaria 2012

R

R

Typical RDF
Autumn School, Sofia
Bulgaria 2012

g(r)

1.0

separation (r)System Properties: Dynamic
Autumn School, Sofia
Bulgaria 2012

Single Correlation Functions:• Mean Square Displacement (Einstein
Relation)• Velocity Autocorrelation (Green-Kubo
relation)

System Properties: Collective
Autumn School, Sofia
Bulgaria 2012

Collective Correlation Functions:• van Hove correlation functions (general,
self, distinct)• Particle density• Intermediate scattering function• Dynamic structure factor

Correlation Function Uses
Autumn School, Sofia
Bulgaria 2012

• Complete description of bulk dynamical
properties• Space-time Fourier Transform of van Hove
function• Elastic properties of materials• Energy dissipation• Sound propagation
Obtained directly from neutron scatteringAutumn School, Sofia
Bulgaria 2012

Part 2
DL_POLY Basics & Algorithms

DL_POLY Supported Molecular Entities
Autumn School, Sofia
Bulgaria 2012

Rigid
molecules

Point ionsand atoms

Flexibly
linked rigid
molecules

Polarisable
ions (core+shell)
Flexiblemolecules

Rigid
bonds

Rigid bondlinked rigid
molecules

Force Field Definitions – I
Autumn School, Sofia
Bulgaria 2012

• particle: rigid ion or atom (charged or not), a core or a shell of apolarisable ion(with or without associated degrees of freedom), a
massless charged site. A particle is a countable object and has a
global ID index.

• site: a particle prototype that serves to defines the chemical &physical nature (topology/connectivity/stoichiometry) of a particle(mass, charge, frozen-ness). Sites are not atoms they areprototypes!
• Intra-molecular interactions: chemical bonds, bond angles,dihedral angles, improper dihedral angles, inversions. Usually, themembers in a unit do not interact via an inter-molecular term.
However, this can be overridden for some interactions. These aredefined by site.

• Inter-molecular interactions: van der Waals, metal (EAM, Gupta,
Finnis-Sinclair, Sutton-Chen), Tersoff, three-body, four-body.
Defined by species.

Force Field Definitions – II
Autumn School, Sofia
Bulgaria 2012

• Electrostatics: Standard Ewald*, Hautman-Klein (2D) Ewald*,
SPM Ewald (3D FFTs), Force-Shifted Coulomb, Reaction Field,
Fennell damped FSC+RF, Distance dependent dielectric constant,
Fuchs correction for non charge neutral MD cells.
• Ion polarisation via Dynamic (Adiabatic) or Relaxed shell model.
• External fields: Electric, Magnetic, Gravitational, Oscillating &Continuous Shear, Containing Sphere, Repulsive Wall.
• Intra-molecular like interactions: tethers, core shells units,
constraint and PMF units, rigid body units. These are also definedby site.
• Potentials: parameterised analytical forms defining theinteractions. These are always spherically symmetric!

• THE CHEMICAL NATURE OF PARTICLES DOES NOT CHANGE
IN SPACE AND TIME!!! *Force Field by Sums
Autumn School, Sofia
Bulgaria 2012

N'
 

 
V( r1 , r2 ,....., rN )   U pair (| ri  rj |) 
i, j
N'

U

Tersoff

i, j,kN'   ri , rj , rk    U3body

1 N' q i q j  4π 0 i, j | ri  rj |N'   ri , rj , rk    U 4body

i, j,k    ri , rj , rk , rn  i, j,k, n

N
 N'
 N'





ε metal   Vpair (| ri  rj |)   F  ρ ij (| ri  rj |)   
 i, ji
 i, j
N bond U bondi bonddihedi dihedi tether  ,angle ra , rb , rc  

N invers   i dihed , ra , rb , rc , rd    U invers   iinvers , ra , rb , rc , rd  i invers

N tetherUii angleN dihed

U

N

angle
 
i bond , ra , rb    U angle

tether

N core-shell 
i tether , rt , rt 0    U core -shelli core-shellN 
icore -shell , | ri  rj |   Φexternali 1ri 

DL_POLY Boundary Conditions
Autumn School, Sofia
Bulgaria 2012

• None (e.g. isolated macromolecules)• Cubic periodic boundaries• Orthorhombic periodic boundaries• Parallelepiped (triclinic) periodic boundaries• Truncated octahedral periodic boundaries*• Rhombic dodecahedral periodic boundaries*• Slabs (i.e. x,y periodic, z non-periodic)

Assumed Parallel Architecture
Autumn School, Sofia
Bulgaria 2012

DL_POLY is designed for homogenious
distributed parallel machines
M0

P0

P4

M4

M1

P1

P5

M5

M2

P2

P6

M6

M3

P3

P7

M7

Replicated Data Strategy – I
Autumn School, Sofia
Bulgaria 2012

A

B

C

D

Initialize

Initialize

Initialize

Initialize

ForcesForcesForcesForcesMotionMotionMotionMotionStatistics

Statistics

Statistics

Statistics

Summary

Summary

Summary

SummaryReplicated Data Strategy – II
Autumn School, Sofia
Bulgaria 2012

• Every processor sees
the full system• No memory distribution(performance overheads
and limitations increase
with increasing systemsize)• Functional/algorithmic
decomposition of theworkload• Cutoff ≤ 0.5 min systemwidth• Extensive global
communications(extensive overheads
increase with increasingsystem size)

Parallel (RD) Verlet List
Autumn School, Sofia
Bulgaria 2012

P0

P0

P0

1,2

1,3

1,4

1,5

1,6

1,7

2,3

2,4

2,5

2,6

2,7

2,8

3,4

3,5

3,6

3,7

3,8

3,9

4,5

4,6

4,7

4,8

4,9

4,10

5,6

5,7

5,8

5,9

5,10

5,11

6,7

6,8

6,9

6,10

6,11

6,12

7,8

7,9

7,10

7,11

7,12

8,9

8,10

8,11

8,12

8,1

9,10

9,11

9,12

9,1

9,2

10,11

10,12

10,1

10,2

10,3

11,12

11,1

11,2

11,3

11,4

12,1

12,2

12,3

12,4

12,5

Distributed
list!

Domain Decomposition MD
Autumn School, Sofia
Bulgaria 2012

A

B

C

D

Linked Cell Lists
Autumn School, Sofia
Bulgaria 2012

•

•

•

Linked lists provide an elegant way to scale short-ranged
two body interactions from O(N2/2) to ≈O(N). The
efficiency increases with increasing link cell partitioning– as a rule of thumb best efficacy is achieved for cubiclike partitioning with number of link-cells per domain ≥ 4
for any dimension.
Linked lists can be used with the same efficiency for 3body (bond-angles) and 4-body (dihedral & improper
dihedral & inversion angles) interactions. For these thelinked cell halo is double-layered and since cutoff3/4-body≤ 0.5*cutoff2-body which makes the partitioning more
effective than that for the 2-body interactions.
The larger the particle density and/or the smaller thecutoff with respect to the domain width, (the larger the
sub-selling and the better the spherical approximation of
the search area), the shorter the Verlet neighbour-list
search.

Linked Cell List Idea
Autumn School, Sofia
Bulgaria 2012

6

12

Cell number
1

2

Head of Chain
List

6

Link
List 1

3

4

3

4

5

6

17

5

10

16

12

10
2

Link
Cell 2

7

8

9

16

17

0

10 11 12 13 14 15 16 17 18 19 20

Atom numberSub-celling of LCs
Autumn School, Sofia
Bulgaria 2012

1 2 3 4 5 6 7• Provides dynamicallyadjustable workload forvariable local density
and VNL speed up of ≈30% (45% theoretically).
• Provides excellent serial
performance, extremely
close to that of BrodeAhlrichs method forconstruction of the Verlet
neighbour-list whensystem sizes are smaller< 5000 particles.

Parallel Force Calculation
Autumn School, Sofia
Bulgaria 2012•

Bonded forces:- Algorithmic decomposition- Interactions managed by bookkeeping arrays,
i.e. explicit bond definition- Shared bookkeeping arrays•

Nonbonded forces:– Distributed Verlet neighbour list (pair forces)– Link cells (3,4-body forces)•

Implementations differ between DL_POLY_4 & C!

Scheme for Distributing Bond Forces
Autumn School, Sofia
Bulgaria 2012

A2
A1

A3

P0

A6

A4

P1

A5

P2

A10

A8
A7

P3

A9

A11

P4

P5

A16

A14

A12
A13

P6

A15

P7

A17

P8

DL_POLY_C and Bonded Forces

P0Localforce
terms
P1Localforce
terms
P2Local
force
terms

Processors

Molecular
force fielddefinition

Global Force Field

Autumn School, Sofia
Bulgaria 2012

DL_POLY_4 and Bonded Forces

P0Localatomicindices
Tricky!

P1Localatomicindices
P2Localatomicindices

Processor Domains

Molecular
force fielddefinition

Global force field

Autumn School, Sofia
Bulgaria 2012

Ensembles and Algorithms
Autumn School, Sofia
Bulgaria 2012

Integration:
Available as velocity Verlet (VV) or leapfrog Verlet (LFV)
generating flavours of the following ensembles• NVE• NVT (Ekin) Evans• NVT Andersen^, Langevin^, Berendsen, Nosé-Hoover, GST• NPT Langevin^, Berendsen, Nosé-Hoover, MartynaTuckerman-Klein^• NT/NPnAT/NPnT Langevin^, Berendsen, Nosé-Hoover,Martyna-Tuckerman-Klein^
Constraints & Rigid Body Solvers:• VV dependent – RATTLE, No_Squish, QSHAKE*• LFV dependent – SHAKE, Euler-Quaternion, QSHAKE*

Integration Algorithms
Autumn School, Sofia
Bulgaria 2012

Essential Requirements:• Computational speed• Low memory demand• Accuracy• Stability (energy conservation, no drift)• Useful property - time reversibility• Extremely useful property – symplecticness= time reversibility + long term stability

Integration: Essential Idea
Autumn School, Sofia
Bulgaria 2012r’ (t+t)v (t)tr (t+t)r (t)f(t)t2/m[r (t), v(t), f(t)][r (t+t), v(t+t), f(t+t)]

Simulation Cycle and Integration Schemes
Autumn School, Sofia
Bulgaria 2012

Setup

Taylor expansion:

t 2 f n
rn1  rn  t vn 
 O t 3
2 m

Leapfrog Verlet (LFV)

Velocity Verlet (VV)0. x i (t ), vi (t  12 t )VV1. 0. x i (t ), vi (t ), f i (t )1. f i (t )  calculated afreshVV1. 1. vi (t  12 t )  vi (t ) 

Set up initialsystemForcesCalculateforcesMotionCalculatemotionStats.

Accumulate
statistical data

Results

Summarise
simulation 

t f i (t )2 mif i (t )mi VV1. 2. x (t  t )  xi (t )  t v (t  1 t )
i
i
2
2
1
3. xi (t  t )  xi (t )  t vi (t  2 t )VV2. 0. f i (t  t )  calculated afresh
2. vi (t  12 t )  vi (t  12 t )  t

VV2.1 . vi (t  t )  vi (t  12 t ) t f i (t  t )2
mi

Integration Algorithms: Leapfrog Verlet
Autumn School, Sofia
Bulgaria 2012

f n-2
rn-2

f n-1
rn-1
vn-3/2 n 1/ 2  n 1/ 2 t  nvi
 vi
 Fi   (t 3 )mi n 1  n
 n 1/ 2
4
ri  ri  tvi  (t )fn
rn

Application in Practicern+1

vn-1/2 vn+1/2

Discrete time

 n 1/ 2  n 1/ 2 t  n
vi
 vi

Fi
mi
 n 1  n
 n 1/ 2
ri  ri  tvi
 n 1/ 2  n 1/ 2
 n vi vivi 2Integration Algorithms: Velocity Verlet
Autumn School, Sofia
Bulgaria 2012

f n-2

f n-1

 n 1  n
 n t 2  nri  ri  tvi Fi   (t 4 )2mi n 1  n t  n  n 1vi  vi ( Fi  Fi )   (t 2 )
2mi
fn

f n+1

rn-2

rn-1

rn

rn+1

vn-2

vn-1

vn

vn+1

Discrete time

Application in Practice n 1/ 2  nt  nvi
 vi 
Fi2mi
 n 1  n
 n 1/ 2
ri  ri  tvi
 n 1  n 1/ 2
t  n 1vi  viFi
2mi

Constraint Solvers
Autumn School, Sofia
Bulgaria 2012t 2 f n  g n
rn1  rn  t vn 
 O t 32mt f n  hn
vn1  vn 1 
 Ot 3 
2
2 m

 

Taylor expansions:
SHAKE

iu


Gij i

io

u
d
 ij
dij

j

o
o
dij
j
o


Gij   G ji  g ij d ij
 2 u 2
ij (d ij  d ij )o ug ij 
2
2td ij  d ij

1

ij



1
1

mi m jRATTLERATTLE_R (SHAKE)G ji

o2 u 2
ij (d ij  d ij )g ij  2  o  utd ij  d ij

ju


vio
i

j


dij

RATTLE_V


vio

o


H ij   H ji  hij  d ij  o ij (vi  v j )  d ij
hij 
t
d ij2

Extended Ensembles in VV casting
Autumn School, Sofia
Bulgaria 2012

Velocity Verlet integration algorithms can be naturally derived from thenon-commutable Liouvile evolution operator by using a second orderSuzuki-Trotter expansion. Thus they are symplectic/true ensembles(with conserved quantities) warranting conservation of the phase-spacevolume, time-reversibility and long term numerical stability…
Examplary VV Expansion of NVE to NVE kin, NVT, NPT & NσT
VV1 :

VV2 :

x i (t ), vi (t ), f i (t )x i (t  t ), vi (t  12 t ), f i (t  t )  afresh

Thermostat(t  t  14 t ): 14 t

Barostat(t  t  12 t ): 12 t

Thermostat(t  14 t  t  12 t ) : 14 t
t f i (t )vi (t  t )  vi (t ) 
2 mi
1
2

: t
1
2

t
xi (t  t )  xi (t ) vi (t  12 t ) : t
2
RATTLE _ R(t  t  t ): tvi (t  t )  vi (t  12 t ) t f i (t  t ) 1: 2 t2
mi

RATTLE _ V (t  12 t  t  t ): t

Thermostat(t  12 t  t  34 t ): 14 t

Barostat(t  12 t  t  t ): 12 t

Thermostat(t  32 t  t  t ): 14 tOther Integration Algorithms
Autumn School, Sofia
Bulgaria 2012

• Gear Predictor-Corrector – generally easily
extendable to any high order of accuracy.
It is used in satellite trajectory
calculations/corrections. However, lacking
long term stability.
• Trotter derived evolution algorithms –generally easily extendable to any highorder of accuracy. Symplectic.

Base Functionality
Autumn School, Sofia
Bulgaria 2012

• Molecular dynamics of polyatomic systems with
options to save the micro evolution trajectory at
regular intervals• Optimisation by conjugate gradients method or zero
Kelvin annealing• Statistics of common thermodynamic properties(temperature, pressure, energy, enthalpy, volume)with options to specify collection intervals and stacksize for production of rolling and final averages• Calculation of RDFs and Z-density profiles• Temperature scaling, velocity re-Gaussing• Force capping in equilibration

DL_POLY_4 Specials
Autumn School, Sofia
Bulgaria 2012

• Radiation damage driven features:• defects analysis• boundary thermostats• volumetric expansion• replay history• variable time step algorithm• Extra ensembles:• Langevin, Andersen, MTK, GST• extensions of NsT to NPnAT and NPnT• Infrequent k-space Ewald evaluation• Direct VdW• Direct Metal• Force shifted VdW• I/O driven features Parallel I/O & netCDF• Extra ReportingAutumn School, Sofia
Bulgaria 2012

Part 3
DL_POLY I/O FilesI/O Files
Autumn School, Sofia
Bulgaria 2012

• Crystallographic(Dynamic) data
• Simulation
Control data
• Molecular /Topological Data
• Tabulated twobody potentials• EAM potential
data• Reference data
for DEFECTS• Restart data

• Finalconfiguration• Simulation
summary data
• Trajectory Data
• Defects data
• MSD & T inst data
• Statistics data
• Best CGM
configuration• RDF data
• Z density data
• Restart data

DL_POLY Units
Autumn School, Sofia
Bulgaria 2012

Internally, DL_POLY uses atomic scale units:• Mass• Charge• Length• Time• Force• Energy- mass of H atom (D) [Daltons]- charge on proton (e)- Angstroms (Å)- picoseconds (ps)- D Å ps-2
- D Å2 ps-2 [10 J mol-1]

pressure is expressed in k-atm for I/O

Acceptable DL_POLY Units
Autumn School, Sofia
Bulgaria 2012

UNITS directive in FIELD file allows to opt for thefollowing energy units• Internal DL_POLY units• Electron-volts• Kilo calories per mol• Kilo Joules per mol• Kelvin per Boltzmann- 10 J mol-1- eV- k-cal mol-1- k-J mol-1- K Boltzmann-1

All interaction MUST have the same energy units!

CONTROL File
Autumn School, Sofia
Bulgaria 2012

•
•
•
•

SIMULATION CONTROLFree Format
Mandatory
Driven by keywords:

keyword [options] {data}e.g.:ensemble NPT Hoover 1.0 8.0

CONFIG [REVCON,CFGMIN] File
Autumn School, Sofia
Bulgaria 2012

• Initial atomic coordinates• Format- Integers(I10)- Reals(F20)- Names(A8)• Mandatory• Units:- Position - Angstroms (Å)- Velocity – Å ps-1
- Force - D Å ps-2
• Construction:- Some kind of GUI essential
for complex systems

FIELD File
Autumn School, Sofia
Bulgaria 2012

Force Field specification• Mandatory• Format:•-

Integers
Reals
Names
Keywords

(I5)
(F12)
(A8)
(A4)

Maps on to CONFIG filestructure• Construction•- Small systems - by hand
- Large systems – nfold orGUI!

TABLE File
Autumn School, Sofia
Bulgaria 2012

Defines non-analytic
pair (vdw) potentials• Format- Integers(I10)- Reals(F15)- Names(A8)• Conditional, activatedby FIELD file option• Potential & Force• NB force (here) is:•
G (r )  r U (r )rTABEAM File
Autumn School, Sofia
Bulgaria 2012

Defines embedded atom potentials• Format- Integers(I10)- Reals(F15)- Names(A8)• Conditional, activated by FIELD file option• Potentials only
• pair, embed & dens keywords for atom types
followed by data records (4 real numbers perrecord)• Individual interpolation arrays•

REVOLD [REVIVE] File
Autumn School, Sofia
Bulgaria 2012•

Provides program restart capability•

File is unformatted (not human readable)•

Contains thermodynamic accumulators, RDF
data, MSD data and other checkpoint data•

REVIVE (output file) ---> REVOLD (input file)OUTPUT File
Autumn School, Sofia
Bulgaria 2012

•
•
•

•

Provides Job Summary (mandatory!)
Formatted to be human readable
Contents:- Summary of input data
- Instantaneous thermodynamic data at selected intervals- Rolling averages of thermodynamic data
- Statistical averages- Final configuration- Radial distribution data
- Estimated mean-square displacements and 3D diffusion
coefficient
Plus:- Timing data, CFG and relaxed shell model iteration data
- Warning & Error reports

STATIS File
Autumn School, Sofia
Bulgaria 2012

• System properties atintervals selected by user
• Optional• Formatted (I10,E14)• Intended use: statisticalanalysis (e.g. error) andplotting vs. time.
• Recommend use with GUI!
• Header:- Title- Units• Data:- Time step, time, #entries- System data

HISTORY File
Autumn School, Sofia
Bulgaria 2012

Configuration data at user
selected intervals• Formatted• Optional
Header:• Title• Data level, cell key, number
Configuration data:• Time step and data keys
• Cell Matrix• Atom name, mass, charge• X,Y,Z coordinates (level 0)• X,Y,Z velocities (level 1)• X,Y,Z forces (level 2)

RDFDAT [ZDNDAT] File
Autumn School, Sofia
Bulgaria 2012

Formatted (A8,I10,E14)• Plotable• Optional• RDFs from pair forces• Header:- Title- No. plots & length of plot
• RDF data:- Atom symbols (2)- Radius (A) & RDF- Repeated……
• ZDNDAT file has same format•

DL_POLY_4 Extra Files
Autumn School, Sofia
Bulgaria 2012•

REFERENCE file– Reference structure to compare against•

DEFECTS file– Trajectory file of vacancies and interstitials
migration•

MSDTMP file– Trajectory like file containing the each particle’sSqrt(MSDmean) and Tmean

Autumn School, Sofia
Bulgaria 2012

Part 5
DL_POLY_4 PerformanceProof of Concept on IBM p575 2005
Autumn School, Sofia
Bulgaria 2012

300,763,000 NaCl with full SPME electrostatics
evaluation on 1024 CPU cores

Start-up timeTimestep time
FFT evaluation

 1 hour 68 seconds
 55 seconds

In theory ,the system can be seen by the eye. Although
you would need a very good microscope – the MD cellsize for this system is 2μm along the side and as the
wavelength of the visible light is 0.5μm so it should betheoretically possible.

Benchmarking BG/L Jülich 2007
Autumn School, Sofia
Bulgaria 2012

Perfect
MD step total
Link cells
van der Waals
Ewald real
Ewald k-space

16000
14000

Speed Gain

12000
10000
8000
6000
4000

14.6 million particle Gd2Zr2O7 system

2000
2000

4000

6000

8000

10000

Processor count

12000

14000

16000

DL_POLY_4 Weak Scaling
Autumn School, Sofia
Bulgaria 2012

Solid Ar (32'000 atoms per CPU)
NaCl (27'000 ions per CPU)
SPC Water (20'736 ions per CPU)

1000

800

Speed Gain

fec
r
pe

n
tio
a
s
eli
l
l
a
33 million atoms
arpt
28 million atoms

600

400

21 million atomsnallelisatio
r
a
p
d
o
o
g

200
max load
max load
max load

0
0

200

400

600

Processor Count

700'000 atoms per 1GB/CPU
220'000 ions per 1GB/CPU
210'000 ions per 1GB/CPU

800

1000

DL_POLY_4 RB v/s CB
on 450,000 particles
Autumn School, Sofia
Bulgaria 2012

Scaling
109

ICE7

steps per second8

ICE7_CB

7
6
5
4

3
2
1
0
0

100

200

300
Np

400

500

600

Benchmarking Main Platforms
Autumn School, Sofia
Bulgaria 2012

9
8

6

-1

Evaluations [s ]7

CRAY XT4 SC
CRAY XT4 DC
CRAY XT3 SC
CRAY XT3 DC
BG/L
BG/P
IBM p575
3GHz Woodcrest DC

5
4
3
2
1

3.8 million particle Gd2Zr2O7 system

0
0

500

1000

Processor count

1500

2000

I/O Solutions in DL_POLY_4
Autumn School, Sofia
Bulgaria 2012

1. Serial read and write (sorted/unsorted) – where only a single
MPI task, the master, handles it all and all the rest communicate inturn to or get broadcasted to while the master completes writing a
configuration of the time evolution.
2. Parallel write via direct access or MPI-I/O (sorted/unsorted) –where ALL / SOME MPI tasks print in the same file in some orderly
manner so (no overlapping occurs using Fortran direct accessprinting. However, it should be noted that the behaviour of thismethod is not defined by the Fortran standard, and in particular wehave experienced problems when disk cache is not coherent with thememory).
3. Parallel read via MPI-I/O or Fortran4. Serial NetCDF read and write using NetCDF libraries formachine-independent data formats of array-based, scientific data(widely used by various scientific communities).

MPI-I/O Write Performance for
216,000 Ions of NaCl on XT5Autumn School, Sofia
Bulgaria 2012

3.09
3.10
3.09
3.10
Cores
I/O Procs Time/s Time/s Mbyte/s Mbyte/s
32
32 143.30
1.27
0.44
49.78
64
64
48.99
0.49
1.29 128.46
128
128
39.59
0.53
1.59 118.11
256
128
68.08
0.43
0.93 147.71
512
256 113.97
1.33
0.55
47.60
1024
256 112.79
1.20
0.56
52.47
2048
512 135.97
0.95
0.46
66.39

MPI-I/O Read Performance for
216,000 Ions of NaCl on XT5Autumn School, Sofia
Bulgaria 20123.10 New
3.10 NewCores
I/O Procs Time/s Time/s Mbyte/s Mbyte/s
32
16
3.71
0.29
17.01 219.76
64
16
3.65
0.30
17.28 211.65
128
32
3.56
0.22
17.74 290.65
256
32
3.71
0.30
16.98 213.08
512
64
3.60
0.48
17.53 130.31
1024
64
3.64
0.71
17.32
88.96
2048
128
3.75
1.28
16.84
49.31

Autumn School, Sofia
Bulgaria 2012

Part 4
Obtaining & Building DL_POLY

DL_POLY Licensing and Support
Autumn School, Sofia
Bulgaria 2012

• Online Licence Facility at http://www.ccp5.ac.uk/DL_POLY/• The licence is• To protect copyright of Daresbury Laboratory• To reserve commercial rights• To provide documentary evidence justifying continued
support by UK Research Councils• It covers only the DL_POLY_4 package• Registered users are entered on the DL_POLY e-mailing list
• Support is available (under CCP5 and EPSRC service levelagreement to CSED) only to UK academic researchers• For the rest of the world there is the online user forum• Last but not least there is a detailed, interactive, selfreferencing PDF (LaTeX) user manualSupply of DL_POLY_4
Autumn School, Sofia
Bulgaria 2012

• Register at http://www.ccp5.ac.uk/DL_POLY/• Registration provides the decryption - procedureand password (sent by e-mail)• Source is supplied by anonymous FTP (in an
invisible directory)• Source is in a tarred, gzipped and encrypted form

• Successful unpacking produces a unix directory
structure• Test and benchmarking data are also available

DL_POLY_Classic Support
Autumn School, Sofia
Bulgaria 2012

• Full documentation of software supplied with source• Support is available through the online user forum or
the CCP5 user community
WWW:

http://www.ccp5.ac.uk/DL_POLY_CLASSIC/
FTP:
ftp://ftp.dl.ac.uk/ccp5/DL_POLY/
FORUM:
http://www.cse.stfc.ac.uk/disco/forums/ubbthreads.php/

Supply of DL_POLY_Classic
Autumn School, Sofia
Bulgaria 2012

• Downloads are available from CCPForge at
http://ccpforge.cse.rl.ac.uk/gf/project/dl_poly_classic/• No registration required – BSD licence• Download source from: CCPForge: Projects: DL_POLY
Classic: Files: dl_poly_classic: dl_poly_classic1.9• Sources is a in tarred and gzipped form

• Successful unpacking produces a unix directory
structure• Test data are also available

DL_POLY Directory Structure
Autumn School, Sofia
Bulgaria 2012

DL_POLY

build

Home of makefiles

sourceDL_POLY source code

executeHome of executable&
Working Directory

javaJava GUI source code

utility

Utility codes

data

Test data

Compiling the DL_POLY code
Autumn School, Sofia
Bulgaria 2012

• Read README.txt file in source and follow
instructions within• Copy Makefile from `build’ to`source’• Use `make machinetype’ to compile e.g

make gfortran• Executable in `execute’ directory• DLPOLY.X is DL_POLY_Classic executable• DLPOLY.Z is DL_POLY_4 executable• Standard executables have been madeavailable on many architectures

DL_POLY_C versus DL_POLY_4
for parallel jobs
Autumn School, Sofia
Bulgaria 2012

1.
2.
3.
4.
5.

Note differences in capabilities (e.g. linked rigid bodies) !!!
Less than 10,000 atoms (if in parallel)? – DL_POLY Classic
More than 30,000 atoms? – DL_POLY_4
Ratio cell_width/rcut < 3 (in any direction)? – DL_POLY_Classic
Less than 500 particles per processor? – DL_POLY_Classic

DL_POLY_Classic
Simple molecules (no SHAKE)• 8 or less, 10,000 atoms• 16 or less, 20,000 atoms• 32 or less, 30,000 atoms
Simple ionics• 16 or less, 10,000 atoms• 64 or less, 20,000 atoms• 128 or less, 30,000 atoms
Molecules (with SHAKE)• 64 max!

DL_POLY_4• Golden Rule 1: No fewer than
3x3x3 link cells per processor (ifin parallel)• Golden Rule 2: No fewer than
500 particles per processor (if inparallel)!

Autumn School, Sofia
Bulgaria 2012

Part 6
DL_POLY_Classic Functionality
W. SmithSpecial Algorithms
Autumn School, Sofia
Bulgaria 2012

• Hyperdynamics– Bias potential dynamics– Temperature accelerated dynamics– Nudged elastic band
• Solvation properties:– Energy decomposition– Spectroscopic solvent shifts– Free energy of solution• Metadynamics

I/O Files
Autumn School, Sofia
Bulgaria 2012

Standard Input

Special Input

Standard Output

Special Output

CONFIGREVOLDOUTPUTHISTORY

FIELD

TABLE

STATIS

RDFDAT

CONTROLTABEAM

REVIVEZDNDAT

REVCONHYPOLD

HYPRESEVENTSCFGBSNnn

Operation Type:CFGTRAnnStandard use

PROnn.XY

Hyperdyn./TAD

SOLVAT

Solvation

FREENG

Hyperdynamics

STEINHARDT

OptimisationZETA

METADYNAMICS

CFGMIN

Solvation Features
Autumn School, Sofia
Bulgaria 2012

• Molecular Solvation Energy• Energy decomposition• Energy distribution functions• Free Energy of Solvation• Mixed Hamiltonian method• Thermodynamic Integration• Solution Spectroscopy• Solvent induced shifts• Solvation relaxation

Solvation Files
Autumn School, Sofia
Bulgaria 2012

• SOLVAT– Breakdown of system energy based onmolecular types– Energies of ground and excited states• FREENG– Energy data for thermodynamic integration

Hyperdynamics
Autumn School, Sofia
Bulgaria 2012

Bias Potential Dynamics

Temperature Accelerated Dynamics



Metadynamics

Bias Potential Dynamics
Autumn School, Sofia
Bulgaria 2012

Original
Potential
Modified
Potential• Construct bias potential
to reduce well depth ofstate A.
• Bias potential is zero at
saddle point.
• Ratios of rates from state
A to states B, C, etc.
preserved:• Suitable bias potential:
Bias
Potential
State A

Bias Potential Dynamics 2
Autumn School, Sofia
Bulgaria 2012f  exp  H  d exp H  df  exp  H    V R   V R  d exp  H    V R   V R  d
N

f

A

N

N

N

f

f

N

N

N

N

N

b

A

N


A

N

b

N

N

b

N

b

 
exp V R  V  R 
 V  R exp V R / exp V R  V  R 
/ exp V R 
if V R   0
k/ exp V R  

f  N exp  Vb R N

Ab

N

b

Now kTST
So

kTST

and

kTST

or

kTST

Ab

*

N

A

*

N

N

N

b

Ab

*

N

Ab
TST

b

N

b

Ab

N

b

Ab

Ab

*

Ab

b

Temperature Accelerated Dynamics
Autumn School, Sofia
Bulgaria 2012

First order reactions:

P dt = k exp(-kt) dt
Lifetime of state: =1/k

log(1/)Hopping probability:

increasing
simulation
time

p1


p2

Arrhenius:

E2

k = A exp(-E*/RT)log(1/) = log A - E*/RT

stop time tend

1/RToo

1/RTh

E1
1/RTl

Temperature Accelerated Dynamics 2
Autumn School, Sofia
Bulgaria 2012

• Simulate system at high T & watch for transitions• When transition found, stop simulation and:– Determine activation energy using nudged elastic
band– Record transition time, save `new’ state
configuration• Restart simulation in original state with new velocities.
• Search for new transitions. Hence build `library’ oftransition data.
• Stop searching after time tend given by:
tend=exp[E2+(Th-Tl)(E2-)/Th]• Commence new search from `first’ low T state.

Nudged Elastic Band
Autumn School, Sofia
Bulgaria 2012

E

A
C0

C2
C1

C3

C4
CN-1

B
CN
R

A

B

• N+1 configs (C0…CN) linearlyinterpolated From A to B
• Connect by spring (stiffness K)• Remove `off tangent’ forces• Minimise all configs subject topresence of spring forces• Resulting path is reaction paththrough saddle point

Kinetic Monte Carlo
Autumn School, Sofia
Bulgaria 2012

pi ; i  1,, N • Simulate set of competing processes• Rate of process p i is ri (make a list).
N• Define sum of ratesR   ri

• Generate random number• Select processi 1

u : 0  u 1

 i 1  i pi :  rj   uR   rj 
 j 1  j 1 R• Advance time
• Repeat!

uRt   log(u) / Rpi

Invoking the Hyperdynamics Options
Autumn School, Sofia
Bulgaria 2012

In the CONTROL file:tadunits kJ
num_block500
num_track10
blackout
1000
catch_radius 3.5
neb_spring 10.0
deltad
6.91
low_temp
40.0
force
0.0025
endtadOR

bpd path
units eVvmin -3.9935E03
ebias -3.5000E03
num_block 300
num_track 10
catch_radius 3.5
neb_spring 1.0
force
0.00025
endbpdHyperdynamics Files
Autumn School, Sofia
Bulgaria 2012

Additional files for TAD and Bias potential
dynamics:• HYPRES/HYPOLD – restart files• EVENTS –program activity report• CFGBSNnn – Basin CONFIG files (new states)• PROnn.XY – Reaction path profiles• CFGTRAnn – Tracking CONFIG files

Subdirectories required in execute directory:
BASINS, PROFILES, TRACKSTAD – DL_POLY Test Case 32
Autumn School, Sofia
Bulgaria 2012

255 L-J Argon atoms FCC crystal + 1 vacancy• Atoms `hop’ into vacancies• Each vacancy has 12 nearest neighbour
atoms• So 12 possible escapes from PE basin• Use TAD to find them!
• Use NEB to find activation energy• Extrapolate to low temperature for low
T rate• Put results into KMC simulation

EVENTS file extract:
Event
nt
BasinsTRA
38500 0 1
TRA55500 0 2
TRA 127500 0 3
TRA 750500 0 4

Nt
E
1 7.28338E+00
1 7.20808E+00
1 7.28160E+00
1 7.19597E+00

Time(ps)
3.82250E+01
5.49650E+01
1.26145E+02
7.47515E+02

Extrap.(ps)4.31244E+07
5.36891E+07
1.41830E+08
7.13444E+08

Stop time(ps)2.04398E+03
2.04398E+03
2.04398E+03
2.04398E+03

BPD – DL_POLY Test Case 33
Autumn School, Sofia
Bulgaria 2012

998 NaCl ions rocksalt crystal + 2 vacancies• Overall neutral system• Ions `hop’ into vacancies• Escapes from PE basin unknown (apriori)• Use BPD to find them!
• Use NEB to find activation energy• Extrapolate hopping time for zerobias
• Put results into KMC simulation

EVENTS file extract:
Event
nt
TRA
4500
TRA 399300
TRA 466500

Basins Nt
E
Time(ps)
Extrap.(ps)0 1 1 6.74301E-01 4.39500E+00 7.34793E+03
1 2 1 1.11127E+00 3.99185E+02 6.45155E+05
2 3 1 6.57466E-01 4.66495E+02 7.53837E+05

NEB Reaction Profiles
Autumn School, Sofia
Bulgaria 2012

Lennard Jones Argon

Sodium Chloride

Metadynamics
Autumn School, Sofia
Bulgaria 2012

Metadynamics is a method devised by Alessandro Laioand Michele Parrinello for accelerating the explorationof a free energy landscape as the function of collective
variables.
Method:• The system potential energy is augmented by a timedependent bias potential consisting of Gaussianfunctions of the collective variables• The longer a simulation remains in a particular freeenergy minimum, the larger the bias potential becomes– thus forcing the system to seek out a new
thermodynamic state.
• The accumulated bias potential provides a descriptionof the free energy surface

Metadynamics in 1D
Autumn School, Sofia
Bulgaria 2012

A. Laio & M. Parrinello, PNAS 99 (2002) 12562

Collective Variables?
Autumn School, Sofia
Bulgaria 2012

A collective variable is a single number that definesN

an atomic structure (i.e. it is a function of r ). Most often
they are called Order Parameters. Particular examples
used in metadynamics are:
• The system potential energy:• Simulation cell vectors:N

U (r )Q

• The Steinhardt order parameters:h  (a, b, c)• Tetrahedral order parameters:Q and   are maximum for particular structures.
Defining the bias potential in terms of order parameters
allows destabilization of particular structural phases.

Metadynamics Formulae
Autumn School, Sofia
Bulgaria 2012

Order parameter vector: s (r )  s1 (r ),, sM (r )
M

N

N

N

pi2
N
M
N
System Hamiltonian: H   2m  U (r )  V [ s (r ), t ]
i 1
i
N

Ng

2
M
M
2

V
[
s
(
r
),
t
]

W
exp

s
(

)

s
(
t
)
/
2

h
Bias Potential:

k



k 1
M

n

W and h are chosen to `fill’ surface at acceptable rate
Force:

VN
f i  iU (r )   i s j (r )j 1 s j
N

MFree Energy Surface: F ( s M )   lim V [ s M (r N ), t ]gt Metadynamics Files
Autumn School, Sofia
Bulgaria 2012

• METADYNAMICS– Data defining the metadynamics hypersurface• STEINHARDT– Defines the Steinhardt order parameters• ZETA– Defines the tetrahedral order parameters

Steinhardt Order Parameters
Autumn School, Sofia
Bulgaria 2012

 4 

1
Q 
Q m

 2  1 m  N C Nfor atom types  and  , with


2 1/ 2Nb

Q m   f c (rb )Ym  b , b 


b 1

where b runs over all N b vectors connecting atom types  and  and
if r  r1 
 1
 1   (r  r )  1f c (r )   cos   1 if r1  r  r2 
 2   (r2  r1 )   0
if r  r2 Tetrahedral Order Parameters
Autumn School, Sofia
Bulgaria 2012

N N N1
2
T 
f c (rij ) f c (rik )(cos jik  1 / 3)

N c N i 1 j i k  j
Where i, j , k run over all N atoms of species
and N c is the number of pairs of atoms linked to
atom  (assuming all atoms are of type  )Ice Nucleation and Growth 1
Autumn School, Sofia
Bulgaria 2012

Bias:
Q4OO, Q6OO, T & PE

0.5ns
D Quigley and PM Rodger, Molec. Sim. 35 (2009) 613

Ice Nucleation and Growth 2
Autumn School, Sofia
Bulgaria 2012

Bias:
Q4OO, Q6OO, T & PE

0.75nsD Quigley and PM Rodger, Molec. Sim. 35 (2009) 613

Ice Nucleation and Growth 3
Autumn School, Sofia
Bulgaria 2012

Bias:
Q4OO, Q6OO, T & PE

1.25nsD Quigley and PM Rodger, Molec. Sim. 35 (2009) 613

Ice Nucleation and Growth 4
Autumn School, Sofia
Bulgaria 2012

Bias:
Q4OO, Q6OO, T & PE

1.5ns
D Quigley and PM Rodger, Molec. Sim. 35 (2009) 613

Conclusions
Autumn School, Sofia
Bulgaria 2012

• DL_POLY Classic is free
• It's very versatile with advanced features• Go get it!

Autumn School, Sofia
Bulgaria 2012

Part 7
The DL_POLY Java GUI
W. SmithGUI Overview
Autumn School, Sofia
Bulgaria 2012

• Java - Free!
• Facilitate use of code
• Selection of options (control of capability)• Construct (model) input files• Control of job submission• Analysis of output• Portable and easily extended by user

Compiling/Editing the GUI
Autumn School, Sofia
Bulgaria 2012

• Edit source in java directory• Edit using vi,emacs, whatever• Compile in java directory:
javac *.javajar -cfm GUI.jar manifesto *.class• Executable is GUI.jar• But.....
****Don't Panic!****The GUI.jar file is provided in the download

Invoking the GUI
Autumn School, Sofia
Bulgaria 2012

• Invoke the GUI from within the execute directory (orequivalent):
java -jar ../java/GUI.jar• Colour scheme options:java -jar ../java/GUI.jar –colourschemewith colourscheme one of:monet, vangoch, picasso, cezanne, mondrian (default
picasso).

The Monitor Window
Autumn School, Sofia
Bulgaria 2012

Menus

Using Menus
Autumn School, Sofia
Bulgaria 2012Show
EditorOptionThe Molecular Viewer
Autumn School, Sofia
Bulgaria 2012

Editor
ButtonGraphics
Buttons

Graphics
Window

The Molecular Editor
Autumn School, Sofia
Bulgaria 2012

Editor
Buttons

EditorWindow

Available Menus
Autumn School, Sofia
Bulgaria 2012

• File - Simple file manipulation, exit etc.
• FileMaker - make input files:– CONTROL, FIELD, CONFIG, TABLE• Execute– Select/store input files, run job• Analysis– Static, dynamic,statistics,viewing,plotting• Information– Licence, Force Field files, disclaimers etc.

A Typical GUI Panel
Autumn School, Sofia
Bulgaria 2012

Buttons

Text BoxesDL_POLY & VMD
Autumn School, Sofia
Bulgaria 2012

• VMD is a free software package for visualising MD
data.
• Website: http://www.ks.uiuc.edu/Research/vmd/• Useful for viewing snapshots and movies.
– A plug in is available for DL_POLY HISTORY files– Otherwise convert HISTORY to XYZ or PDB format

“Hands-On Session”
Autumn School, Sofia
Bulgaria 2012

This will consist of (up to) five components:• A demonstration of the Java GUI• Download & compile DL_POLY classic• Trying some DL_POLY simulations:– prepared exercises, or– creative play

• DL_POLY clinic - what’s up doc?
• Group therapy – all for one and one for all....

Autumn School, Sofia
Bulgaria 2012

Part 9DL_FIELD
C.W. YongDL_FIELD
Autumn School, Sofia
Bulgaria 2012

• AMBER & CHARM to DL_POLY Input Format• OPLS-AA & Drieding to DL_POLY Input Format

xyz
PDB
Protonated

DL_FIELD‘black box’

FIELD

CONFIGDL_POLY & VMD
Autumn School, Sofia
Bulgaria 2012•••

VMD is a free software package for visualising MD data
Website: http://www.ks.uiuc.edu/Research/vmd/
Useful for viewing snapshots and movies.
– A plug in is available for DL_POLY HISTORY files– Otherwise convert HISTORY to XYZ or PDB format

Autumn School, Sofia
Bulgaria 2012

Part 10
DL_MESO
M. Seaton

Autumn School, Sofia
Bulgaria 2012

Part 11
DL_POLY Hands-On
http://www.ccp5.ac.uk/DL_POLY/TUTORIAL/EXERCISES/index.html

The Ewald Summation
Autumn School, Sofia
Bulgaria 2012

Uc 

2

 exp(  k 2 / 4 2 ) N

12V o k  0 q j exp ik  rj  

k2


j 1

 
qn q j  erfc  R  rjn  
4 o R  0 n  j R  rjn
1



N



N

 q 2j4 3 / 2  o j 1

with:

 2
k  1/ 3 (, m, n) V

Smoothed Particle-Mesh Ewald
Autumn School, Sofia
Bulgaria 2012

Ref: Essmann et al., J. Chem. Phys. (1995) 103 8577

The crucial part of the SPME method is the conversionof the Reciprocal Space component of the Ewald sum
into a form suitable for Fast Fourier Transforms (FFT).
Thus:
2





1  exp( k 2 / 4 2 ) N
U recip 
q j exp  ik  rj


2
 
2V o k 0
kj 1

becomes:

1
T
U recip G(k1 , k2 , k3 )Q(k1 , k2 , k3 )2V o k1 ,k2 ,k3

where G and Q are 3D grid arrays (see later)

SPME: Spline Scheme
Autumn School, Sofia
Bulgaria 2012

Central idea - share discrete charges on 3D grid:
Cardinal B-Splines Mn(u) - in 1D:

exp2iu j k / L   b(k )  M n (u j  ) exp2ik / K   1
b(k )  exp( 2i (n  1)k / K )  M n (  1) exp2ik / K   0n
1n!
k
M n (u ) (1)max( u  k ,0) n 1(n  1)! k 0
k!(n  k )!
unuRecursion
M n (u ) 
M n 1 (u ) 
M n 1 (u  1)
n 1
n 1
relationn2

SPME: Building the Arrays
Autumn School, Sofia
Bulgaria 2012

Q ( 1 ,  2 ,  3 ) N q  M (u    n K )M (u    n K )M (u    n K )
j

j 1

n

1j

1

1

1

n

2j

2

2

2

n

3j

3

3

3

n1 , n2 , n3

Is the charge array and QT(k1,k2,k3) its discrete Fourier transform.
GT (k1,k2,k3) is the discrete Fourier Transform of thefunction:exp( k 2 / 4 2 )T*
G(k1 , k2 , k3 ) 
B
(
k
,
k
,
k
)(
Q
(
k
,
k
,
k
))
1
2
3
1
2
3
k2

with

B (k1 , k 2 , k3 )  b1 (k1 ) b2 (k 2 ) b3 (k3 )222

SPME: Comments
Autumn School, Sofia
Bulgaria 2012

• SPME is generally faster then conventional Ewald sum inmost applications. Algorithm scales as O(NlogN)• In DL_POLY_Classic the FFT array is built in pieces oneach processor and made whole by a global sum for the FFT
operation• In DL_POLY_4 the FFT array is built in pieces on eachprocessor and kept that way for the distributed FFToperation (DaFT)• The DaFT FFT `hides’ all the implicit communications

Parallel FFTs - The Basics
Autumn School, Sofia
Bulgaria 2012

FFTs are• Fast (!) - O(VlogV) operations where V is the number ofpoints in the grid
• Global operations - to perform a FFT you need all the
points• This makes it difficult to write an efficient, good scaling
FFTTraditional Parallel FFTs
Autumn School, Sofia
Bulgaria 2012

• Distribute the data by planes• Each processor has a complete set of points in the x and y
directions so can do those Fourier transforms• Redistribute data so that a processor holds all the points in z
• Do the z transforms• Allows efficient implementation of the serial FFTs ( use alibrary routine )• In practice for large enough 3D FFTs can scale reasonably• However, the distribution does not map onto DL_POLY’s large amounts of data redistribution

Daresbury advanced Fourier Transforms
Autumn School, Sofia
Bulgaria 2012

• Takes data distributed as in DL_POLY_4 (Domain
Decomposition) thus avoiding a lot of communication incomparison to traditional 3D FFT routines (due to data
redistribution) when used on more than 8-64 CPUs

• So do a distributed data FFT in the x direction then the yand finally the z directions• Disadvantage is that can not use the library routine for the
1D FFT ( not quite true … )• Scales quite well - e.g. on 512 procs, an 8x8x8 CPU grid, a
1D FFT need only scale to 8 CPU• Totally avoids data redistribution, exists as a stand alone
F90 library.

RD Scheme for long-ranged part of SPME
Autumn School, Sofia
Bulgaria 2012

U. Essmann, L. Perera, M.L. Berkowtz, T. Darden, H. Lee, L.G.
Pedersen, J. Chem. Phys., (1995), 103, 8577
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.

Calculate self interaction correction
Initialise FFT routine (FFT – 3D FFT)
Calculate B-spline coefficients
Convert atomic coordinates to scaled fractional units
Construct B-splines
Construct charge array Q
Calculate FFT of Q array
Construct array G
Calculate FFT of G array
Calculate net Coulombic energy
Calculate atomic forces

DD Scheme for long-ranged part of SPME
Autumn School, Sofia
Bulgaria 2012

U. Essmann, L. Perera, M.L. Berkowtz, T. Darden, H. Lee, L.G.
Pedersen, J. Chem. Phys., 103, 8577 (1995)
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.

Calculate self interaction correction
Initialise FFT routine (FFT – IJB’s DaFT: 3M2 1D FFT)
Calculate B-spline coefficients
Convert atomic coordinates to scaled fractional units
Construct B-splines
Construct partial charge array Q
Calculate FFT of Q array
Construct partial array G
Calculate FFT of G array
Calculate net Coulombic energy
Calculate atomic forces

I.J. Bush, I.T. Todorov, W. Smith, Comp. Phys. Commun., 175, 323(2006)

Choosing Ewald Sum Parameters
Autumn School, Sofia
Bulgaria 2012

Simplest option:ewald precision [] (tolerance)
Advanced option (for anisotropic systems,
particularly with slabs):
ewald [] [kmax1] [kmax2] [kmax3]with = exp[-( rcut)2]/rcut (solve for ) = exp[-(km/2)2]km2 (solve for km)
kmax n=Lnkm/2 , where Ln is the MD box width inn direction