計算科学技術特論ＡParallelization ofmolecular dynamics

2017年 7月 13日
Jaewoon Jung
(RIKEN Advanced Institute for Computational Science)Overview of MD

Molecular Dynamics (MD)
1. Energy/forces are described by classical molecular mechanics force field.
2. Update state according to equations of motion

dri pi

dt
m
dp i
 Fi
dt

pi
ri (t  t )  ri (t )  tmpi (t  t )  pi (t )  Fi t

Equation of motion

Integration

Long time MD trajectory=> Ensemble generation

Long time MD trajectories are important to obtainthermodynamic quantities of target systems.

Potential energy in MD
Etotal  kb (b  b0 )2

O(N)

bonds



Total number of particles

 ka (  0 )2

O(N)angles Vn [1  cos(n   )]O(N)dihedrals

12
6
 




r
r
q
q
0ij
0ij
i j

    ij   2  rij  rij  rij 
j 1 i  j 1  







N 1 NO(N2)Main bottleneck in MD

12
6
 

2
2



r
r
q
q
erfc(

r
)
exp(

k
/
4

)
0ij
i
j
ij
  0ij




2

FFT(Q(k ))




 
  ij  r 
2



rij 
rij
k
 k 0i  j  R   ij 






Real space, O(CN)

Reciprocal space, O(NlogN)4

Non-bonded interaction
1. Non-bond energy calculation is reduced by introducing cutoff
O(N2)

O(N1)

2. The electrostatic energy calculation beyond cutoff will be done in the reciprocal
space with FFT||Real part

Reciprocal part

Self energy

3. Further, it could be reduced by properly distributing over parallel processors, inparticular good domain decomposition scheme.
5Difficulty to perform long time MD simulation
1. One time step length (Δt) is limited to 1-2 fs due to vibrations.
2. On the other hand, biologically meaningful events occur on the time scale ofmilliseconds or longer.
fs

ps

ns

μs

ms

sec

vibrations
Sidechain motions
Mainchain motions
Folding
Protein global motions

How to accelerate MD simulations?
=> Parallelization
SerialParallel

X16?
1cpu
CPU
Core
C

1cpu

16 cpusMPICommGood Parallelization ;1) Small amount ofcomputation in one CPU
2) Small amount ofcommunication time

Parallelization

Shared memory parallelization (OpenMP)
Memory

P1

P2

P3

PP-1

PP

• All processors share data in memory• For efficient parallelization, processors should not access the same memory address.
• It is only available for multi-processors in a physical nodeDistributed memory parallelization (MPI)
M1

M2

M3

MP-1

MP

PP-1

PP

…
P1

P2

P3

• Processors do not share data in memory• We need to send/receive data via communications• For efficient parallelization, the amount of communication data should be minimized

Hybrid parallelization (MPI+OpenMP)
M1

M2

M3

MP-1

MP

PP-1

PP

…
P1

P2

P3

• Combination of shared memory and distributed memory parallelization.
• It is useful for minimizing communicational cost with very large number of processors

SIMD (Single instruction, multiple data)

SIMD

No SIMD
4 times faster• Same operation on multiple data points simultaneously• Usually applicable to common tasks like adjusting graphic image or volume• In most MD programs, SIMD becomes the one of the important topics to increase theperformance

SIMT (Single instruction, multiple threads)• SIMD hardware with MIMD programming modelEx) if then .. end if => SIMD (×)
SIMT (○)• SIMT execution model is usually implement on GPUs and related with GPGPU(General Purpose computing on Graphics Processing Units)• Currently, CUDA allows 32 threads for SIMT (warp size =32)
Grid 1

Block (1,0)

Block
(0,0)

Block
(1,0)

Block
(0,1)

Block
(1,1)

Thread(0,0)

Thread(1,0)

Thread(2,0)

Thread(3,0)

Thread(0,1)

Thread(1,1)

Thread(2,1)

Thread(3,1)

Thread(0,2)

Thread
(1,2)

Thread(2,2)

Thread(3,2)

Parallelization of MD (real space)Parallelization scheme 1 :Replicated data approach
1.
2.

