創薬支援アプリ(ISLiMソフトウェア)講習会 ２０１5年9月1日
理化学研究所 東京事務所,
公益財団法人都市活力研究所

REIN-K講習会
宮下 尚之
近畿大学 生物理工学部
理化学研究所生命システム研究センター

理化学研究所

計算科学研究機構

目次

1. レプリカ交換分子動力学法の基礎
2. REIN-Kプログラムの原理
3. レプリカ交換分子動力学法を使った実際の研究例
4. REIN-Kの使い方

Naoyuki Miyashita, Suyong Re, and Yuji Sugita,“REIN: Replica-Exchange INterface for simulating protein dynamics andfunction.”International Journal of Quantum Chemistry,
115(5), 325–332. (2015).
http://doi.org/10.1002/qua.24785
Volume 115 | Issue 5 | 2015
Issue 5 (March 5, 2015)

Research in Systems Neuroscience
www.q-chem.org

内表紙ですがSpecial Issue: Theoretical Chemistry in Japan
Guest Editors: Hiromi Nakai and Takao Tsuneda

1. レプリカ交換分子動力学法の基礎

タンパク質の立体構造はデータバンク(PDB)に登録されている

年々、登録数が増加	  年

この調子では今年は1万件を超えるかもしれない	  現在、タンパク質の構造は10万件の登録がある。
http://pdbj.orgもしくはhttp://rcsb.org

登録数

PDBJより5タンパク質の立体構造（形）の基本 アンフィンセンのドグマ

wikipediaより1972年ノーベル化学賞
リボヌクレアーゼ分子の研究

リボヌクレアーゼAの実験 ：一度、変性剤（尿素）で立体構造を壊した後変性剤を取り除く（天然条件）

タンパク質の機能が復活する事を発見

タンパク質の天然立体構造はそのアミノ酸の１次配列により一意的に決定づけられるタンパク質の天然構造は全系におけるギブス自由エネルギーが最小の構造に対応する
自由エネルギー的に安定ファネルモデル(Peter Wolyness)

取り得る構造の種類
タンパク質の	  フォールディング

変性構造
天然構造との近さ

自由エネルギー
天然構造

天然構造を見つける	  ＝最適化問題タンパク質の天然構造は全系におけるギブス自由エネルギーが最小の構造に対応する
自由エネルギー的に安定最適化問題・探索問題•

タンパク質の構造予測や、リガンド­タンパク
質間の結合予測、構造変化予測などがしたい
最適解を見つける探索問題•

反応座標に対する自由エネルギーを描きたい
自由エネルギー計算

探索空間が広大

サンプリング法サンプリング• 意味：標本抽出• ある母集団（調査対象全体）からランダ
ムに標本抽出すること• → 全体のごく一部を調べるだけで、大きな母集団の正確な情報がつかめる。

• 非常に便利な方法

サンプリング法
レプリカ交換分子動力学法• 分子動力学法＋• モンテカルロ法

分子動力学法• 力に沿った運動に従った（古典力学に従った、
ある意味決定論的）、時系列データを得る事ができる方法（ダイナミクスが追える）。

• アンサンブル発生装置
決められた条件（アンサンブル）下のデータを得る事ができる
サンプリングデータとして利用できる。

タンパク質の構造予測など

Generalized
ensemble simulation
ここでは、構造揺らぎ・構造配置（ドッキング）・構造変化も含めた広い意味で

天然構造

変性構造

タンパク質の評価関数

力場
タンパク質の評価関数の最適化問題、構造分布

形（構造）
12

タンパク質の評価関数（力場）
ビーズ１個は原子１個バネだと考える1b

クーロン力とLJ

2

χ
共有結合

θ

3

rij

4

5
H (q,r) = K + V
V = Vinternal + Vnonbond

力学Vinternal = ∑ K b (b − b0 )2 + ∑ Kθ (θ − θ 0 )2
bonds

angles+ ∑ K χ [1+ cos(n χ − σ )]dihedrals

Rmin,ij 12 Rmin,ij 6 qi q j
Vnonbond = ∑ (4εij [(
) −(
) ]+
) 電磁気学rij
rij
εD rijnon-bonded
atom pairs

L-J項

クーロン項

13

さまざまな力場と力場パラメータ
力場(FF)

version

CHARMM

C19, C22, C27, C27CMAP, C36 など

Lipid-protein interaction

AMBER

Parm94, Parm99,
Parm99SB, 03 など

Protein, Glycan

GROMACS

43a1, 53a6, 53a7 など

Protein, Lipid

その他

OPLS FF, encad

得意分野

など

•評価関数が微妙に違う•パラメータの決め方が違う•得意不得意がある•最近はどれも結構良いよく使われる MD パッケージMD Packages

HybridSpeed

多次元REMD

on K

多機能
○

○

×

○
GPGPUで速い

×

×

REMD

CHARMM

[1]

×

Slow

AMBER

[2]

×

Medium/Fast

NAMD

[3]

○

Fast

