基本的なGromacsの使い方を勉強してみる


オンラインにGromacs で水の凝集状態をMDするチュートリアルを見つけたので、トレースしてみる。

参照1
参照2
参照3

MDの開始構造の作成

最初のエネルギー最小化＋平衡化の処理

0.最初の準備

これまでにCondaパッケージを利用してCentOS8(GPUなし)上にGromacsをインストールしている

  375  gmx --version

上記コマンドが動くことを確認する。GROMACSの実行により多くのファイルが生成するため、実行ディレクトリ（フォルダ）は専用のものにしておく。

  377  mkdir 200927_gromacs_tutorial
  378  cd 200927_gromacs_tutorial
  380  vim system.top

以下のファイルをvimやemacsなど、Linux上のテキストエディタで作成します。

  • system.top
  • em.mdp
  • nvt.mdp
  • sampling.mdp

1 相互作用パラメータ決め

まず、以下のような内容の system.top ファイルを作成します。vimで作成するText内容

#include "amber03.ff/forcefield.itp"
#include "amber03.ff/spce.itp"

[ system ]

water in water

[ molecules ]

2 初期座標決め

水分子を適当に並べ、初期構造とします。gmx solvate コマンドで作成、4 nm 立方の箱の中に水分子を並べる。

gmx solvate -cs spc216.gro -o init.gro -box 4.0 -p system.top

水分子の座標を記録しているファイルinit.groというファイルが新しくでき、system.top の中身をみて、最終行に"SOL；****"が追加されることを確認。

3 MDの条件決め

以下のような内容で em.mdp というファイルを作成します。vimで作成するText内容

#include "amber03.ff/forcefield.itp"
#include "amber03.ff/spce.itp"

[ system ]

water in water

[ molecules ]

4 MD実行

設定ファイルの「コンパイル」；以上の段階で「初期構造」「力場パラメータ」「MD条件」の３つのファイルを作ったわけですが、これらから"gmx grompp"コマンドでMDを行うプログラムが読み込めるような形式のバイナリファイルを作ります。

gmx grompp -c init.gro -p system.top -f em.mdp -o step1_em.tpr

mdout.mdp と step1_em.tpr の２つのファイルが生成されたことを確認して、さらに以下を実行する。

gmx mdrun -deffnm step1_em

エラーなく終われば、以下のファイルが生成されているはず。

  • step1_em.edr
  • step1_em.log
  • step1_em.trr
  • step1_em.gro
  • step***.pdb (step16c.pdbなど。ない場合もある)

このうち step1_em.groがエネルギー最小化の結果の構造を記録したファイルで、次のMDで使うものになります。

5 結果の解釈

実際にエネルギーはどの程度小さくなっているか。次のコマンドで必要な情報を取出。

gmx energy -f step1_em.edr

6. (ここまで)実行したコマンド

  382  less system.top
  383  gmx solvate -cs spc216.gro -o init.gro -box 4.0 -p system.top
  385  less em.mdp
  387  gmx grompp -c init.gro -p system.top -f em.mdp -o step1_em.tpr
  389  less step1_em.tpr
  390  gmx mdrun -deffnm step1_em
  392  gmx energy -f step1_em.edr

MDの実行(温度を指定しての実行)

1 相互作用パラメータ決め

最初のエネルギー最小化で作成した system.top をそのまま使います。

2 初期座標決め

平衡化の処理で得た最後の構造 step2_nvt.gro が初期構造となります。

3 MDの条件決め

以下の内容で sampling.mdp というファイルを作ります。

integrator  = md-vv
dt          = 0.001
nsteps      = 100000
nstenergy   = 100
nstxout     = 100

coulombtype             = PME

4 MD実行

まずは読み込みを行います。

gmx grompp -c step2_nvt.gro -p system.top -f sampling.mdp -o step3_sampling.tpr

ここでもNOTEが出てきますが、前回同様に無視します。これで step3_sampling.tpr が生成されましたら、いよいよサンプリングMDを実行します。

gmx mdrun -deffnm step3_sampling

通常は3分～10分程度で終わる。今回は貧弱な仮想マシンで実行したので、より時間がかかっている(停止することはなかった)

5．MDの実行(温度を指定しての実行)結果の解釈；動径分布関数

動径分布関数は、原子間の距離を測り、その分布を調べたものになります。これが分かると分子がどれくらい配向しているかなど、さまざまなことが推測できます。

まず、次のコマンドを実行します。

gmx make_ndx -f init.gro

すると、 > 出力が出てくきて、プログラムがこちらの入力を待っている状態となる。ここで

> a OW
> a HW1 HW2
> q

と入力します。すると最後のqで対話モードが終了し、正しく処理されていれば index.ndx というファイルができています。 これで index.ndx ファイルができたのち、

gmx rdf -f step3_sampling.trr -s step3_sampling.tpr -n index.ndx -o rdf_OH.xvg -excl

を実行します。

Available static index groups:
 Group  0 "System" (6495 atoms)
 Group  1 "Water" (6495 atoms)
 Group  2 "SOL" (6495 atoms)
 Group  3 "OW" (2165 atoms)
 Group  4 "HW1_HW2" (4330 atoms)
Specify a selection for option 'ref'
(Reference selection for RDF computation):
(one per line,  for status/groups, 'help' for help)
>

対話モードでどのグループを選ぶかを入力する状態になるので酸素原子と水素原子をを選びます。(あるいは、単に対応するグループ番号の 3 と 4 を入力しても大丈夫)

> OW
> HW1_HW2

このあと Ctrlを押しながらDキーを押すと、軌跡ファイルの読み込みをして動径分布関数の計算に入る。

6．軌跡の可視化

動径分布関数の計算が終わると rdf_OH.xvg というファイルができるので、grace(Xmgrace)で可視化して確認をする(未確認)。軌跡のデータは step3_sampling.trr にあります。 少々重いですので、コピーなどする際には気を付けましょう。 これと最初の座標データ step2_nvt.gro があれば、VMD を使って分子の運動を動画で表示できます。

7．(ここまで)実行したコマンド

  396  vim nvt.mdp
  398  gmx grompp -c step1_em.gro -p system.top -f nvt.mdp -o step2_nvt.tpr
  401  gmx mdrun -deffnm step2_nvt
  402  gmx energy -f step2_nvt.edr -o temperature.xvg
  403  vim sampling.mdp
  404  gmx grompp -c step2_nvt.gro -p system.top -f sampling.mdp -o step3_sampling.tpr
  406  gmx mdrun -deffnm step3_sampling
  407  gmx make_ndx -f init.gro
  410  gmx rdf -f step3_sampling.trr -s step3_sampling.tpr -n index.ndx -o rdf_OH.xvg -excl

(ここまで)　紹介されている操作をトレースできた。環境はほぼ構築したが、grace(Xmgrace)または代替は用意していない。一般ユーザー権限でのCondaパッケージで入れたGromacsでも一応動く(スペックは低い)ことは確認できた。