ext_water.py

 #!/usr/bin/python
import sys, math, re

class MolDat:
    def __init__(self):
        self.molnum   = -1
        self.molname  = ""
        self.atomnum  = -1
        self.atomname = ""
        self.coord    = []
        
def ReadLineGro(line):
    ret          = MolDat()
    ret.molnum   = int(line[0:5])
    ret.molname  = line[5 :10].strip()
    ret.atomname = line[10:15].strip()
    ret.atomnum  = int(line[15:20])
    ret.coord    = [float(line[20:28]), float(line[28:36]), float(line[36:44]) ] 
    return ret

def ReadGroFileAll(rfile, Waters, Lipids):
    for line in rfile[2 : len(rf_lines) - 1]:
        line_dat = ReadLineGro(line)
        if(line_dat.molname != "SOL"):
            Lipids.append(line_dat)
        else:
            Waters.append(line_dat)
    
    if(len(Waters) == 0):
        print "There are no solvent molecules in this system! "
        print "Please check the solvent molecule's name = SOL."
        quit()

def SearchBilayerReg(Lipids, boxlengs):
    ymax = 0.0
    zmax = 0.0
    zmin = boxlengs[2] * 2.0
    for lipid in Lipids:
        pos = lipid.coord
        if(pos[1] > ymax):
            ymax = pos[1]
        if(pos[2] > zmax): 
            zmax = pos[2]
        if(pos[2] < zmin):
            zmin = pos[2]
    return (ymax, zmax, zmin)
    
def RemoveWatersInBilayer(Waters, bound):
    RemWaterId = []
    NewWaters  = []
    for water in Waters:
        pos = water.coord
        if( (pos[1] < bound[0]) and (pos[2] < bound[1]) and (pos[2] > bound[2]) ):
            RemWaterId.append(water.molnum)
    RemWaterId = list(set(RemWaterId))
    for water in Waters:
        molid = water.molnum
        if not(molid in RemWaterId) :
            NewWaters.append(water)
    return NewWaters
            
def ReNumbering(Waters, Lipids):
    mol_n  = Lipids[-1].molnum
    atom_n = Lipids[-1].atomnum
    cur_water_id = -1
    for i in range(len(Waters)):
        w_id = Waters[i].molnum
        if(w_id != cur_water_id):
            cur_water_id = w_id
            mol_n += 1
        atom_n += 1
        Waters[i].molnum  = mol_n
        Waters[i].atomnum = atom_n
    
def Dump2GroFile(wf, Mols):
    for mol in Mols:
        wf.write("%5d%-5s%5s%5d%8.3f%8.3f%8.3f\n" % (mol.molnum, mol.molname, mol.atomname, mol.atomnum, mol.coord[0], mol.coord[1], mol.coord[2]) )
        
if __name__ == "__main__":
    argvs = sys.argv
    if (len(argvs) != 3):
        print "Usage $ python %s input_filename output_filename" % argvs[0]
        quit()
        print "Extract water molecules embedded in bilayer membrane."
        print "Edge line should be placed along X axis.             "

    rf        = open(argvs[1], "r")
    rf_lines  = rf.readlines()
    comment   = rf_lines[0]
    boxlengs  = [float(i) for i in rf_lines[-1].split()]

    Waters    = []
    Lipids    = []
    
    ReadGroFileAll(rf_lines, Waters, Lipids)
    bound     = SearchBilayerReg(Lipids, boxlengs)
    NewWaters = RemoveWatersInBilayer(Waters, bound)
    ReNumbering(NewWaters, Lipids)

    wf        = open(argvs[2], "w")
    wf.write(comment)
    alln = str(len(NewWaters) + len(Lipids))
    wf.write(alln + '\n')
    Dump2GroFile(wf, Lipids)
    Dump2GroFile(wf, NewWaters)
    wf.write("%10.5f%10.5f%10.5f\n"%(boxlengs[0], boxlengs[1], boxlengs[2]))
    
    rf.close()
    wf.close()

Remove water molecules embedded in bilayer membrane.
This script should be used after solvation like this:

gmx solvate -cp bilayer.gro -cs spc216.gro -o bilayerwater.gro -p topol.top