レプリカ交換分子動力学法(Replica exchange MD, REMD)[1]は、拡張アンサンブル法と呼ばれる方法群に属する手法の一つであり、サンプリング効率を向上される方法である。 目的とする統計力学的集団を保つような条件を用いて温度や圧力、傘ポテンシャルなどの外部パラメータをメトロポリス判定を用いて動的にランダムウォークさせることで、サンプリング効率の向上を行う。この統計力学的集団を保つように外部パラメータを変化させるには、主に事前のシミュレーションにより自由エネルギーの推定が必要となる(このランダムウォークを用いる方法にはマルチカノニカル法や焼戻し法などが挙げられる)。しかし、レプリカ交換法は、系のコピー(ここでいうレプリカ)を複数用いて、式に示されるようなメトロポリス判定式を用いることで自由エネルギーの推定を行わずに外部パラメータのランダムウォークを行うことが出来る。

\begin{equation} p=\min\biggl(1,\exp\biggl\{(E_i-E_j)\biggl(\frac{1}{k_BT_i}\frac{1}{k_BT_j}\biggl)\biggl\}\biggl) \end{equation}

以下はREMD法の模式図である。 また、以下は自由エネルギーの多峰構造を持つ系に対してREMD法を用いた場合の状態Xと自由エネルギーの関係の模式図である。

このようにレプリカが高温部分を経由することで、計算機上の時間内では超えることが出来ない自由エネルギーの多峰構造を超えて空間全域を探索することが可能となる。このサンプリングの向上により物理量をより正確に計算することが出来る。

REMD法はNAMDでは既に実装されている。 REMD法のスクリプトはインストールしたNAMDディレクトリのlib/replica/に存在する。 以下、REMD法で計算を行うために必要なスクリプトの解説を行う。

1.replica.namd

レプリカ交換分子動力学法を行うためのマスタースクリプト。REMD法のコードはここにtcl言語を用いて上から順に実行される形で書かれている。 主に設定ファイルとともに読み込むことで実行される。 そのままでは各温度ではエネルギー最小化も平衡化もされていない。エネルギー極小化については温度は関係ないため、別のスクリプトで用意すると良い。 また、ある一つの温度で平衡化しておくと良い。各温度での平衡化は直接replica.namdに書き込む方法が挙げられる。また、何も変更しない場合、温度は指数関数的に与えられるが、外部ファイルから与えたい場合は同じくreplica.namdのソースコードを変更すると良い。上記の2つとも可能にしたプログラムはサンプルスクリプト置き場に配置しておいた。サンプルスクリプトにも書き込んではおいたが、上記サンプルスクリプトを使用する場合には、sortreplicasを用いて並び替える前にcatdcdで平衡化を行った部分を削除する必要があることを注意すべきである。

2.REMDの設定ファイル

lib/replica/example/fold_alanin.confに相当するファイルである。remdのチュートリアルではアラニンの折りたたみのシミュレーションを行っているため、このような名前になっているが、対象のタンパク質に合わせて好きに命名すると良い。
内容はREMD法を実行する上で必要な設定が記述されたファイルである。

#レプリカ数
set num_replicas 8
#温度の最小値。以下の温度の最大値と合わせて指数関数的に分布させる。
set min_temp 300
#温度の最大値
set max_temp 600
#交換判定までに必要な分子動力学計算のステップ数
set steps_per_run 1000
#交換判定を何回行うか
set num_runs 10000
# dcdファイルに書き込む頻度。交換判定何回毎か
set runs_per_frame 10
# dcdファイルに何回書き込む毎にリスタートファイルを書き出すか
set frames_per_restart 10
# 系の設定ファイル。主に温度のパラメータについて記述をしない設定ファイル。
set namd_config_file "alanin_base.namd"
#出力先のディレクトリの指定。もちろん予め作っておく。
set output_root "output/%s/fold_alanin" ; # directories must exist

系の設定などと同じファイルに記述しても良いが、別のファイルに分けておくとわかりやすいため、このようにしていると思われる。

3.系の設定ファイル

lib/replica/example/alanin_base.namdに相当するファイルである。 本wikiのNAMDのチュートリアルで用いた平衡化計算や昇温計算などのスクリプトからrunやminimizationなどの実行用の命令を取り除き、 温度を記述する項目を取り除いたものを使用すると良いだろう。

4.job0.conf

source fold_alanin.conf

# prevent VMD from reading replica.namd by trying command only NAMD has
if { ! [catch numPes] } { source ../replica.namd }

上記のように記述されたファイル。このファイルはREMDのマスタースクリプトと設定ファイルの呼び出しを行うスクリプトである。
主にこのファイルをバッチスクリプトから呼び出すことでジョブの実行が行われる。

5.バッチスクリプトを実行するシェルスクリプト

REMDは並列計算によって行われる。そして、並列コア数は使用するレプリカ数の倍数を用いる必要がある。 例えば以下のように記述したバッチスクリプトを実行することで計算が行われる。

#!/bin/sh
#恐らくこのあたりにPBSなどのジョブ管理システムの情報が記述される

rm -rf /save/users/hogehoge/output
#outputを削除する
../make_output_dirs.sh /save/users/hogehoge/output 40
#outputディレクトリに40個のファイルを作成する。ここにREMDの結果が吐き出される。

charmrun +p200 namd2 +replicas 40 job0.conf +stdout /save/users/hogehoge/output/%d/job0.%d.log
#200コアを用いて40レプリカの計算を行う。

REMDに限らず系のコピーを用いた計算を行いたい場合は+replicasオプションを記述する必要があることに注意する。


参考文献
[1] Y. Sugita and Y. Okamoto: Chem. Phys. Lett. 314 (1999) 141.