第13回 配信講義　計算科学技術特論A（2025）

計算科学技術特論Ａ Parallelization of molecular dynamics 2025年 7⽉ 10⽇ 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 of motion 1 1 𝐩! 𝑡 + ∆𝑡 = 𝐩! 𝑡 − ∆𝑡 + 𝐅! 𝑡 ∆𝑡 2 2 1 +! 𝑡 + ∆𝑡 ∆𝑡 𝐩 2 𝐫! 𝑡 + ∆𝑡 = 𝐫! 𝑡 + 𝑚! Integration MD trajectory => Ensemble generation Long time MD trajectories are important to obtain thermodynamic quantities of target systems.

MD integration (Velocity Verlet case) 1 𝑡 + ∆𝑡 2 𝑡 5 𝑡 + ∆𝑡 2 𝑡 + 2∆𝑡 coordinate update coordinate update coordinate update 𝐪 𝑡 → 𝐪 𝑡 + ∆𝑡 𝐪 𝑡 + ∆𝑡 → 𝐪 𝑡 + 2∆𝑡 𝐪 𝑡 + 2∆𝑡 → 𝐪 𝑡 + 3∆𝑡 𝐩𝑡 1 → 𝐩 𝑡 + ∆𝑡 2 𝐅(𝑡) 3 𝑡 + ∆𝑡 2 𝑡 + ∆𝑡 1 𝐩 𝑡 + ∆𝑡 2 → 𝐩 𝑡 + ∆𝑡 𝐩 𝑡 + ∆𝑡 3 → 𝐩 𝑡 + ∆𝑡 2 𝐅(𝑡 + ∆𝑡) 3 𝐩 𝑡 + ∆𝑡 2 → 𝐩 𝑡 + 2∆𝑡 𝐩 𝑡 + 2∆𝑡 5 → 𝐩 𝑡 + ∆𝑡 2 𝐅(𝑡 + 2∆𝑡) 𝑡 +3∆𝑡 5 𝐩 𝑡 + ∆𝑡 2 → 𝐩 𝑡 + 3∆𝑡 𝐅(𝑡 + 3∆𝑡)

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 of milliseconds or longer. 3. We need to accelerate energy/force calculation for long time MD fs ps ns μs ms sec vibrations Sidechain motions Mainchain motions Folding Protein global motions

