1
2
3
4
5
6
7
8
9
10
11 from rdkit import Geometry
12 from rdkit import Chem
13 import numpy
14 import math
15
16
17
18
19
20
21
23 res = numpy.array([ v1[1]*v2[2] - v1[2]*v2[1],
24 -v1[0]*v2[2] + v1[2]*v2[0],
25 v1[0]*v2[1] - v1[1]*v2[0]],numpy.double)
26 return res
27
29 """
30 Find the IDs of the neighboring atom IDs
31
32 ARGUMENTS:
33 atomId - atom of interest
34 adjMat - adjacency matrix for the compound
35 """
36 nbrs = []
37 for i,nid in enumerate(adjMat[atomId]):
38 if nid >= 1 :
39 nbrs.append(i)
40 return nbrs
41
43
44
45
46 avgVec = 0
47 for nbr in nbrs:
48 nid = nbr.GetIdx()
49 pt = conf.GetAtomPosition(nid)
50 pt -= center
51 pt.Normalize()
52 if (avgVec == 0) :
53 avgVec = pt
54 else :
55 avgVec += pt
56
57 avgVec.Normalize()
58 return avgVec
59
61 """
62 Compute the direction vector for an aromatic feature
63
64 ARGUMENTS:
65 conf - a conformer
66 featAtoms - list of atom IDs that make up the feature
67 featLoc - location of the aromatic feature specified as point3d
68 scale - the size of the direction vector
69 """
70 dirType = 'linear'
71 head = featLoc
72 ats = [conf.GetAtomPosition(x) for x in featAtoms]
73
74 p0 = ats[0]
75 p1 = ats[1]
76 v1 = p0-head
77 v2 = p1-head
78 norm1 = v1.CrossProduct(v2)
79 norm1.Normalize()
80 norm1 *= scale
81
82 norm2 = head-norm1
83 norm1 += head
84 return ( (head,norm1),(head,norm2) ), dirType
85
87 theta = math.pi*theta/180
88 c = math.cos(theta)
89 s = math.sin(theta)
90 t = 1-c
91 X = ax.x
92 Y = ax.y
93 Z = ax.z
94 mat = [ [t*X*X+c, t*X*Y+s*Z, t*X*Z-s*Y],
95 [t*X*Y-s*Z,t*Y*Y+c,t*Y*Z+s*X],
96 [t*X*Z+s*Y,t*Y*Z-s*X,t*Z*Z+c] ]
97 mat = numpy.array(mat)
98 if isinstance(pt,Geometry.Point3D):
99 pt = numpy.array((pt.x,pt.y,pt.z))
100 tmp = numpy.dot(mat,pt)
101 res=Geometry.Point3D(tmp[0],tmp[1],tmp[2])
102 elif isinstance(pt,list) or isinstance(pt,tuple):
103 pts = pt
104 res = []
105 for pt in pts:
106 pt = numpy.array((pt.x,pt.y,pt.z))
107 tmp = numpy.dot(mat,pt)
108 res.append(Geometry.Point3D(tmp[0],tmp[1],tmp[2]))
109 else:
110 res=None
111 return res
112
113
114
115
116
118 """
119 Get the direction vectors for Acceptor of type 2
120
121 This is the acceptor with two adjacent heavy atoms. We will special case a few things here.
122 If the acceptor atom is an oxygen we will assume a sp3 hybridization
123 the acceptor directions (two of them)
124 reflect that configurations. Otherwise the direction vector in plane with the neighboring
125 heavy atoms
126
127 ARGUMENTS:
128 featAtoms - list of atoms that are part of the feature
129 scale - length of the direction vector
130 """
131 assert len(featAtoms) == 1
132 aid = featAtoms[0]
133 cpt = conf.GetAtomPosition(aid)
134
135 mol = conf.GetOwningMol()
136 nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors())
137 hydrogens = []
138 tmp=[]
139 while len(nbrs):
140 nbr = nbrs.pop()
141 if nbr.GetAtomicNum()==1:
142 hydrogens.append(nbr)
143 else:
144 tmp.append(nbr)
145 nbrs = tmp
146 assert len(nbrs) == 2
147
148 bvec = _findAvgVec(conf, cpt, nbrs)
149 bvec *= (-1.0*scale)
150
151 if (mol.GetAtomWithIdx(aid).GetAtomicNum()==8):
152
153
154 v1 = conf.GetAtomPosition(nbrs[0].GetIdx())
155 v1 -= cpt
156 v2 = conf.GetAtomPosition(nbrs[1].GetIdx())
157 v2 -= cpt
158 rotAxis = v1 - v2
159 rotAxis.Normalize()
160 bv1 = ArbAxisRotation(54.5, rotAxis, bvec)
161 bv1 += cpt
162 bv2 = ArbAxisRotation(-54.5, rotAxis, bvec)
163 bv2 += cpt
164 return ((cpt, bv1), (cpt, bv2),), 'linear'
165 else :
166 bvec += cpt
167 return ((cpt, bvec),), 'linear'
168
170 mol = conf.GetOwningMol()
171
172 cpt = conf.GetAtomPosition(aid)
173 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
174 if not _checkPlanarity(conf,cpt,nbrs,tol=0.1):
175 bvec = _findAvgVec(conf, cpt, nbrs)
176 bvec *= (-1.0*scale)
177 bvec += cpt
178 res = ((cpt,bvec),)
179 else:
180 res = ()
181 return res
182
184 """
185 Get the direction vectors for Donor of type 3
186
187 This is a donor with three heavy atoms as neighbors. We will assume
188 a tetrahedral arrangement of these neighbors. So the direction we are seeking
189 is the last fourth arm of the sp3 arrangment
190
191 ARGUMENTS:
192 featAtoms - list of atoms that are part of the feature
193 scale - length of the direction vector
194 """
195 assert len(featAtoms) == 1
196 aid = featAtoms[0]
197
198 tfv = _GetTetrahedralFeatVect(conf,aid,scale)
199 return tfv, 'linear'
200
202 """
203 Get the direction vectors for Donor of type 3
204
205 This is a donor with three heavy atoms as neighbors. We will assume
206 a tetrahedral arrangement of these neighbors. So the direction we are seeking
207 is the last fourth arm of the sp3 arrangment
208
209 ARGUMENTS:
210 featAtoms - list of atoms that are part of the feature
211 scale - length of the direction vector
212 """
213 assert len(featAtoms) == 1
214 aid = featAtoms[0]
215 tfv = _GetTetrahedralFeatVect(conf,aid,scale)
216 return tfv, 'linear'
217
219 hAtoms = []
220 for nid in nbrs:
221 if atomNames[nid] == 'H':
222 hAtoms.append(nid)
223 return hAtoms
224
226 assert len(nbrs) == 3
227 v1 = conf.GetAtomPosition(nbrs[0].GetIdx())
228 v1 -= cpt
229 v2 = conf.GetAtomPosition(nbrs[1].GetIdx())
230 v2 -= cpt
231 v3 = conf.GetAtomPosition(nbrs[2].GetIdx())
232 v3 -= cpt
233 normal = v1.CrossProduct(v2)
234 dotP = abs(v3.DotProduct(normal))
235 if (dotP <= tol) :
236 return 1
237 else :
238 return 0
239
241 """
242 Get the direction vectors for Donor of type 2
243
244 This is a donor with two heavy atoms as neighbors. The atom may are may not have
245 hydrogen on it. Here are the situations with the neighbors that will be considered here
246 1. two heavy atoms and two hydrogens: we will assume a sp3 arrangement here
247 2. two heavy atoms and one hydrogen: this can either be sp2 or sp3
248 3. two heavy atoms and no hydrogens
249
250 ARGUMENTS:
251 featAtoms - list of atoms that are part of the feature
252 scale - length of the direction vector
253 """
254 assert len(featAtoms) == 1
255 aid = featAtoms[0]
256 mol = conf.GetOwningMol()
257
258 cpt = conf.GetAtomPosition(aid)
259
260
261 nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors())
262 assert len(nbrs) >= 2
263
264 hydrogens = []
265 tmp=[]
266 while len(nbrs):
267 nbr = nbrs.pop()
268 if nbr.GetAtomicNum()==1:
269 hydrogens.append(nbr)
270 else:
271 tmp.append(nbr)
272 nbrs = tmp
273
274 if len(nbrs) == 2:
275
276 assert len(hydrogens) == 0
277
278 bvec = _findAvgVec(conf, cpt, nbrs)
279 bvec *= (-1.0*scale)
280 bvec += cpt
281 return ((cpt, bvec),), 'linear'
282 elif len(nbrs) == 3:
283 assert len(hydrogens) == 1
284
285
286
287
288 hid = hydrogens[0].GetIdx()
289 bvec = conf.GetAtomPosition(hid)
290 bvec -= cpt
291 bvec.Normalize()
292 bvec *= scale
293 bvec += cpt
294 if _checkPlanarity(conf, cpt, nbrs):
295
296 return ((cpt, bvec),), 'linear'
297 else :
298
299 ovec = _findAvgVec(conf, cpt, nbrs)
300 ovec *= (-1.0*scale)
301 ovec += cpt
302 return ((cpt, bvec), (cpt, ovec),), 'linear'
303
304 elif len(nbrs) >= 4 :
305
306 res = []
307 for hid in hydrogens:
308 bvec = conf.GetAtomPosition(hid)
309 bvec -= cpt
310 bvec.Normalize()
311 bvec *= scale
312 bvec += cpt
313 res.append((cpt, bvec))
314 return tuple(res), 'linear'
315
317 """
318 Get the direction vectors for Donor of type 1
319
320 This is a donor with one heavy atom. It is not clear where we should we should be putting the
321 direction vector for this. It should probably be a cone. In this case we will just use the
322 direction vector from the donor atom to the heavy atom
323
324 ARGUMENTS:
325
326 featAtoms - list of atoms that are part of the feature
327 scale - length of the direction vector
328 """
329 assert len(featAtoms) == 1
330 aid = featAtoms[0]
331 mol = conf.GetOwningMol()
332 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
333
334
335 hnbr = -1
336 for nbr in nbrs:
337 if nbr.GetAtomicNum()!=1:
338 hnbr = nbr.GetIdx()
339 break
340
341 cpt = conf.GetAtomPosition(aid)
342 v1 = conf.GetAtomPosition(hnbr)
343 v1 -= cpt
344 v1.Normalize()
345 v1 *= (-1.0*scale)
346 v1 += cpt
347 return ((cpt, v1),), 'cone'
348
350 """
351 Get the direction vectors for Acceptor of type 1
352
353 This is a acceptor with one heavy atom neighbor. There are two possibilities we will
354 consider here
355 1. The bond to the heavy atom is a single bond e.g. CO
356 In this case we don't know the exact direction and we just use the inversion of this bond direction
357 and mark this direction as a 'cone'
358 2. The bond to the heavy atom is a double bond e.g. C=O
359 In this case the we have two possible direction except in some special cases e.g. SO2
360 where again we will use bond direction
361
362 ARGUMENTS:
363 featAtoms - list of atoms that are part of the feature
364 scale - length of the direction vector
365 """
366 assert len(featAtoms) == 1
367 aid = featAtoms[0]
368 mol = conf.GetOwningMol()
369 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors()
370
371 cpt = conf.GetAtomPosition(aid)
372
373
374 heavyAt = -1
375 for nbr in nbrs:
376 if nbr.GetAtomicNum()!=1:
377 heavyAt = nbr
378 break
379
380 singleBnd = mol.GetBondBetweenAtoms(aid,heavyAt.GetIdx()).GetBondType() > Chem.BondType.SINGLE
381
382
383 sulfur = heavyAt.GetAtomicNum()==16
384
385 if singleBnd or sulfur:
386 v1 = conf.GetAtomPosition(heavyAt.GetIdx())
387 v1 -= cpt
388 v1.Normalize()
389 v1 *= (-1.0*scale)
390 v1 += cpt
391 return ((cpt, v1),), 'cone'
392 else :
393
394
395
396
397 hvNbrs = heavyAt.GetNeighbors()
398 hvNbr = -1
399 for nbr in hvNbrs:
400 if nbr.GetIdx() != aid:
401 hvNbr = nbr
402 break
403
404 pt1 = conf.GetAtomPosition(hvNbr.GetIdx())
405 v1 = conf.GetAtomPosition(heavyAt.GetIdx())
406 pt1 -= v1
407 v1 -= cpt
408 rotAxis = v1.CrossProduct(pt1)
409 rotAxis.Normalize()
410 bv1 = ArbAxisRotation(120, rotAxis, v1)
411 bv1.Normalize()
412 bv1 *= scale
413 bv1 += cpt
414 bv2 = ArbAxisRotation(-120, rotAxis, v1)
415 bv2.Normalize()
416 bv2 *= scale
417 bv2 += cpt
418 return ((cpt, bv1), (cpt, bv2),), 'linear'
419