rdkit.Chem.Features.FeatDirUtilsRD

1 # $Id: FeatDirUtilsRD.py 997 2009-02-25 06:12:43Z glandrum $ 2 # 3 # Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC 4 # 5 # @@ All Rights Reserved @@ 6 # 7 from rdkit import Geometry 8 from rdkit import Chem 9 import numpy 10 import math 11 12 # BIG NOTE: we are going assume atom IDs starting from 0 instead of 1 13 # for all the functions in this file. This is so that they 14 # are reasonably indepedent of the combicode. However when using 15 # with combicode the caller needs to make sure the atom IDs from combicode 16 # are corrected before feeding them in here. 17

18 -def cross(v1,v2):

19 res = numpy.array([ v1[1]*v2[2] - v1[2]*v2[1], 20 -v1[0]*v2[2] + v1[2]*v2[0], 21 v1[0]*v2[1] - v1[1]*v2[0]],numpy.double) 22 return res

23

24 -def findNeighbors(atomId, adjMat):

25 """ 26 Find the IDs of the neighboring atom IDs 27 28 ARGUMENTS: 29 atomId - atom of interest 30 adjMat - adjacency matrix for the compound 31 """ 32 nbrs = [] 33 for i,nid in enumerate(adjMat[atomId]): 34 if nid >= 1 : 35 nbrs.append(i) 36 return nbrs

37

38 -def _findAvgVec(conf, center, nbrs) :

39 # find the average of the normalized vectors going from the center atoms to the 40 # neighbors 41 # the average vector is also normalized 42 avgVec = 0 43 for nbr in nbrs: 44 nid = nbr.GetIdx() 45 pt = conf.GetAtomPosition(nid) 46 pt -= center 47 pt.Normalize() 48 if (avgVec == 0) : 49 avgVec = pt 50 else : 51 avgVec += pt 52 53 avgVec.Normalize() 54 return avgVec

55

56 -def GetAromaticFeatVects(conf, featAtoms, featLoc, scale=1.5):

57 """ 58 Compute the direction vector for an aromatic feature 59 60 ARGUMENTS: 61 conf - a conformer 62 featAtoms - list of atom IDs that make up the feature 63 featLoc - location of the aromatic feature specified as point3d 64 scale - the size of the direction vector 65 """ 66 dirType = 'linear' 67 head = featLoc 68 ats = [conf.GetAtomPosition(x) for x in featAtoms] 69 70 p0 = ats[0] 71 p1 = ats[1] 72 v1 = p0-head 73 v2 = p1-head 74 norm1 = v1.CrossProduct(v2) 75 norm1.Normalize() 76 norm1 *= scale 77 #norm2 = norm1 78 norm2 = head-norm1 79 norm1 += head 80 return ( (head,norm1),(head,norm2) ), dirType

81

82 -def ArbAxisRotation(theta,ax,pt):

83 theta = math.pi*theta/180 84 c = math.cos(theta) 85 s = math.sin(theta) 86 t = 1-c 87 X = ax.x 88 Y = ax.y 89 Z = ax.z 90 mat = [ [t*X*X+c, t*X*Y+s*Z, t*X*Z-s*Y], 91 [t*X*Y-s*Z,t*Y*Y+c,t*Y*Z+s*X], 92 [t*X*Z+s*Y,t*Y*Z-s*X,t*Z*Z+c] ] 93 mat = numpy.array(mat) 94 if isinstance(pt,Geometry.Point3D): 95 pt = numpy.array((pt.x,pt.y,pt.z)) 96 tmp = numpy.dot(mat,pt) 97 res=Geometry.Point3D(tmp[0],tmp[1],tmp[2]) 98 elif isinstance(pt,list) or isinstance(pt,tuple): 99 pts = pt 100 res = [] 101 for pt in pts: 102 pt = numpy.array((pt.x,pt.y,pt.z)) 103 tmp = numpy.dot(mat,pt) 104 res.append(Geometry.Point3D(tmp[0],tmp[1],tmp[2])) 105 else: 106 res=None 107 return res

