Gromacs

Membrane Simulations

    Running Membrane Simulations

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

    1. Choose a force field for which you have parameters for the protein and lipids.
    2. Insert the protein into the membrane. (For instance, use g_membed on a pre-formed bilayer or do a coarse-grained self-assembly simulation and then convert back to the atomistic representation.)
    3. Solvate the system and add ions to neutralize excess charges and adjust the final ion concentration.
    4. Energy minimize.
    5. Let the membrane adjust to the protein. Typically run MD for ~5-10ns with restraints (1000 kJ/(mol nm2) on all protein heavy atoms.
    6. Equilibrate without restraints.
    7. Run production MD.

    Adding waters with genbox

    When 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

    • copy vdwradii.dat from your $GMXLIB location to the working directory, and edit it to increase the radii of your lipid atoms (between 0.35 and 0.5nm is suggested for carbon) to prevent genbox from seeing interstices large enough for water insertion,
    • a short MD run to get the hydrophobic effect to exclude these waters,
    • editing your structure by hand to delete them (remembering to adjust your atom count for .gro files and to account for any changes in the topology), or
    • use a script someone wrote to remove them.

    Increasing the size of the water segment in the z-direction

    Herein, 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 x/y/z dimensions (around line 41 of the script).

    • run genbox on initial.gro to create solvated.gro
    • cp solvated.gro new_waters.gro
    • use vi to remove everything in new_waters.gro except the new waters (make sure you remove waters that were in initial.gro)
    • use vi to edit keepbyz.sh setting upperz and lowerz variables as you please (but remember you want to keep waters around the head groups of your lipid, so choose wisely) and sol to the name of your solvent molecule
    • ./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 new_waters.gro last_line.gro > new_system.gro
    • editconf -f new_system.gro -o new_system_sequential_numbers.gro

    Here is the script keepbyz.sh (make sure you chmod +x keepbyz.sh):

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

    Page last modified 12:22, 19 Sep 2010 by orbeckst