仓库源文，站点原文

layout: post title: GROMACS如何做之膜模拟 categories:

科 tags:
gmx

原文链接
2016-09-19 21:11:59 翻译: 刘恒江; tcl代码: 刘胜堂; 补充整理: 李继存

运行膜模拟

GROMACS用户在运行双脂质层模拟时往往会遇到各种问题, 尤其是体系中还包含蛋白质的时候. 推荐模拟双脂质层和膜蛋白的用户参考教程KALP-15 in DPPC.

对膜蛋白体系完成一次模拟需要执行下列步骤:

选择一个能用于蛋白质和脂质分子的力场.
将蛋白质插入到膜结构中. (例如, 在预备好的双分子层上使用g_membed命令, 或进行一个粗粒化自组装模拟然后再转换回全原子表示)
向体系中加入溶剂. 并加入离子中和多余的电荷, 确保整个体系最终呈电中性, 同时调整最终的离子浓度
运行能量最小化.
让膜适应蛋白质. 典型的方法是, 在限制蛋白质所有重原子的情况下(力常数可取1000 kJ/mol-nm2)运行约5~10 ns的MD.
去掉限制进行平衡.
运行成品MD.

使用`genbox`命令添加水

当使用genbox命令(5.0以后版本中为gmx solvate)向预备好的双层膜体系中添加水时, 水分子往往会进入到膜的夹层中, 去除这些水分子的解决办法有以下几种:

将vdwradii.dat文件从你的$GMXLIB复制到当前的工作目录下, 调大脂质层分子中原子的范德华半径(建议调整碳原子的半径为0.35~0.5 nm)以阻止genbox命令将水分子加入夹层中;
运行一段短时间的MD, 让脂质层分子的憎水作用将水分子排出;
编辑结构文件, 手动删除双脂层内部的水分子(记住, 不要忘了修改gro文件中的原子数目, 以及top文件中的水分子数目); 或者
借助一些脚本来删除;

删除膜内部的水分子

awk脚本

此脚本可自动计算保留的水分子个数, 更新总原子数目.

<table class="highlighttable"><th colspan="2" style="text-align:left">rmWat.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">usage="\ >>>>>>>>>>>>>>>> rmWat <<<<<<<<<<<<<<<< >>>>>>>>>>>>>>>> Jicun LI <<<<<<<<<<<<<<<< >>>>>>>>>> 2016-09-19 11:24:51 <<<<<<<<< >> Usage: rmWat <File.gro> {Zmin} {Zmax} {Nsit} {Ncol} >> Default: rmWat N/A 2 6 3 0 >> Zmin: 水分子中原子z坐标的最小值 >> Zmax: 水分子中原子z坐标的最大值, 删除 Zmin<z<Zmax 的水分子原子 >> Nsit: 水模型的位点数; 3/4/5 >> Ncol: gro文件中的速度列数, 不含速度为0, 含速度为3"

