生体分子の分子動力学シミュレーション (2) 応用

古明地 勇人, 田島 澄恵, 原口 誠, 高橋 伸幸, 上林 正巳, 長嶋 雲兵


Return

目次


  1. 計算条件の決定
  2. PDB ファイルの加工
  3. マグネシウムイオンのデータ作成
  4. 溶質分子の真空中の分子トポロジー作成
  5. 溶質、イオン、溶媒を含んだ分子トポロジーの作成
  6. エワルド法の条件決め
  7. エネルギー極小化と昇温
  8. 本計算
  9. 結果の解析
  10. おわりに
     参考文献
     付録:PEACHのMPIによる並列化

1 序

本総説では、プログラムPEACH (Program for Energetic Analysis of bioCHemical molecules) ver. 3を使って生体高分子の分子動力学(Molecular Dynamics, MD)シミュレーションを行う方法を、実践的に解説する。筆者らは「生体分子の分子動力学シミュレーション (1) 方法」(古明地ら[1])において、PEACHの中に導入されたMDの方法について詳述した。今回は、前回紹介した方法を使ってMDを行う手順を、具体的に紹介したい。
プログラムPEACH Ver. 3.0は、
        http://www.etl.go.jp/~komeiji/
から入手することができる。以下では、UNIX環境下でPEACHソースプログラムの解凍とコンパイルが完了したものとして、話を進める。
PEACHの使用方法は、ソースプログラムと一緒に配布される英文ドキュメントに詳述してあるが、初心者にとって、ドキュメントを読んですぐにPEACHプログラムを実行するのは、なかなか困難である。そこで、一通りプログラムを動かすための実行例を作成し、ソースプログラムと併せて配布している。

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トラジェクトリーファイルを加工したり
解析するための様々なプログラムの集合。
以下には、この総説で説明するプログラムのみを示す。

個別プログラム名入力ファイルプログラムの機能
ANENtrajene_*エネルギーの解析
GYRCRDpartop回転半径の解析
trajcrd
RMSCRDpartop初期構造からの根平均二乗
restrt変位の解析
trajcrd
TMPVELpartop温度分布の解析
trajvel

MOSBY プログラム(作成、電総研・上野豊[2]): RUNMDのトラジェクトリの
                    アニメーション表示を行う。このプログラムは、PEACH
のモジュールではなく、別途に配布されている。

以上のプログラムモジュールを用いて、MDシミュレーションを行うわけである。以下の章では、MDを始めるまでの準備(2−7章、分子トポロジーファイルと座標ファイルの作成、条件決め、等)、実際のMD計算(8,9章)、結果解析(10章)のそれぞれについて、実行例を説明する。
この総説で説明する peach_3.0/demo/ 下のジョブは、
「ScientificなMD計算を一通り行うための雛形」
として使ってもらうことを想定して作成したため、最終結果を得るまでかなり計算時間がかかる(1,2日)。単にプログラムが正常に動くかどうかだけを確かめたい場合は、
        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による溶液構造も多いが、この二つはX線解析による結晶構造である。
以下、BPTIとDNAの2つの実行例を並行して解説していく。ユーザーは、自分の手で計算を行う前に、結果の検証のために、ディレクトリーごとオリジナルを複製して保存しておくことを薦める。

