GROMACSを用いたProduction Runの方法

GROMACS 2019.2でのProduction Runについての方法を説明します。

## 環境

  - macOSまたは Linux環境。
  - [Gromacsがインストールされている](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)と[GROMACSを用いたNVT, NPT平衡化](https://qiita.com/drafts/24fbb9467787f69fbf1c)の処理を実行し終わった。

## 概要
ここまででトポロジー作成、エネルギー最小化、NVT&NPT平衡化の処理を終わらせていれば、タンパク質のMDシミュレーションの準備はほぼすべて終わりました。あとはひたすらProduction Runをするだけです。

今回説明するのは一般的な球状タンパク質の構造のMDシミュレーションのパラメータ設定例の`mdp`ファイルです。GROMACSにおいて拡張子が`.mdp`のものはMDのパラメータ設定ファイルとして認識されており、基本的にこれ以外のものは受け付けないことを覚えておきましょう。

膜タンパク質を使って膜と一緒にMDをやるときには別のmdpファイルを使ったを試したほうがよいので、それについてはググってください。GromacsのTutorialで有名な[Bevan Lab](http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/)のやつとかがいいと思います。

## 手順

球状タンパク質のProduction Run用のmdpファイル(`mdrunset.mdp`)はこちら。

    title       = hogehoge
    define        = 
    ; Run parameters
    integrator      = md        ; leap-frog integrator
    dt              = 0.002     ; 2 fs
    nsteps          = 50000000  ; 2 fs * 50000000 = 100 ns
    cutoff-scheme   = Verlet
    
    ; Output control
    nstxout     = 5000            ; save coordinates every 10 ps
    nstvout     = 5000            ; save velocities every 10 ps
    nstenergy   = 2500            ; save energies every 5 ps
    nstlog    = 5000              ; update log file every 10 ps
    
    ; Bond parameters
    continuation    = yes               ; Restarting after previous RUN 
    constraint_algorithm = LINCS    ; holonomic constraints 
    constraints      = hbonds     ; Convert the bonds with H-atoms to constraints
    lincs_iter       = 1            ; accuracy of LINCS
    lincs_order      = 4            ; also related to accuracy
    
    ; Langevin Dynamics
    bd_fric = 2.0       ; Brownian dynamics friction coefficient(correspond to gamma_ln in AMBER) 
    ld_seed = 1732  ; random seed
    
    ; Neighborsearching
    ns_type         = grid      ; search neighboring grid cells
    nstlist         = 20          ; 10 fs
    rlist           = 1.0       ; short-range neighborlist cutoff (in nm)
    rcoulomb        = 1.0       ; short-range electrostatic cutoff (in nm)
    rvdw            = 1.0       ; short-range van der Waals cutoff (in nm)
    
    ; Electrostatics
    coulombtype = PME           ; Particle Mesh Ewald for long-range electrostatics
    pme_order   = 4                 ; cubic interpolation
    fourierspacing  = 0.12      ; grid spacing for FFT(nm)
    
    ; Temperature coupling is on
    tcoupl    = V-rescale    ; modified Berendsen thermostat
    tc-grps   = Protein Water_and_ions   ; two coupling groups - more accurate
    tau_t     = 0.1  0.1           ; time constant, in ps
    ref_t     = 300  300       ; reference temperature, one for each group, in K
    
    ; Pressure coupling is on
    pcoupl          = Parrinello-Rahman ; Pressure coupling on in NPT
    pcoupltype      = isotropic     ; uniform scaling of box vectors
    tau_p           = 2.0     ; time constant, in ps
    ref_p           = 1.0     ; reference pressure, in bar
    compressibility = 4.5e-5  ; isothermal compressibility of water, bar^-1
    refcoord_scaling = com
    
    ; Periodic boundary conditions
    pbc         = xyz       ; 3-D PBC
    
    ; Dispersion correction
    DispCorr     = EnerPres ; account for cut-off vdW scheme
    
    ; Velocity generation
    gen_vel         = no       ; Velocity generation
    gen_temp        = 300.0    ; temperature for Maxwell distribution
    gen_seed        = 1732     ; used to initialize random generator for random velocities
    
    ; 以下はdistance restarintsのときに使われる設定
    ; NMR refinement stuff  = 
    ; Distance restraints type: No, Simple or Ensemble = 
    disre                    = simple
    ; Force weighting of pairs in one distance restraint: Equal or Conservative = 
    disre-weighting          = Equal
    ; Use sqrt of the time averaged times the instantaneous violation = 
    disre-mixed              = no
    disre-fc                 = 1000
    disre-tau                = 0
    ; Output frequency for pair distances to energy file = 
    nstdisreout              = 5000

実行スクリプトの例

``` bash
#!/bin/bash

input=1lke

# 同じディレクトリにmd${input}.tprがない場合のみ, gmx gromppコマンドを実行する
# gmx gromppコマンドは平衡化終了後のmd9.tpr, md9.cpt、そして上述のmdrunset.mdpを読み込み、
# GROMACSでのMDを始めるのに必要なmd${input}.tprファイルを生成する。
if [ ! -e md${input}.tpr ];then
    gmx grompp \
        -f mdrunset.mdp \
        -c ../heat/md9.tpr \
        -r ../heat/md9.tpr \
        -t ../heat/md9.cpt \
        -p ../../top/${input}.top \
        -o md${input}.tpr \
        -po md${input}.mdp
fi

# gmx mdrunでproduction runを始める
# 標準出力をstdout.logファイルに、標準エラー出力をstderr.logファイルに書き出す
# 場合によっては -ntmpi 0 も入れておいた方がいいかも
gmx mdrun \
 -v -s md${input} -deffnm ${input} -cpi ${input}.cpt 1> stdout.log 2> stderr.log
```

## 雑な解説
mdpファイルの雑な解説をしていきます。まず、このファイル内では`;`マークによってコメントが可能です。

まずTitleですがこれは適当に設定して構いません。っていうか不要？ `define =　`は上では空欄ですが、`define = -DPOSRES1000`のように設定してあげると、.topファイル内に記述された

    ; Include Position restraint file
    #ifdef POSRES1000
    #include "posre1000.itp"
    #endif

というif文の設定が有効になります。このやり方を使って実はエネルギー最小化や平衡化のときにposition restraintが発動するようにし、重原子を少しずつ動かすようにうまく制御していました。

次にRun parametersのところで、`dt`の値を0.002にするのは現状の技術ではこれがベターだと言われています。これ以上大きくするとLeap-frog法による計算の収束ができなくなってしまうことが多いらしいです。0.002を可能にしているのはSHAKEアルゴリズムとかLINCSアルゴリズムとかのおかげです。値を0.003にすることも可能らしいですが、0.002なら切りがいいところで止めやすいので、こちらを強く勧めます。

`dt`の値と`nsteps`の値を掛けた値がMDで行うシミュレーション時間になります。**GROMACSは最初にProduction runの時間を設定しておき、その時間に到達したら計算を終了するというシステムを採用しています**。例えば、10ナノ秒だけ動かすよう`mdrunset.mdp`に設定して`gmx grompp`コマンドで`.tpr`ファイルを生成したとすると、`gmx mdrun`コマンドでのProduction Runは、一度その時間に到達すると以降`gmx mdrun`を何度動かそうとしてもそれ以上の時間分計算してくれることはありません。
このmdpファイルによる作成をしてしまったあとで、「もうちょっとシミュレーション時間を延長したいな～」ってなるのはよくありますが、それは`tpr`ファイルの設定変更に使われる`gmx convert-tpr`コマンドによって変更することが可能なので、覚えておきましょう。

続いてアウトプットファイルの設定です。GROMACSのMDのトラジェクトリの拡張子は基本的には`.trr`です。`nstxout`と`nstvout`の値はそれぞれそのアウトプット頻度を設定します。`dt`の値を1stepとして5000
stepsごとですと、つまり10 psごとにセーブということになりますね。しかし、MDシミュレーションは10
psごとに全原子の座標データ、速度データをセーブしていても、100 ns（＝100,000 ps）の頃には数十GBとかにも到達します。かといって減らしすぎると変化を追いにくくなるので注意です。
この`.trr`ファイルとは別に、座標データ専用の`.xtc`拡張子ファイルというものを別個に吐き出させることも出来ます。こちらには速度情報は含まれません。これは`nstxtcout`の値によって変更することが可能です。`.trr`ファイルと比べて速度情報が含まれない分、シミュレーション後に座標データのデータ処理を行いたいときは`.xtc`ファイルのほうが速いです。……がまあ私は`.trr`だけあれば十分じゃないかなと思います。

`; Velocity generation`では各原子の初速度の生成を行うかどうかを設定します。もしNPT平衡化などですでに前のシミュレーションを行っている場合、その終了時の速度情報をそのまま流用することができますし、そのほうがオススメです。上の設定例では`gen_vel`をnoにしており、代わりに`continuation = yes`としております。`gen_vel`をyesとしている場合、`gen_temp`と`gen_seed`の値によってランダムに各原子に初速度を与えてMDシミュレーションを開始することになります。`gen_seed`の値は本当に気分でどんな値を入れても構いません。

`; Langevin Dynamics`の部分は、水中で分子がブラウン運動をするのを再現するための設定です。私はAMBERでやっていたときの設定値を2.0 amu/ps<sup>-1</sup>そのまま導入しています。デフォルト値0ですと、自動的に（各原子の質量）/（`tau-t`値）に割り振られるみたいです。

`; Bond parameters`における`constraint_algorithm = LINCS`は、AMBERなど他のMDソフトに比べてGromacsが超速で動く原理の1つを担っています。AMBERではSHAKEアルゴリズムによる拘束を行っていましたが、Gromacsに実装されているLINCSはその3～4倍早く演算を行えるようになります。
原子と原子の伸縮運動のタイムスパンは水素原子みたいな軽い原子が絡んでいると短くなります。そのため、正確に運動方程式を解いていくためにはleap-frog法では本来もっと短い`dt`値が必要となります。しかし、水素の結合伸縮は全体の運動から見るとそこまで重要でないことが多いため、この結合長を固定しておくためのアルゴリズムを導入すると、`dt`値を0.002にしても大丈夫になり、シミュレーション時間を大幅に短縮することが可能となります。
LINCSアルゴリズムは距離拘束専用であり、SHAKEでできていた角度拘束はできません。一応、距離拘束を組み合わせることで擬似的に角度拘束も可能ですが、注意点として、1原子に対してあまりに多い拘束が重なってしまうと使えなくなることがあります（分子の環構造なども含む）。まあ大抵の場合は大丈夫ですけど。

`lincs_order = 4`というのはLINCSアルゴリズム内においてヤコビ行列の逆数に近似する展開を行っていますが、その展開する次数を制御します。

```math
({\bf I}-{\bf A}_n)^{-1} = {\bf I}+{\bf A}_n+{\bf A}_n^2+{\bf A}_n^{3}+\cdots
```

通常のMDシミュレーションだと4次の展開で十分ですが、大きなタイムステップでブラウン運動シミュレーションを見たい場合は8次くらいの展開が必要になります。

`;Temperature coupling`と`;Pressure coupling`にはそれぞれ使いたいアルゴリズムを選んでください。Temperature coupling（系の温度制御）ではNose–Hoover法もひとつ良いやり方です。ただし、Temperature couplingでberendsenを使用するくらいならば、V-rescale法を選択すべきです。こちらは、<b>berendsen法による温度制御では実はcanonical ensembleになっていない</b>という致命的な弱点を克服した新規手法です。
`;Pressure copling`（系の圧力制御）の設定で`pcoupl = berendsen`を使う方法について、これまでの平衡化の時にこちらのアルゴリズムを使用しておりました。<b>実はこちらも正しい canonical ensembleを得られない</b>という致命的弱点があるのですが、圧力を一定に安定化させるのが非常に早い手法であるため、サンプリング前である平衡化のときにはこちらを使って早く収束させようとしていました。繰り返しますが、平衡化のときには`pcoupl = berendsen`を使うべきです。 しかし、本番のProduction runではcanonical ensembleを得るために（＝より現実に即したサンプリングになるように）、`pcoupl = Parrinello-Rahman`を使う方がよいでしょう。こちらのアルゴリズムは`pcoupl = berendsen`よりは遅くなります。また、こちらのアルゴリズムを使う時に初期値が設定値から大きく異なる場合、圧力が安定するまでに大きな圧力値の振動が発生し、場合によっては計算が落ちてしまうことがあります。上述したように、先に`berendsen`でささっと設定値に収束させておき、そこから`Parrinello-Rahman`に切り替えましょう。

あとは各自使いたい設定に書き換えて利用してください。