#はじめに
本シリーズは下記チュートリアルを実施したメモです。
本記事は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の指定を楽にする方法を学びます。