ここではGromacsのチュートリアル1：Lysozyme in Water を日本語化している。初学者向けに全体の流れを重視しているため、詳細は英語版のページを参照してほしい。
シミュレーションの前準備
このチュートリアルでは、水とイオンの入ったボックスにタンパク質（lysozyme）が含まれる系を扱う。チュートリアルではGromacs 2018.xシリーズを使用しているため、手元のバージョンを確認しておく。
タンパク質の構造ファイルを入手
まず、タンパク質の構造ファイルをダウンロードする必要がある。RCSBのサイトからlysozyme（PDBコード：1AKI）を検索し、PDBファイルをダウンロードする。このファイルをVMDで確認すると、タンパク質とその周りの結晶水の構造を見ることができる。今回は結晶水の構造は必要ないので、PDBファイルのHOHを含む行をgrepで除いておく。
grep -v HOH 1aki.pdb > 1AKI_clean.pdb
これで新しく作成したファイルをVMDで確認すると、結晶水が取り除かれている。
Gromacs入力ファイルの作成
次に、タンパク質の構造ファイルからMD計算に使用する以下のファイルを作成する。
- トポロジーファイル（topol.top）
- 位置拘束ファイル（posre.itp）
- Gromacsフォーマットの構造ファイル（1AKI_processed.gro）
これらのファイルはgmx pdb2gmxを使用してPDBファイルから一括で作成できる。-waterは計算で使用する水モデルの種類を指定するオプションである。
gmx pdb2gmx -f 1AKI_clean.pdb -o 1AKI_processed.gro -water spce
コマンドを実行すると使用する力場を聞かれるので、OPLS-AAを選択する。15と入力してEnterキーを押す。
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)
これで必要なファイルが作成されているはずだ。
トポロジーファイルの詳細
ここではgmx pdb2gmxで作成したトポロジーファイル（topol.top）について簡単に紹介する。手元の出力ファイルをテキストエディタなどで開いて確認してほしい。
まず、冒頭のいくつかのコメント行の後、力場パラメータに関する記述がある。今回は力場にOPLS-AAを選択したので、これ以降のパラメータはOPLS-AAで導入されるものとなっている。
#include "oplsaa.ff/forcefield.itp"
次に [ moleculetype ] 、[ atoms ] のセクションがあり、タンパク質とその構成原子の情報が書かれている。特に [ atoms ] では各原子に1から番号が割り振られ、それがタンパク質で何番目のどの残基に属する原子であるか、元素の種類、電荷、質量などが記されている。
（下の例では1番の原子は1番目のLys残基に属する窒素原子であり、電荷-0.3、質量14.0067であることがわかる）
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 LYS rtp LYSH q +2.0
1 opls_287 1 LYS N 1 -0.3 14.0067 ; qtot -0.3
2 opls_290 1 LYS H1 1 0.33 1.008 ; qtot 0.03
3 opls_290 1 LYS H2 1 0.33 1.008 ; qtot 0.36
4 opls_290 1 LYS H3 1 0.33 1.008 ; qtot 0.69
5 opls_293B 1 LYS CA 1 0.25 12.011 ; qtot 0.94
6 opls_140 1 LYS HA 1 0.06 1.008 ; qtot 1
また、続く [ bonds ] 、[ pairs ] 、[ angles ] 、[ dihedrals ] のセクションでは、タンパク質の結合に関する情報が書かれている。ここでどの原子に結合があるか、どの原子で結合角や二面角が生じるかを知ることができる。さらに続けて、位置拘束パラメータ（gmx pdb2gmxで作成したposre.itp）を記述している。ここまででタンパク質のトポロジーの定義が完了する。
最後に、ファイルの残りの部分でタンパク質以外の分子についての定義が行われている。具体的には溶媒（水分子）とイオンである。今回は-waterで水分子のモデルをSPE/Cとしたが、他にもSPC、TIP3P、TIP4Pなどが選択できる。
; Include water topology
#include "oplsaa.ff/spce.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include generic topology for ions
#include "oplsaa.ff/ions.itp"
以上がトポロジーファイルの内容である。ファイル末尾の [ system ] には系の名前、[ molecules ] には系に含まれる分子とその個数が記述されているので、作成した構造と一致するかを確認しておく。
[ system ]
; Name
LYSOZYME
[ molecules ]
; Compound #mols
Protein_A 1
系の構築
Gromacsのトポロジーファイルの内容を理解したところで、系の構築に戻りたいと思う。gmx pdb2gmxで作成したGromacsフォーマットの構造ファイル（.groファイル）には、まだタンパク質の構造しか含まれていない。そこで次に、タンパク質が入るようなシミュレーションボックスを作成し、その中を溶媒で満たす操作を行う。また、系内のタンパク質の電荷を打ち消すためにイオンを追加する。
シミュレーションボックスの作成
シミュレーションボックスの作成はgmx editconfで行う。今回、シミュレーションボックスには周期的境界条件を満たす立方体を用いる。-cでタンパク質をボックス内の中央に配置し、-dでタンパク質とボックス境界の距離を1.0 nm以上空けるようにしている。これにより、周期的境界条件では隣のタンパク質との距離を2.0 nm以上確保することができる。また、-btはボックスタイプを立方体に指定するオプションである。
gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic
溶媒の追加
シミュレーションボックスが定義できたので、続いてgmx solvateにより系を溶媒で満たす。-cpはタンパク質を含む系の構造、-csは溶媒（ここではSPC水モデル）の構造を指定している。SPCの構造ファイル（spc216.gro）はGromacsに標準でインストールされているので、特別にファイルを用意する必要はない。
gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top
これでボックス内を溶媒で満たすことができた。-pで指定したトポロジーファイルを確認すると [ molecules ] に溶媒が追加されていることがわかるだろう。
[ molecules ]
; Compound #mols
Protein_A 1
SOL 10832
イオンの追加
さらに系に対してイオンを追加する。トポロジーファイルを確認するとわかるが、タンパク質の電荷は+8（qtot 8）であるため、イオンを加えて電荷を打ち消す必要がある。イオンの追加にはgmx genionを用いるが、その入力ファイルとして全原子のパラメータを含むファイル形式の.tprファイルが必要となる。.tprファイルの作成は、系の構造ファイルとトポロジーファイルの両方を使用してgmx gromppにより行う。
gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr
この時、gmx gromppはさらに追加で入力ファイルを要求する。-fで指定した.mdpファイル（molecular dynamics parameter file）はMD計算の条件・実行方法などを設定するファイルであり、gmx gromppでバイナリ形式のMD計算実行ファイル（.tprファイル）を作成するために必要となる。.mdpファイルは通常、温度・圧力・シミュレーション時間などの設定を含むが、現時点ではイオンを追加する設定のみが書かれている。ions.mdpファイルをここからコピーして、gmx gromppを実行してほしい。
gmx gromppでions.tprが作成できたら、いよいよイオンを追加する。次のオプションを指定してgmx genionを実行する。
gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral
実行すると溶媒として利用する分子を聞かれるので、13（SOL）を選択する。これにより、溶媒の水分子の一部をイオンに置き換えることができる。-pnameおよび-nnameで追加する陽イオンおよび陰イオンの種類を指定しており、-neutralで系全体の電荷が相殺されるようにイオンを加えている。今回は+8の電荷を打ち消すために8つのClイオンが追加されている。-pで指定したトポロジーファイルを開き、[ molecules ] にイオンが反映されていることを確認しておこう。
[ molecules ]
; Compound #mols
Protein_A 1
SOL 10636
CL 8
以上で系の構築は完了である。VMDで確認してみると、次のような構造が見られるだろう。
まとめ①
今回の内容をまとめる。系の構築のために以下の手順を踏んだ。
- タンパク質のPDBファイルを
gmx pdb2gmxによりGromacsフォーマットに変換した。 - トポロジーファイルの内容を理解した。
- シミュレーションボックスを
gmx editconfで作成し、gmx solvateによりボックス内を溶媒で満たした。 - ボックス内に
gmx genionでイオンを追加した。
次回は作成した系を用いて実際にMD計算をしてみたいと思う。