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

環境

macOSまたは Linux環境。
Gromacsがインストールされている
MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜の記事に沿うように書いてあります。
ここまででGROMACSを用いた簡単なエネルギー最小化とGROMACSを用いたNVT, NPT平衡化の処理を実行し終わった。

概要

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

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

膜タンパク質を使って膜と一緒にMDをやるときには別のmdpファイルを使ったを試したほうがよいので、それについてはググってください。GromacsのTutorialで有名なBevan Labのやつとかがいいと思います。

手順

球状タンパク質の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

実行スクリプトの例

#!/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^-1そのまま導入しています。デフォルト値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アルゴリズム内においてヤコビ行列の逆数に近似する展開を行っていますが、その展開する次数を制御します。

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

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