Development of GENESIS on Fugaku supercomputer and its
application of Spike protein on the surface of SARS-CoV-2 insolution
2020年度第1回計算科学フォーラム
2020/11/24

Jaewoon Jung
Computational Biophysics Research Team,RIKEN Center for Computational Science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
2
2𝐩 𝑡 − ∆𝑡2𝑡𝑡 + ∆𝑡r 𝑡 + ∆𝑡𝐅(𝑡)1𝐩 𝑡 + ∆𝑡
2𝐅(𝑡 + ∆𝑡)

MD software GENESIS• GENESIS has been developed in RIKEN.
• It allows high-performance MD simulations on parallel supercomputers like K,Fugaku, Tsubame, etc.
• It is free software under LGPL license.

GENESIS developers

GENESISの開発の歴史

2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
年度
Start of the development
at RIKEN WakoKobe 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://www.r-ccs.riken.jp/labs/cbrt/

GENESIS 2.0beta: Fugaku Codesign

Continuous effort on the acceleration of GENESIS on Fugaku
FUGAKU Co-design
SystemFrom K to Fugaku

Application

Basic Arch.
DesignManufacturing

Collaboration

11.3 peta-flops

Start-up
Investigation
Development

X 35

GENESIS (2014)

400 peta-flops
GENESIS (2020)

x 125 (dt=2.5fs)
The acceleration was very important for flagship 2020 project.

Speed-up of GENESIS on Intel Xeon Gold 6148 (1 node, 40 cores)
120

Time (s)100Single precision/r-RESPA

80

Optimized kernel

60
40200

Pre-release

GENESIS
1.0

GENESIS
1.1.0

GENESIS
1.2.1

GENESIS
2.0beta

Benchmark calculation of ApoA1 (92, 224 atoms)

OthersLong-range NB
Short-range NB
Pairlist
Integration?

GENESIS2.0 or later

4ps MD

How to optimize GENESIS 2.0beta on Fugaku?
Algorithm1 (X)Algorithm2 (O)do i = 1, 16
do j = 1, 16
f(i,j)
end doend dodo ij = 1, 256
index i and j from ijf(i,j)
end doend do1.

In principle, there is no difference between Algorithm1 and Algorithm2 because of the same operation amount.

2.

On Fugaku, Algorithm2 has better performance than Algorithm1.

3.

The main reason is the operand waiting time for each do loop. On Fugaku, there are non-negligible waiting time before executing calculations in the do loop.

4.

To minimize the waiting time, the most inner do loop length should be long.

5.

If the most inner do loop length is small, it would be better not to vectorize when compiling.

do i = 1, 3f(i)end do!$ocl nosimd
do i = 1, 3f(i)end dof(1:3)!$ocl nosimdf(1:3)
f(1), f(2), f(3)New non-bonded interaction algorithm in GENESIS 2.0beta
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

Two Kinetic energy evaluations in MDIntegration
1
10! 𝑡 + ∆𝑡 = 𝐩! 𝑡 + 𝐅! 𝑡 ∆𝑡𝐩
2
210! 𝑡 + ∆𝑡 ∆𝑡𝐩
2𝐫! 𝑡 + ∆𝑡 = 𝐫! 𝑡 +𝑚!
1
10! 𝑡 + ∆𝑡 + 𝐅! (𝑡 + ∆𝑡)∆𝑡)
𝐩! (𝑡 + ∆𝑡) = 𝐩
2
2

