GROMACSを用いた簡単なエネルギー最小化

このページではGROMACS 2023でのエネルギー最小化（最小点に到達するとは言っていない）についての方法を説明します。

## 環境

  - macOS Sonomaまたは Linux環境。
  - [Gromacs 2023.3がインストールされている](https://qiita.com/Ag_smith/items/2432375399e7d7868138)
  - [MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜](https://qiita.com/Ag_smith/items/2fe20010915ace67ed0a)の記事に沿うように書いてあります。

## 概要

先述の記事で「エネルギー最小化（energy minimization）」が出てくるタイミングがありますが、このエネルギー最小化をするためにはまずGROMACSでトポロジーファイルと初期座標ファイルを作成しておく必要があります。

あらかじめ断っておきますが、**この記事で扱う分子の系はタンパク質です**。簡単な小分子化合物では文字通りの意味での系のエネルギー最小化をすることはまだ容易なのですが、タンパク質という複雑系ではそうではいきません。しかし、それでも「エネルギー最小化」という言葉は使用されます。これは、**分子動力学（MD）法シミュレーションを行うための最初の準備段階として、分子力学法によって系をリラックスさせる**という意味で使われます。

話は戻りまして、タンパク質を含む系のトポロジーファイルと初期座標ファイルが作成できたら、このエネルギー最小化のステップで以下の2つを行います。

1.  200 stepの重原子以外（＝水素原子）のエネルギー最小化
2.  200 stepの全原子のエネルギー最小化

第1段階では水素原子の配向を調整しておき、第2段階で全原子の位置を最適化します。エネルギー最小化の200 stepsの数は自由に変えても構いませんが、だいたいこれくらいで十分です（もしだめだったら大抵の場合は初期構造がむちゃくちゃな構造になっているとか、なにかおかしいことが多いです……）。

GROMACSにおけるエネルギー最小化を含むすべてのMDシミュレーションは、**gmx grompp**コマンドによる前準備と**gmx mdrun**コマンドによる実行の2つから構成されます。

## 実例

以下のようにして `~/Desktop/1LKE/gromacs/minimize` という名前のディレクトリを用意します。

``` bash
cd ~/Desktop/1LKE
mkdir -p gromacs/minimize
cd gromacs/minimize
```

## gmx grompp

> `gmx grompp` (the gromacs preprocessor) reads a molecular topology file, checks the validity of the file, expands the topology from a molecular description to an atomic description. The topology file contains information about molecule types and the number of molecules, the preprocessor copies each molecule as needed. There is no limitation on the number of molecule types. Bonds and bond-angles can be converted into constraints, separately for hydrogens and heavy atoms. Then a coordinate file is read and velocities can be generated from a Maxwellian distribution if requested. `gmx grompp` also reads parameters for `gmx mdrun` (eg. number of MD steps, time step, cut-off), and others such as NEMD parameters, which are corrected so that the net acceleration is zero. Eventually, a binary file is produced that can serve as the sole input file for the MD program.

`gmx grompp`コマンドはGROMACSのMDシミュレーションを実行する前の予備動作となるプログラムです。このコマンドでは主に座標情報とトポロジーファイル、そして**MDの設定ファイル（拡張子.mdp）**を読み込んでそれらに問題がないかをチェックした後で、アウトプットファイル`-o`としてバイナリ形式の拡張子tprのファイルを生成します。エネルギー最小化だけでなく、NVT,　NPT平衡化や本番のProduction Runの前にもこの予備動作が必要となります。

エネルギー最小化ステップにおける設定ファイル（mdpファイル）は次の2つです。 ここに

```
define       = -DPOSRES1000
;
integrator   = steep
nsteps       = 200
;
nstlog       = 1
;
nstlist      = 100
ns_type      = grid
pbc          = xyz
rlist        = 1.0
;
coulombtype  = PME
rcoulomb     = 1.0
;
vdwtype      = cut-off
rvdw         = 1.0
;
optimize_fft = yes
;
constraints  = none
```

という`min1.mdp`と、

```
define       =
;
integrator   = steep
nsteps       = 200
;
nstlog       = 1
;
nstlist      = 100
ns_type      = grid
pbc          = xyz
rlist        = 1.0
;
coulombtype  = PME
rcoulomb     = 1.0
;
vdwtype      = cut-off
rvdw         = 1.0
;
optimize_fft = yes
;
constraints  = none
```

という`min2.mdp`ファイルを用意します。 `min1.mdp`では200 stepの重原子以外のエネルギー最小化を、`min2.mdp`では全体のエネルギー最小化を、それぞれ200 stepずつ行います。

この2つのmdpファイルを用意したら、

    gmx grompp -f min1.mdp -o min1.tpr -p ../../top/1lke.top -c ../../top/1lke.gro -r ../../top/1lke.gro -po md1.out.mdp

というコマンドを実行します。このコマンドの意味は`gmx grompp　-h`のヘルプメッセージに詳細が書かれていますが、`-f`は読み込ませるMDのインプット設定ファイル、`-o`でアウトプットの.tprファイルの指定、`-p`でトポロジーファイルの指定、`-c`で初期座標の指定、`-r`はposition restraintなどを使う時のrestraintの座標が入ったファイル（`-c`と同一のファイルでも良い）、`-po`は`-f`で読み込ませた設定ファイルの確認用ファイル（省略されたデフォルトパラメータも全て記載されている）です。

ここでエラーが発生しなければ`min1.tpr`ファイルが生成されているはずです。エラーが出ているようでしたら、エラーメッセージに応じて該当の箇所を確認してみましょう。

## gmx mdrun

`gmx grompp`がうまくいったら、続いて`gmx mdrun`を実行します。

    gmx mdrun -v -s min1.tpr -o min1.trr -c min1.gro -e min1.edr -g min1.log

これにより、エネルギー最小化が行われている様子が逐一メッセージとして現れるはずです。

    GROMACS:      gmx mdrun, version 2018
    Executable:   /Users/Ag_smith/apps/gmx2018/bin/gmx
    Data prefix:  /Users/Ag_smith/apps/gmx2018
    Working dir:  /Users/Ag_smith/Desktop/1LKE/gromacs/minimize
    Command line:
      gmx mdrun -v -s min1.tpr -o min1.trr -c min1.gro -e min1.edr -g min1.log
    
    Reading file min1.tpr, VERSION 2018 (single precision)
    No option -multi
    Using 1 MPI thread
    Using 4 OpenMP threads
    
    
    Steepest Descents:
       Tolerance (Fmax)   =  1.00000e+01
       Number of steps    =          200
    Step=    0, Dmax= 1.0e-02 nm, Epot= -2.94448e+05 Fmax= 9.98260e+05, atom= 1508
    Step=    1, Dmax= 1.0e-02 nm, Epot= -3.02770e+05 Fmax= 1.18286e+05, atom= 1508
    Step=    2, Dmax= 1.2e-02 nm, Epot= -3.06425e+05 Fmax= 2.11457e+04, atom= 1899
    ..............
    Step=  196, Dmax= 1.7e-02 nm, Epot= -3.80983e+05 Fmax= 1.52034e+04, atom= 647
    Step=  197, Dmax= 2.0e-02 nm, Epot= -3.81036e+05 Fmax= 1.62107e+04, atom= 647
    Step=  199, Dmax= 1.2e-02 nm, Epot= -3.81219e+05 Fmax= 2.65849e+03, atom= 647
    Step=  200, Dmax= 1.5e-02 nm, Epot= -3.81290e+05 Fmax= 2.00573e+04, atom= 647
    
    Energy minimization reached the maximum number of steps before the forces
    reached the requested precision Fmax < 10.
    
    writing lowest energy coordinates.
    
    Steepest Descents did not converge to Fmax < 10 in 201 steps.
    Potential Energy  = -3.8129022e+05
    Maximum force     =  2.0057254e+04 on atom 647
    Norm of force     =  1.6220111e+02
    
    GROMACS reminds you: "All Work and No Play Makes Jack a Dull Boy" (The Shining)

ここで`Epot=`のエネルギー値が徐々に小さくなっているのがわかると思います。これで200ステップ経過すると計算が終了します。正常終了したときには最後にGROMACSの謎名言コレクションから1つランダムにメッセージが選ばれます（今回は映画シャイニングより「All work and no play makes Jack a dull boy.（勉強ばかりしていて遊ばないと子供はバカになる）」）。

## 2段階めのgrompp, mdrun

1段階めのgrompp, mdrunはこれで終わりました。続いて第2段階も同様にやってしまいます。

    gmx grompp -f min2.mdp -po min2.out.mdp -c min1.gro -r ../../top/1lke.gro -p ../../top/1lke.top -o min2.tpr
    gmx mdrun -v -s min2.tpr -o min2.trr -c min2.gro -e min2.edr -g min2.log

問題なく終わっているようでしたら、これでエネルギー最小化ステップは終了です。

## エネルギー最小化のrun.shスクリプト例

上記`min1.mdp`と`min2.mdp`を用意した状態で、

``` bash
#!/bin/sh

# このスクリプトの実行環境は~/Desktop/1LKE/gromacs/minimizeを想定している。
# ~/Desktop/1LKE/top ディレクトリにある1lke.gro, 1lke.topを参照させる。
input="../../top/1lke"

rm -f min1.out.mdp min1.tpr
# gmx grompp (GROMACS PreProcessorコマンド）で、1段階目のmdpファイル、初期座標・トポロジーファイルを読み込み、
# MDシミュレーションを行うために必要なtprファイルを書き出す
gmx grompp \
  -f min1.mdp \
  -po min1.out.mdp \
  -c ${input}.gro \
  -r ${input}.gro \
  -p ${input}.top \
  -o min1.tpr \
  -maxwarn 1 || exit 1
rm -f min1.log min1.gro min1.trr min1.edr
# 上記tprファイルを入力としてMDシミュレーション（エネルギー最小化）を行う
gmx mdrun -v \
  -s min1.tpr \
  -o min1.trr \
  -c min1.gro \
  -e min1.edr \
  -g min1.log || exit 1
rm -f min2.out.mdp min2.tpr
# 2段階目のgrompp。1段階目のgmx mdrunで出力されたmin1.groを初期構造とし、min2.mdpの設定を読み込んで2段階目のエネルギー最小化を行う
gmx grompp \
  -f min2.mdp \
  -po min2.out.mdp \
  -c min1.gro \
  -r ${input}.gro \
  -p ${input}.top \
  -o min2.tpr \
  -maxwarn 1 || exit 1
rm -f min2.log min2.gro min2.trr min2.edr
# 2段階目のエネルギー最小化実行
gmx mdrun -v \
  -s min2.tpr \
  -o min2.trr \
  -c min2.gro \
  -e min2.edr \
  -g min2.log
```

という`run.sh`ファイルを作り、これを実行するだけでエネルギー最小化ができるはず。ちなみに`|| exit 1`は、「処理が正常終了しなかった時にここで打ち切る（そしてコード1を返す）」という意味です。また変数`input`を使って繰り返し書くところを省略しています。シェルスクリプトを書く上ではとても便利なのでこれからも使いましょう。

`grompp`では必ず`-f`（これから行うMDの仕様パラメータを記述した入力ファイル、拡張子`.mdp`）、`-c`（入力座標ファイル）、`-r`（入力restraintファイル）と`-p`（トポロジーファイル）をインプットとして要求し、`-o`と`-po`でアウトプットファイル名を指定します。詳細な仕様は`gmx grompp -h`で確認することができます（他のgmxコマンドでも同様）。

    DESCRIPTION
    （詳しい説明文）
    （中略）
    OPTIONS
    
    Options to specify input files:
    
     -f      [<.mdp>]           (grompp.mdp)
               grompp input file with MD parameters
     -c      [<.gro/.g96/...>]  (conf.gro)
               Structure file: gro g96 pdb brk ent esp tpr
     -r      [<.gro/.g96/...>]  (restraint.gro)  (Opt.)
               Structure file: gro g96 pdb brk ent esp tpr
     -rb     [<.gro/.g96/...>]  (restraint.gro)  (Opt.)
               Structure file: gro g96 pdb brk ent esp tpr
     -n      [<.ndx>]           (index.ndx)      (Opt.)
               Index file
     -p      [<.top>]           (topol.top)
               Topology file
     -t      [<.trr/.cpt/...>]  (traj.trr)       (Opt.)
               Full precision trajectory: trr cpt tng
     -e      [<.edr>]           (ener.edr)       (Opt.)
               Energy file

いくつかの引数指定の部分では、明示されなかった場合には（）のファイル名を自動的にインプットとして認識しようとします。例えば`gmx grompp`のときに必要な`-f`を指定しなかった場合、同じディレクトリ内の`grompp.mdp`ファイルを存在の有無にかかわらず読み込もうとします（なければエラーを返します）。

これが終われば、次は[GROMACSを用いたNVT, NPT平衡化](GROMACSを用いたNVT,_NPT平衡化 "wikilink")に取り掛かります。とはいってもここから先は高性能な計算機やクラスタコンピュータが必要となります。