博文
GROMACS如何做之膜模拟
|
2016-09-19 21:11:59 翻译: 刘恒江; tcl代码: 刘胜堂; 补充整理: 李继存
GROMACS用户在运行双脂质层模拟时往往会遇到各种问题, 尤其是体系中还包含蛋白质的时候. 推荐模拟双脂质层和膜蛋白的用户参考教程KALP-15 in DPPC.
对膜蛋白体系完成一次模拟需要执行下列步骤:
选择一个能用于蛋白质和脂质分子的力场.
将蛋白质插入到膜结构中. (例如, 在预备好的双分子层上使用g_membed命令, 或进行一个粗粒化自组装模拟然后再转换回全原子表示)
向体系中加入溶剂. 并加入离子中和多余的电荷, 确保整个体系最终呈电中性, 同时调整最终的离子浓度
运行能量最小化.
让膜适应蛋白质. 典型的方法是, 在限制蛋白质所有重原子的情况下(力常数可取1000 kJ/mol-nm2)运行约5~10 ns的MD.
去掉限制进行平衡.
运行成品MD.
当使用genbox命令(5.0以后版本中为gmx solvate)向预备好的双层膜体系中添加水时, 水分子往往会进入到膜的夹层中, 去除这些水分子的解决办法有以下几种:
将vdwradii.dat文件从你的$GMXLIB复制到当前的工作目录下, 调大脂质层分子中原子的范德华半径(建议调整碳原子的半径为0.35~0.5 nm)以阻止genbox命令将水分子加入夹层中;
运行一段短时间的MD, 让脂质层分子的憎水作用将水分子排出;
编辑结构文件, 手动删除双脂层内部的水分子(记住, 不要忘了修改gro文件中的原子数目, 以及top文件中的水分子数目); 或者
借助一些脚本来删除;
此脚本可自动计算保留的水分子个数, 更新总原子数目.
| rmWat.bsh | |
|---|---|
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 | 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=$1Zmin=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]=$0YesInc=0
if($(NF-Ncol)<Zmin ||$(NF-Ncol)>Zmax)YesInc=1
for(i=2; i<=Nsit; i++){
getline; Atom[i]=$0if($(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] |
此脚本在VMD中使用, 自动计算所需的z坐标极值.
| rmwat.tcl | |
|---|---|
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 | #remove water molecular near lipids within refered distanceproc rmwat {lipres dz}{set lip [atomselect top "resname $lipres"];set lipbox [measure minmax $lip];set xmin [lindex[lindex$lipbox0]0];set xmax [lindex[lindex$lipbox1]0];set ymin [lindex[lindex$lipbox0]1];set ymax [lindex[lindex$lipbox1]1];set zmin [lindex[lindex$lipbox0]2];set zmax [lindex[lindex$lipbox1]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[$selAget resid]];set rmwatnum [expr$rmwatlen/3];puts"remove water atom is $rmwatlen, and molecular is $rmwatnum"$selBwritepdb rmWat.pdb
set all [atomselect top water]set allwat [expr[llength[$allget resid]]/3];set rmwatnum [expr$allwat-$rmwatnum];puts"Original water molecules number is $allwat"puts"Water molecules left $rmwatnum";} |
假定双层膜体系垂直于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脚本内容如下:
| keepbyz.sh | |
|---|---|
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 | #!/bin/bash# give x.gro as first command line arguementupperz=6.417
lowerz=0.820
sol=SOL
count=0
cat $1 | grep "$sol" | whileread line; dofor first in $line; dobreakdoneif["$count"=3]; thencount=0
ficount=$(expr $count + 1)if["$count" !=1]; thencontinuefil=${#line}m=$(expr $l - 24)# would use -48 if velocities are also in .gro and -24 otherwise. 如果gro文件中包含速度信息, 使用48i=1
for word in ${line:$m}; doif["$i"=1]; thenpopex=$word
elseif["$i"=2]; thenpopey=$word
elseif["$i"=3]; thenpopez=$word
breakfififii=$(expr $i + 1)donenolx=`echo"$popez > $upperz" | bc`nohx=`echo"$popez < $lowerz" | bc`myno=$(expr $nolx + $nohx)if["$myno" !=0]; thenz=${#first}if["$z" !=8]; thensfirst="[[:space:]]$first"elsesfirst=$first
fi`echo grep $sfirst $1`fidone |
脚本中默认使用3位点的水分子模型, 若使用TIP4P水模型, 则应将
修改为:
同理, 使用TIP5P水模型时应修改为5.
为提高效率, 该脚本的作者还提供了一个同样功能的C程序. 所用水模型的位点数通过命令行参数指定, 但你仍然需要修改接近底部的源代码, 指定自己的选择标准, 并使用-24或-48来指定你的gro文件是否包含速度.
| keepbyz.c | |
|---|---|
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 | #include<stdio.h>#include<stdlib.h>#include"string.h"#define LINE_SIZE 1000voidshowUsage(constchar*c){
printf("Usage: %s <solvent.gro> <numAtomsInWaterMolecule>n",c);
printf(" (e.g. 3 for TIP3P, 4 for TIP4P)n");
}
intmain(int argn, char*args[]){
FILE*mf;
char*c;
char linein[LINE_SIZE];
float mx,my,mz;
int count,keep;
int atomInWater;
if(argn!=3){
showUsage(args[0]);
exit(1);
}
if(sscanf(args[2],"%d",&atomInWater)!=1){
printf("error: unable to determine number of atoms in the water modeln");
showUsage(args[0]);
exit(1);
}
mf=fopen(args[1],"r");
if(mf==NULL){
printf("error: unable to open the membrane file %sn",args[1]);
exit(1);
}
count=0;
keep=0;
while(fgets(linein,LINE_SIZE,mf)!=NULL){
if(count==atomInWater){
count=0;
keep=0;
}
count++;
c=&linein[strlen(linein)-24]; // would use -48 if velocities are also in .gro and -24 otherwiseif(sscanf(c,"%f %f %f",&mx,&my,&mz)!=3) continue;
if(count==1){
//Add your selection terms here to set keep=1 when you want to keep that particular solventif(mz<7.0){
keep=1;
}
}
if(keep)printf("%s",linein);
}
fclose(mf);
}
|
膜蛋白模拟(Phil Biggin, Bert de Groot) 幻灯片, Biggin的视频 de Groot的视频
GROMACS膜蛋白模拟教程, 用于演示当模拟嵌入脂质双层中的蛋白质时会遇到那些问题.
组合OPLS-AA力场和Berger脂质, 目的, 方法和测试的详细说明.
膜蛋白与不同力场GAFF, CHARMM BERGER的一些拓扑 Shirley W. I. Siu, Robert Vacha, Pavel Jungwirth, Rainer A. Böckmann: Biomolecular simulations of membranes: Physical properties from different force fields.
Lipidbook提供了膜以及膜蛋白模拟中用到的脂质, 活性剂以及其他分子的力场参数. 相关论文: J. Domański, P. Stansfeld, M.S.P. Sansom, and O. Beckstein. J. Membrane Biol. 236 (2010), 255-258. doi:10.1007/s00232-010-9296-8.
◆本文地址: http://jerkwin.github.io/2016/09/19/GROMACS如何做之膜模拟/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-1003931.html
上一篇：Martini粗粒化力场使用手册：2 常见问题
下一篇：自动更新谷歌host文件的简单程序
