Gromacsチュートリアル1を日本語化②

* 前回：[Gromacsチュートリアル1を日本語化①](https://qiita.com/shalad2/items/5ffc866395966547fd63)
* Gromacsチュートリアル1：[Lysozyme in Water](http://www.mdtutorials.com/gmx/lysozyme/index.html)

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

1. 系の作成
2. エネルギー最小化
3. NVT平衡化
4. NPT平衡化
5. Production run
6. 分析

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

　エネルギー最小化の入力ファイルは前回イオンの追加でも登場した.tprファイルである。`grompp`を使用して系の構造、トポロジー、シミュレーションパラメータをまとめた.tprファイルを作成する。.mdpファイルは[ここ](http://www.mdtutorials.com/gmx/lysozyme/Files/minim.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<sub>pot</sub>と最大力F<sub>max</sub>の値を確認しておく。まず、ポテンシャルエネルギーE<sub>pot</sub>は負の値であり、系の大きさや水分子の数によるが、10<sup>5 </sup>から10<sup>6 </sup>のオーダーとなる。また、最大力F<sub>max</sub>は1000 kJ mol<sup>-1</sup> nm<sup>-1</sup>以下の値となる。これは.mdpファイルにおいて`emtol = 1000.0`（F<sub>max</sub>が1000 kJ mol<sup>-1</sup> nm<sup>-1</sup>より小さくなるまでエネルギー最小化をする）と設定したためである。

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

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

**10**（Potential）と **0** を選択すると（"10 0"としてEnter）、画面にE<sub>pot</sub>の平均値が表示され、potential.xvgにポテンシャルエネルギーが書き込まれる。potential.xvgをグラフにプロットすると、エネルギー最小化のステップに伴ってポテンシャルエネルギーが安定に収束していることを確認できる。
![potential_energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/693380/c9e45190-5ef6-8d5c-b4ac-71ab56cf9285.png)
これにて系のエネルギーが最小化されたので、実際のダイナミクスに取り掛かる。

## NVT平衡化

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

　NVT平衡化は、エネルギー最小化と同様に`grompp`と`mdrun`を用いて行う。.mdpファイルは[ここ](http://www.mdtutorials.com/gmx/lysozyme/Files/nvt.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ファイルのその他の詳細なパラメータについては、ファイル内のコメントまたは[こちら](http://manual.gromacs.org/documentation/2018/user-guide/mdp-options.html)を参照してほしい。

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

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

グラフにプロットすると、系の温度がすぐに300 Kに達し、平衡化の間で安定化していることを確認できる。
![temperature.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/693380/e2c91e1a-9aba-2acf-3f75-e4e425a94d37.png)


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

　同様に`grompp`と`mdrun`を行う。.mdpファイルは[ここ](http://www.mdtutorials.com/gmx/lysozyme/Files/npt.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計算において大きく変動する量であることを注記しておく。
![pressure.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/693380/55de23c6-61da-3523-88dc-ad3290af0984.png)

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

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

![density.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/693380/3e411473-cf43-333e-5f36-d97d34b98dbf.png)

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

## Production Run

　NVT平衡化、NPT平衡化のステップが完了すると、系は目的の温度と圧力で十分に平衡化されている。ここからはようやく、シミュレーションデータを収集するための**本計算**（**Production run**）を行う。
　`-r`で指定していた平衡化での位置拘束を外し、`-t`でNPT平衡化のチェックポイントファイルを指定して`grompp`と`mdrun`を行う。.mdpファイルは[ここ](http://www.mdtutorials.com/gmx/lysozyme/Files/md.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を行い、分析に使用するシミュレーションデータを得た。

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