うみうしの日誌

~(・ ω ・)~nullnull

Gromacsをとりあえず使ってみる

理論をやっている人がどういうモチベーションで分子動力学シミュレーションをするのかちょっと見当がつかないが、実験科学者からするとまずとりあえず使ってみてから考えるというのが1つのアプローチだと思う。私はそうだった。
しかしGromacsはバイナリの配布があるわけでもないので、「とりあえず使ってみる」の段階で既にかなりハードルが高かった。そういうわけで色々試してとりあえず動かしてみることのできた簡単な方法をここで紹介したいと思う。

これで実践できる計算が数学的に大丈夫なのか(近似とか)は私自身は確認していないが、とりあえず使うことがモチベーションなのでそのあたりは後日気にすることにする。

0.Gromacsをインストールする

基本的にLinux上で動かすことを想定しているようだった。しかし、繰り返しになるがとりあえず使うことがモチベーションなので、(人によるだろうが)Windows上で何とか動かしたかった。以下のサイトが非常に参考になった
フリーな分子動力学計算プログラムの一つである GROMACS と連携ソフトウエアを用いて、一般分子(非生体高分子)の1分子モデリングから分子集合体の計算をWindowsノートPCを用いて行う方法を説明する。

『GROMACSと連携ソフトウエアによる分子動力学計算 updated』~国立研究開発法人産業技術総合総合研究所~

というかここを完全に参考にすれば、触って動かすことは可能である。私の場合はC:/Gow/local/gromacs/share/gromacs以下に計算に使うファイルを入れて管理していた。これならGromacsをとりあえずは使わなくなったとき、Gow以下を削除することで関連するファイルをまとめて削除できる。

ただ、ネット上に落ちているGromacsのチュートリアルでは時折Linuxの機能を要求してくる。例えばシミュレーションの結果をグラフにするのにGraceというソフトが便利に使えるらしいが、調べてみるとLinux上で動くソフトだった。
また、LinuxWindowsでは標準出力やエラーをnullに吐き出す記述の仕方などが少し異なるらしく(1やら2やら順番がややこしくてあまり分かっていないが、少なくともLinuxではnullでWindowsではnulである)私がネットに落ちているプログラムを使おうとすると、しばしばエラーが出て動かなかった。

この辺のことを頑張りだすととりあえず使うって感じでなくなるし、多分Linux入れたほうが早いので、操作に慣れてきたらLinuxを入れてその上でGromacsを動かしてみよう(丸投げ)


Linuxでアレコレする場合の参考例】
『CentOS 7 ダウンロード』
『Gromacsのインストールとパフォーマンスの検証 on Docker and CentOS 7』

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

上の「GROMACSと連携ソフトウエアによる分子動力学計算 updated」を完全に参考にした皆さんは、既に4,4'-ビピリジンの水中におけるシミュレーションが出来たと思われる(してなくても良い。私はやってみていない)。

さて、実際に上の手順を実行していなくても流れを一読してもらうとわかると思うが、シミュレーションに必要なファイルを得るまでの手順が案外長い。また、タンパク質について計算する場合はもうちょっと色々平衡化する必要があるらしい(私のモチベーションはペプチド・タンパク質と脂質膜の相互作用の計算だった)。
その辺りを全部気を付けながらMD計算にたどり着くのはしんどかったため、Web上で何か上手くやれる方法がないか探っていたところ、おあつらえ向きのサイトがあった。

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

GUIでシミュレーションの系を作ってくれ、同時に計算に必要なファイルを作ってくれる。Gromacsにも対応しており、ページ左のカラムの「Input Generator」から種々の系に対応したファイルを作製できる。選択できる力場が、多くの場合CHARMm力場になるようだが、気にしなくてもよいだろう。
あまり詳しく覚えていないが、CHARMM24?(多分1世代前。今は36?)ではα-helixを過剰に安定に評価するんだったか、何かあった気がする。AMBER力場やGROMOS力場など色んな力場で、パラメータを変えながら脂質膜中でのα-helixかβ-helix構造の計算をさせ、どの条件で実際のNMRとよく似た結果が得られるかシミュレーション実験した論文があった気がする。気がするが、忘れたし手元に印刷していなかったので何に書いてあったか分からなくなった。
そういうところまで気になる日が来たら自分でちゃんと調べて系をたてればよいのだ。

