Development of GENESIS
on K for large scale
molecular dynamics simulation
Jaewoon Jung(RIKEN Advanced Institute for Computational Science)Importance of molecular dynamics
1. MD is a sampling method : We need ensembles because proteins constantly sample
configuration space

2.

Unlike other sampling method (i.e. Monte Carlo), MD gives dynamics

3. In particular, MD helps us to understand1) Kinetics – chemical reaction (with QM), conformational changes, andprotein folding…
2) Equilibrium Ensemble sampling – thermodynamics (free energy changes,bindings,..)3) Modeling is possible – structure prediction, solvent effects…)

For theses purpose, we need accurate potential from withefficient computational algorithms

Fast non-bonded interaction- New lookup table approachMolecular Mechanics (Force field)
1. The basic functional form of a force field encapsulates both bonded terms relating toatoms that are linked by covalent bonds, and non-bonded (also called "noncovalent")
terms describing the long-range electrostatic and van der Waals forces

2. A force field defines a set of parameters for each type of atom.
∑∑∑∑∑

1cos2

Parameters , , , , , , ,, and are definedaccording to each atom type.
3. The bottleneck of the computation is the non-bond energy calculation due to O(N2)computation steps for N particle.

How to calculate non-bond energy efficiently?
4

Cutoff for non-bond energy calculation
1. Non-bond energy calculation is reduced by introducing cutoff
2

2

O(N2)
4

O(N1)

U

2. The electrostatic energy calculation beyond cutoff will be done in the reciprocal spacewith FFTerfc

U||

4Real part2expGG ⁄4G

4

Reciprocal partcos G ∙4

Self energy

3. Further, it could be reduced by properly distributing over parallel processors, in particular
good domain decomposition scheme.
5Why do we need lookup table?
1. The short-range non-bonded interactions are the main computationalbottleneck in molecular dynamics simulations
2. When we are using PME for electrostatic interactions, short-range nonbonded interaction includes time-consuming inverse square root andcomplementary error functions
3. Such time-consuming operations could be avoided by lookup table approach
example : Short-range non-bonded interaction takes 1.23 sec/step without
using lookup table, and 0.45 sec/step with lookup table for DHFR

How do we change table points ?
(comparison of linear R2 table and linear 1/R2 table)
0.10

0.05

0.05

0.00

0.00

-0.05

-0.05

12

-0.10

121/r -1/r

61/r -1/r

6

0.10

-0.15

-0.10
-0.15

-0.20

-0.20

-0.25

-0.25

-0.30

12

3

4

5

62r

7

8

9

10

-0.30
0.1

0.2

0.3

0.4

0.5

0.6
21/r

0.7

0.8

0.9

1.0

1. Assign more table points at small r2 where energy and gradient change rapidly
2. For large r2 where energy and gradient does not change so much, we can reduce thetable points for faster evaluation.

1E-3

2

Linear R
2
Linear 1/R

1E-4

1E-5

1E-6
0

10

2030

405060708090 100 110

Relative Energy error (kcal/mol)

Relative Energy error (kcal/mol)

Accuracy of single point energy calculation(DHFR)
1E-4

2

Linear R
2
Linear 1/R

1E-5
1E-6
1E-7
1E-8
1E-9

01030

4050

60

708090 100 110

Table densityTable density

Van der Waals∆20Electrostatic/

1. Van der Waals : Linear 1/R2 with table density = 20 is more accurate than linear R2with table density = 100
2. Electrostatic : Linear 1/R2 with table density = 10 is more accurate than linear R2 withtable density = 100

Energy drift (DHFR)-57936-57925

Linear R2 withdensity = 10-57930-57935

0200

400

600

800

Total Energy (kcal/mol)

Total Energy (kcal/mol)-57920-57937-57938-57939-579411000

Linear 1/R2 withdensity = 10-57940

0200-57935-57934-57936-57935

Linear R2 withdensity = 20-57937-579380200400

600

Time (ps)

800

1000

Total Energy (kcal/mol)

Total Energy (kcal/mol)-57933-57936

400

600

800

1000

Time (ps)Time (ps)Linear 1/R2 withdensity = 20-57937-57938-57939-57940

0200

400

600

800

1000

Time (ps)1. Large energy drift for both table density = 10 and 20 using linear R2 table.
2. Slight amount of energy drift for table density = 10 and no energy drift from tabledensity = 20 using linear 1/R2 table.

Energy drift of DHFR according to table density

Energy drift (kcal/mol/ns)

1082

Linear of r
2Linear of 1/r

64
200102030405060708090100 110

Table density

1. Linear R2 table : reliable energy drift from table density = 40.
2. Linear 1/R2 table : reliable energy drift from table density = 20

Non-bonded interaction time of DHFR
(resources : K)

Time per step (sec)

0.60

0.55

0.50

0.45

0.40

2

Linear of r
2Linear of 1/r
0102030405060708090100 110Table density

On K, linear 1/ R2 table shows better performance, due to efficient inverse operationand cache usage.

Why linear 1/ R2 table is efficient for large tabledensity?
Interval

Number ofinteractions

Number of sections inlinear R2 table

Number of sections inlinear 1/R2 table

1<r2<16269245

15

135

16<r2<36730108

20

5

36<r2<1447295834

108