Each processor has a copy of all particle data.
Each processor works only part of the whole works by proper assign in doloops.
my_rank = MPI_Rank
proc = total MPI
do i = 1, N
do j = i+1, N
energy(i,j)
force(i,j)
end doend dodo i = my_rank+1, N, proc
do j = i+1, N
energy(i,j)
force(i,j)
end doend doMPI reduction (energy,force)Hybrid (MPI+OpenMP) parallelization of the
Replicated data approach
1.
2.

Works are distributed over MPI and OpenMP threads.
Parallelization is increased by reducing the number of MPIs involved incommunications.
my_rank = MPI_Rank

my_rank = MPI_Rank
proc = total MPI
do i = my_rank+1, N, proc
do j = i+1, N
energy(i,j)
force(i,j)
end doend doMPI reduction(energy,force)proc = total MPI
nthread = total OMP thread!$omp parallel
id = omp thread id
my_id = my_rank*nthread + id
do i = my_id+1,N,proc*nthread
do j = i+1, N
energy(i,j)
force(i,j)
end doend doOpenmp reduciton!$omp end parallel
MPI reduction (energy,force)

Advantage/Disadvantage of
the Replicated data approach
1. Advantage : easy to implement
2. Disadvantage
1) Parallel efficiency is not gooda) Load imbalanceb) Communication is not reduced by increasing the numberof processors2) Needs a lot of memoryb) Memory usage is independent of the number of processors

Parallelization scheme 2 :Domain decomposition
1.

The simulation space isdivided into subdomains
according to MPI

2.

Domain size is usuallyequal to or greater thanthe cutoff value.

Advantage / Disadvantage ofthe domain decomposition approach
1. Advantage
1) Parallel efficiency is very good compared to the replicateddata method due to small communicational cost2) The amount of memory is reduced by increasing the numberof processors
2. Disadvantage
1) Difficult to implement

Comparison of two parallelization scheme
Computation

Communication

MemoryReplicated data

O(N/P)O(N)

O(N)

Domaindecomposition

O(N/P)

O((N/P)2/3)

O(N/P)

Parallelization of MD(reciprocal space)Smooth particle mesh Ewald method1
E (r)  2 n'

qi q j erfc  ri  r j  nN

i , j 1ri r j n '  rc

ri  r j  nReal part 1

exp( 2k 2 /  2 ) N 2
2S (k ) qi22 V k  0
k i 1

Reciprocal part

Self energy

The structure factor in the reciprocal part is approximated as

S ( k1 , k 2 , k3 )  b1 ( k1 )b2 ( k 2 )b3 ( k3 ) F (Q )(k1 , k 2 , k3 )Using Cardinal B-splines of order nFourier Transform of
chargeIt is important to parallelize the Fast Fourier transform efficiently in PME!!
Ref : U. Essmann et al, J. Chem. Phys. 103, 8577 (1995)

Parallel 3D FFT – slab(1D) decomposition
1. Each processor is assigned a slob of size
N × N × N/P for computing an N × N
× N FFT on P processors.
2. The parallel scheme of FFTW
3. Even though it is easy to implement, thescalability is limited by N, the extent of
the data along a single axis
4. In the case of FFTW, N should be
divisible by P

1D decomposition of 3D FFT

Reference: H. Jagode. Master’s thesis, The University of Edinburgh, 2005

1D decomposition of 3D FFT (continued)
1. Slab decomposition of 3D FFT has three steps
 2D FFT (or two 1D FFT) along the two local dimension
 Global transpose
 1D FFT along third dimension
2. Advantage : The fastest on limited number of processors because it only needs
one global transpose
3. Disadvantage : Maximum parallelization is limited to the length of the largest
axis of the 3D data ( The maximum parallelization can be increased by using a
hybrid method combining 1D decomposition with a thread based
parallelization)

Parallel 3D FFT –2D decomposition

1. Each processor is assigned a slob of size N × N/P × N/Q for computing an N× N × N FFT on P ×Q processors.
2. Current GENESIS adopt this scheme with 1D FFTW

2D decomposition of 3D FFT

Reference: H. Jagode. Master’s thesis, The University of Edinburgh, 2005

2D decomposition of 3D FFT (continued)
1.







2D decomposition of 3D FFT has five steps
1D FFT along the local dimension
Global transpose
1D FFT along the second dimension
Global transpose
1D FFT along the third dimension
Global transpose

