2016年度 ポスト「京」重点課題1

GENESIS講習会（実習）
理化学研究所 生命システム研究センター
分子機能シミュレーション研究チーム
李秀栄、尾嶋拓、新津藍、杉田有治（資料提供：小林千草）

2016年9月8日 FOCUS（神戸）今日の内容n 13:30 – 13:50 講義• 生体分子の分子動力学計算• GENESISの特徴

n 休憩（10分）
n 14:00 – 16:00 実習•••

GENESISのコンパイル
実行例(1) 分子動力学法計算と解析
実行例(2) レプリカ交換分子動力学計算と解析

今日の内容（改定）n 13:30 – 13:45 講義• 生体分子の分子動力学計算• GENESISの特徴

n 13:45 – 14:45 実習（１）••

GENESISのコンパイル
分子動力学法計算と解析

n 休憩（15分）
n 15:00 – 16:00 実習（２）•レプリカ交換分子動力学計算と解析

端末操作と利用ソフト
n FOCUSスパコンへのログイン•••

Xwin Serverを起動
Cygwin64 Terminalを起動
Terminalで「ssh –Y ユーザ名@ﬀ」と入力

n データ転送•

WinSCPを使用

n トラジェクトリーデータの可視化•

VMD1.9.2を使用（※）端末からログアウトすると、ユーザにより生成された全てのデータは消去されます。

実習材料
n 実習で使用する材料は以下のディレクトリにあります../share/Lecture_0908/
GENESIS_binary: (GENESISのバイナリファイル)• spdyn (GENESIS, spdyn)• prst_setup (並列I/Oのセットアップ)• crd_convert (トラジェクトリコンバータ)
Compile: (GENESISのコンパイル)
genesis-1.1.0 (最新のGENESIS code)

Tutorial_1: (Tutorial 1 : タンパク質のMD)• 1_setup• 2_minimiza^on• 3_hea^ng• 4_equilibra^on• 5_produc^on• 6_analysis

Test:（テスト計算用のファイル）• Regression_test/• test_common• test_atdyn• test_spdyn• test_remd …

Tutorial_2: (Tutorial 3: タンパク質のレプリカ交換MD)• 1_setup• 2_minimiza^on• 3_equilibra^on• 4_produc^on• 5_analysis

GENESISのコンパイル (1)n 推奨されるコンパイラ
インテル、富士通コンパイラ、または、gfortran 4.4.7より新しいもの※Analysis toolのコンパイルにはLAPACKが必要になります

GENESIS-1.1.0のソースツリー.:
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)n ../share/Lecture_0908/Compileを自分のディレクトリにコピー
し、GENESISのsrcディレクトリに移動してください% cp -r ../share/Lecture_0908/Compile .
% cd Compile/genesis-1.1.0/src

n コンパイルはmakeを使います% ./conﬁgure --host=FOCUS （もしくは、FCFLAGS=“-axAVX,SSE4.2 –O3 –ip –mkl=parallel” ./conﬁgure）% make installn 出来上がったバイナリはgenesis-1.1.0/binに作成されます••

GENESISはMD、トラジェクトリ変換プログラム、解析プログラム合わせて20個のアプリケーションを作成します
今回は、時間の都合上、既にコンパイル済みのアプリケーションを使います

実行ファイルの場所： ../share/Lecture_0908/GENESIS_binary

FOCUSスパコンシステムの利用について
n FOCUSスパコンシステム概要• フロントエンドサーバ（RHEL 6.6）、演算サーバ（OS：CentOS 6.6）
• 開発環境：Glibc、GNU CC、Intel Parallel Studio XE• ジョブ管理システム：SLURMベース

n 端末利用室に設置の各種端末はインターネットに接続されておりません。このため「FOCUSスパコン利用の手引き」は以下のサイト（FOCUSスパコンネットワーク内）より閲覧してください。
hmp://ﬂ01.j-focus.jp/ablog/usermanual/演算ノードとキュー構成 （hmp://www.j-focus.jp）

