Figure 1. PEACHのプログラムモジュールと入出力ファイル
Return
| peach_3.0/demo/BPTI | （タンパク質の計算例） |
| peach_3.0/demo/DNA | （核酸の計算例） |
の二つのディレクトリに入っている。本総説では、これらの実行例に沿ってMD計算の手続きを解説する。
ディレクトリBPTI/の下にはタンパク質の例としてBPTI (Bovine Pancreatic Trypsin Inhibitor)の入力データと30 ps のMDの実行結果、DNA/の下には核酸の例として7 base pairの二本鎖DNAの入力データと30 psのMDの実行結果が、それぞれ入っている。両者の下はほとんど同じディレクトリ構造なので、使用した順に各ディレクトリーの内容を説明する。なお、以下で「分子トポロジー」という場合は、分子の構成原子の種類、電荷、結合情報など、座標と速度を除いた、分子の全情報のことを示す。
| data/ | 入力用三次元データの作成。 |
| MG_dat/ | マグネシウムイオンのトポロジー作成（DNAのみ）。 |
| mktop_v/ | 溶質部分子（タンパク質か核酸）の |
| トポロジー作成（真空中）。 | |
| mktop_bx/ | 箱形の水中の状態での分子トポロジー作成。 |
| heat/ | エネルギー極小化と温度上昇。 |
| mdprod/ | MDの本計算の実行。 |
| analmd/ | 結果の解析。 |
これらのディレクトリ内のファイルの概要は、名前だけで、内容の見当が付くようになっている。以下、* は任意の文字列を表す。
*.in は、ジョブ制御用のデータやフラッグパラメータが書かれた入力ファイル
*.csh は、実行用のC-shell Scriptファイル
*.out は、診断用の出力ファイル
これら以外は、人間ではなく計算機が読むためのファイル
（ただし、すべてASCII文字列であるので、人間も読む気になれば読める）。
従って、基本的に、この実行例を編集して新しい分子の計算を行うためには、
*.in ファイルを計算したい分子に合わせて編集し、
*.cshを実行し、
*.out ファイルが出力されたら、それを読んで、結果を確認する
という作業を行えばよい。
Figure 1に、PEACHの各プログラムモジュールの機能と相関関係を模式的に示した。以下では、個々のモジュールの機能を簡単に説明する。
| MKDBASプログラム | アミノ酸とヌクレオチドのトポロジーデータベースを作成する。 |
| MKMOLプログラム | 分子のトポロジー（結合情報）を作成する。 |
| MKCORプログラム | PDBファイルから三次元座標を読み込み、 |
| イオンや溶媒の付加などを行い、 | |
| RUNMDプログラムで使用するための | |
| 初期座標ファイルを作成する。 | |
| MKPARAプログラム | 力場パラメータを割り当てて、 RUNMDプログラムで使用する |
| ための最終的なトポロジーファイルを作成する。 | |
| RUNMDプログラム | MKPARAで作成したトポロジーファイルと |
| MKCORで作成した座標ファイルを読み込み、 | |
| MDまたはエネルギー極小化計算を実行する。 | |
| ANALMDプログラム | 作成したMDトラジェクトリーファイルを加工したり |
| 解析するための様々なプログラムの集合。 | |
| 以下には、この総説で説明するプログラムのみを示す。 |
| 個別プログラム名 | 入力ファイル | プログラムの機能 |
| ANEN | trajene_* | エネルギーの解析 |
| GYRCRD | partop | 回転半径の解析 |
| trajcrd | ||
| RMSCRD | partop | 初期構造からの根平均二乗 |
| restrt | 変位の解析 | |
| trajcrd | ||
| TMPVEL | partop | 温度分布の解析 |
| trajvel |
| MOSBY プログラム（作成、電総研・上野豊[2]）： RUNMDのトラジェクトリの | |
| アニメーション表示を行う。このプログラムは、PEACH | |
| のモジュールではなく、別途に配布されている。 | |
以上のプログラムモジュールを用いて、MDシミュレーションを行うわけである。以下の章では、MDを始めるまでの準備（２−７章、分子トポロジーファイルと座標ファイルの作成、条件決め、等）、実際のMD計算（８，９章）、結果解析（10章）のそれぞれについて、実行例を説明する。
この総説で説明する peach_3.0/demo/ 下のジョブは、
「ScientificなMD計算を一通り行うための雛形」
として使ってもらうことを想定して作成したため、最終結果を得るまでかなり計算時間がかかる（１，２日）。単にプログラムが正常に動くかどうかだけを確かめたい場合は、
peach_3.0/test/
の下のファイルを利用したほうがよい。
なお、MDの計算手順と解析方法は、研究目的と対象によってそれぞれ異なる。「水素結合パターン」、「分子内の協同運動」、「フォールディング過程」、「溶媒と溶質の相互作用」、等々、研究したい現象により、シミュレーションの手順は変わり得る。だが、本総説では特定の目的を定めず、最大公約数的なMDの手順を示した。ユーザーは、この実行例を参考にして、各自の目的に合った方法を工夫して欲しい。
Figure 1. PEACHのプログラムモジュールと入出力ファイル
2 計算条件の決定
作業に先立って、MDの計算条件を決める（Table 2）。一番重要なのは力場であるが、ここでは最新のAMBER96を使うことにする。溶媒は、周期境界条件を持った箱型モデルを選んだ（文献[1] 2.2節）。クーロン力の計算は、エワルド法で行う（文献[1]3.5.3節）。
条件を決定した後、Protein Data Bank(PDB)から目的分子の三次元座標ファイルを取ってくる。
URL http://www.rcsb.org/pdb/
から入手可能である。今回はBPTIは6PTI.pdbを、DNAは312D.pdbを使うことにした。最近はNMRによる溶液構造も多いが、この二つはＸ線解析による結晶構造である。
以下、BPTIとDNAの２つの実行例を並行して解説していく。ユーザーは、自分の手で計算を行う前に、結果の検証のために、ディレクトリーごとオリジナルを複製して保存しておくことを薦める。
Table 2. demoディレクトリ中の分子シミュレーションの条件
| BPTI | DNA | |||||
|---|---|---|---|---|---|---|
| 初期構造 | 6PTI.pdb | 312D.pdb | ||||
| 溶媒の箱の一辺の長さ(A) | 立方体46×46×46 | 直方体 39×42×48 | ||||
| 原子数 | 溶質 | 892 | 439 | |||
| イオン | 6 (Cl-) | 6 (Mg2+) | ||||
| 溶媒 | 7941 | 6981 | ||||
| 合計 | 8839 | 7426 | ||||
| （以下の条件は、BPTIとDNAに共通） | ||||||
| 計算機 | COMPAQ Alphastataion 600 au (0.6 Gflops) | |||||
| 力場 | AMBER96（溶質とイオン） | |||||
| Flexible SPC Water （溶媒） | ||||||
| アンサンブル | 能勢のカノニカルアンサンブル（温度一定） | |||||
| 時間積分と刻み幅 | Tuckermanの多時間刻み幅法（RESPA） | |||||
| 堅い力（共有結合、結合角、ねじれ角) | 0.25 fs | |||||
| 中間の力（VDW、エワルド実空間） | 2.00 fs | |||||
| 柔らかい力（エワルド波数空間） | 4.00 fs | |||||
| 非共有結合力の計算 | ファンデルワールス（VDW）力 | 12 Åで切断 | ||||
| クーロン力（エワルド） | 実空間 | 12 Åで切断 | ||||
| 波数空間 | Kmax = 8 | |||||
| (2109ベクトル) | ||||||
| （VDW力とエワルド実空間は、40 fs毎に原子ペアリストを作り、 それを利用して、一つのサブルーチンで計算した） | ||||||
Figure 2. BPTI （左、リボン表示）とDNA（右、スティック表示）の結晶構造。グラフィックスプログラムRASMOL[3]を使用。
3 PDB ファイルの加工（data/ ）
以下の章では、それぞれのディレクトリ内のファイル名とその内容を説明する。「入力ファイル」「出力ファイル」等は、単に「入力」「出力」と略す。また、「ある分子を構成する各原子の三次元座標のデータの集まり」のことを、単に「座標」と呼ぶことにする。
| 入力 | ファイルの内容 | |
| （BPTI） | ||
| 6PTI.pdb | PDBから入手したままの、加工前の座標。 | |
| （DNA） | ||
| 312D.pdb | 〃 | |
| 実行シェルスクリプト | ||
| （BPTI, DNA ） | ||
| pdbrenum.csh | pdbrenumプログラムを実行し、 | |
| 生のPDBファイルを、mkcor用に整形する。 | ||
| 出力 | ||
| （BPTI, DNA） | ||
| pdbrenum.log | 実行ログ。 | |
| in.pdb | mkcorで利用するための、加工後（手作業を含む）の | |
| PDB形式の座標ファイル。 | ||
| その他 | ||
| （BPTI） | ||
| README | ファイルの簡単な説明。 | |
シェルスクリプトを実行させるには、
./pdbrenum.csh
とコマンドを入力すればよい。そうすると、pdbrenum.logとin.pdbのファイルが作成される。このpdbrenum.csh の中では、pdbrenumというプログラムが実行され、PDBファイルの残基の番号を付け替えたり、原子の名前を入れ替えたりして、PEACHで読みやすくなるようなるようにする。だが、このプログラムは万能ではない。現時点では、作成されたin.pdbに対し、以下のような手作業が必要である。
Figure 3. タンパク質（左）のＮ末端と核酸の５’末端（右）
Q: どのように水素の座標を作ればいいのか？
A: AMBER94の場合、H1, H2, H3は、Nに、HAはCAに、H5TはO5'Tに、H5'1とH5'2はC5'に、それぞれ、結合している。そこで、その結合している重い原子(N, CA, O5'T、C5')の座標をコピーして、そこから、だいたい1Å離れた場所に、水素がくっつくように、手作業で座標を作る。つまり、x, y, zのどれかの値を１だけ増せばよい。
「こんないい加減な水素の座標でいいのか？」
と思われるだろうが、あとで水素だけエネルギー極小化して最適化するので（5.4節）、実用上は問題ない。
4 マグネシウムイオンのデータ作成（MG_dat/ ）＊DNAのみ
DNAのMDでは、今回のデモでは、系の電荷を中和するためのカウンターイオンとしてマグネシウムイオン（Mg2+）を利用することにした。しかし、AMBER94(96)の力場では、このイオンは標準データベースには含まれていない。そこで、mkdbasプログラムを使って、Mg2+のトポロジーデータベースを作成する。
| 入力 | ファイルの内容 | |
| mg.in | 分子トポロジー作成用入力データ。 | |
| 実行シェルスクリプト | ||
| mkdbas.csh | mkdbasプログラムを実行する。 | |
| 出力 | ||
| mg.out | 診断出力。 | |
| mg.dat | mkmolで使用するための、Mg2+イオンの分子トポロジー。 | |
5 溶質分子の真空中の分子トポロジー作成（mktop_v/ ）
はじめに、真空中の溶質（BPTIまたはDNA）の分子トポロジーファイルを作る（vはvacuumの略)。一連のジョブは、
./MKTOP_V.csh
と打てば実行されるようになっているが（数分）、このシェルスクリプトの中では、
mkmol.csh
mkcor.csh
mkpara.csh
min.csh
mkpdb.csh
の５つのシェルスクリプトを、順番に実行させるようになっている。それぞれの行っているジョブの内容を以下で説明する。このディレクトリ下で作成されるファイルの中で、最終的に後で利用するものは、
vac.pdb
で、水素を最適化した座標ファイルを、PDB形式に直したものである。これを６章のmktop_bxで利用する。
| 入力 | ファイルの内容 | ||
| molin | ジョブ制御用入力ファイル。 | ||
| この中で、シミュレーションを行う対象の | |||
| タンパク質のアミノ配列や核酸の塩基配列を定義する。 | |||
| ../../../data/dbase/amber94/amb94tot.dbas | Amber94力場のトポロジーデータ | ||
| ベース（末端以外のアミノ酸残基と核酸塩基）。 | |||
| ../../../data/dbase/amber94/amb94nt.dbas | 〃（Ｎ末端のアミノ酸残基）。 | ||
| ../../../data/dbase/amber94/amb94ct.dbas | 〃（Ｃ末端のアミノ酸残基）。 | ||
| 実行シェルスクリプト | |||
| mkmol.csh | mkmolプログラムを実行する。 | ||
| 出力 | |||
| molout | 診断出力。 | ||
| moltop | 作成された分子トポロジー。 | ||
| 5.2節でmkcorの入力となる。 | |||
molinのポイント
| 入力 | ファイルの内容 | |
| corin | ジョブ制御用入力ファイル。 | |
| moltop | mkmolで作った分子トポロジーファイル | |
| ../data/in.pdb | 入力用のPDB形式の座標ファイル。 | |
| 実行シェルスクリプト | ||
| mkcor.csh | mkcorプログラムを実行する。 | |
| 出力 | ||
| corout | 診断出力。 | |
| cortop | 作成された分子トポロジー。mkparaの入力となる。 | |
| corcor | 最終的な（runmdの初期入力となる）座標ファイル。 | |
coroutのポイント
| 入力 | ファイルの内容 | ||
| parin | ジョブ制御用入力 | ||
| cortop | mkcorで作成した分子トポロジー。 | ||
| ../../../data/dbase/amber96/parm.list | Amber96力場のパラメータ | ||
| 実行シェルスクリプト | |||
| mkpara.csh | mkparaプログラムを実行する。 | ||
| 出力 | |||
| parout | 診断出力。 | ||
| partop | 最終的な分子トポロジー。 | ||
mkpara.cshのポイント
| 入力 | ファイルの内容 | |
| minin | ジョブ制御用入力。 | |
| partop | mkparaで作成した分子トポロジー。 | |
| corcor | mkcorで作成した分子座標。 | |
| 実行シェルスクリプト | ||
| min.csh | runmdプログラムを実行する。 | |
| 出力 | ||
| minout | 診断出力。 | |
| minrep | 計算の途中経過を示す、報告ファイル。 | |
| その時点での、ステップとエネルギー等が出力される。 | ||
| atomout | 各原子に関する情報。 | |
| min.rst | 極小化計算された座標のファイル。5.5節で使われる。 | |
mininのポイント
| 入力 | ファイルの内容 | |
| partop | mkparaで作られた分子トポロジーファイル | |
| min.rst | 水素を最適化した後の構造 | |
| 実行シェルスクリプト | ||
| mkpdb.csh | mkpdbプログラムを実行する。 | |
| 出力 | ||
| vac.pdb | 水素を最適化した構造を、PDB形式に変換したもの。 | |
| これを、６章のmktop_bxで利用する。 | ||
6 溶質、イオン、溶媒を含んだ分子トポロジーの作成(mktop_bx/)
このディレクトリでは、5章で作った真空中の溶質の分子トポロジーに対して、箱形の水を発生させる。
MKTOP_BX.csh
を実行すれば、数分で全作業が終わる。このスクリプトの中では
mkmol.csh
mkcor.csh
mkpara.csh
mkpdb.csh
の4つのスクリプトが実行される。ここで作成されたファイルの中で、以後で使うものは
partop 分子トポロジーファイル
corcor 座標ファイル
の二つである。
基本的な作業は、5章の場合と同様であるが、溶媒やイオンを発生させる部分が、少し異なる。5章で説明したファイルについては、ファイルの内容は省略した。
| 入力 | ファイルの内容 | |
| molin | ||
| ../../../data/dbase/amber94/amb94tot.dbas | ||
| ../../../data/dbase/amber94/amb94nt.dbas | ||
| ../../../data/dbase/amber94/amb94ct.dbas | ||
| (DNAのみ | ||
| ../MG_dat/mg.dat | マグネシムイオンのトポロジー。mkdbasの出力。 | |
| molinの中で定義する。） | ||
| 実行シェルスクリプト | ||
| mkmol.csh | ||
| 出力 | ||
| molout | ||
| moltop | ||
molinのポイント
| 入力 | ファイルの内容 | |
| corin | ||
| moltop | ||
| ../mktop_v/vac.pdb | 入力用のPDB形式の座標ファイル。ここでは、 | |
| 5章で作った座標を入力する。 | ||
| ../../../data/dbase/fSPC/amber94/wattop | ||
| flexible SPC waterの分子トポロジー。 | ||
| ../../../data/dbase/fSPC/amber94/watcrd | ||
| 水の座標データ。 | ||
| 実行シェルスクリプト | ||
| mkcor.csh | ||
| 出力 | ||
| corout | ||
| cortop | ||
| corcor | runmdで使うための座標ファイル。 | |
corinのポイント
Figure 6.4. Figure 2に示した、BPTI （左）とDNA（右）の結晶構造に対して、水とイオンを発生させたもの。
Q: 分子トポロジーファイルと座標ファイルができたものの、なぜこんなに作成が面倒なのか？ 作業が複雑なことにメリットがあるのか？
A:今は、GUIがついて、簡単にシミュレーション計算が実行できる市販ソフトウェアが多い。PEACHはそれらに比べて原始的であり、ブラックボックスとして使うのには不向きである。ただし、プログラムの中で一体何が行われているのか各段階毎に確認できるし、ユーザーが習熟すれば細かい条件を自在に設定できるメリットがある。
7 エワルド法の条件決め（cond-ewald/）
| 入力 | |
| なし（入力データは、スクリプト中に記載されている） | |
| 実行シェルスクリプト | |
| vewald.csh | |
| 出力 | |
| vewald.log | |
| vewald.out | |
エワルド法は、箱の大きさや要求される精度により、パラメータを変化させる必要がある。vewaldは、精度とパラメータを自動的に算出するプログラムである。ここでは、一応、精度をε=10-3になる条件を採用した。Kmax=8, Rcut=12 Å, η=4.8 Å。ただし、ここで用いた誤差解析方法はかなり厳しい条件を課しているので、これを厳密に守る必要はない（もっと、甘い条件でも精度は保たれる）。以下のEM計算とMD計算では、この条件を使った。
8 エネルギー極小化と昇温 （heat/）
このディレクトリでは、まず系全体のエネルギーを極小化し、ついで、温度を300 Kまで上昇させる。
./HEAT.csh > & HEAT.log &
とジョブを発行する。計算時間が掛かる（約10時間強）ので、上記のようにバックグラウンドで走らせるとよい。ログは、HEAT.logに出力される。なお、以下で、座標や速度の「トラジェクトリー」と言う表現を使うが、これは、「分子を構成する各原子の座標や速度の、時々刻々の値」を記録したファイルのことを示す。
| 入力 | ファイルの内容 | |
| ../mktop_bx/partop | 分子トポロジー。 | |
| ../mktop_bx/corcor | 初期座標。 | |
| min0.in | ジョブ制御入力（EM用）。 | |
| qd0.in | 〃 （QD用、Quenched Dynamics、後述）。 | |
| md*.in | 〃 （昇温MD用）。 | |
| 実行シェルスクリプト | ||
| HEAT.csh | runmdプログラムを実行する。 | |
| 出力 | ||
| min0.out | EMの診断出力。 | |
| qd0.out | QDの診断出力。 | |
| md*.out | 昇温MDの診断出力。 | |
| atomout | 原子に関する様々な情報のまとめ。 | |
| *.rst | それぞれの段階での最終座標（と速度）ファイル。 | |
| *.crd | 座標のトラジェクトリー（溶媒込み）。 | |
| *.crdv | 〃（溶媒なし）。 | |
| *.vel | 速度のトラジェクトリー（溶媒込み）。 | |
| *.velv | 〃（溶媒なし）。 | |
| *.ene_* | さまざまなエネルギー関係の量の時間変化。 | |
| minrep | EMの経過報告。 | |
| qdrep | QDの経過報告。 | |
| mdrep | MDの経過報告。 | |
HEAT.cshを実行すると、
min0 → qd0 → md2 → md4 → md6 → md8 → md10
の順に、ジョブが実行され、EMから始まって最終的に300 Kまで温度が上がるようにしてある。各過程について、以下で詳しく説明する。
| 結合長、結合角、ねじれ角： | 0.1 - 0.5 fs |
| VDW力、クーロン力（エワルド実空間）： | 0.5 - 2.0 fs |
| エワルド波数空間： | 2.0 - 5.0 fs |
程度の範囲で、精度と速度を考慮して選べば良い。
9 本計算 （mdprod/）
本番のMD計算（production run）は、mdprod/の下で、
./MDPROD.csh > & MDPROD.log &
とジョブを実行した。これも時間が掛かる（約一日）ので、上記のように、バックグラウンドで走らせるとよい。ログは、MDPROD.logに出力される。入出力ファイルは、基本的には8章のheatと同じである。
| 入力 | ファイルの内容 | |
| ../mktop_bx/partop | 分子トポロジー。 | |
| ../heat/md10.rst | 時刻10 psでの座標と速度。 | |
| md.in | 制御用入力。md20, md30ともに共通。 | |
| 実行シェルスクリプト | ||
| MDPROD.csh | runmdプログラムを実行し、md20(10-20 ps)と、 | |
| md30 (20-30ps)のジョブを行う。 | ||
| 出力 | ||
| md*.out | MDの診断出力。 | |
| atomout | 原子に関する様々な情報のまとめ。 | |
| md*.rst | それぞれの段階での最終座標と速度ファイル。 | |
| md*.crd | 座標のトラジェクトリー（溶媒込み）。 | |
| md*.crdv | 〃（溶媒なし）。 | |
| md*.vel | 速度のトラジェクトリー（溶媒込み）。 | |
| md*.velv | 〃（溶媒なし）。 | |
| md*.ene_* | さまざまなエネルギー関係の量の時間変化。 | |
| mdrep | MDの経過報告。 | |
Q: runmdを動かすと何種類もの座標ファイルと速度ファイルが作成されるが、なんのためか？
A: md*.rstは、restart用のファイル。つまり、一つのジョブが終了した時点のシミュレーションでの時刻と各原子の座標と速度が出力される。次のジョブでは、これが入力になる。
一方、md*.crd, crdv, vel, velvには、時々刻々の変化が入っている。これらは、トラジェクトリー解析のためのファイルである。今回のデモでは、*.crdと*.velでは1psおきの座標と速度が、また*.crdvと*.velvでは0.2 psおきの座標と速度が、それぞれ出力されている。
なお、*.crdと*.velには溶媒を含む全原子のデータが出力されるのに対し、*.crdvと*.velvには溶質とイオンのデータしか出力されない。溶媒まで含むとデータが巨大になるが、溶質とイオンだけなら、小さいファイルで済むからである。普通は溶質とイオンだけが解析の主眼となる。しかし、解析の主眼が、「溶媒分子の挙動」ならば、*.crdvと*.velvは作る必要はなく、*.crdと*.velに溶媒分子を含むデータを頻繁に出力させるべきである。
10 結果の解析 （analmd/）
このディレクトリでは、8章と9章で作成した出力の解析を行っている。MDの結果の解析は着眼点により様々であるが（文献[1] 6章）、標準的なものを例として挙げておいた。PEACHでは、解析はanalmdと呼ばれるプログラム群で行う。サブディレクトリの名前は、プログラム名と一致させておいた。
| analmd/ | ||
| anen/ | エネルギーの解析 | |
| tmpvel/ | アミノ酸残基/核酸塩基毎の温度の解析 | |
| rmscrd/ | 初期構造からの根平均二乗変位の解析 | |
| gyrcrd/ | 回転半径の解析 |
analmdに属するプログラム群は、PEACHの本体（mkdbas, mkmol, mkcor, mkpara, runmd）と違って、対話式に実行するようになっている（user friendlyとは言い難いが）。だが、今回の実行例では、バッチ式のシェルスクリプトで実行した。入力パラメータの意味を知りたい方は、対話式に実行して欲しい。
ここで例を挙げた４つのプログラムは、runmdモジュールの出力のどれかを読み込んで、それらをまとめて、さらに何かの量を計算する、という形になっている。つまり、以下の種類のファイルが入力される。
| md*.crd | 座標のトラジェクトリー（溶媒込み） |
| md*.crdv | 〃（溶媒なし） |
| md*.vel | 速度のトラジェクトリー（溶媒込み） |
| md*.velv | 〃（溶媒なし） |
| md*.ene_* | さまざまなエネルギー関係の量の時間変化 |
| 入力 | ||
| ../../heat/md2.ene_tot | 0-2 ps 間のエネルギーの時間変化。 | |
| ../../heat/md4.ene_tot | 2-4 ps 〃 | |
| ../../heat/md6.ene_tot | 4-6 ps 〃 | |
| ../../heat/md8.ene_tot | 6-8 ps 〃 | |
| ../../heat/md10.ene_tot | 8-10 ps 〃 | |
| ../../mdprod/md20.ene_tot | 10-20 ps 〃 | |
| ../../mdprod/md30.ene_tot | 20-30 ps 〃 | |
| 実行シェルスクリプト | ||
| md0_30.csh | anenプログラムを実行する。 | |
| 出力 | ||
| md0_30.log | 実行ログ。エネルギーの平均と標準偏差。 | |
| md0_30.ene | 0-30 ps間のエネルギーの時間変化。 | |
Figure 10.1. エネルギーの時間変化
Figure 10.1に、BPTIとDNAそれぞれの系のエネルギーの時間変化を、全エネルギー（total）、運動エネルギー（kinetic）、ポテンシャルエネルギー（potential）に分けて、示した。昇温過程では運動エネルギーが階段的に上昇し、それに連れてポテンシャルエネルギーも上昇するのがわかる。系全体のエネルギーは、昇温がおわった8 ps以降は大体安定だが、わずかに下がる傾向にある。これは、系がまだ緩和の途中であることを示している。
| 入力 | ファイルの内容 | |
| ../../heat/md2.vel | 0-2 ps 間の速度の時間変化。 | |
| ../../heat/md4.vel | 2-4 ps 〃 | |
| ../../heat/md6.vel | 4-6 ps 〃 | |
| ../../heat/md8.vel | 6-8 ps 〃 | |
| ../../heat/md10.vel | 8-10 ps 〃 | |
| ../../mdprod/md20.vel | 10-20 ps 〃 | |
| ../../mdprod/md30.vel | 20-30 ps 〃 | |
| 実行シェルスクリプト | ||
| md0_30.csh | tmpvel プログラムを実行する。 | |
| 出力 | ||
| md0_30.log | 実行ログ。 | |
| md0_30.tmp | 0-30 ps間の溶質、イオン、溶媒の温度の時間変化、 | |
| 平均と標準偏差。 | ||
| md0_30.tmpres | 0-30 ps間の残基毎の温度の時間変化、 | |
| 平均と標準偏差。 | ||
Figure 10.2. 溶質、イオン、溶媒の温度変化
Figure 10.2は、溶質、イオン、溶媒毎に分けて、温度の時間変化を表示した。温度の揺らぎの大きさは、イオン＞溶質＞溶媒の順であるが、これは、原子数が少ない順に揺らぎが大きくなる、ということである。30 psでは、完全に熱平衡には至らないので、三者の平均温度は完全には一致しないが、数100 psのMDを行ったあと平均をとれば、温度は一致するはずである。ただし、クーロン力に対しカットオフを使った場合は、溶質、イオン、溶媒の間に温度格差が生ずる（文献[1] 6.1.2節）。
| 入力 | ||
| ../../mktop_bx/partop | トポロジーファイル。 | |
| ../../mktop_bx/corcor | 初期構造の座標ファイル。 | |
| ../../heat/md2.crdv | 0-2 ps 間の座標の時間変化。 | |
| ../../heat/md4.crdv | 2-4 ps 〃 | |
| ../../heat/md6.crdv | 4-6 ps 〃 | |
| ../../heat/md8.crdv | 6-8 ps 〃 | |
| ../../heat/md10.crdv | 8-10 ps 〃 | |
| ../../mdprod/md20.crdv | 10-20 ps 〃 | |
| ../../mdprod/md30.crdv | 20-30 ps 〃 | |
| 実行シェルスクリプト | ||
| md0_30.csh | ||
| 出力 | ||
| md0_30.log | 実行ログ。 | |
| md0_30.avg | 全体、溶質、イオン、分子、残基、それぞれの | |
| 0-30 ps間のRMSDの自乗平均。 | ||
| md0_30.tot | 0-30 ps間の全体、溶質、イオンのRMSDの時間変化。 | |
| md0_30.mol | 0-30 ps間の分子毎のRMSDの時間変化。 | |
Figure 10.3. 結晶構造からのずれの時間変化
Figure 10.3にRMSDの時間変化を示す。溶質分子全体（ALL）、主鎖（MAIN）、側鎖（SIDE）に分けて表示した。MAINやSIDEは、Main chain、Side chainの略である。
Q: 水素原子を無視するのはなぜか？
A: PDBの結晶構造には水素の座標は入っていない。よって水素原子のずれをRMSDに考慮しても、「結晶構造との比較」にならないため。
Q: BPTIでは、主鎖のほうが側鎖に比べてRMSDが小さいのに対し、DNAではその逆に なっている。なぜか？
A: BPTIをはじめ、タンパク質では側鎖は溶媒に面していることが多く、動き易いのが普通である。一方、DNAの場合、主鎖のリン酸基は溶媒に露出しているのに対し、ここで側鎖として扱われている塩基は、内部で塩基対水素結合を作っている。よって、側鎖のほうが動きにくく、結果として、初期構造からのずれも小さくなった。
| 入力 | ファイルの内容 | ||
| ../../mktop_bx/partop | 分子トポロジーファイル。 | ||
| ../../heat/md2.crdv | 0-2 ps 間の座標の時間変化。 | ||
| ../../heat/md4.crdv | 2-4 ps 〃 | ||
| ../../heat/md6.crdv | 4-6 ps 〃 | ||
| ../../heat/md8.crdv | 6-8 ps 〃 | ||
| ../../heat/md10.crdv | 8-10 ps 〃 | ||
| ../../mdprod/md20.crdv | 10-20 ps 〃 | ||
| ../../mdprod/md30.crdv | 20-30 ps 〃 | ||
| 実行シェルスクリプト | |||
| md0_30.csh | gyrcrdプログラムを実行する。 | ||
| 出力 | |||
| md0_30.log | ログファイル | ||
| md0_30.gyr | 0-30 ps間の溶質全体の回転半径の時間変化、平均、標準偏差 | ||
| （DNAのみ） | |||
| md0_30.mol | 二本の鎖それぞれの、0-30 ps間の回転半径の時間変化、平均、 | ||
| 標準偏差。 | |||
Figure 10.4に回転半径の時間変化を示した。目盛りを小さく取っているので激しい変化に見えるが、実際は、BPTI、DNAともに、0.5 Å以内の揺らぎに収まっている。ただし、両方とも、初期構造に比べてわずかに広がる傾向にある。これは、結晶構造を初期構造に用いたMDではごく普通の現象で、結晶環境から溶液環境に出たことで、分子が広がる傾向を示している。
Figure 10.4. 回転半径の時間変化
11 おわりに
以上、生体分子のMDについて、例を挙げて解説した。計算が終わるのにそれぞれ30時間以上掛かるが、それでもわずか30 psのトラジェクトリが得られるだけである。MDに詳しい方の中には、runmdモジュールを「遅い」と感じられた方も多いと思う。だが、並列化されたrunmdモジュールを利用すればもっと速く計算できるので、並列計算機をお持ちの方は、そちらも試して頂きたい。詳しくは、付録を参照のこと。今日では、今回例に挙げた程度の小さい分子なら、数100 psから10 ns程度のトラジェクトリを作成して解析するのが普通となっている。
実際にMDを行う場合、今回紹介した以外に、様々な工夫が行われている。たとえば、今回はタンパク質やDNA分子と、溶媒とイオンを、同時に平衡化している。だが、特にDNAの場合は、安定した結果を得るために、まずDNAは固定しておいて、溶媒とイオンをEMとMDで最適化させる、という前処理を行うことも多い。このような工夫は、MDの対象と目的に応じて、各自試みていただきたい。
10章では、MDトラジェクトリーの解析作業の基本的な例を挙げた。PEACHのanalmdとmiscのプログラム群を用いれば、これ以外にも様々な解析を行うことができる。残念ながら、ドキュメントや実行例は完全ではないが、各自試みて下さると幸いである。また、もしPEACH対応の解析プログラムをユーザー自身が作成し、それを寄託して下さるならば、望外の喜びである。
産業技術融合領域研究所の寺倉清之主席研究官のご指導ご支援に感謝いたします。
参考文献
[ 1] 古明地勇人, 上林正巳, 長嶋雲兵, 生体分子の分子動力学シミュレーション(1)方法, J. Chem. Software, 6, 1-36 (2000).
[ 2] Ueno, Y., Asai, K., A new plug-in software architecture applied for a portable molecular structure browser, Proceedings, Intelligent systems for molecular biology, 5, 329-332 (1997).
[ 3] Sayle, R., RASMOL 2.6, Glaxo Research and Development, Hertfordshire, U. K. (1995).
付録：PEACHのMPIによる並列化
本文で紹介したrunmdモジュールによるMD計算の結果は、単一CPUを用いたものである。だが、PEACH ver. 3.0では、runmdモジュールは、MPIライブラリを用いて並列化されている。Table Aに、本文のBPTIのジョブを例に取って、MPIにより並列計算を行ったベンチマークの結果を示す。計算条件はTable 2と同様だが、計算機は、工業技術院情報計算センター(TACC)のIBM RS/6000 SP スカラー計算サーバーを用いた。
Table A. 並列計算に要した計算時間 （Table 2のBPTIの計算と同一条件）MD 1 fsあたりの経過時間を秒で表示。カッコ内は並列化効率 (%)
| CPU数 | 1 | 4 | 8 | 16 |
|---|---|---|---|---|
| エワルド実空間 ＋VDW力 の計算 | 3.3 | 0.83 | 0.42 | 0.22 |
| (100) | (98) | (98) | (91) | |
| エワルド波数空間の計算 | 6.4 | 1.60 | 0.81 | 0.41 |
| (100) | (99) | (99) | (97) | |
| MD全体 | 10.2 | 3.02 | 1.81 | 1.23 |
| (100) | (84) | (70) | (50) |
現在並列化されているのは、クーロン力（エワルド法）とVDW力計算用サブルーチンである。これらの部分は、16台の場合でも90%以上の高い効率で並列化されていることがわかる。MD計算全体では、並列化されていない部分（特にペアリスト作成サブルーチン）が目立つようになるため、台数が増えると効率は落ちるが、それでも16台で50%は保っている。現時点では、10台前後のCPUで計算するのが、効率的だろう。
Return