108 109 110 111 112

113 -def GetAcceptor2FeatVects(conf, featAtoms, scale=1.5):

114 """ 115 Get the direction vectors for Acceptor of type 2 116 117 This is the acceptor with two adjacent heavy atoms. We will special case a few things here. 118 If the acceptor atom is an oxygen we will assume a sp3 hybridization 119 the acceptor directions (two of them) 120 reflect that configurations. Otherwise the direction vector in plane with the neighboring 121 heavy atoms 122 123 ARGUMENTS: 124 featAtoms - list of atoms that are part of the feature 125 scale - length of the direction vector 126 """ 127 assert len(featAtoms) == 1 128 aid = featAtoms[0] 129 cpt = conf.GetAtomPosition(aid) 130 131 mol = conf.GetOwningMol() 132 nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors()) 133 hydrogens = [] 134 tmp=[] 135 while len(nbrs): 136 nbr = nbrs.pop() 137 if nbr.GetAtomicNum()==1: 138 hydrogens.append(nbr) 139 else: 140 tmp.append(nbr) 141 nbrs = tmp 142 assert len(nbrs) == 2 143 144 bvec = _findAvgVec(conf, cpt, nbrs) 145 bvec *= (-1.0*scale) 146 147 if (mol.GetAtomWithIdx(aid).GetAtomicNum()==8): 148 # assume sp3 149 # we will create two vectors by rotating bvec by half the tetrahedral angle in either directions 150 v1 = conf.GetAtomPosition(nbrs[0].GetIdx()) 151 v1 -= cpt 152 v2 = conf.GetAtomPosition(nbrs[1].GetIdx()) 153 v2 -= cpt 154 rotAxis = v1 - v2 155 rotAxis.Normalize() 156 bv1 = ArbAxisRotation(54.5, rotAxis, bvec) 157 bv1 += cpt 158 bv2 = ArbAxisRotation(-54.5, rotAxis, bvec) 159 bv2 += cpt 160 return ((cpt, bv1), (cpt, bv2),), 'linear' 161 else : 162 bvec += cpt 163 return ((cpt, bvec),), 'linear'

164

165 -def _GetTetrahedralFeatVect(conf,aid,scale):

166 mol = conf.GetOwningMol() 167 168 cpt = conf.GetAtomPosition(aid) 169 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors() 170 if not _checkPlanarity(conf,cpt,nbrs,tol=0.1): 171 bvec = _findAvgVec(conf, cpt, nbrs) 172 bvec *= (-1.0*scale) 173 bvec += cpt 174 res = ((cpt,bvec),) 175 else: 176 res = () 177 return res

178

179 -def GetDonor3FeatVects(conf, featAtoms, scale=1.5):

180 """ 181 Get the direction vectors for Donor of type 3 182 183 This is a donor with three heavy atoms as neighbors. We will assume 184 a tetrahedral arrangement of these neighbors. So the direction we are seeking 185 is the last fourth arm of the sp3 arrangment 186 187 ARGUMENTS: 188 featAtoms - list of atoms that are part of the feature 189 scale - length of the direction vector 190 """ 191 assert len(featAtoms) == 1 192 aid = featAtoms[0] 193 194 tfv = _GetTetrahedralFeatVect(conf,aid,scale) 195 return tfv, 'linear'

196

197 -def GetAcceptor3FeatVects(conf, featAtoms, scale=1.5):

198 """ 199 Get the direction vectors for Donor of type 3 200 201 This is a donor with three heavy atoms as neighbors. We will assume 202 a tetrahedral arrangement of these neighbors. So the direction we are seeking 203 is the last fourth arm of the sp3 arrangment 204 205 ARGUMENTS: 206 featAtoms - list of atoms that are part of the feature 207 scale - length of the direction vector 208 """ 209 assert len(featAtoms) == 1 210 aid = featAtoms[0] 211 tfv = _GetTetrahedralFeatVect(conf,aid,scale) 212 return tfv, 'linear'

213

214 -def _findHydAtoms(nbrs, atomNames):

