前節では拡張アンサンブル法の一種であるレプリカ交換分子動力学法(REMD法)について取り扱った。レプリカ交換分子動力学法は様々な系について応用されるなど、サンプリング効率向上の手法として極めてパワフルな方法である。しかしながら系の大きさ(自由度)が大きくなるほど温度変化に対する系のエネルギー変化が大きくなるために、レプリカ間の温度刻みを小さくする必要に迫られ、結果としてレプリカ数が増大する。すなわち大きな自由度の系に対して、REMD法は豊富な計算機資源が必要という欠点がある。そこで、この欠点を解決する方法の一つにReplica exchange with solute tempering(REST)法[1]、またその改良版であるREST2法[2]がある。REST法のレプリカ$i$のポテンシャルエネルギー(ハミルトニアン)は次式である。
\begin{equation} E_i^{REST} = a_iE_{\mathrm{solute - solute}}+b_iE_{\mathrm{solute - solvent}} + c_iE_{\mathrm{solvent - solvent}}, \end{equation} \begin{equation} {\rm where}\qquad \left( a_i, b_i, c_i \right)= \begin{cases} \left(1, \frac{\beta_0 + \beta_i}{2\beta_i}, \frac{\beta_0}{\beta_i} \right) & ({\rm for \ REST})\\ \left(\frac{\beta_0}{\beta_i}, \sqrt{\frac{\beta_0}{\beta_i}}, 1 \right) & ({\rm for \ REST2}) \end{cases} \end{equation}

ここで$\beta_0$は系の逆温度、$\beta_i$はレプリカ$i$の逆温度と呼ばれることがあるが$i\gt0$については物理的に意味を持つ逆温度とは異なり、単にポテンシャルエネルギーを指定するパラメータにすぎないことに注意が必要である。この手法では上式のように系のポテンシャルエネルギーを溶質間、溶質-溶媒間、溶媒-溶媒間に分けてそれぞれの係数をレプリカ交換する。特にREST2法の係数について、溶媒-溶媒間の係数が1であることは、自由度の大きい溶媒分子に起因するレプリカ間のエネルギー差を極めて小さくしてレプリカ数の削減に貢献する。また上式から明らかなように$\beta_i=\beta_0$の時に、ポテンシャルエネルギーはbiasのない通常の形になるが、それ以外の$\beta_i\neq\beta_0$では、biasのかかったポテンシャルエネルギーの形をしており、それらのトラジェクトリは基本的に解析に用いない。拡張アンサンブル法に詳しい人はお気づきかもしれないが、この手法はハミルトニアンレプリカ交換法[3,4]の一種である。

REST2法はNAMDで2018年に実装されている[5]。REST2法のスクリプトはインストールしたNAMDディレクトリのlib/replica/REST2/に存在する。 以下、REST2法で計算を行うために必要なスクリプトの解説を行う。注意として2019年1月以降のナイトメア版のNAMD2.13ではREST2が正しく動くことが確認しているが、少なくともそれ以前にインストールされたものに関してはコーディングミスがありREST2について正しい挙動をしないため確認が必要である。

1. rest2_remd.namd

REST2法を行うためのマスタースクリプト。REST2法のコードはここにtcl言語を用いて上から順に実行される形で書かれている。 主に設定ファイルとともに読み込むことで実行される。 また、何も変更しない場合、温度は指数関数的に与えられるが、外部ファイルから与えたい場合は同じくreplica.namdのソースコードを変更すると良い。

2. REST2の設定ファイル

lib/replica/REST2/init.confに相当するファイルである。

# レプリカ数
set num_replicas 16
# REST2のパラメータとしての最低温度
set min_temp 300 ; # physical temperature, no potential energy rescaling 
# REST2のパラメータとしての最大温度
set max_temp 900 ; # highest temperature where potential energy of selected region is rescaled by 1/3 (300/900)
# 物理的な温度(興味のある温度) 基本的に、min_tempと同じ
set TEMP 300
# 交換判定までに必要な分子動力学計算のステップ数
set steps_per_run 100 ; # 0.2 ps
# 交換判定を何回行うか。runs_per_frame * frames_per_restartで割り切れる数を設定
set num_runs 100 ; # 20000000 = 40 ns
# dcdファイルに書き込む頻度。交換判定何回毎か
set runs_per_frame 10 ; # 5 ps per frame
# dcdファイルに何回書き込む毎にリスタートファイルを書き出すか
set frames_per_restart 10 ; # 1000 ps per restart
# 系の設定ファイル
set namd_config_file "aaqaa3_rest2_base.namd"
#出力先のディレクトリの指定。もちろん予め作っておく。
set output_root "output_spt_aaqaa3/%s/rest2"

3. 系の設定ファイル

lib/replica/REST2/aaqaa3_rest2_base.namdに相当するファイルである。REST2特有の設定部分を抜粋しておく。REST2のパラメータについてさらに細かい設定ができるが次のNAMDのリンク先を参考のこと。

# 次の設定はREST2の設定で記入されているので不要
#outputname          rest2_test
#outputEnergies       1000
#outputPressure       1000
#outputTiming         2000
#temperature 300.0

# REST2を使う
soluteScaling           on
# 溶質フラグを立ててあるファイル名
soluteScalingFile       aaqaa3.spt
# 溶質フラグを立ててある列を指定
soluteScalingCol        O

4. job0.conf

REMD法と同じ役割。ファイル名に注意。

source init.conf

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

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

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

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

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

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

参考文献
[1] P. Liu, B. Kim, R. A. Friesner, and B. J. Berne: PNAS 102 (2005) 13749.
[2] L. Wang, R. A. Friesner, and B. J. Berne: J. Phys. Chem. B 115 (2011) 9431.
[3] Y. Sugita, A. Kitao, and Y. Okamoto: J. Chem. Phys. 113 (2000) 6042.
[4] H. Fukunishi, O. Watanabe, and S. Takada: J. Chem. Phys. 116 (2002) 9058.
[5] S. Jo and W. Jiang: Comp. Phys. Commun. 197 (2015) 304.