AICS公開ソフト講習会 第4回「GENESIS」講義
理研、計算科学研究機構(AICS)
粒子系生物物理研究チーム
2014/06/17

今日の内容• 13:30 – 14:00 講義
 生体分子の分子動力学法
 GENESISの特徴• 14:00 – 14:10 休憩• 14:10 – 16:00 実習
 GENESISのコンパイル 実行例(1) 分子動力学法計算
 実行例(2) 分割I/Oツールを用いた解析
2

分子動力学法(Molecular Dynamics; MD)
粒子間の相互作用力を計算し、ニュートンの運動方程式を解く
2

mid ri (t ) Fi  t   
2dt

NU(r) ：相互作用エネルギーdU (r )
NdrirN：粒子の座標粒子の大きさ（粒度）で全原子モデルと粗視化モデルに分けられる
C原子

GENESISはどちらにも対応

粒子の運動から、構造の安定性、構造ゆらぎ、更にアンサンブルを計算し、自由エネルギー面（結合エネルギーなど）を解析できる
しかし、MD計算はとても時間がかかる…
タンパク質のダイナミクスをミリ秒計算するには、ニュートンの運動方程式を1012回計算する必要がある
例えば1stepに1msかかるなら、約31.7年必要！

3タンパク質のシミュレーション
fs-15

ps-12

s-6

ns-9

分子の振動
タンパク質
タンパク質＋水size

s0timelog10tFolding

二次構造変化

構造緩和

タンパク質の反応

1983

HP362(~104)

BPTI 1(~6x102)
AQP13(~105)

膜タンパク質＋水＋膜分子
超巨大
生体分子＋溶媒

ms-3

BPTI 4, Ubiquitin 5(~2x104)

2000
~2010
Ribosome6(~2x106)

HIV-1 capsid7(~6.5x106)

全原子モデル()は原子数
粗視化モデルFolding of
Globular proteins8

Drug export
from AcrB9, 10

1. Brooks et.al, PNAS,80,pp6571Biol., 9, e1003003 (2013)
6575 (1983)
7. Zhao et al., Nature, 497, pp6432. Duan et al., Science,282,646 (2013)
pp740-3744 (1998)
8. Clementi et al.. JMB, 298, pp
3. de Groot et al., Science,294,937-953 (2000)
pp2353-2357 (2001)
9. Yao et al., Nature Comm. 1:117
4. Shaw et al., Science,330,(2010)
pp341-346 (2010)
10. Yao et al., JACS, 135, pp74745. Piana et al., PNAS,110,7485 (2013)4
pp5915–5920 (2013)
6. Whitford et al., PLOS Comput.

Generalized Ensemble Simulation Systems(GENESIS)
1. 生体系において、効率的で精度の良い自由エネルギー計算
手法の開発が目的
2. 高並列計算 ‐ 京コンピュータなどでの超並列計算が可能（ただし、普通のPCクラスタでも動きます）
3. 巨大な生体系も可能
4. 全原子モデルのみでなく、粗視化モデル等異なる分子モデルへも応用できるアルゴリズムを採用
5. レプリカ交換法（拡張アンサンブル法の一つ）による自由エネルギー計算も可能
5

GENESIS 開発チーム
計算科学研究機構(AICS), 粒子系生物物理研究チーム
Project Leader :杉田 有治

Main developers: Jaewoon Jung
森 貴治(和光, TMS)
岩橋ー小林 千草松永 康佑
GENESIS web site
http://www.riken.jp/TMS2012/cbp/en/research/software/genesis/index.html
6生体分子の相互作用計算
1. 分子動力学法は膨大な回数のニュートン方程式を解く
2. 最も時間がかかる部分は粒子間の相互作用計算の部分である
3. 生体分子の相互作用計算は「力場(force field)」と呼ばれる経験的な関数で記述される
4. force fieldはおおまかに２つの部分に分けられる ‐結合性相互作用：

原子間の共有結合によるもの

非結合性相互作用：長距離相互作用（電荷、van der Waals力)によるもの

粒子数(N) に対してO(N2)で時間がかかるため、最も時間のかかる部分である
結合性相互作用• Bond• Angle• Dihedral

BondAngle
Dihedral

非結合性相互作用• Coulomb• van der Waals

非結合相互作用の高速化が必要
7非結合相互作用の高速化
1. 非結合相互作用はカットオフ長を導入することでO(N2)からO(N1)に変更できる
O(N2)全てのペアを計算

O(N1)カットオフ長以内のペアを計算

2. カットオフ長より長い部分のクーロン相互作用はフーリエ変換(FFT)により逆空間で計算する
erfc
2exp
G ⁄4Ucos G ∙4G44||
G

実空間(近距離)

逆空間
(長距離)

