GENESIS tutorial 2:
Advanced MD (how to speed up GENESIS)
Workshop “Frontiers in ComputationalBiophysics and Biochemistry”
2017/02/28
J. JungContents• 15:30 – 17:00 GENESIS tutorial 2
 General comments in GENESIS How to accelerate MD speed of GENESIS : GPU
 How to accelerate MD speed of GENESIS :Multiple time step integration
 How to accelerate MD speed of GENESIS :Further suggestions

Test system• 92,224 atoms• 92,118 bonds• 74,136 angles• 74,130 dihedral angles

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

dri pidtmdp i Fidtpiri (t  t )  ri (t )  tmpi (t  t )  pi (t )  Fi tEquation 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 using PME
Etotal  kb (b  b0 )2O(N)

bonds



Total number of particles

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

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

126 
22

r
r
qqerfc(r)exp(/4)ki
jij0ij  0ij
2FFT(Q (k ))   ij  r 2
rij rijk k 0i  j  R   ij Real space, O(N)

Reciprocal space, O(NlogN)5

Domain decomposition in GENESIS• In GENESIS, the total system
is subdivided by subdomain
according to MPI processors.
• Subdomain is again
subdivided by cells (basicunit) which has dimensions
larger than half of cutoff
distance.

Comments on system size (1)• The subdomain size should be greater than cutoff distance• Cell size should be greater than half of cutoff distance• In each subdomain, there should be at least two cell in eachdirection (the number of cells in each dimension should bemore than twice of the number of domains).
• You can check this from output file
Setup_Boundary_Cell> Set Variables for Boundary Condition
domains (x,y,z) =222ncells (x,y,z) =888

Comments on system size (1)• The subdomain size should be greater than cutoff distance• Cell size should be greater than half of cutoff distance• In each subdomain, there should be atleast two cell in each direction(the number of cells in each dimension
should be more than twice of thenumber of domains for fast communication).

①rc / 2②②①

Set output
Variables•Setup_Boundary_Cell>
You can check this fromfile for Boundary Condition
domains (x,y,z) =ncells (x,y,z) =416

216210Comments on system size (2)• There should be at least 5 cells in each dimension. If not,program finishes with error message.
Input[ENERGY]
switchdist
cutoffdist
pairlistdist= 9.0= 10.0= 12.0[BOUNDARY]
type = PBC
box_size_x
box_size_ybox_size_z= 37.7= 37.7= 37.7OututSetup_Boundary_Cell> too small boxsize/pairlistdist. shorter
pairlistdist or larger boxsize or less MPI processors should
be used. rank_no =0

You can estimate cell sizes and numbers
 Bx

Bxnx  mod , Px rc 2
 rc 2
nx : number of cells in x dimension
Bx : System size in x dimension
rc : pairlist cutoff distancePx : Subdomain number in x dimension[ENERGY]
switchdist
cutoffdist
pairlistdist[BOUNDARY]
type = PBC
box_size_x
box_size_ybox_size_z

mpirun –np 1 {PATH}/spdyn test1.inp= 10.0= 12.0= 13.5domains (x,y,z) =ncells (x,y,z) =

416

2
16

2
11

mpirun –np 16 {PATH}/spdyn test1.inp= 108.8612= 108.8612= 77.758

domains (x,y,z) =ncells (x,y,z) =

4
16

2
16

2
10

If there are constraints for bonds involvinghydrogens, we should consider larger cell sizes
nx 

 B

Bx
 mod  x , Px rc 2
 rc 2nx No constraints•••Bx
Bx, Px  mod (rc  2) 2 (rc  2) 2Constraints

If there are constraints, only heavy atoms are considered for assigning
subdomain and cells
Hydrogen atoms have the same subdomain and cell assignment as the heavyatoms bonded.
In this case, we should consider pairwise distance up to parilsit cutoff + 2rOH1

HOH2

assigned to cellindex 1

Comments on MPI processors• MPI processors is automatically defined if domain information is
not written in the control input.
• If domain information is written in the control input, thenumber of domains should be identical to the number of MPIprocessors.
Input : test3.inp[BOUNDARY]
type = PBC
domain_x
domain_y
domain_z= 2= 2= 2

execution
mpirun –np 16 {PATH}/spdyn test3.inp

outputSetup_Processor_Number> # of process is not domain_x *domain_y * domain_z rank_no =0