ジョブ管理コマンド （hmp://www.j-focus.jp/user_guide/
ug0000000000.html）

GENESISテストセットn GENESISは、正しくコンパイルされたか確認するために、開発者が実施した計算の入出力ファイルをテストセットとして提供します% cp -r ../share/Lecture_0908/Test .
% cd Test

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

test_remd.py
test_rpath.py
test_rpath_atdyn/
test_rpath_spdyn/
test_spdyn/

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

テストジョブの実施
n SLURM用バッチスクリプト（run.sh）を編集します#!/bin/bash#SBATCH -p c006m#SBATCH -N 1
#SBATCH --ntasks-per-node=8#SBATCH -c 1
#SBATCH -J run#SBATCH -e run.e%J#SBATCH -o run.o%J
module load PrgEnv-intel
module load intel/openmpi165# キュー名# ノード数# ノードあたりのプロセス数# スレッド数# ジョブ名# 標準エラー出力# 標準出力# 環境設定読み込み

export OMP_NUM_THREADS=${SLURM_CPUS_PER_TASK}
bindir=/home1/gleu/share/Lecture_0908/GENESIS_binary./test.py “mpirun -np 8 ${bindir}/spdyn”## 実行ファイルのディレクトリ名# テストジョブの実行

n ジョブを投入します% sbatch run.sh% sjstat# 投入ジョブの確認

出力ファイル（run.oXXX）を開き結果を確認します：Passed: 45/45?

