How to GROMACS(Last update:20070124)

Masaomi Hatakeyama 20041005 the beginning to write

* Contents
* Chap0 はじめに
* Chap1 心の準備
* Chap2 GROMACSインストール
* Chap3 GROMACS概要
* Chap4 GROMACS詳細その1 demo
* Appendix GIFアニメーションの作り方
* Appendix .xvg(xmgrace用)ファイルをgnuplot用ファイルに変換するスクリプト
* Appendix g_rotacf
* Appendix Put in the wall
* Appendix トラジェクトリーファイル(traj.trr)を連結する。
* Appendix 途中でプログラムが止まった時は、、、
* Reference

---
Chap0 はじめに

MDとはMini Disk。ではなくて、Molecular Dynamicsの略で日本語では「分子動力学」
と訳される。私はほとんどMDに関しての知識がない。恥ずかしながら統計力学熱力学
も勉強中である。

とある理由でGROMACSというMD計算ソフトウェアを使わなければならなくなった。
Googleで検索してもまだ日本のユーザーでそれほど広まっていないのか日本語の
ガイドらしきものは見当たらない。そもそも学術関連のソフトウェアは日本で
作られたもの以外英語のものがほとんどのようだが。唯一YONEYAさんの講義のページ
(http://www.nanolc.jst.go.jp/~yoneya/funcPhysChem/MDaction.html)に若干
使い方が載っているのが見つかる程度。

GROMACSの本家サイト(http://www.gromacs.org/)に行くとドキュメント関連
(チュートリアルからマニュアルまで)は結構きっちりそろっている。ただし全部英語。
英語をすらすら読めるならあまり苦にならないのだろうが、私にとってはこれが結構
高く壁としてそびえたつ。

一通りチュートリアルにも目を通して、概要はつかめたもののMDの知識が乏しい私に
はぜんぜん使いこなすことができない。全然動いてくれない。動いてもそれが正しい
のか悪いのか判断すらできない。数学の勉強と同じで、こういうときは分かるところ
まで立ち戻り、一つずつゆっくり解決していくしかない。私の人生経験上遠回りなよ
うでそれが一番早い。ここでは私のGROMACS&MD勉強の苦悩をまとめようと思う。

当然MDの基礎知識がないとGROMACSも使いこなせない(はず)。MDを勉強していくと
おのずと物理(少なくともニュートン力学)と数値計算、さらに統計力学、熱力学の
知識が最低限必要になってくる(はず)。必要性を感じたところで、横道にそれて
一つ一つ打破していく。熱・統計力学を勉強していくと当然それに付随する
数学(いわゆる物理数学。確率統計、微分積分、線形代数もろもろ)の知識が必要と
される。。。。どこまで降りればよいかはわたしの脳みそしだいなのだが、
あがってこれるところまで降りるしかない。。。(涙)

----
Chap1 心の準備

MD用のソフトウェアとしては、他に AMBER や CHARMM といったソフトウェアが有名
(なよう)。OCTA(http://octa.jp/index_jp.html)というのもある。商用では富士通
からMaterials Explorer(http://software.fujitsu.com/jp/winmasphyc/)というの
も出ている。時おり無料講習会をしている。体験版は無料ダウンロードできる。
細かい違いは今のところよくわかっていない。ちなみにMD-GRAPEというのは、理研で
開発されたMD専用の計算機(ボード、ハードウェア)のことで、簡単にいうとMDでよく
使う計算部分をハードウェアに組み込んでしまっている、というもの。通常の
コンピュータの10倍から数100倍の計算能力があるそう。
(参考http://www.peta.co.jp/index_jp.html)

MD用のソフトウェアとしてはGROMACSは後発のようで、参考となる日本語サイトが少な
い。メーリングリストも海外のものしか今のところないようだ。(とりあえず加入して
おいた。参考http://www.gromacs.org/mailing_lists/index.php)
メールは一日に平均数十通届いており、結構活発に議論されている。

GROMACSは並列計算を前提として作られているらしく、この星で最速のMDコードとなる
べく全てにおいて最適化されているらしい。
(参考http://latitude.jp/fah/gromacs.html)また、学術ソフトの中では唯一フリーの
ソフトウエアらしくソースコードもダウンロードして見ることができる。(ぜひとも
そのうちソースコードも読みたい。)(と思うだけでたぶん時間がないといってしない
とも思う。。。)

Googleで「MD」で検索すると古明地さんの「PEACH」と呼ばれるMDソフトウェアの
ページ(http://staff.aist.go.jp/y-komeiji/index_j.html)もよくヒットする。
ここには分子動力学の解説からPEACHの使い方が詳しく掲載さされている。


----
Chap2 GROMACSインストール

* root権限を持っている場合のVine3.2へのインストール
* fftw2.1.5のインストール
* GROMACS3.2.1のインストール
* 並列計算用のインストール

20060410
新しいマシンが入ったのでインストール
OSはVine3.2。ルート権限を持っているととってもインストールは簡単。
rpmファイルを取ってくる。
PCは一台なので並列処理(MPI利用)は考えない。

fftw3-3.0.1-4.i386.rpm.
ftp://ftp.gromacs.org/pub/prerequisite_software/fftw3-3.0.1-4.i386.rpm

gromacs-3.3.1-1.i386.rpm
ftp://ftp.gromacs.org/pub/gromacs/rpm/gromacs-3.3.1-1.i386.rpm

Linuxにログイン。suコマンドでrootになって

rpm -ivh fftw3-3.0.1-4.i386.rpm
rpm -ivh gromacs-3.3.1-1.i386.rpm

とこのふたコマンドでおしまい。

ソースコードからインストールするときは以下を参考に(でもVersion3.2.1のときの話なのでちょっと変わっていると思う)。

----
ここを参考に。
http://www.gromacs.org/installation/index.php

GROMACSをインストールする場合、どうも並列計算用の場合と、
そうでない場合とで若干インストール方法が異なるようだ。

わたしの場合、並列計算はよく分からないので普通にインストール。

まず、

fftw2.1.5

というのをインストールする必要がある。
Linuxのroot権限を持っている場合はrpmを同じページから
ダウンロードしてきて

rpm -ivh hogehoge

とすれば、インストールできるが、わたしはたまたまroot権限がないところに
インストールする必要があったため、ソースコードを

http://www.fftw.org/download.html

からダウンロードしてきた。GROMACSのページを見ると、


fftw-2.1.3-2.i386.rpm.

2.1系のfftwを使っているようなので、あえて最新バージョンの3.0系ではなく、

fftw-2.1.5.tar.gz

をダウンロードしてきた。
ライブラリのバージョン違いで動かなかったらやなので。

ダウンロードしてきたfftw-2.1.5.tar.gzを解凍して、

INSTALLファイル

を見てみる。

自分のマシンでない限り、たいていはroot権限は持っていないはずなので、
ここでは、一般ユーザー用の設定でインストールを進める。
一般ユーザー権限しかない場合には、たいていホームディレクトリ下にインストールする。

インストールページ
http://www.gromacs.org/installation/prerequisites.php

には

./configure --enable-float --enable-type-prefix

としろと書いてある。ここでは、

./configure --prefix=/home/masa --enable-float --enable-type-prefix

として実行した。

gromacs-3.2.1.tar.gzの中のINSTALLファイルによると

----
Read the FFTW installation instructions for details. In short, to
install the single precision library under /usr/local type

./configure --enable-float

and then type

make
make install
Note that in contrast to GROMACS, FFTW defaults to double. Even
if you don't think you'll need it's a good idea to install the
double precision libraries too, once and for all. Clean your
build by issuing 

make distclean

and then type

./configure --enable-type-prefix
make
make install

Your double precision FFTW files will have a "d" prefix.

(It is possible to compile a crippled GROMACS without FFTW, but we
strongly discourage it - check the configure options for details)
----
と書いてある。
要約すると

./configure --enable-float
make
make install
make distclean
./configure --enable-type-prefix
make
make install

と実行しろ。と。

なので、つづいて、

make
make install

とすると勝手にlibディレクトリとincludeディレクトリがホームディレクトリ下
に作成され、そこにfftwのライブラリがインストールされる。

INSTALLファイルには、さらに、

----
* FFTW OR OTHER LIBRARIES IN NON-STANDARD LOCATIONS:

If you install FFTW in your homedirectory or some other place where
it isn't found automatically (not all systems search /usr/local)
by the compiler you should set the environment variables before 
executing configure. Assume we configured and installed FFTW with 
--prefix=/home/erik/fftw. If your shell is tcsh, you set

setenv CPPFLAGS -I/home/erik/fftw/include
setenv LDFLAGS  -L/home/erik/fftw/lib

or, if you are using a bash shell:

export CPPFLAGS=-I/home/erik/fftw/include
export LDFLAGS=-L/home/erik/fftw/lib
----

と書いてある。
つまり、
FFTW用の環境変数を設定しろ、と。
なので

このライブラリをGROMACSから使用できるように環境変数を設定する。
使用しているのはbashなので、.bashrcに以下の記述を加えておく。

export CPPFLAGS=-I/home/masa/usr/local/include
export LDFLAGS=-L/home/masa/usr/local/lib

環境変数を有効にするため、

source .bashrc

と実行しておく。

続いて。
GROMACS本体をインストール。

http://www.gromacs.org/download/gromacs_source.php

から

gromacs-3.2.1.tar.gz

をダウンロードしてきて、解凍。

INSTALLファイル

を読んでみる。

----
   By default, `make install' will install the package's files in
`/usr/local/bin', `/usr/local/man', etc.  You can specify an
installation prefix other than `/usr/local' by giving `configure' the
option `--prefix=PATH'.
----
と書いてあるので、

GROMACSをインストールするディレクトリを作成しておく。

mkdir gromacs

つづいてconfigure。

./configure --prefix=/home/masa/gromacs

で、

make
make install

これでインストールは終了。
ものの5分くらいで完了する。

あとChap3以降の例題を実行するために、

/home/masa/gromacs/share/tutor

ディレクトリをホームディレクトリのどっかにまるごとコピーしてきておく。

cp -r /home/masa/gromacs/share/tutor ~


補足
Fedora Core 2にインストールしたところ、
ngmxがどういうわけかインストールされなかった。

Vine2.9、RedHat9には上手くインストールできた。

あとで出てくるが、Gromacsはデータ解析の出力を

xmgrace(グラフプロットツール)

用に出力するので、xmgraceも合わせてインストールしておくと良い。
自分のノートのVineにはインストールしておいた。
でもどうやってインストールしたかはもう記憶から削除されているらしい。。。(20041130)



20070124 並列計算用のインストール(MPI使用)

ここJAISTには並列用計算機がいくつかある。
Gromacsはシングルプロセスでも速いのだが、並列計算用にコンパイルオプションがあり、
並列計算を使わない手は無い。MPIを使用するが並列計算およびMPIについては別途勉強したが、
ここで詳細は記述しない。MPIって？＿？という人は各自勉強してみよう。
(ここではインストールだけなのでMPI用のコードを記述する知識は要らない。)

並列計算機のうちの一つ SGI Altix3700 がある。共有メモリ型128CPU。
OSにLinuxが入っていて普段からLinuxをつかっているわたしにとっては使い勝手が非常に良い。

恒例のようにまずfftwとgromacsをコンパイルしなければならない。
fftwも並列計算用にコンパイルしなおす。
今現在の最新バージョンはfftw　gromacsは3.3.1だがいままでのプロジェクトを引き継ぎたいので、
あえて低いバージョンのgromacsをいれた。
(*ちなみにgromacsはシミュレーションの途中でヴァージョンを変えてはいけない。というようなことを
  メーリングリストで以前開発者のMarkさんかDavidさんが言っていたような気がする。)

ということで。fftw2.1.5とgromacs3.2.1のmpi版インストール手順。

1.まずfftw2.1.5.tar.gzとgromacs3.2.1.tar.gzをダウンロード
2.両方解凍
3.まずfftwからインストール
　のまえに環境変数の確認
　デフォルトでAltixはtcshになるようなので、.cshrcに以下のように記述。
　このように書けとWebおよびINSTALLファイルに書かれている。
　http://www.gromacs.org/documentation/reference_3.3/gmxfaq.html#noMPIwrapper

　setenv CPPFLAGS -I/home/hatake/include
　setenv LDFLAGS -L/home/hatake/lib
　setenv LIBS -lmpi
　setenv MPICC cc

このオプション指定でわかるように、以前インストールしたときにホームディレクトリ直下に
includeディレクトリとlibディレクトリが置かれてしまっている。
mpi用のコンパイルオプションをつけるとちゃんとライブラリ名も変えてくれるので、ディレクトリが
一緒でも大丈夫な様子。でも念のため今あるライブラリはバックアップしておく。

　./configure --enable-type-prefix --enable-mpi --prefix=/home/hatake
　make
　make install

ちなみにfftwはデフォルトでdouble precision(倍精度64bit？)でコンパイルされるので--enable-float
オプションはいらない(と思う)。Altixは64bitOSなので。ちなみに通常のPCはたいてい32bitだとおもう。
ふつうのWindowsXPとか。
ポイントは--enable-mpiをつけるところ。
single precisionにするときは--enable-floatというオプションをつける。
このようにするとホームディレクトリ(/home/hatake)の下のincludeディレクトリとlibディレクトリの中に
fftw用のライブラリが格納される。

4.で次にgromacsのインストール

　./configure --prefix=/home/hatake/bin/gromacs_mpi --disable-float --enable-mpi --enable-fortran
　make mdrun
　make install-mdrun

gromacsはデフォルトでsingle precisionでコンパイルするので64bitにするには--disable-floatとという
オプションが必要。こちらもポイントは--enable-mpiをつけるとこ。
念のためFortranが使えることを教えておく→　--enable-fortran。

で、なぜか、

　make
　make install

とすると途中でエラーになった。ので、mdrunコマンドだけmpi用にコンパイルする。これでOK。みたい。
トラジェクトリーファイル(.trr)やエネルギーファイル(.edr)はmpi用で無いものと同じ形式で
出力されるので、問題なし。(みたい)。
うまくいくと/home/hatake/bin/gromacs_mpi下のまたいくつか下にmdrun実行ファイルが作成される。
mvでファイル名を書き換えてもいいが、./configure時に--program-prefix=_mpiとかすると、
出力したファイルに_mpiをつけてくれる。
INSTALLファイルには--disable-niceというオプションも書かれているが、これを書いたら
なぜかうまくいかなかった。

5.で無事mdrunコマンドが作成されたら、コマンドパスを設定して、実行する。
　mpi用のmdrunを実行するときはちょっといくつか注意点がある。

gromppしてシミュレーション用ファイルを作るときに　-np CPUノード数　というオプションが必要。
並列して有効なのは大体8個程度までで、16個、32個とCPUを増やして計算速度が倍、倍、となるかというと
そういうわけではないので、たいてい4個程度にしておく。(あまりにCPUを占有しても他に迷惑になるので･･･；)
8個も16個も大差は無い(らしい。(アムダールさん曰く))
でもって

(Altixの場合は)次にtopコマンドを使ってどのCPUがあいているか(idle状態か)を調べる。
具体的にAltixの場合は、top b n 1 とするとワンショットの全CPU情報が表示される。
(他の人がその空きCPUを使う前にすかさず)

　mpirun -np 4 dplace -s1 -c14-18 mdrun_mpi

とかする。-c14-18というのはCPU番号14から18番まで使います。というオプション。
dplaceはAltix特有のコマンドでCPU番号を指定するコマンド。Altixの場合はこれをするほうが早い。
-s1はMPIを使うときのお決まりオプション(らしい･･･)。
通常は(おそらく)

　mpirun -np 4 mdrun_mpi

だけでよいのだろう(と思う･･･；)

ちなみにsingle processで計算した場合と4CPU並列の場合とでは(わたしのブツの場合)

report
(4 cpu)
               NODE (s)   Real (s)      (%)
       Time:   1191.000   1191.000    100.0
                       19:51
               (Mnbf/s)   (MFlops) (ps/NODE hour) (NODE hour/ns)
Performance:     16.594    693.783    302.267      3.308

(1 cpu)
               NODE (s)   Real (s)      (%)
       Time:   3949.877   3965.000     99.6
                       1h05:49
               (Mnbf/s)   (MFlops) (ps/NODE hour) (NODE hour/ns)
Performance:      5.002    208.725     91.142     10.972


大体3倍程度の計算速度の差が確認できた。(計算するブツにもよるだろうが･･･)
微妙に粒子の配置場所は違うが、統計的には同じになるのだろう･･･(と思う･･･；)


----
Chap3 GROMACS概要

* GROMACSはコマンドラインツール
* 各設定ファイル
* 実行の流れ


まず！

* GROMACSはコマンドラインツールであること

GUIではない。コマンドとオプションを覚えて、カタカタタイプして使う。
Windows用のバイナリ(実行ファイル)も用意されているがやはりいろいろな理由から
UNIX(Linux)系のOSで使いたいところ。コマンドの数とそのコマンドに付随するオプションは
数多い。全部を覚えるのはまず無理なので、まずはじめは必要最低限のものを覚えることから
はじめるのが良いと思う。
(参考オンラインリファレンスhttp://www.gromacs.org/documentation/reference_3.2/online.html)

次に！

* 各設定ファイルがいろいろある

GROMACSで使われる代表的なファイルは以下の6つ。(括弧内はファイルの拡張子)

PDBファイル(.pdb)
 Protein Data Bankのフォーマットファイル。
 PDBとはここ->http://pdb.protein.osaka-u.ac.jp/pdb/
 世界中の研究者がたんぱく質の構造データを登録している。
 (くわしことは私も知らない。)

Topologyファイル(.top)
 そのシステムにはどんな分子がいくつかるのかなどの全体的情報。　
 テキストファイル。Vi、Emacs、メモ帳などで閲覧編集可。能

Structureファイル(.gro)
 具体的にどんな分子がどの座標にあるのかというのが書かれている。
 テキストファイル。Vi、Emacs、メモ帳などで閲覧編集可能。

Parameterファイル(.mdp)
 シミュレーションを走らせる時のさまざまなパラメータ(設定値)が書かれている。
 ◎かなり重要!!
 はじめのうちはどのパラメータをどういじったらよいのかわからない。
 テキストファイル。Vi、Emacs、メモ帳などで閲覧編集可能。

Run Inputファイル(.tpr)
 シミュレーション実行に必要な唯一の設定ファイル。(大まかに言うと)上記の
 .top, .gro, .mdp を一つにまとめたファイル。
 シミュレーション実行前に作る。
 バイナリファイル。Vi、Emacs、メモ帳などのテキストエディタでは見れない。

Trajectoryファイル(.trr)
 シミュレーション実行後にできる結果のファイル。これを解析する。
 各分子がどの時間にどの位置にいるかというのが書かれている。
 これもバイナリファイル。
 各種解析用のコマンドが用意されている。

(参考チュートリアルhttp://www.gromacs.org/documentation/reference_3.2/online/getting_started.html#start)

でもって!!

* コマンド実行の流れ

大まかには次のように実行していく。
(参考フローチャートhttp://www.gromacs.org/documentation/reference_3.2/online/flow.html)

1. 構造ファイル(structure file)(.gro)とトポロジーファイル(topology file)(.top)を
   PDBファイル(.pdb)から作る。(pdb2gmxコマンド使用)

   .pdb ----> (pdb2gmx) ----> .gro, .top

2. ボックスの大きさを変更。 (editconfコマンド使用)

   .gro ----> (editconf) ----> .gro

3. 水分子を混ぜる(ペプチドを溶かす)。(genboxコマンド使用)

   .gro, .top ----> (genbox) ----> .gro, .top

4. パラメータファイル(.mdp)を手動で作成。

   ◎ここが肝。

5. シミュレーション実行用の設定ファイル(.tpr)を構造ファイル(.gro)、トポロジー
   ファイル(.top)、パラメータファイル(.mdp)から作る。(gromppコマンド使用)

   .gro, .top, .mdp ----> (grompp) ----> .tpr

6. シミュレーション実行(mdrunコマンド使用)

   .tpr ---> (mdrun) ---> .trr, .gro

7. .trrファイル解析(g_XXXコマンド)

   ◎ここも肝。

という流れ。実際はシミュレーションを行う前に、Energy Minimization(エネルギー
最小化計算)というのとPosition Restrainsといって短めなMD計算を行う。

マニュアルにはこのように書いてある。
(http://www.gromacs.org/documentation/reference_3.2/online/speptide.html)

1. Convert the pdb-file speptide.pdb  to a GROMACS structure file and a GROMACS topology file.
2. Solvate the peptide in water
3. Perform an energy minimization of the peptide in solvent
4. Add ions if necessary (we will omit this step here)
5. Perform a short MD run with position restraints on the peptide
6. Perform full MD without restraints
7. Analysis


なぜこれら(EM, PM)が必要なのかはいまのところ(私自身が)よく分かっていない。
demoプログラムの解説によると、

Energy Minimizationについて
In principle we can start a Molecular Dynamics simulation now. However
it is not very wise to do so, because our system is full of close
contacts. These close contacts are mainly a result of the genbox
program. The added solvent might have some close contacts with the
peptide resulting in very high repulsive energies. If we would start a
Molecular Dynamics (MD) simulation without energy minimisation the
system would not be stable because of these high energies.

The standard procedure to remove these close contacts is
Energy Minimisation (EM). Energy minimisation slightly changes the
coordinates of our system to remove high energies from our system.

要約すると、
genboxコマンドを使ってペプチドを溶かした(水分子を放り込んだ)状態では、
分子同士が極端に近いものがあったりしてエネルギーがものすごく高い状態で
計算がとまってしまう場合がある。
とのこと。それを解決して分子を適当にばらけさせるのがEnergy Minimization
ということらしい。この計算アルゴリズムもさまざまなものがあるだろうが、
それについてはまた別の機会にしよう。とりあえず、ここまでの理解で進めてゆく。

次に
Position Restrained MDについて
Position Restrained MD keeps the peptide fixed and lets all water
molecules equilibrate around the peptide in order to fill holes
etc. which were not filled by the genbox program.

要約すると、
ペプチドの周りに均等に水分子を配置する、処理らしい。
米谷さんの解説によると(http://www.nanolc.jst.go.jp/~yoneya/molsim/MDflow.html)
エネルギ最小化された構造には原子の速度の情報が無いので、適当な温度を設定し、
それに対応したボルツマン分布で原子速度を適当に設定してスタートアップMDをまず
行う。 
とのことらしい。とりあえず、いまのところは概要理解なのでここま進める。


結果の解析はさまざまであろうが、もっとも初期の基本的な操作は、分子の動きを
アニメーションで見ることだろう。そのためのコマンドとして、

ngmx

というコマンドが用意されている。トラジェクトリファイル(.trr)解析用の
コマンドはたいてい g_XXXX というコマンド名になっている。

とりあえず、概要はここまで。大体流れがつかめたと思うので(え？つかめられない？)
次は詳細なコマンド操作とファイル構造について眺めていくことにする。

----
Chap4 GROMACS詳細その1

* demoプログラム

今現在の知識の状態でできるところまで調べながら進める。
題材はチュートリアル
(http://www.gromacs.org/documentation/reference_3.2/online/getting_started.html)
とあわせて現在わたしが取り組み中のCoarse Grain ModelのPEO-PPO-PEOシステムに
ついて、比較を含めながら進める。

Gorase Grain Modelとは簡単に言うと、一つまたは複数の分子を一つの原子と考えて
高分子を粗く抽象化した分子モデルのこと。
こうすることで高分子のシミュレーションを高速に行うことができる。
(その反面精度はさがる)

PPOとは Polypropyleneoxide のことで PEO とは polyethyleneoxide のこと。

ポリエチレンオキサイドの化学式は(C2H4O)nで、
ポリプロピレンオキサイドの化学式は(C3H6O)n。

より具体的に構造を書くと
PEO HO-[-CH2-CH2-O-]n
PPO [-CH-CH2-O-]m
     |
     CH3
PEO [-CH2-CH2-O-]n -H

とつながっている物質。EOの部分が親水性(hydrophilic)で、POの部分が
疎水性(hydorophobic)になっている。非イオン性の界面活性剤(surface-active agents)。
親水性の部分と疎水性の部分をもった分子を両親媒性の分子という。
両親媒性の分子は水溶液中で自己集合(self-assembled)して、膜を形成する。

界面活性剤とは、洗剤の原料にもなっているので有名だが、細胞膜の原料も
界面活性剤である。界面活性剤にはいろいろと種類があり、大きく分類すると、
陽イオン系(カチオンcation)、陰イオン系(アニオンanion)、非イオン系(ノニオン
nonion)とに分類される。人工のもの、天然のものと多数存在する。

はじめ教授と話をしたとき(こちらに来て1週間くらいのとき)何度ボイスレコーダーを
聞いても「POPPO」としか聞こえず、ラボの人に聞いたところPPOとPEOでないかと
教えてくれた。
they are used in making triblock copolymers and have differing solubilities
(as far as hydrophobic and hydrophilic properties are concerned)

恥ずかしながら化学の知識も乏しいので、これまた補いながら進める。
(物理、数学、化学、英語ほぼ全科目!!まさに!!ガチンコ留学生活!!
 唯一の救いはコンピュータがある程度つかえる、ということだけか。。。)
チュートリアルは全部で6つ用意されており、いろいろな物質を解説にしたがって
シミュレーションできるようになっている。設定ファイル等も用意されているので、
各ファイルを比較しつつ勉強しつつ自分のシミュレーションに応用していくのが
良いだろう。

* demoプログラム
* Water分子
* Methanol分子
* Water/Methanol混合
* Ribonuclease S-Peptide
* Protein unfolding

とある。ひとつひとつじっくり見ていこう。
まずは、

----
* demoプログラム

これはCシェルスクリプトで書かれており、コマンドと解説文が書かれている。
これを一通り読むだけでも概要がつかめるが、スクリプトファイルをHackして
コマンドの使い方 & ファイルの構造を学んでみよう！

チュートリアルをみると、
(http://www.gromacs.org/documentation/reference_3.2/online/getting_started.html#start)
Before starting the examples, you have to copy all the neccesary files, 
to your own directory.

tutorディレクトリをまるごとホームディレクトリーにコピーして来い、と書いてある。
わたしの環境では、/usr/local/share/gromacs/tutor/ディレクトリにファイルが
あったのでそれをつかう。

最初のdemoプログラムはgmxdemoディレクトリに入っているので覗いてみる。
gmxdemoディレクトリにはつぎの2つのファイルがある。

cpeptide.pdb	ペプチドのPDBファイル
demo		Cスクリプトファイル

もしCシェルが入っていれば(たいてい入っている)、

./demo

とタイプすれば解説とあわせてシミュレーション(+準備)が始まる。

一つ一つ解説を追いつつ、コマンドと使用(出力)ファイルをじっくり眺めていこう。

まずdemoスクリプトファイルを実行する前に、
cpetide.pdbファイルから。こんな↓ファイルになっている。

ATOM      1  N   LYS     1      24.966  -0.646  22.314  1.00 32.74      1SRN  99
ATOM      2  CA  LYS     1      24.121   0.549  22.271  1.00 32.05      1SRN 100
ATOM      3  C   LYS     1      24.794   1.733  22.943  1.00 31.16      1SRN 101
ATOM      4  O   LYS     1      25.742   1.575  23.764  1.00 31.50      1SRN 102
ATOM      5  CB  LYS     1      22.812   0.323  23.047  1.00 33.09      1SRN 103
ATOM      6  CG  LYS     1      21.763   1.415  22.695  1.00 34.29      1SRN 104
ATOM      7  CD  LYS     1      20.497   1.124  23.561  1.00 34.93      1SRN 105
ATOM      8  CE  LYS     1      20.706   1.659  24.970  1.00 35.35      1SRN 106
ATOM      9  NZ  LYS     1      21.524   0.759  25.825  1.00 35.85      1SRN 107
ATOM     10  N   GLU     2      24.300   2.909  22.632  1.00 29.30      1SRN 108
ATOM     11  CA  GLU     2      24.858   4.145  23.207  1.00 27.38      1SRN 109
ATOM     12  C   GLU     2      24.567   4.201  24.693  1.00 26.12      1SRN 110
ATOM     13  O   GLU     2      23.398   4.051  25.038  1.00 26.39      1SRN 111


PDBのサイト(http://www.protein.osaka-u.ac.jp/pdb/info/faq.html)によると、
PDB ファイルは共通して 
1 桁目から 9 桁目までに項目名、
11 桁目から 72 桁目までにそれについての記述がある。
項目の内容が 2 行以上になる場合は、10 行目に番号が記入される。
73 桁目から 76 桁目までは登録番号であるが、登録番号が 5 桁になる場合は 
77 桁目にも記されている。
77 桁目から 80 桁目までは行番号である。

らしい。(細かいことはすっ飛ばして)
重要なのは先頭が「ATOM」と書いてある行で(これがほとんど)、各原子の配置を
記述している。
詳しくは、
http://www.protein.osaka-u.ac.jp/pdb/info/faq.html
http://www.dna.affrc.go.jp/misc/mm/pdb/
を参考にしてもらうとして、ATOM行で重要なのは、

* 2番目の列に通し番号、
* 3番目の列に原子名、
* 4番目の列に残基名、
* 5番目の列が残基番号、
その後に続く、6番目、7番目、8番目の数値がその原子のそれぞれ
* x座標、y座標、z座標

になる。単位はオングストローム。
とりあえず、これだけしっとけばいいでしょう。いまのところ。たぶん。。。

でもって。
このPDBファイルを3Dで見れるソフトがいくつかある。
MOLMOL(http://129.132.45.141/wuthrich/software/molmol/index.html)
rasmol(http://www.scl.kyoto-u.ac.jp/scl/appli/rasmol.html)

他にもあるようだがどれが良いかなど、調べつつ随時掲載していくとして、
(といいつつもたぶんしないと思うけど)先を進める。インストール方法なども
随時掲載していく(つもり)。(といいつつ絶対しないと思う。。。)

windows版rasmol.exeはただダウンロードしてくるだけで即実行可能なので、
簡単。今回はWindows版を用いて、[File]-[Open]してみて、ペプチドをぐるぐる
回してどんなもんかとりあえず確認してみた。

----
でこれからの手順としては、Chap3 GROMACS概要で書いたように

1. pdb2gmx(.pdbを.groに変換)
2. editconf(.groの修正(箱の大きさを変更))
3. genbox(.gro, .topの修正(箱の中に水分子を入れる))
4. .mdpの作成(手動)
5. grompp(.mdp、.gro、.topから.tprを作成)
6. mdrun(シミュレーション)

と進む。

さて。
では、demoプログラムを実行するとしよう。

./demo

----
Welcome to the GROMACS demo.
This is a script that takes about 10 min to run.
In this demo we will demonstrate how to simulate Molecular
Dynamics (MD) using the GROMACS software package.
The demo will perform a complete molecular dynamics (MD) simulation
of a small peptide in water. The only input file we need to do this
is a pdb file of a small peptide.
If you have any problems or remarks with respect to this demonstration,
please mail to: gromacs@gromacs.org , or check the resources on
our website http://www.gromacs.org .
----
(要約)
ようこそGROMACSデモ。
このデモは約10分くらいかかるよ。
以下略。

はい次。(Enterキーを押す)

----
Before we you can actually start the GROMACS demo, the programs
must be present in your PATH. This might already be the case if
they are linked to /usr/local/bin. If not, follow the instructions
in the getting started section. If GROMACS is not installed 
properly on your computer, contact your system manager.
----
(要約)
環境変数PATHのの設定を確認してねん。

はい次。

----
Before we can start any simulation we need a molecular toplogy
file. This topology file ( .top extension ) is generated by the
program pdb2gmx. The only input file of the pdb2gmx program is the pdb
file of our peptide ( .pdb extension ). 

Because most pdb files do not contain all hydrogen atoms, the pdb2gmx
program will also add them to our peptide. The output file which
contains the structure of the peptide when hydrogen atoms are added is a
gromos structure file ( .gro extension )  
----
ここからが大事そう。

(要約)
pdb2gmxコマンドによってトポロジーファイル(.top拡張子)をPDBファイル(.pdb拡張子)
から作る必要がある。

で、そのコマンドは、、、(Cスクリプトファイルを覗くと)

pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top >& ! output.pdb2gmx

と書いてある。一つずつ見ていこう。

pdb2gmxはコマンド名。これは良いとして、

次の

-fオプションの引数で

${MOL}.pdb

というのは、このスクリプトの中で定義してあるMOLというシェル変数
に.pdbというのをつけたpdbファイルということで、

先頭で、

setenv MOL cpeptide

と定義してあるので、

pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top >& ! output.pdb2gmx

は

pdb2gmx -f cpeptide.pdb -o cpeptide.gro -p cpeptide.top >& ! output.pdb2gmx

と同じことを実行している。ということになる。

コマンドオプションは全部で3つ。

-f 
-o 
-p 

でそれぞれファイルが指定してある。
マニュアルによると
(http://www.gromacs.org/documentation/reference_3.2/online/pdb2gmx.html)
-f　　  eiwit.pdb 	Input 	Generic structure: gro g96 pdb tpr tpb tpa xml
-o 	conf.gro 	Output 	Generic structure: gro g96 pdb xml
-p 	topol.top 	Output 	Topology file
とある。

つまり、
-fオプションは、PDBファイルを入力ファイルとして指定して、
-oオプションは、.groファイルを
-pオプションは、.topファイルを出力ファイルとして指定する。

ということになる。
ちなみにファイル名を指定しなかった場合、自動的に(デフォルトで)
eiwit.pdb、conf.gro、topol.topというファイル名になるようだ。

ここでは、

pdb2gmx -f cpeptide.pdb -o cpeptide.gro -p cpeptide.top >& ! output.pdb2gmx

と全部個別にファイル名を指定してあるので、このコマンドを実行した後には、

cpeptide.groファイルとcpeptide.topファイルができているはず!
最後の

>& ! output.pdb2gmx

というのはコマンドを実行したとき通常画面に出力される情報をすべて
output.pdb2gmxファイルに出力する。というシェルのリダイレクション指定。
>だけの場合は、エラー出力は出力されない。
!マークの意味はよく分からないが、たいした問題ではないと思うので気にせず、
先を進める。

で実行(Enter)してみてlsしてみると、

cpeptide.gro  cpeptide.pdb  cpeptide.top  demo*  output.pdb2gmx  posre.itp

新しく

1. cpeptide.gro
2. cpeptide.top
3. posre.itp
4. output.pdb2gmx

4つできている。(ぉっと、期待してない.itpなるファイルができている?!)

これらを一個一個見ていこう。

まず。
1. cpeptide.gro(通称GROMOSファイル)

Good gRace! Old Maple Actually Chews Slate
  143
    1LYSH     N    1   2.497  -0.065   2.231
    1LYSH    H1    2   2.449  -0.140   2.186
    1LYSH    H2    3   2.583  -0.046   2.184
    1LYSH    H3    4   2.515  -0.089   2.327
    1LYSH    CA    5   2.412   0.055   2.227
    1LYSH    CB    6   2.281   0.032   2.305
    1LYSH    CG    7   2.176   0.141   2.270
    1LYSH    CD    8   2.050   0.112   2.356
    1LYSH    CE    9   2.071   0.166   2.497
    1LYSH    NZ   10   2.152   0.076   2.582
    1LYSH   HZ1   11   2.163   0.116   2.673
    1LYSH   HZ2   12   2.107  -0.013   2.590
    1LYSH   HZ3   13   2.243   0.064   2.541
    1LYSH     C   14   2.479   0.173   2.294
    1LYSH     O   15   2.574   0.157   2.376
    2GLU      N   16   2.430   0.291   2.263
    2GLU      H   17   2.353   0.297   2.199
    2GLU     CA   18   2.486   0.414   2.321
    2GLU     CB   19   2.424   0.535   2.254

何のことはない、PDBファイルとほとんど同じ。並び順などが違うだけ。
でもじっくり見てみると、総原子数が異なる。
PDBファイルで108個、.groファイルで143個。これはH原子の数の違いだというのが
わかる。PDBファイルにはH(水素)原子がついていない。
先頭にはなにやら毎回異なるフレーズ

Good gRace! Old Maple Actually Chews Slate

が出力されるが気にしないでおこう。
2行目には総原子数が書かれている。
最後の行には、
   1.42915   2.16366   1.21831
などという数字が書かれているがこれはボックスのサイズ(x,y,z)をあらわしている。
単位はオングストロームではなくて、nm(ナノメートル)。

ちなみにオングストロームは記述は[Å]こんな
感じで、1Åが0.1nm(ナノメートル)。
せっかくなので、単位の接頭語(?)をまとめておこう。

k	キロ	10の３乗
M	メガ	10の６乗
G	ギガ	10の９乗
T	テラ	10の12乗
P	ペタ	10の15乗
E	エクタ	10の18乗
Z	ゼタ	10の21乗
Y	ヨタ	10の24乗

m	ミリ	10の−３乗
μ	マイクロ10の−６乗
n	ナノ	10の−９乗
p	ピコ	10の−12乗
f	フェムト10の−15乗
a	アト	10の−18乗
z	ゼプト	10の−21乗
y	ヨクト	10の−24乗

最近ナノテクノロジーなどといわれて久しいが、
DNAの太さがだいたい2nmくらいらしい。
ちなみに大腸菌の大きさが1〜2μm(ミクロン)くらい。
水素原子の大きさはだいたい0.25ナノメートル。
炭素原子の大きさはだいたい0.35ナノメートル。
ナノというと、原子数個〜1000個オーダーの大きさの物質
が対象になりそうだ。だいたいタンパク質とかいうとこの
オーダーになりそう。

さて、余談はこのくらいにして本筋に戻って。

続いて。
2. cpeptide.top(トポロジーファイル)

;
;       File 'cpeptide.top' was generated
;       By user: m-hatake (500)
;       On host: localhost.localdomain
;       At date: Sat Oct  2 01:06:41 2004
;
;       This is your topology file
;       "I Live the Life They Wish They Did" (Tricky)
;
; Include forcefield parameters
#include "ffG43a1.itp"

[ moleculetype ]
; Name            nrexcl
Protein             3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  typeB    chargeB      massB
     1         NL      1   LYSH      N      1      0.129    14.0067   ; qtot 0.129
     2          H      1   LYSH     H1      1      0.248      1.008   ; qtot 0.377
     3          H      1   LYSH     H2      1      0.248      1.008   ; qtot 0.625
     4          H      1   LYSH     H3      1      0.248      1.008   ; qtot 0.873
     5        CH1      1   LYSH     CA      1      0.127     13.019   ; qtot 1
     6        CH2      1   LYSH     CB      2          0     14.027   ; qtot 1
     7        CH2      1   LYSH     CG      3          0     14.027   ; qtot 1
     8        CH2      1   LYSH     CD      3          0     14.027   ; qtot 1
(中略)
   141          C     13    MET      C     57       0.27     12.011   ; qtot 2.27
   142         OM     13    MET     O1     57     -0.635    15.9994   ; qtot 1.635
   143         OM     13    MET     O2     57     -0.635    15.9994   ; qtot 1

[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
    1     2     2    gb_2
    1     3     2    gb_2
    1     4     2    gb_2
    1     5     2    gb_20
    5     6     2    gb_26
    5    14     2    gb_26
    6     7     2    gb_26
    7     8     2    gb_26
    8     9     2    gb_26
    9    10     2    gb_20
(中略)
[ pairs ]
;  ai    aj funct            c0            c1            c2            c3
    1     7     1
    1    15     1
    1    16     1
    2     6     1
    2    14     1
    3     6     1
    3    14     1
    4     6     1
    4    14     1
    5     8     1
(中略)
[ angles ]
;  ai    aj    ak funct            c0            c1            c2            c3
    2     1     3     2    ga_9
    2     1     4     2    ga_9
    2     1     5     2    ga_10
    3     1     4     2    ga_9
    3     1     5     2    ga_10
    4     1     5     2    ga_10
    1     5     6     2    ga_12
    1     5    14     2    ga_12
    6     5    14     2    ga_12
    5     6     7     2    ga_14
(中略)
[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3            c4            c5
    2     1     5    14     1    gd_14
    1     5     6     7     1    gd_17
    1     5    14    16     1    gd_20
    5     6     7     8     1    gd_17
    6     7     8     9     1    gd_17
    7     8     9    10     1    gd_17
    8     9    10    11     1    gd_14
    5    14    16    18     1    gd_4
   14    16    18    24     1    gd_19
   16    18    19    20     1    gd_17
   16    18    24    26     1    gd_20
   18    19    20    21     1    gd_17
(中略)
[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            c2            c3
    5     1    14     6     2    gi_2
   14     5    16    15     2    gi_1
   16    14    18    17     2    gi_1
   18    16    24    19     2    gi_2
   20    23    22    21     2    gi_1
   24    18    26    25     2    gi_1
   26    24    28    27     2    gi_1
   28    26    33    29     2    gi_2
   28    32    30    29     2    gi_2
(中略)
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include water topology
#include "spc.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

; Include generic topology for ions
#include "ions.itp"

[ system ]
; Name
Protein

[ molecules ]
; Compound        #mols
Protein             1

ふぅ〜。えらい長い。。。
まとめると。
どうも「;」があるのはコメント行みたいで、
[****]
というのでセクションが分かれているのがわかる。
全部のセクションを列挙してみると

[ moleculetype ]
[ atoms ]
[ bonds ]
[ pairs ]
[ angles ]
[ dihedrals ]
[ dihedrals ]
[ position_restraints ]
[ system ]
[ molecules ]

とある。topology(構造)ファイル、という名前から想像できるように、
きっとこのファイルは分子の構造を記述しているファイルであるに違いない！！
だって.groファイルは分子の種類と位置しか書いてなかったのだから。

一つ一つ焦らず調べてゆこう！！
と思ったけども、マニュアルにChap5で独立して
25ページくらいにわたって解説してあるので、
一筋縄にはいかなそうなので、ここはひとまず後回し。

<後回しList>
1.トポロジーファイルについて

次のファイル。
3. posre.itp

というファイル。チュートリアルにはどこにもないままいきなり登場。
最初のpdb2gmxコマンド実行時にできていたファイル。
マニュアルを見てみると、
(http://www.gromacs.org/documentation/reference_3.2/online/itp.html)
The itp file extension stands for include toplogy. These files are included
in topology files ( with the top extension )

どうも、Topologyファイル(.top)にincludeされるファイルのようだ。
おそるおそるファイルを覗いてみてみると、
----
; In this topology include file, you will find position restraint
; entries for all the heavy atoms in your original pdb file.
; This means that all the protons which were added by pdb2gmx are
; not restrained.

[ position_restraints ]
; atom  type      fx      fy      fz
     1     1  1000  1000  1000
     5     1  1000  1000  1000
     6     1  1000  1000  1000
     7     1  1000  1000  1000
     8     1  1000  1000  1000
     9     1  1000  1000  1000
    10     1  1000  1000  1000
    14     1  1000  1000  1000
    15     1  1000  1000  1000
    16     1  1000  1000  1000
    18     1  1000  1000  1000
(中略)
   133     1  1000  1000  1000
   134     1  1000  1000  1000
   136     1  1000  1000  1000
   137     1  1000  1000  1000
   138     1  1000  1000  1000
   139     1  1000  1000  1000
   140     1  1000  1000  1000
   141     1  1000  1000  1000
   142     1  1000  1000  1000
   143     1  1000  1000  1000

143というのはcpeptide.groファイルの原子数と同じ数。
おそらく一行に一原子の情報が記載されているだろうことが
予想できる。

で、おそらく「;」がコメントで、その内容を見ると
; In this topology include file, you will find position restraint
; entries for all the heavy atoms in your original pdb file.
; This means that all the protons which were added by pdb2gmx are
; not restrained.

よくわからないがとにかくPositionRestraint時に必要なのだろうことが
わかる。どうも喉につかえてる気分だが、このまま進める。でもって
これがcpeputide.topとともに必要なのだろうことも予想できる。

予想通りトポロジーファイルcpeputide.topを見てみると、

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

という記述があった。

で、次に行くことにする。

4. output.pdb2gmx

これは、コマンド実行の際リダイレクション(>&)で指定したおいたファイル名。
なので、結局、

pdb2gmx -f cpeptide.pdb -o cpeptide.gro -p cpeptide.top >& ! output.pdb2gmx

のコマンドを実行した際の出力に相当する。でもってこのファイルの内容は、
demo実行の際に別ウィンドウにも出力されている。

でで、中身を見てみると。
-----

                         :-)  G  R  O  M  A  C  S  (-:

                      GROup of MAchos and Cynical Suckers

                            :-)  VERSION 3.2.1  (-:
(中略)
                               :-)  pdb2gmx  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f   cpeptide.pdb  Input        Generic structure: gro g96 pdb tpr tpb tpa
                                   xml
  -o   cpeptide.gro  Output       Generic structure: gro g96 pdb xml
  -p   cpeptide.top  Output       Topology file
  -i      posre.itp  Output       Include file for topology
  -n      clean.ndx  Output, Opt. Index file
  -q      clean.pdb  Output, Opt. Generic structure: gro g96 pdb xml

(中略)
There are 22 donors and 20 acceptors
There are 23 hydrogen bonds
Will use HISB for residue 12
Making bonds...
Opening library file /usr/local/share/gromacs/top/aminoacids.dat
Number of bonds was 149, now 144
Generating angles, dihedrals and pairs...
Before cleaning: 225 pairs
Before cleaning: 258 dihedrals
There are   74 dihedrals,   66 impropers,  209 angles
           225 pairs,      144 bonds and     0 dummies
Total mass 1547.779 a.m.u.
Total charge 1.000 e
Writing topology

Writing coordinate file...

◎
Select the Force Field:
 0: GROMOS96 43a1 Forcefield (official distribution)
 1: GROMOS96 43b1 Vacuum Forcefield (official distribution)
 2: GROMOS96 43a2 Forcefield (development) (improved alkane dihedrals)
 3: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
 4: Gromacs Forcefield (see manual)
 5: gmx Forcefield with hydrogens for NMR stuff (Do NOT use for new runs)
Looking whether force field file ffG43a1.rtp exists
Reading cpeptide.pdb...
Read 108 atoms
Analyzing pdb file
There are 1 chains and 0 blocks of water and 13 residues with 108 atoms

  chain  #res #atoms
  1 ' '    13    108  

Reading residue database... (ffG43a1)
Processing chain 1 (108 atoms, 13 residues)
Checking for duplicate atoms....
deleting duplicate atom    O  MET  13  pdb nr  108
Now there are 107 atoms
N-terminus: NH3+
C-terminus: COO-
Now there are 13 residues with 143 atoms
		--------- PLEASE NOTE ------------
You have succesfully generated a topology from: cpeptide.pdb.
The select force field and the spc water model are used.
Note that the default mechanism for selecting a force fields has
changed, starting from GROMACS version 3.2.0
		--------- ETON ESAELP ------------

gcq#172: "It Wouldn't Hurt to Wipe Once In a While" (Beavis and Butthead)

----
ふぅ〜。これも長い。
はじめはGROMACSの広告とオプションの説明がずらずら表示されていて、
続いて、分子の情報をずらずら表示している。

でたぶん重要なのが、

Select the Force Field:
 0: GROMOS96 43a1 Forcefield (official distribution)
 1: GROMOS96 43b1 Vacuum Forcefield (official distribution)
 2: GROMOS96 43a2 Forcefield (development) (improved alkane dihedrals)
 3: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
 4: Gromacs Forcefield (see manual)
 5: gmx Forcefield with hydrogens for NMR stuff (Do NOT use for new runs)
Looking whether force field file ffG43a1.rtp exists

の部分で、Forcefieldに何を使うか6択の中から選ぶように指定されている。
今回のケースでは「0」が選択されているようす。これはどこからわかるかというと、

Looking whether force field file ffG43a1.rtp exists

というのからも予想できるし、
シェルスクリプトのコマンドで

echo 0 | pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top >& ! output.pdb2gmx

と先頭に「0」が入力されていることからもわかる。

ForceField(力場)と呼ばれるのは、ごくごく簡単にいうと、
分子間の相互作用を計算するときに、どういうポテンシャル関数とパラメータセット
を使うかという指定で、マニュアルにもChap4に80ページ近くの分量で書かれている。
ポテンシャル関数というのは、分子間力を計算する式のこと。計算式の中には必ず
いくつか定数が入っているのだけども、分子の種類などによって経験的にその係数が
決まっていることが多い。
ここも一筋縄では理解できなそうなので、ひとまずおいておく。

<後回しリスト>
1.Topologyファイル(.top)の構造
2.pdb2gmx実行時のForceFieldの指定


さて。やっとこさ。一つ目のコマンドが終了したので、次のコマンドにうつる。
もう一度復習。

1. pdb2gmx(.pdbを.groに変換)
2. editconf(.groの修正(箱の大きさを変更))
3. genbox(.gro, .topの修正(箱の中に水分子を入れる))
4. .mdpの作成(手動)
5. grompp(.mdp、.gro、.topから.tprを作成)
6. mdrun(シミュレーション)

全工程のうち、「1」を終了したので、次は「2」。

demoを進める。(Enterキー)

----
Because a simulation of a peptide in vacua is a bit unrealistic, we
have to solvate our peptide in a box of water. genbox is the program
we use to do this.

The genbox program reads the peptide structure file and an input file
containing the sizes of the desired water box. The output of genbox is
a gromos structure file of a peptide solvated in a box of water. The
genbox program also changes the topology file ( .top extension ) to
include water. First we will use the program editconf to define the
right boxsize for our system.
----
(要約)
真空中でのシミュレーションは非現実的なので、箱の中に水分子を放り込んで
ペプチドを溶かす。そのためにgenboxというコマンドを使う、と。
genboxコマンドは.groファイルと.topファイルを変更する、と。
で、genboxコマンドを使う前にeditconfコマンドで箱のサイズを変える。

ということで、次はeditoconfコマンドが実行されることが予想される。
実はgenboxコマンドにも箱のサイズを変更するオプションが用意されているようだが、
それは後ほど。
ということでここではdemo通りに進めていく。

(Enterキー)

でなにやら実行されたので、コマンドを除いてみる。

editconf -f ${MOL}.gro -o ${MOL}.gro -d 0.5 >& ! output.genbox
genbox -cp ${MOL}.gro -cs -o ${MOL}_b4em.gro -p ${MOL}.top >>& ! output.genbox

ぉっと、2つも同時に実行されてしまった。
シェル変数を元に戻して

editconf -f cpeputide.gro -o cpeputide.gro -d 0.5 >& ! output.genbox
genbox -cp cpeputide.gro -cs -o cpeputide_b4em.gro -p cpeputide.top >>& ! output.genbox

一つ一つ見てゆこう！

まずeditconfから。

editconf -f cpeputide.gro -o cpeputide.gro -d 0.5 >& ! output.genbox

オプションは全部で三つ。

-f
-o
-d

-fと-oは両方とも.groファイルが指定されているが、-dオプションは、実数値が
指定されている。あとは前回同様リダイレクションで出力情報をファイルに
保存している。

マニュアルを除いてみよう!
(http://www.gromacs.org/documentation/reference_3.2/online/editconf.html)
解説が長い。。。

重要な部分は、、、
editconf converts generic structure format to .gro, .g96 or .pdb.
The box can be modified with options -box, -d and -angles. Both -box and -d
will center the system in the box. 
(要約)
editconfは.gro、.g96、.pdbファイルを変換するコマンドで、
-box、-d、-anglesオプションでボックスを修正できる。らしい。

Option -bt determines the box type: tric is a triclinic box, cubic is a cubic
 box, dodecahedron is a rhombic dodecahedron and octahedron is a truncated
 octahedron. 
(要約)
どうも、オプション-btによってボックスの形状を指定できるよう。

ここでは、立方体(cubic)を前提として話しているのでとくに-btは関係なさげ。

オプション説明のところを見ると、

-f	conf.gro 	Input 	 Generic structure: gro g96 pdb tpr tpb tpa xml
-o 	out.gro 	Output 	 Generic structure: gro g96 pdb xml
-d 	real 	  	Distance between the solute and the box
とあるので、

-f	入力として.groファイルを
-o	変換後の出力として.groファイル名を指定する。

ようだ。
で、

-d	ボックスと溶質(ペプチド)との距離

ここもおそらく単位はnm(ナノメートル)？なのかな？

コマンドでは -d 0.5 としてしてあった。

でもって、出力ファイルoutput.genboxをのぞいてみよう。
先ほどのpdb2gmxの場合と同様にはじめにGROMACSの広告と
オプションの使い方がずらずら表示された後、

#Entries in atommass.dat: 82 vdwradii.dat: 26 dgsolv.dat: 7
Read 143 atoms
Volume: 3.76725 nm^3, corresponds to roughly 1600 electrons
No velocities found
    system size :  1.429  2.163  1.219 (nm)
    center      :  2.469  1.077  2.517 (nm)
    box vectors :  1.429  2.164  1.218 (nm)
    box angles  :  90.00  90.00  90.00 (degrees)
    box volume  :   3.77               (nm^3)
    shift       : -1.254  0.504 -1.407 (nm)
new center      :  1.214  1.581  1.110 (nm)
new box vectors :  2.429  3.163  2.219 (nm)
new box angles  :  90.00  90.00  90.00 (degrees)
new box volume  :  17.05               (nm^3)

とある。
つまり、143原子を読み込んで、ボリュームが3.76725平方立方ナノメートルで、
おおよそ1600個の電子がある。と。
速度(ベクトルのことだろう?)は見つかってない。(まだ速度は指定されていない、
ということだろう。MDを実行する前は分子の速度は0になっている。分子の状態は
X,Y,Z座標のそれぞれの「位置」とそれぞれの向きに対する「速度」、あわせて6個の
状態変数によって決定される。本番MDを実行する前に適当な初速度を与えるために
通常Position Restrained MDという短めのMDをしなくてはならない。)

システムサイズが　	1.429  2.163  1.219 (nm)
中心が			2.469  1.077  2.517 (nm)
角度が			90.00  90.00  90.00 (degrees)
(つまり立方体ということかな。)
ボリューム(体積)が	3.77               (nm^3)

でもってこのシステムを	-1.254  0.504 -1.407 (nm)
だけシフト(移動)させた。と。

すると、あたらしいシステムは
中心が			1.214  1.581  1.110 (nm)
システムサイズが	2.429  3.163  2.219 (nm)
角度が			90.00  90.00  90.00 (degrees)
体積が			17.05               (nm^3)

と変わったと。
なるほでょ。

ためしに、-d 0.0 としてみてやったらどうなるのかな。
じっさいにスクリプトを書き換えてやってみた。
すると出力は次のように変わった。

#Entries in atommass.dat: 82 vdwradii.dat: 26 dgsolv.dat: 7
Read 143 atoms
Volume: 3.76725 nm^3, corresponds to roughly 1600 electrons
No velocities found
    system size :  1.429  2.163  1.219 (nm)
    center      :  2.469  1.077  2.517 (nm)
    box vectors :  1.429  2.164  1.218 (nm)
    box angles  :  90.00  90.00  90.00 (degrees)
    box volume  :   3.77               (nm^3)
(ここまでは同じ)
    shift       : -1.754  0.004 -1.907 (nm)
new center      :  0.714  1.082  0.610 (nm)
new box vectors :  1.429  2.163  1.219 (nm)
new box angles  :  90.00  90.00  90.00 (degrees)
new box volume  :   3.77               (nm^3)


(前回 -d 0.5 の場合)
    shift       : -1.254  0.504 -1.407 (nm)
new center      :  1.214  1.581  1.110 (nm)
new box vectors :  2.429  3.163  2.219 (nm)
new box angles  :  90.00  90.00  90.00 (degrees)
new box volume  :  17.05               (nm^3)

ん!?
ボリュームサイズが変わってない?
なるほど。
0.5 というのは、片辺0.5nm拡大ということで上下左右幅1.0nm広がるということね。
I got it!!

念のため、-d 1.0 の場合もやってみる。

    system size :  1.429  2.163  1.219 (nm)
    box vectors :  1.429  2.164  1.218 (nm)
    box volume  :   3.77               (nm^3)
new box vectors :  3.429  4.163  3.219 (nm)
new box volume  :  45.95               (nm^3)

やっぱり。そのようだ。

次。

genbox -cp cpeputide.gro -cs -o cpeputide_b4em.gro -p cpeputide.top >>& ! output.genbox

オプションは

-cp
-cs
-o
-p

の4つ。恒例のごとくマニュアルを見る。
(http://www.gromacs.org/documentation/reference_3.2/online/genbox.html)

こちらも長いが焦らず読んでみよう。

Genbox can do one of 3 things:

(要約)
genboxコマンドはどうも3つのことができるらしい。
以下の(1)(2)(3)。

1)
Generate a box of solvent. Specify -cs and -box. Or specify -cs and -cp with
 a structure file with a box, but without atoms.

(要約)
まず一つ目。ボックスの中に溶媒(通常、水)を入れる。と。
-csオプションと-boxオプションを指定するか、-csオプションと-cpオプションで構造
ファイルを指定する。かのどちらか。

おそらくここで言っている structure file とはトポロジーファイル(.top)のことだろう。

次。2個目。
2) 
Solvate a solute configuration, eg. a protein, in a bath of solvent molecules. 
Specify -cp (solute) and -cs (solvent). The box specified in the solute
 coordinate file (-cp) is used, unless -box is set, which also centers
 the solute. 
The program editconf has more sophisticated options to change the box and
 center the solute. 
Solvent molecules are removed from the box where the distance between
 any atom of the solute molecule(s) and any atom of the solvent molecule is
 less than the sum of the VanderWaals radii of both atoms. 
A database (vdwradii.dat) of VanderWaals radii is read by the program,
 atoms not in the database are assigned a default distance -vdw.

(要約)
-cpオプションで溶媒、-csオプションで溶質を指定できる。と。
もし-boxオプションを使わなければ、-cpオプションで溶質の座標ファイルをしてしろ、と。
editconfはもっとボックスに関していろいろ設定できるよ。と。
それから、あまりに原子同士が近い場合は、取り除かれるよ。と。

3つ目。
3) 
Insert a number (-nmol) of extra molecules (-ci) at random positions. 
The program iterates until nmol molecules have been inserted in the box. 
To test whether an insertion is successful the same VanderWaals criterium is
 used as for removal of solvent molecules. 
When no appropriately sized holes (holes that can hold an extra molecule) are
 available the program tries for -nmol * -try times before giving up. 
Increase -try if you have several small holes to fill.

(要約)
-ciオプションでその他の分子をランダムポジションで好きな数だけ-nmolで指定して
追加できるよ。と。

まだ解説があるけど、このへんで。
つまりは、溶媒分子を箱にぶち込んで、ペプチドを溶かした状態にする。
ということらしい。

つづいてオプション解説。

-cp 	Input, Opt. 	 Generic structure: gro g96 pdb tpr tpb tpa xml
-cs 	Input, Opt., Lib. 	Generic structure: gro g96 pdb tpr tpb tpa xml
-o 	out.gro 	Output 	Generic structure: gro g96 pdb xml
-p 	topol.top 	In/Out, Opt. 	Topology file

-cp	入力として.groファイル
-cs	入力として溶媒(水)ファイル
-o	出力として.groファイル
-p	入力として.topファイルを受け取り、そのファイルを修正して、
	そのまま同じ名前で.topファイルを出力する。

ということらしい。マニュアルを見ると
(http://www.gromacs.org/documentation/reference_3.2/online/genbox.html)

The default solvent is Simple Point Charge water (SPC), with coordinates from
$GMXLIB/spc216.gro.

とある。つまり、-csオプションで特に指定がなければデフォルトでspc216.groの
水が使われる。と。

ここで、わたしのやるCoarse Grain Modelについて少し。
わたしのブツは、水分子を一つの原子に(仮にWとする)粗く抽象化したものを使う。
のでこのままのspc216.groの水分子だと都合が悪い。
でどうしたかというと、
以下のようなw.groファイルというのをつくった。

----
W, Water molecule in Coarse Grain Model
1
   1W       W 1    .230    .628    .113
   1.86206   1.86206   1.86206
----

でもって、

genbox -cp conf.gro -ci w.gro -nmol 1000 -p topol.top -o conf.gro

とかすると、W分子を1000個、conf.groのランダムな位置に挿入してくれる。
あまりに濃度を高くすると(箱のサイズが小さいと)分子間距離が極端に短くなり、
シミュレーションが動かないことがあった。
そういう場合は、

editconf -f conf.gro -o conf.gro -scale 2.0 2.0 2.0

のようにすれば、箱を広げてくれるそうだが、うまくばらけて配置してくれなかった
ので、また一からやり直した。

-ciオプションは、他の任意の分子を指定して、-nmolで挿入する分子数を指定できる。

ただこのようにやってもtopol.topファイルは更新してくれないようなので
(どういうわけか？)手動でtopol.topの最後の行に

W 1000

とか記述した。
10000分子のW(水分子)を挿入しようと試みたが、FatalErrorとなって8000個くらいで
いつもプログラムが止まってしまった。どうもメモリの容量が少ないととまってしまうらしい。
仕方がないので、水分子8000個のシミュレーションを実行中。
一気に10000個入れようとするダメなので、
いくつかに分割して、たとえば5000個、5000個で分子を挿入する。

genbox -cp conf.gro -ci w.gro -nmol 5000 -p topol.top -o conf.gro
genbox -cp conf.gro -ci w.gro -nmol 5000 -p topol.top -o conf.gro

とすると2回目の時にはまったく同じ場所に水分子が挿入されてしまうので、
乱数シードを変える

-seedオプション

をつけてあげる。デフォルトは1999らしいので、

genbox -cp conf.gro -ci w.gro -nmol 5000 -p topol.top -o conf.gro -seed 1234
genbox -cp conf.gro -ci w.gro -nmol 5000 -p topol.top -o conf.gro -seed 5678

とかなんとかして実行すると10000個の分子を挿入できる。これでもFatalErrorが出るときは
さらに細かく1000個ずつとかやってみる。


さて。
新しくできたファイルを見てみる(ls)と、

#cpeptide.gro.1#
#cpeptide.top.1#
cpeptide.top
cpeptide.gro
cpeptide_b4em.gro
output.genbox

どうも同じ名前のファイル名で出力するときは前あったファイルに##マークがついて
バックアップされるようだ。
なので、

あたらしく、

cpeptide.top
cpeptide.gro
cpeptide_b4em.gro
output.genbox

のファイルができている。
中身を確認してみよう。

どう修正されたのかみるときは、diffコマンドを使おう。

diff #cpeptide.gro.1# cpeptide.gro

3,146c3,146
<     1LYSH     N    1   2.497  -0.065   2.231
<     1LYSH    H1    2   2.449  -0.140   2.186
<     1LYSH    H2    3   2.583  -0.046   2.184
<     1LYSH    H3    4   2.515  -0.089   2.327
(中略)
<    13MET     O1  142   2.790   1.946   2.215
<    13MET     O2  143   3.005   2.017   2.303
<    1.42915   2.16366   1.21831
---
>     1LYSH     N    1   1.243   0.439   0.824
>     1LYSH    H1    2   1.195   0.364   0.779
>     1LYSH    H2    3   1.329   0.458   0.777
>     1LYSH    H3    4   1.261   0.415   0.920
(中略)
>    13MET     O1  142   1.536   2.450   0.808
>    13MET     O2  143   1.751   2.521   0.896
>    2.42900   3.16300   2.21900

.groファイルは全行変わってる。そりゃそうだ。座標変換したんだもの。

次。

diff #cpeptide.top.1# cpeptide.top

916c916
< Protein
---
> Protein in water
920a921
> SOL               488

こちらはいたってシンプル。新しくできたファイルには、

SOL	488

の記述が加わっているだけ。これは水(solvent溶媒)をあらわしているのだろう。

次。
cpeptide_b4em.gro

----
Good gRace! Old Maple Actually Chews Slate
 1607
    1LYSH     N    1   1.243   0.439   0.824
    1LYSH    H1    2   1.195   0.364   0.779
    1LYSH    H2    3   1.329   0.458   0.777
    1LYSH    H3    4   1.261   0.415   0.920
    1LYSH    CA    5   1.158   0.559   0.820
    1LYSH    CB    6   1.027   0.536   0.898
    1LYSH    CG    7   0.922   0.645   0.863
    1LYSH    CD    8   0.796   0.616   0.949
    1LYSH    CE    9   0.817   0.670   1.090
    1LYSH    NZ   10   0.898   0.580   1.175
(中略)
   13MET     O1  142   1.536   2.450   0.808
   13MET     O2  143   1.751   2.521   0.896
   14SOL     OW  144   0.225   0.275   0.996
   14SOL    HW1  145   0.260   0.258   1.088
   14SOL    HW2  146   0.137   0.230   0.984
   15SOL     OW  147   0.019   0.368   0.647
   15SOL    HW1  148  -0.063   0.411   0.686
   15SOL    HW2  149  -0.009   0.295   0.584
(中略)
  500SOL    HW1 1603   1.845   2.179   1.866
  500SOL    HW2 1604   1.968   2.284   1.839
  501SOL     OW 1605   1.901   2.939   2.162
  501SOL    HW1 1606   2.000   2.928   2.153
  501SOL    HW2 1607   1.861   2.853   2.194
   2.42900   3.16300   2.21900
----
恒例のごとく最初の一行はよく分からない文章があるが、無視。
2行目は総原子数。
143行目を境にsolvent(水分子)が付け加えられている。
最後の行は、ボックスのサイズ。

なるほど。

次。
output.genbox

これは、editconfとgenboxコマンド両方の出力内容が書き込まれている。
前半は、editconfの出力なので飛ばして、
genboxの方の出力は、、、
長いので分割しつつ、、、
----
Reading solute configuration
Good gRace! Old Maple Actually Chews Slate
Containing 143 atoms in 13 residues
Initialising van der waals distances...
Reading solvent configuration
"216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984"
solvent configuration contains 648 atoms in 216 residues

ふんふん、ファイルを読み込んだと。

Initialising van der waals distances...
Will generate new solvent configuration of 2x2x2 boxes
Generating configuration
Sorting configuration
Found 1 molecule type:
    SOL (   3 atoms):  1728 residues
Calculating Overlap...
box_margin = 0.315
Removed 2901 atoms that were outside the box
Succesfully made neighbourlist
nri = 11447, nrj = 457532
Checking Protein-Solvent overlap: tested 4083 pairs, removed 231 atoms.
Checking Solvent-Solvent overlap: tested 36887 pairs, removed 588 atoms.
Added 488 molecules
Generated solvent containing 1464 atoms in 488 residues
Writing generated configuration to cpeptide_b4em.gro
Good gRace! Old Maple Actually Chews Slate

ん？ん？residues？removed？added？
なにやら分子数を操作してそうで、、、で、
結果的に

Output configuration contains 1607 atoms in 501 residues
Volume                 :     17.0484 (nm^3)
Density                :     1007.06 (g/l)
Number of SOL molecules:    488   

1607個の原子があって、ボリュームサイズが17立方ナノメートルで、
密度が1007g/lだと。
水分子が488個、と。

なるほど。

Processing topology

Back Off! I just backed up cpeptide.top to ./#cpeptide.top.1#

これはcpeptide.topのファイルがすでに存在していたので、
バックアップとして#cpeptide.top.1#にファイル名を変更した。と。

gcq#83: "Love is Like Moby Dick, Get Chewed and Get Spat Out" (Urban Dance Squad)

いつもなにやら最後にgcq#と偉人先人の有名ワードらしきセンテンスが表示される。
GROMACS開発者のおちゃめな一面(?)らしい。。。

(後の出力省略)
----

ということで、ここまでまとめると。
pdb2gmxコマンドで、.pdbファイルを.groファイルに変換して、
editconfコマンドで、boxのサイズを変更して、
genboxコマンドで水分子を入れて、.groファイル、.topファイルを修正
した。と。

1. pdb2gmx(.pdbを.groに変換)
2. editconf(.groの修正(箱の大きさを変更))
3. genbox(.gro, .topの修正(箱の中に水分子を入れる))
4. .mdpの作成(手動)
5. grompp(.mdp、.gro、.topから.tprを作成)
6. mdrun(シミュレーション)

全工程のうち、1,2,3が終了。次4.

demoを進める。(Enterキー)


----
In principle we can start a Molecular Dynamics simulation now. However
it is not very wise to do so, because our system is full of close
contacts. These close contacts are mainly a result of the genbox
program. The added solvent might have some close contacts with the
peptide resulting in very high repulsive energies. If we would start a
Molecular Dynamics (MD) simulation without energy minimisation the
system would not be stable because of these high energies.

The standard procedure to remove these close contacts is
Energy Minimisation (EM). Energy minimisation slightly changes the
coordinates of our system to remove high energies from our system.  

Before we can start the Energy Minimisation we have to preprocess all
the input files with the GROMACS preprocessor named grompp. grompp
preprocesses the topology file (.top), the structure file (.gro) and a
parameter file (.mdp) resulting in a binary topology file (.tpr
extension). This binary topology file contains all information for a 
simulation (in this case an energy minimisation).
----
(要約)
原理的にはシミュレーションを始められるけど、
このままだと原子同士が近すぎるものがあって(genboxコマンドのせいらしい)、
MDがうまく走らない。まずEnergy Minimisasion(エネルギー最小化)というのをしろ。と。

で、EMをするまえにgromppコマンドをつかって、
トポロジーファイル(.top)、構造ファイル(.gro)、パラメーターファイル(.mdp)
バイナリートポロジーファイル(.tpr)を作る必要がある。と。
このバイナリートポロジーファイルはシミュレーションに必要なすべての
情報を含んでいる。

なるほど。ということで先を進める。

(Enterキー)

なにやらコマンドが実行されたようなので、、、
demoスクリプトファイルを覗いてみる。

----
echo generating energy minimisation parameter file...
cat > em.mdp << _EOF_
title               =  ${MOL}
cpp                 =  /lib/cpp
define              =  -DFLEX_SPC
constraints         =  none
integrator          =  steep
dt                  =  0.002    ; ps !
nsteps              =  100
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
;
;       Energy minimizing stuff
;
emtol               =  1000.0
emstep              =  0.01
_EOF_

----
と、まずある。

echo generating energy minimisation parameter file...

は、ただ単にターミナル画面にgeneratingホゲホゲという文字を表示するだけで、

cat > em.mdp << _EOF_

は以下の行"_EOF_"という文字があるまで、em.mdpファイルに書き出す、
というスクリプトコマンド(というより単にcatコマンドの実行)

その後に、
grompp -f em -c ${MOL}_b4em -p ${MOL} -o ${MOL}_em >& ! output.grompp_em

というコマンドが実行されている。ここでバイナリートポロジーファイル(.tpr)
が作成されているのだろう。

ここも同様にシェル変数${MOL}を具体的に置き換えてみる。

grompp -f em -c cpeptide_b4em -p cpeptide -o cpeptide_em >& ! output.grompp_em

オプションは
-f
-c
-p
-o
の4つ。

マニュアルによると
(http://www.gromacs.org/documentation/reference_3.2/online/grompp.html)

-f 	grompp.mdp 	Input, Opt. 	 grompp input file with MD parameters
-c 	conf.gro 	Input		 Generic structure: gro g96 pdb tpr tpb tpa xml
-p 	topol.top 	Input		 Topology file
-o 	topol.tpr 	Output		Generic run input: tpr tpb tpa xml

とある。つまり、
-f	入力として.mdp(パラメータファイル)
-c	入力として.gro(GROMOSファイル)
-p	入力として.top(トポロジーファイル)を指定して、
-o	出力として.tprファイルを指定する、

らしい。
なるほど。

でもって次は、em.mdpファイルの中身を見ていこう。
とそのまえにそもそもこの.mdpという拡張子のファイルはなんなのか。

Chap3

でも説明したが、Tutorialのページによると
(http://www.gromacs.org/documentation/reference_3.2/online/getting_started.html#mdp)
Molecular Dynamics parameter file (.mdp)
 The Molecular Dynamics Parameter (.mdp) file contains all information about
 the Molecular Dynamics simulation itself e.g. time-step, number of steps,
 temperature, pressure etc. The easiest way of handling such a file is by
 adapting a sample .mdp file. A sample mdp file can be found online.

(要約)
パラメータファイル(.mdp)は、タイムステップ、ステップ数、温度、圧力などの
シミュレーションに関する情報が含まれていて、サンプルファイルを参考に
適当に修正して使ってね。

この設定如何でシミュレーションが決まるといってもいいくらい重要なファイルなので、
ここはゆっくり一つ一つの項目を確認して行こう。

.mdpのマニュアルを見ながら、em.mdpの内容をみてゆく。
(http://www.gromacs.org/documentation/reference_3.2/online/mdp_opt.html)

スクリプトファイルの中で出力されているのは次の項目。
----
title               =  ${MOL}
cpp                 =  /lib/cpp
define              =  -DFLEX_SPC
constraints         =  none
integrator          =  steep
dt                  =  0.002    ; ps !
nsteps              =  100
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
;
;       Energy minimizing stuff
;
emtol               =  1000.0
emstep              =  0.01
-----

「;」はコメント。

一つ一つ見てゆく。

title               =  ${MOL}

this is redundant, so you can type anything you want
ここはなんでも自由に書いてよい。と。

cpp                 =  /lib/cpp

your preprocessor
プリプロセッサ指定。ここも既存の設定のまま。

define              =  -DFLEX_SPC

defines to pass to the preprocessor, default is no defines. 
You can use any defines to control options in your customized topology files.
Options that are already available by default are: 

 -DFLEX_SPC
 Will tell grompp to include FLEX_SPC in stead of SPC into your topology,
 this is necessary to make conjugate gradients or l-bfgs minimization work
 and will allow steepest descent to minimize further.

 -DPOSRE
 Will tell grompp to include posre.itp into your topology, used for position
 restraints.

よくわからないが、Minimization workのときは-DFLEX_SPCオプションで、
position restraintsのときは-DPOSREオプションをつけろ、と？
l-bfgsというのがよくわからない。
-DFLEX_SPCはSPCの代わりにFLEX_SPCを使う(topologyファイル.topに挿入する)そうだが、そもそもFLEX_SPCとSPCというのが何かよくわからない。水分子のことなのだろうが。。。(保留)


constraints         =  none

 Bonds
 constraints:

 none
  No constraints, i.e. bonds are represented by a harmonic or a Morse potential
  (depending on the setting of morse) and angles by a harmonic potential. 
 hbonds
  Only constrain the bonds with H-atoms.
 all-bonds
  Constrain all bonds.
 h-angles
  Constrain all bonds and constrain the angles that involve H-atoms by adding
  bond-constraints.
 all-angles
  Constrain all bonds and constrain all angles by adding bond-constraints.

Bondsコーナーに解説があるということから、
「結合」に関しての何か設定なのだろう。
constraints自体がどうも専門用語のよう。
alcで調べると
強制｛きょうせい｝、制約｛せいやく｝、制限｛せいげん｝すること、しがらみ、
抑制｛よくせい｝、遠慮｛えんりょ｝
などとある。
Google様にお告げを聞いてみる。
「分子動力学 constraint」

ちょっとわからない。
 none
  No constraints, i.e. bonds are represented by a harmonic or a Morse potential
  (depending on the setting of morse) and angles by a harmonic potential. 

noneを指定した場合。constraintsをしない。つまり、、、
結合がharmonicポテンシャルもしくはMorseポテンシャルとharmonicポテンシャル
によって決まる角度とであらわされる。と。
制限なし？
どういう角度でも結合しうる。ということか？
harmonicポテンシャル、MorseポテンシャルというのはForceField(力場)のところで
登場するポテンシャル関数のこと。
詳しくはしらない。

とりあえず、Energy Minimisation時にはnoneにするのがよさそう、というのが予想できる。(ほんとかな〜？)

integrator          =  steep

 Run control
 integrator:
  md
   A leap-frog algorithm for integrating Newton's equations of motion.
  steep
   A steepest descent algorithm for energy minimization. 
   The maximum step size is emstep [nm], the tolerance is emtol [kJ mol-1 nm-1].

Run Controlのセクションに説明があるということは、実行(計算)に関しての
設定だろう。integrateとは積分のことなので、おそらくここは数値計算(数値積分)
のアルゴリズムを設定する場所ではなかろうか。

だいたいsteepとmdの二つが良く使われるようだ。
steepはEnergyMinimization時に、普通のMD計算の時にはmdを指定すればよさげ。
mdを指定するとニュートンの運動方程式に対してリープフロッグというアルゴリズム
が使われると。(そのうち数値積分法についても解説したい。)
steepを指定するとsteepest descent algorithm(最急降下法)というのがつかわれるらしい。
リープフロッグにしても最急降下法にしてもシンプルで有名どころのアルゴリズム
である。

steepを指定したときは、最大のステップサイズをemstepで指定して、
tolerance(許容範囲?誤差?)をemtolで指定する、らしい。。。
ステップサイズの単位[nm]ナノメートルというのがよくわからない。
通常数値積分をする際のステップサイズといったら「時間」がくるのだが、
空間に対して積分している、ということか?
toleranceの単位をみていみると[kJ/(mol×nm)]
つまり[熱量/(分子数×距離)]
なので、1モルの分子・1ナノメートル当たり、の熱量(キロジュール)が
emtolで指定した値以下になるように分子をばらけさせる、ということか。

このemstepとemtolの設定が効いてきそう。

次。

dt                  =  0.002    ; ps !

 dt: (0.001) [ps]
  time step for integration (only makes sense for integrators md, sd and bd)

これはおなじみ、時間のステップサイズ。でもintegratorにmd,sd,bdを使ったときしか
意味がない、と書いてある。steepの場合は必要ないのか？
デフォルトで0.001ピコ秒らしい。
でも今回EnergyMinimisationでsteepが指定してあるにもかかわらず、dtのパラメータも
設定してある。。。。ん〜わからない。
。。

nsteps              =  100

 nsteps: (0)
  maximum number of steps to integrate

こちらもおなじみ。数値積分の計算ステップ数。
つまり今回の例だと、
dt=0.02[ps]
で
nsteps=100
なので
合計2.00[ps]分の計算が行われる。はず。


nstlist             =  10

 Neighbor searching
 nstlist: (10) [steps]
  Frequency to update the neighbor list
 (and the long-range forces, when using twin-range cut-off's). 
 When this is 0, the neighbor list is made only once.

Neighbor Searchingのセクションに説明がある。
直訳すると「お隣さん検索」
Frequencey(頻度) to update(更新) the neighbor list(お隣さんリスト)
お隣さんリストを更新する頻度。
ということらしい。単位は「steps」何ステップごとにお隣さんを更新するのか？
ということらしい。
これが0だとお隣さんリストは一度しか作られない？

おそらくお隣さんリストというのは、ある原子の周りにある原子の何かを
更新するということなのだろう。
ポテンシャルを計算するということなのかな？
ん〜うまくイメージできない。もうすこしMDの勉強をしてからもう一度考えてみよう。
(--ひとまずおいておくリスト---
* Topologyファイル(.top)の構造
* pdb2gmx実行時のForceFieldの指定
* お隣さんリスト更新
)

ns_type             =  grid

 ns_type:
  grid
   Make a grid in the box and only check atoms in neighboring grid cells
   when constructing a new neighbor list every nstlist steps. 
   In large systems grid search is much faster than simple search.
  simple
   Check every atom in the box when constructing a new neighbor list 
   every nstlist steps. 

こちらもNeighbor Searchingのセクションにある項目。
設定値としては、gridかsimpleのどちらか。
gridはボックスをグリッド(格子状)に区切って、隣の格子の原子を調べると。
nstlistステップごとに。
こんなイメージかな？

       |          |         |
       |          |         |
       |          |         |
  -----|----------|---------|--------
       | ○       |     ○  |
       |    ○    |  ○     |
       |          |  ○     |
  -----------------------------------
       | ○       |         |
       |   ○ ○  |  ○ ○  |
       |          |         |
  -----------------------------------
       |          |         |
       |          |         |
       |          |         |

rlist               =  1.0

 rlist: (1) [nm]
  cut-off distance for the short-range neighbor list

カットオフ距離
カットオフというのは簡単に言うと、どこらへんまで原子間相互作用を
考慮するかという距離でこの距離よりも遠い原子との相互作用は計算しない、
ということ。
ただしここで指定してあるのはshort-range(短い距離)のneighbor listに
対して。ということ。

つまりは、long-rangeとshort-rangeとで計算を分けている、ということか。
。。。のちのち解決しよう。

rcoulomb            =  1.0

 Electrostatics and VdW
 rcoulomb: (1) [nm]
  distance for the Coulomb cut-off

Electrostatics(静電気力)とVdW(van der Waal)ファンデルワールス力
の項目にある。
ここでちょっと原子間に働く「力」についてまとめておこう。
bond共有結合力
angle結合角
torsion angle捻れ角
VdWファンデルワールス力
水素結合(力)
Coulomb(Electrostatic)クーロン力(静電相互作用)
このうち、VdW、水素結合、クーロン力は比較的遠い原子間でも相互作用をする「力」
あまりに遠くまで計算しているとたいへんなのである程度のところで相互作用を
断ち切る設定がこのcut-off値。

rvdw                =  1.0

 rvdw: (1) [nm]
  distance for the LJ or Buckingham cut-off

こちらは、Lennard-JonesポテンシャルかBuckinghamポテンシャルのカットオフ。
VdWはこのLennard-Jonesポテンシャルで近似されることが多いらしい。

;
;       Energy minimizing stuff
;
emtol               =  1000.0
emstep              =  0.01

ここがさっきのintegratorでsteepを設定したときの追加設定項目。

 Energy minimization
  emtol: (100.0) [kJ mol-1 nm-1]
   the minimization is converged when the maximum force is smaller than this value 
  emstep: (0.01) [nm]
   initial step-size

ちゃんとEnergyMinimisationの項目に書いてある。
emtolはconverge(収束)する値。単位は先ほどのように[kJ/(mol×nm)]
この値よりも下回れば計算が終わるみたい。
emstepは初期ステップサイズ。これは計算途中である程度変化するみたい。

うまくEnergyMimisationが終わらなければ、この二つの値を変更させてみてみよう。

以上!!
----
長かった。。。。


だんだんとファイル数が多くなってきたが、じっくり前回のファイルリストと比べて
みてみると、、、新しくできたファイルは

cpeptide_em.tpr
em.mdp
mdout.mdp
output.grompp_em

の4つらしい。ここでも期待してないmdout.mdpというファイルが勝手に作成されて
いる。

em.mdpはついさっきみたばっかり、
同じ拡張子のmdout.mdpというのが出てきる。
これはem.mdpに足りない項目をすべて付け足してできたファイルのようだ。
コメントもちゃんとついている。

cpeptide_em.tprは。。。
これはバイナリファイルなので、テキストエディタではみれない。
このファイルの中身を見るコマンドがあり、

gmxdump

というコマンドで見れるらしい。

ためしに見てみる。

gmxdump -s cpeptide_em.tpr

...

長すぎるので、却下。
おそらく、.gro、.top、.mdpの情報が入っているのだろう。。。

output.grompp_em
を見てみると、、、、
恒例のごとくGROMACSの宣伝と使い方が表示された後、

creating statusfile for 1 node...
calling /lib/cpp...
processing topology...
#   G96BONDS:   432
#  G96ANGLES:   836
#      PDIHS:   370
#      IDIHS:   330
#       LJ14:   675
#     SETTLE:   976
Analysing residue names:
There are:   488      OTHER residues
There are:    13    PROTEIN residues
There are:     0        DNA residues
Analysing Protein...
Analysing Other...

とある。
まあ、これはいいだろう。


で次。(Enterキー)

-----------------------------------------------------------------
Now the binary topology file is generated, we can start the energy
minimisation (EM). The program which performs the EM is called
mdrun. In fact all simulations are performed by the same program:
mdrun. 

As the Energy Minimisation is running, watch the output of the
program. The first number ( from left to right ) is the number of the
iteration step. The second number is the step size, which gives an
indication of the change in the system. The third number is the
potential energy of the system. This number starts at a high value and
rapidly drops down, and converges, to a large negative value. 
-----------------------------------------------------------------
(要約)
準備ができたので、EM(エネルギー最小化)をスタートするよん。
コマンドはmdrun。
シミュレーションが実行されいているとき、
ウィンドウに
ステップ数、ステップサイズ、ポテンシャルエネルギー
が順番に表示されいているからみてみてねん。

と。

ウィンドウに注意しつつ、Enterを押す。
たしかにずらずらとなにやら数値が表示されている。

Step=   33, Dmax= 3.1e-03 nm, Epot= -1.70902e+04 Fmax= 1.72036e+03, atom= 1256
Step=   35, Dmax= 1.9e-03 nm, Epot= -1.71066e+04 Fmax= 1.70602e+03, atom= 1256
Step=   37, Dmax= 1.1e-03 nm, Epot= -1.71077e+04 Fmax= 1.70588e+03, atom= 1522
Step=   48, Dmax= 1.3e-06 nm, Epot= -1.71077e+04 Fmax= 1.96611e+03, atom= 130

スクリプトを覗いてみる。
シミュレーション時に実行されたのは、次のコマンド。

mdrun -nice 4 -s ${MOL}_em -o ${MOL}_em -c ${MOL}_b4pr -v >& ! output.mdrun_em

恒例のごとく、シェル変数を置き換えて具体的にしてみると、

mdrun -nice 4 -s cpeptide_em -o cpeptide_em -c cpeptide_b4pr -v >& ! output.mdrun_em

となる。

オプションは
-nice
-s
-o
-c
-v
の5つ。

また恒例のようにマニュアルを見る。
(http://www.gromacs.org/documentation/reference_3.2/online/mdrun.html)

The mdrun program reads the run input file (-s) and distributes the topology
 over nodes if needed. 
The coordinates are passed around, so that computations can begin. 
First a neighborlist is made, then the forces are computed. 
The forces are globally summed, and the velocities and positions are updated. 
If necessary shake is performed to constrain bond lengths and/or bond angles. 
Temperature and Pressure can be controlled using weak coupling to a bath.

オプション-sで入力ファイル(.tpr)を指定して、

mdrun produces at least three output file, plus one log file (-g) per node. 
The trajectory file (-o), contains coordinates, velocities and optionally
 forces. 
The structure file (-c) contains the coordinates and velocities of the 
last step. 
The energy file (-e) contains energies, the temperature, pressure, etc, 
a lot of these things are also printed in the log file of node 0. 
Optionally coordinates can be written to a compressed trajectory file (-x).

えーと、
mdrunは最低3つのファイルを出力する。
-gでログファイル(.log)
-oでトラジェクトリーファイル(.trr)。これには、分子の座標、速度が
書かれている。
-cで構造ファイル(.gro)。最後のステップでの分子の座標と速度が書かれている。
-eでエネルギーファイル(.egr)。これにはエネルギー、温度、圧量などが
書かれている。
-xで圧縮したトラジェクトリーファイルを出力することもできる。と。

なるほど。

これでCスクリプトに記述されていた、
-nice
-s
-o
-c
-v
のうち、

-s
-o
-c

が解決した。と。
で、残り
-nice
-v
は、

-[no]v 	 bool		 no 	 Be loud and noisy
-nice 	 int		 19 	 Set the nicelevel

とある。
まとめると、
-nice	実行のプライオリティー(優先順位度)
-s	入力として.tprファイルを指定
-o	出力として、トラジェクトリーファイル(.trr)を指定
-c	出力として、構造ファイル(.gro)を指定
-v	？？？

最後の-vがよくわからない。be loud and noisy?
直訳すると、「大声でやかましく？」なんじゃそりゃ？
とりあえず、無視。

それぞれ
conf.gro
mdout.mdp
topopl.top
というデフォルトのファイル名を使っているなら、
単にmdrunとタイプするだけOKのようだ。

あたらしくできたファイルを確認。
cpeptide_b4pr.gro
cpeptide_em.trr
ener.edr
md.log
output.mdrun_em

期待通りのファイル出力。ログファイルを指定しないと自動的に
md.logというファイル名のログファイルが出力される。

cpeptide_b4pr.gro
をrasmolで見てみたがはじめと違いがよく分からない。。。
のでパス。
トラジェクトリーファイル.trrも今はパス。
まだこのファイルを解析するわけではないので。
同様にener.edrもパス。

output.mdrun_emとmd.logを見比べると、どうも、
md.logの方が詳しく情報が書かれている様子。
なのでmd.logのほうをチラッと覗いて見る。

またはじめにGROMACS広告があって、
パラメータの情報があって、
どうも、次らへんからが重要そう。。。
各ステップごとのBond(結合)エネルギー、Potentialエネルギー、Kineticエネルギー、
などが出力されている。

           Step           Time         Lambda
              0        0.00000        0.00000

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    3.03610e+02    2.56875e+02    7.71508e+01    6.90061e+00    1.34934e+02
     Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.
    1.46123e+03    1.00916e+04   -1.92894e+04   -6.95719e+03    0.00000e+00
   Total Energy    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00    0.00000e+00

           Step           Time         Lambda
              1        1.00000        0.00000

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.15203e+02    2.01111e+02    7.68892e+01    6.19686e+00    1.31628e+02
     Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.
    1.45763e+03    7.33652e+03   -1.94914e+04   -1.01663e+04    0.00000e+00
   Total Energy    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00    0.00000e+00

           Step           Time         Lambda
              2        2.00000        0.00000

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    7.48740e+01    1.42772e+02    7.66344e+01    8.27204e+00    1.24395e+02
     Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.
    1.45361e+03    5.47839e+03   -1.98514e+04   -1.24925e+04    0.00000e+00
   Total Energy    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00    0.00000e+00

(中略)
           Step           Time         Lambda
             37       37.00000        0.00000

   Energies (kJ/mol)
        G96Bond       G96Angle    Proper Dih.  Improper Dih.          LJ-14
    1.71695e+01    6.82935e+01    7.35121e+01    2.83412e+01    7.16714e+01
     Coulomb-14        LJ (SR)   Coulomb (SR)      Potential    Kinetic En.
    1.45070e+03    2.48046e+03   -2.12978e+04   -1.71077e+04    0.00000e+00
   Total Energy    Temperature Pressure (bar)
    0.00000e+00    0.00000e+00    0.00000e+00

(中略)

Stepsize too small, or no change in energy.
Converged to machine precision,
but not to the requested precision Fmax < 1000

Double precision normally gives you higher accuracy.
You might need to increase your constraint accuracy, or turn
off constraints alltogether (set constraints = none in mdp file)

Steepest Descents converged to machine precision in 49 steps,
but did not reach the requested Fmax < 1000.
Potential Energy  = -1.7107693e+04
Maximum force     =  1.7058767e+03 on atom 1522
Norm of force     =  3.5140055e+04

	M E G A - F L O P S   A C C O U N T I N G

   RF=Reaction-field  Free=Free Energy  SC=Softcore
   T=Tabulated        S=Solvent         W=Water     WW=Water-Water

           Computing:      M-Number      M-Flop's  % Flop's
                   LJ      0.298945      9.267295     1.7
              Coulomb      0.693479     18.723933     3.4
           Coulomb(W)      0.037901      3.069981     0.6
         LJ + Coulomb      0.444072     16.874736     3.1
         LJ + Coul(W)      0.120106     11.049752     2.0
        LJ + Coul(WW)      1.466400    359.268000    65.0
      Innerloop-Iatom      0.979420      9.794200     1.8
             NS-Pairs      4.936512    103.666752    18.8
         Reset In Box      0.078743      0.708687     0.1
              Shift-X      0.078743      0.472458     0.1
               CG-CoM      0.026705      0.774445     0.1
                Bonds      0.007056      0.303408     0.1
               Angles      0.010241      1.669283     0.3
              Propers      0.003626      0.830354     0.2
            Impropers      0.003234      0.672672     0.1
               Virial      0.001323      0.023814     0.0
               Settle      0.047824     15.447152     2.8
          Total                   552.61692   100.0

               NODE (s)   Real (s)      (%)
       Time:      3.000      3.000    100.0
Finished mdrun on node 0 Sun Oct 24 06:57:18 2004

----
と最後に各エネルギーの情報が出力されて終わっている。
途中、「まとめ的」記述があり、

Steepest Descents converged to machine precision in 49 steps,
but did not reach the requested Fmax < 1000.
Potential Energy  = -1.7107693e+04
Maximum force     =  1.7058767e+03 on atom 1522
Norm of force     =  3.5140055e+04

49ステップでEMが終わったと。
1000ステップもする必要なかったよ。というわけですね。ふんふん。
Potentialエネルギーが-1.7107693e+04だったと。

ま、こんなかんじで次に行こう。

(Enterキー)

----
Once all close contacts are removed from the system, we want to do
molecular dynamics of the water molecules, and keep the peptide
fixed. This is called position restrained (PR) MD.

Position Restrained MD keeps the peptide fixed and lets all water
molecules equilibrate around the peptide in order to fill holes
etc. which were not filled by the genbox program.

We are first going to preprocess the input files to generate the
binary topology. The input files are the topology file, the structure
file ( output of the EM ) a paremeter file, and an index file.

By default our system is split into several groups. In this case we
use two of those groups: Protein and SOL(vent). We use these groups to
put position restraints on all the atoms of the peptide.

The parameter file ( .mdp extension ) contains all information about
the PR-MD like: step size, number of steps, temperature, etc. This
Paramter file also tells GROMACS what kind of simulation should be
performed ( like EM, PR-MD and MD etc. )
----

(要約)
次にPosition Restrained MDを行う。
これは、ペプチドを固定したまま、水分子を適当に動かして、
holes(穴)を埋める。

恒例のように.mdpファイルを作成し、gromppコマンドで.tprファイルを作って、
mdrunする。

なるへそ。

次。(Enterキー)

今回新しく

cpeptide_pr.tpr
output.grompp_pr
pr.mdp

3つファイルができた。
cpeptide_pr.tprとoutput.grompp_prはもう見なくてもいいだろう。
ということでpr.mdpをみていこう！

実行コマンドは、

cat > pr.mdp << _EOF_
title               =  ${MOL} position restraining
warnings            =  10
cpp                 =  /lib/cpp
define              =  -DPOSRES
constraints         =  all-bonds
integrator          =  md
dt                  =  0.002    ; ps !
nsteps              =  500      ; total 1.0 ps.
nstcomm             =  1
nstxout             =  10
nstvout             =  1000
nstfout             =  0
nstlog              =  10
nstenergy           =  10
nstlist             =  10
ns_type             =  grid
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  berendsen
tau_t               =  0.1              0.1
tc-grps             =  protein          sol
ref_t               =  300              300
; Pressure coupling is not on
Pcoupl              =  no
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529
_EOF_

と

grompp -f pr -c cpeptide_b4pr -r cpeptide_b4pr -p cpeptide -o cpeptide_pr >& ! output.grompp_pr

の二つ。はじめの方がpr.mdpファイル(パラメータファイル)を出力するコマンドで、
二つ目がmdrun入力ファイル(.tpr)を出力するコマンドgrompp。

gromppコマンドは前回出てきたので省略。
ただ、pr.mdpファイルは前回のem.mdpファイルとちょっと違う。
項目が増えてる。よくよく見比べて、増えた部分のパラメータと同じパラメータを
使っていて値が変わっているところを見てみよう。

1行目、
----
em.mdp
title               =  cpeptide

title               =  cpeptide position restraining




----
* Analysis

シミュレーションが終わる(mdrun実行後)と
md.log
traj.trr
ener.edr
confout.gro
というファイルができる。

それぞれを簡単に説明すると、

md.logは
シミュレーションのログファイル。
実行が上手くいかないときなど見て原因を検討する。

traj.trrは
トラジェクトリーファイル。すべての分子の時間発展の軌道が
含まれている。trajectoryとは「軌道」「軌跡」という意味。
場合によっては結構大き目のファイル。

ener.edrは
エネルギーファイル。エネルギーに関しての情報が含まれているファイル。

confout.groは、
シミュレーション後の構造ファイル(structure file)。
このファイルを元にさらにシミュレーションを実行することも可能。

シミュレーションが終わって、まずはじめにやるのは分子をグラフィカルに
見ることではなかろうか。

シミュレーション後のconfout.groファイルを.pdb形式のファイルに変換して、
rasmol等のPDBファイル閲覧ソフトで見てみる。

.groファイルを.pdbファイルに変換するには、

editconf -f confout.gro -o confout.pdb

というようにする。

-f	が入力.groファイル
-o	が出力.pdbファイル


でも一つ。どのように分子が動いているのかというのをアニメーションで見る
簡易プログラムがGROMACSには付属している。

それが

ngmx

トラジェクトリファイル名がtraj.trrならば、単にngmxとタイプするだけでよい。
これを実行するときはシミュレーション設定ファイル(.tpr)も必要。

ngmxを起動すると、Groupというボックスが現れる。
ここで表示する分子を選択する。Systemという項目を選択すると
すべての分子が表示される。

メニューから[Display]-[Animate]を選択すると画面下部に再生巻戻しボタンが
表示される。そのまま再生ボタンを押すと分子の軌道がアニメーションが見れる
(画像はあまりきれいでないが)。コマ送りに表示することもできる。
[File]-[Quit]を選択するとプログラムが終了する。

さて。次に。

エネルギー状態をみてみてみましょう。

それをするのが、

g_energyコマンド。

デフォルトのener.edrというファイルがあるならば、単に

g_energy

とタイプするだけでOK。
そうでなければ、

g_energy -f [energy file].edr

と-fオプションで指定する。

g_energyコマンドを実行すると以下のような画面が出てくる。

End your selection with 0
   1=          Bond   2=      G96Angle   3=       LJ (SR)   4=  Coulomb (SR)
   5=     Potential   6=   Kinetic En.   7=  Total Energy   8=   Temperature
   9=Pressure (bar)  10=        Vir-XX  11=        Vir-XY  12=        Vir-XZ
  13=        Vir-YX  14=        Vir-YY  15=        Vir-YZ  16=        Vir-ZX
  17=        Vir-ZY  18=        Vir-ZZ  19= Pres-XX (bar)  20= Pres-XY (bar)
  21= Pres-XZ (bar)  22= Pres-YX (bar)  23= Pres-YY (bar)  24= Pres-YZ (bar)
  25= Pres-ZX (bar)  26= Pres-ZY (bar)  27= Pres-ZZ (bar)  28= #Surf*SurfTen
  29=          Mu-X  30=          Mu-Y  31=          Mu-Z  32=Coul-SR:EPE-EPE
  33=    LJ:EPE-EPE  34= Coul-SR:EPE-W  35=      LJ:EPE-W  36=   Coul-SR:W-W
  37=        LJ:W-W  38=         T-EPE  39=           T-W  40=      Lamb-EPE
  41=        Lamb-W

この中から、見たいエネルギー情報の番号を入力して、最後に0を入力すると
energy.xvg
というファイルが出力される。

たとえば、Potential EnergyとKinetic Energyを表示させたい場合は、

5 [Enter]
6 [Enter]
0 [Enter]

と入力する。ちなみにPotential Energyとは位置エネルギーのことでKinetic Energy
とは運動エネルギーのこと。
この.xvgという拡張子のファイルはxmgraceというグラフプロットツールのプロット
ファイル形式でそのまま

xmgrace energy.xvf

ただ。このままでは一番ひだりのデータ(つまり5PotentialEnergy)しかプロットされない。
どうしたら3本同時に表示できるのか、まだよくわからない。
いままでgnuplotを使っていたので、xmgraceはいまいちまだ使い方がよくわからない。日本ではxmgraceよりgnuplotのほうがよく使われているようなので、xmgraceの解説
ページはあまり見当たらない。そのうちxmgraceの使い方もしらべねば。
それかgnuplot用に変換するスクリプトを書いても良いが、いちいち変換するのもなぁ
。。。どうかと思う。

せっかくだから週末当たりにでもxmgraceの使い方でもしらべよう。

参考xmgrace
Grace本家
http://plasma-gate.weizmann.ac.il/Grace/
kensukeさんのページ
http://www.kensuke.org/misc/xmgrace.html


エネルギー選択のところで、

1=          Bond   2=      G96Angle   3=       LJ (SR)   4=  Coulomb (SR)
5=     Potential   6=   Kinetic En.   7=  Total Energy   8=   Temperature
9=Pressure (bar)  10=        Vir-XX  11=        Vir-XY  12=        Vir-XZ

1 Bond	   結合のエネルギー？
2 G96Angle 結合角？
3 LJ	   レナードジョーンズポテンシャル(分子間)？
4 Coulomb  クーロン力
5 Potential位置エネルギー
6 Kinetic  運動エネルギー
7 Total	   全エネルギー
8 Temp	   温度
9 press	   圧量

とここらへんまではなんとなくわかるが、それ以降は何を言ってるのか
よくわからない。とりあえず、保留。


あまりに出力ファイルサイズが大きい場合、出力するデータを間引きすることができる。
-skipオプション

g_energy -skip 10

とかすると、10個ステップ飛ばしでデータを出力してくれる。

xmgraceはいまいち使い方が良く分からない。そこで
xmgrace用のグラフファイルをgnuplot用に書き換えるrubyスクリプトを作った。
(どうも一部動かないオプションがあるようだが、、、いずれ修正予定)
どこかパスの通ったディレクトリ(私は~/binに入れてる)にxvgというファイル名で
実行権を与えて入れておくと便利。

xvg

とするとカレントディレクトリにある.xvgの拡張子がついたファイルを
gnuplotでグラフを表示してくれる。

xvg energy.xvg

とファイルを指定すればそのファイルだけを表示する。

xvg energy.xvg　-o

とすれば、.eps形式で出力

xvg energy.xvg -p 

で、.png形式

xvg energy.xvg　-dd

とすると、gnuplot用のデータを削除せずに残しておける。
複数の出力結果を一つのグラフにまとめたいときなどは、
-ddオプションをつけて、出てきた.pltファイルを修正する。
その場合は、gnuplotのコマンドを参考のこと。
複数のデータを一括して一つのグラフにするようなことは
できるようになっていない。
改良するするつもりもない。

----
#! /usr/bin/ruby

require "options"

=begin
memo
points to modify
1.if you open files with -size, if there are many files,
  you might not see the last graphs.
2.to merge several files to one graph.
=end

class Options
	def usage
		print <<EOF
usage:
 xvg [OPTION] [FILE.xvg]
 no option: show graphs of .xvg in current dir

options:
 -h	Help
 -dd	Don't Delete .dat .plt files
 -o	Output .eps file
 -p     output .png file
 -size  window size
EOF
	end
end
##
# main
#
desktop_size=1000

#
# check option
#
options=Options.new
options.check

#
# check help option 
#
options.help

#
# check window size option
#
window_size=-1
options.values.each_pair do |k,v|
	if k=="-size"
		window_size=v
	end
end

#
# search for .xvg files in current directory
#
filelist=Array.new
options.values.each_pair do |k,v|
	if v=="file"
		filelist << k
		p k
	end
end
if filelist.length==0
	Dir.foreach(".") do |file|
		if file=~/\.xvg/
			filelist << file
		end
	end
end

#
# read .xvg file
#

filelist.each_with_index do |file,file_no|
#file="DE.xvg"
datafile=Array.new
plotfile=Array.new
plotfiledata=""

title=""
linelabel=Array.new
datalineflag=0
open(file,"r").each do |line|
	# replacing for gnuplot script
	if line=~/^\#/
		line=""
	elsif line=~/^\@\s+title\s"(.+)"/
		title=$1
		line=""
	elsif line=~/^\@\ssubtitle\s"(.+)"/
		title+=" "+$1
		line=""
	elsif line=~/^\@\s+xaxis\s+label\s(".+")/
		line="set xlabel "+$1+"\n"
	elsif line=~/^\@\s+yaxis\s+label\s(".+")/
		tmp=$1.gsub(/\\S/,"^{")
		tmp=tmp.gsub(/\\N/,"}")
		line="set ylabel "+tmp+"\n"
	elsif line=~/^\@\ss([0-9])\slegend\s(".+")/
		linelabel << $2
		line=""
	elsif line=~/^\s/
		datafile << line
		line=""
		# initial process
		if datalineflag==0
			plotfiledata="set title \""+title+"\"\n"
			if options.values.keys.index("-o")
				plotfiledata+="set term post enhanced eps\n"
				plotfiledata+="set out \"#{file.gsub(/\.xvg/,".eps")}\"\n"
			end
			if options.values.keys.index("-p")
				plotfiledata+="set term png\n"
				plotfiledata+="set out \"#{file.gsub(/\.xvg/,".png")}\"\n"
			end
			plotfiledata+="plot "
			labelflag=0
			linelabel.each_with_index do |label,i|
                                if labelflag==0
                                       plotfiledata+="\"#{file.gsub(/\.xvg/,".dat")}\" u 1:#{i.to_i+2} t #{label} w l"
			               labelflag=1
                                else
                                        plotfiledata+=",\"#{file.gsub(/\.xvg/,".dat")}\" u 1:#{i.to_i+2} t #{label} w l"
                                end
			end
			if linelabel.size==0
				plotfiledata+="\"#{file.gsub(/\.xvg/,".dat")}\" u 1:2 not w l"
			end
			plotfiledata+="\n"
			datalineflag=1
			line=plotfiledata
		end
	else
		line=""
	end
	plotfile << line
end

=begin
#
#check data
#
plotfile.each do |line|
	print line
end
datafile.each do |line|
	print line
end
exit
=end

#
# output out.dat
#
open(file.gsub(/\.xvg/,".dat"),"w") do |out|
	datafile.each do |line|
		out.print line
	end
end

#
# output out.plt
#
open(file.gsub(/\.xvg/,".plt"),"w") do |out|
	plotfile.each do |line|
		out.print line
	end
end

#
# show graph
#
if window_size==-1
	system("gnuplot -persist #{file.gsub(/\.xvg/,".plt")}")
else
	x=window_size.to_i*file_no%desktop_size
	y=((file_no)*window_size.to_i/desktop_size).to_i*window_size.to_i+(((file_no)*window_size.to_i/desktop_size).to_i*20)
	system("gnuplot -geometry #{window_size}x#{window_size}+#{x}+#{y} -persist #{file.gsub(/\.xvg/,".plt")}")
end
end

#
# delete .dat .plt files
#
ddflag=0
if !options.values.keys.index("-dd")
	system("rm *.dat")
	system("rm *.plt")
end
----

----
mdpファイル(パラメータ)についてはとりあえず保留。

ここから、いきなりジャンプ。
分析(Analysis)用のコマンドに注目。

とりあえず、Getting Startに出ているAnalysis Commandをすべてリストアップ。

make_ndx?

[Water]
g_rdf
g_hbond
g_msd
[Methanol-Water]
g_density
[Ribonuclease S-peptide]
g_energy
g_rms
g_gyrate
g_rama
g_saltbr
my_dssp
[Protein unfolding]
g_sas

まずは、Water分子のシミュレーションの分析から。
http://www.gromacs.org/documentation/reference_3.2/online/water.html
チュートリアルにしたがって、mdrunを終わらせる。




まずは、ここから。
http://www.gromacs.org/documentation/reference_3.2/online/water.html

Getting Startの中の「Water」分子のAnalysis。

ndx file(index.ndx)

動径分布関数radial distribution function
http://homepage.mac.com/mike1336/md/app/j101_150/127_radialDistributionP/radialDistribution.html
動径分布関数とは、１つの分子から別の分子までの距離の分布を表した
関数です。X軸：距離、Y軸：個数


----
Radial Distribution Function
動径分布関数

動径分布関数とは、簡単に言うと、各分子間距離の分布のこと。
ある種類の分子と別のある種類の分子間がどの程度はなれているのか、
くっついているのかが分かる。

動径分布関数を求めるときは、

g_rdf

を用いる。インプットとして、インデックスファイル(index.ndx)か
MD用の入力バイナリファイル(topol.tpr)が必要。
たいていは必要に応じて、index.ndxを作る必要がある。

index.ndxは

make_ndx

で作れる。分子、原子をいったんグループ化してそのグループと他のグループの
分子原子間距離を計算する。

コマンド例
----
make_ndx -f topol.tpr -o index.ndx

1 & a OT1
1 & a OT2
1 & a OT3
2|3|4
1 & a OT4
1 & a OT5
1 & a OT6
6|7|8

q

----
topol.tprというデフォルトのファイル名でも、
-f オプションは省略できない。
デフォルトの出力ファイル名はindex.ndxなので、-oオプションは省略化。

なので、

make_ndx -f topol.tpr

とすればよい。
make_ndxを起動すると、サブコマンド入力状態となる。
どの分子原子をグループ化するかを決める。グループにはそれぞれ
番号がつけられる。

上の例では1がシステム全体
2がOT1全部
3がOT2全部
4がOT3全部
で
5がOT1,OT2,OT3を含む分子全体
となる。

index.ndxは、いわば、分子のグループを指定するようなもので、
ここで指定したものの分子の番号がグループとしてまとめられて、
index.ndxファイルに記述される。

わたしのシステムでの実行例
----
There are:  8050      OTHER residues
There are:     0    PROTEIN residues
There are:     0        DNA residues
Analysing Other...

  0 System              : 10000 atoms
  1 EPE                 :  2000 atoms
  2 W                   :  8000 atoms

 nr : group       !   'name' nr name   'splitch' nr    'l': list residues
 'a': atom        &   'del' nr         'splitres' nr   'h': help
 't': atom type   |   'keep' nr        'splitat' nr    
 'r': residue         'res': group res 'chain' char    'q': save and quit
 "name": group        'case': case sensitive           'v': verbose output   

> a P

Found 1000 atoms with name P

  3 P                   :  1000 atoms
> a E

Found 1000 atoms with name E

  4 E                   :  1000 atoms
> q
----

すでに、
0 System
1 EPE
2 W
のグループができている。
PPO(トポロジーファイル内では、P)とPEO(トポロジーファイル内では、E)と
水分子(W)とのそれぞれの分子間距離が知りたいので、
あとグループPとグループEをつくる。
ちなみにEPEとは、PPOとPEOの分子をまとめたグループ。

これでindex.ndxに5つのセクション、

[ System ]
[ EPE ]
[ W ]
[ P ]
[ E ]

ができ、それぞれの分子に対応する分子番号が並んでいる。


index.ndxが作成できたらg_rdfを実行する。

コマンド例
g_rdf -n index.ndx

ただこれだと、すべての計算時間においての分布が計算されるので、
平衡後の分布を調べたいときは、開始時間を指定する。

g_rdf -n index.ndx -b 10000

とかすると、1000psからの計算になる。

How many groups do you want to calculate the RDF of?

と聞かれるので、
計算したいグループ数を指定。
たとえば、ここで1と選択すると続いて、
計算したい分子原子グループを2つ選択することを要求される。
一つ目に選択したグループが基準のグループとなる。

How many groups do you want to calculate the RDF of?

とのきに、2と入力すると、3つのグループを指定することになる。
たとえば、
W
P
E
と3つ選択すると、
W-P間の分子間距離分布
W-E間の分子間距離分布
が計算される。


だいたいを知るだけで良く、計算時間がえらくかかるときは、

g_rdf -n index.ndx -b 10000 -dt 100

とかすると、100psおきに計算するので速く計算してくれる。
その代わり、計算量が少なくなるので、分布は若干正確じゃなくなるかもしれない。

計算結果はデフォルトで、

rdf.xvg

に出力される。xmgraceで表示するには、
          	
xmgrace rdf.xvg

と入力。xmgraceを使いたくないので、g_energyのときにつくったxvgを
使うときは、

xvg -f rdf.xvg

とする。



----
あまりに分子同士がくっついていると局所的なポテンシャルエネルギーが高すぎて、
EMをしても走らないことがある。
そういう時はボックスを大きくしてDensityを下げる。Densityをあげる場合にも同様の
テクニックが使える。

editconf -f conf.gro -o conf.gro -scale 0.5 0.5 0.5

とすると、箱のサイズをx,y,zそれぞれ50%に縮めて分子の位置を相対的に
ずらしてくれる。

Densityを調べたいときは、

Pcoupl=Berendsen

とかして、圧力を指定して数ステップは知らせた後、実行を止めて、

g_energy

コマンドで14:Density
を選択して、energy.xvgを見てみる。

面倒だが今のところそれしかやり方を知らない。


----
Appendix. GIFアニメーションの作り方

trjconv -f traj.trr -o out.pdb -s topol.tpr -dump 60000
これで、60000ps時のPDBファイルができる。

traj.trrとtopol.tprはデフォルトのファイル名なので、省略できる。
デフォルトでは、PDBファイルが出力ファイルでないので、
-oオプションでPDBファイルを指定する必要がある。
(拡張子を.pdbにしたファイル名にする)

trjconv -o out.pdb -sep -dt 1000
これで1000psごとにout_XX.pdbファイルが出来上がる。
-b、-eを指定しなければ、全時間区間を対象にフレームを作る。


以下、大まかな手順としては、
1. rasmolで各pdbファイルを読み込んで、
2. 適当な色変換、角度変換等の処理＋gifファイル保存のコマンドを
   書いたスクリプトを作る。
3. スクリプトをrasmolに読み込ませ各PDBファイルに対応するGIF画像を作成。
4. convertコマンドでアニメーションGIFを作成。

これらを実行するRubyスクリプトが以下。


----
#! /usr/bin/ruby

# 20050405 update Options.rb

require "Options.rb"

# default values
$trrfile="traj.trr"
$tprfile="topol.tpr"
$rasfile="input.ras"
$outputfile="anime.gif"
$dt=1
$molecular=0
$delay=10
$loop=1

options=Options.new
options.addhelptop(" makeanime [options]\n\n")
options.add("h",0,"  help")
options.add("f",1,"  [file]: input .trr file(default: traj.trr)")
options.add("s",1,"  [file]: input .tpr file(default: topol.tpr)")
options.add("r",1,"  [file]: input .ras file(default: input.ras. forget?, exec with -mr)")
options.add("o",1,"  [file]: output .gif file(default: anime.gif)")
options.add("dt",1," frame time[ps](default: 1 ps)")
options.add("dd",0," don't delete temporary files")
options.add("dl",1," delay time(default 10)")
options.add("lp",1," loop times(default 1)")
options.add("m",1,"  the molecular number(default: 0(system). forget?, exec with -g)")
options.add("g",0,"  just show the molecular group")
options.add("mr",0," make a sample(input.ras) file")
options.check

filelist = Array.new


#
# CAUTION
#
# if there is .gif file in current dir, it doesn't work.
#
#################################################

# options check
if options.value["h"].using==1
	options.help
end
if options.value["f"].using==1
	$trrfile=options.value["f"].args[0]
end

if options.value["s"].using==1
	$tprfile=options.value["s"].args[0]
end

if options.value["o"].using==1
	$outputfile=options.value["o"].args[0]
end

if options.value["dt"].using==1
	$dt=options.value["dt"].args[0]
end

if options.value["m"].using==1
	$molecular=options.value["m"].args[0]
end

if options.value["r"].using==1
	$rasfile=options.value["r"].args[0]
end

if options.value["dl"].using==1
	$delay=options.value["dl"].args[0]
end

if options.value["lp"].using==1
	$loop=options.value["lp"].args[0]
end

# make a ras sample
if options.value["mr"].using==1
	open("input.ras","w") do |ras|
			ras.print "spacefill on\n"
			#ras.print "rotate x 50\n"
			#ras.print "rotate y 50\n"
			#ras.print "rotate z 50\n"
			ras.print "background white\n"
			ras.print "select W\n"
			ras.print "color blue\n"
	end
	exit
end

# show the molecular group
if options.value["g"].using==1
	system("trjconv -f #{$trrfile} -s #{$tprfile} -o out.pdb -sep -dt #{$dt} < /dev/null >& out.dat")
	open("out.dat","r").each do |line|
		if line =~ /^Group/
			print line
		end
	end
	File.unlink("out.dat")
exit
end

# make .pdb files from traj.trr
open("in.dat","w") do |file|
	file.print "#{$molecular}\n"
end

system("trjconv -f #{$trrfile} -s #{$tprfile} -o out.pdb -sep -dt #{$dt} < in.dat")


#look for max value of number
maxval=0
Dir.foreach(".") do |file|
	if file=~/out_(\d+).pdb/
		if maxval<$1.to_i
			maxval=$1.to_i
		end		
	end
end

# rename filenames
maxdig=Math.log10(maxval).to_i

Dir.foreach(".") do |file|
	filename=""
	if file =~ /out_(\d+).pdb/
		if $1.to_i!=0
			num0=maxdig-Math.log10($1.to_i).to_i
		else
			num0=maxdig
		end
		filename="out_"
		num0.times do 
			filename+="0"
		end
		filename+="#{$1.to_i}.pdb"	
		print file," ",filename,"\n"
		system("mv #{file} #{filename}")
	end
	if filename!=""
		filelist << filename
	end
end

# make rasmol script files 
$inputras=Array.new
open("#{$rasfile}","r").each do |line|
	$inputras << line
end

filelist.sort.each do |file|
	#p file
	if file=~/(out_\d+).pdb/
		open("#{$1}.ras","w") do |ras|
			$inputras.each do |line|
				ras.print line
			end
			ras.print "write gif #{$1}.gif\n"
			ras.print "exit\n"
		end
	end
end


# make script for converting pdb to gif
open("convert.sh","w") do |out|
	out.print "#! /bin/sh\n"
	filelist.sort.each do |file|
		if file =~ /(out_\d+).pdb/
			out.print "rasmol -nodisplay -pdb #{$1}.pdb < #{$1}.ras\n"
		end
	end
end

# execute script
system("sh ./convert.sh")

# make script for convert command
open("makeanime.sh","w") do |out|
	out.print "#! /bin/sh\n"
	out.print "convert -loop #{$loop}"
	filelist.sort.each do |file|
		if file =~ /(out_\d+).pdb/
			out.print " -delay #{$delay} #{$1}.gif" 
			#out.print " #{$1}.gif" 
		end
	end
	out.print " anime.gif"
end

# execute convert command
system("sh ./makeanime.sh")


# delete temporary files
if options.value["dd"].using!=1
	File.unlink("makeanime.sh")
	print "delete makeanime.sh\n"
	File.unlink("convert.sh")
	print "delete convert.sh\n"
	system("rm out_*.*")
	print "delete out_*.*\n"
end


----

補足
これを実行させるのは基本的にUNIX上で。
同じディレクトリはRubyのライブラリ検索パスに
以下のOptions.rbを入れておく必要がある。
その他、converコマンドとrasmolコマンドを使っているので、
ImageMagicとRasMolがインストールされていることが前提。

デフォルトの出力ファイルは、
out_0.pdb
out_1.pdb
out_2.pdb
out_3.pdb
...
out_10.pdb
out_11.pdb
...
となるが、これだとGIFアニメ作成時にファイル名の順番が上手く
sortできないので、出力ファイル数にあわせて、例えば、

out_000.pdb
out_001.pdb
out_002.pdb
out_003.pdb
...
out_010.pdb
out_011.pdb
...

のようなファイル名に変換している。

その他メモ
Constant Pressureで実行した場合は画像の大きさ(体積)が変わるので注意
あまりタイムステップを大きくするとアニメにならない。
時間を表示できたらいいな。
convertできるのはおよそ130画像まで？
130以上の画像を使おうとすると、途中でプログラムがとまる。
convert、GIFを使う以外の他のアニメーションの作り方を検討すべきか?!

Options.rb
----
#! /usr/bin/ruby

# 20050408 update print help in order

class Option
	def initialize(name, num, help)
		@name=name
		@arg_num=num
		@args=Array.new
		@using=0	# 0:not use, 1:use
		@help=help
		@helptop=""
		@helplast=""
	end
	attr_accessor:name
	attr_accessor:arg_num
	attr_reader:args
	attr_accessor:using
	attr_reader:help
end

class Options
	def initialize
		@option=Hash.new
		@option_order=Array.new
		@filelist=Array.new
	end
	attr_reader:filelist

	def value
		@option
	end
	def add(op_name, op_arg_num, op_help)
		@option[op_name]=Option.new(op_name, op_arg_num, op_help)
		@option_order << op_name
	end
	def usage
		#print "** you should override this method\n"
		print "usage:\n"
		print @helptop
		print "options:\n"
		@option_order.each do |op_name|
			print " -",op_name,@option[op_name].help,"\n"
		end
		print @helplast
	end
	def addhelptop(str)
		@helptop=str
	end
	def addhelplast(str)
		@helplast=str
	end
	def help
		usage
		exit
	end
	def check
		op_name=""
		ARGV.each do |args|
			if args=~/-([a-z]+)/	# in the case of 'option'
				op_name=$1
				@option[op_name].using=1
			elsif args==""		# in the case of 'no options'
				print "** you should add the '-*' option before args.\n"
				help
			else			# in the case of 'value'
				if op_name==""		# simple check
					print "** you should add the option name before using check method.\n"
					help
				end
				@option[op_name].args << args
			end

		end
		# check arg_num, arg_name
		if @option.length>0
			@option.each do |k,v|
				if v.using==1 && v.args.length!=@option[k].arg_num
					print "** it's wrong number about args. about -#{v.name} option.\n"
					help
				end
			end
		end


	end
	def debug
		@option.each do |k,v|
			print "option:",k,"\tusing:",v.using,"\t\tvalue:"
			v.args.each do |arg|
				print arg," "
			end
			print "\n"
		end
	end

	# option function
	def filecheck
		filelist.each do |list|
			if !File.exist?(list)
				print "#{list} is not found.\n"
				exit
			end
		end
	end

end

if __FILE__==$0
# we check option name, must be "-XXX", and number of arguments.
# we don't check the type of values, integer, double, or string, something.

$input="in.dat"
$output="out.dat"
options=Options.new
options.addhelptop(" ruby Options.rb [options]\n\n")
options.add("h",0,": help")
options.add("f",1," [file]: input .dat file")
options.add("o",1," [file]: output .dat file")
options.check
options.debug
if options.value["f"].using==1
	if options.value["f"].args[0] != /\.dat/
		options.help
	else
		$input=options.value["f"].args[0]
	end
end
if options.value["o"].using==1
	if options.value["o"].args[0] !~ /\.dat/
		options.help
	else
		$output=options.value["o"].args[0]
	end
end
if options.value["h"].using==1
	options.help
end
options.filelist << $input
options.filecheck

open($output,"w") do |out|
open($input,"r").each do |line|
	out.print line
end
end

end

----

----
Appendix
.xvg(xmgrace用)ファイルをgnuplot用ファイルに変換するスクリプト



参考
convertコマンドでGIFアニメーション
http://www.naruto-u.ac.jp/~naosone/imageconvert.html

convert *.gif anime.gif
これでもできるみたい。でもこれだと-delayオプションや-loopオプションが入れられない？


----
Appendix g_rotacf

rotational corelation function
(回転相関関数?)

これはある分子同士を結んだベクトルが時間と共にどう変化していくか、という
ものを時系列に出力するコマンド。出力ファイルは.xvg形式(xmgrace用のファイル)。

相関関数(corelation function)とは、


-d doublet


----
Appendix Put in the wall

通常のPeriodic Boundary Condition(周期的境界条件)を適用した
場合のモデルは、Bulkと呼ばれる。Bulkとは大きなとか大量のとか
という意味があり、左右上下に無限に広がっている空間での
シミュレーションになるからである。
これに対して、ある座標軸一面に不活性の分子を配置した場合、
ごくごく薄い膜状の構造での分子的挙動を観察することができる。
膜表面上での観察をしたいときには、この「壁」に相当する分子を
モデルに組み入れる必要がある。

その解説。
私がいま扱っているのは、PEO-PPO-PEOというtriblock copolymer.

修正する点は、3箇所。
1. topol.topにWallに相当する分子を追加する
2. grompp.mdpにその分子が動かないように設定する
3. conf.groファイルに分子を追加する。

まず、topol.topファイルから。
[ atomtypes ]
の欄に新しい分子を追加する。

--
[ atomtypes ]
;name        mass      charge   ptype         c6           c12
WALL         72.0       0.000   A            0.0           0.0
--

massは分子の重さ。そのほかは良くわからない。
[ moleculetype ] 
のところにも分子の情報を追加。

--
[ moleculetype ] 
; molname       nrexcl    
WALL                0     
 
[ atoms ]  
;id type resnr residu atom cgnr   charge
 1   WALL    1     WALL       P   1      0     72.0
--

[ molecules ]
のところに分子数を記述する

--
[ molecules ] 
; name  number
EPE 50
WALL 3034
--

topol.topについては以上。詳しい数値に関しては、保留。。。

次に、grompp.mdpファイル。

重要なのは、
energygrps               = EPE WALL
; Stuff for the Wall
; No energy for wall-wall interaction
energygrp_excl = WALL WALL
; freeze the wall
freezegrps = WALL
freezedim = Y Y Y

の部分。
energgrpsはg_energyをしたときの出力を得られる様に。
energygrp_excelはよくわからない。
freezegrpsは動かないように設定する。
freezedimはどちらの方向に動かなくするかで、x,y,z全部の座標
を固定するので　全部　Yes Yes Yes　としておく。


今回はPressure Couplingで行った。
tcoupl                   = Berendsen
tc-grps                  = EPE WALL
tau_t                    = 1.0 1.0
ref_t                    = 350 350
Pcoupl                   = Berendsen
Pcoupltype               = Anisotropic
tau_p                    = 2.0 2.0 2.0
compressibility          = 5e-6  5e-6 5e-6 0 0 0
ref_p                    = 1.0 1.0 1.0 1.0 1.0 1.0


で、次にconf.gro

実際はこんな感じ。
--
PEO-2PPO-PEO Copolymers
 5034
    1EPE      E    1  18.647   8.588   9.471 -0.2220  0.1924 -0.0029
    1EPE      E    2  18.548   9.010   9.619 -0.1802  0.0030  0.1563
    1EPE      E    3  18.389   9.201   9.943  0.1363 -0.2364 -0.2098
    1EPE      E    4  18.071   9.250  10.262  0.1240 -0.0410  0.0621
    1EPE      E    5  17.809   9.620  10.405 -0.2436  0.4307 -0.2322
(中略)
 3081WALL     P 5031   4.171  18.500  18.000
 3082WALL     P 5032   4.171  19.000  18.000
 3083WALL     P 5033   4.171  19.500  18.000
 3084WALL     P 5034   4.171  20.000  18.000
  33.83100  20.21700  18.29300
--
手動で挿入するのは困難なので、分子を好きなところに配置できる
rubyスクリプトを作成。
使用するまえに、
元となるconf.gro.bak, topol.top.bakファイルを用意しておく必要がある。
使い方は、スクリプト先頭の各パラメータを参照のこと。
だいたい、
$mol_no=51		# starting molecular number
$atom_no=2001		# starting atom number
$between_mol=0.5	# wall molecular distance each other
$static_coordinate=0	# 0:x 1:y 2:z
$from_box_edge_or_molecule_edge=1	# 0:from box_edge 1:from_molecule_edge
$wall_position=1.0	# from edge
ここらへんを変えれば、つかえるだろう。

これをつかって、分子を挿入する前に、
editconfコマンドで適当にboxの大きさを調節しておいたほうが良いかも。
分子の座標はなるべく「−」(マイナス)でない方が好ましい。
(スクリプト中で分子位置の最大最小を使うため) 
挿入する分子の種類などを変えるときは、コードを書き換える。
挿入後は、

editconf -f conf.gro -o conf.pdb

などをしてrasmolで画像を確認すべし！

--
#! /home/masa/bin/ruby

$boxvec=[0,0,0]
$maxpos=[0.0,0.0,0.0]	# edge position of molecules
$minpos=[0.0,0.0,0.0]
$mol_no=51		# starting molecular number
$atom_no=2001		# starting atom number
$total_mol_no=0		# number of wall molecules


$pos=[0.0,0.0,0.0]

$between_mol=0.5	# wall molecular distance each other
$static_coordinate=0	# 0:x 1:y 2:z
$from_box_edge_or_molecule_edge=1	# 0:from box_edge 1:from_molecule_edge
$wall_position=1.0	# from edge

def putwall(file)
	# which coordinate is static?
	p0=$static_coordinate
	p1=($static_coordinate+1)%3
	p2=($static_coordinate+2)%3

	# set inisial position
	$pos[p1]=0.0
	$pos[p2]=0.0
	if $from_box_edge_or_molecule_edge==0
		$pos[p0]=$boxvec[p0]-$wall_position
	else
		$pos[p0]=$maxpos[p0]+$wall_position
	end
	# put wall	
	while $pos[p2]<$boxvec[p2]
	while $pos[p1]<$boxvec[p1]
		file << format("%5dWALL     P%5d%8.3f%8.3f%8.3f\n",$mol_no,$atom_no,$pos[0],$pos[1],$pos[2])
		$mol_no+=1
		$atom_no+=1
		$pos[p1]+=$between_mol	
		$total_mol_no+=1
	end
	$pos[p1]=0.0
	$pos[p2]+=$between_mol
	end

	# other side
	$pos[p1]=0.0
	$pos[p2]=0.0
	if $from_box_edge_or_molecule_edge==0
		$pos[p0]=$wall_position
	else
		$pos[p0]=$minpos[p0]-$wall_position
	end

	# put wall	
	while $pos[p2]<$boxvec[p2]
	while $pos[p1]<$boxvec[p1]
		file << format("%5dWALL     P%5d%8.3f%8.3f%8.3f\n",$mol_no,$atom_no,$pos[0],$pos[1],$pos[2])
		$mol_no+=1
		$atom_no+=1
		$pos[p1]+=$between_mol	
		$total_mol_no+=1
	end
	$pos[p1]=0.0
	$pos[p2]+=$between_mol
	end


end

$init_flag=0	# 0: called at first 1:after the first
def check_edge(line)
	if line =~ /^\s+[0-9.-]+[a-zA-Z]+\s+[a-zA-Z]+\s+[0-9.-]+\s+([0-9.-]+)\s+([0-9.-]+)\s+([0-9.-]+)\s+/
		if $init_flag==0
			$maxpos[0]=$1.to_f
			$maxpos[1]=$2.to_f
			$maxpos[2]=$3.to_f
			$minpos[0]=$1.to_f
			$minpos[1]=$2.to_f
			$minpos[2]=$3.to_f
			$init_flag=1
		end
		if $maxpos[0] < $1.to_f
			$maxpos[0] = $1.to_f
		end
		if $maxpos[1] < $2.to_f
			$maxpos[1] = $2.to_f
		end
		if $maxpos[2] < $3.to_f
			$maxpos[2] = $3.to_f
		end
		if $minpos[0] > $1.to_f
			$minpos[0] = $1.to_f
		end
		if $minpos[1] > $2.to_f
			$minpos[1] = $2.to_f
		end
		if $minpos[2] > $3.to_f
			$minpos[2] = $3.to_f
		end
	end
	
end

file1=Array.new
file2=Array.new

# check edge
open("conf.gro.bak","r").each do |line|
	file1 << line
	check_edge(line)
end

# modify data
#open("conf.gro.bak","r").each do |line|
file1.each do |line|
	if line =~ /^\s+([.0-9]+)\s+([.0-9]+)\s+([.0-9]+)/
		$boxvec[0]=$1.to_f
		$boxvec[1]=$2.to_f
		$boxvec[2]=$3.to_f

		putwall(file2)
		file2 << line	
	else
		file2 << line
	end
end


# save file
open("conf.gro","w") do |conf|
file2.each do |line|
	if line =~ /^(\s+)(\d+)$/
		conf.print  $1,$2.to_i+$total_mol_no,"\n"		
	else
		conf.print line
	end
end
end

# modify topol.top
open("topol.top","w") do |out|
open("topol.top.bak","r").each do |line|
	out.print line
end
	out.print "WALL ",$total_mol_no,"\n" 
end
print "maxpos"
p $maxpos
print "minpos"
p $minpos
print "Wall molecular size : ",$total_mol_no,"\n"

--
Appendix トラジェクトリーファイル(traj.trr)を連結する。

シミュレーションをしていくと、なんどか継続してシミュレーションをしていくことにな
る。でもって、そのたびに、トラジェクトリーファイル(traj.trr)が出来上がっていく。
アニメーションを作るときなどは、このトラジェクトリーファイルを連結させたい場合が
ある。そんなときのコマンド。

trjcat -f [traj.trr](複数可) -o [出力ファイル.trr] -settime

たいていは、どのtraj.trrもスタートタイムが0[ps]から始まっているので、2番目以降の
ファイルは開始時間をずらす必要がある。-settimeオプションをつけるとどのくらいから
開始するかを逐一聞いてきてくれる。

コマンド実行の様子は以下のよう。

----
$ trjcat -f traj1.trr traj2.trr -o traj.trr -settime
Enter the new start time (ps) for each file.
There are two special options, both disable sorting:

c (continue) - The start time is taken from the end
of the previous file. Use it when your continuation run
restarts with t=0.

l (last) - The time in this file will be changed the
same amount as in the previous. Use it when the time in the
new run continues from the end of the previous one,
since this takes possible overlap into account.

          File             Current start (ps)  New start (ps)
---------------------------------------------------------
                traj1.trr        0.000 ps          
----

この状態で「New start(ps)」の場所に開始時間を入力してもよし、
「c」を押せば自動的に続きの時間を指定してくれる。
でもって、

----
Summary of files and start times used:

          File                Start time
-----------------------------------------
                traj1.trr        0.000 ps
                traj2.trr        Continue from last file

Reading frame       0 time    0.000   
lasttime -12345

Continue writing frames from traj1.trr t=0 ps, frame=0      
Last frame       1000 time 1000.000    ->  frame   1000 time 1000.000 ps     
Reading frame       0 time    0.000   
lasttime 1000
Reading frame       1 time    4.000   
Continue writing frames from traj2.trr t=1004 ps, frame=1001      
Last frame       1000 time 4000.000    ->  frame   2000 time 5000.000 ps     

Last frame written was 2000, time 5000.000488 ps

gcq#16: "Stop Drinking My Beer !" (The Amps)
----

こんな風に結合してくれる。一つ目のファイルは1000[ps]で終わっているで、
次のtraj2.trrファイルは1004[ps]から開始してくれている。

ファイル名が辞書順になっていれば

trjcat -f *.trr -o out.trr -settime

などというふうにワイルドカードも使える。

同様にener.edr(エネルギーファイル)も

eneconv

コマンドで結合できる。オプションはtrjcatと同じ。

eneconv -f [ener.edr](複数可) -o [出力ファイル.edr] -settime

# でもなんでこっちはenecatじゃないのだろう･･･

--
Appendix 途中でプログラムが止まった時は、、、

ハードディスク障害、メモリ障害、などなどでやむなくプログラムが
止まってしまった場合は、tpbconvコマンドを使ってシミュレーションの続きを
実行できる。


tpbconv -s topol.tpr -f traj.trr -e ener.edr -o newtopol.tpr

これを実行すると止まったところから始めてくれる初期設定ファイル

newtopol.tpr

を作ってくれる。
そのまま(gromppコマンドを実行せずに)、mdrunする。
そうすると、続きを実行してくれる。当然md.logやらtraj.trrは新しく作られるので、
traj.trrはあとで、trjcatコマンドで結合してあげる。

ener.edrファイルはeneconvで結合する。

※ただしこのとき-settimeオプションで聞いてくるときに
　'c(continue)'ではなくて'l(last)'を選択するほうが良いみたい。
　もしくは、-settimeオプションをつけづにtrjcat、eneconvを実行する。

つまり、まとめると、
サーバーなどが途中で止まってgrompp.mdpに設定している時間をすべて終わらないうちに
シミュレーションが止まったときは、tpbconvをつかってあらたなtopol.tprをつくって
シミュレーションを続行する。
こうして出来たいくつかのtraj.trr(トラジェクトリーファイル)は、
trjcatで連結するときに、-settimeオプションナシか-settimeオプションで'l'を選ぶ。

ふつうにシミュレーションを実行しおわるとconfout.groが生成される。
この続きをやろうとするときは、confout.groをconf.groに書き直して、
再びgrompp.mdpファイルをチェックして、gromppコマンドであらたなtopol.tprを
作成してシミュレーションを実行する。
こうして実行したときは、シミュレーション時間が0からまたはじまるので、
trjcatでtraj.trrを連結するときは、-settimeオプションで'c'を選択する。

というふうに私は今回理解したのだが、間違っているかもしれないので、
再びtrjcat、eneconvを使うときは、再度自分自身で確かめてほしい。



----
Referene
GROMACS本家		http://www.gromacs.org/
YONEYAさん講義ページ	http://www.nanolc.jst.go.jp/~yoneya/molsim/
Folding@home		http://latitude.jp/fah/
古明地さんPEACH		http://staff.aist.go.jp/y-komeiji/index_j.html



My勉強の原則
1 書いて覚える!
2 使って覚える!
3 じっくり飛ばさず読む!
4 じっくり考える!
5 あせらずゆっくり一つ一つ打破!
6 分かるところまで降りつづけ、
7 降りたら戻る!
8 疲れたら寝る!