○
○
△*
Charmの高速版(GPGPU版もある）

GROMACS

[4]

○

Fast

○
○ GPL (GPGPU版も)△
1CPU性能高い/高速,

DESMOND

[5]

○

Fast

高速 ○

×

×

TINKER

[6]

×

Slow

×
×
コードが読み易い×

MARBLE

[7]

○

Fast

○ (Native)
××
高速、京コンピュータに最適化

GENESIS

[8]

○

Fast

高速、京コンピュータなどに最適化[1] B. R. BROOKS et al., J Comput. Chem 30: 1545–1614, 2009[2] D. A. Pearlman et al., Comput. Phys. Comm. 91, 1-41, 1995
[3] J. C. Phillips et al., J. Comput. Chem. 26:1781-1802, 2005
[4] H.J.C. Berendsen et al., Comput. Phys. Comm. 91:43-56,
1995[5] K. J. Bowers et al., In SC ’06: Proceedings of the 2006 ACM/IEEE conference on Supercomputing, New York, NY,USA, 2006.
ACM Press.
[6] J. W. Ponder and F. M. Richards, J. Comput. Chem., 8,1016-1024, 1987.
[7] M. Ikeguchi, J Comput Chem 25: 529–541, 2004
15
[8] J. Jung et al., (2015)• それぞれのパッケージに特徴・特色• 機能も違う• 計算スピードも違う• 力場も違う（近年、対応できる様になってきた）

分子動力学法@300K

EnergyNativeX

DenatureX
短いシミュレーションだと初期構造に依存正しいアンサンブルが得られない
Conformation

ハミルトニアンレプリカ交換法
Adaptive Tempering法

Wan-Landau法

アンブレラサンプリング

Meta-dynamics法

拡張アンサンブル法
レプリカ交換法

マルチカノニカル法

Tsallis統計を用いた方法
レプリカ交換分子動力学法
Adaptive Biasing法

レプリカ交換法はMCMC法
普通のモンテカルロ法：
静的モンテカルロ法

x1 → x2 →..... → xi → xi+1 →•マルコフ連鎖モンテカルロ（動的モンテカルロ法）•手順：乱数をふるが、一つ手前の情報を用いて、次の状態を得る方法

普通のモンテカルロ法

システムの次元が大きくなると、相対的にターゲットが小さくなり、ダーツがヒットしにくくなる

タンパク質の構造予測：多次元系での構造探索

→マルコフ連鎖モンテカルロ法• どこからか、出発して、だんだん確率の大きい方に移動し、確率が大きくなったら、その近傍を探索する方法がよい• 更に、並列で探索(REMD)

レプリカ交換分子動力学 (REMD)
Replica #4

Not Accept

AcceptAcceptReplica #3• 同じシステムを複数用意（レプリカ）• 異なる温度（パラメータ）を与える

AcceptNot Accept
Replica #2• MDを実行

Accept• あるステップ毎にメトロポリスクライ

Acceptテリアに従って隣合う温度を持つレプ
リカ間で交換評価

Replica #1
AcceptAcceptAccept

時間[REM]	  K.	
  Hukushima,	
  and	  K.	
  Nemoto,	
  J.	
  Phys.	
  Soc.	
  Jpn.	
  65,	
  1604-­‐1608	
  	
  (1996)	
  
[1]	
  Y.	
  Sugita	
  and	  Y.	
  Okamoto,	
  Chem.	
  Phys.	
  Lett.	
  314,	
  141-­‐151,	
  (1999)	
  
[2]	
  Y.	
  Sugita,	
  A.	
  Kitao	
  and	  Y.	
  Okamoto,	
  J.	
  Chem.	
  Phys.	
  113,	
  6042-­‐6051	
  (2000)

6

MCMCの条件• どんな初期分布P から初期状態x を選ん
00

でも、t→ でPt(x)は定常分布P(x)に収束する。

定常マルコフ過程

許されている道路
１日目：1-2, 3-4
２日目：2-3, 4-5

１を出発した人がそれぞれの町にいる確率
Associate probability

保守的な５つの町で商売をしている塩商売人の話。
一本の道でつながっており、関所と門で交通制限されている５つの町の往来を考える。
町には受理率40%の関所がある。町には南北２つの門（関所）があり、１日毎に開門される門が替わる。１番町は北門、５番町は南門のみで１日毎に開け閉めがある。関所破りは厳罰。

5
4
3

1/N

2
1

35日くらいで全ての町を満遍なく回れる。

詳細釣合• 先の条件を満たす条件• 遷移確率R(x|y): 状態yからxを生成
R(x|y)f(y) = R(y|x)f(x)R(x|y)f(y)：事前確率メトロポリス法
同じ確率ならT(x|y)=1/4xxx1 → x2 →..... → xi → xi+1 →

x
x

• 手順
• 候補を選ぶ（y→x) Trial: T(x|y) dxT (x | y) = 1∫• 受理するかどうかを決める
y

受理の確率(Acceptance Prob.): A(x|y)
遷移確率：R(x|y)=A(x|y)T(x|y)
A(x|y)T(x|y)f(y)=A(y|x)T(x|y)f(x)A(x|y)=T(x|y)f(y)=A(y|x)T(x|y)f(x)T (y | x) f (x)A(x | y) =A(y | x)T (x | y) f (y)A(x|y)=r(x|y)A(y|x)
A(x|y)=1 => A(y|x)<1、r(x|y)>1
A(y|x)=1 => A(x|y)<1, r(x|y)<1
この場合A(x|y)=r(x|y)
Acceptance Probability:A(x|y)=min[1,r(x|y)]r(x|y)>1の時は、必ず受理。r(x|y)<1なら、乱数をふって、r(x|y)>randなら受理、
r(x|y)<randなら不採用。

レプリカ交換法の受理確率
今、M個のパラメータで、M個のレプリカがあるとする。全体の分配関数は、