Table 2. demoディレクトリ中の分子シミュレーションの条件
BPTIDNA
初期構造6PTI.pdb312D.pdb
溶媒の箱の一辺の長さ(A)立方体46×46×46直方体 39×42×48
原子数溶質892439
イオン6 (Cl-)6 (Mg2+)
溶媒79416981
合計88397426
(以下の条件は、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.pdbPDBから入手したままの、加工前の座標。
(DNA)
312D.pdb        〃
実行シェルスクリプト
(BPTI, DNA )
pdbrenum.cshpdbrenumプログラムを実行し、
生のPDBファイルを、mkcor用に整形する。
出力
(BPTI, DNA)
pdbrenum.log実行ログ。
in.pdbmkcorで利用するための、加工後(手作業を含む)の
PDB形式の座標ファイル。
その他
(BPTI)
READMEファイルの簡単な説明。

シェルスクリプトを実行させるには、
        ./pdbrenum.csh
とコマンドを入力すればよい。そうすると、pdbrenum.logとin.pdbのファイルが作成される。このpdbrenum.csh の中では、pdbrenumというプログラムが実行され、PDBファイルの残基の番号を付け替えたり、原子の名前を入れ替えたりして、PEACHで読みやすくなるようなるようにする。だが、このプログラムは万能ではない。現時点では、作成されたin.pdbに対し、以下のような手作業が必要である。

  1. 高精度な結晶構造の場合、ひとつの原子に対してA, Bふたつの座標があることがある。ここではとりあえず、Bは消してAのみ残した。
  2. X線結晶解析で得られた構造の場合、普通はPDBに水素原子の座標が入っていない。PEACHのmkcorプログラムは足りない残基や原子を自動発生して補う機能を持っているが、先頭付近の水素(タンパクならH1, H2, H3, HA; DNAならH5T、H5'1, H5'2)は、残念ながら自動発生できない。そこで、人の手で座標を作って付加する必要がある。
  3. AMBER 94と96の力場の場合、DNA5'末端と3'末端の塩基名は、それぞれDT5, DT3等に変更しなければならない。
  4. CYSのうちSSボンドを作るものをCYXに変更する。
  5. HIS(ヒスチジン)がある場合は、HIP (荷電状態によっては、HID, HIE)に変更する
など。


Figure 3. タンパク質(左)のN末端と核酸の5’末端(右)

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のどれかの値を1だけ増せばよい。
「こんないい加減な水素の座標でいいのか?」
と思われるだろうが、あとで水素だけエネルギー極小化して最適化するので(5.4節)、実用上は問題ない。

4 マグネシウムイオンのデータ作成(MG_dat/ )*DNAのみ

DNAのMDでは、今回のデモでは、系の電荷を中和するためのカウンターイオンとしてマグネシウムイオン(Mg2+)を利用することにした。しかし、AMBER94(96)の力場では、このイオンは標準データベースには含まれていない。そこで、mkdbasプログラムを使って、Mg2+のトポロジーデータベースを作成する。

入力ファイルの内容
mg.in分子トポロジー作成用入力データ。
実行シェルスクリプト
mkdbas.cshmkdbasプログラムを実行する。
出力
mg.out診断出力。
mg.datmkmolで使用するための、Mg2+イオンの分子トポロジー。

5 溶質分子の真空中の分子トポロジー作成(mktop_v/ )

はじめに、真空中の溶質(BPTIまたはDNA)の分子トポロジーファイルを作る(vはvacuumの略)。一連のジョブは、
        ./MKTOP_V.csh
と打てば実行されるようになっているが(数分)、このシェルスクリプトの中では、
        mkmol.csh
        mkcor.csh
        mkpara.csh
        min.csh
        mkpdb.csh
の5つのシェルスクリプトを、順番に実行させるようになっている。それぞれの行っているジョブの内容を以下で説明する。このディレクトリ下で作成されるファイルの中で、最終的に後で利用するものは、
        vac.pdb
で、水素を最適化した座標ファイルを、PDB形式に直したものである。これを6章のmktop_bxで利用する。

5. 1 mkmol

溶質分子(タンパク質または核酸)の一次構造(アミノ酸または塩基配列)を入力して、分子トポロジーを作成するモジュールである。

入力ファイルの内容
molinジョブ制御用入力ファイル。
この中で、シミュレーションを行う対象の
タンパク質のアミノ配列や核酸の塩基配列を定義する。
../../../data/dbase/amber94/amb94tot.dbasAmber94力場のトポロジーデータ
ベース(末端以外のアミノ酸残基と核酸塩基)。
../../../data/dbase/amber94/amb94nt.dbas〃(N末端のアミノ酸残基)。
../../../data/dbase/amber94/amb94ct.dbas〃(C末端のアミノ酸残基)。
実行シェルスクリプト
mkmol.cshmkmolプログラムを実行する。
出力
molout診断出力。
moltop作成された分子トポロジー。
5.2節でmkcorの入力となる。

molinのポイント

  1. SS結合のCross link を正しく定義しているかどうか。(今回のデモではBPTIのみ)
    (例:CROSS RES 5 ATOM SG - RES 55 ATOM SGこれは、残基5のSGいう名の原子と、残基55のSGという名の原子の間にCross linkがあることを示している。)
moloutのポイント
  1. 系の全電荷が、整数になっているかどうか(キーワード:TOTAL CHARGE OF THE SYSTEM)。整数でない場合は、何かがおかしい。
  2. プログラムが正常終了したなら、最後に 
            MKMOL: PROGRAM SUCCESSFULLY ENDED
    というメッセージが書かれる。

5. 2 mkcor(BPTI, DNA)

ここでは、moltopと、PDBファイルを読み込み、mkcorプログラムを実行して、座標ファイルを作る

入力ファイルの内容
corinジョブ制御用入力ファイル。
moltopmkmolで作った分子トポロジーファイル
../data/in.pdb入力用のPDB形式の座標ファイル。
実行シェルスクリプト
mkcor.cshmkcorプログラムを実行する。
出力
corout診断出力。
cortop作成された分子トポロジー。mkparaの入力となる。
corcor最終的な(runmdの初期入力となる)座標ファイル。

coroutのポイント

  1. 途中で止まってしまうときは、molinで定義した残基名と、in.pdbでの残基名が一致していないことが多い。
  2. どこかの原子に*INTというフラッグが出ていないか? このフラッグは、「きちんとした座標を読み込めなかったし、それを作成することもできなかった」ということを意味する。もし、出ている場合は、in.pdbをチェック。
  3. プログラムが、正常終了したなら、最後に
            MKCOR: PROGRAM SUCCESSFULLY ENDED
    というメッセージが書かれる。

5. 3 mkpara

以上で作られた分子トポロジーに対し、力場パラメータを割り当てる。

入力ファイルの内容
parinジョブ制御用入力
cortopmkcorで作成した分子トポロジー。
../../../data/dbase/amber96/parm.listAmber96力場のパラメータ
実行シェルスクリプト
mkpara.cshmkparaプログラムを実行する。
出力
parout診断出力。
partop最終的な分子トポロジー。

mkpara.cshのポイント

  1. mkmolので利用したトポロジーデータベースと対応した力場パラメータファイルを使っているか? ここで使ったAmber96の力場パラメータファイルは、Amber94のトポロジーデータベースに対応している。
paroutのポイント
  1. パラメータが足りない、という表示が出たら、parm.listを調べる。
  2. プログラムが、正常終了したなら、最後に
            MKPARA: PROGRAM SUCCESSFULLY ENDED
    というメッセージが書かれる。

5. 4 min

mkcorで作成した座標(corcor)は、水素は内部座標を用いて(あるいは人間が手で)適当に発生したものである。よって、水素だけを、エネルギー極小化計算(Energy Minimization, EM、文献[1] 2.3節)により、最適化しておく。使うモジュールは、runmdである(MDもEMも、同じプログラムモジュールを利用する)。

入力ファイルの内容
mininジョブ制御用入力。
partopmkparaで作成した分子トポロジー。
corcormkcorで作成した分子座標。
実行シェルスクリプト
min.cshrunmdプログラムを実行する。
出力
minout診断出力。
minrep計算の途中経過を示す、報告ファイル。
その時点での、ステップとエネルギー等が出力される。
atomout各原子に関する情報。
min.rst極小化計算された座標のファイル。5.5節で使われる。

mininのポイント

  1. "HYDROGEN" というフラッグを立てること。このフラッグを立てることで、水素だけが最適化される。
minoutのポイント
  1. EM計算の途中でエネルギーが発散していないか?(エネルギーがあるステップで極端に変化していたら要注意)。
  2. プログラムが正常終了したなら、最後に
            ALL FINISHED, BYE-BYE
    というメッセージが書かれる。

5. 5 mkpdb

ここでは、5.4で作成した、水素を最適化した座標をPDB形式に変換する。

入力ファイルの内容
partopmkparaで作られた分子トポロジーファイル
min.rst水素を最適化した後の構造
実行シェルスクリプト
mkpdb.cshmkpdbプログラムを実行する。
出力
vac.pdb水素を最適化した構造を、PDB形式に変換したもの。
これを、6章のmktop_bxで利用する。

6 溶質、イオン、溶媒を含んだ分子トポロジーの作成(mktop_bx/)

このディレクトリでは、5章で作った真空中の溶質の分子トポロジーに対して、箱形の水を発生させる。
        MKTOP_BX.csh
を実行すれば、数分で全作業が終わる。このスクリプトの中では
        mkmol.csh
        mkcor.csh
        mkpara.csh
        mkpdb.csh
の4つのスクリプトが実行される。ここで作成されたファイルの中で、以後で使うものは
        partop    分子トポロジーファイル
        corcor    座標ファイル
の二つである。
基本的な作業は、5章の場合と同様であるが、溶媒やイオンを発生させる部分が、少し異なる。5章で説明したファイルについては、ファイルの内容は省略した。

6. 1 mkmol(BPTI, DNA)

入力ファイルの内容
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のポイント

  1. molin 内で、イオンを定義することが必要。イオンは、5.1で出力されたmoloutのtotal chargeを中和するだけの数を発生させればよい。BPTIの場合は、6つのCl-イオン(CIM)を発生させている。これは、AMBER94のデータベースに標準で入っている。一方、DNAの場合は、6つのMg2+イオン(MG)を発生させているが、これは標準のデータベースにはない。よって、molin内で、mg.datを読むように指定することが必要。
moloutのポイント
  1. 系の全電荷がきちんと中和されているか(TOTAL CHARGE OF THE SYSTEMが、0になっているか)。

6. 2 mkcor

基本的には5.2節と同じだが、ここでは、水を発生させるための設定が必要である。

入力ファイルの内容
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
corcorrunmdで使うための座標ファイル。

corinのポイント

  1. "GENION"フラッグをたてる。これで、イオンの座標を自動発生することができる。
coroutのポイント
  1. イオンの座標ができているか。つまり"GEN"というラベルがついているか。
    別々の例を示すため、BPTIでは立方体の、DNAでは直方体の水を、それぞれ発生させているので注意。
Q: イオンの座標はどういう基準で作られたのか?
A: 溶質+イオン系のクーロンポテンシャルが最小になるような位置を選んで、溶媒分子を置き換える形で発生させている(文献[1] 2.2.2節)。

6. 3 mkpara

入出力とも、5.3節と全く同様なので、省略。
ここで作成されたpartopファイルが、runmdで使用するための最終的なトポロジーファイルである。

6. 4 mkpdb(BPTI,DNA)

partopとcorcorを使って、PDB形式のファイルを作る。このファイルは、MD計算で利用するわけではないが、分子グラフィクスプログラムで分子構造を表示させて、イオンや水などが思ったように発生されているか、チェックするために使う。Figure 6.4に水とイオンを発生させた、最終的な座標ファイル(solv.pdb)をグラフィクスプログラムRASMOL(Sayle [3])で表示した。


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.cshrunmdプログラムを実行する。
出力
min0.outEMの診断出力。
qd0.outQDの診断出力。
md*.out昇温MDの診断出力。
atomout原子に関する様々な情報のまとめ。
*.rstそれぞれの段階での最終座標(と速度)ファイル。
*.crd座標のトラジェクトリー(溶媒込み)。
*.crdv〃(溶媒なし)。
*.vel速度のトラジェクトリー(溶媒込み)。
*.velv〃(溶媒なし)。
*.ene_*さまざまなエネルギー関係の量の時間変化。
minrepEMの経過報告。
qdrepQDの経過報告。
mdrepMDの経過報告。

HEAT.cshを実行すると、
        min0 → qd0 → md2 → md4 → md6 → md8 → md10
の順に、ジョブが実行され、EMから始まって最終的に300 Kまで温度が上がるようにしてある。各過程について、以下で詳しく説明する。

8. 1 min0, qd0

MDの準備としてのEMは、
「初期構造が持っている、不適当な歪みを取り除く」
ため だけ にある。よって、(たとえば基準振動解析を行うときのように)凝る必要はない。今回は、Steepest descent (SD)と、Quenched dynamics (QD)で行った(文献[1] 2.3節)。SDは収束しにくいアルゴリズムなので、収束した/しないに関わらず、ある程度までエネルギーが下がったところで、QDに切り替えてしまう。今回は、200ステップで、次のqd0に切り替えた。QDとは、低温(今回は、0.1 K)でMDを行ってEMの効果を得るための方法。時間積分(文献[1] 4章)の発散を防ぐため時間刻み幅が短い(0.1 fs)。また、系が歪みを持っていた場合に起こる急激な温度変化に対応できるように、Noseの定温アルゴリズム(文献[1] 5.3節)の時間定数も極端に短い(0.1 fs)。さらに、もし時間積分が発散してしまった場合は、仮想熱浴粒子の座標と運動量を0にリセットするようにしてある。

8. 2 md2 - md10

続く昇温過程では、10 ps の間を5段階に分けて10 Kから300 Kまで温度を上げている。昇温過程では、Noseの時間定数は、10 fsと短めにしてある。なお、md2, md4, といったネーミングだが、ここでは、
「ジョブが終了した点でのシミュレーション時刻」
をつけている。たとえば、md2は終了時点で2 psの時刻ということになる。ジョブの名前付けは、それぞれのユーザーが自分のわかりやすいように付ければよい。
今回の実行例での温度の上げ方は多少荒っぽいが、もっと念入りにしたければ「時間を長めに」したり「温度の上げ幅を小刻みに」したりすることで対応できる。
昇温過程でも、それに続く本計算(9章)でも、重心運動は凍結している。また、RESPAによる多重時間刻みを使っている(文献[1] 4.2節)。結合距離、結合角、ねじれ角の力には、0.25 fs、VDW力、エワルド実空間は、2 fs、エワルド逆格子空間は4 fsの時間刻み幅を、それぞれ採用している(Table 2)。
Q: エネルギー極小化と昇温過程において、もっとも重要なポイントは何か?
A: QDを行うことである。これは、説明したとおり、極低温(0.1 - 5 K程度)で、しかも短い時間幅で(0.1 fs程度)でMDを行うことである。わずか数10- 数100ステップ行うだけであるが、これをしないで、いきなり長い時間幅を与えると、温度が急上昇し、時間積分が発散して、ジョブが止まってしまうことが多い。なお、今回は、Noseの定温アルゴリズムでの時定数を、極端に短くすることで、QD温度の上昇を防いだ。もし、これでも温度が急上昇してしまう場合は、Woodcockのスケーリング法で(文献[1] 5.2節)、無理矢理ステップ毎に温度を押さえ込んでしまえばよい。筆者の経験では、QDさえ無事終われば、そのあと時間積分が失敗することは、ほとんどない。
Q: 重心運動を凍結する理由は?
A: 温度一定アルゴリズムを使う場合は、重心運動を凍結しないと、"Flying Ice Cube"問題が生ずる。これは、重心運動にばかり運動エネルギーが分配されてしまい、結果として、系の内部の運動エネルギーが低下して温度が下がってしまう現象(文献[1] 2.5.2節)である。
Q: 多重時間刻み法RESPAの条件は上記でよいのか?
A: ここでは、経験上まず問題ないと思われる条件をとりあえず採用した。ただし、本来は、MDの対象や力場の種類に応じて、様々な時間幅の組み合わせを試すべきである。
目安としては、

結合長、結合角、ねじれ角: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.cshrunmdプログラムを実行し、md20(10-20 ps)と、
md30 (20-30ps)のジョブを行う。
出力
md*.outMDの診断出力。
atomout原子に関する様々な情報のまとめ。
md*.rstそれぞれの段階での最終座標と速度ファイル。
md*.crd座標のトラジェクトリー(溶媒込み)。
md*.crdv〃(溶媒なし)。
md*.vel速度のトラジェクトリー(溶媒込み)。
md*.velv〃(溶媒なし)。
md*.ene_*さまざまなエネルギー関係の量の時間変化。
mdrepMDの経過報告。

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とは言い難いが)。だが、今回の実行例では、バッチ式のシェルスクリプトで実行した。入力パラメータの意味を知りたい方は、対話式に実行して欲しい。
ここで例を挙げた4つのプログラムは、runmdモジュールの出力のどれかを読み込んで、それらをまとめて、さらに何かの量を計算する、という形になっている。つまり、以下の種類のファイルが入力される。

