GROMACS如何做之膜模拟

类别:    标签: gmx   阅读次数:   版权: (CC) BY-NC-SA

运行膜模拟

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

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

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

使用genbox命令添加水

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

删除膜内部的水分子

awk脚本

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

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=$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]

tcl脚本

此脚本在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 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];

	puts "remove water atom is $rmwatlen, and molecular is $rmwatnum"
	$selB writepdb rmWat.pdb

	set all [atomselect top water]
	set allwat [expr [llength [$all get resid]]/3];

	set rmwatnum [expr $allwat-$rmwatnum];
	puts "Original water molecules number is $allwat"
	puts "Water molecules left $rmwatnum";
}

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

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

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

脚本中默认使用3位点的水分子模型, 若使用TIP4P水模型, 则应将

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

修改为:

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

同理, 使用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 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;

	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 model\n");
		showUsage(args[0]);
		exit(1);
	}

	mf=fopen(args[1],"r");
	if(mf==NULL){
		printf("error: unable to open the membrane file %s\n",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 otherwise
		if(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 solvent
			if(mz<7.0){
				keep=1;
			}
		}
		if(keep)printf("%s",linein);
	}
	fclose(mf);
}

外部链接

评论

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 几种颜色映射方案的解析式
后一篇: 自动更新谷歌host文件的简单程序

访问人次(2015年7月 9日起): | 最后更新: 2017-12-13 23:46:45 UTC | 版权所有 © 2008 - 2017 Jerkwin