うみうしの日誌

~(・ ω ・)~nullnull

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

Gromacsにできることの一つに、アンブレラサンプリング法というものがある。詳しい説明はググってもらいたい。
私の雑な理解としては、時間発展を諦める代わりにある反応座標(2集団の重心の距離とか、座標とか)におけるエネルギーが計算できる。このエネルギーというのはあんまり分かっていないが、Potential Mean of Forceの略でPMFというようだ。多分自由エネルギー的な何かなんだと思っている。
これで何が分かるかというと、例えばGromacs初心者が必ず見るだろうJustin A. LemkulさんのGromacs tutorialsでは、Aβの重合体から1分子を引き抜く際のPMF?を計算して、binding energy (ΔGbind)が計算できると言っている。凄い。

そういうわけで私も使ってみたかったが、PCスペックの都合もあり今回粗視化モデルで計算することにしたので、その辺りも含めて記録を残しておくことにする。
ただ、今回の内容はかなり私の雑な理解による説明が多いので、きちんと裏付けをとるべきところが多数あると思う。
正直有識者の方でおかしな点に気づいたら指摘してほしい。独学は大変だ。

0.粗視化モデル(MARTINI)について

気が向いたら少し書く。基本使って覚えるスタイルなので、理論やモデルについてはきちんと説明できるほど理解が及んではいないのだ。


【MARTINIについて】
『Martini Coarse Grain Forcefield for Biomolecules』

1.Gromacs用のファイルを手に入れる

たここに頼る。
CHARMM-GUI provides a web-based graphical user interface to generate various molecular simulation systems and input files to facilitate and standardize the usage of common and advanced simulation techniques. Currently, CHARMM-GUI supports CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM simulation programs mostly based on the CHARMM force fields.

CHARMM-GUI

Martini Makerから、MARTINIを利用した系を作ってくれる。
Membrane Builderでタンパク-脂質膜の系を立てるときとの違いとしては、(何か特殊なことが出来る人は違うかもしれないが)
・基本的に非天然アミノ酸は使えない
・C,N末端の修飾もできず、ChargedかNeutralかだけ選択できる
ことなどがある。私は一部非天然アミノ酸を使いたかったが、素人では太刀打ちできなかったので諦めて性質の似た天然アミノ酸残基に変更した。


【CHARMM-GUIの使い方例(再掲)】
『膜タンパク質のモデリング (CHARMM-GUI, Amber Lipid14 force field)』
香川大学農学部 ケミカルバイオロジー研究室 > 計算機 > の記事。

2.シミュレーションを走らせる

ここからが本番で、CHARMM-GUIで作った系では単に時間発展を計算するだけなので、このままではアンブレラサンプリングは出来ない。そのため、ここで得たmdpファイルをちょこちょこ修正する必要がある。

さて前にも触れたがcharmm-guiからは以下のファイルが得られる。

 topparフォルダ
 いくつかのitpファイル
 system.topとindex.ndx
 step5_assembly.box.pdb
 step6.0_equilibration.mdp~step6.6_equilibration.mdp
 step7_production.mdp
 README

このequilibration.mdpだが、おそらく全てNPT平衡化を行っている、はず。position restraintを少しずつ変えているだけ、なはず。
N(多分粒子数)P(圧力)T(温度)を一定に保つように、なんか上手くsimulation boxの体積を変えながら少し時間発展させて、本格的なMDシミュレーションを走らせる前に系を安定化しているんだと思っている。これをせずにいきなり走らせると、初期構造がおかしいせいで最小エネルギー構造とは別の極小値にはまる可能性があるんだったと思う。僕はAG++さんのとこで勉強したので適当にググって調べるとよい。

で、平衡化の過程は変える必要はないのでそのまま行う。何か参考にしている論文があって、その条件に揃えたいなら適宜mdpファイルを編集する。
変える必要があるのはstep7_production.mdpで、ここに
; Pull code
hogehoge
を追加する必要がある。

Pull codeで何をするかというと、もうほんと雑な理解だが、アンブレラサンプリングではまず指定したグループ(多分index.ndx内で分かれているもののこと。Martini Makerで脂質膜とタンパクの系を立てたなら、ProteinとMembraneとSoluteがあると思う。)を指定した方向に引っ張る計算を走らせる過程が必要なようだ。この過程で、どういう風に引っ張るか、何を引っ張るかとかのルールを、Pull codeで定めている。多分。

Justin A. LemkulさんのGromacs tutorialsで紹介されているPull codeは以下のようになっている。

; Pull code
pull = yes
pull_ngroups = 2
pull_ncoords = 1
pull_group1_name = Chain_B
pull_group2_name = Chain_A
pull_coord1_type = umbrella ; harmonic biasing force
pull_coord1_geometry = distance ; simple distance increase
pull_coord1_groups = 1 2
pull_coord1_dim = N N Y
pull_coord1_rate = 0.01 ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_start = yes ; define initial COM distance > 0
Gromacs tutorials
詳細も全て上のページに書いてあるし、ここと違う使い方がしたい場合は本家Gromacsのマニュアルページに飛ぶと色々書いてある。

