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

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

環境

• macOS Sonomaまたは Linux環境。
• Gromacs 2023.3がインストールされている
• MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜の記事に沿うように書いてあります。

概要

先述の記事で「エネルギー最小化（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 という名前のディレクトリを用意します。

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を用意した状態で、

#!/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平衡化に取り掛かります。とはいってもここから先は高性能な計算機やクラスタコンピュータが必要となります。