Half-time step kinetic energy(𝐩! (𝑡))𝐾%*$$ (𝑡) = /2𝑚!
!&') (𝑁+ 𝑘, 𝑇 = 2 𝐾%*$$ (𝑡) −∆𝑡/ 𝐩! (𝑡) 5 𝐫! (.) 𝑡
6!&'Full-time step kinetic energy))
1
10! 𝑡 − ∆𝑡0! 𝑡 + ∆𝑡𝐩𝐩122𝐾"#$% (𝑡) = /+22𝑚!
2𝑚!
(!&'+ 𝑂 ∆𝑡 0) (𝑁+ 𝑘, 𝑇 = 2 𝐾"#$% (𝑡) +1. Full- and half-time step kinetic energies have the ∆𝑡 ) perturbation terms with opposite signs.
2. The perturbation terms are highly related to the high frequency vibrational motions.

∆𝑡/ 𝐩! (𝑡) 5 𝐫! (.) 𝑡
12!&'+ 𝑂 ∆𝑡 0

Our suggestion of accurate temperature
42𝑁! 𝑘" 𝑇 = 𝐾#$%& (𝑡) + 𝐾&'%% (𝑡) + 𝑂 ∆𝑡 (33
1. By combining full- and half-time step kinetic energies, we can obtain new temperature estimation which is accurate up to thethird power of ∆t.
2. The new temperature estimation does not require further computational cost.
3. With the new temperature estimation, we can enlarge the time step (∆t) up to 3.5 fs for short-range force and 7.0 fs forlong-range force. (In conventional methods, ∆t is limited up to 2.5 fs)

Similarly, we can estimate pressure in a more accurate way:
1X 𝑃 𝑡 = 2𝐾&'%% 𝑡 + 𝑉𝑖𝑟(𝑡)3𝑉However, the error could be negligible with grouppressure up to 2.5 fs

1O 𝑃 𝑡 = 3𝑉 2𝐾#$%& 𝑡 + 𝑉𝑖𝑟(𝑡)Ref: JCP, 148, 164109 (2018)
Ref: JCTC, 15, 84 (2019)

Strong Scaling of GENESIS 2.0beta on Fugaku
4 times better performance than
GENESIS1.0 on K using 1/8 nodes

11.9 ns/day(more than twice better
performance than NAMD)

Spike-protein in solution on Fugaku
Fastest →x 9.04x 6.82x 4.51x 2.96x 1.75Performance between Fugaku and Oakforest-PACS
Spike-protein in solution on Oakforest-PACS

Spike-protein in solution on Fugaku
Fastest → x 9.04
x 6.82x 4.51x 2.96x 1.7564 core is used for each node.

Entry of SARS-CoV-2 into Human Cell
A structural model of SARS-CoV-2• A key structural features of SARS-CoV-2 ⇨ Many spikes on the surface• Binding between Spike protein and a receptor in human cellPurpose of our simulations
6VXXRBDHow conformational changes happen ?
Which amino-acid residues are important ?

6VYB

Conformational changes

In Active form
S proteinDown

What is the role of glycans ?
How Spike protein can bind ACE2 receptor ?

Active formUpTwo purposes:• Prediction of conformational changes between the inactive and active forms.
• Elucidation of functional role of glycans in spike protein

generalized REST (gREST) method (theory)
Replica-Exchange with Solute Tempering (REST2):
T. Terakawa et al., J. Comput. Chem. 32, 1228-1234 (2011),
S. L. C. Moors et al., J. Chem. Theor. Comput. 7, 231-237 (2011),
L. Wang et al., J. Phys. Chem. B 115, 9431-9438 (2011).

𝛽6𝛽6𝐸2345) =𝐸88 +𝛽7
𝛽7REST/REST2'9)

𝐸8< + 𝐸<<In all replicas, solvent issimulated at room temp.

gREST

generalized REST (gREST):M. Kamiya and Y. Sugita, J. Chem. Phys. 149, 072304 (2018)

𝛽6𝛽6𝐸12345 =𝐸 +𝛽7 88
𝛽7𝐸8< + 𝐸<<

Solvent

Solvent-Solvent (vv)
Solute-Solvent (uv)
Protein“solute” temp. isdifferent in each
replica like REMD.

:9;“Solute” in gREST is flexiblydefined.
(ex) A part of a solute
molecule or a part of potential
energy function.

Solute-Solute (uu)LigandSolute = Ligand + Sidechainsat the Binding Sites

Conformational Changes of S-protein in solution
gREST simulations on S-protein• We used only Down form (6VXX) and could predict similar
structures of Up form.
⇨ Only chain B underwent the conformational change.
(N-glycans were not added in the gREST simulations.)

Summary of gREST simulations on Fugaku
Open (MD Snapshot)

Down (Initial Structure)

Up (cryo-EM)

gRESTInactive form

Active form ?

6VYB

Active form

• We predicted large conformational changes from Down to Open forms using Fugaku just for 3 weeks.
• The open form predicted using gREST is similar to Up form by cryo-EM.
• We would like to examine molecular mechanisms for the conformational changes.

Summary• We have developed MD software GENESIS and optimized GENESIS into K, Fugaku,Oakforest, and other supercomputers.
• Conformational changes from Down to Up were predicted using gREST simulations
on Fugaku.
• N-glycans on the surface of S-protein stabilize either Up or Down conformation,
forming different contact and hydrogen-bonded pairs.
• We hope to contribute the development of vaccines or anti-virus drug.

Acknowledgements• RIKEN CPR（理研開拓研究本部）
Dr. Hisham Dokainish
Dr. Takaharu Mori
• RIKEN R-CCS（理研計算科学研究センター）
Dr. Chigusa Kobayashi• Computational Resources
Fugaku, Oakforest-PACS• NIBIOHN （医薬基盤研究所）
Dr. Suyong Re• Computational Resources
Fugaku, Oakforest-PACS• All works are organized by Dr.
Sugita (RIKEN CPR, BDR, R-CCS)