で、CHARMM-GUIではサイトの指示に従って入力していくだけで最終的に自分の走らせたいシミュレーションに必要なファイルを全て用意してくれる。全てのInput Generatorを試したことはなく、「Membrane Builder」と「Martini Maker」しか使ったことはないが、多分大体の場合最後のページの「download.tgz」をクリックしてファイルをダウンロードすればよい。

このCHARMM-GUIの画面も最初は何が何を意味するかは分かりにくいと思うが、使いながら何となく察していってほしい(丸投げ)


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

大学の研究室のWebページで、「この記事は主に自分と研究室の学生向けのメモです.」とあったので直リンクは止めた。
余談だが自分は直リンクを悪とする文化で育ったため直リンクに感情的な抵抗があるが、その辺り正確な引用をしたい気分と衝突するので、こういうときに困る

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

ここまで来たらあとはシミュレーションを走らせるだけだ。
ダウンロードしたファイルは多分charmm-gui.tgzだと思われ、これを解凍するとcharmm-guiというフォルダが得られる。Windowsだと多分デフォで.tgzを解凍できないので、適当なフリーソフトを入れるなどして解凍する。

さてcharmm-guiフォルダの中に入ると色々なファイルがあるが、とりあえず無視してその中のgromacsフォルダに入る。gromacsフォルダの下にはおそらく

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

があるのではないかと思う。

もしLinuxを使っているのであれば、READMEファイルがシミュレーションを走らせるためのシェルスクリプトとなっているので、

chmod +x README
mv README run.csh
./run.csh

とすれば無事シミュレーションが走るはず。完全に下のblogの方を参考にさせていただいた。

gromacsのフォルダ内にある、READMEを見てみると
(中略)
cshスクリプトになっています。

$ chmod +x README
$ mv README run.csh
としてやれば、gromacsがインストールされている環境であれば、

$ ./run.csh
でシミュレーションが走るはずです。

『二重膜の全原子MDシミュレーションを行う。』~tsunekohのblog~


Windowsの場合、ちょっとCygwinなどの環境がなければシェルスクリプトの実行は出来ない。Gowで出来るような気もするが私には出来なかったので、READMEの中身を書き下して直接コマンドプロンプトに書き込んだ。
READMEをメモ帳なりで開いてもらえば分かるが、例えばMartini MakerのBilayer Makerの結果なら以下のようになっており、

#!/bin/sh
gmx grompp -f step6.0_equilibration.mdp -o step6.0_equilibration.tpr -c step5_assembly.box.pdb -p system.top -n index.ndx
gmx mdrun -deffnm step6.0_equilibration

for i in 1 2 3 4 5 6
do
    let j=$i-1        
    gmx grompp -f step6.${i}_equilibration.mdp -o step6.${i}_equilibration.tpr -c step6.${j}_equilibration.gro -p system.top -n index.ndx
    gmx mdrun -deffnm step6.${i}_equilibration
done

gmx grompp -f step7_production.mdp -o step7_production.tpr -c step6.6_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step7_production

単に数字を変えて各mdpファイルにより平衡化を進め、最後にstep7_production.mdpでMD計算をしているだけである。今見たらcshではなくshだった。Linuxでも標準シェルがbashだったりすれば./run.cshは動かなかった気がする。Linux詳しくないし忘れたけど。その辺も適当に試してみてほしい

で、

