前回の記事(AmberTools 16を使ってMD simulationをしてみよう(その2))の続き。
前回、分子動力学シミュレーション(MD)を始めるために必要な、トポロジーファイル(leap.top)と初期座標ファイル(leap.crd)を作りました。MDをやりたい対象の蛋白質を変えたらなかなかこの2つを作るだけでも大変なのですがね。

さて、MDシミュレーションをやるにあたって、やるべきことをおさらいしてみましょう。
  1. LEaPによってトポロジー・初期座標ファイルを作る
  2. 系のエネルギー最小化
  3. 系の平衡化(NVTアンサンブル)
  4. 系の平衡化(NPTアンサンブル)
  5. Production MDの開始(データサンプリングを開始する)


今、1.の部分が終わりました。次にやるべきことはエネルギー最小化と平衡化です。
エネルギー最小化でやるべきことは、系の中に存在する全原子の位置を(ある程度)最適化することです。これをなしにMDをやることはできません。なぜなら、LEaPで生成された初期の水分子や原子の配置は、全然最適化されていません。特に、水分子の位置は氷のように本当に規則正しく並んでしまっています。私達がやりたいのは溶液中・常温でのシミュレーションなので、いきなり氷の状態から溶液のシミュレーションを始めることはできません。また、蛋白質周りの水素結合の形成が全然できていません。水素結合は非常に重要な結合であり、系のエネルギーの安定化にも大きく寄与しています。

系の平衡化のNVTアンサンブルっていうのはそれぞれN(分子数),V(体積),T(温度)が一定の条件のシミュレーションって意味です。Pは圧力です。これらを経て、系を氷の中の状態から徐々に溶液の状態へと変化させ、平衡化させるためのものです。平衡化が終わって、ようやくMDシミュレーションとしてのデータサンプリングを開始することができます。

というわけで、エネルギー最小化をGromacsを使ってやってみようと思います。AmberTools 16でもできるんですがそれはまたその機会に。

さて、以前の記事で、AmberTools 16を用いて力場を作成しました。しかし、AmberToolsで使われているファイルフォーマットとGromacsのそれは当然全然違います。そこで、以前紹介したACPYPE.pyを使ったAMBER→Gromacs topologyの変換を用いて、フォーマットを変換します。

0. 環境

1. ファイル形式の変換

まずはACPYPE.pyを用いてファイル形式を変換します。ここでは、Desktopに以前作ったleap.watディレクトリが存在し、その中にleap.top, leap.crdファイルがあるものとします。そこに、新たにgrotopディレクトリを作成し、その中でGromacs用のファイルを用意することにします。
cd ~/Desktop
mkdir grotop
cd grotop
curl -o acpype.py "http://acpype.googlecode.com/svn/trunk/acpype.py"
python acpype.py -p ../leap.wat/leap.top -x ../leap.wat/leap.crd
これにより、leap_GMX.gro, leap_GMX_top, md.mdp, em.mdpの4つのファイルが生成されます。
ここで、leap_GMX.gro, leap_GMX.topってのは後々のためリネームしておきます。初心者に説明しておくと、ターミナル上でのリネームはmvコマンドを使います。
mv leap_GMX.gro 1c08.gro
mv leap_GMX.top 1c08.top
後の2つのファイルmd.mdp, em.mdpは使わないので削除しておいてください。また、以前の記事で説明したACPYPE.pyの注意点を一読しておいてください。厳密にやろうとするときは.topファイル内の一部を書き換えておいたほうが無難です。ここでは省略します。

2. 位置拘束ファイルの作成

のちのち使うことになる、Gromacs用の位置拘束(position restraint)のファイルを作成します。位置拘束というのは非常に重要で有用なもので、特定の原子に任意の力の重みを付けて、動きにくくすることができます。

上述のエネルギー最小化の操作および平衡化の操作で、はじめに蛋白質の特に重原子(水素原子以外)に重みをつけながらシミュレーションを行うと、水溶媒や水素原子だけが動くことになり、水素結合を効率よく形成させることができるようになります。この目的のために位置拘束がとても使えます。

この位置拘束のファイルを生成するコマンドがGromacsには用意されています。
gmx genrestr -f 1c08.gro -o posre1000.itp -fc 4184 4184 4184
詳しい使い方は例によってgmx genrestr -hとすれば詳しい使い方の説明がでてきます。