3Less cache misses
1. Using Linear 1/R2 table, 90% interactions are dealt with only three sections andincreasing table density does not introduce cache misses for this interval.
2. On the other hand, linear R2 table uses 108 sections for the same interval andcache misses will be severe by increasing table densities.
3. In K, this phenomena highly affect SIMD and this explains higher speed of linear
1/ R2 table even at small table density

GENESIS parallelization scheme- Midpoint cell method

Domain Decomposition
012

3

4

5

6

7

8

9

10

111213 14

15

16

17

18

1920

21 22

2324

25

26

27

28

29 30

31

32

33

34

35

36

37 38

39

40

41

42

43

44

45 46

47

48

49

50

51

52

53 54

55

56

57

58

59

60

61 6263

1. Partitioning space into
fixed-size boxes, withdimension larger than thecutoff distance
2. We need only informationof neighborspace(domain) forcomputation of energies
3. Communication is
reduced by increasing
process number
4. Efficient for goodparallelization andsuitable for large systemwith large processors

Basic scheme of parallelization Midpoint method
1. Two particles interact on a particular box if and only if the midpoint of the segment
connecting them falls within the region of space associated with that box

Each pair of particles separated by a distance less than R (cutoff distance) is connected by a dashed line segment,with “x” at its center lying in the box which will compute the interaction of that pair

2. This scheme applies not only for non-bonded but also bonded interactions.
Ref : KJ. Bowers et al, JCP 124, 184109 (2006)

Midpoint cell method

Midpoint method : interaction between
two particles are decided from the
midpoint position of them.

Midpoint cell method : interaction
between two particles are decided fromthe midpoint cells where each particle
resides.

Small communication, efficient energy/force evaluations

Communication spaceHybrid parallelization in GENESIS
1. The basic is the Domain
decomposition with the midpoint
scheme2. Each node consists of at least two
cells in each direction.
3. Within nodes, thread calculation(OPEN MP) is used forparallelization and communication
is not necessary
4. Only for different nodes, we allow
point to point communication(MPI) for neighboring domains.

Efficient shared memory calculation=>Cell-wise particle data
1. Each cell contains an array with the data of the particles that reside within it
2. This improves the locality of the particle data for operations on individual cells orpairs of cells.

particle data in traditional cell lists

cell-wise arrays of particle data
Ref : P. Gonnet, JCC 33, 76-81 (2012)

How to parallelize by OPENMP?
1. In the case of integrator, every cell indices are divided according to the thread id.
2. As for the non-bond interaction, cell-pair lists are first identified and cell-pair lists are
distributed to each thread

12

3

4

5

6

7

8

9

10

11

12

13

14

15

16

(1,2)
(1,5)
(1,6)
(2,3)
(2,5)
(2,6)
(2,7)
(3,4)
(3,6)
(3,7)
(3,8)
(4,7)
(4,8)
(5,6)
(5,9)
(5,10)(6,7)(6,10)…

thread1
thread2
thread3
thread4
thread1
thread2
thread3
thread4…

Parallelization result(K with 8 OPENMP
Cutoff = 12.0 Å, Particle Mesh Ewald)

Time (sec) per step of 1M System (STMV)Core Number

NAMD

GENESIS256

0.16080.1584512

0.09080.08241024

0.05800.04432048

0.0300

0.024

4096

0.02430.01638192

0.02600.008716384

0.005832768

0.0047

1. Here we used NAMD compiled on K (provided by Dr. Miyashita)
2. Load balance is not considered for NAMD. If load balance option is on, it isfaster than w/o load balance up to 2048 cores, but it becomes slower as we
increase the process numbers
3. Grid size of namd input : 216x216x216, Grid size of GENESIS : 256x256x256

Overall performance for 1.3 M system
3216408204108624

1256

512

1024

2048

4096

8192

16384 32768

Number of cores

4Time (msec/step)

1008060

FFT Forward
FFT Forward communication
FFT Backward
FFT Backward communication
Send/Receive before force
Send/Receive after force

8Simulation time (ns/day)

Simulation time (msec / step)

2002
10.5
0.25
0.125

2565121024

2048

4096

8192 16384 32768

Number of cores

1. It is scalable up to 4,096 nodes (32,768 cores)
2. Not only energy/force evaluation, but also communications for FFT and send/receive
are scalable even though the scaling efficiency is not perfect

16

100
80

860

440
2
20
1
108

2048

4096

8192

16384

32768

Number of cores

65536

0.5131072

FFT Forward
FFT Forward communication
FFT Backward
FFT Backward communication
Send/Receive before force
Send/Receive after force

84Time (msec/step)200

Simulation time (ns/day)

Simulation time (msec / step)

System : 11.2 M atoms (time step = 2fs)

210.5
0.25
0.125

256

512

1024

2048

4096

8192 16384 32768

Number of cores

System : 100 M atoms

Simulation time (ns/day)84

GENESIS on K
NAMD on Blue Gene/Q
NAMD on XT5
NAMD on XK6

2

10.5
16384

32768

65536

131072

Number of cores

262144

Summary•

GENESIS uses newly developed lookup table for more accurate and fastevaluation of non-bonded interaction•

GENESIS uses the midpoint cell method for our domain decomposition
scheme of parallelization.

•Hybrid parallelization (MPI+OPENMP) is considered for efficientparallelization.

•GENESIS is scalable up to the number of cells, in particular, 131,000 coresfor 11.7M systems (more than 16ns/day) and 262,cooo cores for 100M
systems.