Z =Π

M
m=1

確率分布関数はf (X; β ) = ΠM
m=1Z(β m )exp(− β m H (X m ))Z(β m )

詳細釣合の式を思い出す。
R(x|y)f(y) = R(y|x)f(x)T (y | x) f (x)r(x | y) =T (x | y) f (y)β = 1 / kBT

隣り合うパラメータ間の交換は、T(x|y)=T(x|y)

今、H(p,q)=K(p)+U(q), ただしpとK(p)はMDの時
ここではMCの時を仮定。H(x)=U(x)とする。
f (x) exp(− β n H (yn ))exp(− β m H (xm ))=f (y) exp(− β m H (ym ))exp(− β n H (xn ))= exp(− β n (H (yn ) − H (xn )) − β m (H (xm ) − H (ym ))) 尚、MDの時は温度ス
ケーリングを考慮す= exp(− β n (U(y) − U(x)) + β m (U(y) − U(x)))
ること。
= exp[−(β n − β m )(U(y) − U(x))]= exp[−Δ],Δ = (β n − β m )(U(y) − U(x))

A(y|x)=min[1,exp(-Δ)]

例：タンパク質の構造予測
通常のMDでは構造探索あまりできない
通常の分子動力学法（MD）

P(E) 30℃

E
エネルギー空間、局所最小値に捕まる
正しい構造が見つからない29

通常のMDでも高温だといろいろな構
造をサンプルできる。しかし
P(E)

70℃

E
エネルギー空間、局所最小値に捕まらないが知りたい温度の構造が得られない
30

レプリカ交換分子動力学法(REMD)
温度上のランダムウォーク
レプリカ交換法P(E)40℃ 60℃
30℃ 50℃ 70℃

5パラメタ 1 2 3 4 5
5レプリカあるステップ毎にメ
トロポリスクライテ
リアに従って、隣り合う温度を交換

それぞれ異なる温度をもつ
複数のレプリカを用意しMD計算

エネルギー空間、構造空間を幅広くサンプリングできる
統計処理を行い知りたい温度での構造決定
31

E

バリアを乗越える事ができる
温度の異なるMDを行い、ピコ秒毎に温度を交換

@100℃
@90℃
@80℃
@70℃
@60℃
@50℃
@40℃
@30℃

32

多次元レプリカ交換法(MREM)
温度•

温度だけでなく、距離や角度、二面角、その他のパラメー
タ交換にも拡張可能。
たとえばReplica-exchange umbrella sampling(REUS)[2]•

それらを組み合わせて多次元レプリカ交換
Multi-dimensional REMD (MREM) [2]•

ただし、レプリカ数が増える

距離
角度

距離

replica[2] Y. Sugita, A. Kitao and Y. Okamoto, J.
Chem. Phys. 113, 6042-6051 (2000)

京コンピュータの沢山のノードを利用した計算が実現できる。

温度

2. REIN-Kプログラムの原理

REINプログラム

多次元レプリカ交換インターフェースプログラム
Replica-Exchange INterface program for K
REIN-K

REINはバイナリの外部プログラムを制御します
REIN

Control

MD program

REIN：親プログラム
バイナリ外部プログラム(MD)：子プログラム• 子プログラムの実行• 子プログラムの終了判定• 子プログラムのoutputの読込み• エネルギー等の情報取得• レプリカ交換実行• 子プログラムの次のinputを作成

N. Miyashita et al., IJQC, 115(5), 325–332. (2015)

REINは既存のバイナリの外部プログラム(MD)を子プログラムとして複数同時に起動する
２つの仕組みMPI-2 version

MPI-1 version

利用技術

MPI̲COMM̲SPAWNSYSTEM関数

対象コンピュータ

超並列コンピュータ
K, FX10, PCクラスタ

PCやPCクラスタ、小規模並列コン
ピュータ、GPGPU等も含む

制約

MPI-2(動的プロセス)対応, NFS(共有ディスク)対応, 子プログラムは
MPIでコンパイルされている事。

MDはノード内並列のみ対応
MPIを使わない特徴

１ノード分、MD制御プロセスとしMD制御プロセスは各々のレプリカて使用する。レプリカ数 x MDの並の使用するノード上で起動する。
列数 + 1ノード使用する。

計算速度

新しい機能なのでシステムのmpi̲comm̲spawnの対応/整備状
況に依存。大規模系で大規模計算
向き。

速い

レプリカ交換の２つの仕組みを用意
MPI-1 version
一般的なPCクラスタ、MDは1ノード以内
MDはMPI並列でない場合

MPI-2 version
高度なPCクラスタ、MDは1ノード以上も可
MDはMPI並列化もしくはHybrid並列化されていること

N. Miyashita et al., IJQC, 115(5), 325–332. (2015)

MPI-2 version
MPI̲COMM̲SPAWN•

動的プロセス生成•
MPIのプロセスからMPIを実行する•
MPI-2に対応していればどの計算機でも動く•
スケジューラが制限をかけている場合がある•
Master processが必要•
Child Jobはnodeをまたぐ事ができる

制限の無いPCクラスタ

Master
Process

1
1

1

2
2

2

5

1

4

5

2

2

3

5

6

2

3

3

5

6

4

6

6

5

１つのノード内に2つ以上のJobは入らない

5

4

4

2

1

5

5

3

2

6

3

3

2

6

6

Master process: 1 coreChild Job: 12 core x 2 Job生成