Comments on PME grid numbers (1)• GENESIS is currently using FFTE for 1‐dim’l FFT kernel• In FFTE, only the multiple of 2, 3, and 5 are available as PME
grid numbers.
• If you write grid numbers that are not multiple of 2, 3, or 5,
the program automatically increases the grid numbers.
Input : test4.inp[ENERGY]
pme_ngrid_x
pme_ngrid_y
pme_ngrid_z

=
==

99
99
79

execution
mpirun –np 1 {PATH}/spdyn
test4.inp

outputWARNING: PME grid number is different from the input
pme_ngrid(x,y,z)=
100
10080

Comments on deciding PME grid numbers (2)Furthermore, grid numbers should be divisible by some combinationof domain numbers– pme_grid_x : multiple of (2×domain_x)– pme_grid_y : multiple of (domain_y×domain_z)– pme_grid_z : multiple of (domain_x×domain_z) and multiple of(domain_y×domain_z)
Input : test4.inp[ENERGY]
pme_ngrid_x
pme_ngrid_y
pme_ngrid_z

execution
= 99
= 99
= 79

mpirun –np 16 {PATY}/spdyn
test4.inp

outputWARNING: PME grid number is different from the input
pme_ngrid(x,y,z)=
10012080

Why there are restrictions of PME grid numbers•

GENESIS make use of volumetricdecomposition for FFT which requires
one dimensional MPI_Alltoall for threedirections•To perform MPI_Alltoall, some
restrictions exist (divisibility)•

Without restrictions, we can perform
MPI_Alltoallv, but it decreases theperformanceMy suggestions• Write domain information explicitly in the control input andrun with corresponding MPI processors• Given domain information, write PME grid numbers that is
not changed in execution (grid spacing should be less than1.5).
Input : test5.inp[ENERGY]
pme_ngrid_x= 128
pme_ngrid_y= 128
pme_ngrid_z= 96[BOUNDARY]
type
domain_x
domain_y
domain_z= PBC= 4= 2= 2execution
mpirun –np 16 {PATH}/spdyn
test4.inpExercise
1.

Run test1.inp and test2.inp with 16 and 1 MPI processors. Check the domain
size and cell sizes from the outputs.

2.

Run test3.inp with 16 and 8 MPI processors.

3.

Run test4.inp with various MPI processors and check domain numbers andPME grid numbers.

4.

Run test5.inp with 16 MPI processors. Write inputs for 32, 64, and 128 MPI
processors in which PME grid numbers do not change when executed.

Please run mpirun with orte_base_help_aggregate 0(mpirun –mca orte_base_help_aggregate 0 ./spdyn )Contents• 15:30 – 17:00 GENESIS tutorial 2
 General comments in GENESIS How to accelerate MD speed of GENESIS : GPU
 How to accelerate MD speed of GENESIS :Multiple time step integration
 How to accelerate MD speed of GENESIS :Further suggestions

Overview of CPU+GPU calculations
1.

Computation intensive work : GPU• Pairlist• Real space non‐bonded interaction

2.

Communication intensive work or nocomputation intensive work : CPU• Reciprocal space non‐bondedinteraction with FFT• Bonded interaction• Exclusion list

3.

Integration is performed on CPU dueto file I/O.

Compile of GENESIS for GPU usage• CPU+GPU with single precision :– ./configure –enable-single –enable-gpu
with-cuda={path}• CPU+GPU with double precision :– ./Configure –enable-gpu –with-cuda={path}• Only CPU with single precision :– ./configure –enable-single• Only CPU with double precision :– ./configure• Here, we prepared four binaries
CPU single : /home2/data/genesis/bin.CPU.sp/spdyn
CPU double : /home2/data/genesis/bin.CPU.dp/spdyn
CPU+GPU single : /home2/data/genesis/bin.GPU.sp/spdyn
CPU+GPU double : /home2/data/genesis/bin.GPU.dp/spdyn–Comparison of speed(Time for 2000 MD steps)••••

CPU : Intel(R) Xeon(R) CPU E5‐2680 v3 @ 2.50GHz, 24 cores
GPU : NVIDIA GeForce GTX 1080
8 MPIs with 3 OpenMPs
Test input file : test6.inp
Condition

Time (sec)

CPU+GPU (double)

40.685

CPU+GPU (single)

26.307

CPU only (double)

109.52CPU only (single)