md*.crd座標のトラジェクトリー(溶媒込み)
md*.crdv〃(溶媒なし)
md*.vel速度のトラジェクトリー(溶媒込み)
md*.velv〃(溶媒なし)
md*.ene_*        さまざまなエネルギー関係の量の時間変化

10. 1 anen/ エネルギーの解析

anenは、エネルギーの時間変化のファイルをいくつも読み込んで、それらをまとめ、平均と標準偏差を計算するプログラムである。

入力
../../heat/md2.ene_tot0-2 ps 間のエネルギーの時間変化。
../../heat/md4.ene_tot2-4 ps 〃
../../heat/md6.ene_tot4-6 ps 〃
../../heat/md8.ene_tot6-8 ps 〃
../../heat/md10.ene_tot8-10 ps 〃
../../mdprod/md20.ene_tot10-20 ps 〃
../../mdprod/md30.ene_tot    20-30 ps 〃
実行シェルスクリプト
md0_30.cshanenプログラムを実行する。
出力
md0_30.log実行ログ。エネルギーの平均と標準偏差。
md0_30.ene0-30 ps間のエネルギーの時間変化。


Figure 10.1. エネルギーの時間変化

Figure 10.1に、BPTIとDNAそれぞれの系のエネルギーの時間変化を、全エネルギー(total)、運動エネルギー(kinetic)、ポテンシャルエネルギー(potential)に分けて、示した。昇温過程では運動エネルギーが階段的に上昇し、それに連れてポテンシャルエネルギーも上昇するのがわかる。系全体のエネルギーは、昇温がおわった8 ps以降は大体安定だが、わずかに下がる傾向にある。これは、系がまだ緩和の途中であることを示している。

