マルチスケール分⼦動⼒学ソフトウェアGENESISの並列化
Supercomputing Japan 2024
2024/03/12 タワーホール船堀

BDR

Yuji Sugita
RIKEN Center for Computational Science
RIKEN Center for Biosystems Dynamics Research
RIKEN Cluster for Pioneering Research
RIKEN iTHEMS

Protein Dynamics and Functions• Protein structures are determined by X-ray crystallography, NMR , andcryo-EM at atomic resolutions.
• Proteins are flexible for their functions in a living cell .
Enzyme reaction by Adenylate Kinase (Adk)

ATP + AMP

2ADP

Conformational changes of Adk
Open-state

NMP

LID

CORE

Closed-state

LID

NMP

CORE
Y. Matsunaga, et al., PLoS Comput. Biol. 8, e1002555 (2012).

Molecular dynamics (MD) simulation• MD is useful to predict molecular motions at the atomic resolutions.
• In classical MD, newton‘s equation of motion is solved iteratively.

𝐹 = 𝑚𝑎• Time step △ 𝑡 should be very short to reproduce fast vibrations
in molecules. (10-15 sec). Many iterations are required to detect
protein dynamics.
1
𝐩! 𝑡 + ∆𝑡 ∆𝑡
2r 𝑡
𝐫! 𝑡 + ∆𝑡 = 𝐫! 𝑡 +𝑚!
1
1𝐩! 𝑡 + ∆𝑡 = 𝐩! 𝑡 − ∆𝑡 + 𝐅! 𝑡 ∆𝑡
1
22𝐩 𝑡 − ∆𝑡
2𝑡𝑡 + ∆𝑡r 𝑡 + ∆𝑡
1
𝐩 𝑡 + ∆𝑡2𝐅(𝑡)𝐅(𝑡 + ∆𝑡)

Molecular Force Field Function
V ( r ) = ∑ kb ( ri,i+1 − r0 )2

Bond term

bonds+ ∑ ka (θ i,i+1,i+2 − θ 0 )2Angle termtaken from PumMa HP

angles+ ∑ kt ⎡⎣1 + cos ( nφi,i+1,i+2,i+3 − ω )⎤⎦torsions⎡⎛ σ ⎞ 12 ⎛ σ ⎞ 6 ⎤+ ∑ ∑ ε ⎢⎜ ⎟ − ⎜ ⎟ ⎥rij ⎠ ⎥⎢⎝ rij ⎠⎝i j>i⎣⎦qi q j+∑ ∑i j>i 4πε 0 rij

Dihedral term

vdW term(solved with cutoffapproximation)Electrostatic term

O(N2)CHARMM, AMBER, OPLS, andother Force Fields are available
for various molecules.
4Jung, Mori et al., WIREs Comput. Mol. Science. (2015); Kobayashi, Jung et al., JCC. (2017)

Development of GENESIS
Initial Developers• We have developed GENESIS software since 2019.
• GENESIS can show good performances on K andFugaku.
• GENESIS is freely available for everyone under
LGPLv3 license.

2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029

Start of the development
at RIKEN Wako

Ver. 1.0Kobe teams started

https://www.r-ccs.riken.jp/labs/cbrt/

Ver. 1.2
Ver. 1.1

Ver. 1.3

Ver. 1.6

Ver. 1.4
Ver. 1.5

Ver. 1.7

Ver. 2.0
Ver 2.1

Latest Version
Ver. 2.0beta

YearGENESIS is freely available
https://www.r-ccs.riken.jp/labs/cbrt/

Click

Click

How can we accelerate MD
simulations to understand the slow
conformational changes of proteins(if we can use supercomputers) ?

RIKEN Center for Computational Science
The Supercomputer K
NodeCPU

SPARC64™ VIIIfx 2GHz

Performance

128 GF(16 GF x 8 cores)

Memory16 GB

Number of Nodes per Rack

864

Total Number of Nodes

82,944

NetworkTofu Interconnect(6D Mesh/Torus)

Peak Performance

10.62 PF

Total Memory

1.26 PB

File IO

Fujitsu Exabyte File System(FEFS)

Storage

30 PB

Top500No.1 for 2 consecutive times since Jun. 2011

Graph500 No.1 in June 2014 & for 6 consecutive times
since Jul. 2015
No.2 for 4 consecutive times since Nov. 2014
HPCG
No.1 for 3 consecutive times since Nov. 2016
The K computer strikes a balance between performance incalculation, memory, and communication.
Shared use of the K computer has ended on August 16, 2019.

Domain Decomposition Scheme

Divide a whole space into subdomains and cells →Greatly reduce the memory

J. Jung, T. Mori and Y. Sugita, J. Comp. Chem. 35, 1064 (2014)

