1
2
3
4
5
6
7
8
9 """ Import all RDKit chemistry modules
10
11 """
12 import rdBase
13 import RDConfig
14 import numpy.oldnumeric as Numeric
15 import DataStructs
16 from Geometry import rdGeometry
17 from Chem import *
18 from rdPartialCharges import *
19 from rdDepictor import *
20 from rdForceFieldHelpers import *
21 from Chem.ChemicalFeatures import *
22 from rdDistGeom import *
23 from rdMolAlign import *
24 from rdMolTransforms import *
25 from rdShapeHelpers import *
26 from rdChemReactions import *
27 import ForceField
28 Mol.Compute2DCoords = Compute2DCoords
29 Mol.ComputeGasteigerCharges = ComputeGasteigerCharges
30 import numpy
31
48
53
73
86
87 -def GetBestRMS(ref,probe,refConfId=-1,probeConfId=-1,maps=None):
88 if not maps:
89 query = RemoveHs(probe)
90 matches = ref.GetSubstructMatches(query)
91 if not matches:
92 raise ValueError,'mol %s does not match mol %s'%(ref.GetProp('_Name'),
93 probe.GetProp('_Name'))
94 maps = []
95 for match in matches:
96 t=[]
97 for j,idx in enumerate(match):
98 t.append((j,idx))
99 maps.append(t)
100 bestRMS=1000.
101 for amap in maps:
102 rms=AlignMol(probe,ref,probeConfId,refConfId,atomMap=amap)
103 if rms<bestRMS:
104 bestRMS=rms
105 return bestRMS
106
108 """ Returns a generator for the virtual library defined by
109 a reaction and a sequence of sidechain sets
110
111 >>> import Chem
112 >>> from Chem import AllChem
113 >>> s1=[Chem.MolFromSmiles(x) for x in ('NC','NCC')]
114 >>> s2=[Chem.MolFromSmiles(x) for x in ('OC=O','OC(=O)C')]
115 >>> rxn = AllChem.ReactionFromSmarts('[O:2]=[C:1][OH].[N:3]>>[O:2]=[C:1][N:3]')
116 >>> r = AllChem.EnumerateLibraryFromReaction(rxn,[s2,s1])
117 >>> [Chem.MolToSmiles(x[0]) for x in list(r)]
118 ['CNC=O', 'CCNC=O', 'CNC(=O)C', 'CCNC(=O)C']
119
120 Note that this is all done in a lazy manner, so "infinitely" large libraries can
121 be done without worrying about running out of memory. Your patience will run out first:
122
123 Define a set of 10000 amines:
124 >>> amines = (Chem.MolFromSmiles('N'+'C'*x) for x in range(10000))
125
126 ... a set of 10000 acids
127 >>> acids = (Chem.MolFromSmiles('OC(=O)'+'C'*x) for x in range(10000))
128
129 ... now the virtual library (1e8 compounds in principle):
130 >>> r = AllChem.EnumerateLibraryFromReaction(rxn,[acids,amines])
131
132 ... look at the first 4 compounds:
133 >>> [Chem.MolToSmiles(r.next()[0]) for x in range(4)]
134 ['NC=O', 'CNC=O', 'CCNC=O', 'CCCNC=O']
135
136
137 """
138 if len(sidechainSets) != reaction.GetNumReactantTemplates():
139 raise ValueError,'%d sidechains provided, %d required'%(len(sidechainSets),reaction.getNumReactantTemplates())
140
141 def _combiEnumerator(items,depth=0):
142 for item in items[depth]:
143 if depth+1 < len(items):
144 v = _combiEnumerator(items,depth+1)
145 for entry in v:
146 l=[item]
147 l.extend(entry)
148 yield l
149 else:
150 yield [item]
151 for chains in _combiEnumerator(sidechainSets):
152 prodSets = reaction.RunReactants(chains)
153 for prods in prodSets:
154 yield prods
155
156
157
159 match = mol.GetSubstructMatch(core)
160 if not match:
161 raise ValueError,"molecule doesn't match the core"
162 coordMap={}
163 coreConf = core.GetConformer()
164 for i,idxI in enumerate(match):
165 corePtI = coreConf.GetAtomPosition(i)
166 coordMap[idxI]=corePtI
167
168 ci = EmbedMolecule(mol,coordMap=coordMap,randomSeed=randomseed)
169 if ci<0:
170 logger.error('could not embed molecule %s, no coordinates generated.'%mol.GetProp('_Name'))
171
172 algMap=[]
173 for i,itm in enumerate(match):
174 algMap.append((itm,i))
175
176 if not useTethers:
177
178 ff = UFFGetMoleculeForceField(mol,confId=0)
179 for i,idxI in enumerate(match):
180 for j in range(i+1,len(match)):
181 idxJ = match[j]
182 d = coordMap[idxI].Distance(coordMap[idxJ])
183 ff.AddDistanceConstraint(idxI,idxJ,d,d,100.)
184 ff.Initialize()
185 n=4
186 more=ff.Minimize()
187 while more and n:
188 more=ff.Minimize()
189 n-=1
190
191 ssd =AlignMol(mol,core,atomMap=algMap)
192 else:
193
194 ssd = AlignMol(mol,core,atomMap=algMap)
195 ff = UFFGetMoleculeForceField(mol,confId=0)
196 conf = core.GetConformer()
197 for i in range(core.GetNumAtoms()):
198 p =conf.GetAtomPosition(i)
199 pIdx=ff.AddExtraPoint(p.x,p.y,p.z,fixed=True)-1
200 ff.AddDistanceConstraint(pIdx,match[i],0,0,100.)
201 ff.Initialize()
202 n=4
203 more=ff.Minimize(energyTol=1e-4,forceTol=1e-3)
204 while more and n:
205 more=ff.Minimize(energyTol=1e-4,forceTol=1e-3)
206 n-=1
207
208 ssd = AlignMol(mol,core,atomMap=algMap)
209 print numpy.sqrt(ssd/len(algMap))
210 return mol
211
212
213
214
215
216
218 import doctest,sys
219 return doctest.testmod(sys.modules["__main__"])
220
221
222 if __name__ == '__main__':
223 import sys
224 failed,tried = _test()
225 sys.exit(failed)
226