はじめに
以下のチュートリアルを実施したメモ。GROMACSでエタノールの溶媒和自由エネルギーを求める。
下記記事ではメタンの溶媒和自由エネルギーを求めた
実施
前回も似たような内容を行ったが、こちらのほうは分子がエタノールなので、クーロン相互作用も含める必要がある。ただ、このチュートリアルでは計算資源節約のためにVDWとクーロン相互作用を同時に小さくしていく近道をしている。実際にはVDWとクーロンで分けるほうがよい（たしかクーロンを小さくして0にしたあと、VDWを小さくしていく）ということ。
チュートリアルのデータをダウンロードする。
$ wget http://www.gromacs.org/@api/deki/files/261/=gromacs-free-energy-tutorial.tgz
$ tar xzvf \=gromacs-free-energy-tutorial.tgz
$ cd gromacs-free-energy-tutorial/
ethanol0.gro および topol.topがある。ethanolの周りに水分子を置いていく。
# エタノール分子をdodecahedron型のセルに入れる
$ gmx editconf -f ethanol.gro -o box.gro -bt dodecahedron -d 1
# セルに水分子を設置する
$ gmx solvate -cp box.gro -cs -o solvated.gro -p topol.top
エネルギー最小化を実施する。
$ gmx grompp -f em.mdp -c solvated.gro -o em.tpr
$ gmx mdrun -v -deffnm em
なお、コア数が多過ぎてジョブが流れない場合はmdrunのオプションで-nt 4などして並列数を指定するとよい。
続いて系を平衡化する。
$ gmx grompp -f equil.mdp -c em.gro -o equil.tpr
$ gmx mdrun -deffnm equil -v
シェルを実行して各ラムダの値のディレクトリとmdpファイルを生成する。
$ ./mklambdas.sh run.mdp topol.top equil.gro
lambda_00 - lambda_06 のディレクトリが生成されます。
$ grep init-lambda lambda_*/grompp.mdp
lambda_00/grompp.mdp:init-lambda-state = 0
lambda_01/grompp.mdp:init-lambda-state = 1
lambda_02/grompp.mdp:init-lambda-state = 2
lambda_03/grompp.mdp:init-lambda-state = 3
lambda_04/grompp.mdp:init-lambda-state = 4
lambda_05/grompp.mdp:init-lambda-state = 5
lambda_06/grompp.mdp:init-lambda-state = 6
芸がなくてすみませんが、ともかく以下のシェルで各ディレクトリのジョブを実行できます。
各ジョブ200psと短いです。
# !/bin/bash
cd lambda_00
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_01
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_02
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_03
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_04
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_05
gmx grompp
gmx mdrun >& run.log
cd ..
cd lambda_06
gmx grompp
gmx mdrun >& run.log
cd ..
色々オプションが足りない気がするが、topol.topなどデフォルトのファイル名を使っているので、省略されている。
BAR法により自由エネルギーを求める。
$ gmx bar -b 100 -f lambda_*/dhdl.xvg
Final results in kJ/mol:
point 0 - 1, DG -20.49 +/- 0.10
point 1 - 2, DG -14.59 +/- 0.06
point 2 - 3, DG 4.58 +/- 0.12
point 3 - 4, DG 18.32 +/- 1.53
point 4 - 5, DG 4.03 +/- 0.35
point 5 - 6, DG -6.79 +/- 0.69
total 0 - 6, DG -14.95 +/- 2.02
ちなみに-b 100は100psから解析を始める、つまり最初の100psは無視するということ。
チュートリアルでは-20.6 kJ/mol ということだが、私の環境では-14.95となった。
ちなみに実験値は-20.9 kJ/molだそう。