NPTアンサンブルでの体積決定のためのシミュレーションで体積を決定したので、遂に本計算(サンプリング)を行う。これまでと同様にmdinファイルを以下のようにmain.inとして作成する。この段階で、アウトプットファイルをどのくらいの頻度で書き出せばよいかを指定する。この書き出し頻度は少ないと解析に不十分であるし、多過ぎても容量を圧迫してしまう。そのため、本計算の前に、どのような情報が必要かを十分検討して決定することをおすすめする。また、SHAKEアルゴリズムによってボンドの動きを拘束し、Nonbonded相互作用にカットオフ距離を設けていることにも注意されたい。

Production
 &cntrl
  imin=0,
  ntx=5,
  irest=1,
  nstlim=5000000, # 10 nsの分子動力学計算を行う
  dt=0.002, # 2 fsの時間刻みで計算する
  ntf=2,# ntf=2, ntc=2によってSHAKEアルゴリズムを使用する
  ntc=2,
  temp0=300.0,
  ntpr=1000, # mdout形式のアウトプットファイルを指定されたステップ毎に書き出す
  ntwx=500, # mdcrd形式のアウトプットファイルを指定されたステップ毎に書き出す
  ntwr=1000000 # リスタートファイルを2 ns毎に作成する
  cut=12.0, # Nonbonded相互作用のカットオフ距離
  ntb=1,# 体積一定条件下で周期的境界条件を用いる
  ntp=0,
  ntt=3,
  gamma_ln=5.0,
  ig=-1,
 /
 $wt
 type='END',
 &end

今回はpmemdによって以下のように実行する。ここではnpt.rstから分子動力学計算をリスタートする。そのため、体積の値はNPTアンサンブルでの体積決定のためのシミュレーションにおける計算の最後の値を用いることになる。今回は10 nsの計算であるが、昨今の計算化学分野では数百ns〜数十μsの計算が多い。

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

結局、初期構造からMDシミュレーションまでを行うことで、以下のように構造が変化したことが確認できる。 構造変化はポテンシャルエネルギーの変化に関連しており、potential.shのようなスクリプトによって系のポテンシャルエネルギーの時間発展を可視化できる。

grep Etot main.out | awk '{print $9}' > potential.dat
head -n 5000 potential.dat > potential_cut.dat

Pythonを使用して可視化した図を以下に示す(勿論、GnuplotやRを用いてもよい)。 シミュレーション中の構造変化は、Root Mean Square Deviation (RMSD)二次構造二面角といった指標を計算することで、定量的に議論することが出来る。シミュレーションから熱平衡での構造集団を得たい場合には、いくつかの指標を確認して、判断することになるだろう。


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