さて、これをやると対話式に続けて入力を求められます。
Reading structure file
Select group to position restrain
Group 0 ( System) has 90766 elements
Group 1 ( Protein) has 5336 elements
Group 2 ( Protein-H) has 2733 elements
...
Group 14 ( Cl-) has 20 elements
Group 15 ( Water) has 85395 elements
Group 16 ( SOL) has 85395 elements
Group 17 ( non-Water) has 5371 elements
Group 18 ( Water_and_ions) has 85430 elements
Select a group:
ここで、Group 2のProtein-Hを選択します。これは、蛋白質を構成する原子のうち、水素原子を除いたものという意味です。今回は重原子だけを固定したいので、これを選択します。
そこで、2を入力してEnterキーを押すと、posre1000.itpが生成されます。これを開くと、以下のようになっています。
[ position_restraints ]
; i funct fcx fcy fcz
2 1 4184 4184 4184
5 1 4184 4184 4184
6 1 4184 4184 4184
7 1 4184 4184 4184
9 1 4184 4184 4184
これらの左端の値は入力のgroファイル(座標ファイル)の番号に対応、2番めのカラムは1で通常固定(詳しくはマニュアル)、3番目〜5番目はそれぞれx, y, z方向への位置拘束をかける力の値(kJ/mol nm^2)です。ちなみに10 kcal/mol A^2 = 4184 kJ/mol nm^2です。特に10 kcal/mol A^2にする意味は無いので適当な値を代入して構わないのですが、初めは大きな値を、それから徐々に半減させていくとよいでしょう。

ということでこのファイルを-fc の値を1000 1000 1000, 500 500 500, 200..., 100..., 50..., 20..., 10..., 0...としてそれぞれposre1000, 500, 200,,,0までつくります。シェルスクリプトの知識があればforやforeachなどでさくっと作れるでしょう。例えばbash系だったら
#!/bin/sh

array=(4184 2092 837 418 209 84 42 0)
array2=(1000 500 200 100 50 20 10 0)

