昇温シミュレーションでは系の温度を挙げて目標となる速度分布となるようにした。しかし、昇温シミュレーションは体積一定条件下で行われていたため、系の分子の密度が1気圧における密度になっていない。そのため、昇温した後にNPTアンサンブル、つまり、圧力一定条件下で数ns〜数十ns程度のシミュレーションを実行し、体積を調節する。短いNPTアンサンブルでのシミュレーションでは、以下のようなmdinファイルをnpt.inとして作成する。

Decide the size of simulation box
 &cntrl
  imin=0, 
  ntx=5, # ntx=5, irest=1を指定することで前回の分子動力学計算の構造と速度を用いてリスタートする
  irest=1,
  nstlim=1000000, # 2 nsの圧力一定シミュレーションを実行する
  dt=0.002,
  ntf=2,
  ntc=2,
  temp0=300.0,
  ntpr=1000,
  ntwx=500,
  cut=12.0,
  ntb=2, # 圧力一定条件下で周期的境界条件を用いる
  ntp=1, # 等方的に圧力制御を行う
  ntt=3,
  gamma_ln=5.0,
  ig=-1,
 /
 $wt
 type='END',
 &end

次に、これまでのチュートリアルと同様にsanderもしくはpmemdを用いて計算を実行する。この時、分子の座標情報は昇温シミュレーションで得たheat.rstから読み込む。ここではpmemdを用いて、以下のように計算を実行する。

mpirun -np 並列CPUコア数 $AMBERHOME/bin/pmemd -O -i npt.in -o npt.out -x npt.mdcrd -p test.prmtop -c heat.rst -r npt.rst -inf npt.info

実際に、簡単に体積の変化を確認してみる。以下のようなBashスクリプトをvolume.shとして作成し、shコマンドで実行する。

grep EKCMT npt.out | awk '{print $9}' > volume.dat
head -n -9 volume.dat > volume_cut.dat # 最後の9行は体積の時間発展以外の値が格納されてしまうので削除する

実際にPythonを使ってvolume_cut.datを描写した図を以下に示す。この図から体積の値が時間経過によって収束していっている様子が確認できる。描写には、Pythonの他にもRやGnuplotを使ってもよいだろう。 体積の時間発展 本チュートリアルでは本計算をNVTアンサンブルで実行することを想定しているが、読者がNPTアンサンブルでシミュレーションを行いたい場合もあるだろう。その場合には、nstlimの値をより大きな値に変更し、シミュレーション後半部分を本計算(サンプリング)としてもよいだろう。


参考文献
http://ambermd.org/doc12/Amber16.pdf