Molecular Surfaces Using IsoSurfaces

[sargis]

There are a number of programmes for computing Molecular Surfaces; MSMS being one of the best that I use regularly in MGLTools. However, for PyRx, I needed something with permissive free software licence, and since, as far as I know, there is none available, I decided to write one myself. I remember reading A General Mesh Generation Tool for Finite Element Analysis of Bio-Structures; so I just needed to code the first formula from this article, with blobbyness parameter, and use IsoSurface Module from MayaVi, which provides convenient access to appropriate functionality in VTK. Here is the complete code:

Code: Select all

    def DisplayMolecularSurface(self, molecule):
        bloby = -1
        spacing = 10       
        bounds = molecule.assembly.GetBounds()
        dimX = int(abs(bounds[1]-bounds[0])) + spacing
        dimY = int(abs(bounds[3]-bounds[2])) + spacing
        dimZ = int(abs(bounds[5]-bounds[4])) + spacing
        origX =   int(bounds[0]) - spacing/2
        origY =   int(bounds[2]) - spacing/2
        origZ =   int(bounds[4]) - spacing/2               
        img = tvtk.tvtk.tvtk.ImageData(origin=[origX,origY,origZ], spacing=[1,1,1], dimensions=[dimX, dimY, dimZ] )
        scalars = zeros((dimX, dimY, dimZ))
        for atom in molecule.allAtoms:
            coords = atom._coords[0]
            intX = int(coords[0])
            intY = int(coords[1])
            intZ = int(coords[2])
            iX = int(abs(origX - coords[0]))
            iY = int(abs(origY - coords[1]))
            iZ = int(abs(origZ - coords[2]))
            p = 2
            for i in range(-p, p+1):               
                for j in range(-p, p+1):
                    for k in range(-p, p+1):                                   
                        e = (intX + i- coords[0])*(intX + i - coords[0]) +\
                               (intY + j - coords[1])*(intY + j - coords[1])+\
                               (intZ +k - coords[2])*(intZ +k - coords[2])                   
                        scalars[iX+i, iY+j, iZ+k] += math.exp(bloby*(e/(atom.vdwRadius*atom.vdwRadius)-1))
        s = scalars.transpose().copy()
        img.point_data.scalars = ravel(s)
        img.point_data.scalars.name = molecule.name
        src = VTKDataSource(data=img, name=molecule.name+"-grid")
        self.engine.add_source(src)
        mayavi = self.engine
        iso = IsoSurface()
        mayavi.add_module(iso)


I committed this code to CVS under PyRxDIST/PyRx/mayaviEngine.py.