1月15日開催!現年収非公開で企業からスカウトをもらってみませんか?PR

転職ドラフトでリアルな市場価値を測る。レジュメをもとに、企業から年収とミッションが提示されます。

6
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

GROMACSを用いたNVT, NPT平衡化

Last updated at Posted at 2018-06-02

このページではGROMACSでのNVT, NPT平衡化操作についての方法を説明します。
(2023年7月27日追記:md1.mdpファイル内にgen_tempの設定がなかったことを修正)
(2023年11月7日追記:各mdpファイルのpcouplを最新の方式であるc-rescaleに修正)

環境

概要

ここまででタンパク質と水分子の位置の最適化によるエネルギー最小化を行いました。しかしながら別にエネルギー最小化と言っても本当の最小値に至ることはほとんどありません。なぜなら系の自由度がとても高い上に、最急降下法(integrator = steepで指定していた)ではその収束のアルゴリズムはそれほどよいものではありません(代わりに、実装しやすいしやっている処理がわかりやすいというメリットがあります)。

とはいっても、別に最小値に行かせる必要はなく、これから行う平衡化処理をするうえで適度にエネルギーが低くなっていればいいな、という程度のものです。

さて、次にやることはNVT平衡化です。粒子数$N$、体積$V$、温度$T$の値が一定のまま行う平衡化です。ここで、系の圧力$P$は変わります。これはカノニカルアンサンブルで、熱浴(heat bath)との間でエネルギーを自由にやりとりできる系(閉鎖系)に対応する統計分布に対応します。この状態でまず温度を300 Kに十分に平衡化させる操作を行います。だいたい100 psくらいあれば十分と言われています。

実際の手順

私達がやるべきことは、

  1. 以上のNVT平衡化を行うための設定ファイルmdpファイルを作ること
  2. それをgmx gromppでGROMACSが読み込める形にすること
  3. 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がかかっているため、タンパク質に対する影響は小さいと思いますが、水やイオンについてはよろしくなかったかも。)

ここでは圧力制御のためのpcouplnoに設定されていることに気をつけてください。また、gen_vel, gen_temp = 100gen_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時間くらいはかかってしまいます。クラスタマシンやスパコンの利用をした方がいいですね。

6
8
5

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
6
8

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?