京コンピュータ, FX10

6

5

6

4

1

4

3

１

1

京コンピュータ, FX10
44nodeのMPI/OpenMP化

Master process: 1 coreChild Job: 3 core x 6 Job生成
3

4

京コンピュータ, FX10でMaster

例：4 core x 8 nodeのシステム

3

1

Child Jobはノードをまたぐ事ができる

1

1

2

2

2

2

1

1

2

2

2

2

1

1

1

1

2

2

1

1

1

1

2

2

MPI-1 version
system 関数を用いたレプリカの制御
•

•
•

•

•

MPIのプロセスからSystem関数を実行
大抵のPCクラスタで動く
しばしばスケジューラが制限をかけている場合がある
Master processは必要ないが、Child Jobの実行プロセスの１
coreが制御にも使われる。
現状ではMPIプロセスを生成できない（nodeをまたぐ事ができない）
例：4 core x 8 nodeのシステム

Child Jog: 2 core x 6 Job生成
1

1

3

3

2

2

4

4

5

5

6

6

Child Jog: 3 core x 6 Job生成

Child Jog: 4 core x 6 Job生成

Replica-exchange Interface Program (REIN)
多次元レプリカ交換シミュレーションができる
Control

REIN

MD program
MDのコード

パッケージ

MD

言語

MDの制御

REIN version
1.0.4

NAMD
MARBLEFORTRAN95,
03/MPI/
OpenMPbinary run[mpi_comm_spawn orsystem function ]

修正

MREM

K対応

○

○

no

MARBLE: Optimized on K.
NAMD: Currently flat MPI version are available on K.

REIN:
MachineCompiler

K computer (sparc), FX10 (sparc), PC cluster (x86_64), RICC(富士通のスパコン)
Fortran (Intel, fujitsu, (gfortran))

MPI-2 (Open MPI, Fujitsu MPI)/OpenMP

多次元レプリカ交換（MREM）：
REMD
軸

REUS

温度

距離

角度

二面角

Atom

Group

Atom

Group

Atom

Group

NAMD

○

○

○

○

○

○

○

MARBLE

○

○

○

○

×

○

×

京コンピュータとSCLS FX10

京コンピュータ

SCLS FX10#Core/node

8 core16 core#Node88,128 node

48 nodeTotal core

705,024 core

768 core

File system

Staging/Shared

Shared (home)

CPU

2.0GHz(sparc64VIIIfx)

1.6GHz(sparc64IXfx)

性能

10PF

10.1TFGPGPUを搭載するPCクラスターやPCクラスター
TSUBAME	  2.5	
  	
  	
  	
  	
  5.7	  PFLOPS	  	
  @Tokyo	  Ins8tute	
  of	
  Technology	  東工大web pageより1442$node$x$($12$cores$+$3$GPGPU$(K20X)$)
NAMD2.10$nighty$version$4/24

京コンピュータのディスク

ステージングによるデータ転送

Frontend shared HDD
計算ノード shared HDD
X 通信不可計算ノード local HDD
動的プロセス生成では、shared HDDのみ利用可(REINで利用)
SCLS FX10ではステージングは無い

REIN プラットフォーム（計算科学）
モデル

REINはbinary形式のプログラムを京コンピュータやPCクラスタ(GPGPU対応)上で複数起動機能

任意のパラメータを変えて並列に子プログラムを起動
し、ある周期毎に任意のパラメータの交換ができる

京コンピュータで動くnamd（とmarble）に対応
対応モジュール
PCクラスタ(GPGPU対応可)で動くnamdに対応(対応子プログラム)
インターフェースモジュールの開発により他にも対応
REMDの対応機能

現時点では、パラメータの交換対象は、温度とアンブ
レラポテンシャル

今後

将来的には対応するMDの種類と、レプリカ交換計算の種類を増やします。

利用目的・対応状況（分子シミュレーション）
機能

REINは既存の分子動力学(MD)プログラムを利用して、多次元レプ
リカ交換(MREM)シミュレーションを行う為のプログラム対応機種

REINは京コンピュータ, FX10とPCクラスタに対応, TSUBAME
でも利用可

対応スケジュ
ーラ

富士通のスケジューラ、Grid Engine

計算対象は、生体高分子の構造予測, 分子間の結合自由エネルギー
想定する研究計算、構造変化揺らぎの予測です。
現時点ではMDはNAMD2を利用します
利用するMD(MARBLEもKでは対応はしていますが、現時点ではNAMD2でのプログラム
利用を推奨します。)
レプリカ交換の対応機能

現状、温度レプリカ交換、距離、角度、二面角のレプリカ交換ア
ンブレラサンプリングに対応

今後

将来的には対応するMDの種類と、レプリカ交換計算の種類、対応
機種・スケジューラを増やします。

REIN-Kを用いた
Penta-Alaninの構造予測の例

N. Miyashita et al., IJQC, 115(5), 325–332. (2015)

レプリカパラメータ：温度とEnd-to-endの距離

T-REMD

温度

距離

REUS

MREMD(2D)

298-500K[32 replica]

300K

300-350K[4 replica]

ー

4-17Å
[12 replica]

4-17Å[8 replica]

1.0 kcal/mol/Å2

1.0 kcal/mol/Å2

Harmonic
restraining
potential

レプリカ数

32

12

32

Exchange
Trials

10,000[Total 640ns]

10,000[Total 240ns]

