|
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
|