本記事では、GROMACSを用いてタンパク質の分子動力学（MD）シミュレーションを行い、安定な構造を評価・抽出する方法について解説します。この構造はin silico screeningに使うことができます。UCFS Chimeraでうまくタンパク質が安定化できなかった場合などに用いることができます。ぜひトライしてみてください！
本記事は以前紹介した以下のGROMACSの記事とほぼ同じですが、各種パラメータと最後の解析方法を加えています。以前に紹介したものはGROMACSのチュートリアルをフォローしているので、ご希望に応じてご参照ください。
【GROMACS】GROMACSを用いたMD simulation【in silico創薬】 – LabCode
動作検証済み環境
Mac M1, Sequoia 15.3.2
自宅でできるin silico創薬の技術書を販売中
新薬探索を試したい方必読！
ITエンジニアである著者の視点から、wetな研究者からもdryの創薬研究をわかりやすく身近に感じられるように解説しています
MD シミュレーションとは？
分子動力学シミュレーション（MDシミュレーション）は、分子の動きや相互作用を計算機上で再現する手法です。分子の構造と相互作用を記述する力場（potential function）を使用し、時間の経過とともに分子の運動を予測します。これにより、物質の性質や相互作用の理解、材料設計、生体分子の研究などに応用されます。
GROMACSは、MDシミュレーションを実行するためのソフトウェアの一つです。GROMACSは高性能な計算を行うための最適化が施されており、広く使われている信頼性の高いツールです。GROMACSを使うことで、分子のダイナミクスや相互作用をシミュレーションすることができます。
この記事では実際にGROMACSを使って、自分のパソコンでMD シミュレーションを行ない、タンパク質が安定化した構造をPDBファイルとして保存します。
GROMACSの環境構築（Homebrew)
GROMACSを使うにあたり、Homebrewをまずターミナルからインストールします。
Homebrewは、macOS向けのパッケージ管理システムです。パッケージ管理システムは、ソフトウェアのインストールやアップデートを簡単に行えるツールであり、開発者やユーザーにとって便利です。
Homebrewを使用すると、ターミナルを介してコマンドを実行することで、様々なソフトウェアパッケージをインストールできます。これにより、開発ツール、プログラミング言語、データベース、ライブラリなど、さまざまなソフトウェアを簡単に管理できます。
Homebrewを使用すると、Mac上でGROMACSをインストールすることができます。Homebrewを使ってGROMACSをインストールすると、依存関係の解決やバージョン管理などが自動的に行われます。
まずターミナルで
/bin/bash -c "$(curl -fsSL <https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh>)"
と打ち、Homebrewをインストールします。
続いて
brew install gromacs
と打ち、gromacsをインストールします。
gmx
と打ち、以下のような記述が返ってきます。
:-) GROMACS - gmx, 2025.1-Homebrew (-:
Executable: /opt/homebrew/bin/../Cellar/gromacs/2025.1/bin/gmx
Data prefix: /opt/homebrew/bin/../Cellar/gromacs/2025.1
Working dir: /Users/kshinba
Command line:
gmx
SYNOPSIS
gmx [-[no]h] [-[no]quiet] [-[no]version] [-[no]copyright] [-nice <int>]
[-[no]backup]
OPTIONS
Other options:
-[no]h (no)
Print help and quit
-[no]quiet (no)
Do not print common startup info or quotes
-[no]version (no)
Print extended version information and quit
-[no]copyright (no)
Print copyright information on startup
-nice <int> (19)
Set the nicelevel (default depends on command)
-[no]backup (yes)
Write backups if output files exist
Additional help is available on the following topics:
commands List of available commands
selections Selection syntax and usage
To access the help, use 'gmx help <topic>'.
For help on a command, use 'gmx help <command>'.
GROMACS reminds you: "If you want to save your child from polio, you can pray or you can inoculate... choose science." (Carl Sagan)これが出てくるときちんとHomebrewとGROMACSのインストールが完了しています。
GROMACSの環境構築（ソースから)
ソースからもGROMACSをインストールできます。GROMACSはCPUのみの環境と、GPUを活用する環境の両方で動作します。それぞれの環境に合わせたインストール方法を説明します。
1. 必要なライブラリのインストール
# 必要なパッケージのインストール
sudo apt-get update
sudo apt-get install -y build-essential cmake wget
build-essential：C/C++の開発に必要なコンパイラやライブラリなど（例：gcc,makeなど）をまとめたパッケージ。cmake：C/C++などのビルド（コンパイルやリンク）を設定・自動化するためのツール。wget：インターネットからファイルをダウンロードするためのコマンドラインツール。
2. GROMACSのインストール（CPU版）
# GROMACS 2025.1のダウンロード
wget <ftp://ftp.gromacs.org/gromacs/gromacs-2025.1.tar.gz>
# アーカイブの展開
tar -xvzf gromacs-2025.1.tar.gz
cd gromacs-2025.1
mkdir build
cd build
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DGMX_MPI=ON -DGMX_OPENMP=ON -DGMX_GPU=OFF -DCMAKE_INSTALL_PREFIX=/usr/local/gromacs
make -j$(sysctl -n hw.logicalcpu)
sudo make install
wget：GROMACSのソースコードをダウンロード。tar -xvzf：ダウンロードしたファイルを解凍。mkdir build：ビルド用のディレクトリを作成。cmake ..：CMakeを使用してコンパイルオプションを設定。make -j$(sysctl -n hw.logicalcpu)：利用可能なCPUコア数を自動的に検出し並列コンパイル。sudo make install：GROMACSをインストール。
GROMACSのビルド時に使用するCMakeオプションについて、それぞれの引数の意味と役割を以下にまとめます。
DGMX_BUILD_OWN_FFTW=ON- 説明: GROMACSは高速フーリエ変換（FFT）を多用します。このオプションを有効にすると、GROMACS自身がFFTWライブラリをソースからダウンロードし、ビルドします。
- 使用目的: システムに適切なFFTWライブラリがインストールされていない場合や、GROMACSと最適な互換性を持つFFTWを使用したい場合に有効です。
- 参考: GROMACS 2024.4 インストールガイド
DGMX_MPI=ON- 説明: GROMACSでMPI（Message Passing Interface）を使用した並列計算を有効にします。
- 使用目的: 複数の計算ノード間での並列計算を行いたい場合に設定します。ただし、シングルノード内での並列計算には不要です。ArchWiki
- 注意: このオプションを有効にするには、OpenMPIやMPICHなどのMPIライブラリが事前にインストールされている必要があります。
- 参考: ArchWikiのGROMACSページArchWiki
DGMX_OPENMP=ON- 説明: GROMACSでOpenMPを使用したマルチスレッド並列計算を有効にします。
- 使用目的: 単一の計算ノード内で、複数のCPUコアを利用して計算速度を向上させる場合に設定します。
- 参考: GROMACS 2024.4 インストールガイド
DGMX_GPU=OFF- 説明: GROMACSでのGPU（Graphics Processing Unit）サポートを無効にします。
- 使用目的: GPUを使用しない、もしくはGPUが搭載されていない環境でのビルド時に設定します。
- 補足: GPUを使用する場合は、対応するCUDAやOpenCLの環境が必要であり、
DGMX_GPU=CUDAやDGMX_GPU=OpenCLと設定します。 - 参考: GROMACS 2024.4 インストールガイド
DCMAKE_INSTALL_PREFIX=/usr/local/gromacs- 説明: GROMACSをインストールするディレクトリを指定します。
- 使用目的: デフォルト以外の場所にGROMACSをインストールしたい場合に設定します。
- 補足: インストール先のディレクトリに書き込み権限が必要です。
- 参考: GROMACSをVisual Studioでビルドする話Tumblr
make -j$(sysctl -n hw.logicalcpu)：利用可能なCPUコア数を自動的に検出し並列コンパイル。sudo make install：GROMACSをインストール。
環境変数を設定します。
echo 'source /usr/local/gromacs/bin/GMXRC' >> ~/.bash_profile
source ~/.bash_profile
source /usr/local/gromacs/bin/GMXRC：GROMACSの環境変数を適用。
動作確認を行います。
gmx --version
3. GPU版を使う場合
以下の記事によると、-DGMX_GPU=OpenCLでGPUが使えるそうです。（私の搭載しているPCはGPUがないので、確認は各自の責任でお願いします。）
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DGMX_MPI=ON -DGMX_OPENMP=ON -DGMX_GPU=OpenCL -DCMAKE_INSTALL_PREFIX=/usr/local/gromacs
参考記事
Mac Studio M2 Ultra に GROMACS 2023.4 を Install – Qiita
Installation guide – GROMACS 2024.4 documentation
以下のコマンドでGPU support: OpenCL と表示されれば、GPUを認識しています。
gmx --version本記事の流れ
- タンパク質の準備: Protein Data Bank（PDB）から糖尿病の標的タンパク質DPP IVの構造データをダウンロードします。この時、結晶の水があるので、除きます。
- PDBからGROMACSトポロジーへの変換: pdb2gmxユーティリティを使用して、タンパク質のトポロジーファイルを作成します。このファイルは、分子の結合、角度、二面角等の情報を含むGROMACSのフォーマットになります。
- ボックス作成と溶媒の追加: editconfユーティリティを使用して、シミュレーションボックスの形とサイズを設定します。またsolvateユーティリティを使用して、シミュレーションボックスに水分子を追加します。
- イオンの追加: システム全体の電荷を中和するためにイオンを追加します。
- エネルギー最小化: 分子の初期配置が不適切な場合、不自然に高い力が生じる可能性があります。これを解決するためにエネルギー最小化（または最適化）ステップを実行します。
- 平衡化（NVTとNPT）: これらはシステムを平衡状態に持っていくためのステップです。NVTは一定の体積と温度で、NPTは一定の圧力と温度でシミュレーションを実行します。
- MDシミュレーション: ここで、最終的な分子動力学シミュレーションを行います。このステップで生成されたデータは後で分析に使用します。
- 可視化：ここでPymolを用いて、MD simulationを可視化していきます。
- RMSDグラフの作成：MDシミュレーションの安定性を評価するために、ルート平均二乗偏差（RMSD）を時間経過とともにプロットします。これにより、シミュレーション中のタンパク質の構造的安定性を確認できます。
- 安定構造の抽出：RMSDが収束した後の代表的な安定構造を抽出するために、
gmx trjconvを使用して最終的な構造を取得します。
1. タンパク質準備
UCFS chimeraで行います。最近ではChimera Xというバージョンが上がっています。
タンパク質準備については過去の記事を参考にしてください。この記事ではタンパク質は糖尿病の標的タンパク質DPPIV（PDBID: 2I78)を使います。
インシリコ創薬の分子ドッキングで必要なタンパク質の準備を論文に沿ってやってみる。【in silico創薬】
拙著
エンジニアとWET研究者のためのコンピュータ創薬入門 ~コンピュータ創薬を体験して、新しい技術を習得しよう~
自宅でできるin silico創薬【新薬探索編】 ~初心者でもできる”おうち創薬”~
では本記事と同タンパク質でタンパク質の準備をしております。よろしければご参照ください。
エネルギー最小化以外の項目を行い、protein.pdbとして保存してください。
2. PDBからGROMACSトポロジーへの変換
PDBからGromacs用のファイルへ変換していきます。ここでは後の低分子-タンパク質とのMDシミュレーションとの力場、水分子のモデルをそろえています。力場はAMBER99SB-ILDN、水分子モデルはTIP3Pで設定しておきます。
ターミナルで、protein.pdbのあるディレクトリに行き、次のコードをターミナルで実行していください。
gmx pdb2gmx -f protein.pdb -o processed.gro -ignh
gmx pdb2gmx- GROMACSの
pdb2gmxコマンドを実行する。 - タンパク質のPDBファイルを処理し、シミュレーション用のトポロジーを作成する。
- GROMACSの
f protein.pdb- 入力ファイルとして
protein.pdbを指定する。 - 解析対象となるタンパク質の3D構造が含まれるPDBファイル。
- 入力ファイルとして
o processed.gro- 出力ファイルとして
processed.groを作成する。 .groファイルはGROMACSの標準的な構造ファイルで、座標情報が含まれる。
- 出力ファイルとして
ignh- PDBファイルに含まれる水素原子を無視する。
- GROMACSが新たに適切な水素原子を追加するため、シミュレーションの安定性が向上する。
以下のように出てくるので、６を指定してください。
:-) GROMACS - gmx pdb2gmx, 2024 (-:
Executable: /usr/local/gromacs/bin/gmx
Data prefix: /usr/local/gromacs
Working dir: /home/shizuku/mizusama/DPP4
Command line:
gmx pdb2gmx -f protein.pdb -o processed.gro -ignh
Select the Force Field:
From '/usr/local/gromacs/share/gromacs/top':
1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
以下概要の説明です。
| 番号 | 力場名 | 特徴 |
|---|---|---|
| 1 | AMBER03 | AMBERファミリーの初期の力場（2003年）。 |
| 2 | AMBER94 | AMBERシリーズの最初の力場（1995年）。 |
| 3 | AMBER96 | AMBER94の改良版（1996年）。 |
| 4 | AMBER99 | AMBER96をベースに改良（2000年）。 |
| 5 | AMBER99SB | AMBER99の改良版で、タンパク質の二次構造予測を改善（2006年）。 |
| 6 | AMBER99SB-ILDN | AMBER99SBに側鎖の動き（Ile, Leu, Asp, Asn）を改良し、より現実的な動きを再現（2010年）。一般的なタンパク質MDシミュレーションに最適。 |
| 7 | AMBERGS | タンパク質のフォールディング研究向け（2002年）。 |
| 8 | CHARMM27 | タンパク質・DNA・脂質を扱える力場。膜タンパク質の解析に最適。 |
| 9 | GROMOS96 43a1 | GROMOS96シリーズの最初の力場。タンパク質・ペプチド向け。 |
| 10 | GROMOS96 43a2 | 43a1の改良版で、アルカンの二面角を改善。有機分子にも適用可能。 |
| 11 | GROMOS96 45a3 | 43a2の発展版で、溶液中のタンパク質の精度を向上。 |
| 12 | GROMOS96 53a5 | 45a3の改良版で、タンパク質とリガンドの相互作用の精度を向上。 |
| 13 | GROMOS96 53a6 | 53a5の改良版で、極性アミノ酸の扱いを改善し、水和エネルギーの計算精度を向上。 |
| 14 | GROMOS96 54a7 | 最新のGROMOS力場で、溶媒との相互作用や二次構造の安定性を改善。 |
| 15 | OPLS-AA/L | 全原子モデル（All-Atom）で、溶媒中のタンパク質シミュレーションに最適（2001年）。 |
AMBER99SB-ILDNは、タンパク質の構造やダイナミクスを高精度に再現できる力場です。 特に、**イソロイシン（Ile）、ロイシン（Leu）、アスパラギン酸（Asp）、アスパラギン（Asn）**の側鎖の動きを改善しており、以下の研究分野で広く用いられています。
- タンパク質のコンフォメーション変化
- 溶液中や膜環境での構造変化の解析に利用される。
- フォールディングシミュレーション
- タンパク質がどのように折りたたまれるかを再現し、変異体の影響を調べる。
- 酵素のダイナミクス解析
- 触媒作用時の構造変化や基質との相互作用を解析する。
- 分子ドッキングやMDシミュレーション
- 創薬研究で、リガンド-タンパク質結合の安定性評価に用いられる。
次は以下のように出るので、1を指定してください。
Using the Amber99sb-ildn force field in directory amber99sb-ildn.ff
Opening force field file /usr/local/gromacs/share/gromacs/top/amber99sb-ildn.ff/watermodels.dat
Select the Water Model:
1: TIP3P TIP 3-point, recommended
2: TIP4P TIP 4-point
3: TIP4P-Ew TIP 4-point optimized with Ewald
4: TIP5P TIP 5-point (see <https://gitlab.com/gromacs/gromacs/-/issues/1348> for issues)
5: SPC simple point charge
6: SPC/E extended simple point charge
7: None
| 番号 | 水モデル | 特徴 |
|---|---|---|
| 1 | TIP3P | 3点モデル、広く使われる標準モデル、AMBER推奨。 |
| 2 | TIP4P | 4点モデル、TIP3Pより物理特性が改善。 |
| 3 | TIP4P-Ew | Ewald法に最適化、精度向上。 |
| 4 | TIP5P | 5点モデル、水素結合ネットワークが改善。 |
| 5 | SPC | 古典的3点モデル、計算負荷が低い。 |
| 6 | SPC/E | SPCの改良版、水の誘電特性が向上。 |
| 7 | None | 水を含まないシミュレーション（真空中）。 |
TIP3PはAMBER力場で推奨されており、最も一般的に使用される水モデルです。計算負荷が低く、安定した結果を得やすいのが特徴です。特にこだわりがなければTIP3P（1番）を選択するのが無難です。
3. ボックス作成と溶媒の追加
gmx editconf -f processed.gro -o boxed.gro -c -d 1.0 -bt cubic
分子を中央に配置し、シミュレーションボックスを作成します。
f processed.gro→ 入力ファイル（processed.gro）を指定。o boxed.gro→ 出力ファイル（boxed.gro）を作成。c→ 分子をボックスの中央に配置。d 1.0→ 分子の周囲に1.0 nmの間隔を確保。bt cubic→ **立方体（cubic）**のボックス形状を指定。
gmx solvate -cp boxed.gro -cs spc216.gro -o solvated.gro -p topol.top
cp boxed.gro→ 入力ファイル（boxed.gro）を指定。cs spc216.gro→ SPC216（水のクラスター）を使用。o solvated.gro→ 出力ファイル（水を加えたsolvated.gro）。p topol.top→ トポロジーファイル（topol.top）を更新し、水分子の情報を追加。
4. イオンの追加（中性化）
gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr
エネルギー最小化の前に、シミュレーション用のバイナリ入力ファイル（.tpr）を作成する。
f ions.mdp→ パラメータファイル（ions.mdp）を指定。c solvated.gro→ 入力構造ファイル（溶媒を加えたsolvated.gro）を指定。p topol.top→ トポロジーファイル（topol.top）を指定。o ions.tpr→ 出力ファイル（イオン追加用のions.tpr）を作成。
ions.mdpファイルを空ファイルで作成し、内容は以下のコピペで大丈夫です。こちらはGromacsチュートリアルから入手できるのもので、パラメータはほぼ同じに設定しています。
ions.mdp
; LINES STARTING WITH ';' ARE COMMENTS
title = Minimization ; Title of run
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 10.0 kJ/mol
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.0 ; Cut-off for making neighbor list (short range forces)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; long range electrostatic cut-off
rvdw = 1.0 ; long range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
続いて、以下を実行してください。
gmx genion -s ions.tpr -o ionized.gro -p topol.top -pname NA -nname CL -neutral
水の一部をNa⁺（ナトリウム）とCl⁻（塩化物イオン）に置換し、システムを中性化します。
s ions.tpr→ 入力ファイル（ions.tpr）を指定。o ionized.gro→ 出力ファイル（イオンを追加したionized.gro）を作成。p topol.top→ トポロジーファイルを更新し、イオン情報を反映。pname NA→ 追加する陽イオンの名前を「NA（ナトリウム）」に指定。nname CL→ 追加する陰イオンの名前を「CL（塩化物イオン）」に指定。neutral→ システムの総電荷が0になるようにイオンを追加。
5. エネルギー最小化
gmx grompp -f minim.mdp -c ionized.gro -p topol.top -o em.tpr -maxwarn 1
エネルギー最小化のためのバイナリ入力ファイル（.tpr）を作成する。
f minim.mdp→ エネルギー最小化のパラメータファイル（minim.mdp）を指定。c ionized.gro→ イオンを追加した構造ファイル（ionized.gro）を使用。p topol.top→ トポロジーファイル（topol.top）を指定。o em.tpr→ エネルギー最小化用のバイナリファイル（em.tpr）を作成。maxwarn 1→ 警告を1件まで許容し、強制的に実行。（ただし、多用は非推奨）
minim.mdpファイルは以下のコピペで大丈夫です。こちらはGromacsチュートリアルから入手できるのもので、パラメータはほぼ同じに設定しています。
minim.mdp
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 5000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
続いて以下を実行してください。
gmx mdrun -v -deffnm em
エネルギー最小化を実行し、分子の構造を安定化させます。
v→ 詳細な実行ログを表示（verboseモード）。deffnm em→ すべての出力ファイルの接頭辞を「em」に設定。
6. 平衡化（NVT & NPT）
エネルギー最小化（EM）が完了したら、NVT（定積・等温）平衡化とNPT（定圧・等温）平衡化を実行し、シミュレーションを安定させます。
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
このコードはエネルギー最小化後の構造（em.gro）を用いて、NVTシミュレーション用の.tprファイルを作成します。
f nvt.mdp→ NVTシミュレーションのパラメータファイル（nvt.mdp）を指定。c em.gro→ エネルギー最小化後の構造ファイルを使用。r em.gro→ 拘束（position restraints）を適用する基準構造。p topol.top→ トポロジーファイルを指定。o nvt.tpr→ NVT用のバイナリファイル（nvt.tpr）を作成。
nvt.mdpファイルは以下のコピペで大丈夫です。こちらはGromacsチュートリアルから入手できるのもので、パラメータはほぼ同じに設定しています。
nvt.mdp
title = protein_nvt
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
gmx mdrun -deffnm nvt -v
NVTシミュレーションを実行し、システムを目標温度（例: 300K）に安定化させます。
deffnm nvt→ 出力ファイルの接頭辞を「nvt」に設定。v→ 実行状況を詳細に表示（verboseモード）。
終了したら、以下を実行してください。
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
NVT平衡化後の構造（nvt.gro）を用いて、NPTシミュレーション用の.tprファイルを作成。
f npt.mdp→ NPTシミュレーションのパラメータファイル（npt.mdp）を指定。c nvt.gro→ NVT平衡化後の構造ファイルを使用。r nvt.gro→ 拘束（position restraints）を適用する基準構造。t nvt.cpt→ NVTのチェックポイント（温度・速度情報）を継続。p topol.top→ トポロジーファイルを指定。o npt.tpr→ NPT用のバイナリファイル（npt.tpr）を作成。
npt.mdpファイルは以下のコピペで大丈夫です。こちらはGromacsチュートリアルから入手できるのもので、パラメータはほぼ同じに設定しています。
npt.mdp
title = protein_npt
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off
gmx mdrun -deffnm npt -v
NPTシミュレーションを実行し、システムの圧力を目標値（例: 1 atm）に安定化する。
deffnm npt→ 出力ファイルの接頭辞を「npt」に設定。v→ 実行状況を詳細に表示（verboseモード）。
7. MDシミュレーション
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr
NPT平衡化後のデータをもとに、本番MDシミュレーション用の.tprファイルを作成。
f md.mdp→ MDシミュレーションのパラメータファイル（md.mdp）を指定。c npt.gro→ NPT平衡化後の構造ファイル（npt.gro）を使用。t npt.cpt→ NPTのチェックポイント（npt.cpt）を引き継ぎ、温度・速度情報を保持。p topol.top→ トポロジーファイルを指定。o md_0_1.tpr→ 本番MDシミュレーション用のバイナリファイルを作成。
md.mdpファイルは以下のコピペで大丈夫です。こちらはGromacsチュートリアルから入手できるのもので、パラメータはほぼ同じに設定しています。
md.mdp
title = protein_MD
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 2500000 ; 2 * 2500000 = 5000 ps (1 ns)
dt = 0.002 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
gmx mdrun -v -deffnm md_0_1
本番MDシミュレーションを開始し、分子の動きを解析します。
v→ 実行状況を詳細に表示（verboseモード）。deffnm md_0_1→ 出力ファイルの接頭辞を「md_0_1」に設定。
解析が終了すると、以下のファイルが作成しています。
| ファイル名 | 内容 |
|---|---|
| md_0_1.tpr | 本番MDシミュレーション用の入力ファイル。 |
| md_0_1.gro | MDシミュレーション後の最終構造。 |
| md_0_1.trr | 全原子の座標・速度・力の情報（フル精度）。 |
| md_0_1.xtc | 圧縮版のトラジェクトリ（軽量で解析用）。 |
| md_0_1.edr | エネルギー情報（温度、圧力、エネルギー）。 |
| md_0_1.log | 実行ログ（MDの詳細情報）。 |
| md_0_1.cpt | チェックポイント（中断・再開に使用）。 |
8. 可視化
終了したら、md_0_1.groとmd_0_1.xtcファイルが出ていると思います。まずmd_0_1.gro をpymolで開き、開いた状態で、md_0_1.xtc をドラックアンドドロップすると、次のような画面が出てきます。この状態になると、loadを教えてください。
Hide→watersから水を消し、右下のGlobal Framesにある再生ボタンを押すと、MD simulationが開始できます。
以下のように再生できます。
9. RMSD解析
論文に乗せるのはRMSD解析したグラフなので、RMSDのグラフを以下で作成します。
gmx rms -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg
MDシミュレーション中のタンパク質の構造変化（RMSD）を解析します。
s md_0_1.tpr→ リファレンス構造（基準となる構造データ）を含む.tprファイルを指定。f md_0_1.xtc→ トラジェクトリファイル（シミュレーション中の座標データ）を指定。o rmsd.xvg→ RMSDデータを格納する出力ファイル（rmsd.xvg）を作成。
以下のように出力されます。
4を選択してください。
:-) GROMACS - gmx rms, 2024 (-:
Executable: /usr/local/gromacs/bin/gmx
Data prefix: /usr/local/gromacs
Working dir: /home/shizuku/mizusama/DPP4
Command line:
gmx rms -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg
Reading file md_0_1.tpr, VERSION 2024 (single precision)
Reading file md_0_1.tpr, VERSION 2024 (single precision)
Select group for least squares fit
Group 0 ( System) has 327482 elements
Group 1 ( Protein) has 23242 elements
Group 2 ( Protein-H) has 11898 elements
Group 3 ( C-alpha) has 1452 elements
Group 4 ( Backbone) has 4356 elements
Group 5 ( MainChain) has 5810 elements
Group 6 ( MainChain+Cb) has 7182 elements
Group 7 ( MainChain+H) has 7214 elements
Group 8 ( SideChain) has 16028 elements
Group 9 ( SideChain-H) has 6088 elements
Group 10 ( Prot-Masses) has 23242 elements
Group 11 ( non-Protein) has 304240 elements
Group 12 ( Water) has 304212 elements
Group 13 ( SOL) has 304212 elements
Group 14 ( non-Water) has 23270 elements
Group 15 ( Ion) has 28 elements
Group 16 ( Water_and_ions) has 304240 elements
Select a group:
次に以下で出てくるので4を選択してください。
Select group for RMSD calculation
Group 0 ( System) has 327482 elements
Group 1 ( Protein) has 23242 elements
Group 2 ( Protein-H) has 11898 elements
Group 3 ( C-alpha) has 1452 elements
Group 4 ( Backbone) has 4356 elements
Group 5 ( MainChain) has 5810 elements
Group 6 ( MainChain+Cb) has 7182 elements
Group 7 ( MainChain+H) has 7214 elements
Group 8 ( SideChain) has 16028 elements
Group 9 ( SideChain-H) has 6088 elements
Group 10 ( Prot-Masses) has 23242 elements
Group 11 ( non-Protein) has 304240 elements
Group 12 ( Water) has 304212 elements
Group 13 ( SOL) has 304212 elements
Group 14 ( non-Water) has 23270 elements
Group 15 ( Ion) has 28 elements
Group 16 ( Water_and_ions) has 304240 elements
Select a group:
作成されるrmsd.xvgのファイルを見ると、以下のようになっています。
１行目の数字が時間単位がTime (ps)、2行目がRMSDとなっています。RMSD（nm)となっているので、1/10にして、Åの単位に直し、エクセルなどで、グラフを作成してください。
# This file was created Sun Mar 16 09:24:15 2025
# Created by:
# :-) GROMACS - gmx rms, 2024 (-:
#
# Executable: /usr/local/gromacs/bin/gmx
# Data prefix: /usr/local/gromacs
# Working dir: /home/shizuku/mizusama/DPP4
# Command line:
# gmx rms -s md_0_1.tpr -f md_0_1.xtc -o rmsd.xvg
# gmx rms is part of G R O M A C S:
#
# GRoups of Organic Molecules in ACtion for Science
#
@ title "RMSD"
@ xaxis label "Time (ps)"
@ yaxis label "RMSD (nm)"
@TYPE xy
@ subtitle "Backbone after lsq fit to Backbone"
0.0000000 0.0004978
10.0000000 0.0943481
20.0000000 0.0992674
30.0000000 0.1152738
(省略）
4940.0000000 0.2464143
4950.0000000 0.2461264
4960.0000000 0.2485593
4970.0000000 0.2470721
4980.0000000 0.2526056
4990.0000000 0.2768429
5000.0000000 0.2775707
グラフを書くと以下のような感じでした。1000 psほどで、RMSDが0.02Åになり、タンパク質が安定したことがわかります。また、4500 ps付近でさらにRMSDが変化しているので、さらにsimulaitonの長さを長くして、安定しているかどうか見てもよいとは思います！
10. 安定な構造の抽出
1000 ps以降は安定なタンパク質構造だと仮定して、2000 ps付近のpdbファイルを抽出していきます。
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o frame_2ns.pdb -dump 2000
シミュレーションの座標データから、2ns時点（2000 ps）の構造を抽出し、PDBファイルとして保存する。
s md_0_1.tpr→ 構造情報を含むトップファイル（md_0_1.tpr）を指定。f md_0_1.xtc→ シミュレーションのトラジェクトリ（md_0_1.xtc）を指定。o frame_2ns.pdb→ 出力ファイル（2000 ps時点の構造をframe_2ns.pdbに保存）。dump 2000→ 2000 ps（=2 ns）時点のフレームを抽出。
以下の通りに出てくるので、1を選択してください。
:-) GROMACS - gmx trjconv, 2024 (-:
Executable: /usr/local/gromacs/bin/gmx
Data prefix: /usr/local/gromacs
Working dir: /home/shizuku/mizusama/DPP4
Command line:
gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o frame_2ns.pdb -dump 2000
Note that major changes are planned in future for trjconv, to improve usability and utility.
Will write pdb: Protein data bank file
Reading file md_0_1.tpr, VERSION 2024 (single precision)
Reading file md_0_1.tpr, VERSION 2024 (single precision)
Select group for output
Group 0 ( System) has 327482 elements
Group 1 ( Protein) has 23242 elements
Group 2 ( Protein-H) has 11898 elements
Group 3 ( C-alpha) has 1452 elements
Group 4 ( Backbone) has 4356 elements
Group 5 ( MainChain) has 5810 elements
Group 6 ( MainChain+Cb) has 7182 elements
Group 7 ( MainChain+H) has 7214 elements
Group 8 ( SideChain) has 16028 elements
Group 9 ( SideChain-H) has 6088 elements
Group 10 ( Prot-Masses) has 23242 elements
Group 11 ( non-Protein) has 304240 elements
Group 12 ( Water) has 304212 elements
Group 13 ( SOL) has 304212 elements
Group 14 ( non-Water) has 23270 elements
Group 15 ( Ion) has 28 elements
Group 16 ( Water_and_ions) has 304240 elements
Select a group:
こちらで作成したpdbファイルを基に、続くin silicoスクリーニングをしてください。
最後に
GROMACSを用いたタンパク質の安定化について、環境構築からRMSD解析、安定構造の抽出まで解説しました。AMBER99SB-ILDNフォースフィールドを使用し、より精密な分子動力学シミュレーションを行いましょう！
