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