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

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

過去回はこちら↓

#使用するファイル

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

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

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

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

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

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

# 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部分は自分で考える形式になっていて正解は書かれていません。私が書いたインプットは下記です)

# 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

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

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

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

co

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

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