1
2
3
4
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
18
19
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
42
54
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
85
86
87 if __name__=='__main__':
88 from rdkit.Chem import AllChem,ChemicalFeatures
89 from rdkit.Chem.PyMol import MolViewer
90
91
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