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 -def cross(v1,v2):
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
28 -def findNeighbors(atomId, adjMat):
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
42 -def _findAvgVec(conf, center, nbrs) :
43 # find the average of the normalized vectors going from the center atoms to the 44 # neighbors 45 # the average vector is also normalized 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
60 -def GetAromaticFeatVects(conf, featAtoms, featLoc, scale=1.5):
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 #norm2 = norm1 82 norm2 = head-norm1 83 norm1 += head 84 return ( (head,norm1),(head,norm2) ), dirType
85
86 -def ArbAxisRotation(theta,ax,pt):
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
117 -def GetAcceptor2FeatVects(conf, featAtoms, scale=1.5):
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 # assume sp3 153 # we will create two vectors by rotating bvec by half the tetrahedral angle in either directions 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
169 -def _GetTetrahedralFeatVect(conf,aid,scale):
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
183 -def GetDonor3FeatVects(conf, featAtoms, scale=1.5):
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
201 -def GetAcceptor3FeatVects(conf, featAtoms, scale=1.5):
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
218 -def _findHydAtoms(nbrs, atomNames):
219 hAtoms = [] 220 for nid in nbrs: 221 if atomNames[nid] == 'H': 222 hAtoms.append(nid) 223 return hAtoms
224
225 -def _checkPlanarity(conf, cpt, nbrs, tol=1.0e-3):
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
240 -def GetDonor2FeatVects(conf, featAtoms, scale=1.5) :
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 # find the two atoms that are neighbors of this atoms 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 # there should be no hydrogens in this case 276 assert len(hydrogens) == 0 277 # in this case the direction is the opposite of the average vector of the two neighbors 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 # this is a little more tricky we have to check if the hydrogen is in the plane of the 285 # two heavy atoms (i.e. sp2 arrangement) or out of plane (sp3 arrangement) 286 287 # one of the directions will be from hydrogen atom to the heavy atom 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 # only the hydrogen atom direction needs to be used 296 return ((cpt, bvec),), 'linear' 297 else : 298 # we have a non-planar configuration - we will assume sp3 and compute a second direction vector 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 # in this case we should have two or more hydrogens we will simple use there directions 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
316 -def GetDonor1FeatVects(conf, featAtoms, scale=1.5) :
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 # find the neighboring heavy atom 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
349 -def GetAcceptor1FeatVects(conf, featAtoms, scale=1.5) :
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 # find the adjacent heavy atom 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 # special scale - if the heavy atom is a sulfur (we should proabably check phosphorous as well) 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 # ok in this case we will assume that 394 # heavy atom is sp2 hybridized and the direction vectors (two of them) 395 # are in the same plane, we will find this plane by looking for one 396 # of the neighbors of the heavy atom 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