10. 2 tmpvel/ 温度分布の解析

tmpvelは、速度の時間変化のファイルをいくつも読み込んで、それらをまとめ、さらに残基の重心運動の温度と内部運動の温度を計算するプログラムである(文献[1] 6.1.1節)。

入力ファイルの内容
../../heat/md2.vel0-2 ps 間の速度の時間変化。
../../heat/md4.vel2-4 ps 〃
../../heat/md6.vel4-6 ps 〃
../../heat/md8.vel6-8 ps 〃
../../heat/md10.vel8-10 ps 〃
../../mdprod/md20.vel10-20 ps 〃
../../mdprod/md30.vel    20-30 ps 〃
実行シェルスクリプト
md0_30.cshtmpvel プログラムを実行する。
出力
md0_30.log実行ログ。
md0_30.tmp0-30 ps間の溶質、イオン、溶媒の温度の時間変化、
平均と標準偏差。
md0_30.tmpres0-30 ps間の残基毎の温度の時間変化、
平均と標準偏差。


Figure 10.2. 溶質、イオン、溶媒の温度変化

Figure 10.2は、溶質、イオン、溶媒毎に分けて、温度の時間変化を表示した。温度の揺らぎの大きさは、イオン>溶質>溶媒の順であるが、これは、原子数が少ない順に揺らぎが大きくなる、ということである。30 psでは、完全に熱平衡には至らないので、三者の平均温度は完全には一致しないが、数100 psのMDを行ったあと平均をとれば、温度は一致するはずである。ただし、クーロン力に対しカットオフを使った場合は、溶質、イオン、溶媒の間に温度格差が生ずる(文献[1] 6.1.2節)。

