AICS公開ソフト講習会 第16回
「GENESIS」 実習
理研、計算科学研究機構(AICS)
粒子系生物物理研究チーム
2017/01/13

今日の内容• 13:00 – 17:00 実習
 GENESISのコンパイル 実行例(1) 分子動力学法計算
 GENESISのトラブルシューティング 実行例(2) レプリカ交換分子動力学法計算と自由エネルギー計算

実習の材料
今回の実習の材料は以下のディレクトリにあります
最初にこのファイルを自分のホームディレクトリで展開してください/data/acourse/Lecture_0113/tutorial_data.tar% tar xvf /data/acourse/Lecture_0113/tutorial_data.tar
Compile: (GENESISのコンパイル)
genesis‐1.1.1 (最新のGENESIS code)
tests‐1.1.1 (コンパイルテスト）
env: (今回のチュートリル用の設定ファイル）
Tutorial_1: (Tutorial 1 : タンパク質MD)• 1_setup• 2_minimization• 3_heating• 4_equilibration••

5_production
6_analysis

Tutorial_2: (Tutorial 3: タンパク質のレプリカ交換MD)• 1_setup• 2_minimize_pre‐equi• 3_equilibrate• 4_production• 5_analysis

今日のチュートリアルの内容は以下からもダウンロードできます
Tutorial_1
http://www.aics.riken.jp/labs/cbrt/tutorial/basic_md_tutorials/tutorial‐1‐1/
Tutorial_2
http://www.aics.riken.jp/labs/cbrt/tutorial/remd‐tutorials/tutorial‐2‐2/

GENESISのコンパイル (1)GENESISのsrcディレクトリに移動してください% cd Compile/genesis-1.1.1/src
GENESISのソースツリー.:
COPYINGREADME
src/
bin/

簡単な説明
ソースコードバイナリ（コンパイル後）./src:
GENESIS_VERSION
Makefile.am
Makefile.in
aclocal.m4
analysis/
解析ツールのコード
atdyn/
atdynのコード
bootstrap
cleanupconfig.h.in
configureconfigure.ac

depcompinstall-sh
lib/
共通コード
Missing
spdyn/
spdynのコード./src/analysis:
Makefile.am
Makefile.in
libana/analysis共通コード
pcaana/
PCA解析
rpath_generator/ rpath生成コード
rstcnv/
rst_convertのコード
trjana/trj_analysisのコード
trjcnv/
crd/pcrd/remd_convertのコード

GENESISのコンパイル (2)推奨されるコンパイラ
Intel , Fujitsu compiler
gfortran 4.4.7より新しいもの

GENESISのコンパイルにはlapack, blasが必要です

1. コンパイルはautoconfとmakeを使います% ./configure --host=fx10% make% make install

2. 出来上がったバイナリはgenesis‐1.1.1/binに作成されます

GENESISはMD, トラジェクトリ変換プログラム、解析プログラム合わせて9つのア
プリケーションを作成します
今回は、時間の都合で、すでにコンパイル済みのアプリケーションを使いますGENESISのコンパイル (3)その他の環境でのコンパイルの仕方
0. configureはヘルプで使用可能なオプション一覧を表示できます% ./configure --help
1. (クロスコンパイルではない）通常のコンパイル% ./configure% make install
2. Mixed precision(混合精度）での計算の仕方% ./configure --enable-single% make install
3. GPU&混合精度での計算の仕方% ./configure --enable-gpu --enable-single% make install
4. デバッグモードでのコンパイル% ./configure --enable-debug=3% make install

GENESISのテスト
GENESISが正しくコンパイルされている事を確認するために、開発者が提供する入出力ファイルによるテストを行うことを推奨します
こちらも、ソースコードと同様にGENESISのウェブサイトからダウンロードが可能です% cd Compile/tests-1.1.1/regression_test

GENESIS‐1.1.0テストセットのツリー./regression_test:
test_remd.py
test_rpath.pybuild/
計算用の初期ファイル
charmm.py
test_rpath_atdyn/
genesis.py
test_rpath_spdyn/
param/
パラメータファイル
test_spdyn/
run.shバッチジョブ用のスクリプト
test.py
test_atdyn/
atdynのテスト
test_common/
共通テスト
test_nonstrict.py
test_parallel_IO/ パラレルIOのテスト
test_remd/
レプリカ交換のテスト
test_remd.csh

rpath(atdyn)のテスト
rpath(spdyn)のテスト
spdynのテスト

GENESISのテスト (2)コンパイルテストの方法
FX10環境では、今回こちらで用意したrun.shというスクリプトを使います（省略）通常のワークステーションでのやり方は以下です% cd regression_tests% export OMP_NUM_THREADS=1% ./test.py “mpirun -np 8 ${PATH_GENESIS}/atdyn”% ./test.py “mpirun -np 8 ${PATH_GENESIS}/spdyn”

プロセス数は必ず８にしてください。スレッド数はプロセス数以下が速度上、推奨です
GPUを用いた際(‐‐enable‐gpu)のテスト方法% export OMP_NUM_THREADS=1% ./test.py “mpirun -np 8 ${PATH_GENESIS}/spdyn” gpu

REMD, R‐PATHのテスト方法% export OMP_NUM_THREADS=1% ./test_remd.py “mpirun -np 8 ${PATH_GENESIS}/spdyn”% ./test_remd.py “mpirun -np 8 ${PATH_GENESIS}/atdyn”% ./test_rpath.py “mpirun -np 8 ${PATH_GENESIS}/spdyn”% ./test_rpath.py “mpirun -np 8 ${PATH_GENESIS}/atdyn”

注意：倍精度・混合精度は自動で判別されますTutorial 1: タンパク質のMD
GENESIS マニュアルのTutorial 1 Tutorial_1/
水中のBPTI(牛膵臓トリプシン阻害酵素）のMD計算

http://www.aics.riken.jp/labs/cbrt/tutorial/basic_md_tutorials/tutorial‐1‐1/

1. Simulation systemの作成 (1_setup/)
GENESISは、現在インプット用の分子ファイルを作成する部分は存在せず、
CHARMMやVMDで作成された分子ファイル(top/par)を読み込むことができます
今回は時間の都合により、予め作成した分子ファイルを用います（分子ファイルの作成方法に関しては、詳しくはTutorialをご覧ください）Tutorial 1: 構造最適化(2_minimization/)
初期構造では、原子間の距離が近すぎたり、不自然な原子結合を持つようなことがあります。不安定なシミュレーションを避けるため、本計算を行う前に“平衡化”と呼ばれる予備計算を十分に行います。

構造最適化は、その第１段階でエネルギー値が下がる方向に粒子を動かします
GENESISでは現在、最急降下法(Steepest Decent)を使うことができます
Tutorial_1/2_minimization/
run.inp : GENESIS入力ファイル
run.sh : fx10用のジョブサブミットスクリプト
output/ : 計算結果

構造最適化の入力ファイル[INPUT]
topfile = ../1_setup/top_all27_prot_lipid.rtf # topology fileparfile = ../1_setup/par_all27_prot_lipid.prm # parameter filepsffile = ../1_setup/ionize.psf# protein structure file
pdbfile = ../1_setup/ionize.pdb# PDB filereffile = ../1_setup/ionize.pdb# reference for restraints[OUTPUT]
dcdfile = run.dcd # DCD trajectory filerstfile = run.rst # restart file

Coulomb相互作用でPMEを用いる[ENERGY]
electrostatic = PME# [CUTOFF,PME]
switchdist= 10.0 # switch distancecutoffdist= 12.0 # cutoff distancepairlistdist = 13.5 # pair‐list distancecontact_check = YES
pme_ngrid_x = 72
# grid size_x in [PME]pme_ngrid_y = 80
# grid size_y in [PME]pme_ngrid_z = 72
# grid size_z in [PME]# restart output period[BOUNDARY]
type= PBC # [PBC,NOBC]
box_size_x = 70.8250 # box size (x) in [PBC]
box_size_y = 83.2579 # box size (y) in [PBC]
box_size_z = 69.0930 # box size (z) in [PBC][SELECTION]
group1= an:CA | an:C | an:O | an:N # index of restraint group 1[RESTRAINTS]
nfunctions = 1 # number of functions
function1 = POSI # restraint function type
位置拘束
constant1 = 10.0 # force constantselect_index1 = 1 # restrained group どの粒子に拘束をか

原子間距離を確認し「Bad contact」や不自然な距離をレポートし、さらに力の値に制限をかける（最適化・平衡化のみで使用してください）[MINIMIZE]
nsteps= 1000 # number of steps
eneout_period = 100 # energy output period
crdout_period = 100 # coordinates output period

ステップ数rstout_period = 1000

けるのかは[SELECTION]で選ぶ（[SELECTION]に関してはManual 11、または巻末を参照）

長距離相互作用で用いるFFTのグリッド数(2,3,5の倍数かつ、並列数によって最適値が変わる、
一般的にはグリッドサイズが1Å程度になるようにする。詳しくはマニュアルの5.3.1を参照)

構造最適化の実行
計算の実行% cd ~/Lecture_0113/Tutorial_1/2_minimization% pjsub run.sh

Jobの状況を調べるには% pjstat

Jobが終わったら以下のファイルが出力されているはずですrun.dcd :err.o :
run.out :run.rst :dcd形式のトラジェクトリファイル (binary)スクリプトのエラーファイル
GENESIS出力ファイル (ascii)
GENESISリスタートファイル (binary)

出力ファイル(run.out)*******************************************************
GENESIS SPDYN****
A Molecular Dynamics Simulator with**
Spatial Decomposition Scheme****
Developed by RIKEN AICS*******************************************************

出力ファイルの見方
GENESISの出力ファイルは7段階([STEP 0]‐[STEP 6])に分かれている[STEP 0] : 計算環境の確認[STEP 1] : 入力パラメータの確認[STEP 2] : 並列数（プロセス数、スレッド数）の確認[STEP 3] : 分子・エネルギー関数情報の確認[STEP 4] : 初期座標のエネルギー計算結果[STEP 5] : シミュレーション計算結果
INFO:
STEP POTENTIAL_ENE
RMSG
BOND
ANGLE UREY‐BRADLEYDIHEDRALIMPROPERCMAP
VDWAALS
ELECT RESTRAINT_TOTAL‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐INFO:
0 ‐101597.1585
30.2479 11247.1385
2939.9532
74.9561260.8373
62.6065‐72.0093 11798.1692 ‐127908.8099
0.0000
INFO:

100 ‐114456.2091

5.45973977.1511

2353.8183

51.5709

256.4578

22.6346‐74.3462[STEP 6] : 終端処理と計測時間
Output_Time> Averaged timer profile (Min, Max)
total time = 137.159
setup
=
5.843
dynamics = 131.316
energy = 122.896
integrator =0.838
pairlist =5.774 (
5.629,
5.963)
energybond=
0.657 (
0.587,
0.708)
angle=
0.717 (
0.696,
0.751)
dihedral =1.275 (
1.236,
1.295)

nonbond= 109.986 ( 109.658, 110.141)
pme real = 83.219 ( 81.500, 84.919)
pme recip = 26.077 ( 24.516, 27.847)
restraint =0.221 (
0.178,
0.251)
integrator
constraint =0.000 (
0.000,
0.000)
update=0.000 (
0.000,
0.000)
comm_coord =0.000 (
0.000,
0.000)
comm_force =0.000 (
0.000,
0.000)
comm_migrate =0.000 (
0.000,
0.000)

9730.9841 ‐130775.5114

1.0317

Tutorial 1:温度上昇シミュレーション (3_heating/)
構造最適化で、ある程度の安定構造にした後、システムの温度をゆっくり変化させ、本計算で用いる温度まで上昇させる。
GENESISでは’Annealing’ オプションを使ってゆっくり温度を変化させるようなシミュレー
ションができます
Tutorial_1/3_heating/run.inp : GENESIS入力ファイル
run.sh : fx10用のジョブサブミットスクリプト
output/ : 計算結果

入力ファイル抜粋[INPUT]
topfile = ../1_setup/top_all27_prot_lipid.rtf # topology file(skip)
rstfile = ../2_minimization/run.rst# restart file
[DYNAMICS]
LEAP: Leapfrog
integrator = LEAP # [LEAP,VVER]
nsteps= 5000 # number of MD steps VVER: Velocity Verlet
timestep= 0.002 # timestep (ps)(skip)
50 stepに1回、3Kずつ
annealing = YES# simulated annealing
温度を上げるanneal_period = 50
# annealing period
dtemperature = 3# temperature change at annealing (K)[CONSTRAINTS]
rigid_bond = YES# constraints all bonds# involving hydrogen[ENSEMBLE]
ensemble= NVT# [NVE,NVT,NPT]
tpcontrol = LANGEVIN # thermostat
temperature = 0.1# initial temperature (K)チュートリアルでは時間の関係でこのステップ数にしていますが、実際はもっと長い計算を行います温度上昇シミュレーションの実行
計算の実行% cd ~/Tutorial_1/3_heating% pjsub run.sh

計算が終わったら以下のファイルが出力されているはずですrun.dcd : dcd形式のトラジェクトリファイル (binary)
run.err : スクリプトのエラーファイル
run.out : GENESIS出力ファイル (ascii)
run.rst : GENESISリスタートファイル (binary)

出力ファイル抜粋[STEP 5]
INFO:
STEP
TIME
TOTAL_ENE POTENTIAL_ENE KINETIC_ENE
RMSG
BOND
ANGLE UREY‐
BRADLEYDIHEDRALIMPROPERCMAP
VDWAALS
ELECT RESTRAINT_TOTAL TEMPERATURE
VOLUME‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ ‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐(..skip..)
Simulated_Annealing_Leapfrog> Anneal temperature from 213.100 to 216.100
INFO:
3650
7.3000 ‐121015.8486 ‐134668.2442 13652.3956
14.8398
129.1861
316.5600
34.9546301.9510
20.3742‐76.0000 17444.0571 ‐152886.2023
46.8752
184.8766 407423.5098
Simulated_Annealing_Leapfrog> Anneal temperature from 216.100 to 219.100
INFO:
3700
7.4000 ‐120524.6696 ‐134365.7978 13841.1281
14.8065
135.7281
319.1936
42.6740291.7295
14.9113‐82.0820 17460.5756 ‐152587.8837
39.3559187.4323 407423.5098

計算結果の確認
 温度の変化を確認します• GENESIS出力ファイルから各ステップでの温度を抽出し、gnuplotでプ
ロットします (今回は省略)% sh temp.sh% gnuplot
gnuplot> plot "temp.dat" w lp

実行結果
$ cat temp.sh#!/bin/bash
grep "^INFO" output/run.out | tail -n+2 | awk '{print $3, $17}' >temp.dat

Tutorial 1:平衡化(4_equilibration/)
システムを計算させたい温度に上げた後、本計算と同じ条件の計算を行い、システムを平衡化させる (研究を行う上では、非常に重要なステップ)
Tutorial_1/4_equilibration/
run.inp : GENESIS入力ファイル
run.sh : fx10用のジョブサブミットスクリプト
output/ : 計算結果

入力ファイル抜粋[OUTPUT]
dcdfile = run.dcd # DCD trajectory file

トラジェクトリ書き出しを行い、温度・エネルギーだけでなく、中の分子も平衡に達していることも確認する[ENERGY][CONSTRAINTS][RESTRAINS (今回はない）][ENSEMBLE]
これらのパラメータは基本的に本計算と同じものを使う(3)と同様にrun.shを実行してください[ENSEMBLE]
ensemble = NPT# [NVE,NVT,NPT]
tpcontrol = LANGEVIN # thermostat and barostat
temperature = 300.0 # initial temperature (K)pressure = 1.0# target pressure (atm)

チュートリアルでは時間の関係でこのステップ数にしていますが、実際はもっと長い計算を行います（温度・エネルギーの安定性だけでなく、分子の観測したい性質が安定することも確認します）

計算結果の確認
 体積の変化を確認します• GENESIS出力ファイルから各ステップでの体積を抽出し、gnuplotでプ
ロットします(今回は省略)$ sh vol.sh$ gnuplot
gnuplot> plot "vol.dat" w lp

実行結果
$ cat vol.sh#!/bin/bash
grep "^INFO" output/run.out | tail -n+2 | awk '{print $3, $17}' > vol.dat

Tutorial 1:本計算(5_production/)
実際にデータを取るための計算
Tutorial_1/5_production/
run.inp : GENESIS入力ファイル
run.sh : fx10用のジョブサブミットスクリプト
output/ : 計算結果

入力ファイル抜粋

トラジェクトリ書き出しを行うことで、研究に用いる様々な性質を計算する[OUTPUT]
dcdfile = run.dcd # DCD trajectory file

チュートリアルでは時間の関係でこのステップ数にしていますが、実際はもっと長い計算を行います。（講義の資料の4P参照）(3)と同様にrun.shを実行してくださいTutorial 1:解析(6_analysis/) 求めたトラジェクトリから様々な解析を行います• GENESISは、得られたトジェクトリーを必要に応じて加工し、基本的な解析を行う独自のツールを提供します（例）./6_analysis:
1_RMSD
2_DIST
3_RMSF
4_PCA# RMSDの計算# 距離の計算# RMSFの計算# PCA解析• 今回は一例として、「trjana/rmsd_analysis」ツールを用いて、システムの中からCa原子のみ抜き出しそのRMSD(根二乗平均変位:基準の座標からどれだけずれたかを見る)を計算します./1_RMSD:
run.inp
run.sh# GENESISの入力ファイル# バッチスクリプト

解析の入力ファイル
 Trajectory、selection、fittingで解析対象と内容を指定します[INPUT]reffile = ../../1_setup/ionize.pdb # PDB file
[OUTPUT]
rmsfile= run.rms

trj_natom=0# (0:uses reference PDB atom count)

Caのみ選択[SELECTION]
group1= an:CA # selection group 1# RMSD file

[FITTING]
fitting_method = TR+ROT # method[TRAJECTORY]# atom group
trjfile1= ../../5_production/output/run.dcd # trajectory fitting_atom = 1
mass_weight = YES# mass‐weight is applied
file
md_step1= 500000 # number of MD steps[OPTION]
mdout_period1 = 500# MD output period
基準構造に対する重ana_period1 = 500# analysis period
check_only = NO
# (YES/NO)
ね合わせを行う
repeat1=1
analysis_atom = 1# atom group(TR+ROT: 並進と回転
trj_format = DCD# (PDB/DCD)
複数trjfileを扱う場
で重ね合わせする)
trj_type= COOR# (COOR/COOR+BOX)合に指定簡略化

YES：本計算しない run.shの内容は基本的にanalysisツールを実行するものですrmsd_analysis run.inp > run.out’trj_analysis’を用いる事で、距離・角度・二面角などを計算することも可能

GENESISのトラブルシューティング• GENESISが異常終了する
1. “contact_check” optionをYES, nstep=20で計算
MDステップに入る前の段階で構造のチェックを行います
ログに”too short distance”などのメッセージが表示された場
合は、構造が不安定ですので構造最適化・平衡化を行ってください
近日中にcontact_check機能を高めたversionの公開予定

2. ./configure ‐‐enable‐debug=3を用いて再コンパイル
セル内のアトム数などを確認しながら計算を行います
配列外アクセスが起きた際にはエラーを返しますので、セル数の最大値を変更するなどしてください

3. GENESISのForumに投稿してくださいTutorial 2: タンパク質のレプリカ交換MD
Tutorial_2/

アンブレラサンプリング
レプリカ交換MD計算

http://www.aics.riken.jp/labs/cbrt/tutorial/remd‐tutorials/tutorial‐2‐2/

レプリカ交換MD
 オリジナルのコピー（＝レプリカ）を多数作成し、温度、ハミルトニ
アン等の物理条件を変え、お互いを相互作用させることなくサンプリング
 特定のタイミングで隣接レプリカの物理条件を詳細釣り合わせに従い交換 物理条件の交換により、より効率の良いサンプリングを可能
 交換のタイミング以外では他のレプリカと相互作用しないため、超並列計算において高効率を達成

Relative free energy

Sugita and Okamoto, Chem. Phys. Lett., 314:141‐151 (1999)
Sugita et al., J. Chem. Phys., 113:6042‐6051 (2000)

レプリカ交換MDの手順
実際の研究においてレプリカ数と温度の幅はシステムの原子数、システム自体の性質によって変わります

Replica 14
Replica 4
Replica 3
Replica 2

13.76
…
4.56

…
…
…
exchanged
exchanged

3.64
2.72
1.80

Setup
step1

Minimization Equilibration
step2.1
Step2.2‐2.4

Equilibration
step3

rejected
Production
step4

今回は時間と使用可能なノード数の都合により、8レプリカ、かつ、ステップ数を短くして行います。

チュートリアルの内容
1. Simulation systemの作成 (1_setup/)
2. 構造最適化と平衡化 (2_minimize_pre‐equil/)
こちらはレプリカではなく単独MDとして走らせます

この手順はチュートリア
ル１と同様です
今回は時間の都合により、予め作成したデータを使います

3. 各レプリカの平衡化 (3_equilibrate/)
各レプリカ上のシステムを計算させたい距離に変化させ平衡化(レプリカの交換なし）

4. 本計算(4_production/)
レプリカ交換を行う

5. 解析(5_analysis/)
レプリカ交換を行う

Tutorial 2: 平衡化 (3_equilibrate/)[INPUT]
topfile = ../1_setup/top_all36_prot.rtf# topology fileparfile = ../1_setup/par_all36_prot.prm# parameter filestrfile = ../1_setup/toppar_water_ions.str # stream filepsffile = ../1_setup/wbox.psf# protein structure file
pdbfile = ../1_setup/wbox.pdb# PDB filerstfile = ../2_minimize_pre‐equi/step2.4.rst # restart file
[OUTPUT]
logfile = step3_rep{}.log
dcdfile = step3_rep{}.dcd
rstfile = step3_rep{}.rst[REMD]
dimension=1exchange_period = 0
iseed= 3141592
type1= RESTRAINT#nreplica1= 14
nreplica1=8cyclic_params1 = NO
rest_function1 = 1# log file of each replica# DCD trajectory file# restart file

レプリカ交換MDの指定セクション
レプリカの更新頻度(交換をしない)

交換されるrestraint関数の指定[ENERGY]
forcefield= CHARMM # CHARMM
force field（今日は時間と利用ノード数
electrostatic = PME# Particle Mesh Ewald methodの関係で8に変えています）
switchdist= 10.0 # switch distancecutoffdist= 12.0 # cutoff distancepairlistdist = 13.5 # pair‐list distance
vdw_force_switch = YES# force switch option for van der Waals
pme_nspline = 4# order of B‐spline in [PME]pme_max_spacing = 1.0# max grid spacing[DYNAMICS]
integrator = LEAP # [LEAP,VVER]#nsteps= 50000 # number of MD steps in REMD

nsteps= 5000 # number of MD steps in REMD
timestep= 0.002 # timestep (ps)eneout_period = 500 # energy output period
crdout_period = 500 # coordinates output period#rstout_period = 50000 # restart output period
rstout_period = 5000 # restart output period
stoptr_period = 500 # remove translational and rotational motions period
nbupdate_period = 10 # nonbond update period

単独のリスタートファイルからスター[CONSTRAINTS]
ト = YESrigid_bond# constraints all bonds involving hydrogen

レプリカ毎に出力ファイルを生成[ENSEMBLE]
ensemble = NVT# NVT ensemble（※）{}内はレプリカを示す数字で展
tpcontrol= LANGEVIN # Langevin thermostat
temperature= 300.0# target temperature (K)開される
gamma_t= 1.0# thermostat friction (ps‐1) in [LANGEVIN](e.g. step3_rep1.dcd)[BOUNDARY]
type
= PBC

# Periodic boundary condition[SELECTION]
group1= an:OY and resno:1 # restraint group 1
group2= an:HNT and resno:3 # restraint group 2[RESTRAINTS]
nfunctions = 1function1 = DIST
#constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2#reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 8.24 9.16 10.08 11.00 11.92
12.84 13.76
constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2
reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 8.24
select_index1 = 1 2

Tutorial 2: 本計算 (4_production/)[INPUT]
topfile = ../1_setup/top_all36_prot.rtf # topology fileparfile = ../1_setup/par_all36_prot.prm # parameter filestrfile = ../1_setup/toppar_water_ions.str # stream filepsffile = ../1_setup/wbox.psf# protein structure file
pdbfile = ../1_setup/wbox.pdb# PDB filerstfile = ../3_equilibrate/step3_rep{}.rst # restart file
[OUTPUT]
logfile = step4_rep{}.log
dcdfile = step4_rep{}.dcd
rstfile = step4_rep{}.rst
remfile = step4_rep{}.rem[REMD]
dimension=1exchange_period = 1000
type1
= RESTRAINT#nreplica1= 14
nreplica1=8cyclic_params1 = NO
rest_function1 = 1# log file of each replica# DCD trajectory file# restart file# parameter index file

レプリカ交換MDの指定セクション
レプリカの更新頻度

timestep= 0.002 # timestep (ps)#eneout_period = 50 # energy output period#crdout_period = 50 # coordinates output period
eneout_period = 500 # energy output period
crdout_period = 500 # coordinates output period#rstout_period = 2000000 # restart output period
rstout_period = 5000 # restart output period
nbupdate_period = 10 # nonbond update period[CONSTRAINTS]
レプリカ毎のリスタートファイルからスタート
rigid_bond = YES# constraints all bonds involving hydrogen
レプリカの数だけ出力ファイルを生成[ENSEMBLE]
ensemble = NVT# NVT ensemble
tpcontrol = LANGEVIN # Langevin thermostat
temperature = 300.0# target temperature (K)gamma_t= 1.0# thermostat friction (ps‐1) in [LANGEVIN][BOUNDARY]
type
= PBC

# Periodic boundary condition

交換されるrestraint関数の指定 [SELECTION]
group1= an:OY and resno:1 # restraint group 1[ENERGY]
forcefield= CHARMM # CHARMM force field
electrostatic = PME# Particle Mesh Ewald method
switchdist= 10.0 # switch distancecutoffdist= 12.0 # cutoff distancepairlistdist = 13.5 # pair‐list distance
vdw_force_switch = YES# force switch option for van der Waals
pme_nspline = 4# order of B‐spline in [PME]pme_max_spacing = 1.0# max grid spacing[DYNAMICS]
integrator = LEAP # [LEAP,VVER]#nsteps= 2000000 # number of MD steps in REMD
nsteps= 5000 # number of MD steps in REMD

group2= an:HNT and resno:3 # restraint group 2[RESTRAINTS]
nfunctions = 1function1 = DIST
#constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2#reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 8.24 9.16 10.08 11.00 11.92
12.84 13.76
constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2
reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 8.24
select_index1 = 1 2

レプリカ交換MDの実行
計算が終わったら以下のファイルが出力されているはずです
remd_run[1‐8].dcd
remd_run.out
remd_run[1‐8].log
remd_run[1‐8].dcd
remd_run[1‐8].rst
remd_run[1‐8].rem:各レプリカのdcd形式のトラジェクトリファイル (binary): 全体のアウトプットファイル: 各レプリカのエネルギーファイル: レプリカ交換出力ファイル (ascii): GENESISリスタートファイル (binary): レプリカ情報ファイル

各レプリカがどの距離（パラメータ)を計算しているのかを出力
トラジェクトリの解析に重要！
レプリカ交換出力ファイル

レプリカ情報ファイル

REMD> Step: 13000 Dimension: 1 ExchangePattern: 2
Replica ExchangeTrial
AcceptanceRatio Before
After1
1> 0 N0/01
1
2
3> 2 A2/
7
3
2
3
2> 3 A2/
7
2
3
4
7> 6 A3/
7
7
6
A: 交換成功(Accepted)
56> 7 A3/
7
6
7
R: 交換失敗(Rejected)
64> 5 R3/
7
4
4
7
5> 4 R3/
7
5
5
8
8> 0 N0/0
8
8Parameter :RepIDtoParmID:
ParmIDtoRepID:

1

2
1
1

3
2
2

6
3
3

7
6
6

4
7
7

5
4
4

8
5
5

8
8

0
1000
2000
3000

1
1
2
3

Step数

距離（パラ
メータ）の
番号

Tutorial 2: 解析 (5_analysis/)5_anaylsisではいくつかの解析のための例がありますが、今日はremd_convertを用
いて、求めた各レプリカのトラジェクトリからパラメータ毎のトラジェクトリファイルへと変換し、距離を求めた後にWHAMを用いて自由エネルギー面を求めます
1_convert_dcd
2_sort_dcd
3_calc_ratio
4_plot_index
5_calc_dist
6_calc_wham
今回は時間の都合により、予め作成したデータを使い、流れを説明します。
詳しくはチュートリアルのページをご覧ください
http://www.aics.riken.jp/labs/cbrt/tutorial/remd‐tutorials/tutorial‐2‐2/解析 (5_analysis/)
1. Coordinateファイルの変換 (1_convert_dcd) : crd_convert[INPUT]reffilepsffile= ../../1_setup/wbox.pdb= ../../1_setup/wbox.psf[OUTPUT]
pdbfile= step4_rep1.pdb
trjfile= step4_rep1.dcd[TRAJECTORY]#trjfile1= ../../4_production/step4_rep1.dcd#md_step1= 2000000# number of MD steps#mdout_period1 = 50
# MD output period#ana_period1 = 50
# analysis period
trjfile1= ../../4_production/saved_data/step4_rep1.dcd
md_step1= 40000# number of MD steps
mdout_period1 = 500# MD output period
ana_period1 = 500# analysis period
repeat1=1
trj_format = DCD# (PDB/DCD)

trj_type= COOR+BOX# (COOR/COOR+BOX)
trj_natom = 0# (0:uses reference PDB atom count)[SELECTION]
group1= an:N or an:CA or an:C or an:O
group2= sid:PROA[FITTING]
fitting_method = TR+ROT# methodfitting_atom = 2# atom groupmass_weight = YES# mass‐weight is applied

座標のfitting指定
セクション[OPTION]
check_only = NO
# (YES/NO)
trjout_format = DCD# (PDB/DCD)
trjout_type = COOR+BOX# (COOR/COOR+BOX)
trjout_atom = 2# atom grouppbc_correct = NO
# (NO/MOLECULE)

出力させる原子を選ぶ

2. Coordinateファイルを各パラメータ毎に配置し直す (2_sort_dcd) : remd_convert[FITTING][INPUT]
fitting_method = NO
# method
reffile= ../1_convert_dcd/step4_rep1.pdb # PDB filedcdfile= ../1_convert_dcd/step4_rep{}.dcd # DCD file
[OPTION]#remfile= ../../4_production/step4_rep{}.rem # REMD parameter ID file
配置の仕方を指定# (YES/NO)
remfile= ../../4_production/saved_data/step4_rep{}.rem # REMD parameter ID check_only = NOconvert_type = PARAMETER# (REPLICA/PARAMETER)
file
dcd/remファイルが必須
convert_ids =# (empty = all)(example: 1 2 5‐10)#dcd_md_period = 50
# input dcdfile MD period[OUTPUT]
dcd_md_period = 500# input dcdfile MD period
pdbfile= step4_par.pdb # PDB file
trjout_format = DCD# (PDB/DCD)
trjfile= step4_par{}.dcd # trajectory file
trjout_type = COOR+BOX# (COOR/COOR+BOX)
trjout_atom = 1[SELECTION]
pbc_correct = NO
# (NO/MOLECULE)
group1= all# selection group 1

解析 (5_analysis/)
5. Distanceの計算 (5_calc_dist) : trj_analysis[INPUT]reffile= ../../2_sort_dcd/step4_par.pdb[OUTPUT]
disfile= parameter1.dis# PDB file# distance file

[TRAJECTORY]
trjfile1 = ../../2_sort_dcd/step4_par1.dcd # trajectory file#md_step1 = 2000000# number of MD steps#mdout_period1 = 50
# MD output period#ana_period1 = 50
# analysis period
md_step1 = 40000# number of MD steps
mdout_period1 = 500# MD output period

ana_period1 = 500# analysis period
repeat1=1
trj_format = DCD# (PDB/DCD)
trj_type = COOR+BOX# (COOR/COOR+BOX)
trj_natom = 0# (0:uses reference PDB atom count)[OPTION]
check_only = NOdistance1 = PROA:1:ALA:OY PROA:3:ALA:HNT

それぞれのパラメータ毎に計算を行う
6. 自由エネルギーの計算 (6_calc_wham) : wham_analysis[INPUT]
cvfile= ../5_calc_dist/parameter/parameter{}.dis[OUTPUT]
pmffile= output.pmf[WHAM]
dimension = 1
nblocks=1
temperature = 300.0
rest_function1 = 1
tolerance = 1E‐08#grids1= 0.0 15.0 301# potential of mean force filegrids1= 0.0 8.0 201[RESTRAINTS]
function1 = DIST
#constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 ¥#1.2 1.2 1.2 1.2 1.2 1.2 1.2#reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 ¥#8.24 9.16 10.08 11.00 11.92 12.84 13.76
select_index1 = 1 2
constant1 = 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2
reference1 = 1.80 2.72 3.64 4.56 5.48 6.40 7.32 8.24
is_periodic1 = NOAppendix

GENESISの原子・分子選択”selection”（GENESIS User Guide 1.1.0）