Tutorial 1: タンパク質のMD (1)n 計算内容水中のBPTI(牛膵臓トリプシン阻害酵素）のMD計算
参照：hmp://www.aics.riken.jp/labs/cbrt/tutorial/
Basic molecular dynamics simula^ons

n 計算データ

BPTI% cp –r ../share/Lecture_0908/Tutorial_1 .
% cd Tutorial_1

Tutorial_1の中身.:
1_setup/2_minimization/
3_heating/
4_equilibration/
5_production/
6_analysis/# セットアップ
# 最適化
# 昇温
# 平衡化
# プロダクション# 解析システムのセットアプ（1_setup）
n GENESISは、現在インプット用の分子ファイルを作成する部分は存在せず、CHARMMやVMDで作成された分子ファイル(top/par)を読み込むことができます• 今回は時間の都合により、予め作成した分子ファイル（ionize.pdb、
ionize.psf）を用います• 分子ファイル作成の詳細は、HPをご覧くださいn 初期構造の確認• WinSCPを使って端末にionize.pdbとionize.psfをダウンロードしますØ ホスト名：「ﬀ」• VMDを起動し、ionize.psf、ionize.pdbの順で読み込みますØ 「File」→「New Molecule…」
Ø Filename：の「Browse…」で「ionize.psf」を選択して「load」
Ø Load ﬁles for：が「0: ionize.psf」になっている状態で、Filename：「Browse…」で「ionize.pdb」を選択して「load」

構造最適化（2_minimizaCon）
n 初期構造では、原子間の距離が近すぎたり、不自然な原子結合を持つようなことがあります。不安定なシミュレーションを避けるため、本計算を行う前に「平衡化」と呼ばれる予備
計算を十分に行います。
n 構造最適化は、その第１段階でエネルギー値が下がる方向に粒子を動かします。GENESISでは現在、最急降下法(Steepest Decent)を使うことができます./2_minimization:
run.inp
run.shoutput# GENESISの入力ファイル# FOCUS用のバッチスクリプト# 計算結果構造最適化の入力ファイル[INPUT]
topﬁle = ../1_setup/top_all27_prot_lipid.ry # topology ﬁleparﬁle = ../1_setup/par_all27_prot_lipid.prm # parameter ﬁle
psﬃle = ../1_setup/ionize.psf# protein structure ﬁle
pdbﬁle = ../1_setup/ionize.pdb# PDB ﬁlereﬃle = ../1_setup/ionize.pdb# reference for restraints[OUTPUT]
dcdﬁle = run.dcd # DCD trajectory ﬁlersyile = run.rst # restart ﬁleCoulomb相互作用でPMEを用いる[ENERGY]
electrostaCc = PME# [CUTOFF,PME]
switchdist= 10.0 # switch distancecutoﬀdist= 12.0 # cutoﬀ distancepairlistdist = 13.5 # pair-list distancepme_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]contact_check = YES

rstout_period = 1000

# 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]
nfunc^ons = 1 # number of func^ons
funcCon1 = POSI # restraint func^on type
位置拘束
constant1 = 10.0 # force constantselect_index1 = 1 # restrained group どの粒子に拘束をかけるのかは[SELECTION]で選ぶ（[SELECTION]に関しては次ページ参照）[MINIMIZE]
nsteps= 1000 # number of steps
eneout_period = 100 # energy output period
crdout_period = 100 # coordinates output period

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

ステップ数

原子間距離を確認し「Bad contact」や不自然な距離をレポートする (本計算には使用しない）

that is defined by mol_name below.
mole_nameN molecule starting-residue ending-residue(‘N’ must be replaced by a positive integer」
value.) Define molecule by starting-residue and
GENESISの「selecCon
ending-residue. Those residues are defined by[segment id]:residue number:residue name
Table. Available keywords and operators in group.
expression
an:nameai:number[-[number]]atno:number[-[number]]rnam:namerno:number[-[number]]mname:namesegid:IDhydrogen
heavy
all
and
or
not
()

meaning
atom nameatom index*1
atom number*1
residue nameresidue number
molecule namesegment index
hydrogen atoms
heavy atoms
all atoms
conjunction
logical add
negationassemble

examplean:CA
ai:1-5
atno:6rnam:GLY
rno:1-5mname:molAsegid:PROA

other available expression
atomname, atom_nameatomindex, atomidx
atomnoresiduename, resname
residueno, resno
moleculename, molname
segmentid, sidhydrogenatom
heavyatom*&|!

（GENESIS User Guide 1.1.0）42構造最適化の実行
n ジョブの投入% cd 2_minimization% sbatch run.sh% sjstatsqueuesでキューのノード実行状況を調べて、
run.shを適宜編集してください# 投入ジョブの確認

n ジョブ終了後、以下のファイルが出力されているはずですrun.dcd
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*******************************************************[STEP0] Architecture and Compiler Information（※）
計算しなおす場合は、
dcdとrstファイルを消して下さい出力ファイルの見方
n 出力ファイルは6段階に分かれている[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 -101787.3074
30.2325 11248.6855
2939.6563
74.9556260.8343
62.6065-72.0093 11784.2874-128086.3236
0.0000
INFO:
100 -114610.3616-130925.2074
1.03165.42093975.9892

2354.4594

51.4906

256.4631

22.6729-74.3336

9727.0726[STEP 6] : 終端処理と計測時間
Output_Time> Averaged ^mer proﬁle (Min, Max)
total ^me = 95.076
setup
=
2.334
dynamics = 92.742
energy = 86.844
integrator =0.319
pairlist =5.858 (
5.740,
5.991)
energybond=
0.271 (
0.261,
0.293)
angle=
0.276 (
0.252,
0.300)
dihedral =0.086 (
0.057,
0.119)

nonbond= 83.908 ( 83.701, 84.087)
pme real = 72.354 ( 70.931, 73.537)
pme recip = 11.545 ( 10.218, 13.149)
restraint =0.011 (
0.005,
0.018)
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)

計算結果の確認
n ポテンシャルエネルギーの変化を確認します• GENESIS出力ファイルから各ステップでのエネルギー値を抽出し、gnuplotでプロットします% sh pot.sh% gnuplot
gnuplot> plot “pot.dat” w lp

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

昇温計算 (3_heaCng/)n 構造最適化で、ある程度の安定構造にした後、システムの温度をゆっくり変化させ、本計算で用いる温度まで上昇させる。
n GENESISでは’Annealing’ オプションを使ってゆっくり温度を変化させるようなシミュレーションができます./3_heating:run.inp
run.shoutput# GENESISの入力ファイル# FOCUS用のバッチスクリプト# 計算結果昇温計算の入力ファイル抜粋[INPUT]
topﬁle = ../1_setup/top_all27_prot_lipid.ry # topology ﬁle(skip)
rsZile = ../2_minimizaCon/run.rst# restart ﬁle[DYNAMICS]
LEAP: Leapfrog
integrator = LEAP # [LEAP,VVER]
nsteps= 5000 # number of MD steps VVER: Velocity Verlet^mestep= 0.002 # ^mestep (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# ini^al temperature (K)チュートリアルでは時間の関係でこのステップ数にしていますが、実際はもっと長い計算を行います（※）実際の計算では、計算システムに応じて、標的に温度になるようゆっくり上昇させる

昇温計算の実行
n ジョブの投入% cd 3_heating% sbatch run.sh% sjstat# 投入ジョブの確認

n ジョブ終了後、以下のファイルが出力されているはずですrun.dcd
run.out
run.rst# dcd形式のトラジェクトリファイル (binary)# GENESIS出力ファイル (ascii)# GENESISリスタートファイル (binary)

例：出力ファイルの[step 5]部分
INFO:
STEP
TIME
TOTAL_ENE POTENTIAL_ENE KINETIC_ENE
RMSG
BOND
ANGLE
DIHEDRALIMPROPERVDWAALS
ELECT UREY-BRADLEYCMAPRESTRAINT
BOXX
BOXY
BOXZ
VOLUME TEMPERATUREPRESSXXPRESSYYPRESSZZVIRIALPRESSURE(..skip..)
Simulated_Annealing_Leapfrog> Anneal temperature from 210.100 to 213.100
INFO:
3600
7.2000 -121579.0827 -135018.8050 13439.7223
14.8579
126.2918
326.4705
296.5269
18.2507 17470.6626 -153272.5717
43.5707-75.9996
47.993170.8250
83.2579
69.0930 407423.5098
181.9966 -1507.9456 -1629.7733 -1447.1411 -18040.6915 -1528.2867
Simulated_Annealing_Leapfrog> Anneal temperature from 213.100 to 216.100
INFO:
3650
7.3000 -121202.4067 -134760.3492 13557.9425
14.8412133.4514
322.5176
297.5756
18.2225 17256.3516 -152798.7766
47.1838-80.8475
43.972570.8250
83.2579
69.0930 407423.5098
183.5975 -1553.3990 -1690.9955 -1543.6580 -18521.9527 -1596.0175

計算結果の確認
n 温度の変化を確認します• 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

平衡化(4_equilibraCon/)n 本計算と同じ条件の計算を行い、システムを平衡化させる(研究を行う上では、非常に重要なステップ)./4_equilibration:
run.inp
run.shoutput

入力ファイル抜粋[OUTPUT]
dcdﬁle = run.dcd # DCD trajectory ﬁle# GENESISの入力ファイル# FOCUS用のバッチスクリプト# 計算結果トラジェクトリ書き出しを行い、温度・エネルギーだけでなく、中の分子も平衡に達していることも確認する[ENERGY][CONSTRAINTS][RESTRAINS (今回はない）][ENSEMBLE]
これらのパラメータは基本的に本計算と同じものを使う[ENSEMBLE]
ensemble = NPT# [NVE,NVT,NPT]
tpcontrol = LANGEVIN # thermostat and barostat
temperature = 300.0 # ini^al temperature (K)pressure = 1.0# target pressure (atm)

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

計算結果の確認
n 体積の変化を確認します• 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本計算(5_producCon/)n 実際にデータを取るための計算./5_production:
run.inp
run.shoutput

入力ファイル抜粋# GENESISの入力ファイル# FOCUS用のバッチスクリプト# 計算結果トラジェクトリ書き出しを行うことで、研究に用いる様々な性質を計算する[OUTPUT]
dcdﬁle = run.dcd # DCD trajectory ﬁle

チュートリアルでは時間の関係でこのステップ数にしていますが、実際はもっと長い計算を行います。

解析(6_analysis/)n 求めたトラジェクトリから様々な解析を行います• 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の入力ファイル# FOCUS用のバッチスクリプト

解析の入力ファイル
n Trajectory、selec^on、ﬁ•ngで解析対象と内容を指定します
参照：hmp://www.aics.riken.jp/labs/cbrt/tutorial/
Trajectory analysis tools[INPUT]reﬃle = ../../1_setup/ionize.pdb # PDB ﬁle[OUTPUT]
rmsﬁle= run.rms# RMSD ﬁletrj_natom=0# (0:uses reference PDB atom count)

Caのみ選択[SELECTION]
group1= an:CA # selec^on group 1[FITTING][TRAJECTORY]ﬁ•ng_method = TR+ROT # method
trjﬁle1= ../../5_produc^on/output/run.dcd # trajectory ﬁ•ng_atom = 1# atom groupﬁlemass_weight = YES# mass-weight is applied
md_step1= 500000 # number of MD steps
mdout_period1 = 500# MD output period[OPTION]
基準構造に対する重ana_period1 = 500# analysis period
check_only = NO
# (YES/NO)
ね合わせを行う
repeat1=1
analysis_atom = 1# atom group(TR+ROT: 並進と回転
trj_format = DCD# (PDB/DCD)
複数trjﬁleを扱う場
で重ね合わせする)
trj_type= COOR# (COOR/COOR+BOX)合に指定簡略化

YES：本計算しない

n run.shの実行部分を以下の様に変更して実行します${bindir}/rmsd_analysis run.inp > run.out

ドメインとセルの決め方
ドメイン：

、セル：
• MPIプロセス数に合わせてドメインに分割：

ドメイン数＝プロセス数• セル幅 d は、
２つ隣まで
考慮

d > (pairlistdist + buffer)/2を満たしつつ、セル数がドメイン数で割り切れるようにセル幅を決定する。
• ２つ隣のセルまで相互作用を考慮する
ため、周期境界条件の場合、最低でも5 x 5 x 5 のセルが必要。

Tutorial 2: タンパク質のレプリカ交換MD (1)n 計算内容水中のTrialanineのレプリカ交換MD計算
参照：hmp://www.aics.riken.jp/labs/cbrt/tutorial/
Replica-exchange molecular dynamics simula^ons

Trialaninen 計算データ% cp –r ../share/Lecture_0908/Tutorial_2 .
% cd Tutorial_2

Tutorial_2の中身.:
1_setup/2_minimize_pre-equil/
3_equilibration/
4_production/
5_analysis/# セットアップ# 最適化と初期平衡化# 平衡化# プロダクション# 解析レプリカ交換MD法
n エネルギー極小状態に捕われず幅広い配位空間を探索
温度の異なる複数のレプリカ（コピー）を用意し独立にシミュレーションを実行、定期的にレプリカの温度を交換することでサンプリング効率を格段に向上(Y. Sugita and Y. Okamoto. 1999)
Replica 4

Replica 3

Replica 2

Replica 1

MD

MD

Temperatureexchange

MD

MD

Temperatureexchange
Rejectedレプリカ交換MD計算の手順n 実際の研究においてレプリカ数と温度の幅はシステムの原
子数、システム自体の性質によって変わります
レプリカ数や温度幅の推定に役立つサイト
Temperature generator for REMD-simula^ons：hmp://folding.bmc.uu.se/remd/
Patriksson and van der Spoel, Phys. Chem. Chem. Phys., 2008, 10:2073-2077

レプリカ交換MD計算スキーム（実習では、Replica数は20、温度範
囲は300K〜351K）

Replica 20
Replica 4
Replica 3
Replica 2

351K
…
303K

…
…
…
exchanged exchanged

302K
301K
300K

Setup
step1

Minimize & pre-equilibrate
step2

Equilibrate
step3

rejected
Produc^on
step4

各ステップの説明•

Simula^on systemの作成 (1_setup/)•

構造最適化と平衡化（2_minimize_pre-equil/）
–
–
–
–

•

Step2.1：構造最適化
Step2.2：昇温計算（NVT、T=300K）（拘束あり）
Step2.3：平衡化（NPT、T=300K、P=1atm）（拘束あり）
Step2.4：平衡化（NPT、T=300K、P=1atm）

単独MD、
Tutorial_1と同様

各レプリカの平衡化(3_equilibrate/)– 各レプリカのシステムを該当温度で平衡化(レプリカの交換なし）•

各レプリカの平衡化と本計算(4_produc^on/)•

解析（5_analysis）

レプリカ交換MD平衡化計算の入力[INPUT]
topﬁle = ../1_setup/top_all36_prot.ry# topology ﬁleparﬁle = ../1_setup/par_all36_prot.prm# parameter ﬁlestrﬁle = ../1_setup/toppar_water_ions.str # stream ﬁle
psﬃle = ../1_setup/wbox.psf# protein structure ﬁle
pdbﬁle = ../1_setup/wbox.pdb# PDB ﬁlersyile = ../2_minimize_pre-equi/step2.4.rst # restart ﬁle[OUTPUT]
dcdﬁle = step3_rep{}.dcd # DCD trajectory ﬁlersyile = step3_rep{}.rst# restart ﬁleremﬁle = step3_rep{}.rem# replica exchange ID ﬁlelogﬁle = step3_rep{}.log# log ﬁle of each replica[ENERGY]
forceﬁeld= CHARMM # CHARMM force ﬁeld
electrosta^c = PME# Par^cle Mesh Ewald method
switchdist= 10.0 # switch distancecutoﬀdist= 12.0 # cutoﬀ distancepairlistdist = 13.5 # pair-list distance
vdw_force_switch = YES# force switch op^on for van der Waals
pme_nspline = 4# order of B-spline in [PME]pme_max_spacing = 1.0# max grid spacing[DYNAMICS]
integrator= LEAP # leapfrog Verlet integrator
nsteps= 1000 # number of MD steps^mestep= 0.002 # ^mestep (ps)eneout_period = 100 # energy output period
crdout_period = 100 # coordinates output period
rstout_period = 1000 # restart output period
nbupdate_period = 10 # nonbond update period[CONSTRAINTS]
rigid_bond= YES# constraints all bonds involving hydrogen[ENSEMBLE]
ensemble= NVT# NPT ensemble
tpcontrol= LANGEVIN # thermostat and barostat
temperature = 300.00 # ini^al and target temperature (K)[BOUNDARY]
type
= PBC

# periodic boundary condi^on[REMD]
dimension=1exchange_period = 0

レプリカ交換MDの指定セクション# number of parameter types# NO exchange for equilibra^on

レプリカの更新頻度

type1= temperature # T-REMD
nreplica1= 20# number of replicas (交換をしない)
parameters1 = 300.00 302.53 305.09 307.65 310.24 \
312.85 315.47 318.12 320.78 323.46 \326.16 328.87 331.61 334.37 337.14 \339.94 342.75 345.59 348.44 351.26

単独のリスタートファイルからスタート
レプリカ毎に出力ファイルを生成（※）{}内はレプリカを示す数字で展開される(e.g. step3_rep1.dcd)

レプリカ交換MD本計算の入力[INPUT]
topﬁle = ../1_setup/top_all36_prot.ry# topology ﬁleparﬁle = ../1_setup/par_all36_prot.prm# parameter ﬁlestrﬁle = ../1_setup/toppar_water_ions.str # stream ﬁle
psﬃle = ../1_setup/wbox.psf# protein structure ﬁle
pdbﬁle = ../1_setup/wbox.pdb# PDB ﬁlersyile = ../3_equilibrate/step3_rep{}.rst # restart ﬁle[OUTPUT]
dcdﬁle = step4_rep{}.dcd# DCD trajectory ﬁlersyile = step4_rep{}.rst# restart ﬁleremﬁle = step4_rep{}.rem# replica exchange ID ﬁlelogﬁle = step4_rep{}.log# log ﬁle of each replica[ENERGY]
forceﬁeld= CHARMM # CHARMM force ﬁeld
electrosta^c = PME# Par^cle Mesh Ewald method
switchdist= 10.0 # switch distancecutoﬀdist= 12.0 # cutoﬀ distancepairlistdist = 13.5 # pair-list distance
vdw_force_switch = YES# force switch op^on for van der Waals
pme_nspline = 4# order of B-spline in [PME]pme_max_spacing = 1.0# max grid spacing[DYNAMICS]
integrator= LEAP # leapfrog Verlet integrator
nsteps= 2500000 # number of MD steps^mestep= 0.002 # ^mestep (ps)eneout_period = 500 # energy output period
crdout_period = 500 # coordinates 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# NPT ensemble
tpcontrol= LANGEVIN # thermostat and barostat
temperature = 300.00 # ini^al and target temperature (K)[BOUNDARY]
type
= PBC

# periodic boundary condi^on

レプリカ交換MDの指定セクション
レプリカの次元数[REMD]
dimension=1# number of parameter types
exchange_period = 2500
# NO exchange for equilibra^on
レプリカの更新頻度

type1= temperature # T-REMD
nreplica1= 20# number of replicas
parameters1 = 300.00 302.53 305.09 307.65 310.24 \
312.85 315.47 318.12 320.78 323.46 \326.16 328.87 331.61 334.37 337.14 \339.94 342.75 345.59 348.44 351.26

レプリカ毎のリスタートファイルからスタート
レプリカの数だけ出力ファイルを生成（※）{}内はレプリカを示す数字で展開される(e.g. step3_rep.dcd)

レプリカ交換MD計算の実行（3_equilibrate/）n 今回は、時間の関係上、平衡化計算のみを実施していただ
きます% cd 3_equilibrate% sbatch step3.sh% sjstat# 投入ジョブの確認

n バッチスクリプト（step3.sh）の編集#!/bin/bash#SBATCH -p c006m#SBATCH -N 2
#SBATCH --ntasks-per-node=10#SBATCH -c 1
#SBATCH -J run#SBATCH -e run.e%J#SBATCH -o run.o%J
module load PrgEnv-intel
module load intel/openmpi165# キュー名# ノード数# ノードあたりのプロセス数# スレッド数
squeuesコマンドでキューの空き状況を確認し、適宜変更して下さい※今の指定は、1プロセス1レプリカの計算で、合計20コアで流れます（2ノード（10コア/ノード））

export OMP_NUM_THREADS=${SLURM_CPUS_PER_TASK}
bindir=/home1/giuo/share/Lecture_0908/GENESIS_binary
mpirun -np 20 ${bindir}/spdyn step3.inp > step3.log

レプリカ交換MD計算の出力ファイル
n 各レプリカ毎に出力ファイルが生成されます
•
•
•
•
•

step3.log
step3_rep[1-20].log
step3_rep[1-20].dcd
step3_rep[1-20].rst
step3_rep[1-20].rem：全体の出力ファイル：各レプリカのエネルギーファイル：各レプリカのdcd形式のトラジェクトリファイル (binary)：各レプリカのGENESISリスタートファイル (binary)：レプリカ情報ファイル
各レプリカが何番の温度を計算しているのかを出力
トラジェクトリの解析に重要！

レプリカ交換MD出力ファイルの抜粋
レプリカ交換出力ファイル

Step数

R: 交換失敗(Rejected)
A: 交換成功(Accepted)
REMD> Step:
Replica
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

2500
Dimension:
ExchangeTrial
1 >
0
N
2 >
3
A
3 >
2
A
4 >
5
A
5 >
4
A
6 >
7
A
7 >
6
A
8 >
9
R
9 >
8
R
10 >
11
R
11 >
10
R
12 >
13
A
13 >
12
A
14 >
15
A
15 >
14
A
16 >
17
A
17 >
16
A
18 >
19
R
19 >
18
R
20 >
0
N

Parameter:
323.460
326.160
RepIDtoParmID:
10
11
ParmIDtoRepID:
10
11

300.000
331.610113
1
13305.090328.870312
3
12

レプリカ情報ファイル

1

交換前後の温度

ExchangePattern:
AcceptanceRatio
0 /
0
1 /
1
1 /
1
1 /
1
1 /
1
1 /
1
1 /
1
0 /
1
0 /
1
0 /
1
0 /
1
1 /
1
1 /
1
1 /
1
1 /
1
1 /
1
1 /
1
0 /
1
0 /
1
0 /
0

2
Before300.000
302.530305.090307.650
310.240
312.850315.470318.120320.780
323.460
326.160
328.870331.610334.370
337.140
339.940
342.750345.590348.440351.260

302.530
337.140215
2
15307.650342.750417
4
17310.240
334.370
5
14
5
14

After
300.000
305.090302.530310.240
307.650315.470
312.850
318.120320.780
323.460
326.160
331.610328.870
337.140
334.370
342.750
339.940
345.590348.440351.260
315.470339.940
7
16
7
16

0
2500
5000
7500
10000
12500
15000
17500
20000
22500
25000

1
1
1
1
1
1
2
2
1
1
2

温度の番号

312.850345.590618
6
18

318.120348.440819
8
19320.780351.260920
920解析(5_analysis/)n remd_convertを用いて、求めた各レプリカのトラジェクトリから温度毎のトラジェクトリファイルへと変換します% cd 5_analysis

5_analysisの中身.:
1_calc_ratio/
2_plot_index/
3_plot_potential/
4_sort_dcd/
5_PCA/# 交換効率の計算# 各温度におけるレプリカIDの時間変化# 各レプリカのポテンシャルエネルギーの時間変化# 温度毎のトラジェクトリファイルの生成# PCA解析解析の入力ファイル
n remd_convertを用いて、求めた各レプリカのトラジェクトリから温度毎のトラジェクトリファイルへと変換します[INPUT]
psﬃle = ../../1_setup/wbox.psf
reﬃle = ../../1_setup/wbox.pdb
dcdﬁle = ../../4_produc^on/step4_rep{}.dcd
remﬁle = ../../4_produc^on/step4_rep{}.rem[OUTPUT]
pdbﬁle = trialanine.pdb
trjﬁle = remd_paramID{}_trialanine.dcd[SELECTION]
group1 = atomno:1-42
group2 = resno:2 and (an:N or an:CA or an:C or an:O)[FITTING]ﬁ•ng_method = TR+ROT
ﬁ•ng_atom = 2
mass_weight = YES

レプリカの番号からパラメータ（今回は温度）の空間への変換[OPTION]
check_only = NOconvert_type = PARAMETER
convert_ids = 1# (empty = all)dcd_md_period = 500
trjout_format = DCD
trjout_type = COOR+BOX
温度1 (300K) のみの出力
trjout_atom = 1(デフォルトまたは空欄で
pbc_correct = NO

全てのパラメータを出力)（※）convert_ids = 1 3 5とすると1 (300K) , 5 (304K),
6 (305K)のファイルを出力