101.284

Contents• 15:30 – 17:00 GENESIS tutorial 2
 General comments in GENESIS How to accelerate MD speed of GENESIS : GPU How to accelerate MD speed of GENESIS :Multiple time step integration How to accelerate MD speed of GENESIS :Further suggestions

Multiple time step integration (1)• Real space bonded interaction : small amount of computation• Real space nonbonded interaction : large amount ofcomputation, well scalable by increasing processors (mainbottleneck when using small number of processors)• Reciprocal space nonbonded interaction : medium amount ofcomputation, not well scalable by increasing processors (mainbottleneck when using large number of processors)• Purpose of multiple time step integration : Not to calculate thecommunication‐intensive reciprocal space nonbonded
interactions every step to increase the performance

Multiple time step integration (2)Conventional
Update velocity
t F (t ) tv i  t  short   v i  t    i
2 
2mi

Update coordinate t ri  t  t   ri  t   t  v i  t  2 Force calculation Fi (t  t )MTS
If slow force,
Update velocity

vi  t   vi  t  tlong Fi ,long (t )mi2tshort Fi ,short (t ) tshort vtvt
Update velocity i i2 2mi

Update coordinate tri  t  tshort   ri  t   tshort  v i  t  short 2 Force calculation Fi ,short (t  t )If slow force, calculate Fi ,long (t  t )Update velocity

Update velocity t  t F (t  t )v i  t  t   v i  t     i
2  2miF(t  tshort ) t tvi  t  tshort   vi  t  short   short  i ,short2 2mi

If slow force,Ft(t )
vi  t   vi  t   long  i ,longmiUpdate velocity
2

Multiple time step integration with GPU
1. In the multiple time step
integrator, we do not perform
reciprocal space interaction
every step.
2. If reciprocal space interaction
is not necessary, we assign
subset of real spaceinteraction on CPU tomaximize the performance.
3. Integration is performed on
CPU only.

Performance result with Multiple time step
integration
Input : test7.inp

Input : test6.inp[DYNAMICS]
integrator = VVER

[DYNAMICS]
integrator= VRES
nsteps= 2000elec_long_period = 2
thermostat_period = 2
barostat_period= 2

Timing for 2000 MD steps
Condition

Test6.inp

Test7.inp

CPU+GPU (double)

40.68539.153

CPU+GPU (single)

26.30722.749CPU only (double)

109.5299.520CPU only (single)

101.2894.691

Speed up due to multiple time step becomes more critical as we increase the number ofprocessors.

Contents• 15:30 – 17:00 GENESIS tutorial 2
 General comments in GENESIS How to accelerate MD speed of GENESIS : GPU
 How to accelerate MD speed of GENESIS :Multiple time step integration
 How to accelerate MD speed of GENESIS :Further suggestions

Energy output•••

If energy output is required, the program calculates both energy and force
If energy output is not required, the program calculates only force.
Please do not print energy output frequently.
Input : test6.inp[DYNAMICS]
eneout_period = 100

Input : test8.inp[DYNAMICS]
eneout_period = 1

Timing for 2000 MD steps
Condition

Test6.inp

Test8.inp

CPU+GPU (double)

40.68548.032

CPU+GPU (single)

26.30729.897

CPU only (double)

109.52137.89CPU only (single)

101.28120.73

Decision of MPI and OpenMP processors (1)• MPI– Distributed memory parallelization– Communication is necessary• OpenMP– Shared memory parallelization– Communication is not necessary• MPI/OpenMP hybrid parallelization in GENESIS– Domain decomposition by MPI– In each domain, force decomposition by OpenMP

Decision of MPI and OpenMP processors (2)• Our suggestions– For small number of processors, flat MPI is preferred– For multiple nodes, OpenMP threads could be helpful– More than 16 OpenMP threads is not recommended.

Timing for 2000 MD steps
Condition

MPI=16, OMP=1

MPI=4, OMP=4CPU+GPU (double)

43.76550.496

CPU+GPU (single)

30.44736.406

CPU only (double)

151.34158.89CPU only (single)

145.84146.85

Exercise1.

Compare the performance by running test6.inp (conventional integration) andtest7.inp (multiple time step integration)

2.

Compare the performance by running test6.inp (energy output every 100steps) and test8.inp (energy output every step)

3.

Compare the performance by running test9.inp with various combinations ofMPIs and OpenMPs