215 hAtoms = [] 216 for nid in nbrs: 217 if atomNames[nid] == 'H': 218 hAtoms.append(nid) 219 return hAtoms

220

221 -def _checkPlanarity(conf, cpt, nbrs, tol=1.0e-3):

222 assert len(nbrs) == 3 223 v1 = conf.GetAtomPosition(nbrs[0].GetIdx()) 224 v1 -= cpt 225 v2 = conf.GetAtomPosition(nbrs[1].GetIdx()) 226 v2 -= cpt 227 v3 = conf.GetAtomPosition(nbrs[2].GetIdx()) 228 v3 -= cpt 229 normal = v1.CrossProduct(v2) 230 dotP = abs(v3.DotProduct(normal)) 231 if (dotP <= tol) : 232 return 1 233 else : 234 return 0

235

236 -def GetDonor2FeatVects(conf, featAtoms, scale=1.5) :

237 """ 238 Get the direction vectors for Donor of type 2 239 240 This is a donor with two heavy atoms as neighbors. The atom may are may not have 241 hydrogen on it. Here are the situations with the neighbors that will be considered here 242 1. two heavy atoms and two hydrogens: we will assume a sp3 arrangement here 243 2. two heavy atoms and one hydrogen: this can either be sp2 or sp3 244 3. two heavy atoms and no hydrogens 245 246 ARGUMENTS: 247 featAtoms - list of atoms that are part of the feature 248 scale - length of the direction vector 249 """ 250 assert len(featAtoms) == 1 251 aid = featAtoms[0] 252 mol = conf.GetOwningMol() 253 254 cpt = conf.GetAtomPosition(aid) 255 256 # find the two atoms that are neighbors of this atoms 257 nbrs = list(mol.GetAtomWithIdx(aid).GetNeighbors()) 258 assert len(nbrs) >= 2 259 260 hydrogens = [] 261 tmp=[] 262 while len(nbrs): 263 nbr = nbrs.pop() 264 if nbr.GetAtomicNum()==1: 265 hydrogens.append(nbr) 266 else: 267 tmp.append(nbr) 268 nbrs = tmp 269 270 if len(nbrs) == 2: 271 # there should be no hydrogens in this case 272 assert len(hydrogens) == 0 273 # in this case the direction is the opposite of the average vector of the two neighbors 274 bvec = _findAvgVec(conf, cpt, nbrs) 275 bvec *= (-1.0*scale) 276 bvec += cpt 277 return ((cpt, bvec),), 'linear' 278 elif len(nbrs) == 3: 279 assert len(hydrogens) == 1 280 # this is a little more tricky we have to check if the hydrogen is in the plane of the 281 # two heavy atoms (i.e. sp2 arrangement) or out of plane (sp3 arrangement) 282 283 # one of the directions will be from hydrogen atom to the heavy atom 284 hid = hydrogens[0].GetIdx() 285 bvec = conf.GetAtomPosition(hid) 286 bvec -= cpt 287 bvec.Normalize() 288 bvec *= scale 289 bvec += cpt 290 if _checkPlanarity(conf, cpt, nbrs): 291 # only the hydrogen atom direction needs to be used 292 return ((cpt, bvec),), 'linear' 293 else : 294 # we have a non-planar configuration - we will assume sp3 and compute a second direction vector 295 ovec = _findAvgVec(conf, cpt, nbrs) 296 ovec *= (-1.0*scale) 297 ovec += cpt 298 return ((cpt, bvec), (cpt, ovec),), 'linear' 299 300 elif len(nbrs) >= 4 : 301 # in this case we should have two or more hydrogens we will simple use there directions 302 res = [] 303 for hid in hydrogenss: 304 bvec = numpy.array(coords[3*(hid):3*hid+3]) 305 bvec -= cpt 306 bvec /= (numpy.sqrt(numpy.dot(bvec, bvec))) 307 bvec *= scale 308 bvec += cpt 309 res.append((cpt, bvec)) 310 return tuple(res), 'linear'

311

312 -def GetDonor1FeatVects(conf, featAtoms, scale=1.5) :