Midpoint Cell Method

Midpoint Method:

Interaction is computed at the cell where themidpoint of two interacting particles is located.
(¢) Low communication cost(✕) Need to decide the midpoint every time step.
Kevin J. Bowers et al. , J. Chem. Phys. 124, 184109 (2006).

Midpoint Cell Method:

Interaction is computed at the midpoint cell foreach cell pair.
(¢) Low communication cost(¢) Midpoint cell is decided only once.

J. Jung, T. Mori and Y. Sugita, J. Comp. Chem. 35, 1064 (2014)

Subdomains and cells
(a)

012

2. A whole system is divided into subdomains
and cells.

bac

dn−1
1n−1

Subdomain assigned by MPICell

1. The basic domain decomposition scheme is
the midpoint cell method.

Boundary cell

3. Information of the particles residing at each
subdomain (plus surrounding cells) is stored ateach MPI process.
4. Each subdomain is further divided into cells
for OpenMP thread parallelization.
This method is suitable to Hybrid parallelization.

Intracellular Environments• Concentration of proteins and RNA
in the E. coli cell is 300-400mg/ml.
• 20-40% of the volume inside of thecell is occupied by macromolecules.
• This condition has been called as
macromolecular crowding, which
can affect structure, dynamics, andfunctions of macromolecules in thecell.
Zimmerman and Minton, Ann. Rev. Biophys. Biomol. Struct. (1993) 22: 27-65.

Illustration by Goodsell for E. coli cellZhou, H.X., Rivas, G., Minton, A.G., Ann. Rev. Biophys. (2008) 37, 375-397.

Go to see more arts in his homepage: https://ccsb.scripps.edu/goodsell/Feig, Harada, Mori, Yu, Takahashi, Sugita : J. Mol. Graph. Model 58, 1-9 (2015).

Atomistic Modeling of the Cytoplasm
Homology models of protein structures

AlphaFold and RosettaFold werenot available in 2015.

M. Feig(Michigan State University)

Disordered regions were notincluded in the models.

Metabolites in M. Genitalium

Fully hydrated models

・# of atoms: 103,708,785
・# of proteins: 1238
・# of RNAs: 284
・# of metabolites: 41,006
・# of ions: 214,000
・# of water: 26,263,505

Yu et al. eLife 5, e19274 (2016).

The cytoplasm in M. genitalium

The first all-atom
MD calculation of
the bacterialcytoplasm(> 100 M atoms)

Isseki Yu (RIKEN→ Maebashi Inst. Tech.)

M. Feig(Michigan State University)

Destabilization in Crowded Environments
PDHA destabilization in Crowded Environments

Volume Exclusion Effect in crowdedenvironments

Blue: available, Red: excluded (Minton, 2001)
Free space for large molecules is very limited.

This is inconsistent with the excluded volume effect.
Yu et al. eLife 5, e19274 (2016).

In the crowded condition, compact structures
are preferred.
Zimmerman and Minton, Ann. Rev. Biophys. Biomol. Struct. (1993) 22: 27-65.

ATP controls solubility
ATP distributions on the Protein Surfaces

Interactions between proteins and ATPs areimportant.
Yu et al. eLife 5, e19274 (2016).

A hidden function of ATP in the cells(hydrotrope in crowded conditions)Patel et al. Science 356, 753-756 (2017)

RIKEN Center for Computational Science
The Supercomputer Fugaku
NodeCPU

Armv8.2-A SVE 512bit

Performance

3.072TF（DP）
6.144TF（SP）
48 (+2) coreMemory32 GB

Number of Nodes per Rack384 node x 396 rack
192 node x 36 rack

Total Number of Nodes

158,976

NetworkTofu Interconnect D

Peak Performance

488 PF（DP）
977 PF（SP）

Total Memory

4.85 PB

File IOLLIOOSRed Hat Enterprise Linux 8

Graph500

No. 1 in June and Nov. 2020
No. 1 in June and Nov. 2020

HPL-AI

No. 1 in June and Nov. 2020

HPCG

No. 1 in June and Nov. 2020

Top500

Fugaku is now the fourth fastest supercomputer in the world.

Optimization of GENESIS on Fugaku
New non-bonded interaction algorithm in GENESIS

do ijcel = 1, cell_pair
obtain icel and jcel
do i=1,natom(icel)
do k=1,neighbor(i,ijcel)
j = list(k,ijcel)
interaction end do
end doend dodo icel = 1, cell
do i=1,natom(icel)
do k=1,neighbor(i,icel)
list(k,i,icel)
interaction end do
end doend doApoA1 on one MPI processor

Single/Mixed/Double precision•