8,000[Total 512ns]

2 ps for 1 exchange
力場：charmm c36 CMAP
MD: NAMD2.9, MPI-2 version on K-computer

T-REMD

REUS

MREMD

N. Miyashita et al., IJQC, 115(5), 325–332. (2015)

MREMDの収束性は非常に良い
WHAM @ 300K

赤：T-REMD
緑：REUS
青：MREMD

REIN-Kを利用した生体分子構造予測シミュレーション

糖鎖の構造予測Re, S., NM et al., 2011. Biophys
Journal, 101(10), pp.L44–L46

糖鎖の構造予測２
Nishima, W., NM et al., 2012. Journal Of Physical
Chemistry B, 116(29), pp.8504–8512.

TOM20の機能解析
Komuro, Y., NM et al., 2013. Journal Of Physical
Chemistry B, 117(10), pp.2864–2871.

膜タンパク質の構造予測

さまざまな生体分子の構造予測に利用できる52

3. レプリカ交換分子動力学法を用いた実際の研究例

膜タンパク質の二量体構造予測

アミロイド前駆体タンパク質の二量体構造予測
インターフェースの変化がAβペプチド生成に影響
Naoyuki Miyashita, John E Straub, D Thirumalai, and Yuji Sugita,
Transmembrane structures of amyloid precursor protein dimer predicted byreplica-exchange molecular dynamics simulations. Journal of the American
Chemical Society, 131(10), pp.3438–3439 (2009)32 replica

20 ns

30 ns

REINの応用例

糖鎖分子の構造予測
Suyong Re, Naoyuki Miyashita, Yoshiki Yamaguchi, Yuji Sugita, Structural Diversity and Changes in
Conformational Equilibria of Biantennary Complex-Type N-Glycans in Water Revealed by Replica-Exchange
Molecular Dynamics Simulation,

李博士(理研)Biophys. J. 101, Pages L44–L46 (2011)

REINでのREMD計算で糖鎖（N-glycan）の構造予測を行った。分子に付いた糖鎖の揺らぎは分子間の認識に重要である。実験の予想どおり
bisecting GlcNAcは揺らぎにおおきな影響を及ぼしている事がわかった。

90 replica

ゆらぐ

抑制

bisect無し

bisectありカポジウィルス由来タンパク質の免疫システムへの攻撃実験との共同研究
Mizuho Kajikawa, Pai-Chi Li, Eiji Goto, Naoyuki Miyashita, Masami Aoki-Kawasumi, Mari Mito-Yoshida,
Mika Ikegaya, Yuji Sugitate, and Satoshi Ishido,
The intertransmembrane region of Kaposi's sarcoma-associated herpesvirus modulator of immune
recognition 2 contributes to B7-2 downregulation. Journal of Virology, 86(9), pp.5288–5296 (2012)

梶川博士(昭和薬大)Pai-chi Lee博士（理研）

32 replica, 450 replica

4. REIN-Kの使い方

REINパッケージの構成

•

rein

•

supporting software
input̲builder
batch̲builder
remd̲convert３つのサポーティングソフトウェア

名称

機能

理由

input̲builder

REINのインプットを自動生成する

柔軟性の為、数多くのインプ
ットファイルが必要であるがそれらを一括で自動生成する

batch̲builder

バッチキュー用のスクリプ
ト生成と、REINのインプ
ットの修正を自動で行なう

レプリカ数やノード数、スレ
ッド数、ステージングの設定などの設定は非常に複雑な為

レプリカ毎やパラメータ毎
remd̲convertにトラジェクトリーファイ
ルを収集する

レプリカが多くなると、デー
タ収集が非常に複雑になる為サポーティングプログラムのしくみ
input_builder← ib.inp$ input_builder -h ctrl

インプットのテンプレート$ input_builder -h ctrl all 全てのオプションが表示される$ input_builder ib.inp

実行

テンプレートファイルを読んで、reinのインプットを作成• namdやmarble, rein.inpのテンプレートが入っている• 必要に応じてテンプレートを直す• 追加も同じ形式（形式）で作成できる

基本的には３ステップで実行
1) Input_builder ! rein.inp

reinのインプットを自動作成#input_builder ib.inp
2) Batch_builder ! rein.sh

ステージングなどの煩雑なバッチスクリプトを自動生成#batch_builder bb.inp
3) あとはバッチスクリプトを投入するだけ。
# pjsub rein.sh

Restart#batch_builder bb.inp restart（restartと付けるだけで自動生成）

remd_converter ! トラジェクトリをまとめる#remd_converter rc.inp
61