History of Biological MD simulations (1977 to current) ps µs ns ms relaxation vibration folding Secondary structure change protein Reaction of protein 1977 BPTI (~6x102) Membrane protein HP36 (~104) AQP1 (~105) BPTI, Ubiquitin (~2x104) Ribosome (~2x106) 2013 ANTON (Gordon Bell Prize 2009) 2000 HIV-1 capsid (~6.5x107) Macromolecule complex Mycoplasma genitalium (~1.0x108) Cellular environment ~2019 KNL size (# of atoms) time s Chromatin (~1.0x109) Sars-Cov-2 (~3.0x108) HIV-1 capsid (~6.5x107) Chromatophore vesicle (~1.4x108) We need to extend simulation time & size efficient parallelization 6

Potential energy in MD (all-atom case)
𝐸"#$%&$'()
= ( 𝑘- 𝑏 − 𝑏. /

bond

O(N) Total number

angle

O(N)

of particles

*#&+,

+ ( 𝑘1 𝜃 − 𝜃. /
(&0)%,

+

(

𝑘5 1 + 𝑐𝑜𝑠 𝑛𝜑 − 𝛿

dihedral angle

O(N)

𝑘7 𝜔 − 𝜔. /

improper dihedral angle

O(N)

+'2%+3()4

+

(
'6"3#"%3

+

(
&#&8*#&+%+
"('3,

:!;
𝜀!9

6'&
𝑅!9
𝑟!9

</

6'&
𝑅!9
−
𝑟!9

Main bottleneck in MD

=

+

𝑞! 𝑞9
𝑟!9

non-bonded
(van der Waals and electrostatic)O(N2)
7

Practical evaluation of non-bonded interaction
1. We evaluate non-bonded interaction with periodic boundary condition with infinite
images.
:!;
( 𝜀!9
!9

6'&
𝑅!9

𝑟!9

</

−

6'&
𝑅!9

𝑟!9

=

𝑞! 𝑞9
:!;
+
→ ( 𝜀!9
𝑟!9
!9,𝐧

6'&
𝑅!9

𝑟!9,𝐧

</

−

6'&
𝑅!9

𝑟!9,𝐧

=

+

𝑞! 𝑞9
𝑟!9,𝐧

Practical evaluation of non-bonded interaction
2. Truncation of van der Waals interaction with 𝑉 𝑟 = 0 with 𝑟 > 𝑅+
:!;
( 𝜀!9
!9,𝐧

6'&
𝑅!9

𝑟!9,𝐧

</

−

6'&
𝑅!9

𝑟!9,𝐧

=

+

𝑞! 𝑞9
→ (
𝑟!9,𝐧

!9,𝐧
@!" AB#

i-th
particle

:!;
𝜀!9

6'&
𝑅!9

𝑟!9,𝐧

</

−

6'&
𝑅!9

𝑟!9,𝐧

=

+(
!9

𝑞! 𝑞9
𝑟!9,𝐧

Range of j-th
particle that
interacts with i-th
particle

Practical evaluation of non-bonded interaction
3. Truncation of real-space electrostatic interaction by decomposing the electrostatic
interaction into real- and reciprocal-space

(

:!;
𝜀!9

!9,𝐧
@!" AB#

→ (
!9,𝐧
@!" AB#

:!;
𝜀!9

6'&
𝑅!9

</

𝑟!9,𝐧
6'&
𝑅!9

𝑟!9,𝐧

−

6'&
𝑅!9

=

𝑟!9,𝐧

</

−

6'&
𝑅!9

𝑟!9,𝐧

real-space

+(
!9
=

𝑞! 𝑞9
𝑟!9,𝐧
G

𝑞! 𝑞9 erfc(𝛼𝑟!9,𝐧 )
1
exp( − 𝜋 / 𝐤 / /𝛼 / )
𝛼
/
+
+
(
𝑆(𝐤)
−
( 𝑞!/
/
𝑟!9,𝐧
2𝜋𝑉
𝐤
𝜋
𝐤D𝟎

reciprocal-space

!F<

self term

Parallelization is one of the solutions to accelerate MD simulations Serial Parallel X16? 1cpu CPU Core C 1 CPU MPI Comm 16 CPUs Good Parallelization ; 1) Small amount of computation in one process 2) Small communicational cost

Parallelization

Shared memory parallelization (OpenMP) C Memory #pragma omp parallel Fortran !$omp parallel P1 P2 P3 PP-1 PP • 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) M M M 1 2 3 MP-1 M P … P1 P2 P3 call mpi_init call mpi_comm_rank PP PP call_mpi_comm_size -1 • 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) call mpi_init M1 M2 M3 MP-1 MP call_mpi_comm_size … P1 P2 P3 call mpi_comm_rank . . . PP-1 PP !$omp parallel do do i = 1, N .. • Combination of shared memory and distributed memory parallelization. • It is useful for minimizing communicational cost with very large number of processors end do !$omp end parallel do

Low level parallelization: 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 the performance

Low level parallelization in GPU: SIMT (Single instruction, multiple threads) • Many threads execute the same instruction simultaneously. • Each thread can work on different data. Figure is from NVIDIA’s Volta architecture whitepaper

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. do i = 1, N do j = i+1, N energy(i,j) force(i,j) end do end do 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 do end do MPI reduction (energy,force)

Parallelization scheme 1: Replicated data approach atom indices Computational cost 1 2 3 4 5 6 7 8 MPI rank 0 MPI rank 1 MPI rank 2 MPI rank 3 1. Perfect load balance is not guaranteed in this parallelization scheme 2. Reduction with communication prevents the larger number of MPIs.

Hybrid (MPI+OpenMP) parallelization of the Replicated data approach 1. Works are distributed over MPI and OpenMP threads. 2. Parallelization is increased by reducing the number of MPIs involved in communications. do i = my_rank+1, N, proc do j = i+1, N energy(i,j) force(i,j) end do end do MPI reduction !$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 do end do Openmp reduciton !$omp end parallel do i = my_rank+1,N,proc !omp parallel do do j = i+1, N or energy(i,j) force(i,j) end do !$omp end parallel do end do MPI reduction MPI reduction my_rank: MPI rank, proc: total number of MPIs, nthread: total number of OMP threads

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 of processes • 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 into subdomains according to MPI (different colors for different MPIs). 2. Each MPI only considers the corresponding subdomain. 3. MPI communications only among neighboring processes. Communications between processors

Parallelization scheme 2: Domain decomposition 1. For the interaction between different subdomains, it is necessary to have the data of the buffer region of each subdomain. Buffer Subdomain 𝒓𝒄 2. The size of the buffer region is dependent on the cutoff values. 3. The interaction between particles in different subdomains should be considered very carefully. 4. The size of the subdomain and buffer region is decreased by increasing 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

• https://en.wikipedia.org/wiki/Hilbert_curve

Comparison of two parallelization scheme Computation Communicatio n Memory Replicated data O(N/P) O(N) O(N) Domain decomposition O(N/P) O((N/P)2/3) O(N/P) N: system size P: number of processes

MPI Parallelization of MD (reciprocal space)

Electrostatic interaction with particle mesh Ewald method G 𝑞! 𝑞9 erfc(𝛼𝑟!9,𝐧 ) 1 exp( − 𝜋 / 𝐤 / /𝛼 / ) 𝛼 / ( + ( 𝑆(𝐤) − ( 𝑞!/ / 𝑟!9,𝐧 2𝜋𝑉 𝐤 𝜋 !9,𝐧 @!" AB# 𝐤D𝟎 real-space !F< reciprocal-space self energy The structure factor in the reciprocal part is approximated as 𝑆(𝑘< , 𝑘/ , 𝑘H ) = 𝑏< (𝑘< )𝑏/ (𝑘/ )𝑏H (𝑘H )𝐹(𝑄)(𝑘< , 𝑘/ , 𝑘H ) Fourier Transform of charge It is important to parallelize the Fast Fourier transform efficiently in PME!!

Overall procedure of reciprocal space calculation 𝑸 (charge on real space) " 𝑸(charge on grid) FFT Energy calculation " 𝚯 " 𝑭𝑭𝑻(𝑸)× " 𝑭𝑭𝑻(𝑸) Inverse FFT " 𝚯) " 𝑰𝑭𝑭𝑻(𝑭𝑭𝑻(𝑸)× Force calculation " 𝚯) " from 𝑰𝑭𝑭𝑻(𝑸×

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 D2 B1 B2 B3 B4 C1 C2 C3 C4 D1 D2 D3 D4 A3 B3 C3 D3 A4 B4 C4 D4 MPI_alltoall MPI_alltoall is same as matrix transpose!!

Parallel 3D FFT – slab (1D) decomposition 1. Each process is assigned a slab of size 𝑁×𝑁× 𝑁⁄𝑃 for computing FFT of 𝑁×𝑁×𝑁 grids on 𝑃 processes. 2. The scalability is limited by 𝑁 3. 𝑁 should be divisible by 𝑃

Parallel 3D FFT – slab (1D) decomposition No communication MPI_Alltoall communication among all preocessors

Parallel 3D FFT – slab (1D) decomposition 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

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

Parallel 3D FFT – 3D decomposition • More communications than existing FFT • MPI_Alltoall communications only in one dimensional space • Reduce communicational cost for large number of processes Ref) J. Jung et al. Comput. Phys. Comm. 200, 57-65 (2016)

Case Study (1) : Parallelization in GENESIS MD software on the Fugaku supercomputer

MD software GENESIS • GENESIS has been developed in RIKEN (main director: Dr. Yuji Sugita). • It allows high-performance MD simulations on parallel supercomputers like Fugaku, LUMI, etc. • It is free software under LGPL license. GENESISの開発の歴史 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 Start of the development at RIKEN Wako Kobe teams started Ver. 1.0 Ver. 1.4 Ver. 1.2 Ver. 2.0 Optimized for Fugaku Ver. 1.1 Ver. 1.3 Ver. 1.5 Latest release https://mdgenesis.org Ver. 2.1 年度

• https://mdgenesis.org

Domain decomposition in GENESIS SPDYN 𝒓𝒄 /𝟐 1. The simulation space is divided into subdomain according to the number of 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 neighboring processes. 5. In each subdomain, OpenMP parallelization is used.

Cell-wise particle data in GENESIS 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 or pairs of cells. particle data in traditional cell lists cell-wise arrays of particle Ref : P. Gonnet, JCC 33, 76 (2012) data J. Jung et al. JCC, 35, 1064 (2014)

Midpoint cell method • For every cell pair, midpoint cell is decided. • If midpoint cell is not uniquely decided, we decide it by considering load balance. • For every particle pairs, interaction subdomain is decided from the midpoint cell. • With this scheme, we can only communicate data of each subdomain only considering the adjacent cells of each subdomain. Ref : J. Jung et al. JCC, 35, 1064 (2014)

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 𝑟# /2 Local sub-domain Boundary

Hybrid Parallelization of GENESIS (real-space) 5 6 7 9 11 1 2 10 3 4 12 8 13 14 15 16 Unit cell Subdomain of each MPI OpenMP parallelization thread1 Cell pair Midpoint cell? (1,1) Yes (2,2) Yes … … (1,2) Yes thread4 (1,3) Yes thread1 … … (1,5) Yes (1,6) No (1,7) No (1,8) Yes (1,9) No (1,10) Yes … … thread2 thread2 thread3 thread4 The parallelization scheme is changed on Fugaku for higher speed!!

Overall procedure of reciprocal-space interaction 𝑸 (charge in real space) PME_Pre " 𝑸(charge grid data) FFT ") E(𝑸 " 𝚯 " 𝑸× ") 𝑭𝑭𝑻(𝑸 Inverse FFT " 𝚯) " 𝑰𝑭𝑭𝑻(𝑸× PME_Post Force from " 𝚯) " 𝑰𝑭𝑭𝑻(𝑸×

Charge grid data evaluation in GENESIS 1.X grids 1.Charge values is transferred to the neighbor grids. 2.The amount of neighbor depends on the spline order. Real-space 𝐷Q 𝐷R 𝐷S 𝐷T 𝐷U 𝐷V 𝐷W 𝐷X Generated charge grid Accepted charge grid Reciprocal-space "Q 𝐷 "R 𝐷 "S 𝐷 "T 𝐷 "U 𝐷 "V 𝐷 "W 𝐷 "X 𝐷 𝐷Y 𝐷QZ 𝐷QQ 𝐷QR "Y 𝐷 "QZ 𝐷 "QQ 𝐷 "QR 𝐷 𝐷QS 𝐷QT 𝐷QU 𝐷QV "QS 𝐷 "QT 𝐷 "QU 𝐷 "QV 𝐷 In this scheme, we can avoid communication in charge grid data generation.

Problem of charge grid data in GENESIS 1.X Reciprocal-space Real-space 𝐷Q 𝐷R 𝐷S 𝐷U 𝐷V 𝐷T 𝐷W 𝐷X ? "Q 𝐷 "R 𝐷 "S 𝐷 "T 𝐷 "U 𝐷 "V 𝐷 "W 𝐷 "X 𝐷 𝐷Y 𝐷QZ 𝐷QQ 𝐷QR "Y 𝐷 "QZ 𝐷 "QQ 𝐷 "QR 𝐷 𝐷QS 𝐷QT 𝐷QU 𝐷QV "QS 𝐷 "QT 𝐷 "QU 𝐷 "QV 𝐷 • We assume that charge grids in reciprocal-space is obtained from the real-space subdomain and its adjacent cells. • If the spline order and grid spacing is large (e.g. spline order 8 with grid spacing 2Å), the charge data of adjacent cell is not sufficient to obtain the charge grid data in the subdomain.

New scheme of charge grid data for Fugaku (1) Real-space Reciprocal-space • Charge grid data is obtained from the charge values in the subdomain only. • In this way, we don’t have to consider charge grid data loss even for large spline order.

New scheme of charge grid data for Fugaku (2) Real-space Reciprocal-space 1. In the new scheme, we generate charge grid data only from the subdomain. 2. Charge grid data generated out of the subdomain is sent to next subdomain with communication. 3. It accelerates the speed by reducing operation and higher order spline can be used. 4. The communication used is not global, so it does not lose the performance significantly.

Conventional MD (strong scaling) 4 times better performance than GENESIS1.0 on K using 1/8 nodes 11.9 ns/day Ref : J. Jung et al. JCC, 42, 231 (2021)

Case Study (2) : Parallelization of coarse-grained MD simulations

Problem of large scale CG 180x180x3000 Å3, ~20,000 particles More than 100 times larger number of particles => It will take very long time for simulation. < 20,000 particles We can consider large scale CG MD with parallelization technique.

Problem of large-scale CG simulation with parallelization IDP proteins with CG Low density Macromolecular crowding (AAMD) High density Particle density is not uniform in CG MD simulations => difficult to parallelize Particle density is almost uniform due to water molecules

New developments for large scale CG MD simulations • Different subdomain size with similar number of particles. • With this scheme, we can avoid load imbalance problem in CG. • We implement this scheme into CGDYN engine in GENESIS software Ref: Jung et al. Nat. Commun. 15, 3370 (2024)

Dynamic load balancing in CGDYN In CGDYN, we can apply load balancing during MD simulations [DYNAMICS] . . . lbupdate_period = 1000000 Jung et al. Nat. Commun. 15, 3370 (2024)

CGDYN improves the performance significantly SPDYN-like CGDYN Jung et al. Nat. Commun. 15, 3370 (2024)

Effect of dynamic domain decomposition (on Fugaku) • Total simulation time is 107 steps. • We used 256 nodes for the system shown as the right figure. • Tlb is the dynamic domain decomposition period (in the input, it is written as lbupdate_period). Jung et al. Nat. Commun. 15, 3370 (2024)

CG MD simulation of multiple droplet fusions (1) Jung et al. Nat. Commun. 15, 3370 (2024)

CG MD simulation of multiple droplet fusions (2) Jung et al. Nat. Commun. 15, 3370 (2024)

Case Study (3) : GPU Parallelization by GROMACS

Cluster for efficient GPU calculation • Particle-particle interaction => Group-group interaction • The group size depends on the warp size of GPU. Páll et al,. J. Chem. Phys. 153, 134110 (2020)

CPU-GPU heterogeneous offloading Páll et al,. J. Chem. Phys. 153, 134110 (2020)

Case Study (4) : GPU Parallelization by GENESIS

GENESIS GPU kernel calculation Jung et al. JCTC, 12, 4947-4958 (2016)

GENESIS GPU kernel calculation (SIMT) • Assignment of GPU threads according to the GPU hardware architecture is important. Diego et al. HPCAsia 25, https://doi.org/10.1145/3712031.3712036

GENESIS CPU-GPU offloading • GPU: Computationally intensive part • CPU: Communicationally intensive part Jung et al. JCTC, 12, 4947-4958 (2016)

Summary • Parallelization: MPI and OpenMP • MPI: Distributed memory parallelization • OpenMP: Shared memory parallelization • Hybrid: both MPI and OpenMMP • Parallelization of MD: mainly by domain decomposition with parallelization of non-bonded interaction • Key issues in parallelization: minimizing communication cost to maximize the parallel efficiency, hardware aware optimization.