Single precision• 32 bit• Fast evaluation is possible by SIMD.
• Less accurate than double precision•

Double precision• 64 bit• More accurate than single precision•We can define precision levels in GENESIS 2.0 according to our purposes• Single precision: All real numbers except energy/constraints are described by single precision• Mixed precision: Integration by double precision and force calculation by single precision• Double precision: All numbers are described by double precision

Picture is from wikipedia.

Multiple kernels for nonbonded energy calculations
User Friendly Multiple Kernels for Nonbonded Energy Calculations

GENESIS Multi-Kernel (Version 2.0)
Optimal code is generated, when users are compilingIntel CPU••••CPU + GPU

K computerFugaku

GENESIS 2.0 (the latest version) includes multiple kernels for nonbonded energy calculations.
Users need to select a proper option when they compile GENESIS on each computational platform.
Then, optimal code is generated automatically.
Precision is also needed to select when compile GENESIS on each platform. We recommend
single_precision for intel compiler and mixed_precision for Fugaku, considering the performance as well asaccuracy.

Hybrid GPU+CPU application in GENESIS
start
Control &Enefunc
Integration (VV1)

GPUY
Pairlist

Update?
NForce (real)

Send coordinates
Force (bonding)
Force (reciprocal)

1. Computation intensive work : GPU• Pairlist• Real space non-bonded interaction
2. Communication intensive work : CPU• Reciprocal space non-bonded interaction with FFT• Bonded interaction• Exclusion list
3. Integration is performed on CPU due to file I/O.

Gather forces
Integration (VV2)
Finish?
Yend

NJung et al. JCTC 12, 4947-4958 (2016).

Our suggestion of HMR ratio
33

3

10103

1012 12103123HMR

10383101212833

31012 13321112211123 (2.5)

312112
11HMR(our suggestion)

10123331032

1112 1210122103103

3

2 or 3 (2.5)
11 or 10 (10.5)

283 (2.5)
2 or 3 (2.5)
13 or 12 (12.5)

1.

HMR ratio is 2 for XH1 type or ring case2.
3.

For XH2 or XH3, we assign HMR ratio 3 (CHARMM) or 2.5 (AMBER)
It should be noted that this does not confirm the reliability of 5 fs all the time. To avoid shake error all the time, 3.5 fs
with RESPA MTS could be better solution.

Jung et al. JCTC 17, 5312-5321 (2021)

21AccelerationofourMDsoftware,GENESIS
Kcomputer0.5

ster512(f)Performance (ns/day)

8PDYNputer

128512

2048

8192

GENESIS version 1 on K

32768

GENESIS version 2 on Fugaku

J. Jung (RIKEN)

Crowding (103.7 million)

SPDYN

4

100 million

21 billion

10.5131072 16384
2048

K computer
32768
4096

65536
32768
8192 131072
16384 262144

Numberof Nodes
Numberof CPU
cores
Jung, Mori et al., WIREs Comput. Mol. Science. 5, 310-323 (2015)

Jung et al. JCC 42, 231-241 (2021).

How accurate MD simulation with HMR?
Stable integration with a large time step

MD trajectories of a bond distance between nitrogen andcarbon atoms in the pyrrole ring (“d” shown in Fig. 4(a))
of the ligand binding domain of AMPA receptor. The
black line starts from 210 steps before the constraint error
using HMR. In contrast, HMR_xh1 gives a stabletrajectory of the distance, starting from the same structure.
(Inset) Two structures of distorted indole rings.

Diffusive motions of proteins and ligands

Translational mean squared displacements (left) and rotationalcorrelational function (right) of (a) PP1 and (b) CI2. (c) Lateral
diffusion coefficients of DPPC in 160DPPC systems in sevensimulations.
Jung et al. JCTC2417, 5312-5321 (2021)

Jung et al. JCTC 17, 5312-5321 (2021)

How accurate MD simulation with HMR?
Membrane protein dynamics (a: Calcium pump, b: Sec translocon)

Free-energy landscape of Chignolin in water25

MD simulations of SARS-CoV-2 S-protein on Fugaku

Interactions between atoms→ Force (F) is applied to atoms→ Atoms move

Spike Protein

Sugar chains
Water

Simulation on Fugaku

26

Mori et al. Biophys. J. 120, 1060-1071 (2021).

27

What can we study in biology(using a fast MD program onsupercomputers)?

Liquid-liquid phase separation (LLPS)

Banani et al. 2017 Nature Rev. Mol. Cell Biol.

Alberti et al. 2019 Cell

Names of LLPS in biology:•••

Biomolecular Condensates
Membraneless Organelles
Liquid Droplets

Multivalency: the General Principle of LLPS
Weak transient interactions between multivalent proteins drive LLPS.