[[ $# -lt 1 ]] && { echo "$usage"; exit; }

File=$1 Zmin=2; [[ $# -ge 2 ]] && { Zmin=$2; } Zmax=6; [[ $# -ge 3 ]] && { Zmax=$3; } Nsit=3; [[ $# -ge 4 ]] && { Nsit=$4; } Ncol=0 [[ $# -ge 5 ]] && { Ncol=$5; }

Fout=${File%.gro}~rmWat.gro

# 获取脂质分子的原子数目, 输出其坐标 Natm=$(awk ' BEGIN{Natm=0} $1!~/SOL/ && NF>4 { Natm++; print >"_Lip.gro"; next } END { print Natm} ' $File)

# 获取删除水后的总原子数, 输出水的坐标 Ntot=$(awk -v Ntot=$Natm -v Zmin=$Zmin -v Zmax=$Zmax \ -v Nsit=$Nsit -v Ncol=$Ncol ' $1~/SOL/ { Atom[1]=$0 YesInc=0 if($(NF-Ncol)<Zmin || $(NF-Ncol)>Zmax) YesInc=1 for(i=2; i<=Nsit; i++) { getline; Atom[i]=$0 if($(NF-Ncol)<Zmin || $(NF-Ncol)>Zmax) YesInc=1 } if(YesInc) { Ntot += Nsit for(i=1; i<=Nsit; i++) print Atom[i] >"_Wat.gro" } } END{print Ntot} ' $File)

# 组合文件 head -n 1 $File > $Fout echo $Ntot >>$Fout cat _Lip.gro _Wat.gro >>$Fout tail -n 1 $File >>$Fout

# 使用gmx工具重新设定分子, 原子编号 gmx editconf -f $Fout -o $Fout

# 删除临时文件 rm -rf _Lip.gro _Wat.gro

echo '保留水分子数目: ' $[($Ntot-$Natm)/$Nsit] </pre></div> </td></tr></table>

tcl脚本

此脚本在VMD中使用, 自动计算所需的z坐标极值.

<table class="highlighttable"><th colspan="2" style="text-align:left">rmwat.tcl</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">#remove water molecular near lipids within refered distance proc rmwat {lipres dz} { set lip [atomselect top "resname $lipres"]; set lipbox [measure minmax $lip]; set xmin [lindex [lindex $lipbox 0] 0]; set xmax [lindex [lindex $lipbox 1] 0]; set ymin [lindex [lindex $lipbox 0] 1]; set ymax [lindex [lindex $lipbox 1] 1]; set zmin [lindex [lindex $lipbox 0] 2]; set zmax [lindex [lindex $lipbox 1] 2]; set zmin [expr $zmin+$dz]; set zmax [expr $zmax-$dz]; set selA [atomselect top "same residue as water and z>$zmin and z<$zmax"] set selB [atomselect top "all not {same residue as water and z>$zmin and z<$zmax}"] puts "all not {same residue as water and z>$zmin and z<$zmax}" set rmwatlen [llength [$selA get resid]]; set rmwatnum [expr $rmwatlen/3];

<span style="color: #AA22FF">puts</span> <span style="color: #BB4444">&quot;remove water atom is $rmwatlen, and molecular is $rmwatnum&quot;</span>
<span style="color: #B8860B">$selB</span> <span style="color: #B8860B">writepdb</span> rmWat.pdb

<span style="color: #AA22FF; font-weight: bold">set</span> all <span style="color: #AA22FF; font-weight: bold">[</span><span style="color: #B8860B">atomselect</span> top water<span style="color: #AA22FF; font-weight: bold">]</span>
<span style="color: #AA22FF; font-weight: bold">set</span> allwat <span style="color: #AA22FF; font-weight: bold">[expr</span> <span style="color: #AA22FF; font-weight: bold">[</span><span style="color: #AA22FF">llength</span> <span style="color: #AA22FF; font-weight: bold">[</span><span style="color: #B8860B">$all</span> <span style="color: #B8860B">get</span> resid<span style="color: #AA22FF; font-weight: bold">]]</span><span style="color: #666666">/3</span><span style="color: #AA22FF; font-weight: bold">];</span>

<span style="color: #AA22FF; font-weight: bold">set</span> rmwatnum <span style="color: #AA22FF; font-weight: bold">[expr</span> <span style="color: #B8860B">$allwat-$rmwatnum</span><span style="color: #AA22FF; font-weight: bold">];</span>
<span style="color: #AA22FF">puts</span> <span style="color: #BB4444">&quot;Original water molecules number is $allwat&quot;</span>
<span style="color: #AA22FF">puts</span> <span style="color: #BB4444">&quot;Water molecules left $rmwatnum&quot;</span><span style="color: #AA22FF; font-weight: bold">;</span>

} </pre></div> </td></tr></table>

增大水相z方向的尺寸(此shell脚本效率较低, 只有存档意义, 不建议使用)

假定双层膜体系垂直于z轴延展, 其法向沿z轴, Chirs Neale提供了一个方法用以去除不需要的水分子. 具体操作如下(如果需要, 可以修改第41行使脚本适用于x/y/z方向延展的体系):

对初始双层膜构型initial.gro运行genbox命令, 得到solvated.gro;
cp solvate.gro new_waters.gro;
使用文本编辑器删除new_waters.gro中除第一步新添加的水之外的内容(若初始构型initial.gro中含有水分子, 这些水分子也要删去). 在vim中删除多行可以用指令ndd, n为行数;
修改脚本keepbyz.sh中的upperz和lowerz参数为需要的值, 并将sol参数改为溶剂分子的名称(upperz和lowerz分别为双层膜上下边界的z轴坐标, 注意脂质分子的头基周围需要保留一定数目的水分子, 所以请选择合适的值. sol对应溶剂分子的名称, 默认为SOL)
执行chmod +x keepbyz.sh确保脚本可执行, 然后运行./keepbyz.sh new_waters.gro > keep_these_waters.gro;
运行tail -1 initial.gro > last_line.gro获取盒子信息;
head -$(expr $(cat initial.gro | wc -l | awk '{print $1}') - 1 ) initial.gro > not_last_line.gro;
cat not_last_line.gro keep_these_waters.gro last_line.gro > new_system.gro; 手动修改new_system.gro中的原子数目, 使其与体系的总原子数目相等;
editconf -f new_system.gro -o new_system_sequential_numbers.gro; (5.0以后版本中使用gmx editconf); 修改拓扑文件, 确保与gro中的分子数一致.

keepbyz.sh脚本内容如下:

<table class="highlighttable"><th colspan="2" style="text-align:left">keepbyz.sh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">#!/bin/bash # give x.gro as first command line arguement upperz=6.417 lowerz=0.820 sol=SOL count=0 cat $1 | grep "$sol" | while read line; do for first in $line; do break done if [ "$count" = 3 ]; then count=0 fi count=$(expr $count + 1) if [ "$count" != 1 ]; then continue fi l=${#line} m=$(expr $l - 24) # would use -48 if velocities are also in .gro and -24 otherwise. 如果gro文件中包含速度信息, 使用48 i=1 for word in ${line:$m}; do if [ "$i" = 1 ]; then popex=$word else if [ "$i" = 2 ]; then popey=$word else if [ "$i" = 3 ]; then popez=$word break fi fi fi i=$(expr $i + 1) done nolx=`echo "$popez > $upperz" | bc` nohx=`echo "$popez < $lowerz" | bc` myno=$(expr $nolx + $nohx) if [ "$myno" != 0 ]; then z=${#first} if [ "$z" != 8 ]; then sfirst="[[:space:]]$first" else sfirst=$first fi `echo grep $sfirst $1` fi done </pre></div> </td></tr></table> 脚本中默认使用3位点的水分子模型, 若使用TIP4P水模型, 则应将

if [ "$count" = 3 ]; then
    count=0
fi

修改为:

if [ "$count" = 4 ]; then
    count=0
fi

同理, 使用TIP5P水模型时应修改为5.

为提高效率, 该脚本的作者还提供了一个同样功能的C程序. 所用水模型的位点数通过命令行参数指定, 但你仍然需要修改接近底部的源代码, 指定自己的选择标准, 并使用-24或-48来指定你的gro文件是否包含速度.

<table class="highlighttable"><th colspan="2" style="text-align:left">keepbyz.c</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%">#include <stdio.h> #include <stdlib.h> #include "string.h"

#define LINE_SIZE 1000

void showUsage(const char *c){ printf("Usage: %s <solvent.gro> <numAtomsInWaterMolecule>\n",c); printf(" (e.g. 3 for TIP3P, 4 for TIP4P)\n"); }

int main(int argn, char *args[]){ FILE *mf; char *c; char linein[LINE_SIZE]; float mx,my,mz; int count,keep; int atomInWater;

<span style="color: #AA22FF; font-weight: bold">if</span>(argn<span style="color: #666666">!=3</span>){
    showUsage(args[<span style="color: #666666">0</span>]);
    exit(<span style="color: #666666">1</span>);
}
<span style="color: #AA22FF; font-weight: bold">if</span>(sscanf(args[<span style="color: #666666">2</span>],<span style="color: #BB4444">&quot;%d&quot;</span>,<span style="color: #666666">&amp;</span>atomInWater)<span style="color: #666666">!=1</span>){
    printf(<span style="color: #BB4444">&quot;error: unable to determine number of atoms in the water model</span><span style="color: #BB6622; font-weight: bold">\n</span><span style="color: #BB4444">&quot;</span>);
    showUsage(args[<span style="color: #666666">0</span>]);
    exit(<span style="color: #666666">1</span>);
}

mf<span style="color: #666666">=</span>fopen(args[<span style="color: #666666">1</span>],<span style="color: #BB4444">&quot;r&quot;</span>);
<span style="color: #AA22FF; font-weight: bold">if</span>(mf<span style="color: #666666">==</span><span style="color: #AA22FF">NULL</span>){
    printf(<span style="color: #BB4444">&quot;error: unable to open the membrane file %s</span><span style="color: #BB6622; font-weight: bold">\n</span><span style="color: #BB4444">&quot;</span>,args[<span style="color: #666666">1</span>]);
    exit(<span style="color: #666666">1</span>);
}

count<span style="color: #666666">=0</span>;
keep<span style="color: #666666">=0</span>;
<span style="color: #AA22FF; font-weight: bold">while</span>(fgets(linein,LINE_SIZE,mf)<span style="color: #666666">!=</span><span style="color: #AA22FF">NULL</span>){
    <span style="color: #AA22FF; font-weight: bold">if</span>(count<span style="color: #666666">==</span>atomInWater){
        count<span style="color: #666666">=0</span>;
        keep<span style="color: #666666">=0</span>;
    }
    count<span style="color: #666666">++</span>;
    c<span style="color: #666666">=&amp;</span>linein[strlen(linein)<span style="color: #666666">-24</span>];           <span style="color: #008800; font-style: italic">// would use -48 if velocities are also in .gro and -24 otherwise</span>
    <span style="color: #AA22FF; font-weight: bold">if</span>(sscanf(c,<span style="color: #BB4444">&quot;%f %f %f&quot;</span>,<span style="color: #666666">&amp;</span>mx,<span style="color: #666666">&amp;</span>my,<span style="color: #666666">&amp;</span>mz)<span style="color: #666666">!=3</span>) <span style="color: #AA22FF; font-weight: bold">continue</span>;
    <span style="color: #AA22FF; font-weight: bold">if</span>(count<span style="color: #666666">==1</span>){
        <span style="color: #008800; font-style: italic">//Add your selection terms here to set keep=1 when you want to keep that particular solvent</span>
        <span style="color: #AA22FF; font-weight: bold">if</span>(mz<span style="color: #666666">&lt;7.0</span>){
            keep<span style="color: #666666">=1</span>;
        }
    }
    <span style="color: #AA22FF; font-weight: bold">if</span>(keep)printf(<span style="color: #BB4444">&quot;%s&quot;</span>,linein);
}
fclose(mf);