gmx grompp -f step6.0_equilibration.mdp -o step6.0_equilibration.tpr -c step5_assembly.box.pdb -p system.top -n index.ndx
gmx mdrun -deffnm step6.0_equilibration -v
gmx grompp -f step6.1_equilibration.mdp -o step6.1_equilibration.tpr -c step6.0_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.1_equilibration -v
gmx grompp -f step6.2_equilibration.mdp -o step6.2_equilibration.tpr -c step6.1_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.2_equilibration -v
gmx grompp -f step6.3_equilibration.mdp -o step6.3_equilibration.tpr -c step6.2_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.3_equilibration -v
gmx grompp -f step6.4_equilibration.mdp -o step6.4_equilibration.tpr -c step6.3_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.4_equilibration -v
gmx grompp -f step6.5_equilibration.mdp -o step6.5_equilibration.tpr -c step6.4_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.5_equilibration -v
gmx grompp -f step6.6_equilibration.mdp -o step6.6_equilibration.tpr -c step6.5_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step6.6_equilibration -v
gmx grompp -f step7_production.mdp -o step7_production.tpr -c step6.6_equilibration.gro -p system.top -n index.ndx
gmx mdrun -deffnm step7_production -v
gmx trjconv -s step7_production.tpr -f step7_production.xtc -o step7-mod.xtc -pbc mol
0

私は単純に書き下した上の内容を適当なテキストファイルに保存しておき、全部マウスでコピーしコマンドプロンプトに張り付けていた。そうすると改行ごとにコマンドが実行され、./run.cshと実質同じことが出来る。




少しだけ改変したので補足する。
-v:これをつけると、計算終了予定時刻を表示してくれる。シミュレーションは系の大きさにもよるがかなり時間がかかるため、これがないとどこまで計算が進んでいるのか判断しづらく、地味に困る(ログが定期的に書き出されるはずなのでそれを見れば判断できる)

gmx trjconv -s step7_production.tpr -f step7_production.xtc -o step7-mod.xtc -pbc mol
0

これは『GROMACSと連携ソフトウエアによる分子動力学計算 updated』の内容からもらってきた。
gmx trjconv~pbc molの操作を、どの分子に対して行うかを「0」で選んでいる。一応対話型で選択でき、0は全原子に対して操作を行う。確か1でProteinだった。あとは忘れた。

3.その他

これでとりあえず走らせることはできたはず。上記のコードはGromacsのVersion 5以降の話で、Version 4までの場合は"gmx hogehoge"の形式ではなく"g_hogehoge"の形のコマンドとなる。『GROMACSと連携ソフトウエアによる分子動力学計算 updated』からGromacsをとってきている場合はGromacs5.0.7なので、上記の通り実行できる。

ただ、CHARMM-GUIで提供されるmdpファイルの形式がどうやらversion 4の頃のものに対応しているらしく(CHARMM-GUIのサイトにはVersion 5に対応したと書いてあった気がするが、定かではない)、コマンドを実行していると途中でいくつか注意が表示される(この計算方法は古いよ!とか、このパラメータ名古いバージョンのだから書き換えるで!とか)
おそらく実行する分には問題ないはずなので気にせず推し進めればよい。ダメな場合は何がダメか結構丁寧に表示してくれたはずなので、それに従ってmdpその他のファイルを訂正して計算しなおすといい。

また脂質膜にタンパク質やペプチドを撒いた形でシミュレーションを走らせると、たまに平衡化の段階で構造が変化しすぎてしまい、エラーとなって計算が止まることがある。この場合、CHARMM-GUI曰く、平衡化を1fsにしても良いらしい(ただしそれはレアなケースだとも言っている)。あるいは、それが正しい計算をもたらすかは多分慎重に検討するべきだが、itpファイルの[position_restraints]をいじるのも1つの手段だと思われる。

シミュレーション的には多分計算後の解析の方が重要だと思われるが、今回はとりあえず動かすことだけ目指してきたので、データの解析などは後日またまとめることにする。


【[position_restraints]変更の際の参考】
位置拘束というのは非常に重要で有用なもので、特定の原子に任意の力の重みを付けて、動きにくくすることができます。

上述のエネルギー最小化の操作および平衡化の操作で、はじめに蛋白質の特に重原子(水素原子以外)に重みをつけながらシミュレーションを行うと、水溶媒や水素原子だけが動くことになり、水素結合を効率よく形成させることができるようになります。この目的のために位置拘束がとても使えます。

『Gromacsを使ったエネルギー最小化』~Ag++~

この方のblogには日頃より大変お世話になっております