インストールREINパッケージの構成
2014/02/26時点でversion 1.0.4です。
REIN---- README.md: README file|-- AUTHORS.md: authors file|-- note.md: release notes|-- update.md: update|-- COPYING: GPL licence file|-- src/: main code of REIN|-- tools/: supporting software for REIN|-- examples/: For tutorial`-- manual/manual.pdf : manual

情報：
質問など：

@rein_devel http://twitter.com/rein_devel
yukimya+rein@gmail.com

REINのコンパイル１
京コンピュータ、FX10$ cd src$ ./conﬁgure.sh sparc k
京用のNAMDのバイナリは提供しますのでご連絡下さい。

PCクラスタなど(no spawn版)$ cd src$ ./conﬁgure.sh x86̲64 gfortran nospawn
NAMDのhomepageからDLしたバイナリがそのまま使えます。

PCクラスタなど(spawn版)$ cd src$ ./conﬁgure.sh x86̲64 ifort
NAMDをMPI用にコンパイルして下さい。

REINのコンパイル２
インストール$ make$ make installsrc/と同じ場所にあるbin/にreinという実行ファイルが格納されます

clean$ make clean
objectファイルなどが消去されます。

distclean$ make distclean
objectファイルだけでなくbin/に入っているreinも消去されます。

サポーティングソフトウェアのコンパイル
京コンピュータ、FX10$ cd tools$ ./conﬁgure.sh x86̲64 gfortran京, fx10のフロントエンドはx86̲64マシンなのでPCクラスタなど(no spawn版)$ cd src$ ./conﬁgure.sh x86̲64 ifort
ifortも使えます。gfortranも可。

インストール$ make$ make installbin/にinput̲builder, batch̲builder, remd̲converter, templates/
がインストールされます。clean, distcleanはreinと同じ。

実行の為の準備

REIN実行前の準備１

事前に、通常の分子動力学シミュレーションを行う
時と全く同じ準備をする。

1. NAMD: VMDなどでpdbファイルからpsfファイルを作成。（REINの初期設定の力場はC27+CMAP）
MARBLE: molxでmdatとcrdファイルを作成。
2. MDプログラムを用いてMinimizationする。
3. MDプログラムを用いて平衡化計算を行う。

REIN実行前の準備２
REMD/MREMシミュレーションのプログラムとイ
ンプットの準備をします。
ここではrun/というディレクトリ内でシミュレーションを行うとする。
ここでは、最も単純な方法を述べる。

1. rein/bin/にある全てのソフトウェアとtemplates/をrun/にコピーする。
2. 平衡化で得られた計算結果(coord, vel, xsc, xst)とinputのpdb, psfを、initial.XXX (XXXはそれぞれの拡張子)という名前に変えてrun/へコピーする。
3. parameterファイルもrun/へコピーする

REIN実行前の準備３
input̲builderを使ってreinのinputファイルと実行用バッチファイルを作成する。
1. $ cd run/
2. $ ./input̲builder -h ctrl > ib.inp
3. $ vi ib.inp
レプリカの情報とPMEメッシュの修
正
その他はbatch̲builderでも修正可能
4. $ ./input̲builder ib.inp > ib.out
5. rep.1/, rep.2/, .....というディレクトリと、ﬁle̲*.dというファイルが生成。rein.inpのひな形が生成。

例: Input ﬁle of Input builder
REINのインプットは複雑な計算にも対応している。
簡単にREINのインプットを自動生成する。
温度 - 距離
2D-REMD# control parameters in input_builder
rein_input
spawn

=

= rein.inp
yes[REMD]# rein input file nameex)

１次元目：この場合、温度

axis001= temperature# axis (TEMPERATURE/HARMONIC)

num_replicas001= 4# # of replicas for the axis

range001= 300.0 360.0# lowest, highest value

division001= exp# (EQUIV/EXP/CYCLE)

axis002= harmonic

num_replicas002= 4

range002= 2.0 8.0

division002= equiv

distance002= ADP CR 1-1 ADP CL 1-1

spring_const002= 2.0

Alanine Dipeptide

REIN/NAMD, NVT

２次元目：この場合、原子間距離

MD steps: 1000 steps (2ps)
Exchange: 100 exchanges

等間隔

Number of Replicas:
Temperature : 4 replicas
Distance: 4 replicas
Total: 4 x 4 = 16 replicas[MD]
md_program= namd# md program (NAMD/MARBLE)

md_exe_file= ./namd2# md exe filemd_temperature= 300.0# md temperature

md_steps= 1000# md steps / exchangesystem_size_x= 23.4# system size --automatically create pme grids

system_size_y= 23.2#system_size_z= 23.9#システムサイズを入れると自動的にPMEのメッシュの数を計算する
71

Range of Temperature: 300-360K
Range of Distance: 2.0-8.0Å$ input̲builder ib.inp
実演REINのインプット：サポーティングソフトウェアで自動生成
sample/ ------ rein.inp||--- file_rank.d(For REIN)

レプリカとパラメータの対応表|--- file_value.d(For REIN)

パラメータの具体的な値|--- file_on.d(For REIN)計算に用いるレプリカ（input_builderのoutputをそのまま使って下さい）|--- file_axis_ex.d(For REIN)(For REIN)

軸の順序（同上）|--- file_coef.d(For REIN)

インプットに対する温度比|--- file_ex.d(For REIN)

交換ルール

|
|--- initial.pdb(NAMD)

テンプレート構造|--- initial.psf(NAMD)

トポロジーファイル
パラメータファイル|--- par_all27_prot_lipid.prm (NAMD)|--- initial.inp(NAMD modified)

namdインプットファイル。パラメータにタグを付けている。

|--- initial.colvar [for REUS] (NAMD modified)

namdインプットファイル。束縛情報。同上||--- rep.1/ ----- rep.1.0.coor (NAMD)

初期構造

|

|--- rep.1.0.vel (NAMD)

初期速度

|

`--- rep.1.0.xsc (NAMD)

初期環境（体積、圧力）|--- rep.2/ ----- rep.2.0.coor (NAMD)||--- rep.2.0.vel (NAMD)|`--- rep.2.0.xsc (NAMD)|--- rep.3/ ----- rep.3.0.coor (NAMD)||--- rep.3.0.vel (NAMD)|`--- rep.3.0.xsc (NAMD)`--- rep.4/ ----- rep.4.0.coor (NAMD)|--- rep.4.0.vel (NAMD)`--- rep.4.0.xsc (NAMD)

