GENESISでのMD計算ハンズオン
岩橋ー小林 千草
理化学研究所 計算科学研究センター
2021/09/28 GENESIS講習会 −Cygnusを用いたハンズオン−
1

内容• MD計算を高速に動かすために• REMD計算のデモ2

本ハンズオンの目的• 今回のハンズオンはCygnusでのGENESISによるMD計算を高速に動かすためのやり方を学びます (MDの本計算を想定)• MD/REMD計算の入力ファイルの作成方法や基本的なMD計算の手順は対象外です• ハンズオンで取り扱わなかった件に関しては、GENESISのウェブページをご覧いただくか、ご質問をフリーディスカッションで受け付けます
Usage

https://www.r-ccs.riken.jp/labs/cbrt/

Tutorial & Samples

3

ハンズオン資料の取得方法
1.
2.

ディレクトリからtutorial_GENESIS_2021Sep.tar.bz2ファイルをダウンロードし、cygnusのワーク領域へ置いてください
cygnusにログインしていただき、ワーク領域で展開します% cd /work/EDU2/<ユーザ名>% tar xvhf tutorial_GENESIS_2021Sep.tar.bz2% cd tutorial_GENESIS_2021Sep% ls
1_MD/ 2_REMD/

4

Hands-on 1

MD計算を高速に動かすために5

資料について今回の資料はGENESIS websiteのTutorial 3.3に基づいています
https://www.r-ccs.riken.jp/labs/cbrt/tutorials2019/tutorial-3-3/

MDの入力ファイルの作成方法、MD計算の準備（平衡化）などの手順は、上のウェブサイトに詳しく記載されています
水中のタンパク質の等温等圧(NPT)条件下でのMD計算

水溶液中のProtein G6

Step 1: 基本のMD
% cd
1_MD
% ls
1_VVER
3_RESPA_2.5fs
2_VVER_2.5fs 4_RESPA_2.5fs_tb10% cd 1_VVER% ls
INP run.sh sample

5_Benchmark
scripts

setup
topparvi等エディタでINPファイルを開いてください

7

Step 1: GENESIS control file
GENESISのインプットは[]で囲まれたセクションに分けられており、各セクション毎に決められたパラメータが存在します[INPUT] 入力ファイル[OUTPUT] 出力ファイル[ENERGY] 力場、エネルギー計算のパラメータ[MINIMIZE] or [DYNAMICS] エネルギー最小法かMD計算か (排他的)[REMD] or [GAMD] or [RPATH] or [ALCHEMY]各計算手法でのパラメータ (排他的)[CONSTRAINT] 分子内拘束[RESTRAINTS] restraint計算[SELECTION] restraintなどでの原子の選択[ENSEMBLE] 温度、圧力の情報[BOUNDARY] 周期境界条件の有無、Boxサイズなど
8

Step 1: Control file
[INPUT]
topfile = ../toppar/top_all36_prot.rtf
parfile = ../toppar/par_all36m_prot.prm
strfile = ../toppar/toppar_water_ions.str
psffile = ../setup/ionized.psf
pdbfile = ../setup/ionized.pdb
rstfile = ../setup/eq3.rst[ENERGY]
forcefield= CHARMM Coulomb相互作用でelectrostatic= PME
PMEを用いる
switchdist= 10.0cutoffdist= 12.0pairlistdist= 13.5
vdw_force_switch = YES
pme_nspline= 4
pme_max_spacing = 1.2[DYNAMICS]
integrator
nstepstimestepeneout_period
nbupdate_period=====[CONSTRAINTS]
rigid_bond= YES[ENSEMBLE]
ensemble
tpcontrol
temperaturepressure

= NPT
NPT計算
= BUSSI= 298.15= 1.0[BOUNDARY]
type= PBC

VVER
50000.002
10010Velocity Verlet法
ステップ数
ステップ幅（単位はps)

