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

2019年 7月 11日
Jaewoon Jung
(RIKEN Center 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

Equation ofmotion

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 )2bonds

 ka (  0 )2O(N)Total number of particles

O(N)angles Vn [1  cos(n   )]O(N)dihedrals126 rrqq
0ij
0iji j    ij   2  rij  rij  rij j 1 i  j 1  





N 1 NO(N2)Main bottleneck in MD

126
 
22

r
r
qq
erfc(r)exp(k/4)0ijijij  0ij2
FFT(Q(k ))   ij  r 2rij rijk 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 spacewith FFT||Real part

Reciprocal part

Self energy

3. Further, it could be reduced by properly distributing over parallel processes, in particular
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.

fsps

nsμs

ms

sec

vibrations
Sidechain motions
Mainchain motions
Folding
Protein global motions

How to accelerate MD simulations?
=> Parallelization
SerialParallel

X16?
1cpu
CPUCoreC

1cpu

16 cpus
MPICommGood Parallelization ;1) Small amount ofcomputation in one process
2) Small communicational cost

Parallelization

Shared memory parallelization (OpenMP)

MemoryP1

P2

P3PP-1PP• All processes share data in memory• For efficient parallelization, processes 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

MPPP-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 combined with multithreading• SIMT execution model is usually implemented 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)

MPI Parallelization of MD(non-bonded interaction in real space)Parallelization scheme 1: Replicated data approach

1.

Each process has a copy of all particle data.

2.

Each process works only part of the whole works by proper assign in do loops.
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)Parallelization scheme 1: Replicated data approach
atom indices
12

3

45678

MPI rank 0

Computational cost

MPI rank 1
MPI rank 2
MPI rank 3

Perfect load balance is not guaranteed in this parallelization schemeHybrid (MPI+OpenMP) parallelization of theReplicated 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
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)

my_rank = MPI_Rank
proc = total MPI
nthread = total OMP thread!$omp parallelid = 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)Pros and Cons of the Replicated data approach
1.

Pros : easy to implement

2.

Cons
•

Parallel efficiency is not good•No perfect load balance•Communication cost is not reduced by increasing the number ofprocesses•

We can parallelize only for energy calculation (with MPI,parallelization of integration is not so much efficient)•Needs a lot of memory•

Usage of global data

Parallelization scheme 2: Domain decomposition
1.

The simulation space is divided intosubdomains according to MPI(different colors for different MPIs).

2.

Each MPI only considers thecorresponding subdomain.

3.

MPI communications only amongneighboring processes.

Communications
between processors

Parallelization scheme 2: Domain decomposition

BufferSubdomain

1.

For the interaction between
different subdomains, it isnecessary to have the data of thebuffer region of each subdomain.

2.

The size of the buffer region isdependent on the cutoff values.

3.

The interaction between particles indifferent subdomains should beconsidered very carefully.

4.

The size of the subdomain andbuffer region is decreased byincreasing the number of processes.

𝒄Pros and Cons of the domain decomposition approach
1. Pros••

Good parallel efficiency•

Reduced computational cost by increasing the number of processes•

We can easily parallelize not only energy but also integration

Availability of a huge system•

Data size is decreased by increasing the number of processes

2. Cons•

Implementation is not easy•

Domain decomposition scheme is highly depend on the potential energy
type, cutoff and so on•

Good performance cannot be obtained for nonuniform particle distributions

Special treatment for nonuniform distributions
Tree method: Subdomain size is adjusted to have the same number of particles

Hilbert space filling curve: a map that relates multi-dimensional space to one-dimensional curve

The above figure is from Wikipedia
https://en.wikipedia.org/wiki/Hilbert_curve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)

N: system sizeP: number of processes

MPI Parallelization of MD(reciprocal space)Smooth particle mesh Ewald method
1E (r)  2 n'

qi q j erfc  ri  r j  nN

i , j 1ri r j n '  rcri  r j  nReal part 1exp( 2k 2 /  2 ) N 2
2S (k ) qi22 V k  0k
 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 ofcharge

It is important to parallelize the Fast Fourier transform efficiently in PME!!
Ref : U. Essmann et al, J. Chem. Phys. 103, 8577 (1995)

Overall procedure of reciprocal space calculation(charge on real space) PME_Pre(charge on grid)
FFT

E( )
InverseFFT
PME_Post

Force calculationfromSimple note of MPI_alltoall communications
Proc1 Proc2

Proc3 Proc4

Proc1 Proc2

Proc3 Proc4

A1

B1

C1

D1

A1

A2

A3

A4

A2

B2

C2

D2B1B2

B3

B4C1C2C3C4

D1

D2

D3

D4A3

B3C3

D3A4B4C4D4MPI_alltoall

MPI_alltoall is same as matrix transpose!!

Parallel 3D FFT – slab (1D) decomposition

1. Each process is assigned a slob of sizefor computing FFT ofgrids on

Proc 1
Proc 2

processes.

2. The scalability is limited by

Proc 3
Proc 43.

should be divisible byParallel 3D FFT – slab (1D) decomposition

No communication

MPI_Alltoall communication
among all preocessors

Parallel 3D FFT – slab (1D) decomposition (continued)
1. Slab decomposition of 3D FFT has three steps•

2D FFT (or two 1D FFT) along the two local dimension•Global transpose (communication)•

1D FFT along third dimension

2. Pros•

fast using small number of processes

3. Cons•

Limitation of the number of processes

