うみうしの日誌

~(・ ω ・)~nullnull

アンブレラサンプリングとPMF計算(2)

アンブレラサンプリングとPMF計算の続き
とりあえず出来たのだが、だいぶ前になったので忘れてきている。実行したのも研究室のパソコンなのでインプットファイルが手元になく、少し嘘を書くかもしれない。間違いに気づいたら随時修正する。
今回も引き続きJustin A. LemkulさんのGromacs tutorialsを完全に追う形で進めて、個人的に躓いたところをメモする。

ちょっと解説とメモ
前回までで反応座標が決定された。(という言い方で正しいかは分からないが。)
実際ちゃんとやってみるまでアンブレラサンプリングがどんなことをしているのか分かりづらかったので、メモも兼ねて自分の理解をここで少し簡単に説明してみる。

アンブレラサンプリングは、最終的に反応座標に沿って自由エネルギー(具体的に何のかはまだ分かっていない。ギブズの?なんか他の?)を得る方法である。自由エネルギー地形という表現をしている人もいた。
実際どういう風に計算するかというと、まず前回までにやった通り反応座標のある地点(z軸方向に引っ張ったならならz=1とか1.5とか。英語だとよくwindowとか表記されている)に目的の分子がいる構造のファイル(.groと.tpr。他にも何か必要だったかな?とりあえず特に指定していなくても必要なファイルは出力されているはず)を得る。この時点でクラッシュしたり止まっていたりする場合、前回の.mdpで設定したpull_coord1_rateが速すぎるか、dtが大きすぎるか何かだと思われる。
で、この時点では無理やり引っ張っていっているので、各反応座標にいる分子が最安定な状態になっていない=その反応座標の求めたい自由エネルギーの状態に落ち着いていない(という言い方であってるかは知らない)
雑なイメージは下図。

f:id:CHfour:20171212221143p:plain
青:求めたい自由エネルギー 赤:前回までで得た構造のファイルの自由エネルギー

この赤い線を青い線まで落とすときの計算手法がアンブレラサンプリング、であっているはず。ここまでの一連の計算とその後の処理もあわせてアンブレラサンプリングと言うのかも。
赤い線を青い線までどうやって落とすかと言えば、50個ぐらいwindowをとって、各windowについてその座標から分子が逃げ出さないようにしながらMDを走らせて安定構造を得る、のだと思う。Justin A. LemkulさんのGromacs tutorialsの下の図はそれを多分表している。

f:id:CHfour:20171212214719j:plain
アンブレラサンプリングの概念図 GROMACS Tutorial > Umbrella Sampling > Step Five: Generating Configurationsより

その後WHAMという方法で、分子が逃げ出さないようにかけていた拘束?重みづけ?の影響を除きながら青い自由エネルギー地形を得ることになる。


3.アンブレラサンプリング

そういうわけで、

  1. 前回MDシミュレーションを走らせて得た.tprファイルと.xtc(または.trr)ファイルから、各windowの初期構造となる.groファイルを得る
  2. 各windowについてNPT平衡化する
  3. 各windowについてMDシミュレーションを走らせる

ことが必要になる。

1.groファイルを得る
Justin A. LemkulさんのGromacs tutorialsに便利なプログラムがおいてあるので適宜活用する。
まず

gmx trjconv -s pull.tpr -f traj_comp.xtc -o conf.gro -sep

と打ち込んで(pull.tprとtraj_comp.xtcは自分が前回までで得たファイル名を指定)、xtcファイルをばらしてgroファイルを得る。Justinさんのところでは501個のファイルが得られると書いてあるが、これは前回自分が走らせた時のMDシミュレーションの長さとどれぐらいの頻度でlogを書き出したかによって変わる。
長さは(dtと)nstepsのところで、頻度は(確か)nsthogehogeのところで何ステップごとに書き出すかを変更できる。あまり頻度が高いと多分計算が重くなるが、頻度が低いとこの後で詰むのでそこそこの頻度で書き出すようにする。そこそこってどれぐらいだよって思われるだろうが、これは今のところ僕にはこの後作業してみて詰んだら頻度を上げてやり直す、以上の方法は分からない

