GROMACSを用いたNVT, NPT平衡化

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

## 環境

  - macOSまたは Linux環境。
  - [Gromacs 2023.2 がインストールされている](https://qiita.com/Ag_smith/items/2432375399e7d7868138)
  - [MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜](https://qiita.com/Ag_smith/items/2fe20010915ace67ed0a)の記事に沿うように書いてあります。
  - ここまでで[GROMACSを用いた簡単なエネルギー最小化](https://qiita.com/Ag_smith/items/3a438adb0dab4cfa06eb)の処理を実行し終わった。

## 概要

ここまででタンパク質と水分子の位置の最適化による[エネルギー最小化](https://qiita.com/Ag_smith/items/3a438adb0dab4cfa06eb)を行いました。しかしながら別にエネルギー最小化と言っても本当の最小値に至ることはほとんどありません。なぜなら系の自由度がとても高い上に、最急降下法（`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がかかっているため、タンパク質に対する影響は小さいと思いますが、水やイオンについてはよろしくなかったかも。）

ここでは圧力制御のための`pcoupl`は`no`に設定されていることに気をつけてください。また、`gen_vel`, `gen_temp = 100`と`gen_seed = -1`によって各原子に温度100 Kに対応したBoltzmann Distributionに基づくランダムな初速を与えることで熱運動を開始させています。タイムステップ<code>dt</code>は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`ディレクトリに作ります。

```bash
cd ~/Desktop/1LKE/gromacs/
mkdir heat
```

ここの`heat`ディレクトリに`md1`~`md9.mdp`ファイルと、次の`run.sh`を作成します。

``` bash
#!/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時間くらいはかかってしまいます。クラスタマシンやスパコンの利用をした方がいいですね。