Self energy(定数)3. 1,2に加えて、非結合相互作用の高速化、高並列化は近距離項、長距離項の双方で更に行われる
8

GENESISの高速化・高並列化
最も時間のかかる非結合相互作用計算を高速化するため
GENESISでは下のアルゴリズムを新規に開発、組み込む• Inverse lookup table approach(Jung et al., J. Comput. Chem., 34:2412–2420, (2013))• Midpoint cell methods(Jung et al., J. Comput. Chem., 35:1064–1072, (2014))• Parallelization of FFT(Jung et al., (in submitted))(近距離相互作用計算の高速化）(近・長距離相互作用計算の高並列化）(長距離相互作用計算の高並列化）

上の３つに加えて、トラジェクトリなどの書き出し・読み込みを高速化するため並列I/Oスキームを独自に持つ
9

Inverse lookup table法
Jung et al., J. Comput. Chem., 34:2412–2420, (2013)

Lookup table法とは近距離相互作用は距離の関数であるため、カットオフ長までの距離の代表点での相互作用の値を計算、メモリ(table)に記憶しておく。
MDステップ中では、実際の距離に近いtableの値から内挿して求めるメモリに記録されている相互作用値

r

従来の方法

GENESISの方法

Jung et al., J. Comput. Chem., 34:2412–2420, (2013)より

従来の方法ではr2の線形関数、３次関数で内挿していた物を、GENESISでは1/r2の線形関数で内挿し、高速で精度の良い計算を可能とする
10

Midpoint cell 法

Jung et al., J. Comput. Chem.,
35:1064–1072, (2014)

ドメイン分割法（ほとんどのMDプログラムで採用）• 系全体をカットオフ長より長い一辺をもつセルで区切る

075663• エネルギー計算は隣接のboxのみを考えればよい(黄色のセル内の原子の相互作用は黄と緑のセル内の原子のみをカウントすればよい）• 通信回数が減少されるため、並列度が良くなる• ドメイン分割法で異なるセル間の相互作用計算をどのCPUコアに割り振るかが並列度に重要な問題になる
Midpoint cell法• 従来のmidpoint法では、２つの原子の中間地点のセルを受け持つコアが計算する• それぞれの原子が存在する二つのセル（８, 10とする）の中間のセル（7,11のどちらか）を受け持つコアが計算を担当する

Midpoint法1

Midpoint cell法

Jung et al., J. Comput. Chem.,
35:1064–1072, (2014)より• 通信回数が減少されるため、並列度が良くなる
1: KJ. Bowers et al, JCP 124, 184109 (2006)

11

System :(time step = 2fs)
京コンピュータでの計算(unpublished)
12

GENESISの内容物（1）MDプログラム• ATDYN (ATomic decomposition DYNamics simulator): atomic decompositionを使用
 レプリカ交換法が計算可能
 粗視化モデル（C‐GO, all atom GO)が計算可能
 わかり易いコードでユーザによる開発可能(SPDYNでの開発の前にATDYNでテストをすることも可能）
今回の実習では• SPDYN (SPatial decomposition DYNamics simulator):
こちらを使用 domain decompositionを使用
 超並列・高速なルーチン (Midpoint cell/3次元分割FFT/並列 I/O）
特徴

ATDYN

SPDYN

システムの分割法

原子分割

ドメイン分割New lookup table○○

レプリカ交換法○

×

粗視化モデル○×

3次元分割FFT

×

○

並列I/O

×

○

13

GENESISの内容物（2）ATDYN/SPDYN共通に計算可能な物• 最適化
 Steepest Decent法• Integrator
 Leapfrog
 Velocity Verlet• アンサンブル
 NVE NVT
LangevinBerendsen
 NPT
Langevin Piston(Isotropy of Simulation box:
Isotropic, Semi‐iso, An‐iso,XY‐fixed)• 拘束計算(constraint) SHAKE (Leapfrog)
 RATTLE (Velocity Verlet)
 SETTLE• FFT (PME) FFTE FFTW• Restraint functions
 Position
 Bond Angle Dihedral angle

14

GENESISの内容物（3）

その他の主なツール• trj_analysis
 トラジェクトリを解析するツール 距離・角度・二面角などが計算可能• crd_convert トラジェクトリを変換するツール トラジェクトリから一部分のみを抜き出して、新しいトラジェクトリを作成(例えば水分子を抜くとか、C原子のみにするなど）
 同時にfittingやRMSDなども計算可能• prst_convert
 GENESISのrestartファイルから並列I/O用のrestartファイルに交換• pcrd_convert
 並列I/O用計算で出された複数のトラジェクトリファイルを変換するツール
 crd_convertと同じ機能が使える
ダウンロードはGENESIS web siteから
http://www.riken.jp/TMS2012/cbp/en/research/software/genesis/index.html
15