Parallel 3D FFT –2D decomposition
MPI_Alltoall
communications among 3
processes of the same color

MPI_Alltoall
communications among 3
processes of the same color

Each process is assigned a data of size
grids onprocesses.

for computing a FFT of

Parallel 3D FFT –2D decomposition (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

2. The global transpose requires communication only between subgroups of all nodes
3. Cons: Slower than 1D decomposition for a small number of processes
4. Pros : Maximum parallelization is increased

Pseudo code of reciprocal space calculation with
2D decomposition of 3D FFT! 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 fft(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 fft(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 fft(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 fft(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 fft(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 fft(work_local)! FFT (F(X))
Q(my_rank) = work_local
end doend do!compute force

Parallel 3D FFT – 3D decomposition•

More communications thanexisting FFT•

MPI_Alltoall communications
only in one dimensional space•

Reduce communicational costfor large number of processesRef) J. Jung et al. Comput. Phys. Comm.
200, 57-65 (2016)

Case Study:
Parallelization in GENESIS MD softwareDomain decomposition in GENESIS SPDYN𝒄

1.

The simulation space is divided into
subdomain according to the numberof MPI processes.

2.

Each subdomain is further divided
into the unit domain (named cell).

3.

Particle data are grouped together
in each cell.

4.

Communications are considered
only between neighboringprocesses.

5.

In each subdomain, OpenMPparallelization is used.

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)Midpoint cell method•

For every cell pair, midpoint cell isdecided.

•If midpoint cell is not uniquelydecided, we decide it by consideringload balance.

•For every particle pairs, interaction
subdomain is decided from themidpoint cell.

•With this scheme, we can only
communicate data of each
subdomain only considering theadjacent cells of each subdomain.

Subdomain structure with midpoint cell method
1. Each MPI has local subdomain which is surrounded by boundary cells.
2. At least, the local subdomain has at least two cells in each direction

Local sub-domain
Boundary

Communication pattern (example of coordinates)
1. Each subdomain can evaluate interactions by obtaining the data of boundary region fromother processes
2. Using the communication pattern shown below, we can minimize the frequency ofcommunication.

① Receive from the upper process in x

③dimension② Receive from the lower process in x
dimension②①③ Receive from the upper process in ydimension④④ Receive from the lower process in ydimension

Communication pattern (example of forces)
1. Communication pattern of forces is opposite to the communication pattern of coordinates.
2. Each process send the force data in boundary cells.

① Send to the lower process in ydimension②② Send to the upper process in ydimension③④③ Send to the lower process in x
dimension①④ Send to the upper process in x
dimension

OpenMP parallelization in GENESIS
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

1

2

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 of FFT in GENESISFFT in GENESIS (2 dimensional view)

GENESIS(Identical domain
decomposition
between two space)

NAMD, Gromacs(Different domaindecomposition)

Parallelization based on CPU/GPU calculation

Computationally expensive part: GPU
Communicationally expensive part: CPU

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

z1x2

y2z2xNyN zNIn GENESIS, it is expressed as coord_pbc(1:3,1:natom,1:ncell)• Structure of Array (SoA)x1 x2x3x4zN-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)GENESIS developments on KNL
Simple pair-list scheme with better performance•

Pair-list scheme in GENESIS 1.X does not show good performance due to noncontiguous memory access (Figure (a))

GENESIS developments on KNL
Memoryusage

Number ofinteractions

CPU time(pair-list)

CPU time(non-bonded)

CPU time(total)

Algorithm 1 (AoS)

305.88 MiB47364770

9.95

36.42

46.37

Algorithm 1 (SoA)

305.88 MiB

47364770

6.60

33.71

40.31

Algorithm 2 (AoS)

48.4 MiB

122404436

4.91

102.46107.37

Algorithm 2 (SoA)

48.4 MiB

122404436

4.15

27.86

32.01

• Algorithm 1: pair-list algorithm used in GENESIS 1.X• Algorithm 2: New pair-list scheme• Target system : ApoA1 (about 90,000 atoms, 2000 time steps)Benchmark performance of GENESISGENESIS 1.0 performance on K (1)atom size:
macro molecule size:
metabolites size:
ion size:water size:
PBC size:

11737298
216 (43 type)4212 (76 type)23049 (Na+, Cl-, Mg2+, K+)
2944143
480 x 480 x 480(Å3)

Simulation time (ns/day)

0.06

Time per step (sec)

••
•
•
•
•

0.04

0.02

0.01
0.008
0.006
8192

16384

32768

65536

131072 262144

Number of cores

3224
16
128
64
3
2

8192

16384

32768

65536

131072 262144

Number of cores

GENESIS 1.0 performance on K (2)100nm = 1/10000mm

Simulation time (ns/day)Time per step (sec)

0.3
0.2

0.1
0.09
0.08
0.070.06
0.05
0.04
0.03
0.02
16384

32768

65536

131072

Number of cores

26214410987654321
0.9
0.8
0.7
0.6
0.5

16384

32768

65536

131072

Number of cores

262144

GENESIS performance on KNL (Oak-foreast)GENESIS performance on KNL (LANL Trinity)Efficient parallelization enable us to perform theMD simulations of the world largest system

Summary• Parallelization: MPI and OpenMP– MPI: Distributed memory parallelization– OpenMP: Shared memory parallelization• Parallelization of MD: mainly by domain decomposition with FFTparallelization• Key issues in parallelization: minimizing communication cost to maximize
the parallel efficiency