10. 3 rmscrd/ 初期構造からの根平均二乗変位の解析

rmscrdは、初期構造からのMD構造の違い(根平均二乗変位, RMSD)を比較するプログラムである(文献[1] 6.2.1節)。

入力
../../mktop_bx/partopトポロジーファイル。
../../mktop_bx/corcor初期構造の座標ファイル。
../../heat/md2.crdv0-2 ps 間の座標の時間変化。
../../heat/md4.crdv2-4 ps 〃
../../heat/md6.crdv4-6 ps 〃
../../heat/md8.crdv6-8 ps 〃
../../heat/md10.crdv8-10 ps 〃
../../mdprod/md20.crdv10-20 ps 〃
../../mdprod/md30.crdv    20-30 ps 〃
実行シェルスクリプト
md0_30.csh
出力
md0_30.log実行ログ。
md0_30.avg全体、溶質、イオン、分子、残基、それぞれの
0-30 ps間のRMSDの自乗平均。
md0_30.tot0-30 ps間の全体、溶質、イオンのRMSDの時間変化。
md0_30.mol0-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の場合、主鎖のリン酸基は溶媒に露出しているのに対し、ここで側鎖として扱われている塩基は、内部で塩基対水素結合を作っている。よって、側鎖のほうが動きにくく、結果として、初期構造からのずれも小さくなった。