Shin et al. 2017 Science

Computational Biophysics Research Team (Molecular Models)

TL: Yuji Sugita

Multi-scale molecular models in GENESIS

Hybrid QM/MM MD

Atomistic MDCoarse-Grained MDProtein

1 bead / amino acid

S. Takada (Kyoto U.)DNA

3 beads / nucleotide
Yagi et al., JCTC (2019)→ Electronic structures→ Enzyme reactions→ Material designs→ Vibration and spectroscopy

Yu et al., eLife (2016)→ Membrane protein→ Protein/DNA complex→ Cellular environments
M. Feig (Michigan State U.)

Tan et al. PLoS Comp. Biol. (2022)→ Slow dynamics→ Macromolecular interaction→ Data-driven simulations→ Cellular-scale simulations

Structure-based Coarse-Grained Models
for Biomolecules
Residue-level coarse-graining: ~10 atoms / CG particle

Residue-level CG models: protein-DNA systems•~500 particles

Protein: AICG2+
W. Li et al. 2014, PNAS.

•DNA: 3SPN.2C~2,000 particles

1 bead / amino acid

G. Freeman et al. 2014, JCP.

•RNA: Go-like~1,000 particles

N. Hori et al. 2012, JCTC.

3 beads / nucleotide: Phosphate, Sugar, BaseC.Tan et al. 2016, JACS.
K. Kamagata et al. 2020, NAR; 2021, NAR.

Protein-DNA: PWMcos

https://www.cafemol.org/

Prof. S. Takada(Kyoto U.)

C. Tan(Kyoto U. → RIKEN)

C.Tan et al. 2018, JCTC.
C.Tan et al. 2020, PNAS.
C.Tan et al. 2018, JCTC.

These CG models, which were originated from the structure-based
Go-model, have been developed and used in the CafeMol program byProf. Shoji Takada’s group.
We (RIKEN) are collaborating with the Takada group (Kyoto U) todevelop CG MD simulations in GENESIS software.

CG models in GENESIS MD software
GENESIS supports CG models of protein, RNA, DNA and Lipid
C. Tan (RIKEN)

D. Ugarte La Torre(RIKEN)← New, for lipids!

Summary and Perspectives
Macromolecular Crowding EffectSignal Transduction Pathway
Multi-scale models(AA, CG, QM/MM)
Enhanced Sampling(REMD, REUS, gREST,GaMD, etc.)
Okamoto, Hibino, Sako. BBA – General Subjects 1864
(2020) 129358.

Protein Condensates
Larger biologicalsystems in
atomistic or
residue-level MD
simulations
Zimmerman and Minton, Ann. Rev. Biophys. Biomol. Struct. (1993) 22: 27-65.
Zhou, H.X., Rivas, G., Minton, A.G., Ann. Rev. Biophys. (2008) 37, 375-397.

Shin et al., Science 357, 1253 (2017)
Wang et al., Sci. Adv. 9, eadf6205 (2023)

Acknowledgement
Collaborators

RIKEN:• Dr. Jaewoon Jung (RIKEN CPR/R-CCS)• Dr. Chigusa Kobayashi (RIKEN/R-CCS)• Dr. Cheng Tan (RIKEN/R-CCS)• Dr. Diego Ugarte La Torre (RIKEN/R-CCS)• Dr. Ai Niitsu (RIKEN/CPR)• Dr. Isseki Yu (RIKEN/CPR→Maebashi Inst. Tech)• Dr. Kiyoshi Yagi (RIKEN/CPR)• Dr. Takaharu Mori (RIKEN/CPR→Tokyo U. Sci)
For Bacterial cytoplasm• Prof. Michael Feig (Michigan State University)
For CafeMol• Prof. Shoji Takada (Kyoto University)

RIKEN Sugita Lab
BDR @ Kobe (2022)

Research grants:
•••••

•

MEXT “Priority Issue 1 on post-K computer”
MEXT “Program for Promoting Researches on theSupercomputer Fugaku”
MEXT Grant-in-Aid for Scientific Research (S)(19H05645)
MEXT Grant-in-Aid for Transformative Research Area
(A), Cross-Scale Biology (21H05249)
RIKEN Pioneering Projects (Dynamic StructuralBiology/Glycolipidologue Initiative/Biology of
Intracellular Environments)
Research grants from RIKEN R-CCS and BDR

R-CCS @ Kobe (2022)

Computer resources:• RIKEN HOKUSAI GreatWave, HOKUSAI
BigWaterfall, K computer, and Fugaku
supercomputer• HPCI supercomputer resources(hp170254, hp180201, hp180274,hp190097, hp190181,hp200129,
hp200135, hp210107, hp210177,
hp220087).

CPR @ Wako (2023)