1
2
3
4
5
6 from __future__ import print_function
7 from rdkit import Chem,Geometry
8 from rdkit.Chem import AllChem
9 from rdkit.Chem.Subshape import SubshapeObjects
10 from rdkit.Chem.Subshape import BuilderUtils
11 from rdkit.six.moves import cPickle
12 import time
13
14
19
20
22 gridDims=(20,15,10)
23 gridSpacing=0.5
24 winRad=3.0
25 nbrCount=7
26 terminalPtRadScale=0.75
27 fraction=0.25
28 stepSize=1.0
29 featFactory=None
30
43
55
57 if conf and skelFromConf:
58 pts = BuilderUtils.FindTerminalPtsFromConformer(conf,self.winRad,self.nbrCount)
59 else:
60 pts = BuilderUtils.FindTerminalPtsFromShape(shape,self.winRad,self.fraction)
61 pts = BuilderUtils.ClusterTerminalPts(pts,self.winRad,self.terminalPtRadScale)
62 BuilderUtils.ExpandTerminalPts(shape,pts,self.winRad)
63 if len(pts)<3:
64 raise ValueError('only found %d terminals, need at least 3'%len(pts))
65
66 if not terminalPtsOnly:
67 pts = BuilderUtils.AppendSkeletonPoints(shape.grid,pts,self.winRad,self.stepSize)
68 for i,pt in enumerate(pts):
69 BuilderUtils.CalculateDirectionsAtPoint(pt,shape.grid,self.winRad)
70 if conf and self.featFactory:
71 BuilderUtils.AssignMolFeatsToPoints(pts,conf.GetOwningMol(),self.featFactory,self.winRad)
72 shape.skelPts=pts
73
86
87
88 if __name__=='__main__':
89 from rdkit.Chem import AllChem,ChemicalFeatures
90 from rdkit.Chem.PyMol import MolViewer
91
92
93 if 1:
94 cmpd = Chem.MolFromSmiles('C1=CC=C1C#CC1=CC=C1')
95 cmpd = Chem.AddHs(cmpd)
96 AllChem.EmbedMolecule(cmpd)
97 AllChem.UFFOptimizeMolecule(cmpd)
98 AllChem.CanonicalizeMol(cmpd)
99 print(Chem.MolToMolBlock(cmpd), file=file('testmol.mol','w+'))
100 else:
101 cmpd = Chem.MolFromMolFile('testmol.mol')
102 builder=SubshapeBuilder()
103 if 1:
104 shape=builder.GenerateSubshapeShape(cmpd)
105 v = MolViewer()
106 if 1:
107 import tempfile
108 tmpFile = tempfile.mktemp('.grd')
109 v.server.deleteAll()
110 Geometry.WriteGridToFile(shape.grid,tmpFile)
111 time.sleep(1)
112 v.ShowMol(cmpd,name='testMol',showOnly=True)
113 v.server.loadSurface(tmpFile,'testGrid','',2.5)
114 v.server.resetCGO('*')
115
116 cPickle.dump(shape,file('subshape.pkl','w+'))
117 for i,pt in enumerate(shape.skelPts):
118 v.server.sphere(tuple(pt.location),.5,(1,0,1),'Pt-%d'%i)
119 if not hasattr(pt,'shapeDirs'): continue
120 momBeg = pt.location-pt.shapeDirs[0]
121 momEnd = pt.location+pt.shapeDirs[0]
122 v.server.cylinder(tuple(momBeg),tuple(momEnd),.1,(1,0,1),'v-%d'%i)
123