自動生成して、必要な部分を修正すると良いREIN実行前の準備４

batch̲builderを使って計算機環境に合わせたreinのinputの書き換え、と実行用バッチファイルを作成する。
1. $ ./batch̲builder -h ctrl > bb.inp
2. $ vi bb.inp MPIやthreadに関する情報を記載する。
スケジューラに関する情報を記載する。
3. $ ./batch̲builder bb.inp > bb.out
4. rein.inpが完成し、rein.shが生成

例: input ﬁle of Batch builder
レプリカの数だけ、ファイル転送する必要があるが、その際のステージングやFTLなど面倒な記述を自動で作成する。
REIN: 16 core (2 node)
MD: 4 node x 16 replica = 64 node(512 core)# control parameters in batch_builder
rein_input= rein.inp

spawn

= yes

chain_job= no# REIN input file name

Total node: 66 nodeThread: 8 (hybrid MPI)

[BATCH]
batchfile_name= rein.sh# batch filename

job_scheduler= K# job scheduler type ([K]/GE/RICC)

queue_name= small

num_mpiproc_rein= 2# number of process for rein

num_threads_md= 8# number of threads in md

num_mpiproc_md= 4# # of mpi process for md

cpu_time= 01:00:00# cpu time

num_core_cpu= 8# number of core in cpu, for k = 8, for fx10 = 16, ...

replica_exchange= yes# ([YES]/NO): replica exchange / umbrella sampling

num_of_exchange_steps= 100# number of exchange steps# queue name, k:([SMALL]/LARGE),fx10: REGULAR..

[REMD]

Exchange steps :100 exchanges

現在、K (fx10) , RICC, SGE タ
イプのスケジューラに対応.
新規スケジューラもテンプレートの追加で対応できる.
$ batch̲builder bb.inp$ batch̲builder bb.inp restart74

REINの実行

REIN実行
1. $ pjsub rein.sh
2. $ ./batch̲builder bb.inp restart > bb.out2
REIN終了後、run/内に結果データの入ったSTGOUT̲DIR/が作成される。
batch̲builderでrestartをかけると全ての結果ファイルがresults/の中に格納される。
以上を繰返す。

実行後：トラジェクトリーをreplica単位かパラメー
タ単位にまとめ直す
1. $ ./remd̲convert -h ctrl > rc.inp
2. $ vi rc.inp
ファイルのあるディレクトリを指定
3. $ ./remd̲convert rc.inp > rc.out

例: input ﬁle of remd̲convert
解析の前に、結果のトラジェクトリをレプリカ毎かパラメータ毎にまとめると便利# control parameters in remd_convert

REINを3回リスタートした場
合の例[INPUT]
psffilereffiletrj_format= initial.psf= initial.pdb= DCD

directory001
directory002
directory003= ./results/result0000001_0000010= ./results/result0000011_0000020= ./results/result0000021_0000030

repeat001
trj_filename001
energy_trj001= 3= rep.{rep_no}/rep.{rep_no}.{step}.dcd= replica.trj

energy_column=# (DCD/MARBLE/PDB)

レプリカでまとめている。

# column of output energy, ex. = 1 2-4, empty: all column[OUTPUT]
out_pdbfile
trj_filename
trj_format= out.pdb= out.{type}.{number}.{extension}= DCD# (DCD/MARBLE/AMBER/TINKER/PDB/GROMACS)[CONVERT]
type
rank

= REPLICA=# (REPLICA/CONDITION)# ex. = 1 2-8, empty: output all rank

[SELECTION]# select_atom# select_atom# mole_name001# mole_name002= all= molname:protein | molname:lipid= protein P1:1:TYR P1:5:MET= lipid
OLEO:PCGL:OLEO

77$ remd̲convert rc.inp

REIN連続投入
京、FX10の場合：tools/auto̲submitter/auto̲submitter.plを用いる。

••$ vi auto̲submitter.pl$ perl auto̲submitter.pl

PCクラスタの場合、２通り••

auto̲submitter.plの利用
chain̲job = yes (bb.inp)とすると
rein.shの最後に１行加わる。

...
${MPIEXEC} ${EXEC}./batch_builder bb.inp restart...
${MPIEXEC} ${EXEC}./batch_builder bb.inp restart${MPIEXEC} ${EXEC}./batch_builder bb.inp restartOutput format (replicat.trj)#========= EXCHANGE_STEP = 1
1
1
290.14150
22309.04970
3
3314.56290
44313.8316055304.91160
6
6309.8476077299.44280
88288.62440#========= EXCHANGE_STEP = 21
2
306.76840
21
317.18100
34310.83750
43
284.34980
56294.21380
6
5336.6426078296.55280
87299.56080-1624.33470-1677.37840-1658.91060-1657.91540-1642.47860-1635.80170-1685.76920-1683.47140

370.15570
394.27840
401.31190
400.37900
388.99900
395.29630
382.02210
368.22020-1994.49040-2071.65680-2060.22250-2058.29450-2031.47760-2031.09800-2067.79130-2051.69170-1664.28540-1651.98420-1629.92330-1655.03960-1680.86470-1606.88210-1698.13390-1662.84080

