Package rdkit :: Package Chem :: Package Pharm3D :: Module EmbedLib
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.Pharm3D.EmbedLib

   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 __future__ import print_function 
  12  from rdkit import RDConfig 
  13   
  14  import sys,time,math 
  15  from rdkit.ML.Data import Stats 
  16  import rdkit.DistanceGeometry as DG 
  17  from rdkit import Chem 
  18  import numpy 
  19  from rdkit.Chem import rdDistGeom as MolDG 
  20  from rdkit.Chem import ChemicalFeatures 
  21  from rdkit.Chem import ChemicalForceFields 
  22  import Pharmacophore,ExcludedVolume 
  23  from rdkit import Geometry 
  24  _times = {} 
  25   
  26  from rdkit import RDLogger as logging 
  27  logger = logging.logger() 
  28  defaultFeatLength=2.0 
  29   
30 -def GetAtomHeavyNeighbors(atom):
31 """ returns a list of the heavy-atom neighbors of the 32 atom passed in: 33 34 >>> m = Chem.MolFromSmiles('CCO') 35 >>> l = GetAtomHeavyNeighbors(m.GetAtomWithIdx(0)) 36 >>> len(l) 37 1 38 >>> isinstance(l[0],Chem.Atom) 39 True 40 >>> l[0].GetIdx() 41 1 42 43 >>> l = GetAtomHeavyNeighbors(m.GetAtomWithIdx(1)) 44 >>> len(l) 45 2 46 >>> l[0].GetIdx() 47 0 48 >>> l[1].GetIdx() 49 2 50 51 """ 52 res=[] 53 for nbr in atom.GetNeighbors(): 54 if nbr.GetAtomicNum() != 1: 55 res.append(nbr) 56 return res
57
58 -def ReplaceGroup(match,bounds,slop=0.01,useDirs=False,dirLength=defaultFeatLength):
59 """ Adds an entry at the end of the bounds matrix for a point at 60 the center of a multi-point feature 61 62 returns a 2-tuple: 63 new bounds mat 64 index of point added 65 66 >>> boundsMat = numpy.array([[0.0,2.0,2.0],[1.0,0.0,2.0],[1.0,1.0,0.0]]) 67 >>> match = [0,1,2] 68 >>> bm,idx = ReplaceGroup(match,boundsMat,slop=0.0) 69 70 the index is at the end: 71 >>> idx == 3 72 True 73 74 and the matrix is one bigger: 75 >>> bm.shape == (4, 4) 76 True 77 78 but the original bounds mat is not altered: 79 >>> boundsMat.shape == (3, 3) 80 True 81 82 83 We make the assumption that the points of the 84 feature form a regular polygon, are listed in order 85 (i.e. pt 0 is a neighbor to pt 1 and pt N-1) 86 and that the replacement point goes at the center: 87 >>> print(', '.join(['%.3f'%x for x in bm[-1]])) 88 0.577, 0.577, 0.577, 0.000 89 >>> print(', '.join(['%.3f'%x for x in bm[:,-1]])) 90 1.155, 1.155, 1.155, 0.000 91 92 The slop argument (default = 0.01) is fractional: 93 >>> bm,idx = ReplaceGroup(match,boundsMat) 94 >>> print(', '.join(['%.3f'%x for x in bm[-1]])) 95 0.572, 0.572, 0.572, 0.000 96 >>> print(', '.join(['%.3f'%x for x in bm[:,-1]])) 97 1.166, 1.166, 1.166, 0.000 98 99 100 """ 101 maxVal = -1000.0 102 minVal = 1e8 103 nPts = len(match) 104 for i in range(nPts): 105 idx0 = match[i] 106 if i<nPts-1: 107 idx1 = match[i+1] 108 else: 109 idx1 = match[0] 110 if idx1<idx0: 111 idx0,idx1 = idx1,idx0 112 minVal = min(minVal,bounds[idx1,idx0]) 113 maxVal = max(maxVal,bounds[idx0,idx1]) 114 maxVal *= (1+slop) 115 minVal *= (1-slop) 116 117 scaleFact = 1.0/(2.0*math.sin(math.pi/nPts)) 118 minVal *= scaleFact 119 maxVal *= scaleFact 120 121 replaceIdx = bounds.shape[0] 122 if not useDirs: 123 bm = numpy.zeros((bounds.shape[0]+1,bounds.shape[1]+1),numpy.float) 124 else: 125 bm = numpy.zeros((bounds.shape[0]+2,bounds.shape[1]+2),numpy.float) 126 bm[0:bounds.shape[0],0:bounds.shape[1]]=bounds 127 bm[:replaceIdx,replaceIdx]=1000. 128 129 if useDirs: 130 bm[:replaceIdx+1,replaceIdx+1]=1000. 131 # set the feature - direction point bounds: 132 bm[replaceIdx,replaceIdx+1]=dirLength+slop 133 bm[replaceIdx+1,replaceIdx]=dirLength-slop 134 135 for idx1 in match: 136 bm[idx1,replaceIdx]=maxVal 137 bm[replaceIdx,idx1]=minVal 138 if useDirs: 139 # set the point - direction point bounds: 140 bm[idx1,replaceIdx+1] = numpy.sqrt(bm[replaceIdx,replaceIdx+1]**2+maxVal**2) 141 bm[replaceIdx+1,idx1] = numpy.sqrt(bm[replaceIdx+1,replaceIdx]**2+minVal**2) 142 return bm,replaceIdx
143
144 -def EmbedMol(mol,bm,atomMatch=None,weight=2.0,randomSeed=-1, 145 excludedVolumes=None):
146 """ Generates an embedding for a molecule based on a bounds matrix and adds 147 a conformer (id 0) to the molecule 148 149 if the optional argument atomMatch is provided, it will be used to provide 150 supplemental weights for the embedding routine (used in the optimization 151 phase to ensure that the resulting geometry really does satisfy the 152 pharmacophore). 153 154 if the excludedVolumes is provided, it should be a sequence of 155 ExcludedVolume objects 156 157 >>> m = Chem.MolFromSmiles('c1ccccc1C') 158 >>> bounds = MolDG.GetMoleculeBoundsMatrix(m) 159 >>> bounds.shape == (7, 7) 160 True 161 >>> m.GetNumConformers() 162 0 163 >>> EmbedMol(m,bounds,randomSeed=23) 164 >>> m.GetNumConformers() 165 1 166 167 168 """ 169 nAts = mol.GetNumAtoms() 170 weights=[] 171 if(atomMatch): 172 for i in range(len(atomMatch)): 173 for j in range(i+1,len(atomMatch)): 174 weights.append((i,j,weight)) 175 if(excludedVolumes): 176 for vol in excludedVolumes: 177 idx = vol.index 178 # excluded volumes affect every other atom: 179 for i in range(nAts): 180 weights.append((i,idx,weight)) 181 coords = DG.EmbedBoundsMatrix(bm,weights=weights,numZeroFail=1,randomSeed=randomSeed) 182 #for row in coords: 183 # print(', '.join(['%.2f'%x for x in row])) 184 185 conf = Chem.Conformer(nAts) 186 conf.SetId(0) 187 for i in range(nAts): 188 conf.SetAtomPosition(i,list(coords[i])) 189 if excludedVolumes: 190 for vol in excludedVolumes: 191 vol.pos = numpy.array(coords[vol.index]) 192 193 #print(' % 7.4f % 7.4f % 7.4f Ar 0 0 0 0 0 0 0 0 0 0 0 0'%tuple(coords[-1]), file=sys.stderr) 194 mol.AddConformer(conf)
195
196 -def AddExcludedVolumes(bm,excludedVolumes,smoothIt=True):
197 """ Adds a set of excluded volumes to the bounds matrix 198 and returns the new matrix 199 200 excludedVolumes is a list of ExcludedVolume objects 201 202 203 >>> boundsMat = numpy.array([[0.0,2.0,2.0],[1.0,0.0,2.0],[1.0,1.0,0.0]]) 204 >>> ev1 = ExcludedVolume.ExcludedVolume(([(0,),0.5,1.0],),exclusionDist=1.5) 205 >>> bm = AddExcludedVolumes(boundsMat,(ev1,)) 206 207 the results matrix is one bigger: 208 >>> bm.shape == (4, 4) 209 True 210 211 and the original bounds mat is not altered: 212 >>> boundsMat.shape == (3, 3) 213 True 214 215 >>> print(', '.join(['%.3f'%x for x in bm[-1]])) 216 0.500, 1.500, 1.500, 0.000 217 >>> print(', '.join(['%.3f'%x for x in bm[:,-1]])) 218 1.000, 3.000, 3.000, 0.000 219 220 """ 221 oDim = bm.shape[0] 222 dim = oDim+len(excludedVolumes) 223 res = numpy.zeros((dim,dim),numpy.float) 224 res[:oDim,:oDim] = bm 225 for i,vol in enumerate(excludedVolumes): 226 bmIdx = oDim+i 227 vol.index = bmIdx 228 229 # set values to all the atoms: 230 res[bmIdx,:bmIdx] = vol.exclusionDist 231 res[:bmIdx,bmIdx] = 1000.0 232 233 # set values to our defining features: 234 for indices,minV,maxV in vol.featInfo: 235 for index in indices: 236 try: 237 res[bmIdx,index] = minV 238 res[index,bmIdx] = maxV 239 except IndexError: 240 logger.error('BAD INDEX: res[%d,%d], shape is %s'%(bmIdx,index,str(res.shape))) 241 raise IndexError 242 243 # set values to other excluded volumes: 244 for j in range(bmIdx+1,dim): 245 res[bmIdx,j:dim] = 0 246 res[j:dim,bmIdx] = 1000 247 248 if smoothIt: DG.DoTriangleSmoothing(res) 249 return res
250
251 -def UpdatePharmacophoreBounds(bm,atomMatch,pcophore,useDirs=False, 252 dirLength=defaultFeatLength, 253 mol=None):
254 """ loops over a distance bounds matrix and replaces the elements 255 that are altered by a pharmacophore 256 257 **NOTE** this returns the resulting bounds matrix, but it may also 258 alter the input matrix 259 260 atomMatch is a sequence of sequences containing atom indices 261 for each of the pharmacophore's features. 262 263 >>> feats = [ChemicalFeatures.FreeChemicalFeature('HBondAcceptor', 'HAcceptor1', Geometry.Point3D(0.0, 0.0, 0.0)), 264 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 265 ... ] 266 >>> pcophore=Pharmacophore.Pharmacophore(feats) 267 >>> pcophore.setLowerBound(0,1, 1.0) 268 >>> pcophore.setUpperBound(0,1, 2.0) 269 270 >>> boundsMat = numpy.array([[0.0,3.0,3.0],[2.0,0.0,3.0],[2.0,2.0,0.0]]) 271 >>> atomMatch = ((0,),(1,)) 272 >>> bm = UpdatePharmacophoreBounds(boundsMat,atomMatch,pcophore) 273 274 275 In this case, there are no multi-atom features, so the result matrix 276 is the same as the input: 277 >>> bm is boundsMat 278 True 279 280 this means, of course, that the input boundsMat is altered: 281 >>> print(', '.join(['%.3f'%x for x in boundsMat[0]])) 282 0.000, 2.000, 3.000 283 >>> print(', '.join(['%.3f'%x for x in boundsMat[1]])) 284 1.000, 0.000, 3.000 285 >>> print(', '.join(['%.3f'%x for x in boundsMat[2]])) 286 2.000, 2.000, 0.000 287 288 """ 289 replaceMap = {} 290 for i,matchI in enumerate(atomMatch): 291 if len(matchI)>1: 292 bm,replaceIdx = ReplaceGroup(matchI,bm,useDirs=useDirs) 293 replaceMap[i] = replaceIdx 294 295 for i,matchI in enumerate(atomMatch): 296 mi = replaceMap.get(i,matchI[0]) 297 for j in range(i+1,len(atomMatch)): 298 mj = replaceMap.get(j,atomMatch[j][0]) 299 if mi<mj: 300 idx0,idx1 = mi,mj 301 else: 302 idx0,idx1 = mj,mi 303 bm[idx0,idx1] = pcophore.getUpperBound(i,j) 304 bm[idx1,idx0] = pcophore.getLowerBound(i,j) 305 306 return bm
307 308 309
310 -def EmbedPharmacophore(mol,atomMatch,pcophore,randomSeed=-1,count=10,smoothFirst=True, 311 silent=False,bounds=None,excludedVolumes=None,targetNumber=-1, 312 useDirs=False):
313 """ Generates one or more embeddings for a molecule that satisfy a pharmacophore 314 315 atomMatch is a sequence of sequences containing atom indices 316 for each of the pharmacophore's features. 317 318 - count: is the maximum number of attempts to make a generating an embedding 319 - smoothFirst: toggles triangle smoothing of the molecular bounds matix 320 - bounds: if provided, should be the molecular bounds matrix. If this isn't 321 provided, the matrix will be generated. 322 - targetNumber: if this number is positive, it provides a maximum number 323 of embeddings to generate (i.e. we'll have count attempts to generate 324 targetNumber embeddings). 325 326 returns: a 3 tuple: 327 1) the molecular bounds matrix adjusted for the pharmacophore 328 2) a list of embeddings (molecules with a single conformer) 329 3) the number of failed attempts at embedding 330 331 332 >>> m = Chem.MolFromSmiles('OCCN') 333 >>> feats = [ChemicalFeatures.FreeChemicalFeature('HBondAcceptor', 'HAcceptor1', Geometry.Point3D(0.0, 0.0, 0.0)), 334 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 335 ... ] 336 >>> pcophore=Pharmacophore.Pharmacophore(feats) 337 >>> pcophore.setLowerBound(0,1, 2.5) 338 >>> pcophore.setUpperBound(0,1, 3.5) 339 >>> atomMatch = ((0,),(3,)) 340 341 >>> bm,embeds,nFail = EmbedPharmacophore(m,atomMatch,pcophore,randomSeed=23,silent=1) 342 >>> len(embeds) 343 10 344 >>> nFail 345 0 346 347 Set up a case that can't succeed: 348 >>> pcophore=Pharmacophore.Pharmacophore(feats) 349 >>> pcophore.setLowerBound(0,1, 2.0) 350 >>> pcophore.setUpperBound(0,1, 2.1) 351 >>> atomMatch = ((0,),(3,)) 352 353 >>> bm,embeds,nFail = EmbedPharmacophore(m,atomMatch,pcophore,randomSeed=23,silent=1) 354 >>> len(embeds) 355 0 356 >>> nFail 357 10 358 359 """ 360 global _times 361 if not hasattr(mol,'_chiralCenters'): 362 mol._chiralCenters = Chem.FindMolChiralCenters(mol) 363 364 if bounds is None: 365 bounds = MolDG.GetMoleculeBoundsMatrix(mol) 366 if smoothFirst: DG.DoTriangleSmoothing(bounds) 367 368 bm = bounds.copy() 369 #print '------------' 370 #print 'initial' 371 #for row in bm: 372 # print ' ',' '.join(['% 4.2f'%x for x in row]) 373 #print '------------' 374 375 bm = UpdatePharmacophoreBounds(bm,atomMatch,pcophore,useDirs=useDirs,mol=mol) 376 377 if excludedVolumes: 378 bm = AddExcludedVolumes(bm,excludedVolumes,smoothIt=False) 379 380 if not DG.DoTriangleSmoothing(bm): 381 raise ValueError("could not smooth bounds matrix") 382 383 #print '------------' 384 #print 'post replace and smooth' 385 #for row in bm: 386 # print ' ',' '.join(['% 4.2f'%x for x in row]) 387 #print '------------' 388 389 390 if targetNumber<=0: 391 targetNumber=count 392 nFailed = 0 393 res = [] 394 for i in range(count): 395 tmpM = bm[:,:] 396 m2 = Chem.Mol(mol) 397 t1 = time.time() 398 try: 399 if randomSeed<=0: 400 seed = i*10+1 401 else: 402 seed = i*10+randomSeed 403 EmbedMol(m2,tmpM,atomMatch,randomSeed=seed, 404 excludedVolumes=excludedVolumes) 405 except ValueError: 406 if not silent: 407 logger.info('Embed failed') 408 nFailed += 1 409 else: 410 t2 = time.time() 411 _times['embed'] = _times.get('embed',0)+t2-t1 412 keepIt=True 413 for idx,stereo in mol._chiralCenters: 414 if stereo in ('R','S'): 415 vol = ComputeChiralVolume(m2,idx) 416 if (stereo=='R' and vol>=0) or \ 417 (stereo=='S' and vol<=0): 418 keepIt=False 419 break 420 if keepIt: 421 res.append(m2) 422 else: 423 logger.debug('Removed embedding due to chiral constraints.') 424 if len(res)==targetNumber: break 425 return bm,res,nFailed
426
427 -def isNaN(v):
428 """ provides an OS independent way of detecting NaNs 429 This is intended to be used with values returned from the C++ 430 side of things. 431 432 We can't actually test this from Python (which traps 433 zero division errors), but it would work something like 434 this if we could: 435 >>> isNaN(0) 436 False 437 438 #>>> isNan(1/0) 439 #True 440 441 """ 442 if v!=v and sys.platform=='win32': 443 return True 444 elif v==0 and v==1 and sys.platform!='win32': 445 return True 446 return False
447 448
449 -def OptimizeMol(mol,bm,atomMatches=None,excludedVolumes=None, 450 forceConstant=1200.0, 451 maxPasses=5,verbose=False):
452 """ carries out a UFF optimization for a molecule optionally subject 453 to the constraints in a bounds matrix 454 455 - atomMatches, if provided, is a sequence of sequences 456 - forceConstant is the force constant of the spring used to enforce 457 the constraints 458 459 returns a 2-tuple: 460 1) the energy of the initial conformation 461 2) the energy post-embedding 462 NOTE that these energies include the energies of the constraints 463 464 465 466 >>> m = Chem.MolFromSmiles('OCCN') 467 >>> feats = [ChemicalFeatures.FreeChemicalFeature('HBondAcceptor', 'HAcceptor1', Geometry.Point3D(0.0, 0.0, 0.0)), 468 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 469 ... ] 470 >>> pcophore=Pharmacophore.Pharmacophore(feats) 471 >>> pcophore.setLowerBound(0,1, 2.5) 472 >>> pcophore.setUpperBound(0,1, 2.8) 473 >>> atomMatch = ((0,),(3,)) 474 >>> bm,embeds,nFail = EmbedPharmacophore(m,atomMatch,pcophore,randomSeed=23,silent=1) 475 >>> len(embeds) 476 10 477 >>> testM = embeds[0] 478 479 Do the optimization: 480 >>> e1,e2 = OptimizeMol(testM,bm,atomMatches=atomMatch) 481 482 Optimizing should have lowered the energy: 483 >>> e2 < e1 484 True 485 486 Check the constrained distance: 487 >>> conf = testM.GetConformer(0) 488 >>> p0 = conf.GetAtomPosition(0) 489 >>> p3 = conf.GetAtomPosition(3) 490 >>> d03 = p0.Distance(p3) 491 >>> d03 >= pcophore.getLowerBound(0,1)-.01 492 True 493 >>> d03 <= pcophore.getUpperBound(0,1)+.01 494 True 495 496 If we optimize without the distance constraints (provided via the atomMatches 497 argument) we're not guaranteed to get the same results, particularly in a case 498 like the current one where the pharmcophore brings the atoms uncomfortably 499 close together: 500 >>> testM = embeds[1] 501 >>> e1,e2 = OptimizeMol(testM,bm) 502 >>> e2 < e1 503 True 504 >>> conf = testM.GetConformer(0) 505 >>> p0 = conf.GetAtomPosition(0) 506 >>> p3 = conf.GetAtomPosition(3) 507 >>> d03 = p0.Distance(p3) 508 >>> d03 >= pcophore.getLowerBound(0,1)-.01 509 True 510 >>> d03 <= pcophore.getUpperBound(0,1)+.01 511 False 512 513 """ 514 try: 515 ff = ChemicalForceFields.UFFGetMoleculeForceField(mol) 516 except Exception: 517 logger.info('Problems building molecular forcefield',exc_info=True) 518 return -1.0,-1.0 519 520 weights=[] 521 if(atomMatches): 522 for k in range(len(atomMatches)): 523 for i in atomMatches[k]: 524 for l in range(k+1,len(atomMatches)): 525 for j in atomMatches[l]: 526 weights.append((i,j)) 527 for i,j in weights: 528 if j<i: 529 i,j = j,i 530 minV = bm[j,i] 531 maxV = bm[i,j] 532 ff.AddDistanceConstraint(i,j,minV,maxV,forceConstant) 533 if excludedVolumes: 534 nAts = mol.GetNumAtoms() 535 conf = mol.GetConformer() 536 idx = nAts 537 for exVol in excludedVolumes: 538 assert exVol.pos is not None 539 logger.debug('ff.AddExtraPoint(%.4f,%.4f,%.4f)'%(exVol.pos[0],exVol.pos[1], 540 exVol.pos[2])) 541 ff.AddExtraPoint(exVol.pos[0],exVol.pos[1],exVol.pos[2],True) 542 indices = [] 543 for localIndices,foo,bar in exVol.featInfo: 544 indices += list(localIndices) 545 for i in range(nAts): 546 v = numpy.array(conf.GetAtomPosition(i))-numpy.array(exVol.pos) 547 d = numpy.sqrt(numpy.dot(v,v)) 548 if i not in indices: 549 if d<5.0: 550 logger.debug('ff.AddDistanceConstraint(%d,%d,%.3f,%d,%.0f)'%(i,idx,exVol.exclusionDist,1000,forceConstant)) 551 ff.AddDistanceConstraint(i,idx,exVol.exclusionDist,1000, 552 forceConstant) 553 554 else: 555 logger.debug('ff.AddDistanceConstraint(%d,%d,%.3f,%.3f,%.0f)'%(i,idx,bm[exVol.index,i],bm[i,exVol.index],forceConstant)) 556 ff.AddDistanceConstraint(i,idx,bm[exVol.index,i],bm[i,exVol.index], 557 forceConstant) 558 idx += 1 559 ff.Initialize() 560 e1 = ff.CalcEnergy() 561 if isNaN(e1): 562 raise ValueError('bogus energy') 563 564 if verbose: 565 print(Chem.MolToMolBlock(mol)) 566 for i,vol in enumerate(excludedVolumes): 567 pos = ff.GetExtraPointPos(i) 568 print(' % 7.4f % 7.4f % 7.4f As 0 0 0 0 0 0 0 0 0 0 0 0'%tuple(pos), file=sys.stderr) 569 needsMore=ff.Minimize() 570 nPasses=0 571 while needsMore and nPasses<maxPasses: 572 needsMore=ff.Minimize() 573 nPasses+=1 574 e2 = ff.CalcEnergy() 575 if isNaN(e2): 576 raise ValueError('bogus energy') 577 578 if verbose: 579 print('--------') 580 print(Chem.MolToMolBlock(mol)) 581 for i,vol in enumerate(excludedVolumes): 582 pos = ff.GetExtraPointPos(i) 583 print(' % 7.4f % 7.4f % 7.4f Sb 0 0 0 0 0 0 0 0 0 0 0 0'%tuple(pos), file=sys.stderr) 584 ff = None 585 return e1,e2
586
587 -def EmbedOne(mol,name,match,pcophore,count=1,silent=0,**kwargs):
588 """ generates statistics for a molecule's embeddings 589 590 Four energies are computed for each embedding: 591 1) E1: the energy (with constraints) of the initial embedding 592 2) E2: the energy (with constraints) of the optimized embedding 593 3) E3: the energy (no constraints) the geometry for E2 594 4) E4: the energy (no constraints) of the optimized free-molecule 595 (starting from the E3 geometry) 596 597 Returns a 9-tuple: 598 1) the mean value of E1 599 2) the sample standard deviation of E1 600 3) the mean value of E2 601 4) the sample standard deviation of E2 602 5) the mean value of E3 603 6) the sample standard deviation of E3 604 7) the mean value of E4 605 8) the sample standard deviation of E4 606 9) The number of embeddings that failed 607 608 """ 609 global _times 610 atomMatch = [list(x.GetAtomIds()) for x in match] 611 bm,ms,nFailed = EmbedPharmacophore(mol,atomMatch,pcophore,count=count, 612 silent=silent,**kwargs) 613 e1s = [] 614 e2s = [] 615 e3s = [] 616 e4s = [] 617 d12s = [] 618 d23s = [] 619 d34s = [] 620 for m in ms: 621 t1 = time.time() 622 try: 623 e1,e2 = OptimizeMol(m,bm,atomMatch) 624 except ValueError: 625 pass 626 else: 627 t2 = time.time() 628 _times['opt1'] = _times.get('opt1',0)+t2-t1 629 630 e1s.append(e1) 631 e2s.append(e2) 632 633 d12s.append(e1-e2) 634 t1 = time.time() 635 try: 636 e3,e4 = OptimizeMol(m,bm) 637 except ValueError: 638 pass 639 else: 640 t2 = time.time() 641 _times['opt2'] = _times.get('opt2',0)+t2-t1 642 e3s.append(e3) 643 e4s.append(e4) 644 d23s.append(e2-e3) 645 d34s.append(e3-e4) 646 count += 1 647 try: 648 e1,e1d = Stats.MeanAndDev(e1s) 649 except Exception: 650 e1 = -1.0 651 e1d=-1.0 652 try: 653 e2,e2d = Stats.MeanAndDev(e2s) 654 except Exception: 655 e2 = -1.0 656 e2d=-1.0 657 try: 658 e3,e3d = Stats.MeanAndDev(e3s) 659 except Exception: 660 e3 = -1.0 661 e3d=-1.0 662 663 try: 664 e4,e4d = Stats.MeanAndDev(e4s) 665 except Exception: 666 e4 = -1.0 667 e4d=-1.0 668 if not silent: 669 print('%s(%d): %.2f(%.2f) -> %.2f(%.2f) : %.2f(%.2f) -> %.2f(%.2f)' % 670 (name,nFailed,e1,e1d,e2,e2d,e3,e3d,e4,e4d)) 671 return e1,e1d,e2,e2d,e3,e3d,e4,e4d,nFailed
672
673 -def MatchPharmacophoreToMol(mol, featFactory, pcophore):
674 """ generates a list of all possible mappings of a pharmacophore to a molecule 675 676 Returns a 2-tuple: 677 1) a boolean indicating whether or not all features were found 678 2) a list, numFeatures long, of sequences of features 679 680 681 >>> import os.path 682 >>> dataDir = os.path.join(RDConfig.RDCodeDir,'Chem/Pharm3D/test_data') 683 >>> featFactory = ChemicalFeatures.BuildFeatureFactory(os.path.join(dataDir,'BaseFeatures.fdef')) 684 >>> activeFeats = [ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(0.0, 0.0, 0.0)), 685 ... ChemicalFeatures.FreeChemicalFeature('Donor',Geometry.Point3D(0.0, 0.0, 0.0))] 686 >>> pcophore= Pharmacophore.Pharmacophore(activeFeats) 687 >>> m = Chem.MolFromSmiles('FCCN') 688 >>> match,mList = MatchPharmacophoreToMol(m,featFactory,pcophore) 689 >>> match 690 True 691 692 Two feature types: 693 >>> len(mList) 694 2 695 696 The first feature type, Acceptor, has two matches: 697 >>> len(mList[0]) 698 2 699 >>> mList[0][0].GetAtomIds() 700 (0,) 701 >>> mList[0][1].GetAtomIds() 702 (3,) 703 704 The first feature type, Donor, has a single match: 705 >>> len(mList[1]) 706 1 707 >>> mList[1][0].GetAtomIds() 708 (3,) 709 710 """ 711 return MatchFeatsToMol(mol, featFactory, pcophore.getFeatures())
712
713 -def _getFeatDict(mol,featFactory,features):
714 """ **INTERNAL USE ONLY** 715 716 >>> import os.path 717 >>> dataDir = os.path.join(RDConfig.RDCodeDir,'Chem/Pharm3D/test_data') 718 >>> featFactory = ChemicalFeatures.BuildFeatureFactory(os.path.join(dataDir,'BaseFeatures.fdef')) 719 >>> activeFeats = [ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(0.0, 0.0, 0.0)), 720 ... ChemicalFeatures.FreeChemicalFeature('Donor',Geometry.Point3D(0.0, 0.0, 0.0))] 721 >>> m = Chem.MolFromSmiles('FCCN') 722 >>> d =_getFeatDict(m,featFactory,activeFeats) 723 >>> sorted(list(d.keys())) 724 ['Acceptor', 'Donor'] 725 >>> donors = d['Donor'] 726 >>> len(donors) 727 1 728 >>> donors[0].GetAtomIds() 729 (3,) 730 >>> acceptors = d['Acceptor'] 731 >>> len(acceptors) 732 2 733 >>> acceptors[0].GetAtomIds() 734 (0,) 735 >>> acceptors[1].GetAtomIds() 736 (3,) 737 738 """ 739 molFeats = {} 740 for feat in features: 741 family = feat.GetFamily() 742 if not family in molFeats: 743 matches = featFactory.GetFeaturesForMol(mol,includeOnly=family) 744 molFeats[family] = matches 745 return molFeats
746
747 -def MatchFeatsToMol(mol, featFactory, features):
748 """ generates a list of all possible mappings of each feature to a molecule 749 750 Returns a 2-tuple: 751 1) a boolean indicating whether or not all features were found 752 2) a list, numFeatures long, of sequences of features 753 754 755 >>> import os.path 756 >>> dataDir = os.path.join(RDConfig.RDCodeDir,'Chem/Pharm3D/test_data') 757 >>> featFactory = ChemicalFeatures.BuildFeatureFactory(os.path.join(dataDir,'BaseFeatures.fdef')) 758 >>> activeFeats = [ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(0.0, 0.0, 0.0)), 759 ... ChemicalFeatures.FreeChemicalFeature('Donor',Geometry.Point3D(0.0, 0.0, 0.0))] 760 >>> m = Chem.MolFromSmiles('FCCN') 761 >>> match,mList = MatchFeatsToMol(m,featFactory,activeFeats) 762 >>> match 763 True 764 765 Two feature types: 766 >>> len(mList) 767 2 768 769 The first feature type, Acceptor, has two matches: 770 >>> len(mList[0]) 771 2 772 >>> mList[0][0].GetAtomIds() 773 (0,) 774 >>> mList[0][1].GetAtomIds() 775 (3,) 776 777 The first feature type, Donor, has a single match: 778 >>> len(mList[1]) 779 1 780 >>> mList[1][0].GetAtomIds() 781 (3,) 782 783 """ 784 molFeats = _getFeatDict(mol,featFactory,features) 785 res = [] 786 for feat in features: 787 matches = molFeats.get(feat.GetFamily(),[]) 788 if len(matches) == 0 : 789 return False, None 790 res.append(matches) 791 return True, res
792
793 -def CombiEnum(sequence):
794 """ This generator takes a sequence of sequences as an argument and 795 provides all combinations of the elements of the subsequences: 796 797 >>> gen = CombiEnum(((1,2),(10,20))) 798 >>> next(gen) 799 [1, 10] 800 >>> next(gen) 801 [1, 20] 802 803 >>> [x for x in CombiEnum(((1,2),(10,20)))] 804 [[1, 10], [1, 20], [2, 10], [2, 20]] 805 806 >>> [x for x in CombiEnum(((1,2),(10,20),(100,200)))] 807 [[1, 10, 100], [1, 10, 200], [1, 20, 100], [1, 20, 200], [2, 10, 100], [2, 10, 200], [2, 20, 100], [2, 20, 200]] 808 809 """ 810 if not len(sequence): 811 yield [] 812 elif len(sequence)==1: 813 for entry in sequence[0]: 814 yield [entry] 815 else: 816 for entry in sequence[0]: 817 for subVal in CombiEnum(sequence[1:]): 818 yield [entry]+subVal
819
820 -def DownsampleBoundsMatrix(bm,indices,maxThresh=4.0):
821 """ removes rows from a bounds matrix that are 822 that are greater than a threshold value away from a set of 823 other points 824 825 returns the modfied bounds matrix 826 827 The goal of this function is to remove rows from the bounds matrix 828 that correspond to atoms that are likely to be quite far from 829 the pharmacophore we're interested in. Because the bounds smoothing 830 we eventually have to do is N^3, this can be a big win 831 832 >>> boundsMat = numpy.array([[0.0,3.0,4.0],[2.0,0.0,3.0],[2.0,2.0,0.0]]) 833 >>> bm = DownsampleBoundsMatrix(boundsMat,(0,),3.5) 834 >>> bm.shape == (2, 2) 835 True 836 837 we don't touch the input matrix: 838 >>> boundsMat.shape == (3, 3) 839 True 840 841 >>> print(', '.join(['%.3f'%x for x in bm[0]])) 842 0.000, 3.000 843 >>> print(', '.join(['%.3f'%x for x in bm[1]])) 844 2.000, 0.000 845 846 if the threshold is high enough, we don't do anything: 847 >>> boundsMat = numpy.array([[0.0,4.0,3.0],[2.0,0.0,3.0],[2.0,2.0,0.0]]) 848 >>> bm = DownsampleBoundsMatrix(boundsMat,(0,),5.0) 849 >>> bm.shape == (3, 3) 850 True 851 852 If there's a max value that's close enough to *any* of the indices 853 we pass in, we'll keep it: 854 >>> boundsMat = numpy.array([[0.0,4.0,3.0],[2.0,0.0,3.0],[2.0,2.0,0.0]]) 855 >>> bm = DownsampleBoundsMatrix(boundsMat,(0,1),3.5) 856 >>> bm.shape == (3, 3) 857 True 858 859 """ 860 nPts = bm.shape[0] 861 k = numpy.zeros(nPts,numpy.int0) 862 for idx in indices: k[idx]=1 863 for i in indices: 864 row = bm[i] 865 for j in range(i+1,nPts): 866 if not k[j] and row[j]<maxThresh: 867 k[j]=1 868 keep = numpy.nonzero(k)[0] 869 bm2 = numpy.zeros((len(keep),len(keep)),numpy.float) 870 for i,idx in enumerate(keep): 871 row = bm[idx] 872 bm2[i] = numpy.take(row,keep) 873 return bm2
874
875 -def CoarseScreenPharmacophore(atomMatch,bounds,pcophore,verbose=False):
876 """ 877 878 >>> feats = [ChemicalFeatures.FreeChemicalFeature('HBondAcceptor', 'HAcceptor1', Geometry.Point3D(0.0, 0.0, 0.0)), 879 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 880 ... ChemicalFeatures.FreeChemicalFeature('Aromatic', 'Aromatic1', Geometry.Point3D(5.12, 0.908, 0.0)), 881 ... ] 882 >>> pcophore=Pharmacophore.Pharmacophore(feats) 883 >>> pcophore.setLowerBound(0,1, 1.1) 884 >>> pcophore.setUpperBound(0,1, 1.9) 885 >>> pcophore.setLowerBound(0,2, 2.1) 886 >>> pcophore.setUpperBound(0,2, 2.9) 887 >>> pcophore.setLowerBound(1,2, 2.1) 888 >>> pcophore.setUpperBound(1,2, 3.9) 889 890 >>> bounds = numpy.array([[0,2,3],[1,0,4],[2,3,0]],numpy.float) 891 >>> CoarseScreenPharmacophore(((0,),(1,)),bounds,pcophore) 892 True 893 894 >>> CoarseScreenPharmacophore(((0,),(2,)),bounds,pcophore) 895 False 896 897 >>> CoarseScreenPharmacophore(((1,),(2,)),bounds,pcophore) 898 False 899 900 >>> CoarseScreenPharmacophore(((0,),(1,),(2,)),bounds,pcophore) 901 True 902 903 >>> CoarseScreenPharmacophore(((1,),(0,),(2,)),bounds,pcophore) 904 False 905 906 >>> CoarseScreenPharmacophore(((2,),(1,),(0,)),bounds,pcophore) 907 False 908 909 # we ignore the point locations here and just use their definitions: 910 >>> feats = [ChemicalFeatures.FreeChemicalFeature('HBondAcceptor', 'HAcceptor1', Geometry.Point3D(0.0, 0.0, 0.0)), 911 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 912 ... ChemicalFeatures.FreeChemicalFeature('Aromatic', 'Aromatic1', Geometry.Point3D(5.12, 0.908, 0.0)), 913 ... ChemicalFeatures.FreeChemicalFeature('HBondDonor', 'HDonor1', Geometry.Point3D(2.65, 0.0, 0.0)), 914 ... ] 915 >>> pcophore=Pharmacophore.Pharmacophore(feats) 916 >>> pcophore.setLowerBound(0,1, 2.1) 917 >>> pcophore.setUpperBound(0,1, 2.9) 918 >>> pcophore.setLowerBound(0,2, 2.1) 919 >>> pcophore.setUpperBound(0,2, 2.9) 920 >>> pcophore.setLowerBound(0,3, 2.1) 921 >>> pcophore.setUpperBound(0,3, 2.9) 922 >>> pcophore.setLowerBound(1,2, 1.1) 923 >>> pcophore.setUpperBound(1,2, 1.9) 924 >>> pcophore.setLowerBound(1,3, 1.1) 925 >>> pcophore.setUpperBound(1,3, 1.9) 926 >>> pcophore.setLowerBound(2,3, 1.1) 927 >>> pcophore.setUpperBound(2,3, 1.9) 928 >>> bounds = numpy.array([[0,3,3,3],[2,0,2,2],[2,1,0,2],[2,1,1,0]],numpy.float) 929 930 >>> CoarseScreenPharmacophore(((0,),(1,),(2,),(3,)),bounds,pcophore) 931 True 932 933 >>> CoarseScreenPharmacophore(((0,),(1,),(3,),(2,)),bounds,pcophore) 934 True 935 936 >>> CoarseScreenPharmacophore(((1,),(0,),(3,),(2,)),bounds,pcophore) 937 False 938 939 """ 940 for k in range(len(atomMatch)): 941 if len(atomMatch[k])==1: 942 for l in range(k+1,len(atomMatch)): 943 if len(atomMatch[l])==1: 944 idx0 = atomMatch[k][0] 945 idx1 = atomMatch[l][0] 946 if idx1<idx0: 947 idx0,idx1=idx1,idx0 948 if bounds[idx1,idx0] >= pcophore.getUpperBound(k, l) or \ 949 bounds[idx0,idx1] <= pcophore.getLowerBound(k, l) : 950 if verbose: 951 print('\t (%d,%d) [%d,%d] fail'%(idx1,idx0,k,l)) 952 print('\t %f,%f - %f,%f' % 953 (bounds[idx1,idx0],pcophore.getUpperBound(k,l), 954 bounds[idx0,idx1],pcophore.getLowerBound(k,l))) 955 #logger.debug('\t >%s'%str(atomMatch)) 956 #logger.debug() 957 #logger.debug('\t %f,%f - %f,%f'%(bounds[idx1,idx0],pcophore.getUpperBound(k,l), 958 # bounds[idx0,idx1],pcophore.getLowerBound(k,l))) 959 return False 960 return True
961
962 -def Check2DBounds(atomMatch,mol,pcophore):
963 """ checks to see if a particular mapping of features onto 964 a molecule satisfies a pharmacophore's 2D restrictions 965 966 967 >>> activeFeats = [ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(0.0, 0.0, 0.0)), 968 ... ChemicalFeatures.FreeChemicalFeature('Donor',Geometry.Point3D(0.0, 0.0, 0.0))] 969 >>> pcophore= Pharmacophore.Pharmacophore(activeFeats) 970 >>> pcophore.setUpperBound2D(0,1,3) 971 >>> m = Chem.MolFromSmiles('FCC(N)CN') 972 >>> Check2DBounds(((0,),(3,)),m,pcophore) 973 True 974 >>> Check2DBounds(((0,),(5,)),m,pcophore) 975 False 976 977 """ 978 dm = Chem.GetDistanceMatrix(mol,False,False,False) 979 nFeats = len(atomMatch) 980 for i in range(nFeats): 981 for j in range(i+1,nFeats): 982 lowerB = pcophore._boundsMat2D[j,i] #lowerB = pcophore.getLowerBound2D(i,j) 983 upperB = pcophore._boundsMat2D[i,j] #upperB = pcophore.getUpperBound2D(i,j) 984 dij=10000 985 for atomI in atomMatch[i]: 986 for atomJ in atomMatch[j]: 987 try: 988 dij = min(dij,dm[atomI,atomJ]) 989 except IndexError: 990 print('bad indices:',atomI,atomJ) 991 print(' shape:',dm.shape) 992 print(' match:',atomMatch) 993 print(' mol:') 994 print(Chem.MolToMolBlock(mol)) 995 raise IndexError 996 if dij<lowerB or dij>upperB: 997 return False 998 return True
999
1000 -def _checkMatch(match,mol,bounds,pcophore,use2DLimits):
1001 """ **INTERNAL USE ONLY** 1002 1003 checks whether a particular atom match can be satisfied by 1004 a molecule 1005 1006 """ 1007 atomMatch = ChemicalFeatures.GetAtomMatch(match) 1008 if not atomMatch: 1009 return None 1010 elif use2DLimits: 1011 if not Check2DBounds(atomMatch,mol,pcophore): 1012 return None 1013 if not CoarseScreenPharmacophore(atomMatch,bounds,pcophore): 1014 return None 1015 return atomMatch
1016
1017 -def ConstrainedEnum(matches,mol,pcophore,bounds,use2DLimits=False, 1018 index=0,soFar=[]):
1019 """ Enumerates the list of atom mappings a molecule 1020 has to a particular pharmacophore. 1021 We do check distance bounds here. 1022 1023 1024 """ 1025 nMatches = len(matches) 1026 if index>=nMatches: 1027 yield soFar,[] 1028 elif index==nMatches-1: 1029 for entry in matches[index]: 1030 nextStep = soFar+[entry] 1031 if index != 0: 1032 atomMatch = _checkMatch(nextStep,mol,bounds,pcophore,use2DLimits) 1033 else: 1034 atomMatch = ChemicalFeatures.GetAtomMatch(nextStep) 1035 if atomMatch: 1036 yield soFar+[entry],atomMatch 1037 else: 1038 for entry in matches[index]: 1039 nextStep = soFar+[entry] 1040 if index != 0: 1041 atomMatch = _checkMatch(nextStep,mol,bounds,pcophore,use2DLimits) 1042 if not atomMatch: 1043 continue 1044 for val in ConstrainedEnum(matches,mol,pcophore,bounds,use2DLimits=use2DLimits, 1045 index=index+1,soFar=nextStep): 1046 if val: 1047 yield val
1048
1049 -def MatchPharmacophore(matches,bounds,pcophore,useDownsampling=False, 1050 use2DLimits=False,mol=None,excludedVolumes=None, 1051 useDirs=False):
1052 """ 1053 1054 if use2DLimits is set, the molecule must also be provided and topological 1055 distances will also be used to filter out matches 1056 1057 """ 1058 for match,atomMatch in ConstrainedEnum(matches,mol,pcophore,bounds, 1059 use2DLimits=use2DLimits): 1060 bm = bounds.copy() 1061 bm = UpdatePharmacophoreBounds(bm,atomMatch,pcophore,useDirs=useDirs,mol=mol); 1062 1063 if excludedVolumes: 1064 localEvs = [] 1065 for eV in excludedVolumes: 1066 featInfo = [] 1067 for i,entry in enumerate(atomMatch): 1068 info = list(eV.featInfo[i]) 1069 info[0] = entry 1070 featInfo.append(info) 1071 localEvs.append(ExcludedVolume.ExcludedVolume(featInfo,eV.index, 1072 eV.exclusionDist)) 1073 bm = AddExcludedVolumes(bm,localEvs,smoothIt=False) 1074 1075 sz = bm.shape[0] 1076 if useDownsampling: 1077 indices = [] 1078 for entry in atomMatch: 1079 indices.extend(entry) 1080 if excludedVolumes: 1081 for vol in localEvs: 1082 indices.append(vol.index) 1083 bm = DownsampleBoundsMatrix(bm,indices) 1084 if DG.DoTriangleSmoothing(bm): 1085 return 0,bm,match,(sz,bm.shape[0]) 1086 1087 return 1,None,None,None
1088
1089 -def GetAllPharmacophoreMatches(matches,bounds,pcophore,useDownsampling=0, 1090 progressCallback=None, 1091 use2DLimits=False,mol=None, 1092 verbose=False):
1093 res = [] 1094 nDone = 0 1095 for match in CombiEnum(matches): 1096 atomMatch = ChemicalFeatures.GetAtomMatch(match) 1097 if atomMatch and use2DLimits and mol: 1098 pass2D = Check2DBounds(atomMatch,mol,pcophore) 1099 if verbose: 1100 print('..',atomMatch) 1101 print(' ..Pass2d:',pass2D) 1102 else: 1103 pass2D = True 1104 if atomMatch and pass2D and \ 1105 CoarseScreenPharmacophore(atomMatch,bounds,pcophore,verbose=verbose): 1106 if verbose: 1107 print(' ..CoarseScreen: Pass') 1108 1109 bm = bounds.copy() 1110 if verbose: 1111 print('pre update:') 1112 for row in bm: 1113 print(' ',' '.join(['% 4.2f'%x for x in row])) 1114 bm = UpdatePharmacophoreBounds(bm,atomMatch,pcophore); 1115 sz = bm.shape[0] 1116 if verbose: 1117 print('pre downsample:') 1118 for row in bm: 1119 print(' ',' '.join(['% 4.2f'%x for x in row])) 1120 1121 if useDownsampling: 1122 indices = [] 1123 for entry in atomMatch: 1124 indices += list(entry) 1125 bm = DownsampleBoundsMatrix(bm,indices) 1126 if verbose: 1127 print('post downsample:') 1128 for row in bm: 1129 print(' ',' '.join(['% 4.2f'%x for x in row])) 1130 1131 if DG.DoTriangleSmoothing(bm): 1132 res.append(match) 1133 elif verbose: 1134 print('cannot smooth') 1135 nDone+=1 1136 if progressCallback: 1137 progressCallback(nDone) 1138 return res
1139 1140
1141 -def ComputeChiralVolume(mol,centerIdx,confId=-1):
1142 """ Computes the chiral volume of an atom 1143 1144 We're using the chiral volume formula from Figure 7 of 1145 Blaney and Dixon, Rev. Comp. Chem. V, 299-335 (1994) 1146 1147 >>> import os.path 1148 >>> dataDir = os.path.join(RDConfig.RDCodeDir,'Chem/Pharm3D/test_data') 1149 1150 R configuration atoms give negative volumes: 1151 >>> mol = Chem.MolFromMolFile(os.path.join(dataDir,'mol-r.mol')) 1152 >>> Chem.AssignStereochemistry(mol) 1153 >>> mol.GetAtomWithIdx(1).GetProp('_CIPCode') 1154 'R' 1155 >>> ComputeChiralVolume(mol,1) < 0 1156 True 1157 1158 S configuration atoms give positive volumes: 1159 >>> mol = Chem.MolFromMolFile(os.path.join(dataDir,'mol-s.mol')) 1160 >>> Chem.AssignStereochemistry(mol) 1161 >>> mol.GetAtomWithIdx(1).GetProp('_CIPCode') 1162 'S' 1163 >>> ComputeChiralVolume(mol,1) > 0 1164 True 1165 1166 Non-chiral (or non-specified) atoms give zero volume: 1167 >>> ComputeChiralVolume(mol,0) == 0.0 1168 True 1169 1170 We also work on 3-coordinate atoms (with implicit Hs): 1171 >>> mol = Chem.MolFromMolFile(os.path.join(dataDir,'mol-r-3.mol')) 1172 >>> Chem.AssignStereochemistry(mol) 1173 >>> mol.GetAtomWithIdx(1).GetProp('_CIPCode') 1174 'R' 1175 >>> ComputeChiralVolume(mol,1)<0 1176 True 1177 1178 >>> mol = Chem.MolFromMolFile(os.path.join(dataDir,'mol-s-3.mol')) 1179 >>> Chem.AssignStereochemistry(mol) 1180 >>> mol.GetAtomWithIdx(1).GetProp('_CIPCode') 1181 'S' 1182 >>> ComputeChiralVolume(mol,1)>0 1183 True 1184 1185 1186 1187 """ 1188 conf = mol.GetConformer(confId) 1189 Chem.AssignStereochemistry(mol) 1190 center = mol.GetAtomWithIdx(centerIdx) 1191 if not center.HasProp('_CIPCode'): 1192 return 0.0 1193 1194 nbrs = center.GetNeighbors() 1195 nbrRanks = [] 1196 for nbr in nbrs: 1197 rank = int(nbr.GetProp('_CIPRank')) 1198 pos = conf.GetAtomPosition(nbr.GetIdx()) 1199 nbrRanks.append((rank,pos)) 1200 1201 # if we only have three neighbors (i.e. the determining H isn't present) 1202 # then use the central atom as the fourth point: 1203 if len(nbrRanks)==3: 1204 nbrRanks.append((-1,conf.GetAtomPosition(centerIdx))) 1205 nbrRanks.sort() 1206 1207 ps = [x[1] for x in nbrRanks] 1208 v1 = ps[0]-ps[3] 1209 v2 = ps[1]-ps[3] 1210 v3 = ps[2]-ps[3] 1211 1212 res = v1.DotProduct(v2.CrossProduct(v3)) 1213 return res
1214 1215 1216 1217 1218 1219 1220 #------------------------------------ 1221 # 1222 # doctest boilerplate 1223 #
1224 -def _test():
1225 import doctest,sys 1226 return doctest.testmod(sys.modules["__main__"])
1227 1228 if __name__ == '__main__': 1229 import sys 1230 failed,tried = _test() 1231 sys.exit(failed) 1232