313 """ 314 Get the direction vectors for Donor of type 1 315 316 This is a donor with one heavy atom. It is not clear where we should we should be putting the 317 direction vector for this. It should probably be a cone. In this case we will just use the 318 direction vector from the donor atom to the heavy atom 319 320 ARGUMENTS: 321 322 featAtoms - list of atoms that are part of the feature 323 scale - length of the direction vector 324 """ 325 assert len(featAtoms) == 1 326 aid = featAtoms[0] 327 mol = conf.GetOwningMol() 328 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors() 329 330 # find the neighboring heavy atom 331 hnbr = -1 332 for nbr in nbrs: 333 if nbr.GetAtomicNum()!=1: 334 hnbr = nbr.GetIdx() 335 break 336 337 cpt = conf.GetAtomPosition(aid) 338 v1 = conf.GetAtomPosition(hnbr) 339 v1 -= cpt 340 v1.Normalize() 341 v1 *= (-1.0*scale) 342 v1 += cpt 343 return ((cpt, v1),), 'cone'

344

345 -def GetAcceptor1FeatVects(conf, featAtoms, scale=1.5) :

346 """ 347 Get the direction vectors for Acceptor of type 1 348 349 This is a acceptor with one heavy atom neighbor. There are two possibilities we will 350 consider here 351 1. The bond to the heavy atom is a single bond e.g. CO 352 In this case we don't know the exact direction and we just use the inversion of this bond direction 353 and mark this direction as a 'cone' 354 2. The bond to the heavy atom is a double bond e.g. C=O 355 In this case the we have two possible direction except in some special cases e.g. SO2 356 where again we will use bond direction 357 358 ARGUMENTS: 359 featAtoms - list of atoms that are part of the feature 360 scale - length of the direction vector 361 """ 362 assert len(featAtoms) == 1 363 aid = featAtoms[0] 364 mol = conf.GetOwningMol() 365 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors() 366 367 cpt = conf.GetAtomPosition(aid) 368 369 # find the adjacent heavy atom 370 heavyAt = -1 371 for nbr in nbrs: 372 if nbr.GetAtomicNum()!=1: 373 heavyAt = nbr 374 break 375 376 singleBnd = mol.GetBondBetweenAtoms(aid,heavyAt.GetIdx()).GetBondType() > Chem.BondType.SINGLE 377 378 # special scale - if the heavy atom is a sulfur (we should proabably check phosphorous as well) 379 sulfur = heavyAt.GetAtomicNum()==16 380 381 if singleBnd or sulfur: 382 v1 = conf.GetAtomPosition(heavyAt.GetIdx()) 383 v1 -= cpt 384 v1.Normalize() 385 v1 *= (-1.0*scale) 386 v1 += cpt 387 return ((cpt, v1),), 'cone' 388 else : 389 # ok in this case we will assume that 390 # heavy atom is sp2 hybridized and the direction vectors (two of them) 391 # are in the same plane, we will find this plane by looking for one 392 # of the neighbors of the heavy atom 393 hvNbrs = heavyAt.GetNeighbors() 394 hvNbr = -1 395 for nbr in hvNbrs: 396 if nbr.GetIdx() != aid: 397 hvNbr = nbr 398 break 399 400 pt1 = conf.GetAtomPosition(hvNbr.GetIdx()) 401 v1 = conf.GetAtomPosition(heavyAt.GetIdx()) 402 pt1 -= v1 403 v1 -= cpt 404 rotAxis = v1.CrossProduct(pt1) 405 rotAxis.Normalize() 406 bv1 = ArbAxisRotation(120, rotAxis, v1) 407 bv1.Normalize() 408 bv1 *= scale 409 bv1 += cpt 410 bv2 = ArbAxisRotation(-120, rotAxis, v1) 411 bv2.Normalize() 412 bv2 *= scale 413 bv2 += cpt 414 return ((cpt, bv1), (cpt, bv2),), 'linear'

415

Source Code for Module rdkit.Chem.Features.FeatDirUtilsRD