Package rdkit :: Package Chem :: Package Subshape :: Module SubshapeBuilder

[hide private]
[frames] | no frames]

# Source Code for Module rdkit.Chem.Subshape.SubshapeBuilder

      1  # $Id$
      2  #
      3  # Copyright (C) 2007 by Greg Landrum
      4  #  All rights reserved
      5  #
      6  from rdkit import Chem,Geometry
      7  from rdkit.Chem import AllChem
      8  from rdkit.Chem.Subshape import SubshapeObjects
      9  from rdkit.Chem.Subshape import BuilderUtils
     10  from rdkit.six.moves import cPickle
     11  import time
     12
     13  #-----------------------------------------------------------------------------
     14 -class SubshapeCombineOperations(object):
     15    UNION=0
     16    SUM=1
     17    INTERSECT=2
     18
     19  #-----------------------------------------------------------------------------
     20 -class SubshapeBuilder(object):
     21    gridDims=(20,15,10)
     22    gridSpacing=0.5
     23    winRad=3.0
     24    nbrCount=7
     25    terminalPtRadScale=0.75
     26    fraction=0.25
     27    stepSize=1.0
     28    featFactory=None
     29
     30 -  def SampleSubshape(self,subshape1,newSpacing):
     31      ogrid=subshape1.grid
     32      rgrid = Geometry.UniformGrid3D(self.gridDims[0],self.gridDims[1],self.gridDims[2],
     33                                        newSpacing)
     34      for idx in range(rgrid.GetSize()):
     35        l = rgrid.GetGridPointLoc(idx)
     36        v = ogrid.GetValPoint(l)
     37        rgrid.SetVal(idx,v)
     38
     39      res = SubshapeObjects.ShapeWithSkeleton()
     40      res.grid = rgrid
     41      return res;
     42
     43 -  def GenerateSubshapeShape(self,cmpd,confId=-1,addSkeleton=True,**kwargs):
     44      shape = SubshapeObjects.ShapeWithSkeleton()
     45      shape.grid=Geometry.UniformGrid3D(self.gridDims[0],self.gridDims[1],self.gridDims[2],
     46                                        self.gridSpacing)
     47      AllChem.EncodeShape(cmpd,shape.grid,ignoreHs=False,confId=confId)
     48      if addSkeleton:
     49        conf = cmpd.GetConformer(confId)
     50        self.GenerateSubshapeSkeleton(shape,conf,kwargs)
     51      return shape
     52 -  def __call__(self,cmpd,**kwargs):
     53      return self.GenerateSubshapeShape(cmpd,**kwargs)
     54
     55 -  def GenerateSubshapeSkeleton(self,shape,conf=None,terminalPtsOnly=False,skelFromConf=True):
     56      if conf and skelFromConf:
     57        pts = BuilderUtils.FindTerminalPtsFromConformer(conf,self.winRad,self.nbrCount)
     58      else:
     59        pts = BuilderUtils.FindTerminalPtsFromShape(shape,self.winRad,self.fraction)
     60      pts = BuilderUtils.ClusterTerminalPts(pts,self.winRad,self.terminalPtRadScale)
     61      BuilderUtils.ExpandTerminalPts(shape,pts,self.winRad)
     62      if len(pts)<3:
     63        raise ValueError('only found %d terminals, need at least 3'%len(pts))
     64
     65      if not terminalPtsOnly:
     66        pts = BuilderUtils.AppendSkeletonPoints(shape.grid,pts,self.winRad,self.stepSize)
     67      for i,pt in enumerate(pts):
     68        BuilderUtils.CalculateDirectionsAtPoint(pt,shape.grid,self.winRad)
     69      if conf and self.featFactory:
     70        BuilderUtils.AssignMolFeatsToPoints(pts,conf.GetOwningMol(),self.featFactory,self.winRad)
     71      shape.skelPts=pts
     72
     73 -  def CombineSubshapes(self,subshape1,subshape2,operation=SubshapeCombineOperations.UNION):
     74      import copy
     75      cs = copy.deepcopy(subshape1)
     76      if operation==SubshapeCombineOperations.UNION:
     77        cs.grid |= subshape2.grid
     78      elif operation==SubshapeCombineOperations.SUM:
     79        cs.grid +=  subshape2.grid
     80      elif operation==SubshapeCombineOperations.INTERSECT:
     81        cs.grid &=  subshape2.grid
     82      else:
     83        raise ValueError('bad combination operation')
     84      return cs
     85
     86
     87  if __name__=='__main__':
     88    from rdkit.Chem import AllChem,ChemicalFeatures
     89    from rdkit.Chem.PyMol import MolViewer
     90    #cmpd = Chem.MolFromSmiles('CCCc1cc(C(=O)O)ccc1')
     91    #cmpd = Chem.AddHs(cmpd)
     92    if 1:
     93      cmpd = Chem.MolFromSmiles('C1=CC=C1C#CC1=CC=C1')
     94      cmpd = Chem.AddHs(cmpd)
     95      AllChem.EmbedMolecule(cmpd)
     96      AllChem.UFFOptimizeMolecule(cmpd)
     97      AllChem.CanonicalizeMol(cmpd)
     98      print >>file('testmol.mol','w+'),Chem.MolToMolBlock(cmpd)
     99    else:
    100      cmpd = Chem.MolFromMolFile('testmol.mol')
    101    builder=SubshapeBuilder()
    102    if 1:
    103      shape=builder.GenerateSubshapeShape(cmpd)
    104    v = MolViewer()
    105    if 1:
    106      import tempfile
    107      tmpFile = tempfile.mktemp('.grd')
    108      v.server.deleteAll()
    109      Geometry.WriteGridToFile(shape.grid,tmpFile)
    110      time.sleep(1)
    111      v.ShowMol(cmpd,name='testMol',showOnly=True)
    112      v.server.loadSurface(tmpFile,'testGrid','',2.5)
    113    v.server.resetCGO('*')
    114
    115    cPickle.dump(shape,file('subshape.pkl','w+'))
    116    for i,pt in enumerate(shape.skelPts):
    117      v.server.sphere(tuple(pt.location),.5,(1,0,1),'Pt-%d'%i)
    118      if not hasattr(pt,'shapeDirs'): continue
    119      momBeg = pt.location-pt.shapeDirs[0]
    120      momEnd = pt.location+pt.shapeDirs[0]
    121      v.server.cylinder(tuple(momBeg),tuple(momEnd),.1,(1,0,1),'v-%d'%i)
    122