GROMACS如何做之膜模拟

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]

#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";
}

#!/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

#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);
}