10. 4 gyrcrd/ 回転半径の解析

gyrcrdは、分子の回転半径を計算するプログラムである(文献[1] 6.2.3節)。

入力ファイルの内容
../../mktop_bx/partop分子トポロジーファイル。
../../heat/md2.crdv0-2 ps 間の座標の時間変化。
../../heat/md4.crdv2-4 ps 〃
../../heat/md6.crdv4-6 ps 〃
../../heat/md8.crdv6-8 ps 〃
../../heat/md10.crdv8-10 ps 〃
../../mdprod/md20.crdv10-20 ps 〃
../../mdprod/md30.crdv    20-30 ps 〃
実行シェルスクリプト
md0_30.cshgyrcrdプログラムを実行する。
出力
md0_30.logログファイル
md0_30.gyr0-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数14816
エワルド実空間
    +VDW力 の計算
3.30.830.420.22
(100)(98)(98)(91)
エワルド波数空間の計算6.41.600.810.41
(100)(99)(99)(97)
MD全体10.23.021.811.23
(100)(84)(70)(50)

現在並列化されているのは、クーロン力(エワルド法)とVDW力計算用サブルーチンである。これらの部分は、16台の場合でも90%以上の高い効率で並列化されていることがわかる。MD計算全体では、並列化されていない部分(特にペアリスト作成サブルーチン)が目立つようになるため、台数が増えると効率は落ちるが、それでも16台で50%は保っている。現時点では、10台前後のCPUで計算するのが、効率的だろう。


Return