周期境界条件を利用
9Step 1: 基本のMDの動かし方
vi等エディタでrun.shファイルを開いてください#!/bin/bash -f##PBS -A EDU2#PBS -q genesis-b#PBS -l elapstim_req=00:10:00 要求時間
#PBS -b 1 ノード数
#PBS –l cpunum_job=24#PBS -T openmpi#PBS -v OMP_NUM_THREADS=3 スレッド数#PBS -v NQSV_MPI_VER=gdr/4.0.3/intel19.1.3-cuda11.2.1
module load load genesis/1.5.1-openmpi-4.0.3-intel19.1.3-cuda11.2.1_single
cd ${PBS_O_WORKDIR}
mpirun ${NQSII_MPIOPTS} –np 8 –npernode 8 -bind-to socket spdyn INP > LOG

プロセス数

run.shファイルを実行します（２０秒程度）% qsub run.sh
10

Step 1: 出力ファイルの見方
問題なく実行された場合はLOGという出力ファイルが出力されます
GENESISの出力ファイルは7段階([STEP 0]-[STEP 6])に分かれています[STEP 0] : 計算環境の確認[STEP 1] : 入力パラメータの確認[STEP 2] : 並列数（プロセス数、スレッド数）の確認[STEP 3] : 分子・エネルギー関数情報の確認[STEP 4] : 初期座標のエネルギー計算結果[STEP 5] : シミュレーション計算結果[STEP 6] : 終端処理と経過時間
11

Step 1: 計算時間の確認の仕方[STEP 6] : 経過時間を確認ください‘total time’はsetupも含めた全体の経過時間です
計算時間の見積もりを行うためには”dynamics”で示された時間を利用してください
経過時間は基本的に最大値が使われます。
Performance (ns/day)=Output_Time> Averaged timer profile (Min, Max)
単位＝秒
total time
=23.005
setup
=
1.164
dynamics=21.841 MDのメインループの経過時間
energy=13.200
integrator =7.475
pairlist=
1.191 (1.118,1.260)
energybond=0.023 (0.017,0.030)
angle=0.084 (0.049,0.126)
dihedral=0.228 (0.119,0.349)
nonbond=11.725 (11.193,12.128)
pme real=11.724 (11.192,12.127)
pme recip=
8.354 (8.341,8.370)
各プロセスでのsolvation=0.000 (0.000,0.000)
polar=0.000 (0.000,0.000) 最大値
non-polar=0.000 (0.000,0.000)
restraint=0.000 (0.000,0.000)
qmmm=0.000 (0.000,0.000)
integrator
constraint=
1.331 (1.298,1.397)
update=0.594 (0.585,0.609)
comm_coord=0.796 (0.725,0.903)(skip)

24*3600*0.01/${dynamics}

最大値と最小値の差はプロセス間のインバランスとなるので、時間が遅いときは確認が必要です12Step 2: 時間幅の変更% cd ../2_VVER_2.5fs% lsrun.sh sample% cp ../1_VVER/INP .

GENESISのIntegratorは、2.5fsでの計算が可能です
Jung et al. J. Chem. Phys. 148, 164109 (2018)
Jung et al. J. Chem. Theory Comput. 15, 84-94 (2018)vi等エディタでINPファイルを開いて、編集してください[DYNAMICS]
integrator
nstepstimestepeneout_period% qsub run.sh

====

VVER
4000
0.002580

ステップ数 (5000 → 4000)
ステップ幅 (0.002 → 0.0025)
エネルギー表示回数 (100 → 80)

13

Step 3: RESPAの利用(1)% cd ../3_RESPA_2.5fs% ls
run.sh sample% cp ../2_VVER_2.5fs/INP .

GENESISのIntegratorは、RESPA(Tuckerman et al., J. Chem. Phys. 97, 1990 (1992))の計算が可能です
GENESISでのRESPAについてJung et al. J. Chem. Phys. 148, 164109 (2018)

Fast motion

VVER

Step=1

Slow motion
Step=2OthersLf: Fast motion
Ls: Slow motion

Step=3

Step=4

TimeRESPA
Step=1

Step=2

Step=3

Step=4

Speed up

Fast motion : 結合性相互作用、非結合性相互作用（短距離）
Slow motion : 非結合性相互作用（長距離）、温度・圧力計算