私の理解を簡単に説明すると
pull:yesかnoだけ。デフォだと多分no。yesにしないと始まらないのでyesにする。たまにGromacsのメーリングリストでここをumbrellaにしてるために上手くいっていないだろう人を見かけるので、人のファイルを参照する場合(特に、「こういうファイル使っているけどうまく走らないんだ」と言っている人のファイルを間借りする場合)は気を付けてみる。
pull_ngroups:引っ張るグループの数?デフォでは1のようだ。ここで2になっているのは、Tutorial内では5分子のAβが会合したもののうちChainAだけを引っ張って外したいようだが、ただ引っ張ると重合体全体が引っ張られてしまって引き離せないので、隣のChainBを参照グループとして使用する?ためらしい。実際Tutorial内ではtopol_B.itpを別に用意し、#ifdef POSRES_BでChainBに重みづけをしているようだ。多分。
pull_ncoords:反応座標の数。1以外のやり方はわからない。
pull_group1_name:ndxファイル内のグループ名で指定する。
pull_coord1_type:umbrella以外知らない。色々あるらしい。
pull_coord1_geometry:引っ張る方向の指定方法の指定(ややこしい。)distanceなら、指定した2つのグループの重心を近づける方向に引っ張るはず。これは重心が重なったら止まる。それはそう。
私はそれが少し都合が悪かったので、direction-periodicを使った。
pull_coord1_dim:引っ張る方向の指定。普通に中学高校の数学の座標と同じで、(X,Y,Z)についてYesかNoを指定する。上の例の場合だとZ軸方向だけYesなので、Z軸方向だけに引っ張り力が働いているはず。多分。
pull_coord1_rate:よくわからないが、引っ張る速度を指定している。マイナスにするとちゃんと逆方向に引っ張れる。0.01でも速度としては速めらしいが、Justin A. Lemkulさんらは0.005と0.0001でもやってみて、時間-エネルギー(多分引っ張る力の)グラフを書いたら全て概形が同じだったので、0.01で問題ないとした、と言っていた。多分実際使うときもそうするのがよい。私は計算コストが惜しかったので確認なく-0.002を使った。-0.01は速すぎたようだったが検証はしてないし大丈夫だったかも。
pull_coord1_k:分からない。手を触れないことにした。おそらく参照グループを抑える力か、引っ張る力かのどっちか
pull_coord1_start:よく分からない。yesだと良いっぽかったのでyesにした。

さて、今回私がやりたかった内容はCPPが脂質膜を通過するときのエネルギーだったのだが、非対称膜を作ったのでCOMが重なったところでCPPが止まるのは正直不都合だった。参考にしていた論文ではpull_coord1_geometry=distanceのまま計算をしているようだったが、私はdirection-periodicにした。
ここが結構やっても大丈夫だったのか怪しいところ。pull_coord1_geometryにはdistanceの他に、単純に指定した軸方向に引っ張るdirectionというモードがある。このモードのとき、シミュレーションボックスの長さの半分以上の距離を分子が動くとシミュレーションは止まってしまう。多分何かしらまずいことが起きることが分かっているのだろう。この場合、例えば溶媒の量を増やしてシミュレーションボックスを長くし、動かしたい距離の2倍以上の長さになるようにするのがよりらしい。だがそんなことすると計算コストが上がるのであまりやりたくなかった。
そのときGromacsさんは代わりにdirection-periodicを使用することをコマンドライン上で勧めてきた。これなら箱の半分以上の長さを動いても良いらしい。
ただdirection-periodicにも制限があり、引っ張りたい方向についてはシミュレーションボックスの長さが固定されている必要がある。これは何かというと、系の圧力を一定に保つ際、上でも触れたがシミュレーションボックスの体積を少しずつ変えているのだ。これを、引っ張りたい方向については行わないようにする必要がある。
Justinさんがメーリングリスト上で答えていたところによると、NPT平衡化を諦めてNVT平衡化だけにする?か、Pressure couplingの内容を変更し、pcoupltypeをsemiisotropicに、そしてcompressibilityを"(何かの値) 0"にする必要がある。semiisotropicでは体積を変更する際に、XYZ方向に等しく変化を加えるのではなく、XYとZ方向に分けて変化を加える。で、compressibilityではXY方向についてとZ方向についての2つの値が指定されているので、2つ目のZ方向についての値を0にすることでZ方向には長さが変化しなくなる。

私は後者の方法をとることにしたが、果たしてそれで正しい結果の出る計算になっているのかは分からない。正しくはもっと長い箱を用意するべきなのだろう。

ここまで出来たら、以前書いた通り(というか引用させていただいた通り)READMEファイルをcshかshにして実行可能にし、プログラムを実行する。
この段階はまだアンブレラサンプリングではなく、その前の反応座標を決定する段階なので、ここからが(計算コスト的に)長い。

長くなったし続きはまたまとめることにする。
実際ここから先はまだやっている途中なので、どうしたら上手くいくのか分かっていないところもあるので。