Chem.Features.FeatDirUtilsRD

1 # $Id: FeatDirUtilsRD.py 52 2006-09-15 16:44:38Z glandrum $ 2 # 3 # Copyright (C) 2004-2006 Rational Discovery LLC 4 # 5 # @@ All Rights Reserved @@ 6 # 7 import Geometry 8 import Chem 9 from Numeric import * 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 = 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]],Float) 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 = array(mat) 94 if isinstance(pt,Geometry.Point3D): 95 pt = array((pt.x,pt.y,pt.z)) 96 tmp = matrixmultiply(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 = array((pt.x,pt.y,pt.z)) 103 tmp = matrixmultiply(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 = mol.GetAtomWithIdx(aid).GetNeighbors() 133 assert len(nbrs) == 2 134 135 bvec = _findAvgVec(conf, cpt, nbrs) 136 bvec *= (-1.0*scale) 137 138 if (mol.GetAtomWithIdx(aid).GetAtomicNum()==8): 139 # assume sp3 140 # we will create two vectors by rotating bvec by half the tetrahedral angle in either directions 141 v1 = conf.GetAtomPosition(nbrs[0].GetIdx()) 142 v1 -= cpt 143 v2 = conf.GetAtomPosition(nbrs[1].GetIdx()) 144 v2 -= cpt 145 rotAxis = v1 - v2 146 rotAxis.Normalize() 147 bv1 = ArbAxisRotation(54.5, rotAxis, bvec) 148 bv1 += cpt 149 bv2 = ArbAxisRotation(-54.5, rotAxis, bvec) 150 bv2 += cpt 151 return ((cpt, bv1), (cpt, bv2),), 'linear' 152 else : 153 bvec += cpt 154 return ((cpt, bvec),), 'linear'

155

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

157 mol = conf.GetOwningMol() 158 159 cpt = conf.GetAtomPosition(aid) 160 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors() 161 if not _checkPlanarity(conf,cpt,nbrs,tol=0.1): 162 bvec = _findAvgVec(conf, cpt, nbrs) 163 bvec *= (-1.0*scale) 164 bvec += cpt 165 res = ((cpt,bvec),) 166 else: 167 res = () 168 return res

169

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

171 """ 172 Get the direction vectors for Donor of type 3 173 174 This is a donor with three heavy atoms as neighbors. We will assume 175 a tetrahedral arrangement of these neighbors. So the direction we are seeking 176 is the last fourth arm of the sp3 arrangment 177 178 ARGUMENTS: 179 featAtoms - list of atoms that are part of the feature 180 scale - length of the direction vector 181 """ 182 assert len(featAtoms) == 1 183 aid = featAtoms[0] 184 185 tfv = _GetTetrahedralFeatVect(conf,aid,scale) 186 return tfv, 'linear'

187

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

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

204

205 -def _findHydAtoms(nbrs, atomNames):

206 hAtoms = [] 207 for nid in nbrs: 208 if atomNames[nid] == 'H': 209 hAtoms.append(nid) 210 return hAtoms

211

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

213 assert len(nbrs) == 3 214 v1 = conf.GetAtomPosition(nbrs[0].GetIdx()) 215 v1 -= cpt 216 v2 = conf.GetAtomPosition(nbrs[1].GetIdx()) 217 v2 -= cpt 218 v3 = conf.GetAtomPosition(nbrs[2].GetIdx()) 219 v3 -= cpt 220 normal = v1.CrossProduct(v2) 221 dotP = abs(v3.DotProduct(normal)) 222 if (dotP <= tol) : 223 return 1 224 else : 225 return 0

226

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

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

297

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

299 """ 300 Get the direction vectors for Donor of type 1 301 302 This is a donor with one heavy atom. It is not clear where we should we should be putting the 303 direction vector for this. It should probably be a cone. In this case we will just use the 304 direction vector from the donor atom to the heavy atom 305 306 ARGUMENTS: 307 308 featAtoms - list of atoms that are part of the feature 309 scale - length of the direction vector 310 """ 311 assert len(featAtoms) == 1 312 aid = featAtoms[0] 313 mol = conf.GetOwningMol() 314 nbrs = mol.GetAtomWithIdx(aid).GetNeighbors() 315 316 # find the neighboring heavy atom 317 hnbr = -1 318 for nbr in nbrs: 319 if nbr.GetAtomicNum()>1: 320 hnbr = nbr.GetIdx() 321 break 322 323 cpt = conf.GetAtomPosition(aid) 324 v1 = conf.GetAtomPosition(hnbr) 325 v1 -= cpt 326 v1.Normalize() 327 v1 *= (-1.0*scale) 328 v1 += cpt 329 return ((cpt, v1),), 'cone'

330

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

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

401

Source Code for Module Chem.Features.FeatDirUtilsRD