PLUMED入門 その３ 簡単なCVを計算する

#はじめに
本シリーズは下記チュートリアルを実施したメモです。

https://www.plumed.org/doc-v2.6/user-doc/html/master-_i_s_d_d-1.html

本記事はPLUMED入門の第３回です

過去回はこちら↓

https://qiita.com/tacoma/items/8b445de0a716477843f8

https://qiita.com/tacoma/items/0da07516dab6a7e2da2d


#使用するファイル

チュートリアルで使用するファイルはこちら↓

https://github.com/maxbonomi/plumed-nest-data/raw/master/projects-2020/master-ISDD-1.tar.gz

```bash
tar xzvf master-ISDD-1.tar.gz
```
で解凍する

中身は、

* GB1_native.pdb：GB1の天然構造
* traj-whole.xtc：xtcフォーマットのトラジェクトリー。PBCで分子が切断される問題が解消されている
* traj-broken.xtc：こちらも同じくトラジェクトリーだが、GROMACSが出力したものでありPBCで分子が切断されている


# Exercise 1 : 簡単なCVを計算、出力する
計算済みのトラジェクトリーからCVを計算するために

* PLUMEDのインプットファイル（plumed.dat）を作成する
* PLUMED driverユーティリティを実行する

VMDでトラジェクトリーを可視化することもできる

```
vmd GB1_native.pdb traj-whole.xtc
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/f730c462-3daa-48d4-79d2-0b9721a9bd35.png)

GB1プロテインがフォールディングしていく


PLUMEDのインプットファイルを用意しよう
以下を計算したい

* すべてのCA原子（Cα原子）を使った慣性半径（gyration of radius）
* すべてのCA原子間のcontactを計算する。contactの定義は距離R_0が8.0Å以下になることとする

下記のサンプルをテンプレートとして使うと良い
FILLとなっている箇所は、書き換える必要がある

```python
# Compute gyration radius on CA atoms:
r: GYRATION ATOMS=__FILL__

# Compute number of contacts between CA atoms
co: COORDINATION GROUPA=__FILL__ R_0=__FILL__

# Print the two collective variables on COLVAR file every step
PRINT ARG=__FILL__ FILE=COLVAR STRIDE=__FILL__
```
(注：チュートリアル上ではFILL部分は自分で考える形式になっていて正解は書かれていません。私が書いたインプットは下記です)

```python
# define ca atom group
ca: GROUP ATOMS=2,10,17,29,38,46,54,62,70,74,83,90,98,107,111,120,127,134,141,150,155,162,170,175,180,187,192,201,210,217,228,237,246,258,263,271,279,287,291,298,306,310,319,333,340,352,360,368,373,380,389,396,407,414,421

# Compute gyration radius on CA atoms:
r: GYRATION ATOMS=ca

# Compute number of contacts between CA atoms
co: COORDINATION GROUPA=ca R_0=0.8

# Print the two collective variables on COLVAR file every step
PRINT ARG=r,co FILE=COLVAR STRIDE=1

```
（なお、下記サイトを参考にしました）

https://www.plumed.org/doc-v2.6/user-doc/html/_g_r_o_u_p.html

Cα原子のインデックスを求めるのに、GROMACSのインデックスファイルを使いました。

```bash
gmx make_ndx -f GB1_native.pdb
> q
```

index.ndxが出力されます。
この中のC-alphaグループがCα原子のインデックスです。

```
[ C-alpha ]
   2   10   17   29   38   46   54   62   70   74   83   90   98  107  111
 120  127  134  141  150  155  162  170  175  180  187  192  201  210  217
 228  237  246  258  263  271  279  287  291  298  306  310  319  333  340
 352  360  368  373  380  389  396  407  414  421
```

さて、PLUMEDを実行してみましょう。

```
> plumed driver --plumed plumed.dat --mf_xtc traj-broken.xtc
```

ログの途中に以下のような文面があり、PBCでbrokenした分子も修復されていることがわかります。

```
broken molecules will be rebuilt assuming atoms are in the proper order
```
以下の記述は何でしょうか？

```
between two groups of 55 and 0 atoms
```

これはGROUPBというオプションがあり、もし何も指定しないとGROUPAの中での組み合わせを計算するということだそうです。


出力された"COLVAR"ファイルを見ると、rおよびcoが出力されていることが分かる

```
#! FIELDS time r co
 0.000000 2.457214 162.290544
 1.000000 2.302935 161.701507
 2.000000 2.351218 160.458395
 3.000000 2.412862 141.860172
 4.000000 2.546138 143.713653
 5.000000 2.274152 161.323902
 6.000000 2.099441 173.566495
 7.000000 2.037371 181.014834
 8.000000 2.000826 179.489130
 9.000000 2.113100 154.312681
```

gnuplotで確認してみよう

r
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/67060ced-f667-6c9d-db19-137b985f4c34.png)

co
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/788299/b87fa880-5479-caea-537e-a92022e950fe.png)

余裕があれば、もう一度コマンドを実行して何が起こるか確かめよう。（前の出力ファイルのバックアップが作成される）
また、traj-whole.xtcを実行してみよう。結果は同じになるはずだ。なぜか？（brokenのときにもプログラムが修正してくれるから）

次回はMOLINFOショートカットを使って、CVの指定を楽にする方法を学びます。