Source Code for Module MolKit.hydrogenBondBuilder

  1  ############################################################################ 
  2  # 
  3  # Author:  Ruth Huey 
  4  # 
  5  # Copyright: M. Sanner TSRI 2005 
  6  # 
  7  ############################################################################# 
  8   
  9  # 
 10  # $Header: /opt/cvs/python/packages/share1.5/MolKit/hydrogenBondBuilder.py,v 1.4 2007/02/20 21:20:28 rhuey Exp $ 
 11  # 
 12  # $Id: hydrogenBondBuilder.py,v 1.4 2007/02/20 21:20:28 rhuey Exp $ 
 13  # 
 14   
 15  """ 
 16  This module implements the HydrogenBondBuilder class which builds hydrogen 
 17  bonds between appropriate atoms. 
 18   
 19  """ 
 20   
 21  #need these in order to typeAtoms 
 22  from PyBabel.atomTypes import AtomHybridization 
 23  from PyBabel.bo import BondOrder 
 24   
 25  from MolKit.molecule import Atom, AtomSet, HydrogenBond 
 26  from MolKit.distanceSelector import DistanceSelector 
 27   
 28  #needed for some math.. 
 29  import Numeric, math 
 30   
 31  #set up lists of types for donors/acceptors 
 32  sp2Donors = ['Nam', 'Ng+', 'Npl'] 
 33  sp3Donors = ['N3+', 'O3','S3'] 
 34  allDonors = ['Nam', 'Ng+', 'Npl', 'N3+', 'O3','S3'] 
 35  #added Nam because of bdna cytidine N3 
 36  #NB: Npl cannot be an acceptor if the sum of its bonds' bondOrder>2 
 37  sp2Acceptors = ['O2', 'O-', 'Npl', 'Nam'] 
 38  sp3Acceptors = ['S3', 'O3'] 
 39  allAcceptors = ['O2', 'O-', 'Npl', 'Nam', 'S3', 'O3'] 
 40  #???? 
 41  #for n in allAcceptors: 
 42  #    allAcceptors.append('a'+n) 
 43   
 44   
 45 -def dist(c1, c2): 
 46      d = Numeric.array(c2) - Numeric.array(c1)  
 47      ans = math.sqrt(Numeric.sum(d*d)) 
 48      return round(ans, 3) 
 49   
 50   
 51 -def getAngle(ac, hat, don ): 
 52      #acCoords = getTransformedCoords(ac) 
 53      #hCoords = getTransformedCoords(hat) 
 54      #dCoords = getTransformedCoords(don) 
 55      acCoords = ac.coords 
 56      hCoords = hat.coords 
 57      dCoords = don.coords 
 58      pt1 = Numeric.array(acCoords, 'f') 
 59      pt2 = Numeric.array(hCoords, 'f') 
 60      pt3 = Numeric.array(dCoords, 'f') 
 61      #pt1 = Numeric.array(ac.coords, 'f') 
 62      #pt2 = Numeric.array(hat.coords, 'f') 
 63      #pt3 = Numeric.array(don.coords, 'f') 
 64      v1 = Numeric.array(pt1 - pt2) 
 65      v2 = Numeric.array(pt3 - pt2) 
 66      dist1 = math.sqrt(Numeric.sum(v1*v1)) 
 67      dist2 = math.sqrt(Numeric.sum(v2*v2)) 
 68      sca = Numeric.dot(v1, v2)/(dist1*dist2) 
 69      if sca>1.0: 
 70          sca = 1.0 
 71      elif sca<-1.0: 
 72          sca = -1.0 
 73      ang =  math.acos(sca)*180./math.pi 
 74      return round(ang, 5) 
 75   
 76   
 77 -def applyTransformation(pt, mat): 
 78      pth = [pt[0], pt[1], pt[2], 1.0] 
 79      return Numeric.matrixmultiply(mat, pth)[:3] 
 80   
 81   
 82 -def getTransformedCoords(atom): 
 83      # when there is no viewer, the geomContainer is None 
 84      if atom.top.geomContainer is None: 
 85          return atom.coords 
 86      g = atom.top.geomContainer.geoms['master'] 
 87      c = applyTransformation(atom.coords, g.GetMatrix(g)) 
 88      return  c.astype('f') 
 89   
 90   
 91  distCutoff=2.25     #H-Acc Distance cutoff 
 92  distCutoff2=3.00    #Donor-Acc Distance cutoff 
 93  d2min=120           #sp2-donor-hydrogen-acceptor angle  
 94  d2max=180           # 
 95  d3min=120           #sp3-donor-hydrogen-acceptor angle 
 96  d3max=180           # 
 97  a2min=110           #sp2-donor-acceptor-acceptorN angle,  
 98  a2max=150           #    where acceptorN is 'neighbor' atom bonded to acceptor 
 99  a3min=100           #sp3-donor-acceptor-acceptorN angle 