14Step 3: RESPAの利用(2)vi等エディタでINPファイルを開いて、編集してください[DYNAMICS]
ntegrator=VRESnsteps=
4000
timestep= 0.0025
eneout_period=80
nbupdate_period=10
elec_long_period =2thermostat_period =2barostat_period=2% qsub run.sh

RESPAの利用 (VVER → VRES)

長距離非結合性相互作用計算の頻度(2 steps, 5fs毎)
温度計算の頻度(2 steps, 5fs毎)
圧力計算の頻度(2 steps, 5fs毎)Step 4: 温度・圧力計算の頻度を落とすRESPA% cd ../4_RESPA_2.5fs_tb10% lsrun.sh sample% cp ../3_RESPA_2.5fs/INP .

GENESISの温度・圧力計算の頻度は10(25fs)へ落とすことが可能です[DYNAMICS]
ntegrator=(skip)elec_long_period =thermostat_period =barostat_period=

VRES
210 温度計算の頻度(2 → 10)10 圧力計算の頻度(2 → 10)% qsub run.sh
16計算時間の比較% cd ../% grep “ dynamics” */LOG | sort –g –r –k 4 > timer.log

timer.logから、インプットを変更することで計算時間が短縮していることをご確認ください (sortコマンドによって計算時間が短いものが最後に来ます）
1_VVER/LOG:
dynamics=
21.841
2_VVER_2.5fs/LOG:
dynamics=
17.445
3_RESPA_2.5fs/LOG:
dynamics=
16.335
4_RESPA_2.5fs_tb10/LOG:
dynamics=
15.404
Tutorial 3.3で使われている入力パラメータ(Integrator, dt,…)は4_RESPA_2.5fs_tb10と同じものです17Step 5: ベンチマーク計算
効果的に研究を行うためには、ベンチマークを計測し効率的なコア数・並列パラメータを決定することが大事ですベンチマークの計算によって必要とされる資源量（時間、ノード時間）を知ることもできます
Cygnusは1ノード24コア持つため、プロセスとスレッドの組み合わせは(processes, threads)=(1,24), (2,12), (3, 8), (4, 6), (6, 4), (8, 3), (12, 2), (24, 1)
が理論上考えられます
数字を変えて試してみましょう （異なる名前のスクリプトを作成ください）#!/bin/bash -f#(skip)#PBS -b 1 ノード数（固定）#PBS -T openmpi#PBS -v OMP_NUM_THREADS=3 スレッド数(skip)
mpirun ${NQSII_MPIOPTS} -np 8 –npernode 8 -bind-to socket spdyn INP

全プロセス(MPI)数

ノード当たりのプロセス数>

p8t3.log18

Step 5: 計算が失敗する理由 (1)
いくつかの組み合わせでは、ログが以下の部分で止まります
Define_Molecule> Summary of molecules
num_atoms=24552 num_bonds=
num_angles=
9438 num_dihedrals=
num_impropers=
141 num_cmap_terms =num_residues=
7986 num_molecules=
num_segments=
3 num_deg_freedom =total_charge=-0.000

16628
2269
54
7931
73656

Setup_Restart_Pre> Coordinates and velocities were replaced

その場合は(script名).eXXXXX (XはJOBID)をご覧になってください19Case 1:

Step 5: 計算が失敗する理由 (2)Setup_Processor_Number> MPI Process number can notbe defined, please set them manualy
rank_no =2

ドメイン セル

GENESISのドメイン分割は、全体のMPIプロセスを各方向(X, Y, Z)に振り分けますが、自動で振り分けることができない場合にこのエラーが出ます (特に少ないプロセス数の時によく出ます）[BOUNDARY]
type
domain_x
domain_y
domain_z= PBC= 2= 2 各方向でのMPIプロセスの数= 3 domain_x * domain_y * domain_z =# of MPIs

ただし、domain_[x,y,z] = 1の場合は(domain_x, domain_y, domain_z) = (2, 1, 1) or (2, 2, 1) or (1, 1, 1)以外は使えません
つまり、3や6のMPIプロセスは計算できません
20

Step 5: 計算が失敗する理由 (3)Case 2:

Setup_Processor_Number> Cannot define domains and cells. Smaller or adjusted MPI
processors, or shorter pairlistdist, or larger boxsize should be used (see "Chapter:
Trouble shooting" in the user manual). rank_no =13

セルのサイズはドメイン数とペアリストの長さで決まります。適切なドメイン分割ができない場合にはこのようなエラーが出ます。MPIの数を減らすか、セットアップを見直してください。
ペアリスト長=13.5Åの場合、最短のセル長(Rcell, min)、実際のセル長(Rcell, i )、セル数(Ncell, i)はR Rbuff 13.5  2.6   8.05
22
N cell ,i  int( Bi / Rcell ,min )  mod  int( Bi / Rcell ,min ), domain i  i  x, y, z

R cell ,min 

pairlist

Rbuff= 2.6 (NPT & (const || water_model))Rbuff= 2.0 ((const || water_model))Rcell ,i  Bi / N cell ,iNcell, i >= 5、かつ、 Ncell, i > 2*domaini である必要があります
Rcell, i >> Rcell, min の場合にはメモリ使用量が大きくなりすぎたり、シミュレーションが遅くなります[STEP3]のSetup_Boundary_Cell>にこのシミュレーションでのドメインとセル数が出力されます
Setup_Boundary_Cell> Set Variables for Boundary Condition
domains (x,y,z) =22
1ncells (x,y,z) =66721

Step 5: ベンチマーク計算
計算ができたものは以下の５セットになるはずです(processes, threads)=(1,24), (2,12), (4, 6), (8, 3), (12, 2)
grep “ dynamics” p*log | sort –g –r –k 4 > timer.log

計算時間がプロセス、スレッド数の組み合わせによって大きく変わっていることをご確認ください
p1t24.log:
p2t12.log:
p12t2.log:
p4t6.log:
p8t3.log:

dynamics
dynamics
dynamics
dynamics
dynamics=====

39.60922.06217.31315.89215.404同じコア数でも２倍以上の差プロセス、スレッド数の組み合わせは、システムのサイズ、計算条件で大きく変化します
まずは5000-10000ステップで最適な組み合わせを探してください22

Hands-on 2

REMD計算のデモンストレーション
23

資料について今回の資料はGENESIS websiteのTutorial 10.1に基づいています
https://www.r-ccs.riken.jp/labs/cbrt/tutorials2019/tutorial-10-1/
初期構造作成などは、上のウェブサイトに詳しく記載されています

水中のペプチドの等温(NVT)条件下でのREMD計算

水溶液中のペプチド

24

ファイル
% cd ../../2_REMD% ls
1_REMD_DEMO toppar setup% cd 1_REMD_DEMO% ls
INP run.sh samplevi等エディタでINPファイルを開いてください25

コントロールファイル（抜粋）[OUTPUT]
dcdfile = step6_rep{}.dcd
rstfile = step6_rep{}.rst
remfile = step6_rep{}.rem
logfile = step6_rep{}.log(skip)

入出力ファイルは
header.N.sufの様にレプリカ番号が入る
それらの番号部分は｛｝で示す[REMD]
dimensionexchange_period
type1
nreplica1
parameters1
310.24 ¥= 1= 2000

REMDの変数（次元）の数

レプリカの交換頻度= temperature レプリカの種類
レプリカ数= 8= 300.00 302.53 305.09 307.65[DYNAMICS]
integrator= VRES
312.85 315.47 318.12
ステップ数はレプリカの交換頻度の倍数
nsteps= 4000
各レプリカが持つ値
timestep= 0.0025
継続行は”¥”でつなぐ（ただし、¥の後
eneout_period= 500に空白などをいれないこと、次行の先
出力頻度はレプリカの交換頻度の約数
crdout_period= 500
頭は一つ以上の空白を入れること）
rstout_period= 2000stoptr_period= 10
注意:デモ利用のため、
elec_long_period = 2
10.1のパラメータ(step6)からレプリカ数など
thermostat_period = 10
barostat_period= 10

変更しています26

スクリプトファイル
vi等エディタでrun.shファイルを開いてください
今回のデモでは1レプリカ当たり、4 プロセス * 3 スレッドで計算を行います
全体のプロセス数 = レプリカ数 (8レプリカ) * 4 プロセス = 32 プロセス
32 プロセス * 3 スレッド = 96 コア → 4ノード 利用#!/bin/bash -f##PBS -A EDU2#PBS -q genesis#PBS -l elapstim_req=00:10:00#PBS -b 4#PBS -T openmpi#PBS -v OMP_NUM_THREADS=3#PBS -v NQSV_MPI_VER=gdr/4.0.3/intel19.1.3-cuda11.2.1
module load genesis/1.5.1-openmpi-4.0.3-intel19.1.3-cuda11.2.1_single
cd ${PBS_O_WORKDIR}
mpirun ${NQSII_MPIOPTS} -np 32 -npernode 8 -bind-to socket –x PATH spdyn INP>

LOG27出力ファイルの見方 (1)
問題なく実行された場合はLOGという出力ファイル、更に各レプリカ毎のエネルギー、トラジェクトリなどが別ファイルで出力されます
REMDの場合も出力ファイルは7段階([STEP 0]-[STEP 6])に分かれています
並列度の情報はMDと同様に[STEP2]ですが、レプリカ毎の情報もあります
Setup_Mpi_Remd> Summary of Setup MPI
number of MPI processes=
number of MPI processes in one replica =number of OpenMP threads=
total number of CPU cores=

324396

レプリカの情報が実際にどう使われているかは[STEP3] に表示されます
Setup_Remd> Replica information
ParmsetID
1 =(skip)ParametersDim
1 =1300.000

302.530305.090307.650

310.240 …

28出力ファイルの見方 (2)
問題なく実行された場合はLOGという出力ファイル、更に各レプリカ毎のエネルギー、トラジェクトリなどが別ファイルで出力されます
REMDの場合も出力ファイルは7段階([STEP 0]-[STEP 6])に分かれています
並列度の情報はMDと同様に[STEP2]ですが、レプリカ毎の情報もあります
302.530K
307.650K
Rejected
Accepted
Replica 2
Replica 4
300.000K
305.090K
Replica 1
Replica 3
REMD> Step:
Replica
123(skip)

4000
Dimension:
ExchangeTrial
1 >
2
R
2 >
1
R
3 >
4A

Parameter:
RepIDtoParmID:
ParmIDtoRepID:

1

3番と４番の交換

R: rejected
A: accepted300.000
11302.530
22

ExchangePattern:
AcceptanceRatio
0 /
10 /
11 /
11Before
300.000
302.530305.090After
300.000
302.530307.650

各レプリカが担当する
変数(温度)の変化
307.650
44305.090
33312.850
66310.240
55318.120
88315.470
7729

最後に
GENESIS開発チームへの質問、バグ報告等はGENESIS forumをお使いください
Forum

日本語専用のフォーラムもあります

今日はありがとうございました！
30

謝辞
•

GENESIS開発
開発プロジェクトリーダー: 杉田 有治 博士 (理化学研究所）
開発者: （順不同、特に1.5.x、1.6.x）
Jaewoon Jung 博士 (理化学研究所）
八木 清博士 (理化学研究所）
森 貴治 博士 (理化学研究所）
尾嶋 拓博士 (理化学研究所）
伊東 真吾 博士 (理化学研究所）
神谷 基司
松永 康佑
松岳 大輔
優乙石博士
博士
博士
博士

(理化学研究所、現・分子科学研究所）(理化学研究所、現・埼玉大学）(理化学研究所、現・RIST）(理化学研究所、現・前橋工科大学）
31

講習会当日にあげられた質問補足(2021/09/28)••••

Q: RESPA計算で elec_long_period/thermo_period/baro_periodを設定しない場合はそれぞれいくつですか？
A: 1です。
Q: gREST+REUSは計算できますか？インプット例はありますか？
A: 可能です。GENESIS manualの15.4.4 gREST(p. 95)にインプット例があります。
Q: REMD計算のログにある”N”はなんでしょうか？
A: 交換の試行を行わないレプリカです。
Q: LAMMPSで扱うポテンシャル関数は計算可能ですか？
A: 多体ポテンシャルは計算できません。(講習会ではN-bodyと言っていましたが、Many-bodyのつもりでした。)

32