はじめに

Martiniを使って脂質二重膜のシミュレーションを行うチュートリアルを（途中まで）実施する。
環境としては

GROMACS ver5以上
VMD

がインストール済みでパスも通っていること。
また、GROMACSの使用についての知識があることが望ましい。

脂質2重膜の形成

bilayer-lipidome-tutorial-GMX5_2017Aug04.tgzを落とす。

$ tar -xzvf bilayer-lipidome-tutorial-GMX5_2017Aug04.tgz
$ cd  bilayer-lipidome-tutorial/

DPPCの二重膜を自己組織化するところから始める。
まずは128のDPPC分子をボックス内にランダムに配置する。

$ cd spontaneous-assembly/initial_assembly
$ gmx insert-molecules -ci DPPC-em.gro -box 7.5 7.5 7.5 -nmol 128 -radius 0.21 -try 500 -o 128_noW.gro

radiusでCG粒子の径を指定している。noWとは水無しの意味である。

近接による力を緩和するため、この状態でエネルギー最小化を行う。
DPPCのトポロジーファイルを下記リンクから入手する。

DPPCのitpの内容をコピペしてファイル名"martini_v2.0_DPPC_01.itp"に保存する。

DPPCの各CG粒子を定義するitpファイルも必要である。
下記リンクから"martini_v2.1.itp"を保存する。

エネルギー最小化計算のtprファイルを作成してから、MD計算を実行する。

$ gmx grompp -f minimization.mdp -c 128_noW.gro -p dppc.top -o dppc-min-init.tpr
$ gmx mdrun -deffnm dppc-min-init -v -c 128_noW.gro

gmx solvateで脂質1分子あたりCG水粒子を6つを追加する。
（CG水粒子は水分子を4つまとめたものであるので24分子分に相当する）

$ gmx solvate -cp 128_minimized.gro -cs water.gro -o waterbox.gro -maxsol 768 -radius 0.21

waterbox.groが水和された構造である。
ここから再びエネルギー最小化計算を行う。
水粒子が768粒子追加されたのでdppc.topを編集する。
”W 768”の前にあるセミコロンを削除し、コメントアウトを解除する。
つまり、下記のようになる。

[ molecules ]
DPPC 128
W 768

$ gmx grompp -f minimization.mdp -c waterbox.gro -p dppc.top -o dppc-min-solvent.tpr 
$ gmx mdrun -deffnm dppc-min-solvent -v -c minimized.gro

いよいよ自己組織化MDを行う。
約30nsの計算で十分である。

$ gmx grompp -f martini_md.mdp -c minimized.gro -p dppc.top -o dppc-md.tpr
$ gmx mdrun -deffnm dppc-md -v

手元の計算機では4コアを使って約3分で完了した。

               Core t (s)   Wall t (s)        (%)
       Time:     1227.409      153.426      800.0
                 (ns/day)    (hour/ns)
Performance:    16894.135        0.001

となっており、1日で16マイクロ秒の計算ができるらしい。全原子と比べると笑ってしまうくらい早い。

結果を可視化しよう。

$ gmx view -f dppc-md.xtc -s dppc-md.tpr

これはGROMACSのツールとして用意されているビューアーだが、私は使ったことがない。viewerを含む形でコンパイルもしていない。

pymolを使う場合は、trjconvで結合情報を含んだPDBファイルを書き出すことで可視化が可能である。

$ gmx trjconv -s dppc-md.tpr -f dppc-md.xtc -o dppc-md.pdb -pbc whole -conect
$ pymol dppc-md.pdb

VMDを使う場合は、そのためのスクリプトcg_bonds.tcl が既にディレクトリに用意されている。

$ chmod +x do_vmd.sh
$ ./do_vmd.sh

を行うとVMDが起動する。do_vmd.shの中身としてはpbc=wholeでトラジェクトリを変換し可視化ファイルviz.vmdを読み込むということをやっている。

適当な方向に複製すると脂質二重膜っぽくなる。

このシミュレーション、二重膜を形成するように少し水を少なめにしているそうで、gmx solvate の-maxsolオプションで水分子を増量してみるのもよいかもということ。

二重膜の解析

さて、次のステップのために二重膜がz軸に沿って形成されているか（つまりxy平面に広がっているか）をチェックする。
そうなっていなければgmx editconfで回転する。

私はx軸に沿って形成されていたのでy軸方向に90度回転。

gmx editconf -f dppc-md.gro -rotate 0 90 0 -o dppc-md-rotate.gro

先ほどはxyzで等方的な圧力カップリングを用いた。現在は膜を形成したので圧力は等方的でない。そのため、xy方向とz方向で圧力カップリングを分離してMD計算を行い、その後温度を341Kに上げて30nsの計算を行う。

$ mkdir dspc-eq
$ cp dppc-md.gro (私の場合はdppc-md-rotate.gro) dspc-eq/.
$ cp martini_md.mdp dspc-eq/.
$ cp martini*.itp dspc-eq/.
$ cp dppc.top dspc-eq/.
$ cd dspc-eq

後はシミュレーションを実行する。
"martini_md.mdp"を下記のように変更して圧力カップリングの設定を変更する。

Pcoupltype               = semiisotropic  ; semiisotropic
tau_p                    = 6.0       ; 12.0  ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility          = 3e-4   3e-4    ; 3e-4
ref_p                    = 1.0    1.0    ; 1.0

計算を実施する。

gmx grompp -f martini_md.mdp -c dppc-md-rotate.gro -p dppc.top -o dppc-md-semi.tpr
gmx mdrun -deffnm dppc-md-semi

計算が終わったら続いて温度を上げて30nsのMD計算を行う。
mdpファイルは下記のように変更し、martini_md_ext.mdpに保存した。

ref_t                    = 341 341
Pcoupl                   = parrinello-rahman
Pcoupltype               = semiisotropic  ; semiisotropic
tau_p                    = 12.0         ; 12.0  ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility          = 3e-4   3e-4    ; 3e-4
ref_p                    = 1.0    1.0    ; 1.0

gen_vel                  = no

計算を実行する。

gmx grompp -f martini_md_ext.mdp -c dppc-md-semi.gro -p dppc.top -o md.tpr
gmx mdrun -deffnm md

脂質1分子あたりの面積

脂質1分子あたりの面積を求める。
gmx energyでBox-XとBox-Yを選択してXY平面の面積を求める。

$ gmx energy -f md.edr -o box-xy.xvg

     > 10 11 0 [press Enter]

私のバージョンでは11と12だった。
さて面積を求めて脂質分子の数で割りる。
Tutorial中ではxmgraceを使うやり方が書いてあったが、エクセルでもpythonでも好きな方法でやれば良いと思う。

だいたい0.65nm2くらい。

二重膜の厚み

$ gmx make_ndx -f md.gro
    > a P*
    > q

$ gmx density -f md.xtc -s md.tpr -b 15000 -n index.ndx -o p-density.xvg -xvg no
    > 4

$ xmgrace p-density.xvg

ボックスをまたいでいるので変な表示になっている。

ここで終わります。
この後もいろいろやられていまして、物性値として水平方向への拡散や配向の計算が行われていました。
また、脂質のヘッドグループを変えて計算を行い、脂質の種類を変えると面積や厚みがどう変わるのかも検討しているようです。

VMDでの可視化について

他の計算でVMDで可視化する際はvmdでgroファイル等を開いた後、コンソールで

source ./cg_bonds.tcl
cg_bonds -top dppc.top

などしてcg_bonds.tclを読み込んだ後に、topファイルを読み込む。

すると、結合を示すことができる。