100  a3max=150           #    where acceptorN is 'neighbor' atom bonded to acceptor 
101   
102   
103 -class HydrogenBondBuilder: 
104      """ 
105      object which can build hydrogen bonds between atoms according 
106      to their coords and atom type 
107      """ 
108   
109 -    def __init__(self, distCutoff=distCutoff, distCutoff2=distCutoff2, 
110                          d2min=d2min, d2max=d2max, d3min=d3min, d3max=d3max, 
111                          a2min=a2min, a2max=a2max, a3min=a3min, a3max=a3max, 
112                          donorTypes=allDonors, acceptorTypes=allAcceptors, 
113                          distOnly=False): 
114          d = self.paramDict = {} 
115          d['distCutoff'] = distCutoff 
116          d['distCutoff2'] = distCutoff2 
117          d['d2min'] = d2min 
118          d['d2max'] = d2max 
119          d['d3min'] = d3min 
120          d['d3max'] = d3max 
121          d['a2min'] = a2min 
122          d['a2max'] = a2max 
123          d['a3min'] = a3min 
124          d['a3max'] = a3max 
125          d['donorTypes'] = donorTypes 
126          d['acceptorTypes'] = acceptorTypes 
127          d['distOnly'] = distOnly 
128          self.distSelector = DistanceSelector(return_dist=0) 
129   
130   
131 -    def check_babel_types(self, ats): 
132          try: 
133              ats.babel_type 
134              ats._bndtyped 
135          except AttributeError: 
136              babel = AtomHybridization() 
137              bond_orderer = BondOrder() 
138              tops = ats.top.uniq() 
139              for mol in tops: 
140                  babel.assignHybridization(mol.allAtoms) 
141                  bond_orderer.assignBondOrder(mol.allAtoms, mol.allAtoms.bonds[0]) 
142                  mol.allAtoms._bndtyped = 1 
143   
144   
145 -    def reset(self, ats): 
146          tops = ats.top.uniq() 
147          for mol in tops: 
148              for a in mol.allAtoms: 
149                  if hasattr(a, 'hbonds'): 
150                      for item in a.hbonds: 
151                          del item 
152                      delattr(a, 'hbonds') 
153   
154   
155 -    def build(self, group1, group2=None, reset=True, paramDict=None): 
156          """atDict <- build(group1, group2, reset, paramDict=None, **kw): 
157              group1: atoms  
158              group2: atoms  
159              reset: remove all previous hbonds, default True! 
160   
161              paramDict: a dictionary with these keys and default values 
162              distCutoff: 2.25  hydrogen--acceptor distance 
163              distCutoff2: 3.00 donor... acceptor distance 
164              d2min: 120 <min theta for sp2 hybridized donors> 
165              d2max: 180 <max theta for sp2 hybridized donors> 
166              d3min: 120 <min theta for sp3 hybridized donors> 
167              d3max: 170 <max theta for sp3 hybridized donors> 
168   
169              a2min: 120 <min phi for sp2 hybridized donors> 
170              a2max: 150 <max phi for sp2 hybridized donors> 
171              a3min: 100 <min phi for sp3 hybridized donors> 
172              a3max: 150 <max phi for sp3 hybridized donors> 
173              @@FIX THIS: these do not seem to be input here 
174              donorTypes = allDonors 
175              acceptorTypes = allAcceptors 
176             """ 
177          #setup parameter dictionary 
178          if paramDict is None: 
179              paramDict = self.paramDict 
180   
181          #setup group2  
182          if group2 is None: 
183              group2 = group1 
184   
185          #process each group 
186          #group1 
187          if group1.__class__!=Atom: 
188              group1 = group1.findType(Atom) 
189              #print "now group1=", group1.full_name() 
190              if not len(group1): 
191                  return "ERROR"  #@@ OR WHAT? 
192          self.check_babel_types(group1) 
193   
194          #group2 
195          if group2.__class__!=Atom: 
196              group2 = group2.findType(Atom) 
197              #print "now group2=", group2.full_name() 
198              if not len(group2): 
199                  return "ERROR"  #@@ OR WHAT? 
200          self.check_babel_types(group2) 
201   
202          if reset: 
203              #do optional reset: remove all prior hbonds  
204              self.reset(group1) 
205              self.reset(group2) 
206   
207          #print "group1=", len(group1) 
208          #print "group2=", len(group2) 
209   
210          #buildHbonds 
211          atDict = {} 
212          dict1 = self.buildD(group1, paramDict) 
213   
214          #@@what is this doing here??? 
215          #sp2 hybridized atoms 
216          #dAts2 = ats.get(lambda x, l=sp2: x.babel_type in l) 
217   
218          if group1==group2: 
219              #only one step: 
220              atD1 = self.process(dict1, dict1, paramDict) 
221              #print 'len(atD1).keys()=', len(atD1.keys()) 
222          else: 
223              # two steps: 
224              # 1: group1 donors v group2 acceptors  
225              # 2: group1 acceptors vs group2 donors 
226              dict2 = self.buildD(group2, paramDict) 
227              atD1 = self.process(dict1, dict2, paramDict) 
228              atD2 = self.process(dict2, dict1, paramDict) 
229              #print 'len(atD1).keys()=', len(atD1.keys()) 
230              #print 'len(atD2).keys()=', len(atD2.keys()) 
231          #if called with 1 atom could get tuple of two empty dictionaries 
232          if type(atD1)==type(atDict): 
233              atDict.update(atD1) 
234              if group1!=group2 and type(atD2)==type(atDict): 
235                  atDict.update(atD2) 
236          #@@DESCRIBE atDict 
237          #this dict has atomSets as keys and as values(?check that) 
238          return atDict 
239   
240   
241 -    def checkForPossibleH(self, ats, blen): 
242          #@@FIX THIS: WHAT IS THE POINT OF THIS??? 
243          #check that if at has all bonds, at least one is to a hydrogen 
244          # have to do this by element?? 
245          probAts = AtomSet(ats.get(lambda x, blen=blen: len(x.bonds)==blen)) 
246          #probOAts = ats.get(lambda x, blen=blen: len(x.bonds)==blen) 
247          #probSAts = ats.get(lambda x, blen=blen: len(x.bonds)==blen) 
248          if probAts: 
249              rAts = AtomSet([]) 
250              for at in probAts: 
251                  if not len(at.findHydrogens()): 
252                      rAts.append(at) 
253              if len(rAts): 
254                  ats =  ats.subtract(rAts) 
255          return ats 
256       
257   
258 -    def getHBDonors(self, ats, donorList): 
259          #getHBDonors  
260          sp2 = [] 
261          sp3 = [] 
262          for item in sp2Donors: 
263              if item in donorList: sp2.append(item) 
264          for item in sp3Donors: 
265              if item in donorList: sp3.append(item) 
266           
267          dAts2 = ats.get(lambda x, l=sp2: x.babel_type in l) 
268          if not dAts2: dAts2=AtomSet([]) 
269          else:  
270              dAts2 = self.checkForPossibleH(dAts2, 3) 
271           
272          dAts3 = ats.get(lambda x, l=sp3: x.babel_type in l) 
273          if not dAts3: dAts3=AtomSet([]) 
274          else: dAts3 = self.checkForPossibleH(dAts3, 4) 
275          return dAts2, dAts3 
276   
277   
278 -    def filterAcceptors(self, accAts): 
279          ntypes = ['Npl', 'Nam'] 
280          npls = accAts.get(lambda x, ntypes=ntypes: x.babel_type=='Npl') 
281          nams = accAts.get(lambda x, ntypes=ntypes: x.babel_type=='Nam') 
282          #nAts = accAts.get(lambda x, ntypes=ntypes: x.babel_type in ntypes) 
283          restAts = accAts.get(lambda x, ntypes=ntypes: x.babel_type not in ntypes) 
284          if not restAts: restAts = AtomSet([]) 
285          #if nAts: 
286          if npls: 
287              #for at in nAts: 
288              for at in npls: 
289                  s = 0 
290                  for b in at.bonds: 
291                      if b.bondOrder=='aromatic': 
292                          s = s + 2 
293                      else: s = s + b.bondOrder 
294                  #if s<3: 
295                  #apparently this is wrong 
296                  if s<4: 
297                      restAts.append(at) 
298          if nams: 
299              #for at in nAts: 
300              for at in nams: 
301                  s = 0 
302                  for b in at.bonds: 
303                      if b.bondOrder=='aromatic': 
304                          s = s + 2 
305                      else: s = s + b.bondOrder 
306                      #s = s + b.bondOrder 
307                  if s<3: 
308                      restAts.append(at) 
309          return restAts 
310                   
311   
312 -    def getHBAcceptors(self, ats, acceptorList): 
313          #print "getHBAcceptors: acceptorList=", acceptorList 
314          #getHBAcceptors  
315          sp2 = [] 
316          sp3 = [] 
317          for item in sp2Acceptors: 
318              if item in acceptorList: sp2.append(item) 
319          for item in sp3Acceptors: 
320              if item in acceptorList: sp3.append(item) 
321   
322          dAts2 = AtomSet(ats.get(lambda x, l=sp2: x.babel_type in l)) 
323          if dAts2:  
324              dAts2 = self.filterAcceptors(dAts2) 
325          dAts3 = AtomSet(ats.get(lambda x, l=sp3: x.babel_type in l)) 
326          return dAts2, dAts3 
327         
328   
329 -    def buildD(self, ats, paramDict=None): 
330          if paramDict is None: 
331              paramDict = self.paramDict 
332          #these are from the __call__ method of vf.buildHBonds 
333          if not paramDict.has_key('distCutoff'): 
334              paramDict['distCutoff'] = 2.25 
335          if not paramDict.has_key('distCutoff2'): 
336              paramDict['distCutoff2'] = 3.00 
337          if not paramDict.has_key('d2min'): 
338              paramDict['d2min'] = 120. 
339          if not paramDict.has_key('d2max'): 
340              paramDict['d2max'] = 180. 
341          if not paramDict.has_key('d3min'): 
342              paramDict['d3min'] = 120. 
343          if not paramDict.has_key('d3max'): 
344              paramDict['d3max'] = 170. 
345          if not paramDict.has_key('a2min'): 
346              paramDict['a2min'] = 130. 
347          if not paramDict.has_key('a2max'): 
348              paramDict['a2max'] = 170. 
349          if not paramDict.has_key('a3min'): 
350              paramDict['a3min'] = 120. 
351          if not paramDict.has_key('a3max'): 
352              paramDict['a3max'] = 170. 
353          if not paramDict.has_key('distOnly'): 
354              paramDict['distOnly'] = 0 
355          if not paramDict.has_key('donorTypes'): 
356              paramDict['donorTypes'] = allDonors 
357          if not paramDict.has_key('acceptorTypes'): 
358              paramDict['acceptorTypes'] = allAcceptors 
359   
360          d = {} 
361          donorTypes = paramDict['donorTypes'] 
362          donor2Ats, donor3Ats = self.getHBDonors(ats, donorTypes) 
363          d23 = donor2Ats + donor3Ats 
364          #hAts = ats.get(lambda x, d23=d23: x.element=='H' \ 
365                      #and x.bonds[0].neighborAtom(x) in d23) 
366          hydrogen_atoms = ats.get(lambda x: x.element=='H' and len(x.bonds)) 
367          #hAts = AtomSet(ats.get(lambda x, donorTypes=donorTypes: x.element=='H' \ 
368          hAts = AtomSet(hydrogen_atoms.get(lambda x, donorTypes=donorTypes:  
369                      x.bonds[0].atom1.babel_type in donorTypes\ 
370                      or x.bonds[0].atom2.babel_type in donorTypes)) 
371          d['hAts'] = hAts 
372          d['donor2Ats'] = donor2Ats 
373          d['donor3Ats'] = donor3Ats 
374       
375          acceptorTypes = paramDict['acceptorTypes'] 
376          #print "about to call getHBAcceptors with acceptorTypes=", acceptorTypes 
377          acceptor2Ats, acceptor3Ats = self.getHBAcceptors(ats, acceptorTypes) 
378          d['acceptor2Ats'] = acceptor2Ats 
379          d['acceptor3Ats'] = acceptor3Ats 
380          if acceptor2Ats: 
381              acceptorAts = acceptor2Ats 
382              if acceptor3Ats: 
383                  acceptorAts = acceptorAts + acceptor3Ats 
384          elif acceptor3Ats: 
385              acceptorAts = acceptor3Ats 
386          else: 
387              #CHECK THIS: should it be None or AtomSet([]) 
388              acceptorAts = None 
389          d['acceptorAts'] = acceptorAts 
390          return d 
391               
392   
393 -    def getMat(self, ats): 
394          pass 
395          #tops = ats.top.uniq() 
396          #if len(tops)>1:  
397          #    self.warningMsg('transformation mat=None:>1 mol in atomset!') 
398          #    return None 
399          #g = tops[0].geomContainer.geoms['master'] 
400          #return g.GetMatrix(g) 
401   
402 -    def process(self, dict1, dict2, paramDict): 
403          #hAts are keys, aceptorAts are checks 
404          hAts = dict1['hAts'] 
405          tAts = hAts 
406          dist = paramDict['distCutoff'] 
407          distOnly = paramDict['distOnly'] 
408   
409          if not hAts: 
410              #then use donors and a different distance 
411              tAts = dict1['donor2Ats'] + dict1['donor3Ats'] 
412              dist = paramDict['distCutoff2'] 
413               
414          acceptorAts = dict2['acceptorAts'] 
415          #print "acceptorAts=", acceptorAts 
416          if not acceptorAts or not tAts: #6/14/2004 
417              return {}, {} 
418   
419          #call distanceSelector on two groups of atoms with dist 
420          #keyMat = self.getMat(tAts) 
421          #checkMat = self.getMat(acceptorAts) 
422          atDict = self.distSelector.select(tAts, acceptorAts, dist) 
423          #    keyMat=keyMat, checkMat=checkMat) 
424          #atDict = self.distSelector.select(tAts, acceptorAts, dist) 
425          #first remove bonded angles 
426          atDict = self.removeNeighbors(atDict) 
427   
428          donor2Ats = dict1['donor2Ats'] 
429          donor3Ats = dict1['donor3Ats'] 
430          acceptor2Ats = dict2['acceptor2Ats'] 
431          acceptor3Ats = dict2['acceptor3Ats'] 
432   
433          if distOnly:  
434              #need to build hbonds and return dictionary 
435              self.makeBonds(atDict, donor2Ats, donor3Ats, \ 
436                      acceptor2Ats, acceptor3Ats, paramDict) 
437              return atDict 
438   
439          badAtDict = self.filterBasedOnAngs(atDict, donor2Ats, donor3Ats, \ 
440                      acceptor2Ats, acceptor3Ats, paramDict) 
441          atDict = self.removeBadAts(atDict, badAtDict) 
442          if atDict is None: 
443              atDict = {} 
444          return atDict 
445   
446   
447 -    def makeBonds(self, pD, d2Ats, d3Ats, a2Ats, a3ats, paramDict): 
448          for k in pD.keys(): 
449              if k.element=='H': 
450                  if hasattr(k, 'hbonds') and len(k.hbonds): 
451                      continue 
452                  d = k.bonds[0].atom1 
453                  if id(d)==id(k): d = k.bonds[0].atom2 
454                  #d = k.bonds[0].neighborAtom(k) 
455                  h = k 
456              else:  
457                  d = k 
458                  h = None 
459              #pD[k] is a list of close-enough ats 
460              for ac in pD[k]: 
461                  if ac==d: continue 
462                  dSp2 = d in d2Ats 
463                  aSp2 = ac in a2Ats 
464                  if dSp2: 
465                      if aSp2: typ = 22 
466                      else: typ = 23 
467                  elif aSp2: typ = 32 
468                  else: typ = 33 
469                  #THEY could be already bonded 
470                  alreadyBonded = 0 
471                  if hasattr(d, 'hbonds') and hasattr(ac,'hbonds'): 
472                      for hb in d.hbonds: 
473                          if hb.donAt==ac or hb.accAt==ac: 
474                              alreadyBonded = 1 
475                               
476                  if not alreadyBonded: 
477                      newHB = HydrogenBond(d, ac, h, typ=typ) 
478                      if not hasattr(ac, 'hbonds'): 
479                          ac.hbonds=[] 
480                      if not hasattr(d, 'hbonds'): 
481                          d.hbonds=[] 
482                      ac.hbonds.append(newHB) 
483                      d.hbonds.append(newHB) 
484                      if h is not None: 
485                          #hydrogens can have only 1 hbond 
486                          h.hbonds = [newHB] 
487   
488   
489 -    def filterBasedOnAngs(self, pD, d2Ats, d3Ats, a2Ats, a3ats, paramDict): 
490          badAtDict = {} 
491          d2max = paramDict['d2max'] 
492          d2min = paramDict['d2min'] 
493          d3max = paramDict['d3max'] 
494          d3min = paramDict['d3min'] 
495          #NEED these parameters 
496          a2max = paramDict['a2max'] 
497          a2min = paramDict['a2min'] 
498          a3max = paramDict['a3max'] 
499          a3min = paramDict['a3min'] 
500          #NB now pD keys could be hydrogens OR donors 
501          for k in pD.keys(): 
502              if k.element=='H': 
503                  d = k.bonds[0].atom1 
504                  if id(d)==id(k): d = k.bonds[0].atom2 
505                  #d = k.bonds[0].neighborAtom(k) 
506                  h = k 
507              else:  
508                  d = k 
509                  h = None 
510              badAts = AtomSet([]) 
511              ct = 0 
512              for ac in pD[k]: 
513                  if h is not None: 
514                      ang = getAngle(ac, h, d) 
515                  else: 
516                      acN = ac.bonds[0].atom1 
517                      if id(acN) == id(ac): acN = ac.bonds[0].atom2 
518                      #acN = ac.bonds[0].neighborAtom(ac) 
519                      ang = getAngle(d, ac, acN) 
520                  #print 'ang=', ang 
521                  dSp2 = d in d2Ats 
522                  aSp2 = ac in a2Ats 
523                  #these limits could be adjustable 
524                  if h is not None: 
525                      if dSp2: 
526                          upperLim = d2max 
527                          lowerLim = d2min 
528                          #upperLim = 170 
529                          #lowerLim = 130 
530                      else: 
531                          upperLim = d3max 
532                          lowerLim = d3min 
533                          #upperLim = 180 
534                          #lowerLim = 120 
535                  else: 
536                      #if there is no hydrogen use d-ac-acN angles 
537                      if dSp2: 
538                          upperLim = a2max 
539                          lowerLim = a2min 
540                          #upperLim = 150 
541                          #lowerLim = 110 
542                      else: 
543                          upperLim = a3max 
544                          lowerLim = a3min 
545                          #upperLim = 150 
546                          #lowerLim = 100 
547                  if ang>lowerLim and ang <upperLim: 
548                      #AT THIS MOMENT BUILD HYDROGEN BOND: 
549                      if dSp2: 
550                          if aSp2: typ = 22 
551                          else: typ = 23 
552                      elif aSp2: typ = 32 
553                      else: typ = 33 
554                      #THEY could be already bonded 
555                      alreadyBonded = 0 
556                      if hasattr(d, 'hbonds') and hasattr(ac,'hbonds'): 
557                          for hb in d.hbonds: 
558                              if hb.donAt==ac or hb.accAt==ac: 
559                                  alreadyBonded = 1 
560                      if not alreadyBonded: 
561                          newHB = HydrogenBond(d, ac, h, theta=ang, typ=typ) 
562                          if not hasattr(ac, 'hbonds'): 
563                              ac.hbonds=[] 
564                          if not hasattr(d, 'hbonds'): 
565                              d.hbonds=[] 
566                          ac.hbonds.append(newHB) 
567                          d.hbonds.append(newHB) 
568                          if h is not None: 
569                              #hydrogens can have only 1 hbond 
570                              h.hbonds = [newHB] 
571                          #    newHB.hlen = dist 
572                          #else: 
573                          #    newHB.dlen = dist 
574                  else: 
575                      badAts.append(ac) 
576                  ct = ct + 1 
577              badAtDict[k] = badAts 
578          return badAtDict 
579   
580   
581 -    def removeBadAts(self, atDict, badAtDict): 
582          #clean-up function called after filtering on angles 
583          badKeys= badAtDict.keys() 
584          for at in atDict.keys(): 
585              if at not in badKeys: 
586                  continue 
587              if not len(badAtDict[at]): 
588                  continue 
589              closeAts = atDict[at] 
590              badAts = badAtDict[at] 
591              goodAts = [] 
592              for i in range(len(closeAts)): 
593                  cAt = closeAts[i] 
594                  if cAt not in badAts: 
595                      goodAts.append(cAt) 
596              if len(goodAts): 
597                  atDict[at] = goodAts 
598              else: 
599                  del atDict[at] 
600          return atDict 
601   
602   
603 -    def removeNeighbors(self, atDict): 
604          #filter out at-itself and at-bondedat up to 1:4 
605          #NB keys could be hydrogens OR donors 
606          for at in atDict.keys(): 
607              closeAts = atDict[at] 
608              bondedAts = AtomSet([]) 
609              for b in at.bonds: 
610                  ###at2 = b.neighborAtom(at) 
611                  at2 = b.atom1 
612                  if id(at2)==id(at): at2 = b.atom2 
613                  bondedAts.append(at2) 
614                  #9/13 remove this: 
615                  ##also remove 1-3 
616                  for b2 in at2.bonds: 
617                      at3 = b2.atom1 
618                      if id(at3)==id(at2): at3 = b.atom2 
619                      #at3 = b2.neighborAtom(at2) 
620                      if id(at3)!=id(at): 
621                          bondedAts.append(at3) 
622                      #for b3 in at3.bonds: 
623                          #at4 = b2.neighborAtom(at3) 
624                          #if at4!=at and at4!=at2: 
625                              #bondedAts.append(at4) 
626              bondedAts = bondedAts.uniq() 
627              goodAts = [] 
628              for i in range(len(closeAts)): 
629                  cAt = closeAts[i] 
630                  if cAt not in bondedAts: 
631                      goodAts.append(cAt) 
632              if len(goodAts): 
633                  atDict[at] = goodAts 
634              else: 
635                  del atDict[at] 
636          return atDict 
637   
638   
639 -    def getDonors(self, nodes, paramDict): 
640          donorList = paramDict['donorTypes'] 
641          #print 'donorList=', donorList 
642          # currently this is a set of hydrogens 
643          hats = AtomSet(nodes.get(lambda x: x.element=='H')) 
644          #hats are optional: if none, process donors 
645          # if there are hats: dAts are all atoms bonded to all hydrogens 
646          if hats: 
647              dAts = AtomSet([]) 
648              for at in hats: 
649                  for b in at.bonds: 
650                      at2 = b.atom1 
651                      if id(at2)==id(at): at2 = b.atom2 
652                      dAts.append(at2) 
653                      #dAts.append(b.neighborAtom(at)) 
654          else: 
655              dAts = nodes 
656          #get the sp2 hybridized possible donors which are all ns 
657          sp2 = [] 
658          for t in ['Nam', 'Ng+', 'Npl']: 
659              if t in donorList: 
660                  sp2.append(t) 
661          #ntypes = ['Nam', 'Ng+', 'Npl'] 
662   
663          sp2DAts = None 
664          if len(sp2): 
665              sp2DAts = AtomSet(dAts.get(lambda x, sp2=sp2: x.babel_type in sp2)) 
666   
667          hsp2 = AtomSet([]) 
668          if sp2DAts: 
669              if hats: 
670                  hsp2 = AtomSet(hats.get(lambda x, sp2DAts=sp2DAts:x.bonds[0].atom1 \ 
671                          in sp2DAts or x.bonds[0].atom2 in sp2DAts)) 
672          if sp2DAts: 
673              #remove any sp2 N atoms which already have 3 bonds not to hydrogens 
674              n2Dons = AtomSet(sp2DAts.get(lambda x: x.element=='N')) 
675              if n2Dons: 
676                  n2Dons.bl=0 
677                  for at in n2Dons: 
678                      for b in at.bonds: 
679                          if type(b.bondOrder)==type(2): 
680                              at.bl = at.bl + b.bondOrder 
681                          else: 
682                              at.bl = at.bl + 2 
683                          #allow that there might already be a hydrogen 
684                      nH = at.findHydrogens() 
685                      at.bl = at.bl - len(nH) 
686                  badAts = AtomSet(n2Dons.get(lambda x: x.bl>2)) 
687                  if badAts: 
688                      sp2DAts = sp2DAts - badAts 
689                  delattr(n2Dons,'bl') 
690          #get the sp3 hybridized possible donors 
691          sp3 = [] 
692          for t in ['N3+', 'S3', 'O3']: 
693              if t in donorList: 
694                  sp3.append(t) 
695          n3DAts = None 
696          if 'N3+' in sp3: 
697              n3DAts = AtomSet(dAts.get(lambda x: x.babel_type=='N3+')) 
698          o3DAts = None 
699          if 'O3' in sp3: 
700              o3DAts = AtomSet(dAts.get(lambda x: x.babel_type=='O3')) 
701          if o3DAts: 
702              #remove any O3 atoms which already have 2 bonds not to hydrogens 
703              badO3s = AtomSet([]) 
704              for at in o3DAts: 
705                  if len(at.bonds)<2: continue 
706                  if len(at.findHydrogens()): continue 
707                  else: 
708                      badO3s.append(at) 
709              if len(badO3s): 
710                  o3DAts = o3DAts - badO3s 
711          s3DAts = None 
712          if 'S3' in sp3: 
713              s3DAts = AtomSet(dAts.get(lambda x: x.babel_type=='S3')) 
714          sp3DAts = AtomSet([]) 
715          for item in [n3DAts, o3DAts, s3DAts]: 
716              if item: 
717                  sp3DAts = sp3DAts + item 
718          hsp3 = AtomSet([]) 
719          if sp3DAts: 
720              if hats: 
721                  hsp3 = AtomSet(hats.get(lambda x, sp3DAts=sp3DAts:x.bonds[0].atom1 \ 
722                      in sp3DAts or x.bonds[0].atom2 in sp3DAts)) 
723          hsp = hsp2 + hsp3 
724          #print 'hsp=', hsp.name 
725          #print 'sp2DAts=', sp2DAts.name 
726          #print 'sp3DAts=', sp3DAts.name 
727          return hsp, sp2DAts, sp3DAts 
728   
729 -    def getAcceptors(self, nodes, paramDict): 
730          acceptorList = paramDict['acceptorTypes'] 
731          #print 'acceptorList=', acceptorList 
732   
733          sp2 = [] 
734          for t in ['Npl', 'Nam']: 
735              if t in acceptorList: sp2.append(t) 
736          n2Accs = None 
737          if 'Npl' in sp2: 
738              n2Accs = AtomSet(nodes.get(lambda x: x.babel_type=='Npl')) 
739          if 'Nam' in sp2: 
740              n2Accs2 = AtomSet(nodes.get(lambda x: x.babel_type=='Nam')) 
741              if n2Accs2: 
742                  if n2Accs: 
743                      n2Accs = n2Accs+n2Accs2 
744                  else: 
745                      n2Accs = n2Accs2 
746          if n2Accs is None:  
747              n2Accs = AtomSet([]) 
748   
749          o_sp2 = [] 
750          for t in ['O2', 'O-']: 
751              if t in acceptorList: sp2.append(t) 
752   
753          o2Accs = None 
754          if 'O2' in o_sp2: 
755              o2Accs = AtomSet(nodes.get(lambda x: x.babel_type=='O2')) 
756          if 'O-' in sp2: 
757              o2Accs2 = AtomSet(nodes.get(lambda x: x.babel_type=='O-')) 
758              if o2Accs2: 
759                  if o2Accs: 
760                      o2Accs = o2Accs+o2Accs2 
761                  else: 
762                      o2Accs = o2Accs2 
763          if o2Accs is None:  
764              o2Accs = AtomSet([]) 
765   
766           
767          o3Accs = None 
768          if 'O3' in acceptorList: 
769              o3Accs = AtomSet(nodes.get(lambda x: x.babel_type=='O3')) 
770          if o3Accs is None: o3Accs = AtomSet([]) 
771   
772          s3Accs = None 
773          if 'S3' in acceptorList: 
774              s3Accs = AtomSet(nodes.get(lambda x: x.babel_type=='S3')) 
775          if s3Accs is None: s3Accs = AtomSet([]) 
776   
777          ret2Ats = AtomSet([]) 
778          for item in [n2Accs, o2Accs]: 
779              ret2Ats = ret2Ats + item 
780   
781          ret3Ats = AtomSet([]) 
782          for item in [s3Accs, o3Accs]: 
783              ret3Ats = ret3Ats + item 
784          if ret2Ats: print 'ret2Ats=', ret2Ats.name 
785          else: print 'no ret2Ats' 
786          if ret3Ats: print 'ret3Ats=', ret3Ats.name 
787          else: print 'no ret3Ats' 
788          return ret2Ats, ret3Ats 
789