このページではGROMACSでのNVT, NPT平衡化操作についての方法を説明します。
(2023年7月27日追記:md1.mdp
ファイル内にgen_temp
の設定がなかったことを修正)
(2023年11月7日追記:各mdp
ファイルのpcouplを最新の方式であるc-rescale
に修正)
環境
- macOSまたは Linux環境。
- Gromacs 2023.2 がインストールされている
- MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜の記事に沿うように書いてあります。
- ここまででGROMACSを用いた簡単なエネルギー最小化の処理を実行し終わった。
概要
ここまででタンパク質と水分子の位置の最適化によるエネルギー最小化を行いました。しかしながら別にエネルギー最小化と言っても本当の最小値に至ることはほとんどありません。なぜなら系の自由度がとても高い上に、最急降下法(integrator = steep
で指定していた)ではその収束のアルゴリズムはそれほどよいものではありません(代わりに、実装しやすいしやっている処理がわかりやすいというメリットがあります)。
とはいっても、別に最小値に行かせる必要はなく、これから行う平衡化処理をするうえで適度にエネルギーが低くなっていればいいな、という程度のものです。
さて、次にやることはNVT平衡化です。粒子数$N$、体積$V$、温度$T$の値が一定のまま行う平衡化です。ここで、系の圧力$P$は変わります。これはカノニカルアンサンブルで、熱浴(heat bath)との間でエネルギーを自由にやりとりできる系(閉鎖系)に対応する統計分布に対応します。この状態でまず温度を300 Kに十分に平衡化させる操作を行います。だいたい100 psくらいあれば十分と言われています。
実際の手順
私達がやるべきことは、
- 以上のNVT平衡化を行うための設定ファイル
mdp
ファイルを作ること - それを
gmx grompp
でGROMACSが読み込める形にすること - GROMACSを用いて実行することです。
そのmdpファイル(md1.mdp
)はこちらとなります。
; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)
define = -DPOSRES1000
; RUN CONTROL PARAMETERS
integrator = md
dt = 0.002
nsteps = 100000
;
nstxout = 500
nstlog = 500
nstenergy = 500
;
nstlist = 5
ns_type = grid
pbc = xyz
rlist = 0.8
;
coulombtype = PME
rcoulomb = 0.8
;
vdwtype = cut-off
rvdw = 0.8
dispcorr = enerpres
;
optimize_fft = yes
;
tcoupl = v-rescale
tc_grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300.0 300.0
;
pcoupl = no
pcoupltype = isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = all
;
gen_vel = yes
gen_temp = 100
gen_seed = -1
;
continuation= no; Restarting after previous RUN
constraints = hbonds
constraint_algorithm = LINCS
annealing = single single
annealing_npoints = 2 2
annealing_time = 0 200 0 200
annealing_temp = 100 300 100 300
(2023年7月27日追記:100 KのBoltzmann Distributionに基づいて初速を生成するように修正。 これまではgen_temp
デフォルトの300 Kが自動選択されており、300 Kで初速が生成された後に0 Kまで急激に冷やされてる形になっており、望ましい計算になっていませんでした。これまでもタンパク質にはposition restraintsがかかっているため、タンパク質に対する影響は小さいと思いますが、水やイオンについてはよろしくなかったかも。)
ここでは圧力制御のためのpcoupl
はno
に設定されていることに気をつけてください。また、gen_vel
, gen_temp = 100
とgen_seed = -1
によって各原子に温度100 Kに対応したBoltzmann Distributionに基づくランダムな初速を与えることで熱運動を開始させています。タイムステップdt
は0.002 ps, それを100000 steps行うので、つまり200 ps間シミュレーションを行います。他のパラメータの意味については説明書を見てください。
続いて、そのままNPT平衡化を行うためのファイル(md2
~md9.mdp
)を作ります。
; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)
define = -DPOSRES1000 (-DPOSRES500, -DPOSRES200, -DPOSRES100, -DPOSRES50,-DPOSRES20,-DPOSRES10,-DPOSRES0 )
; RUN CONTROL PARAMETERS
integrator = md
dt = 0.002
nsteps = 50000
;
nstxout = 500
nstlog = 500
nstenergy = 500
;
nstlist = 5
ns_type = grid
pbc = xyz
rlist = 0.8
;
coulombtype = PME
rcoulomb = 0.8
;
vdwtype = cut-off
rvdw = 0.8
dispcorr = enerpres
;
optimize_fft = yes
;
tcoupl = v-rescale
tc_grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300.0 300.0
;
pcoupl = c-rescale
pcoupltype = isotropic
tau_p = 1.0
compressibility = 4.5e-5
ref_p = 1.0
refcoord_scaling = all
;
gen_vel = no
;
continuation= yes; Restarting after previous RUN
constraints = hbonds
constraint_algorithm = LINCS
(2023年11月7日: pcouplの設定値を、最近使用可能になったおすすめ設定のc-rescale
に修正)
ここで、md2.mdp
, md3.mdp
, md4.mdp
, ... , md9.mdp
を作るときに、-DPOSRES
の値をそれぞれ上に対応するように0に近づくよう1つ1つ変えていってください。こうすることで、以前1lke.top
の中に追記したposition restraint
がそれぞれ有効になります。
これら9つのmdpファイルを以下の~/Desktop/1LKE/gromacs/heat
ディレクトリに作ります。
cd ~/Desktop/1LKE/gromacs/
mkdir heat
ここのheat
ディレクトリにmd1
~md9.mdp
ファイルと、次のrun.sh
を作成します。
#!/bin/sh
input="1lke"
prev_tpr="../minimize/min2.gro"
prev_cpt="../minimize/min2.trr"
for (( i=1;i<10; i++));do
if [ $i -ne 1 ];then
j=$((i-1))
prev_tpr=md$j.tpr
prev_cpt=md$j.cpt
fi
#まだmd${i}ステップが終わっていなければ(md${i}.groが存在していないなら)実行
if [ ! -f md${i}.gro ];then
rm -f md$i.out.mdp md$i.tpr
gmx grompp \
-f md$i.mdp \
-c $prev_tpr \
-r $prev_tpr \
-p ../../top/${input}.top \
-n ../../top/index.ndx \
-t $prev_cpt \
-o md$i.tpr \
-po md$i.out.mdp || exit $?
rm -f md$i.trr md$i.cpt md$i.gro md$i.edr md$i.log
gmx mdrun \
-v -s md$i -deffnm md$i -cpi md$i.cpt
#すでに${i}ステップが終わっていれば(md${i}.groが存在していれば)次のサイクルへ
else
echo "md${i}.gro present. Start next cycle."
fi
done
Gromacs-2018からはgrompp
中に-r $prev_tpr
を入れる必要ができたみたいです。
以上10のファイルを置いて、
chmod +x run.sh
./run.sh
とすればこれらの操作が一気に行われ、NVT平衡化、NPT平衡化が一気に終わるはずです。
なんですけど、これらの処理は非常に重く、また計算時間はお使いのCPUのコア数に反比例します。最新のMacBook Proは4コアですが、これでも1時間くらいはかかってしまいます。クラスタマシンやスパコンの利用をした方がいいですね。