上のチュートリアルに従い、BAR法でメタンが水に解ける際の自由エネルギー変化を計算したので、その手順をメモする。
- 座標ファイルをmethane_water.gro, とポロジーファイルをtopol.topとして保存
- 各プロセスのmdpファイルを保存
- λを変更したmdpファイルを生成するためのスクリプト、write_mdp.plとして保存、Jobを連続して実施するためのスクリプトjob.shを保存
write_mdp.plを実行する
perl write_mdp.pl em_steep.mdp
perl write_mdp.pl nvt.mdp
perl write_mdp.pl npt.mdp
perl write_mdp.pl md.mdp
λを変更したmdpファイルが生成されているので、これをMDPというディレクトリに移しておく
mkdir MDP
mv *.mdp MDP
job.shの必要な箇所を変更する
私の場合は以下を変更した（自分が使うGROMACSのbinの箇所を指定する）
GMX=/usr/local/gromacs/2020.6_GPU/bin
job.shを実行する
chmod u+x job.sh
./job.sh
あとは最後までジョブを実行してくれる。ちなみにλ=0（メタンと水が100%相互作用している）からλ＝1（メタンと水が全く相互作用しない）までλを0.05ずつ変更しながら計算を行っている。ラムダの値は21個で、各Production MDが1nsなので、合計ではざっと20nsの計算となる
md xvg
cp Lambda_*/Production_MD/md*.xvg xvg
cd xvg/
gmx bar -f md*.xvg -o -oi
出力は以下。
Final results in kJ/mol:
point 0 - 1, DG 0.19 +/- 0.01
point 1 - 2, DG 0.17 +/- 0.01
point 2 - 3, DG 0.14 +/- 0.00
point 3 - 4, DG 0.10 +/- 0.01
point 4 - 5, DG 0.06 +/- 0.01
point 5 - 6, DG 0.01 +/- 0.01
point 6 - 7, DG -0.05 +/- 0.01
point 7 - 8, DG -0.15 +/- 0.01
point 8 - 9, DG -0.24 +/- 0.02
point 9 - 10, DG -0.37 +/- 0.01
point 10 - 11, DG -0.55 +/- 0.01
point 11 - 12, DG -0.78 +/- 0.02
point 12 - 13, DG -1.12 +/- 0.02
point 13 - 14, DG -1.43 +/- 0.01
point 14 - 15, DG -1.49 +/- 0.03
point 15 - 16, DG -1.32 +/- 0.02
point 16 - 17, DG -1.03 +/- 0.01
point 17 - 18, DG -0.70 +/- 0.01
point 18 - 19, DG -0.39 +/- 0.00
point 19 - 20, DG -0.12 +/- 0.01
total 0 - 20, DG -9.09 +/- 0.06
となり、自由エネルギー変化は-9.09 +/- 0.06 kJ/mol ということになります。