391.36800
404.65210
396.55920
362.76680
375.35110
429.48070
378.33510
382.17260-2055.65340-2056.63630-2026.48240-2017.80640-2056.21580-2036.36280-2076.46900-2045.01340

Replica#, condition#, Temperature, ( Distance etc.. ), Total, Kinetic,Potential, Extra#========= EXCHANGE̲STEP = 1111322.72880
28333.08190
3
6355.33950
4
5321.22060
52351.72040
6
9371.11030
7
7284.53280
8
4276.94870
93
334.14340#========= EXCHANGE̲STEP = 126.61182
7.132126.55514
6.67915
6.10475
7.01426
7.15645-1643.34560-1575.02380-1523.72050-1597.99230-1554.42320-1485.12630-1680.79330

411.72990
424.93810
453.33390
409.80570
448.71660
473.45380
363.00030-2055.07550-1999.96190-1977.05440-2007.79800-2003.13980-1958.58020-2043.79360

6.52260
7.480200.09900
0.05150
4.81970
7.943107.386406.584866.26521-1677.11400-1528.19240

353.32460
426.29240-2030.43860-1954.48480

0.086205.3308079

remd_convertで得られたOutputのformat (out.xxx.#.enetrj)
1
2
3
4
5
6
7
8
9
10

1
1
1
1
1
1
1
1
1
1

1
2
3
4
4
3
3
4
4
4

297.75270
317.02580
310.59180
315.72080
319.46520
321.01360
300.87790
298.17360
299.28800
295.37950-1671.67350-1639.88920-1617.55910-1601.63910-1626.64460-1633.93310-1678.91540-1661.42740-1656.02570-1657.20690

379.86590
404.45410
396.24580
402.78920
407.56620
409.54160
383.85290
380.40290
381.82470
376.83830-2051.53940-2044.34340-2013.80490-2004.42820-2034.21090-2043.47480-2062.76830-2041.83030-2037.85040-2034.04520#, Replica#, condition#, Temperature, ( Distance etc.. ), Total, Kinetic,Potential, Extra

1
2
3
4
5
6
7
8
9
10

1
1
1
1
1
1
1
1
1
1

1
1
1
1
1
2
2
2
2
1

286.98240
312.47330
289.37010
296.36820
289.88810
314.75100
330.75970
332.70610
330.66870
308.05460

6.37033
6.44694
6.00928
6.238496.28432
6.02690
6.421436.61133
5.60459
6.21432-1717.17430-1653.73810-1643.07150-1679.83210-1653.50920-1618.84470-1606.04790-1573.92340-1552.42150-1668.87810

80366.12540
398.64610
369.17160
378.09960
369.83240
401.55200
421.97540
424.45860
421.85930
393.00890-2083.29970-2052.38420-2012.24300-2057.93170-2023.34170-2020.39670-2028.02330-1998.38200-1974.28080-2061.88700

5.67950
5.94070
4.52790
5.24390
5.39340
4.581105.853106.52090
3.39190
5.16590

まとめレプリカ交換インターフェースプログラム
Replica-exchange interface program for K
REIN-K
京コンピュータで大規模な多次元レプリカ交換シミュレーシ

プログラム開ョンを実施する事を想定したプログラム
発の目的
様々なMDプログラムを利用したレプリカ交換をサポート
ファンド

次世代計算科学研究開発プログラム 次世代生命体統合シミュレーション研究開発プロジェクト(ISLIM)にて作成

その他ファン
理化学研究所生命システム研究センター(QBiC)
ド
SCLS（HPCI戦略プログラム 戦略分野1）
講習会
RIKEN情報基盤センター
ダウンロードバイオグリッドセンター関西、バイオスーパーコンピューティ
サイトなど
ング研究会

REINの情報
AuthorNaoyuki Miyashita and Yuji Sugita

ライセンス

ライセンス：GPL version 3

引用

Naoyuki Miyashita, Suyong Re, and Yuji Sugita,“REIN: Replica-Exchange INterface for simulating protein
dynamics and function.”International Journal of Quantum Chemistry,
115(5), 325‒332. (2015).
Naoyuki Miyashita, Yuji Sugita, M-2: Replica-exchangeinterface program (REIN), http://www.islim.org/islimdl̲e.html

Main developer& サポート近畿大学

生物理工学部

宮下尚之

e-mailyukimya+rein@gmail.com

Twitter@rein̲devel http://twitter.com/rein̲devel

謝辞

•
•
•
•
•
•
•
•

杉田有治 主任研究員, 理研, 理研QBiC, 理研AICS
池口満徳 准教授, 横浜市大
李秀栄 研究員, 理研
二島渉 協力研究員, 理研, Pai-chi Li 協力研究員, 理研
高瀬規男, 磯子ソフト
HPCIヘルプデスク
松田元彦 研究員, 理研AICS, 丸山直也 チームリーダー, 理研AICS
本研究開発は以下の支援を受けて行われたものである。

•

文部科学省 最先端・高性能汎用スーパーコンピューターの開発利用「次世代生命体統合シミュレーションソフトウェアの研究開発」(RIKEN SCRP ISLiM)

•
•
•
•
•

理研

生命システム研究センター

理研

計算科学研究機構

JSPS科研費
理研

京コンピュータ

若手研究（B）(2012-2014) Number 24700299.

計算機センター RIKEN Integrated Cluster of Clusters (RICC)

JSPS科研費

基盤研究（C）(2014-2016) 開発中
84