ここで得たすべてのgroファイルについてシミュレーションを走らせる必要はなく、十分近い間隔のファイルだけを抜き出して初期構造とすればよい。JustinさんはTutorialとして反応座標が(大体)0.2nm変わるごとに抜き出している。
この反応座標とはJustinさんの場合は多分AβのChainAとChainBの重心間の距離のはず。そのため、各conf1~500.groのファイルについてChainAとChainBの重心をそれぞれ計算し、重心間距離を求め、それが0.2nm変わっていたgroファイルを見つける必要がある。はっきり言って手作業では不可能なので、プログラムを書いて実行する。
Justinさんのところではそのプログラム(スクリプト?)が公開されていて、Perlで実行できる。Perlwindowsでも無料で利用できるパッケージというかなんというかググってもらえば見つかるやつがあるので(Strawberry Perlとか)、入れておく。

ここで公開されているdistance.plは重心間距離を計算するので、前回のpull_codeでpull_coord1_geometry = distanceを指定した場合はそのまま使えるが、directionやdirection-periodicの場合は少し書き直して使う必要がある。具体的には

for (my $i=0; $i<=500; $i++) {
    print "Processing configuration $i...\n";
    system("gmx distance -s pull.tpr -f conf${i}.gro -n index.ndx -oall dist${i}.xvg -select \'com of group \"Chain_A\" plus com of group \"Chain_B\"\' &>/dev/null");
}

のところを編集し、z軸方向の変位を得るようにする。編集するのだが、流石に人様が書いたスクリプトを勝手に書き換えたものを再配布するのはギルティなので、ここは申し訳ないが各自でなんとか頑張ってほしい。GROMACS公式のgmx distanceのページを読んであれこれしたら何とかなる。なったはず。gmx trajを使っていた。自分で編集したけど、何をしているのか今の自分には分からない。

これが出来たら、summary_distances.datというファイルが作成されていて、各0~500.groにおける目的の分子の反応座標(JustinさんのならChainAとChainBの重心間距離。書き換えていたらz軸座標)が一覧になっている。この反応座標を確認し、大体0.2ずつ値が変わっていくようにconfhoge.groをいくつか選ぶ。
ここで前回分子を引っ張ったときのmdpファイルの設定でlogの書き出しの頻度が低かったりpull_rateが速い場合、1つ次のgroファイルに行っただけで0.3とか0.4とか反応座標が変化してしまったりする。とりあえず一番変位が小さいファイルを選んで次に進んでもよいが、おそらく最後まではいけないのでpullingからやり直すほうがよい。これが初めてならとりあえず詰まるとこまでやってみるのもあり。ただこのあとかなり時間のかかる計算をするので、ここで引き返すほうが無難。

2.NPT平衡化
初期構造となるgroファイルが選べたら、

gmx grompp -f npt_umbrella.mdp -c conf0.gro -p topol.top -n index.ndx -o npt0.tpr
...
gmx grompp -f npt_umbrella.mdp -c conf450.gro -p topol.top -n index.ndx -o npt22.tpr

を実行。上はJustinさんのところの例で、confhoge.groの数字とnpthoge.tprのところは自分が選んだファイルに依存して変える。mdpファイルは基本的に前回のときのNPT平衡化のものでよい。変更点として末尾にpull_codeを加え、pull_rateを0.0にする。あとCHARMM-GUIで最初のファイルを作成していた場合topファイルの名前がsystem.topになっているので、topol.topに名前を変更しておくとよい。

これができたら

gmx mdrun -deffnm npt0
...
gmx mdrun -deffnm npt22

を実行する。これでNPT平衡化が出来る。
僕の場合CHARMM-GUIで得たstep6.0~6.6_equilibration.mdpのどこからどこまでがNPT平衡化にあたるのか自信がなかったため、integrator = mdになるstep6.2から6.6まで行った。溶媒の平衡化は改めてやる必要ないだろうし、こんなにたくさんの計算は必要はなかったかも。

3.アンブレラサンプリング
tprファイルを作る。

gmx grompp -f md_umbrella.mdp -c npt0.gro -t npt0.cpt -p topol.top -n index.ndx -o umbrella0.tpr
...
gmx grompp -f md_umbrella.mdp -c npt22.gro -t npt22.cpt -p topol.top -n index.ndx -o umbrella22.tpr

