Package rdkit :: Package Chem :: Package Features :: Module FeatDirUtilsRD
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.Features.FeatDirUtilsRD

  1  # $Id$ 
  2  # 
  3  # Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC 
  4  # 
  5  #   @@ All Rights Reserved @@ 
  6  #  This file is part of the RDKit. 
  7  #  The contents are covered by the terms of the BSD license 
  8  #  which is included in the file license.txt, found at the root 
  9  #  of the RDKit source tree. 
 10  # 
 11  from rdkit import Geometry 
 12  from rdkit import Chem 
 13  import numpy 
 14  import math 
 15   
 16  # BIG NOTE: we are going assume atom IDs starting from 0 instead of 1 
 17  # for all the functions in this file. This is so that they 
 18  # are reasonably indepedent of the combicode. However when using 
 19  # with combicode the caller needs to make sure the atom IDs from combicode 
 20  # are corrected before feeding them in here. 
 21   
 22   
23 -def cross(v1, v2):
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
29 -def findNeighbors(atomId, adjMat):
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
44 -def _findAvgVec(conf, center, nbrs):
45 # find the average of the normalized vectors going from the center atoms to the 46 # neighbors 47 # the average vector is also normalized 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
63 -def GetAromaticFeatVects(conf, featAtoms, featLoc, scale=1.5):
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 #norm2 = norm1 85 norm2 = head - norm1 86 norm1 += head 87 return ((head, norm1), (head, norm2)), dirType
88 89
90 -def ArbAxisRotation(theta, ax, pt):
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
118 -def GetAcceptor2FeatVects(conf, featAtoms, scale=1.5):
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 # assume sp3 154 # we will create two vectors by rotating bvec by half the tetrahedral angle in either directions 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
172 -def _GetTetrahedralFeatVect(conf, aid, scale):
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
187 -def GetDonor3FeatVects(conf, featAtoms, scale=1.5):
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
206 -def GetAcceptor3FeatVects(conf, featAtoms, scale=1.5):
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
224 -def _findHydAtoms(nbrs, atomNames):
225 hAtoms = [] 226 for nid in nbrs: 227 if atomNames[nid] == 'H': 228 hAtoms.append(nid) 229 return hAtoms
230 231
232 -def _checkPlanarity(conf, cpt, nbrs, tol=1.0e-3):
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
248 -def GetDonor2FeatVects(conf, featAtoms, scale=1.5):
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 # find the two atoms that are neighbors of this atoms 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 # there should be no hydrogens in this case 284 assert len(hydrogens) == 0 285 # in this case the direction is the opposite of the average vector of the two neighbors 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 # this is a little more tricky we have to check if the hydrogen is in the plane of the 293 # two heavy atoms (i.e. sp2 arrangement) or out of plane (sp3 arrangement) 294 295 # one of the directions will be from hydrogen atom to the heavy atom 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 # only the hydrogen atom direction needs to be used 304 return ((cpt, bvec), ), 'linear' 305 else: 306 # we have a non-planar configuration - we will assume sp3 and compute a second direction vector 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 # in this case we should have two or more hydrogens we will simple use there directions 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
326 -def GetDonor1FeatVects(conf, featAtoms, scale=1.5):
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 # find the neighboring heavy atom 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
360 -def GetAcceptor1FeatVects(conf, featAtoms, scale=1.5):
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 # find the adjacent heavy atom 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 # special scale - if the heavy atom is a sulfur (we should proabably check phosphorous as well) 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 # ok in this case we will assume that 405 # heavy atom is sp2 hybridized and the direction vectors (two of them) 406 # are in the same plane, we will find this plane by looking for one 407 # of the neighbors of the heavy atom 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