2. The global transpose requires communication only between subgroups of all nodes
3. Disadvantage : Slower than 1D decomposition for a number of processors possible with
1D decomposition
4. Advantage : Maximum parallelization is increased
5. Program with this scheme Parallel FFT package by Steve Plimpoton (Using MPI_Send and MPI_Irecv)[1] FFTE by Daisuke Takahashi (Using MPI_AlltoAll)[2] P3DFFT by Dmitry Pekurovsky (Using MPI_Alltoallv)[3][1] http://www.sandia.gov/~sjplimp/docs/fft/README.html[2] http://www.ffte.jp[3] http://www.sdsc.edu/us/resources/p3dfft.php

2D decomposition of 3D FFT (pseudo-code)! compute Q factor
do i = 1, natom/P
compute Q_orig
end docall mpi_alltoall(Q_orig, Q_new, …)
accumulate Q from Q_new!FFT : F(Q)do iz = 1, zgrid(local)
do iy = 1, ygrid(local)
work_local = Q(my_rank)
call fftw(work_local)
Q(my_rank) = work_local
end doend docall mpi_alltoall(Q, Q_new,…)
do iz = 1, zgrid(local)
do ix = 1, xgrid(local)
work_local = Q_new(my_rank)
call fftw(work_local)
Q(my_rank) = work_local
end doend docall mpi_alltoall(Q,Q_new,..)

do iy = 1, ygrid(local)
do ix = 1, xgrid(local)
work_local = Q(my_rank)
call fftw(work_local)
Q(my_rank) = work_local
end doend dodo iy = 1, ygrid(local)
do ix = 1, xgrid(local)
work_local = Q(my_rank)
call fftw(work_local)
Q(my_rank) = work_local
end doend docall mpi_alltoall(Q,Q_new,..)! compute energy and virial
do iz = 1, zgrid(local)
do iz = 1, zgrid
do ix = 1, xgrid(local)
do iy = 1, ygrid(local)
work_local = Q(my_rank)
do ix = 1, xgrid(local)
call fftw(work_local)
energy = energy + sum(Th*Q)
Q(my_rank) = work_local
virial = viral + ..
end doend doend doend docall mpi_alltoall(Q,Q_new)
end dodo iz = 1, zgrid(local)
do iy = 1, ygrid(local)! X=F_1(Th)*F_1(Q)work_local = Q(my_rank)
call fftw(work_local)! FFT (F(X))
Q(my_rank) = work_local
end doend docompute force

Parallelization scheme of existing MD
programs

Gromacs
1. Gromacs makes use of multiple parallelism scheme

•
•
•
•

SIMD register : cluster interaction kernel or bonded interaction
OpenMP multithreading : inside domain decomposition
Nonbed interactions by GPU if GPU cards exist
MPI for each single domain
Ref : M. J. Abraham et al, Software X, 1-2, 19-25 (2015)

Gromacs
2. Gromacs makes use of the 8-th

shell scheme as the domain
decomposition scheme

Rc

Coordinates in zones 1 to 7 arecommunicated to the corner cell 0

8th shell scheme
Ref : B. Hess et al. J. Chem. Theor. Comput. 4, 435 (2008)

Gromacs
3. Gromacs assigns kernels according to available SIMD / SIMT widths

Ref : M. J. Abraham et al, Software X, 1-2, 19-25 (2015)

NAMD
1.
2.
3.

NAMD is based on the Charmm++ parallel program system and runtime
library.
Subdomains named patch are decided according to MPIForces are calculated by independent compute objects

Ref :L. Kale et al. J. Comput. Phys. 151, 283 ((1999)
J. C. Phillips et al. J. Comput. Chem. 26, 1781
(2005)Desmond
1. Midpoint method•

Two particles interact on a particular box if and only if the midpoint of thesegment connecting them falls within the region of space associated with thatbox

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

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

Desmond
2. FFT scheme••

Not distribute the data (not to use all-to-all communications)
Parallelization of one dimensional FFT by performing directly the butterfly
communications

Ref : KJ. Bowers et al, Proceedings of the 2006 ACM/IEEE conference on Supercomputing. ACM, 2006

GENESIS
1. Midpoint cell method• Partitioning space into fixed size
boxes, with dimension larger thanthe cutoff distance.
• We need only information ofneighbor space(domain) forcomputation of energies.
• Communication is reduced byincreasing process number .
• Efficient for good parallelizationand suitable for large system withmassively parallel
supercomputers.

GENESIS
2. Volumetric decomposition
FFT• More communications thanexisting FFT•

MPI_Alltoall
communications only in one
dimensional space•

Reduce communicational
cost for large number ofprocessors

Ref) J. Jung et al. Comput. Phys. Comm.
200, 57-65 (2016)

