MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜

MDシミュレーションのチュートリアル記事を書いてみました。
今回題材として選んだタンパク質はPDB: 1LKEです。

![1lke1.png](https://qiita-image-store.s3.amazonaws.com/0/224327/80cde3f7-1103-e5d8-48dc-eb7c511934d9.png)

この構造のポイントは以下の3つです。

  - リガンドとしてDigoxigeninが結合している
  - TYR-117とGLY-123の間に構造データが存在しないミッシングループが存在し、SER-166の後ろにもいくつか残基が続く（参考：PDB: 1LNM）
  - Cys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170でもジスルフィド結合が形成されている（参考：PDB: 1LNM）

これらにどう対処すればよいかを知っていればほとんどのタンパク質MDシミュレーションに対応できるからです。
（2018年9月13日）　GitHub（ https://github.com/YoshitakaMo/MDsimulations ）に使用した全ファイルと完成済みファイル例を置きました（minimizationまで）。参考にしてください。
では行ってみましょう。

## 環境

  - macOSまたはLinux環境推奨。Windows 10の場合は[Windows Subsystem for Linux（通称WSL）をインストール](https://qiita.com/Aruneko/items/c79810b0b015bebf30bb)すれば、Ubuntuのような環境を構築することができる。**GPUがあると計算が最大10倍速くらいになる**。
  - [AmberTools 22（AMBER20があればそれも）がインストールされている](https://qiita.com/Ag_smith/items/d2e86dda17f190a3dd0f)
  - [Gromacs 2023.2がインストールされている](https://qiita.com/Ag_smith/items/2432375399e7d7868138)
  - Modellerがインストールされている（ミッシング残基が存在する場合）

環境設定についてはそれぞれのリンクを読んで下さい。Modellerについてはすぐ下で述べます。

## Modellerのインストール(macOSの場合)

macOSの場合は[ターミナルからHomebrewをインストール](http://brew.sh/index_ja.html )してある状態で、以下のコマンドを打ちます。

```bash
brew tap salilab/salilab
brew install modeller
```

これだけでインストールできます。

ただし、使用開始する前にファイル `/usr/local/opt/modeller/modlib/modeller/config.py` にLicense Keyを入れておく必要がある(XXはバージョン。現在は10.1)。vimやVS Codeなどでこのファイルを開いて中身を以下のように書き換える。

```
install_dir = r'/usr/local/Cellar/modeller/10.1'
license = r'MODELIRANJE'
```

これで完了です。

Windows 10、WSLまたはUbuntuを使っている方は、ターミナルを開いてModeller公式ページの[Ubuntu用のインストール方法](https://salilab.org/modeller/10.1/release.html#deb)を試してみてください。

## PDBファイルのダウンロード

適当な作業ディレクトリを作成します。ここでは`~/Desktop/1LKE`とします。まずはそのディレクトリを作り、そこにRCSB PDBのウェブサイトからPDBファイルのPDB: 1LKEをダウンロードしてきます。

```bash
# ディレクトリの作成
mkdir -p ~/Desktop/1LKE/PDB
# ディレクトリへの移動
cd ~/Desktop/1LKE/PDB
# PDB:1LKEのファイルをcurlコマンドでダウンロード
curl -o 1LKE.pdb "https://files.rcsb.org/download/1LKE.pdb"
```

このPyMOLなどの可視化ソフトで構造を見ておきましょう。この構造で確認しておくポイントは以下の3つです。

  - リガンドとしてDigoxigeninが結合している
  - TYR-117とGLY-123の間に構造データが存在しないミッシングループが存在し、SER-166の後ろにもいくつか残基が続く（参考：[PDB: 1LNM](http://www.rcsb.org/structure/1LNM)）
  - Cys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170でもジスルフィド結合が形成されている（参考：[PDB: 1LNM](http://www.rcsb.org/structure/1LNM)）

このMDシミュレーションのチュートリアルではこれらの部分の設定方法も含めながら進んでいきます。

## AmberToolsを使った1LKEのトポロジーファイル作成

### トポロジーファイルとは

MDシミュレーションの基本的な考えは、タンパク質などの生体分子を構成する各原子が剛体球とばねで繋がれているとするモデルを大前提としています。これに加え、3原子間の角度についても同様な球とばねモデルを適用し（例として水分子のH-O-H角度が約104.5°となるように制御するなど）、さらに4原子間の二面角についても設定し、分子の平面性を制御したりしています。また、例えば同じ炭素原子であってもsp3炭素とsp2炭素の振る舞いが違うように、各原子について個別に適切なパラメータを設定する必要があります。

幸いなことに、タンパク質のシミュレーションではタンパク質を構成するアミノ酸はおよそ20種類に限られ、かつパラメータもアミノ酸ごとに同じものを使いまわしても大きな問題にはならないことが判明していますので、これらのパラメータセットは**力場（forcefield）**としてライブラリ化されており、いくつかの信頼と実績のある力場が世界中で広く使われています。この力場を使ってタンパク質のMDシミュレーションの系に対して作成したファイルのことをトポロジーファイルと言います。

トポロジーファイルはタンパク質のMDシミュレーションをする上で必須のファイルであり、パラメータは動きと精度に大きく関わるので非常に重要なものとなります。

### Step 1-1. PDBファイルの原子の座標部分だけ取り出す

AMBERとAmberToolsに付属の**tleap**というモジュールを使うと、PDBファイルからトポロジーファイルを作成することができます。さっそくこれを使ってみましょう。

まず`cp 1LKE.pdb 1LKE_truncated.pdb`で先程ダウンロードしたファイルのバックアップを確保しておきます。続いてこの`1LKE_truncated.pdb`ファイルをテキストエディタで見てみます。すると、462行目から1793行目まで、ATOM、TER、HETATMで始まる文字列があるのがわかります。この部分だけを確保し、残りは消してしまいましょう。

    ATOM      1  N   ASP A   5       2.533 -40.650  52.271  1.00 69.29           N
    ATOM      2  CA  ASP A   5       1.130 -40.378  52.680  1.00 70.86           C
    ATOM      3  C   ASP A   5       1.018 -39.026  53.360  1.00 66.11           C
    ATOM      4  O   ASP A   5       1.950 -38.566  54.021  1.00 68.69           O
    ATOM      5  CB  ASP A   5       0.598 -41.480  53.609  1.00 74.53           C
    （中略）
    HETATM 1330  O   HOH A 556      28.620 -15.856  57.423  1.00 53.96           O
    HETATM 1331  O   HOH A 557       4.205 -24.295  47.080  1.00 39.36           O
    HETATM 1332  O   HOH A 558      27.529 -13.304  59.699  1.00 58.19           O
    END

として切り出しておきます。これを`1LKE_truncated.pdb`として保存します。

### Step 1-2. Modellerを使ったミッシング残基の補完（存在する場合のみ）

このステップの動作は、上に引き続き`~/Desktop/1LKE/PDB`を想定しています。

今回のタンパク質1LKEには一部の残基が見えていないミッシング領域が存在します。
RCSB PDBのウェブサイトでタンパク質構造1LKEの[sequenceタブ](https://www.rcsb.org/sequence/1LKE)を見ると、TYR-117とGLY-123の間にはDEDKKのアミノ酸が存在することがわかります。

![28391.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/881e1a68-514a-db94-4ad8-9a83bf2e7eca.png)

この画像で**UNMODELED**の行に注目します。この行は座標情報が存在していないアミノ酸残基の場所を示しています。UNMODELEDの行の118-122の灰色の四角にマウスカーソルを載せると情報が表示されます。

![734719.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/f41554de-ca65-fbb0-b514-c84f7fe15688.png)

118番から122番までのアミノ酸残基、DEDKKがミッシング（UNMODELED）であることがわかります。この箇所は結晶中でさえもゆらぎが大きいためか、この相当する部位の電子雲がほとんど見えず、その座標部分がPDBファイル内に記載されていません。しかし、このDEDKK配列は必ず存在するはずであり、またMDシミュレーションをする場合には当然ながらこれらを含んだ構造にしている必要があります。

またジスルフィド結合（**DISULFIDE BRIDGE**）についても注目します。上と同じSequenceタブでマウスをDISULFIDE BRIDGEの部分に載せると、CYS-8とCYS-115の間でジスルフィド結合が形成されていることがわかります。

![1922.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/3e110a83-89a4-8ba0-9d8f-96b55c403979.png)

**ジスルフィド結合はタンパク質構造を保つ上でとても重要な働きをするので、必ず入れる必要があります**。

ここで、やや発展的な内容ですが、ジスルフィド結合が他にも存在していないか調べてみましょう。RCSB PDBのStructureタブに戻って下の方にはLiteratureという情報欄があります。そこには、PDB: 1LKEに関係するタンパク質がいくつか示されています。

![12821.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/1b9d680d-2d9d-8097-1892-832b7696be0a.png)

ここにPDB: 1KXO, 1LNMのことが書かれています。こちらのいずれかのPDBデータについてSequenceタブを調べてみると、そちらの構造にはCYS-42, CYS-170間にもう1つのジスルフィド結合が存在することがわかります。（参考：http://www.rcsb.org/sequence/1LNM）

![94829.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/2f6b7ecf-6856-5651-07f5-a7a771a84b8e.png)

**アミノ酸配列が同じであればタンパク質構造の概形は同じ**なので、こちらでもう1つのジスルフィド結合が見えたのであれば、1LKEでもこのジスルフィド結合が存在している可能性は極めて高いと言えます。そしてジスルフィド結合はタンパク質構造を保つ上でとても重要なので（2回め）、この情報を入れた上でMDシミュレーションをしたいということになります。ただし、このような座標情報がないアミノ酸残基の座標情報を人間が手で適切に補完してあげるのは困難です。そこで機械的に行いたいわけですが、こうしたミッシング領域の座標の補完を行うとき、よく用いられるのが[Modeller](https://salilab.org/modeller/)というソフトです。ModellerはSali Labが開発しているPythonベースのタンパク質構造編集ソフトウェアです。使いこなせれば非常に強力にタンパク質構造を編集できます。

今回はModellerで以下の3つの設定を踏まえたモデリングを行います。

  - TYR-117とGLY-123の間のミッシング領域についてDEDKKのアミノ酸ループ構造を補完する。
  - さらに、SER-166の後ろにもいくつか残基が続くので、EAACKVNのタンパク質構造を補完する（今回はとりあえずCYS-170で形成されるジスルフィド結合周辺の構造を補完するのが狙いです）。
  - PDBファイル上のCys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170（EAACKVN配列中のC）でもジスルフィド結合が形成されている。

Modellerの実行には実行バイナリ`mod10.1`（10.1はバージョン番号）にPATHが通っていることの他に（`which mod10.1`とコマンドを入力してみてPATHが表示されれば成功）、Modellerに読み込ませるPythonスクリプトの`align.py`と配列アラインメントファイルの`alignment.ali`の3つが必要です。

では、まずは以下のような`align.py`ファイルを用意します。

```python:align.py
#!/usr/bin/env python3

import modeller
import modeller.environ
import modeller.util.logger
import modeller.automodel
import modeller.parallel

# 並列計算用のコア数を設定(2021年3月17日追加)
# https://salilab.org/modeller/9.15/manual/node457.html
ncpus = 8

j = modeller.parallel.job(host="localhost")
for i in range(ncpus):
    j.append(modeller.parallel.LocalWorker())

modeller.util.logger.log.verbose()
env = modeller.environ()
# HETATMの座標を保存するようにフラグを追加
env.io.hetatm = True
# 水分子の座標を保存するようにフラグを追加
env.io.water = True

# directories for input atom files
env.io.atom_files_directory = ['.', '1LKE_truncated.pdb']

a = modeller.automodel.AutoModel(env, alnfile = 'alignment.ali', knowns = '1LKE', sequence = '1LKE_fill')
a.starting_model= 1
a.ending_model  = 8
a.md_level       = modeller.automodel.refine.fast

# 並列計算の実行（2021年3月17日追加）
a.use_parallel_job(j)
a.make()
```

現行のmodellerでは、マルチコアCPUを利用して並列に構造モデリングを行う機能が実行されています。これを利用することでたくさんのモデリングを行うときに大幅に計算時間を短縮することができます。

残基番号が4つずれていることに注意してください。これはModellerでは以下の`alignment.ali`で定義した構造情報を持つアミノ酸残基を必ず1番めに再設定するという性質があるためです（言い換えれば、アウトプットされるPDBの残基番号を基準としてここの数字を設定します）。PDB: 1LKEの構造情報はChain Aの残基番号5から始まっています。
ミッシング領域に含まれるCys-170（完成後のモデル構造上では166）についてもSS結合を定義しています。

次にアラインメントファイルの`alignment.ali`です。

```shell:alignment.ali
>P1;1LKE
structureX:1LKE_truncated:5:A::A:undefined:undefined:-1.00:-1.00
DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWAEYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSK
IGKIYHSYTIGGVTQEGVFNVLSTDNKNYIIGYFCSY-----GHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKL
VYSDFS-------...........................................................*
>P1;1LKE_fill
sequence:::::::::
DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWAEYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSK
IGKIYHSYTIGGVTQEGVFNVLSTDNKNYIIGYFCSYDEDKKGHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKL
VYSDFSEAACKVN...........................................................*
```

上の方が構造情報を持つテンプレート配列（PDBの構造と同じに揃える）、下のほうが補完した配列です。必ず配列を対応させておく必要があり、ミッシングに相当する配列は`-`記号で埋めます。また、`align.py`で`env.io.hetatm = True`、`env.io.water = True`を定義していた場合、この配列情報中の`.`記号でそれぞれの残基を代用することができます（この後の計算中のModeller中間ファイル`*.ini`では`BLK`として処理しているようです。）。上の配列ファイルでの`.`の数は、Digoxigeninが1つ分、水分子が58つで計59個あります。

これの2つのファイルを用意したら、`cd ~/Desktop/1LKE/PDB`でそのディレクトリを移動し、`1LKE_truncated.pdb`と`align.py`, `alignment.ali`の3つのファイルが同じディレクトリに存在することを確認した上で、ターミナル上で

```shell
$ python3 align.py > align.log
```

を実行します。

並列計算の数値によりますが、およそ30秒くらいで計算は終了し、8つのモデリング構造がアウトプットされます。

`align.log`の最後の方に書かれている構造スコアを見てみます。この値が低いほどよい構造とされます。

```shell
>> Summary of successfully produced models:
Filename                          molpdf
----------------------------------------
1LKE_fill.B99990001.pdb       1138.31165
1LKE_fill.B99990002.pdb       1333.03198
1LKE_fill.B99990003.pdb       1128.88806
1LKE_fill.B99990004.pdb       1168.13550
1LKE_fill.B99990005.pdb       1195.39551
1LKE_fill.B99990006.pdb       1367.64893
1LKE_fill.B99990007.pdb       1110.27405
1LKE_fill.B99990008.pdb       1391.83472
```

モデリングにはある程度の乱数を利用しているため、この中で最も良かった構造の番号は毎回異なります。

値が最もよかったものを選んでPyMOLなどで構造を見てみましょう。またここで、構造を見てみたときに注意してほしい点があります。

![1lke3.png](https://qiita-image-store.s3.amazonaws.com/0/224327/4701693d-9a44-093a-6a62-4fb8d58120ac.png)

このようにモデリングしたループがSS結合の中に入ってしまっている構造がいくつかあります。この構造は明らかに良くないので選ばないようにしましょう。

また、水分子とDigoxigeninがモデル後の構造にも保存されているでしょうか？大丈夫であれば、ここでは最もスコアの小さかった`1LKE_fill.B99990007.pdb`を選んで話を続けます。

Modellerで作られたこのファイルを開いてみると、残基番号が1番からに振り直されています。ここで（必要というわけではないですが、話をわかりやすくするため）、Gromacsの`gmx editconf`コマンドを使って残基番号を改めて5番からに振り直します。

```shell
# gmx editconfコマンドで、残基番号を5から振り直したものを1LKEmodel.pdbに出力する
gmx editconf -f 1LKE_fill.B99990007.pdb -o 1LKEmodel.pdb -resnr 5
```

入力`-f`を先程のモデル構造ファイルに、`-o`に出力ファイル名を指定し、renumberingを5番からに指定します。こうして得られた`1LKEmodel.pdb`ファイルを使って次のStepに進みます。

### Step 1-3. His残基の水素の位置の設定

ヒスチジン（His）残基は[図のようなpKa曲線](http://bio.wikia.com/wiki/Histidine)を持っており、pH = 6.0以下で側鎖のイミダゾール基に水素を2つもつ状態になり、それ以上になると水素を1つ、δ位かε位のどちらか1つについた状態になります。MDシミュレーションは通常pH 6.0〜9.0程度の溶液を想定して行うことが多いため、水素原子が1つついた状態のヒスチジンを取り扱うことがほとんどになります。

さて、イミダゾール基のδ位かε位のどちらか1つに水素原子がつくわけですが、水素原子はとても小さいのでタンパク質結晶構造解析であってもなかなか水素原子の位置を捉えることはできません。しかしどちらに設定するべきかというのは、場合によっては重要な要素になります。そこで通常は、**タンパク質結晶構造中で、水素結合が形成されていそうな原子の方に水素原子を付加し、水素結合が形成されてそうでなければ適当にどちらか一方（通常ε位）に水素原子を設定する**という手法を取ります。

今回使用するAMBERのtleapモジュールの場合、読み込むPDBファイルのHIS残基名を`HID`に変更するとδ-protonatedのヒスチジン、`HIE`にするとε-protonatedのヒスチジンとして認識してくれるようになります。

では`1LKE_model.pdb`についてやってみます。このタンパク質内には5つのHis残基があります。このうち注目するべきはHis-86です。

![1lke2.png](https://qiita-image-store.s3.amazonaws.com/0/224327/fabfeaa9-3f4c-3ec1-d2be-94897d3bdfc1.png)

His-86をみると、δ位に水素原子があったほうがdigoxigeninと水素結合を形成していい感じになりそうなのがわかります。His-35はε位に水素があったほうが水素結合を形成しやすそうなのがわかります（δ位ってのはイミダゾールの五角形の根元から近い方の窒素原子がδ、遠い方がεです）。他のヒスチジン残基はどっちでもよさそうです。なので、His-86だけは残基名を`HID`にしてδ-protonatedということにしておきましょう。

先程の`1LKEmodel.pdb`ファイルのHIS 86あたりの残基名をHIDにします。

```
ATOM    654  N   HID    86       4.953 -19.957  36.923  1.00 79.07           N
ATOM    655  CA  HID    86       5.853 -20.098  35.812  1.00 79.07           C
ATOM    656  ND1 HID    86       9.513 -20.269  37.045  1.00 79.07           N
ATOM    657  CG  HID    86       8.205 -19.890  36.843  1.00 79.07           C
ATOM    658  CB  HID    86       7.209 -20.754  36.129  1.00 79.07           C
ATOM    659  NE2 HID    86       9.292 -18.249  37.946  1.00 79.07           N
ATOM    660  CD2 HID    86       8.086 -18.655  37.401  1.00 79.07           C
ATOM    661  CE1 HID    86      10.118 -110.12  37.707  1.00 79.07           C
ATOM    662  C   HID    86       5.229 -21.009  34.812  1.00 79.07           C
ATOM    663  O   HID    86       4.616 -22.013  35.174  1.00 79.07           O
```

ここのHIS残基をすべてHIDに変更して保存したファイルを`1LKE_trun_hid.pdb`としておきましょう。あとは自動的にε-protonatedになります。

### Step 1-4. SS結合の設定（存在する場合のみ）

今回のケースのように、SS結合が存在する場合はもういくつか設定する必要があります。AmberToolsの場合、SS結合を形成させたい場合は以下の2つを行う必要があります。

  - SS結合を形成する2つのCYS残基の名前を**CYX**に変更する
  - SS結合を形成させる2残基のCYXのSGを結合させる設定を書く

1つめは上のHIS→HIDと同じように名前をCYS→CYXに変えるだけです。

```
ATOM     15  C   ALA A   7      -0.641 -33.155  52.030  1.00 49.97           C
ATOM     16  O   ALA A   7       0.142 -33.209  52.966  1.00 42.23           O
ATOM     17  CB  ALA A   7      -3.095 -33.395  52.296  1.00 58.78           C
ATOM     18  N   CYX A   8      -0.532 -32.270  51.035  1.00 46.12           N
ATOM     19  CA  CYX A   8       0.517 -31.252  51.058  1.00 43.44           C
ATOM     20  C   CYX A   8       0.439 -30.449  52.312  1.00 46.59           C
ATOM     21  O   CYX A   8      -0.638 -29.981  52.691  1.00 51.88           O
ATOM     22  CB  CYX A   8       0.410 -30.329  49.847  1.00 38.64           C
ATOM     23  SG  CYX A   8       0.925 -31.183  48.413  1.00 36.24           S
ATOM     24  N   PRO A   9       1.577 -30.296  52.978  1.00 45.14           N
ATOM     25  CA  PRO A   9       1.646 -29.492  54.192  1.00 47.94           C
ATOM     26  C   PRO A   9       1.342 -28.044  53.866  1.00 53.73           C
（中略）
ATOM    880  CE1 PHE A 114       3.892 -29.444  41.508  1.00 29.36           C
ATOM    881  CE2 PHE A 114       3.175 -27.133  41.524  1.00 28.15           C
ATOM    882  CZ  PHE A 114       2.961 -28.488  41.165  1.00 23.55           C
ATOM    883  N   CYX A 115       4.434 -28.614  45.678  1.00 20.99           N
ATOM    884  CA  CYX A 115       3.916 -29.898  46.202  1.00 20.88           C
ATOM    885  C   CYX A 115       2.614 -30.206  45.489  1.00 31.70           C
ATOM    886  O   CYX A 115       1.789 -29.310  45.288  1.00 24.52           O
ATOM    887  CB  CYX A 115       3.632 -29.806  47.726  1.00 26.63           C
ATOM    888  SG  CYX A 115       2.960 -31.364  48.464  1.00 28.36           S
ATOM    889  N   SER A 116       2.467 -31.439  45.049  1.00 32.36           N
ATOM    890  CA  SER A 116       1.156 -31.884  44.551  1.00 48.05           C
```

今回、このファイルにはCYSが8, 42, 115, 170にあるのですが、いずれもCYXに変更します。**115, 170のCYSもCYXにします（大切なことなので2回言いました）**。 そしてこのファイルを`1LKE_trun_hid_cyx.pdb`として保存しておきます。

2つめのSS結合の設定ですが、これについては以下のStep 2. tleapで紹介します。

### Step 1-5. リガンドのパラメータファイルの用意（存在する場合のみ）

タンパク質は化合物と結合して機能を発揮するタイプのものが非常にたくさんあります。1LKEにもリガンド（Digoxigenin）が結合しています。しかし、20種類の"カノニカル"な（正統な）アミノ酸と違って、一般の化合物の形状は非常に多くの形状を持つので、その都度それについて適切な力場や電荷を用意してあげなければなりません。

この化合物のパラメータファイルの用意自体は非常に長くなってしまうため、作成方法は別の記事を参照してください。ここではすでに完成されたファイル`dog.mol2`が与えられているものとして次へ進みます。

```shell:dog.mol2
@<TRIPOS>MOLECULE
DOG
   62    66     1     0     0
SMALL
resp
@<TRIPOS>ATOM
      1 O23          7.2640     0.6350     0.3660 o          1 DOG     -0.603790
      2 C23          6.1540     0.2200     0.1420 c          1 DOG      0.913720
      3 C22          4.9650     0.9130    -0.3940 ce         1 DOG     -0.552430
      4 O21          5.7960    -1.0930     0.3690 os         1 DOG     -0.543340
      5 C21          4.4240    -1.2860     0.0110 c3         1 DOG      0.253200
      6 C20          3.9340     0.0590    -0.4810 c2         1 DOG      0.088750
      7 C17          2.5540     0.3660    -0.9980 c3         1 DOG     -0.002670
      8 C16          2.0410    -0.6750    -2.0480 c3         1 DOG     -0.263190
      9 C13          1.4010     0.4740     0.0890 c3         1 DOG      0.635860
     10 C18          1.9120     0.6990     1.5210 c3         1 DOG     -0.416180
     11 C12          0.4900     1.6560    -0.3400 c3         1 DOG      0.078140
     12 O12          1.2290     2.8590    -0.1190 oh         1 DOG     -0.703970
     13 C14          0.5950    -0.8510    -0.0860 c3         1 DOG     -0.051480
     14 C15          0.6220    -1.0940    -1.6100 c3         1 DOG      0.108190
     15 O14          1.4070    -1.8660     0.5410 oh         1 DOG     -0.712600
     16 C8          -0.8080    -0.8380     0.5810 c3         1 DOG      0.042480
     17 C7          -1.5700    -2.1580     0.3480 c3         1 DOG     -0.040690
     18 C9          -1.6660     0.3920     0.1730 c3         1 DOG      0.040630
     19 C11         -0.8550     1.6840     0.3920 c3         1 DOG     -0.475450
     20 C10         -3.0760     0.4040     0.8680 c3         1 DOG      0.163710
     21 C19         -2.9500     0.6320     2.3920 c3         1 DOG     -0.573120
     22 C5          -3.8090    -0.9540     0.6040 c3         1 DOG      0.141680
     23 C6          -2.9460    -2.1570     1.0260 c3         1 DOG     -0.331830
     24 C4          -4.2990    -1.1120    -0.8530 c3         1 DOG      0.042980
     25 C3          -5.1640     0.0520    -1.3340 c3         1 DOG      0.272340
     26 O32         -6.3830     0.0000    -0.5830 oh         1 DOG     -0.710140
     27 C2          -4.4320     1.3850    -1.1280 c3         1 DOG     -0.171740
     28 C1          -3.9800     1.5490     0.3280 c3         1 DOG     -0.056550
     29 H1           4.9820     1.9620    -0.6620 ha         1 DOG      0.218680
     30 H2           3.8680    -1.6470     0.8800 h1         1 DOG      0.060400
     31 H3           4.3750    -2.0640    -0.7600 h1         1 DOG      0.060400
     32 H4           2.6320     1.3530    -1.4600 hc         1 DOG      0.057230
     33 H5           2.6940    -1.5500    -2.0730 hc         1 DOG      0.070310
     34 H6           2.0490    -0.2540    -3.0590 hc         1 DOG      0.070310
     35 H7           1.0840     0.8150     2.2280 hc         1 DOG      0.103440
     36 H8           2.5110    -0.1450     1.8640 hc         1 DOG      0.103440
     37 H9           2.5160     1.6060     1.5680 hc         1 DOG      0.103440
     38 H10          0.2860     1.5600    -1.4200 h1         1 DOG      0.036850
     39 H11          0.6880     3.5980    -0.4370 ho         1 DOG      0.444350
     40 H12          0.4090    -2.1410    -1.8550 hc         1 DOG     -0.017340
     41 H13         -0.1430    -0.4960    -2.1160 hc         1 DOG     -0.017340
     42 H14          1.0310    -2.7280     0.3040 ho         1 DOG      0.443700
     43 H15         -0.6050    -0.7670     1.6580 hc         1 DOG      0.058100
     44 H16         -0.9860    -3.0040     0.7360 hc         1 DOG      0.031780
     45 H17         -1.6900    -2.3330    -0.7290 hc         1 DOG      0.031780
     46 H18         -1.8490     0.3110    -0.9080 hc         1 DOG      0.078700
     47 H19         -1.4190     2.5530     0.0280 hc         1 DOG      0.151810
     48 H20         -0.6740     1.8570     1.4580 hc         1 DOG      0.151810
     49 H21         -3.9280     0.5080     2.8710 hc         1 DOG      0.137930
     50 H22         -2.2570    -0.0600     2.8790 hc         1 DOG      0.137930
     51 H23         -2.6080     1.6490     2.6160 hc         1 DOG      0.137930
     52 H24         -4.7110    -0.9490     1.2300 hc         1 DOG      0.041280
     53 H25         -2.8100    -2.1520     2.1150 hc         1 DOG      0.095400
     54 H26         -3.4760    -3.0890     0.7900 hc         1 DOG      0.095400
     55 H27         -4.8810    -2.0380    -0.9330 hc         1 DOG     -0.016620
     56 H28         -3.4560    -1.2110    -1.5480 hc         1 DOG     -0.016620
     57 H29         -5.3880    -0.0840    -2.4050 h1         1 DOG      0.002960
     58 H30         -6.9330     0.7500    -0.8610 ho         1 DOG      0.433910
     59 H31         -5.0920     2.2190    -1.4090 hc         1 DOG      0.028770
     60 H32         -3.5790     1.4310    -1.8190 hc         1 DOG      0.028770
     61 H33         -3.4880     2.5190     0.4630 hc         1 DOG      0.039300
     62 H34         -4.8800     1.5760     0.9550 hc         1 DOG      0.039300
@<TRIPOS>BOND
     1     1     2 2
     2     2     3 1
     3     2     4 1
     4     3     6 2
     5     3    29 1
     6     4     5 1
     7     5     6 1
     8     5    30 1
     9     5    31 1
    10     6     7 1
    11     7     8 1
    12     7     9 1
    13     7    32 1
    14     8    14 1
    15     8    33 1
    16     8    34 1
    17     9    10 1
    18     9    11 1
    19     9    13 1
    20    10    35 1
    21    10    36 1
    22    10    37 1
    23    11    12 1
    24    11    19 1
    25    11    38 1
    26    12    39 1
    27    13    14 1
    28    13    15 1
    29    13    16 1
    30    14    40 1
    31    14    41 1
    32    15    42 1
    33    16    17 1
    34    16    18 1
    35    16    43 1
    36    17    23 1
    37    17    44 1
    38    17    45 1
    39    18    19 1
    40    18    20 1
    41    18    46 1
    42    19    47 1
    43    19    48 1
    44    20    21 1
    45    20    22 1
    46    20    28 1
    47    21    49 1
    48    21    50 1
    49    21    51 1
    50    22    23 1
    51    22    24 1
    52    22    52 1
    53    23    53 1
    54    23    54 1
    55    24    25 1
    56    24    55 1
    57    24    56 1
    58    25    26 1
    59    25    27 1
    60    25    57 1
    61    26    58 1
    62    27    28 1
    63    27    59 1
    64    27    60 1
    65    28    61 1
    66    28    62 1
@<TRIPOS>SUBSTRUCTURE
     1 DOG         1 TEMP              0 ****  ****    0 ROOT
```

上では表現できませんが、**最後に空行をはさんでおきます。**これを`dog.mol2`として`PDB`ディレクトリに入れて保存します。

参考： [antechamberを使ったMD用小分子化合物の電荷パラメータの作り方](https://qiita.com/Ag_smith/items/430e9efb32a855d4c511)

### Step 2. tleapへのインプットファイル

続いては`tleap`を用いたトポロジーファイル作成をするために、topディレクトリを作成します。

```bash
# top（トポロジー）のディレクトリ作成と移動
mkdir -p ~/Desktop/1LKE/top
cd ~/Desktop/1LKE/top
```

ここに以下の内容で`leap.in`ファイルを作成します。Step 1-4で紹介したジスルフィド結合の設定はこのファイルの中で行います（`bond mol~ mol~`の部分に対応します）。

(2023年2月8日、小分子向けの汎用力場をGAFF2に変更)

```:leap.in
#AMBER の力場パラメータff14SBを読み込む
source leaprc.protein.ff14SB #for AmberTools22
source leaprc.gaff2 #小分子向けの汎用力場「GAFF2」を読み込む
source leaprc.water.tip3p #for AmberTools22
#追加のイオンの力場の導入
loadAmberParams frcmod.ionsjc_tip3p
#.mol2ファイルの読み込み
DOG = loadMol2 ../PDB/dog.mol2
#pdbを"mol"として読み込む
mol = loadPDB ../PDB/1LKE_trun_hid_cyx.pdb
#SS-bondを設定
bond mol.8.SG mol.115.SG
bond mol.42.SG mol.170.SG
#イオンの追加。末尾の数だけ対応するイオンを配置する。0を指定すると系全体の電荷を0にするようにイオンを入れるいう設定になる。
addIons2 mol Na+ 20
addIons2 mol Cl- 0
#直方体型の周期境界の溶媒和ボックスを形成する。入力PDBファイルのサイズに対し、最低10A余裕を持って形成する。
solvateBox mol TIP3PBOX 10.0
#最後に、"mol"という溶媒和ボックスの系の電荷情報を表示する。0.00000になっていることが理想。
charge mol
#溶媒和された系のトポロジー・初期座標をleap.parm7, leap.rst7としてそれぞれ保存
saveAmberParm mol leap.parm7 leap.rst7
#溶媒和された系のPDBファイルをleap.pdbとして保存
savePDB mol leap.pdb
quit
```

このファイルを作成したら、このディレクトリで`tleap -f leap.in`のコマンドを実行します。しかし、次のようなメッセージが出るはずです。

    Total unperturbed charge:  -0.000000
    Total perturbed charge:    -0.000000
    Checking Unit.
    FATAL:  Atom .R<ASN 173>.A<OXT 15> does not have a type.
    Failed to generate parameters
    Parameter file was not saved.
    Writing pdb file: leap.pdb
       printing CRYST1 record to PDB file with box info
     Converting N-terminal residue name to PDB format: NASP -> ASP
            Quit

`Parameter file was not saved.`と表示されるとトポロジーファイルの作成に失敗しています。まずはここで`Total perturbed charge: -0.000000`となっていることを確認することと、そして`FATAL`が**表示されていない**ことを確認します（上の例では表示されています）。この例では読み込んだPDBファイルのASN 173に謎のOXT15という原子が存在することが原因で失敗しているようです。`XXX does not have a type.`ってのは、力場に対応する原子タイプが存在しないために発生するエラーメッセージです。

読み込んでいるPDBファイルを編集し、問題が起きている箇所の原子を削っておきましょう。削った箇所には区切れを意味する`TER`を書いておくと良いです。

    ATOM   1335  ND2 ASN   173      30.509 -41.532  33.732  1.00 51.92           N
    ATOM   1336  C   ASN   173      27.959 -39.594  35.439  1.00 51.92           C
    ATOM   1337  O   ASN   173      27.577 -39.210  36.576  1.00 51.92           O
    TER
    HETATM 1339  O23 DOG   174      15.551 -21.450  37.436  1.00105.39           O
    HETATM 1340  C23 DOG   174      14.900 -22.454  37.433  1.00105.39           C
    HETATM 1341  C22 DOG   174      13.769 -22.807  36.606  1.00105.39           C

こうして再び`tleap -f leap.in`コマンドを実行します。

    Total unperturbed charge:   0.000000
    Total perturbed charge:     0.000000
    Checking Unit.
    Building topology.
    Building atom parameters.
    Building bond parameters.
    Building angle parameters.
    Building proper torsion parameters.
    Building improper torsion parameters.
     total 556 improper torsions applied
    Building H-Bond parameters.
    Incorporating Non-Bonded adjustments.
    Not Marking per-residue atom chain types.
    Marking per-residue atom chain types.
      (Residues lacking connect0/connect1 -
       these don't have chain types marked:

            res     total affected

            CASN    1
            DOG     1
            NASP    1
            WAT     8527
      )
     (no restraints)
    Writing pdb file: leap.pdb
       printing CRYST1 record to PDB file with box info

    ~/apps/amber18/bin/teLeap: Warning!
     Converting N-terminal residue name to PDB format: NASP -> ASP

    ~/apps/amber18/bin/teLeap: Warning!
     Converting C-terminal residue name to PDB format: CASN -> ASN
            Quit

    Exiting LEaP: Errors = 0; Warnings = 2; Notes = 0.

こうして`Building topology.`と出力され、0バイトではない`leap.parm7`と`leap.rst7`ファイルが作成されていればトポロジーファイルの作成ステップは完了です。一応`leap.log`の中にWARNINGなどがないか確認しておき、あってもそれが大きな問題でなさそうかを確認しておきましょう。

さて、ここで作成された`leap.pdb`がどのようになっているか、分子描画ソフトウェアPyMOLなどで開いてみましょう。するとこのようになっているはずです。

![1lke4.png](https://qiita-image-store.s3.amazonaws.com/0/224327/364fb5b2-8ee4-128e-59ca-7b84445ef613.png)

図のように、タンパク質の周りに水分子、Naイオン、Clイオンが散りばめられている様子が見えます。特に、タンパク質の周りにたくさん存在する水分子は、ある種の見えない箱みたいなものに閉じ込められているかのように映るはずです。これは、今のleap.inの中で行っていたsolvateBoxによって、周期境界条件の溶媒和ボックスを生成したためです。

なぜこんな水ボックスを生成する必要があるかというと……現実のタンパク質1分子に対して水分子の占める体積は非常に大きいものになっていてほしいはずなのですが、コンピュータの計算の制約上、それをすべて表現するというのは困難です。そこで、無限ループする箱の中にタンパク質と水分子を入れるということで、擬似的に無限を表現することで一応の解決を得ています。図にするとこんな感じ。

![Pbc1.png](https://qiita-image-store.s3.amazonaws.com/0/224327/0e15f5dd-0524-f055-aac0-af3703a8bdfb.png)

この図で水色のボックスを先程生成したタンパク質と水分子の入った基準ボックスとしたとき、$x, y, z$方向に基準ボックスのコピーが無限に広がっているというイメージです。そしてボックス間で原子・分子の移動は可能となっており、ボックスの境界を超えた原子は見かけ上ボックスの反対側からまた流入するというような動きを見せます。また、静電相互作用のような長距離の力も適切に作用します。したがって、あまり小さな周期境界ボックスにタンパク質を入れてしまうと、各ボックスに存在するタンパク質のコピー同士間の静電相互作用が強まってしまいます。そこで上の`leap.in`では10 Åの余裕を持って周期境界条件ボックスを作るような設定にしていました（本当はもうちょっと大きな値にしておいた方がいいような気がしますがここではチュートリアルなので勘弁してください）。つまり、隣り合うタンパク質とは20 Å以上常時離れていることになります。この影響を小さくするためにはボックスを大きく取ればいいのですが、そうすると今度は水分子の個数が加速度的に増えていきますのでシミュレーションの計算速度が落ちてしまうことになり、そこはトレードオフとなってしまいます。

この周期境界ボックスのサイズは`leap.pdb`の1行目にかかれてあります。このファイルをテキストエディタで開いてみますと、

    CRYST1   71.492   72.959   68.898  90.00  90.00  90.00 P 1           1

という記述があり、$x=71.492, y=72.959, z=68.898$（Å）のサイズになっていることがわかります。


### Step 3. GROMACS形式への変換

Step 2でAMBER形式のトポロジーファイルが完成しました。このままAMBERでのMDシミュレーションを開始することができるのですが、AMBERでのMD計算はGROMACSに比べ計算が2〜3倍遅いので、私はGROMACS形式にトポロジーファイルを変換してからGROMACSでのシミュレーションを使うことを好んでいます（最近はともにGPU対応になり、あまり差がなくなりました）。

[GROMACS/ACPYPE.pyを用いたAMBER→Gromacs形式へのtopologyファイル変換](https://qiita.com/Ag_smith/items/6fd0afa8effc119b0f38)のページを参考にして、`~/apps`ディレクトリに`acpype.py`, `groconvert.sh`, `gen_posre.pl`の3つのファイルを置き、以下のコマンドで変換します。

```bash
# macOSのみ、gsed, gawkコマンドのインストールを行っておく (デフォルトのsed, awkは流派が異なり使えない機能があるため)
brew install gnu-sed gawk
~/apps/groconvert.sh -i leap -o 1lke -r
```

これで`1lke.gro`と`1lke.top`ファイルを作成します。また、`posreXX.itp`というファイルがいくつか作成されます。これは後のMDシミュレーションで使う、各重原子（＝水素以外の原子）に対しての位置拘束(position restraint)ファイルです。

`1lke.gro`ファイルを開いてみると、中はこのようになっています。

    leap_GMX.gro created by acpype (Rev: 0) on Fri Apr  6 11:52:56 2018
     28284
        1  ASP    N    1   2.444   1.727   4.364
        1  ASP   H1    2   2.350   1.760   4.349
        1  ASP   H2    3   2.442   1.627   4.377
        1  ASP   H3    4   2.504   1.768   4.294
    （中略）
     8732  WAT   H228281   1.518   1.755   0.175
     8733  WAT    O28282   1.421   1.671   1.497
     8733  WAT   H128283   1.500   1.622   1.477
     8733  WAT   H228284   1.450   1.739   1.558
        7.14920     7.29594     6.88980

GROMACSの座標ファイル形式である`.gro`の記述はシンプルで、1行目はコメント行（何を書いてもよい）、**2行目は系の中に存在する原子数（超重要）**、3行目以降は原子の記述（残基番号、残基名、atom type、原子の通し番号、x, y, z座標（**単位はÅではなくnm**））です。2行目の原子数は必ず実際にファイル内に存在する原子の数と一致させておく必要があります。3行目以降の記述は文字の列認識が固定されており、一文字でもずれると正しく認識してくれないことに注意です。ちなみに残基番号と原子の通し番号については実は飾りで、どんな数字を付けていてもエラーにはならず、上から順に認識します。そのため100000を超える原子を含むgroファイルであっても問題なく作動します。最終行は周期境界ボックスのサイズ指定であり、この指定は必須です。先程`leap.pdb`で見た値から1/10となっていますが、これはGROMACS形式がnm表記で統一していることと、AMBER形式ではÅ表記で統一しているために1/10ずれていることに起因します。なので問題ありません。

### Step 4. gromacs indexファイルの作成

後のMDシミュレーションのために、インデックスファイル（通常`index.ndx`）を作成します。MDシミュレーション内に含まれる原子は、タンパク質に属するものやリガンド、水もイオンも含まれます。また「水＋イオン」とか「（水＋イオン）以外」というグループを作りたいこともあります。このインデックスファイルを作っておくと、「リガンドだけの座標データを切り出して取得」「温度の制御を水イオングループとそれ以外に別々に設定」などということが可能になります。このチュートリアル方式では、MDシミュレーションを始めたら早速このインデックスファイルを要求されることになるので今のうちに作っておきます。

インデックスファイルを作るコマンドは`gmx make_ndx`です。

```shell
$ gmx make_ndx -f 1lke.gro
```

すると、続いて追加入力を求められます。

    There are:  8527      Water residues
    Analysing Protein...
    Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
    Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...

      0 System              : 28284 atoms
      1 Protein             :  2605 atoms
      2 Protein-H           :  1338 atoms
      3 C-alpha             :   169 atoms
      4 Backbone            :   507 atoms
      5 MainChain           :   677 atoms
      6 MainChain+Cb        :   830 atoms
      7 MainChain+H         :   842 atoms
      8 SideChain           :  1763 atoms
      9 SideChain-H         :   661 atoms
     10 Prot-Masses         :  2605 atoms
     11 non-Protein         : 25679 atoms
     12 Other               :    62 atoms
     13 DOG                 :    62 atoms
     14 Na+                 :    20 atoms
     15 Cl-                 :    16 atoms
     16 Ion                 :    36 atoms
     17 DOG                 :    62 atoms
     18 Na+                 :    20 atoms
     19 Cl-                 :    16 atoms
     20 Water               : 25581 atoms
     21 SOL                 : 25581 atoms
     22 non-Water           :  2703 atoms
     23 Water_and_ions      : 25617 atoms

     nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
     'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
     't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
     'r': residue              'res' nr         'chain' char
     "name": group             'case': case sensitive           'q': save and quit
     'ri': residue index

    >

ここでは`>`を求められるたびに、以下のように追加入力を行います。`q`を押すと終了です。

```
 23 Water_and_ions      : 25617 atoms

 nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
 'r': residue              'res' nr         'chain' char
 "name": group             'case': case sensitive           'q': save and quit
 'ri': residue index

> !23

Copied index group 23 'Water_and_ions'
Complemented group: 2667 atoms

 24 !Water_and_ions     :  2667 atoms

> name 24 nowation


> q
```

一度目は`!23`を入力します。これは`23: Water_and_ions`に含まれる原子の補集合を表す`!`をつけた設定です。これを指定すると上記のようなComplemented　group 2667 atomsが`24 !Water_and_ions : 2667　atoms`として出てきます（番号は系によって異なります）。この24番のグループを改名したい場合は`name 24　nowation`とします。最後に`q`を入力して終了です。

長かったですが、これでトポロジーファイルの生成は終了です。操作に慣れてくればここまでで10〜20分くらいで終わります。

この時点でのこの作業ディレクトリ`~/Desktop/1LKE/top`のファイル構成はこのようになっているはずです。`ls -lt`コマンドで表示してみますと

```shell
$ ls -lt ~/Desktop/1LKE/top
total 25512
-rw-r--r--  1 Agsmith  staff   833741  4  6 11:52 index.ndx
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre0.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre10.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre100.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre1000.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre20.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre200.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre50.itp
-rw-r--r--  1 Agsmith  staff    42407  4  6 11:52 posre500.itp
-rw-r--r--  1 Agsmith  staff  2056681  4  6 11:52 1lke.top
-rw-r--r--  1 Agsmith  staff  1272891  4  6 11:52 1lke.gro
-rw-r--r--  1 Agsmith  staff  1032459  4  6 11:50 leap.rst7
-rw-r--r--  1 Agsmith  staff   257958  4  6 11:50 leap.log
-rw-r--r--  1 Agsmith  staff  1955061  4  6 11:50 leap.pdb
-rw-r--r--  1 Agsmith  staff  5274680  4  6 11:50 leap.parm7
-rw-r--r--  1 Agsmith  staff     1294  4  6 11:45 leap.in
```

## GROMACSを使ったエネルギー最小化と平衡化

ここからは水溶液中に存在するタンパク質を物理演算で動かしていくことになります。しかし、これまでの操作で作り出した初期座標では水溶液の様子を再現しているとは言えません。その初期座標ファイル`1lke.gro`をPyMOLなどで見てもらえばわかりますが、このファイル内の水分子は規則正しく配置されたままですし、何より原子に速度情報が与えられていません。そこで、まずは系のエネルギーをある程度小さくしてから（＝エネルギー最小化）、各原子に（温度に対するBoltzmann-distributionに従う範囲で）ランダムに速度を与え、徐々に系全体を設定温度まで引き上げていく（＝平衡化）という予備動作を行うことで水溶液の状態に近づけていくという操作を行います。流れは

1.  [GROMACSを用いたエネルギー最小化](https://qiita.com/Ag_smith/items/3a438adb0dab4cfa06eb)
2.  [GROMACSを用いたNVT, NPT平衡化](https://qiita.com/Ag_smith/items/24fbb9467787f69fbf1c)

となり、後のProduction Runというステップにつなげていきます。

エネルギー最小化と平衡化の操作はほとんどすべてのMDシミュレーションにおいて使いまわすことができます。

さて、ここからは計算に時間がとてもかかることになります。しかし近年のGPUを使った計算手法の発展（GPGPU: GPUの演算資源を画像処理以外の目的に応用する技術のこと）のおかげで、GROMACSやAMBER、NAMDなどのMDシミュレーションのプログラムはGPUを使って計算を大幅に加速させることができます。その結果、GPUの取り付けてある計算機とそうでないマシンでは計算速度に雲泥の差が出てくるようになりました。以降の操作はGPUの有無にかかわらず可能ですが、GPUのあるマシンで行うことを推奨します。

上の最小化の操作を `~/Desktop/1LKE/gromacs/minimize` で行うものとし、ここに`min1.mdp`, `min2.mdp`とエネルギー最小化実行スクリプト`run.sh`を置き、NVT, NPT平衡化の操作は `~/Desktop/1LKE/gromacs/heat`で行うものとし、ここに`md1.mdp ~ md9.mdp`と平衡化実行スクリプト`run.sh`を置いて実行します。

## Production Runの実行

平衡化が終わったら、あとはもう時間とハードディスク容量が許す限りいくらでもMDシミュレーションを行うことができます。この平衡化以降のMDシミュレーションのことを**Production Run**と呼びます。実行ディレクトリを`~/Desktop/1LKE/gromacs/pr`に移します。

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

GROMACSの場合、Production Runの設定ファイル`mdrunset.mdp`は以下のようになります。

```:mdrunset.mdp
define            =
; Run parameters
integrator      = md            ; leap-frog integrator
dt              = 0.002     ; 2 fs
nsteps          = 50000000  ; 2 fs * 50000000 = 100 ns
cutoff-scheme   = Verlet
; Output control
nstxout         = 5000          ; save coordinates every 10 ps
nstvout         = 5000            ; save velocities every 10 ps
nstenergy       = 2500              ; save energies every 5 ps
nstlog          = 5000        ; update log file every 10 ps
; Bond parameters
continuation    = yes           ; Restarting after previous RUN
constraint_algorithm = LINCS    ; holonomic constraints
constraints          = hbonds  ; Convert the bonds with H-atoms to constraints
lincs_iter           = 1            ; accuracy of LINCS
lincs_order          = 4              ; also related to accuracy
; Langevin Dynamics
bd_fric = 2.0           ; Brownian dynamics friction coefficient(correspond to gamma_ln in AMBER)
ld_seed = 1732  ; random seed
; Neighborsearching
ns_type         = grid          ; search neighboring grid cells
nstlist         = 20              ; 10 fs
rlist           = 1.0       ; short-range neighborlist cutoff (in nm)
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.12      ; grid spacing for FFT(nm)
; Temperature coupling is on
tcoupl    = V-rescale   ; modified Berendsen thermostat
tc-grps   = nowation Water_and_ions   ; 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         = Parrinello-Rahman      ; 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
; Dispersion correction
DispCorr     = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel         = no            ; Velocity generation
gen_temp        = 300.0         ; temperature for Maxwell distribution
gen_seed        = 1732          ; used to initialize random generator for random velocities
```

特に重要なのは`dt = 0.002`と`nsteps`です。この値は結合長を拘束するSHAKEアルゴリズムやLINCSアルゴリズム（GROMACSでは`constraint_algorithm`で設定を切り替えられる）を利用している場合に、1ステップあたり2フェムト秒（`dt = 0.002`）後まで時間をワープさせても計算値が発散せずにシミュレーションを進めることが可能になっています。この繰り返しを規定するのが`nsteps`の値であり、掛け算した値がシミュレーション時間ということになります。現在、AMBER, NAMD, GROMACS、どのソフトウェアの設定をみても`dt = 0.002`としてあることが多いです（一応0.003までやってもいい方法が開発されているのですけど、これをやると計算の解析のときに微妙に切りの悪い数字になりがちなので結局使われることは少ないです）。

GROMACSは初めに`nsteps`で規定しておいたMDシミュレーション分の時間に達するまで計算し続けます。一応、後からでも`gmx convert-tpr`コマンドでこの時間を変更することが可能なのですが、ふつうは最初の時点でどれくらいのシミュレーション時間をしようか、よく考えてから設定しましょう。今回は100 nsのシミュレーションを行うことにし、`nsteps =　50000000`（5000万回ステップ計算と`dt = 0.002`の設定にしておきましょう。

温度制御と圧力制御はそれぞれ`V-rescale`法と`Parrinello-Rahman`法を使っています。これの詳細はものすごく物理学チックな話なので詳細は省きますが、圧力制御のParrinello-Rahman法を使うときは、すでに圧力がある程度落ち着いた状態から始めなければならないという点に気をつけておきましょう。詳細は[GROMACS/Production Runの方法](https://qiita.com/Ag_smith/items/1ebfd1accd4d3984ead2) を読んで下さい。

`mdrunset.mdp`ファイルを用意したら、この実行ファイル`run.sh`を用意してコンピュータ上（NVIDIA GPUがあるマシンを強く推奨）で実行します。

```bash
#!/bin/sh

input=1lke

# if md${input}.tpr doesn't exist, grompp will make the tpr file.
if [ ! -e md${input}.tpr ];then
    gmx grompp \
        -f mdrunset.mdp \
        -c ../heat/md9.tpr \
        -r ../heat/md9.tpr \
        -t ../heat/md9.cpt \
        -p ../../top/${input}.top \
        -n ../../top/index.ndx \
        -o md${input}.tpr \
        -po md${input}.mdp
fi

# gmx mdrun
gmx mdrun \
 -v -s md${input} -deffnm ${input} -cpi ${input}.cpt
```

例によって`gmx grompp`で先にMDシミュレーションの実行バイナリファイル`md1lke.tpr`を作り出します。今回は`-c ../heat/md9.tpr`と`-t ../heat/md9.cpt`で座標情報と速度情報をheatの最終ステップから読み込む、つまりheatの計算直後からデータを引き継いで計算を実行すると設定しています。またここではif文を使っており、このファイルが同じディレクトリ上にない場合には`gmx grompp`を実行する方式です。

これができたら続けて`gmx mdrun`でProduction Runを実行します。

```bash
# 作成したrun.shに実行許可のパーミッションを付与する
chmod +x run.sh
# 上で実行許可してあれば、そのプログラムは ./hoge という記法でターミナルから実行可能となる。
./run.sh
```

です。

この系を構成する原子数はおよそ28000ほどで、計算時間はNVIDIA GPUの1つ、GTX 1080Tiを積んだマシンでおよそ8〜10時間程度かかります。このような計算時間が長い場合には、計算途中でもしマシンがなんらかの理由でシャットダウンさせられてしまったり、もしくはユーザーの都合で計算を打ち切るというような事がありえます。そんなとき、計算を再スタートさせるときにまたProduction Runを一からやり直すというのはとても大変なことです。そこで、たいていのMDソフトウェアには計算の再スタートを円滑にするための**チェックポイント機能**が備わっています。Gromacsのデフォルト設定では15分に一度、`$input.cpt`という名前のチェックポイントファイルを作り出してくれます。このファイルにはその計算ステップ時点での時間情報、座標情報と速度情報を保存してくれているので、これさえ定期的に作り出していればもし計算が異常終了してもここから再スタートすることができるという仕組みです。

以上でProduction Runの説明は終了です。気の済むまでProduction Runをやりましょう。

ちなみにtprファイルを作成してしまった後で、「やっぱりMDのProduction Runの時間を変更したいな……」ってときは、 `gmx convert-tpr` コマンドでその時間を変更することができます。

```bash
# 上で作成したmd1lke.tprファイルの設定時間を300,000 ps (=300 ns)に変更したものをhoge.tprに書き出す
gmx convert-tpr -s md1lke.tpr -o hoge.tpr -until 300000
# hoge.tprをmd1lke.tprにリネーム（上書き）。
mv hoge.tpr md1lke.tpr
```

こうして作り変えた `md1lke.tpr` を作って再び上記 `run.sh` スクリプトを実行すれば300 nsまでのシミュレーションが行われるようになります。

## AmberToolsを使ったトラジェクトリファイルの解析

Production Runで得られたトラジェクトリファイル`1lke.trr`を解析するには再びAmberToolsを使います。特に、AmberToolsの`cpptraj`というモジュールが大活躍します。

`cpptraj`というモジュールはトラジェクトリの処理や各種座標情報の解析に用います。詳しくはAmberToolsのマニュアルのcpptrajの章を読んでください。amber22のマニュアルはamber22ディレクトリの下の`amber22_src/doc/Amber22.pdf`にあります。2023年5月以降はAmberTools23が登場したことで`Amber23.pdf`という名前になっています。

```bash
# /path/to/はamber22_srcがあるところに読みかえる。
# openコマンドはmacOSのみ
$ open /path/to/amber22_src/doc/Amber23.pdf
```

これでマニュアルが開きます。

上記の`cpptraj`コマンドはbashやpythonのようなインタプリタ言語式でスクリプトを記述します。基本的には`trajin`コマンドで

```shell:trajfix.in
## 1フレーム目だけを取り出す準備
# trajinで処理させたいトラジェクトリ名を指定。その後の1 1 1は、[最初のフレーム] [最後フレーム] [サンプリング間隔]いう意味。ここでは最初の1枚めのフレームだけ取り出す。
trajin 1lke.trr 1 1 1

# reference構造の指定
reference ../../top/leap.rst7
# 残基番号で1から170に対応する残基を、同じBoundary box内に位置するようにする（ボックスの外に行ってしまっても途中でループさせない）
unwrap :1-170
# 残基番号1-170のCα原子の重心を原点上にセンタリングするよう、各フレームのトラジェクトリを平行移動させる
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.dat @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
#以上の処理を施したものを1lke_init.pdbという名前で出力する
trajout 1lke_init.pdb pdb nobox
# 1フレーム目の取り出し実行
go
# これまでの入力情報をリセット
clear trajin

## 2段階目のtrajin処理の開始準備
trajin 1lke.trr 1 last 1
# もし複数のトラジェクトリファイルに分割している場合はこれに続けてtrajinで追加する
# trajin 1lke_2.trr 1 last 1
# trajin 1lke_3.trr 1 last 1
unwrap :1-170
trajout 1lkenj.trr
## 2段階目のtrajin処理実行
go
# これまでの入力情報をリセット
clear trajin

## 3段階目のtrajin処理の開始準備。
trajin 1lkenj.trr 1 last 10
unwrap :1-170
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.dat @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-

#以上の処理を施したものを1lkenj.trrという名前で出力する
trajout 1lkenj2.trr
## 3段階目のtrajin処理実行
go
```

このファイル`trajfix.in`を`~/Desktop/1LKE/gromacs/pr`ディレクトリに用意した後、以下のコマンドで実行します。`-p`にはtleapのときに作成したparm7ファイルを指定します。

```bash
cpptraj -p ../../top/leap.parm7 -i trajfix.in
```

これで`1lkenj.trr`と`1lkenj2.trr`という名前のファイルが出力されます。`1lkenj.trr`には周期境界条件を修正した状態、`1lkenj2.trr`には周期境界条件の修正と水・イオン分子の削除を行い、それを10ステップ刻みで出力した状態のトラジェクトリデータが含まれます。また、`1lke_init.pdb`ファイルにはMD開始時点の構造情報が入っています。以降の可視化では`1lkenj2.trr`の方を用います。

ちなみに`1lkenj.trr`の処理は必要でない場合もありますが、あらかじめ一度`unwrap`かつステップ刻みを1ずつに設定して出力しておいた方が、第2段階以降のトラジェクトリ変換に失敗しにくくなるので、あえてこのようにしています。

## PyMOLでトラジェクトリを表示する

タンパク質構造可視化ソフトウェアであるPyMOLを使ってMDシミュレーションの結果を可視化します。インストール方法は[macOS/Ubuntu 22.04へのオープンソース版PyMOL 2.5のインストール方法](https://qiita.com/Ag_smith/items/58e917710c4eddab46ee)この記事を参考にして下さい。

`1lke_init.pdb`と`1lkenj2.trr`というファイルを使って、手持ちのPCのPyMOL上でトラジェクトリファイルを表示させることで、タンパク質の動きを様々な角度から見ることができます。ここではPyMOLを用います。

まずPyMOLで`1lke_init.pdb`のファイルを開いた後、さらにその構造オブジェクトの上に`1lkenj2.trr`のトラジェクトリファイルをドラッグ＆ドロップで載せることでロードすることができます。コマンドで行う場合、PyMOLのコマンドウィンドウのところに`load 1lkenj2.trr, 1lke_init, 1`と入力することで可能です。

![Screenshot 2023-11-08 1.03.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/946fdebd-e654-a028-4844-e3bc008b0b9b.png)

## cpptrajで距離・角度・二面角と残基ゆらぎの値などを取得する

`cpptraj`モジュールはトラジェクトリについての情報を解析し、ファイルに出力させることができます。よく出力に用いられるのはトラジェクトリ中の2原子間の距離とか3原子間の角度、4原子間の二面角などです。

以下がそのサンプルコードです。HIE31のNE2原子とDOG170のO14原子の間の距離とDOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角を、トラジェクトリのフレームごとに出力します。また、`atomicfluct`で各残基のゆらぎ(**RMSF**)を計算させることができます。

```bash:analysis.in
trajin 1lkenj2.trr 1 last 1
# HIE31のNE2原子とDOG170のO14原子の間の距離をfoo.txtに出力
distance dist1 :31@NE2 :170@O14 out foo.txt
# MDシミュレーション中の各残基のゆらぎ（RMSF）を計算するatomicfluct
atomicfluct out rmsf.txt @CA byres
# DOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角の値をfoo.txtに出力
dihedral dih1 :170@C15 :170@C14 :170@O14 :31@NE2 out foo.txt
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
go
```

ところで、この`1lkenj2.trr`はすでに水・イオンを抜いた状態のトラジェクトリになっているため、その状態に対応するトポロジーファイル(.parm7)を用意してあげる必要があります。これはAmberToolsの`ante-MMBSA.py`で以下のようにして用意することができます。

```bash
cd ~/Desktop/1LKE/top/ # leap.parm7があるtopディレクトリに移動
ante-MMPBSA.py -p leap.parm7 -s :WAT,Na+,Cl- -c leap_stripped.parm7
```

これを行うと

```
Stripping :WAT,Na+,Cl- (solvent) from original topology, output is leap_stripped.parm7
Done stripping solvent!
```

と表示され、同ディレクトリに`leap_stripped.parm7`というファイルが生成されます。これは`ante-MMPBSA.py`のときに`-s`で指定した残基名を抜いた状態のparm7ファイルとなっています。

これらを用いて、再度`~/Desktop/1LKE/gromacs/pr`ディレクトリで上記`analysis.in`ファイルとともに

```bash
cpptraj -p ../../top/leap_stripped.parm7 -i analysis.in
```

を実行すればOKです。

ちなみに、先述のトラジェクトリを書き出すコード`trajfix.in`の中に組み込んでおいてそれにまとめることもできます。この場合は`leap_stripped.parm7`ではなく`leap.parm7`を用います。

```bash:trajfix.in
## 1フレーム目だけを取り出す準備
# trajinで処理させたいトラジェクトリ名を指定。その後の1 1 1は、[最初のフレーム] [最後フレーム] [サンプリング間隔]いう意味。ここでは最初の1枚めのフレームだけ取り出す。
trajin 1lke.trr 1 1 1

# reference構造の指定
reference ../../top/leap.rst7
# 残基番号で1から170に対応する残基を、同じBoundary box内に位置するようにする（ボックスの外に行ってしまっても途中でループさせない）
unwrap :1-170
# 残基番号1-170のCα原子の重心を原点上にセンタリングするよう、各フレームのトラジェクトリを平行移動させる
center :1-170@CA mass origin
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
#以上の処理を施したものを1lke_init.pdbという名前で出力する
trajout 1lke_init.pdb pdb nobox
# 1フレーム目の取り出し実行
go
# これまでの入力情報をリセット
clear trajin

## 2段階目のtrajin処理の開始準備
trajin 1lke.trr 1 last 1
# もし複数のトラジェクトリファイルに分割している場合はこれに続けてtrajinで追加する
# trajin 1lke_2.trr 1 last 1
# trajin 1lke_3.trr 1 last 1
unwrap :1-170
trajout 1lkenj.trr
## 2段階目のtrajin処理実行
go
# これまでの入力情報をリセット
clear trajin

## 3段階目のtrajin処理の開始準備。
## Gromacsは1フレーム目に初期情報が入っているので2フレームめ抽出すると良い。
trajin 1lkenj.trr 2 last 10
unwrap :1-170
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.txt @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
# HIE31のNE2原子とDOG170のO14原子の間の距離をoutput.txtに出力
distance dist1 :31@NE2 :170@O14 out output.txt
# MDシミュレーション中の各残基のゆらぎ（RMSF）を計算するatomicfluct
atomicfluct out rmsf.txt @CA byres
# DOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角の値をoutput.txtに出力
dihedral dih1 :170@C15 :170@C14 :170@O14 :31@NE2 out output.txt

#以上の処理を施したものを1lkenj.trrという名前で出力する
trajout 1lkenj2.trr
## 3段階目のtrajin処理実行
go
```

これを用いて同様に`cpptraj -p ../../top/leap.parm7 -i trajfix.in`コマンドを実行すれば、距離・二面角情報の書き出しとトラジェクトリの書き出しが同時に行われます。

出力された`output.txt`を見てみましょう。

```:output.txt
#Frame          dist1         dih1
       1       2.9622       0.4116
       2       2.9025      -4.1638
       3       2.8128     -30.7028
       4       2.8926     -17.0340
       5       2.9170     -20.7868
       6       3.0112      -4.1119
...
```

このように2つの値が並べて1000フレーム分（100ナノ秒）出力されています。これとPyMOLで`1lkenj2.trr`を表示してみたときの結果を表示してみると正しく値を取り出せていることも確認できるでしょう。

![Screenshot 2023-11-08 1.00.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/7ae27447-7eee-cab9-b650-a71fd196fdd0.png)


また、`rmsd.txt`の方もこうなっています。

```:rmsd.txt
#Frame     RMSD_00002
       1       0.0000
       2       1.2040
       3       1.2095
       4       1.2948
       5       1.3648
       6       1.6311
       7       1.8407
       8       1.8810
       9       1.8093
      10       1.7317
...
```

これは後ほどRMSDの時間発展として使うことができます。

RMSFの値は`rmsf.txt`に出力されます。

```:rmsf.txt
#Res        AtomicFlx
   1.000       3.9846
   2.000       3.0131
   3.000       2.1024
   4.000       1.4635
   5.000       2.1263
   6.000       1.8951
   7.000       1.4176
...
```

## matplotlibで結果を図にする

`cpptraj`で得られる角度情報をmatplotlibでプロットしてみましょう。この技術をマスターすれば、論文に載せられる図も簡単に作成することができます。

例として以下の`rmsd.txt`ファイルのRMSDの時間変化を図にすることを考えます。

```:rmsd.txt
#Frame     RMSD_00002
       1       0.0000
       2       1.2040
       3       1.2095
       4       1.2948
       5       1.3648
       6       1.6311
       7       1.8407
       8       1.8810
       9       1.8093
      10       1.7317
...
```

以下の`plot.py`でプロットすることができます。当然python3にmatplotlibとnumpyがインストールされていることが条件です。まだインストールしていない方はターミナルから`pip`でインストールしておきましょう。

```bash
$ python3.11 -m pip install numpy matplotlib
```

```python:plot.py
#!/usr/bin/env python3.11

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams.update(
    {
        "text.usetex": False,
        "font.family": "sans-serif",
        "font.sans-serif": "Arial",
        "font.size": 18,
        "lines.color": "b",
        "lines.linewidth": 1.0,
    }
)

# 入力ファイル名
file = "rmsd.txt"

# スペース区切りのテキストファイルからデータセットを生成
# "#"で始まる行はコメント行として無視する。
# 1列目をX軸、2列目をY軸とみなしている。
x, y = np.loadtxt(file, comments="#", unpack=True)
# 1000フレーム分のデータを100nsに置き換えるため、x軸の値に0.1をかける
x = 0.1 * x

# plot
fig = plt.figure()
fig1 = fig.add_subplot(111)
fig1.plot(x, y, label="rmsd")

# x軸の設定
fig1.set_xlim(0, 100)  # 表示範囲
fig1.set_xlabel("Simulation time [ns]")  # ラベル
# y軸の設定
fig1.set_ylim(0, 10)  # 表示範囲
fig1.set_yticks([0, 2, 4, 6, 8, 10])
fig1.set_ylabel(
    r"RMSD [$\AA$]",
)
fig1.legend()
plt.show()

# save
fig.savefig("rmsd.png", bbox_inches="tight", pad_inches=0.05)
```

これを実行すると`rmsd.png`という画像ファイルが出力されます。

![rmsd.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/178a9025-1571-2f92-5c71-652e227cb019.png)

ちなみに`"text.usetex": False,`を`True`に、fig.savefig時に`.pdf`拡張子を指定するとベクター画像形式になります。後でこのファイルをIllustratorなどで開いて目盛りなど見た目を再調整したいときはこちらが便利でしょう。

ここから見た目の変更などを行いたい場合は、matplotlibのチュートリアルとかAPIリファレンスとかを読み、上で何をしているかを理解して思い通りのフォーマットを作ってあげましょう。

また、RMSFを同様にしてプロットしてみます。$x$軸は残基番号、$y$軸がRMSF[Å]になります。

```python
#!/usr/bin/env python3.11

# %%
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams.update(
    {
        "text.usetex": False,
        "font.family": "sans-serif",
        "font.sans-serif": "Arial",
        "font.size": 18,
        "lines.color": "b",
        "lines.linewidth": 1.0,
    }
)

# 入力ファイル名
file = "rmsf.txt"

# スペース区切りのテキストファイルからデータセットを生成
# "#"で始まる行はコメント行として無視する。
x, y = np.loadtxt(file, comments="#", unpack=True)

# plot
fig = plt.figure()
fig1 = fig.add_subplot(111)
fig1.plot(x, y, label="rmsf")

# x軸の設定
fig1.set_xlim(0, 170)  # 表示範囲
# np.arrangeで0から170まで20刻みの配列を作成
fig1.set_xticks(np.arange(0, 170, step=20))
fig1.set_xlabel("Residue Number")  # ラベル
# y軸の設定
fig1.set_ylim(0, 10)  # 表示範囲
fig1.set_yticks([0, 2, 4, 6, 8, 10])
fig1.set_ylabel(
    r"RMSF [$\AA$]",
)
fig1.legend()
plt.show()

# save
fig.savefig("rmsf.png", bbox_inches="tight", pad_inches=0.05)
```

結果はこんな感じ。

![rmsf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/224327/077eed80-1eab-03e0-45c4-923e367638d9.png)

N末端とC末端の残基のゆらぎがそれぞれ大きいことは、PyMOLで表示させたトラジェクトリの動きと対応していることでしょう。また、アミノ酸Asp114からLys118までの"DEDKK"の周辺でRMSFの値が大きくなっていることがわかります。これは結晶構造中で**ミッシング残基**になっていた範囲の残基と一致しており、この領域が本質的にゆらぎが大きく、そのために結晶構造中でも残基の座標が決定できなかったことと一致しています。