例によって上はJustinさんのところの例。なぜcptファイルが要るのかは今の僕には理解できない。mdpファイルは前回引っ張った時のstep7_production.mdpとほぼ同じでよい。pull_rateを0.0にし、適宜dt、nstepsを変更する。
で、各umbrella0~22について以下を実行。

gmx mdrun -deffnm umbrella0 -pf pullf-umbrella0.xvg -px pullx-umbrella0.xvg
...
gmx mdrun -deffnm umbrella22 -pf pullf-umbrella22.xvg -px pullx-umbrella22.xvg

ここがかなり時間がかかる。window数だけMDシミュレーションを走らせることになるので。僕がやった時は12000個の粒子を含む系で、50windowについて1windowあたり30ns(めちゃくちゃ短い)走らせたため、1日180ns計算出来るパソコンだったが1週間強かかった。これを2種類の分子についてやってみたので合計3週間かかった。CPU温度が上がりすぎてkernelエラーが頻発するなど色々苦労したが、そこはマシンスペック次第。

ここまでの作業を全て自動化したPythonスクリプトがJustinさんのところで公開されている。tar.gzか何かの形式で、これを解凍すると
caught-output.txt
frame-hoge_run-umbrella.shがたくさん
README
run-umbrella.sh
setupUmbrella.py
summary_distances.dat
が入っていると思う。自分はPythonの知識がほとんどないのでsetupUmbrella.pyがどういう実装なのかは全然理解出来なかった。caught-output.txt、frame-hoge_run-umbrella.sh、run-umbrella.shはこのスクリプトで生成されるファイルなので、ここに入っているファイルはあくまでJustinさんの例のファイルになる。
自分で使う場合は自分のsummary_distances.datを用いてsetupUmbrella.pyを実行し、生成されたrun-umbrella.shを実行すればよいはず。
が、pull_coord1_geometry = directionを指定した場合はこのままは使えない。僕にはこのスクリプトを改変する技術がなかったので、泣く泣く上記のように手動で全て行った。

4.重み付きヒストグラムの解析(gmx wham)

上のシミュレーションが終わったら、あとは自由エネルギー地形を得るだけ。
基本的にJustinさんのところの日本語訳でしかないが、まず新規ファイルを作成し

umbrella0.tpr
umbrella1.tpr
...
umbrella22.tpr

とだけ記入してtpr-files.datとして保存する。同じようにファイル名の一覧を記したファイルについてpullf-files.dat(欲しければpullx-files.datも)を用意する。

で、以下を実行

gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal

ここでなんか上手くいかないことがある。
まず一つはいつまでも平衡化?の段階が終わらないことがある。多分系が不安定だってことになるので、よくわからないがダメな気がする。気がするが、toleranceか何か(何だったか忘れたがコマンドラインか端末上に示される)を引数として少し大きい値を与えて終わらせてしまうことが出来る。
もう一つはwindow間が離れすぎていて、自由エネルギー地形を復元できないことがある。これはこのままではどうしようもないので、離れすぎていたwindow間にwindowを足してやる。端末上でどこの間に隙間があるかは教えてくれた気がする。そうでなくても、gmx whamを実行した時点でhisto.xvgが生成されているはずで、そのヒストグラムを確認すれば分かる。LINUXならxmgrace -nxy histo.xvgでgraceがグラフを描いてくれる。(-nxyはなくてもよい)
Justinさんのところの例なら下のような感じ。

f:id:CHfour:20171213002600j:plain
histo.xvg GROMACS Tutorial > Umbrella Sampling > Step Seven: Data Analysisより

この図の山が重なっていない部分があるとgmx whamは上手くいかないので、離れている山があったらその間の座標にあたるconfhoge.groファイルをとってきてNPT平衡化、MDrunを行う。もし間にあたる座標がない場合は、前回の引っ張るところからやり直しになる。ヤバい。

ここまでが全て上手くいけば、無事Justinさんのところのような自由エネルギー地形のグラフが得られたはず。自分が参考にした論文ではこれと同じことを、前回の一番最初のところから14回シミュレーションしてその平均と標準誤差を取ってグラフを出していた。
そこまでしなければならないのかは不明。

振り返ってみるとかなりGromacs Tutorialのままの内容なので、引用の範囲を超えてるかな……
後でまずいと思ったら消すかも。