FFT in GENESIS (2 dimensional view)
GENESIS(Identical domain
decomposition
between two space)

NAMD, Gromacs(Different domain
decomposition)

SIMD in GENESIS (developer version)• Array of Structure (AoS)
x1 y1

z1

x2

y2

z2

xN

yN zN

In GENESIS, it is expressed as coord_pbc(1:3,1:natom,1:ncell)• Structure of Array (SoA)
x1 x2

x3

x4

zN-2 zN-1 zN

In updated GENESIS source code for KNL, it will be expressed ascoord_pbc(1:natom,1:3,1:ncell)• For Haswell/Broadwell/KNL machines, SoA shows better performance
than AoS due to efficient vectorization (SIMD)

Further improvements in
Parallelization of MD :MPMD (Multiple Program andMultiple Data)

Multiple program/multiple data (MPMD) scheme
1. A subset of nodes are dedicated to real space and another set of nodes are dedicated toreciprocal space interactions.
2. Representative program : Gromacs, CHARMM.
3. Usually different domain decomposition of real space is different from that of reciprocal
space => reciprocal space nodes is used only for calculation (not integration).

MPMD scheme of Gromacs : 3D domain decomposition in real space combined with 2D
pencil decomposition in reciprocal space (ref : Bioinformatics, btt055 (2013)

What is the problem of existing MPMD scheme?
1. Different domain decomposition
between real and reciprocal space=> Large amount ofcommunication
2. Integration of reciprocal space is
not possible => RESPA cannot be
implemented.
3. Communication pattern between
real and reciprocal nodes isdependent on the number ofprocessors

To overcome the problem, wesuggest a new MPMD scheme
where real and reciprocal spacehave same domain decomposition.

MPMD scheme with multiple time step integrator

1. Both real and reciprocal space nodes are involved in integrator (it is possible by same
domain decomposition between two spaces).
2. Real space non-bonded interactions are divided into two subspaces.
3. When reciprocal space interaction is not necessary, real and reciprocal space nodes
are assigned to the evaluation of subspace 1 and subspace 2, respectively.

Flowchart of MPMD scheme

Conventional integrator

Multiple time
step integrator

Benchmark performance

Gromacs

Ref : S. Pall et al., International Conference on Exascale Applications and Software, EASC 2014

NAMD

NAMD 2.11 on TACC Stampede(PME every 3 steps)(From NAMD webpage)

NAMD 2.10(PME every 3 steps)(In proceedings of the 2014 International
Conference for High Performance
Computing, Networking, Storage, andAnalysis (SC14))

GENESIS 1.1

GENESIS 1.1 on PC clusters(1M system)

GENESIS 1.1 on K
(28M system)

GENESIS MPMD performance
ApoA1

STMV

ACRB

8_STMV

27_STMV

1. MPMD is not suitable for small number of processors.
2. MPMD could be a good solution for very large number of processors.
3. MPMD even increases the available number of processors.

GENESIS performance on KNL
32

Performance (ns/day)

16

Performance (ns/day)

8

GENESIS (VVER)
GENESIS (RESPA,2)
GENESIS (RESPA,4)

8
4
2
1
1

2

4

8

16

Number of Nodes

1 STMV (1M)

32

64

128

GENESIS (VVER)
GENESIS (RESPA,2)
GENESIS (RESPA,4)

4

2

1

0.5

64

128

256

512

Number of Nodes

64 STMV (66M)1024

Summary
1. Parallelization : Distributed memory (MPI), shared memory (OpenMP),hybrid (MPI+OpenMP), SIMD, SIMT, etc.
2. Real space parallelization => Mainly domain decomposition scheme
3. Reciprocal space parallelization => Parallelization of FFT
4. Further improvements : MPMD, GPU, and so on.