for ((i=0; i < ${#array[@]}; i++));do
gmx genrestr -f 1c08.gro -o posre${array2[$i]}.itp -fc ${array[${i}]} ${array[${i}]} ${array[${i}]}<<EOF
2
EOF
done

みたいな感じですかねえ。

3. .topファイルへの設定追記

1c08.topには力場・電荷パラメータが入っています。初心者のうちはむやみにいじると死にます。
が、慣れてくると、割と書き換えやすく、わかりやすい記述になっていることに気付くでしょう。

ここでは必要な記述だけにしておきます。

まず操作を始める前にバックアップをとっておきます。これは大事。
cp 1c08.top 1c08.top.bak
次に、以下の赤色の記述をファイル1c08.topの一番下の方に追記します。
  1848   1852   1850   1851      1   180.00   4.60240   2 ;      C-    CA-     N-     H
1852 1864 1862 1863 1 180.00 43.93200 2 ; CA- N- C- O
1862 1866 1864 1865 1 180.00 4.60240 2 ; C- CH3- N- H

; Include Position restraint file
#ifdef POSRES1000
#include "posre1000.itp"
#endif
#ifdef POSRES500
#include "posre500.itp"
#endif
#ifdef POSRES200
#include "posre200.itp"
#endif
#ifdef POSRES100
#include "posre100.itp"
#endif
#ifdef POSRES50
#include "posre50.itp"
#endif
#ifdef POSRES20
#include "posre20.itp"
#endif
#ifdef POSRES10
#include "posre10.itp"
#endif
#ifdef POSRES0
#include "posre0.itp"
#endif


[ moleculetype ]
; molname nrexcl
Na+ 1

[ atoms ]
Gromacsのtopファイルは大体の場合、かぎかっこ[ hogehoge ]を用いてそれぞれのセクションを管理しています。今回挿入するべきところは、蛋白質の力場の二面角[ dihedrals ]が記述されている箇所が終わった直後です。次のからは[ moleculetype ]は蛋白質ではなくイオンのパラメータ管理に入っています。
ここに以下の赤色の記述を入れます。#ifdefは条件分岐であり、直後に定義された任意の配列が呼び出されれば、その値が有効になり、それ以外は無効にするというものです。#include では.itpファイルを呼び出しています。先ほどgmx genrestrによって.itpファイルを生成していましたが、#includeを使うことで.topファイル中に.itp形式の指令を代入することが簡単にできます。

4. .topファイルの細かい修正

AMBERのff14SBとGromacsのtopologyファイルは完全に互換、というわけではないので、いくつか細かい修正を入れる必要があります。
以前の記事にも書いたのですが、Gromacsでは.topファイル内における[ atomtype ]にて、原子種の指定の1文字目に数字を置くことは禁止されています。なので、これを一括で変更します。
sed -i -e "s/ 2C/C2C/g" 1c08.top
sed -i -e "s/ 3C/C3C/g" 1c08.top
のように、原子種2C3Cの1文字目をなんらかの方法でアルファベットにします。スペースを挟んでsedしているのがミソです。こうしていないと後々
-------------------------------------------------------
Program grompp, VERSION 4.6.6
Source code file: /users/agsmith/apps/gmx4.6.6/src/kernel/toppush.c, line: 1336

Fatal error:
Atomtype 2C not found
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
みたいなエラーを喰らいます。


もう1つ、イオンのatom typeに互換性がないので、名前を統一しておきます。
1c08.topを再び開き、一番下の方に存在する2つの行をそれぞれ変更します。
[ atoms ]
; id_ at type res nr residu name at name cg nr charge mass
1 IPNa+ 1 NA+ NA+ 1 1 22.9898

[ atoms ]
; id_ at type res nr residu name at name cg nr charge mass
1 IMCl- 1 CL- CL- 1 -1 35.45300




さて、これでようやく最小化に移ります。

5. エネルギー最小化

やりますのは200 stepの重原子以外のエネルギー最小化、そののちに200 stepの全体のエネルギー最小化です。

別のディレクトリを用意します。
cd ~/Desktop
mkdir minimize
cd minimize
ここに
define       = -DPOSRES1000
;
integrator = steep
nsteps = 200
;
nstlog = 1
;
nstlist = 100
ns_type = grid
pbc = xyz
rlist = 1.0
;
coulombtype = PME
rcoulomb = 1.0
;
vdwtype = cut-off
rvdw = 1.0
;
optimize_fft = yes
;
constraints = none
というmin1.mdpと、
define       =
;
integrator = steep
nsteps = 200
;
nstlog = 1
;
nstlist = 100
ns_type = grid
pbc = xyz
rlist = 1.0
;
coulombtype = PME
rcoulomb = 1.0
;
vdwtype = cut-off
rvdw = 1.0
;
optimize_fft = yes
;
constraints = none
というmin2.mdpファイルを用意します。
min1.mdpでは200 stepの重原子以外のエネルギー最小化を、min2.mdp200 stepの全体のエネルギー最小化を行います。

これを入力したら、いざ
gmx grompp -f min1.mdp -o min1.tpr -p ../grotop/1c08.top -c ../grotop/1c08.gro -po md1.out.mdp
と入力してみます。

1c08.top, 1c08.groの設定がうまく行っていましたら、ここでエラーが出ずにmin1.tprファイルが生成されているはずです。エラーが出ているようでしたら、それに応じて該当の箇所を確認してみましょう。ここはよくエラーが登場する処理です。逆に、これがうまく行っていたら残りは簡単です。

ここでsegmentation faultと出てエラーが出ずに終了するパターンも有ります。その場合は前回の記事を参照してください。

これがうまくできましたら、もう後は簡単です。
#!/bin/sh

rm -f min1.mdp.out min1.tpr
gmx grompp \
-f min1.mdp \
-po min1.out.mdp \
-c ../grotop/1c08.gro \
-r ../grotop/1c08.gro \
-p ../grotop/1c08.top \
-o min1.tpr
rm -f min1.log min1.gro min1.trr min1.edr
gmx mdrun -v \
-s min1.tpr \
-o min1.trr \
-c min1.gro \
-e min1.edr \
-g min1.log
rm -f min2.mdp.out min2.tpr
gmx grompp \
-f min2.mdp \
-po min2.out.mdp \
-c min1.gro \
-r ../grotop/1c08.gro \
-p ../grotop/1c08.top \
-o min2.tpr
rm -f min2.log min2.gro min2.trr min2.edr
gmx mdrun -v \
-s min2.tpr \
-o min2.trr \
-c min2.gro \
-e min2.edr \
-g min2.log
というrun.shファイルを作り、これを実行するだけでエネルギー最小化ができるはずです。お疲れ様でした。

これが終われば、次はGromacsを使ったNPT, NVT平衡化に取り掛かります。とはいってもここから先はきちんとした計算機・クラスタコンピュータが欲しいですね。