Replica exchange molecular dynamics
レプリカ交換分子動力学(Replica Exchange MD, REMD)[1]は、通常のMDシミュレーションのサンプリング効率を改善するためのアルゴリズムである。この方法では、ある温度範囲で異なる温度を持つレプリカを用意し、それぞれで独立なMDシミュレーションを実行する。そして、ある一定時間間隔で、レプリカの温度もしくは原子座標を交換する。結果として、REMDシミュレーションの各レプリカは様々な温度空間を経験し、サンプリング効率が上昇する。REMDシミュレーションを実行するために、サンプリングの段階で以下のようなmdinファイルをremd.mdinとして用意する。
Sampling by REMD &cntrl irest=0, ntx=1, nstlim=500, # レプリカ交換の頻度 numexchg=1000, # シミュレーション中のレプリカ交換の回数 dt=0.002, ntt=3, gamma_ln=1.0, temp0=XXXXX, # この段階では温度は指定しない ig=RANDOM_NUMBER, ntc=2, ntf=2, nscm=1000, ntb=0, igb=5, cut=999.0, rgbmax=999.0, ntpr=100, ntwx=1000, ntwr=100000, / &wt TYPE='END'
このmdinファイルでは、各レプリカの総シミュレーション時間は$nstlim\times numexchg$によって設定される。REMDシミュレーションでは複数の温度を持ったファイルが必要である。それらを用意するために、以下のようなスクリプトをsetup_remd_input.shとして作成する。
#!/bin/bash -f if [ -f groupfile ]; then rm groupfile fi nrep=`wc temperatures.dat | awk '{print $1}'` echo $nrep count=0 for TEMP in `cat temperatures.dat` do let COUNT+=1 REP=`printf "%03d" $COUNT` echo "TEMPERATURE: $TEMP K ==> FILE: remd.mdin.$REP" sed "s/XXXXX/$TEMP/g" remd.mdin > temp sed "s/RANDOM_NUMBER/$RANDOM/g" temp > remd.mdin.$REP echo "-O -rem 0 -i remd.mdin.$REP -o remd.mdout.$REP -c npt.rst -r remd.rst.$REP -x remd.mdcrd.$REP -inf remd.mdinfo.$REP -p test.prmtop" >> remd.groupfile rm -f temp done echo "#" >> groupfile echo "N REPLICAS = $nrep" echo " Done."
このスクリプトはremd.mdinとtemperatures.datを読み込んで、複数のmdinファイルとremd.groupfileを生成する。temperatures.datは例えば、以下のように作成する。
269.5 300.0 334.0 371.8 413.9 460.7 512.9 570.9
この時、最高温度は目的に合わせて十分に大きくとる必要がある。温度間隔は大き過ぎると交換が行われないので注意が必要である。大体の目安としては20~50%程度の確率で交換されていればよいだろう。温度は指数関数的に分布させる方法が一般的である。全レプリカ数が$n_{\rm tot}$で、最高温度が$T_{\rm max}$、最低温度が$T_{\min}$であるような時には、$n+1$番目のレプリカ温度は \begin{equation} T_{n+1} = T_{\rm min}\exp({\rm A}n),\ \ {\rm A} = \frac{\log(T_{\rm max}/T_{\rm min})}{n_{\rm tot}-1} \end{equation} で与えられる。setup_remd_input.shを実行すると、複数のmdinファイルとremd.groupfileが用意されていることを確認しよう。
mpirun -np 並列CPUコア数 $AMBERHOME/bin/pmemd -ng 8 -groupfile remd.groupfile
ここで、-ngではレプリカの数を指定することに注意されたい。シミュレーションが終了したら、最初にレプリカの交換でそれぞれのレプリカがどのような温度空間を経験しているかを確認しよう。例えば以下のようなスクリプトをtemp.shとして作成し、実行することで、datファイルに各レプリカの温度の時間発展を書き出すことができる。
for i in {1..8}; do grep NSTEP remd.mdout.00$i | awk '{print $6, $9}' > temp$i.dat done
また、各mdcrdファイルは連続なトラジェクトリーになっている。通常、ある温度の一定条件での構造集団を解析したいことが多い。その場合には、以下のようなスクリプトによって温度一定条件でのトラジェクトリーを再構成することが出来る(extract_traj.sh)。
cpptraj << EOF parm test.prmtop trajin remd.mdcrd.001 remdtraj remdtrajtemp 300.0 # 300.0 Kにおけるトラジェクトリーを再構成する trajout ファイル名.mdcrd EOF
参考文献
[1] Y. Sugita and Y. Okamoto: Chem. Phys. Lett. 314 (1999) 141.
[2] http://ambermd.org/doc12/Amber16.pdf