 |
|
|
|
Membrane SimulationsRunning Membrane SimulationsUsers frequently encounter problems when running simulations of lipid bilayers, especially when a protein is involved. Users seeking to simulate membrane proteins may find this tutorial useful. One protocol for the simulation of membrane proteins consists of the following steps:
Adding waters with genboxWhen generating waters around a pre-formed lipid membrane with genbox you may find that water molecules get introduced into interstices in the membrane. There are several approaches to removing these, including
Increasing the size of the water segment in the z-directionHerein, the z-direction is assumed to be the one normal to the membrane (bi)layer. This method was suggested by Chris Neale here. The script could be modified to get rid of new waters in any
Here is the script #!/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
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
That script assumes that you use a 3 atom water molecule. If you use TIP4P then you would want if [ "$count" = 3 ]; then
count=0
fi
to be changed to: if [ "$count" = 4 ]; then
count=0
fi
and analogously for TIP5P. Alternatively, use this simple C program, written by the same author, that was created because the bash script takes so long to complete. The number of atoms in your water model is specified on the command line, but you still need to edit the source near the bottom and craft your particular inclusion specifications and also to choose between -24 and -48 based on whether your .gro has velocities or not #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);
}
External material
|