前回：Gromacsチュートリアル1を日本語化①
Gromacsチュートリアル1：Lysozyme in Water

　前回は水とイオンの中にタンパク質（lysozyme）が含まれる系を作成した。今回はいよいよ、実際にMD計算を行っていく。
　はじめに、MD計算の手順を軽く紹介する。MD計算では系を作成した後、エネルギー最小化、平衡化を経てようやく本計算（Production run）を行うことになる。今回は次の手順2から手順5までを扱う。

系の作成
エネルギー最小化
NVT平衡化
NPT平衡化
Production run
分析

エネルギー最小化

　系が作成できたら、まずはエネルギー最小化（Energy minimization）を行う。これは、系に立体的な衝突や不適切なジオメトリがないように構造緩和を行う操作である。

　エネルギー最小化の入力ファイルは前回イオンの追加でも登場した.tprファイルである。gromppを使用して系の構造、トポロジー、シミュレーションパラメータをまとめた.tprファイルを作成する。.mdpファイルはここからコピーする。

gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr

前回のイオンの追加では.tprファイルをgenionに渡したが、エネルギー最小化ではmdrunを用いる。-deffnmで入力ファイル名と出力ファイル名を定義している。また、-vでエネルギー最小化のステップを画面表示するように指定している。

gmx mdrun -v -deffnm em

mdrunを実行すると、以下のファイルが得られているだろう。

em.log：エネルギー最小化のログファイル
em.edr：バイナリのエネルギーファイル
em.trr：バイナリのトラジェクトリファイル
em.gro：エネルギー最小化が完了した構造ファイル

　エネルギー最小化が成功しているかを判断するために、ポテンシャルエネルギーE_potと最大力F_maxの値を確認しておく。まず、ポテンシャルエネルギーE_potは負の値であり、系の大きさや水分子の数によるが、10⁵から10⁶のオーダーとなる。また、最大力F_maxは1000 kJ mol^-1 nm^-1以下の値となる。これは.mdpファイルにおいてemtol = 1000.0（F_maxが1000 kJ mol^-1 nm^-1より小さくなるまでエネルギー最小化をする）と設定したためである。

　ここで、エネルギー最小化におけるポテンシャルエネルギーの変化を調べてみる。energyを用いてem.edrからポテンシャルエネルギーのデータを取り出す。

gmx energy -f em.edr -o potential.xvg

10（Potential）と 0 を選択すると（"10 0"としてEnter）、画面にE_potの平均値が表示され、potential.xvgにポテンシャルエネルギーが書き込まれる。potential.xvgをグラフにプロットすると、エネルギー最小化のステップに伴ってポテンシャルエネルギーが安定に収束していることを確認できる。

これにて系のエネルギーが最小化されたので、実際のダイナミクスに取り掛かる。

NVT平衡化

　次のステップはNVT平衡化である。実際のダイナミクスを行うためには、タンパク質周りの溶媒とイオンを平衡化しなければならない。ここではNVTアンサンブルを適用して、300 Kで100 psの平衡化を行う。

　NVT平衡化は、エネルギー最小化と同様にgromppとmdrunを用いて行う。.mdpファイルはここからコピーして使用する。

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

gmx mdrun -deffnm nvt

-rで.groファイルを指定することで、その座標に前回作成したposre.itpを適用し、タンパク質の重たい原子に対して位置拘束のパラメータを付与する。これによって、タンパク質の構造変化に関する変数を追加することなく溶媒の平衡化を行うことができる。
また.mdpファイルでは、タイムステップdt = 0.002 (ps)、ステップ数nsteps = 50000とすることで100 psの平衡化を設定し、gen_vel = yesで各原子にダイナミクス開始時の速度を与えている。さらにtcoupl = V-rescaleおよびpcoupl = noとすることで、圧力制御を行わずに温度制御のみを適用してNVTアンサンブルを実現している。.mdpファイルのその他の詳細なパラメータについては、ファイル内のコメントまたはこちらを参照してほしい。

　NVT平衡化が完了したら、系の温度制御が正しく行われていることを確認する。energyで 16（Temperature）と 0 を選択してtemperature.xvgを作成する。

gmx energy -f nvt.edr -o temperature.xvg

グラフにプロットすると、系の温度がすぐに300 Kに達し、平衡化の間で安定化していることを確認できる。

NPT平衡化

　NVT平衡化により温度が安定したら、続いてNPT平衡化により圧力制御を行う。NPTアンサンブルを用いて、300 K、1 barで100 psの平衡化を行う。

　同様にgromppとmdrunを行う。.mdpファイルはここからコピーして使用する。

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

gmx mdrun -deffnm npt

ここでgromppにおいて、NVT平衡化の時には指定しなかった-tを付け加えている。-tはNVT平衡化からシミュレーションを継続するためのチェックポイントファイル（.cptファイル）を指定しており、NVT平衡化により生成した速度を引き継いでNPT平衡化を行っている。そのため、.mdpファイルにはcontinuation = yesおよびgen_vel = no と設定されている。また、pcoupl = Parrinello-Rahmanの設定によって圧力制御が適用されている。

　NPT平衡化が完了したら、今度は系の圧力制御が正しく行われていることを確認する。energyで 18（Pressure）と 0 を選択してpressure.xvgを作成する。

gmx energy -f npt.edr -o pressure.xvg

グラフにプロットすると、圧力の変動はかなり大きいことが確認できる。基準圧力が1 barであるのに対して、その平均値は7.5±160.5 barとなっている。このように、圧力はMD計算において大きく変動する量であることを注記しておく。

　同様にして、密度変化も確認する。energyで 24（Density）と 0 を選択してdensity.xvgをグラフにプロットする。

gmx energy -f npt.edr -o density.xvg

密度の平均値は1019±3 kg m^-3であり、実験値の1000 kg m^-3やSPC/E水モデルの1008 kg m^-3と近い値になっている。密度の値は100 psの間非常に安定しており、これで系は圧力と密度に関して平衡化されたと判断することができる。

Production Run

　NVT平衡化、NPT平衡化のステップが完了すると、系は目的の温度と圧力で十分に平衡化されている。ここからはようやく、シミュレーションデータを収集するための本計算（Production run）を行う。
　-rで指定していた平衡化での位置拘束を外し、-tでNPT平衡化のチェックポイントファイルを指定してgromppとmdrunを行う。.mdpファイルはここからコピーして使用する。

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr

gmx mdrun -deffnm md_0_1

　計算が終了すると、構造ファイルと1 nsのトラジェクトリファイルが出力される。以降はこれらのファイルを用いて、様々な分析を行うことになる。

まとめ②

　今回はMD計算の一連の流れを見てきた。

エネルギー最小化のステップで構造緩和を行った。
NVT平衡化を行い、温度制御がなされていることを確認した。
NPT平衡化を行い、圧力制御がなされていることを確認した。
Production runを行い、分析に使用するシミュレーションデータを得た。

次回は得られたデータを用いた分析の例を示していく。