Source Code for Module PyBabel.atomTypes

  1  ############################################################################# 
  2  # 
  3  # Author: Michel F. SANNER 
  4  # Reimplemented from Babel v1.6 from Pat Walters and Math Stahl 
  5  # 
  6  # Copyright: M. Sanner TSRI 2000 
  7  # 
  8  ############################################################################# 
  9  # 
 10  # $Header: /opt/cvs/python/packages/share1.5/PyBabel/atomTypes.py,v 1.3 2006/12/01 19:17:38 sargis Exp $ 
 11  # 
 12  # $Id: atomTypes.py,v 1.3 2006/12/01 19:17:38 sargis Exp $ 
 13  # 
 14   
 15  """ 
 16  This file implements the AtomHybridization class that can be used to assign 
 17  atom types. 
 18   
 19  example: 
 20       
 21        >>> atype = AtomHybridization() 
 22        >>> atype.assignHybridization(atoms) 
 23   
 24        atoms has to be a list of atom objects 
 25        Atom: 
 26            a.element : atom's chemical element symbol (string) 
 27            a.coords : 3-sequence of floats 
 28            a.bonds : list of Bond objects 
 29        Bond: 
 30            b.atom1 : instance of Atom 
 31            b.atom2 : instance of Atom 
 32   
 33        after completion each atom has the following additional members: 
 34            babel_type: string 
 35            babel_atomic_number: int 
 36            babel_organic 
 37         
 38  reimplmentation of Babel1.6 in Python by Michel Sanner April 2000 
 39  Original code by W. Patrick Walters and Matthew T. Stahl  
 40  """ 
 41   
 42   
 43   
 44  import string 
 45   
 46  from babelAtomTypes import babel_types 
 47  from babelElements import babel_elements 
 48  from util import * 
 49   
 50   
 51  # for valence_three 
 52  SP3_MAX     = 114.8 
 53  MAY_BE_SP2  = 122.0 
 54  SP_MIN      = 160.0 
 55   
 56  # for valence_one 
 57  V1_C1_C1_CUTOFF = 1.22 
 58  V1_C2_C_CUTOFF  = 1.41 
 59  V1_C2_N_CUTOFF  = 1.37 
 60   
 61  V1_N1_C1_CUTOFF = 1.20  
 62  V1_N3_C_CUTOFF  = 1.38 
 63  V1_N3_N3_CUTOFF = 1.43 
 64  V1_N3_N2_CUTOFF = 1.41 
 65   
 66  V1_O2_C2_CUTOFF = 1.30 
 67  V1_O2_AS_CUTOFF = 1.685 
 68   
 69  V1_S2_C2_CUTOFF = 1.76 
 70  V1_S2_AS_CUTOFF = 2.11 
 71   
 72  V2_C3_C_CUTOFF  = 1.53 
 73  V2_C3_N_CUTOFF  = 1.46 
 74  V2_C3_O_CUTOFF  = 1.44 
 75   
 76  V2_N2_C_CUTOFF  = 1.38 
 77  V2_N2_N_CUTOFF  = 1.32 
 78   
 79  V2_C2_C_CUTOFF  = 1.42 
 80  V2_C2_N_CUTOFF  = 1.41 
 81  GEN_C3_C_CUTOFF = 1.45 
 82   
 83   
 84 -class AtomHybridization: 
 85   
 86   
 87 -    def __init__(self): 
 88          """constructor""" 
 89          self.atoms = None 
 90           
 91           
 92 -    def get_atomic_number(self, name): 
 93          """return the element number for a given name or raises a 
 94          ValueError exception if the element is not known""" 
 95          _name = string.upper(name[0]) 
 96          if len(name)>1: 
 97              if not name[1] in string.digits: 
 98                  _name = _name + string.lower(name[1]) 
 99          if _name in babel_elements.keys(): 
100              return babel_elements[_name]['num']  
101          else: 
102              raise ValueError( "Could not find atomic number for %s %s"% \ 
103                                (name,_name) ) 
104   
105   
106 -    def assignHybridization(self, atoms): 
107          """atoms is a list of objects of type Atom having the following 
108          members: 
109          Atom: 
110              a.element : atom's chemical element symbol (string) 
111              a.coords : 3-sequence of floats 
112              a.bonds : list of Bond objects 
113          Bond: 
114              b.atom1 : instance of Atom 
115              b.atom2 : instance of Atom 
116   
117          after completion each atom has the following additional members: 
118          babel_type: string 
119          babel_atomic_number: int 
120          babel_organic 
121          """ 
122           
123          self.atoms = atoms 
124          for a in self.atoms: 
125              a.babel_type = a.element 
126              a._babel_redo = 0 
127              a.babel_atomic_number = self.get_atomic_number(a.babel_type) 
128   
129              if a.babel_type[0] in [ 'C', 'H', 'O', 'N', 'S', 'P' ]: 
130                  a.babel_organic=1 
131              else: a.babel_organic = 0 
132           
133          self.phase1() 
134          self.valence_four() 
135          self.valence_three() 
136          self.valence_two() 
137          self.valence_one() 
138   
139          self.phase4() 
140          self.phase5()  
141          self.phase6() 
142          self.check_for_amides() 
143   
144          # cleanup 
145          for a in self.atoms: 
146              delattr(a,'_babel_redo') 
147   
148          delattr(self,'atoms') 
149   
150   
151 -    def count_heavy_atoms(self, atom): 
152          count = 0 
153          for b in atom.bonds: 
154              bonded_atom = b.atom1 
155              if bonded_atom==atom: bonded_atom=b.atom2 
156              if bonded_atom.babel_type[0] == 'H': count = count + 1 
157          return len(atom.bonds) - count 
158   
159   
160 -    def count_free_ox(self, atom): 
161          free_O_count=0 
162          for b in atom.bonds: 
163              bonded_atom = b.atom1 
164              if bonded_atom==atom: bonded_atom=b.atom2 
165              if bonded_atom.babel_type[0] == 'O' and \ 
166                 self.count_heavy_atoms(bonded_atom) == 1: 
167                  free_O_count = free_O_count+1 
168          return free_O_count 
169   
170           
171 -    def phase1(self): 
172          for a in self.atoms: 
173              if a.babel_type[0] == 'H': 
174                  a.babel_type = 'H' 
175   
176                  if len(a.bonds): 
177                      k = a.bonds[0].atom1 
178                      if k==a: k = a.bonds[0].atom2 
179                      if k.babel_type[0] == 'C': 
180                          a.babel_type = 'HC' 
181   
182   
183 -    def valence_four(self): 
184   
185        for a in self.atoms: 
186          if len(a.bonds) == 4 and a.babel_organic: 
187   
188              if a.babel_type[0] == 'C': 
189                  if a.babel_type=='C': a.babel_type = "C3" 
190   
191              elif a.babel_type[0] == 'N': 
192                  if self.count_free_ox(a) >= 1: a.babel_type = "Nox" 
193                  else: a.babel_type = "N3+" 
194   
195              elif a.babel_type[0] == 'P': 
196                  if len(a.babel_type) == 1: 
197                      count = self.count_free_ox(a) 
198                      if count >= 2: a.babel_type = "Pac" 
199                      elif count == 1: a.babel_type = "Pox" 
200                      else: a.babel_type = "P3+" 
201   
202              elif a.babel_type[0] == 'S': 
203                  if a.babel_type=='S': 
204                      count = self.count_free_ox(a) 
205                      if count >= 3: a.babel_type = "Sac" 
206                      elif count >= 1: a.babel_type = "Sox" 
207                      else: a.babel_type = "S" 
208   
209              elif a.babel_type[0] == 'B': 
210                  count = self.count_free_ox(a) 
211                  if count >= 3: a.babel_type = "Bac" 
212                  if count >= 1: a.babel_type = "Box" 
213                  else: a.babel_type = "B" 
214   
215   
216 -    def valence_three(self): 
217   
218        for a in self.atoms: 
219            if len(a.bonds) == 3 and a.babel_organic: 
220   
221                k = a.bonds[0].atom1 
222                if k==a: k = a.bonds[0].atom2 
223                l = a.bonds[1].atom1 
224                if l==a: l = a.bonds[1].atom2 
225                m = a.bonds[2].atom1 
226                if m==a: m = a.bonds[2].atom2 
227   
228                angle1 = bond_angle(k.coords, a.coords, l.coords) 
229                angle2 = bond_angle(k.coords, a.coords, m.coords) 
230                angle3 = bond_angle(l.coords, a.coords, m.coords) 
231                avg_angle = (angle1 + angle2 + angle3)/3 
232   
233                if a.babel_type[0] =='C': 
234                    if avg_angle < SP3_MAX: a.babel_type="C3" 
235                    elif self.count_free_ox(a) >= 2: a.babel_type="Cac" 
236                    else: a.babel_type = "C2" 
237   
238                elif a.babel_type[0] =='N': 
239                    if avg_angle < SP3_MAX: a.babel_type = "N3" 
240                    elif self.count_free_ox(a) >= 2: a.babel_type = "Ntr" 
241                    else: a.babel_type = "Npl" 
242   
243                elif a.babel_type[0] =='B': 
244                    if self.count_free_ox(a) >= 1: a.babel_type = "Box" 
245                    else: a.babel_type = "B" 
246   
247                elif a.babel_type =='S': 
248                    if self.count_free_ox(a) >= 1: a.babel_type = "Sox" 
249                    else: a.babel_type = "S3+"       
250   
251   
252 -    def valence_two(self): 
253   
254        for a in self.atoms: 
255            if len(a.bonds) == 2 and a.babel_organic: 
256   
257                k = a.bonds[0].atom1 
258                if k==a: k = a.bonds[0].atom2 
259                l = a.bonds[1].atom1 
260                if l==a: l = a.bonds[1].atom2 
261   
262                angle1 = bond_angle(k.coords, a.coords, l.coords) 
263   
264                if a.babel_type[0] == 'C': 
265                    if a.babel_type =="C": 
266                        if angle1 < SP3_MAX: 
267                            a.babel_type = "C3" 
268                            a._babel_redo = 1 
269                        elif angle1 < SP_MIN: 
270                            a.babel_type = "C2" 
271                            if angle1 < MAY_BE_SP2: 
272                                a._babel_redo = 3 
273                        else: a.babel_type = "C1" 
274   
275                elif a.babel_type[0] == 'N': 
276                    if angle1 <= SP3_MAX: 
277                        a.babel_type = "N3" 
278                        a._babel_redo = 2 
279                    elif angle1 <= SP_MIN: a.babel_type = "Npl" 
280                    else: a.babel_type = "N1" 
281   
282                elif a.babel_type[0] == 'O': 
283                  a.babel_type = "O3" 
284   
285                elif a.babel_type[0] == 'S': 
286                    if a.babel_type =="S": a.babel_type = "S3" 
287   
288   
289 -    def valence_one(self): 
290   
291        for a in self.atoms: 
292   
293            if len(a.bonds) == 1 and a.babel_organic: 
294                k = a.bonds[0].atom1 
295                if k==a: k = a.bonds[0].atom2 
296                bond_length = distance(a.coords, k.coords) 
297   
298                if a.babel_type[0] == 'C': 
299                    if a.babel_type=="C": 
300                        if k.babel_type[:2]=='C1' and \ 
301                           bond_length <= V1_C1_C1_CUTOFF: 
302                            a.babel_type = "C1" 
303                        elif k.babel_type[0] == "C" and \ 
304                             bond_length <= V1_C2_C_CUTOFF: 
305                            a.babel_type = "C2" 
306                        else: a.babel_type = "C3" 
307   
308                    if k.babel_type[0]=="N": 
309                        if bond_length <= V1_C2_N_CUTOFF: a.babel_type = "C2" 
310                        else: a.babel_type = "C3" 
311   
312                if a.babel_type[0] == 'N': 
313                     if a.babel_type=="N": 
314                        if k.babel_type[:2]=='C1' and \ 
315                           bond_length <= V1_N1_C1_CUTOFF: 
316                            a.babel_type = "N1" 
317                        elif (k.babel_type[:2] == "C2" or \ 
318                              k.babel_type[:2] == "C3") \ 
319                             and bond_length > V1_N3_C_CUTOFF: 
320                            a.babel_type = "N3" 
321                        elif a.babel_type[:2]== "N3" and \ 
322                             bond_length > V1_N3_N3_CUTOFF: 
323                            a.babel_type = "N3" 
324                        elif a.babel_type[:2]== "Npl" and \ 
325                             bond_length > V1_N3_N2_CUTOFF: 
326                            a.babel_type = "N3" 
327                        else: 
328                            a.babel_type = "Npl" 
329   
330                if a.babel_type[0] == 'O': 
331                    if a.babel_type=="O": 
332                        if k.babel_type in ["Cac", "Pac", "Sac", "Ntr"]: 
333                            a.babel_type = "O-" 
334                        elif k.babel_type in ["Nox", "Pox", "Sox"]: 
335                            a.babel_type = "O2" 
336                        elif k.babel_type[0] =="C" and \ 
337                             bond_length <= V1_O2_C2_CUTOFF: 
338                            a.babel_type = "O2" 
339                            k.babel_type = "C2" 
340                            k._babel_redo = 0 
341                        elif k.babel_type=="As" and \ 
342                             bond_length <= V1_O2_AS_CUTOFF: 
343                            a.babel_type = "O2" 
344                        else: a.babel_type = "O3" 
345   
346                if a.babel_type[0] == 'S': 
347                    if a.babel_type=="S": 
348                        if k.babel_type[0] =="P": a.babel_type = "S2" 
349                        elif k.babel_type[0]=="C" and \ 
350                             bond_length <= V1_S2_C2_CUTOFF: 
351                            a.babel_type = "S2" 
352                            k.babel_type = "C2" 
353                            a._babel_redo = 0 
354                        elif k.babel_type=="As" and \ 
355                             bond_length <= V1_S2_AS_CUTOFF: 
356                            a.babel_type = "S2" 
357                        else: a.babel_type = "S3" 
358   
359   
360 -    def phase4(self): 
361   
362        for a in self.atoms: 
363            if a._babel_redo==1: 
364                for b in a.bonds: 
365                    j = b.atom1 
366                    if j==a: j = b.atom2 
367                    bond_length = distance(a.coords, j.coords) 
368                    if bond_length <= V2_C2_C_CUTOFF and j.babel_type[0] == 'C': 
369                        a.babel_type = "C2" 
370                    elif bond_length <= V2_C2_N_CUTOFF and j.babel_type[0] =='N': 
371                        a.babel_type = "C2" 
372   
373                for b in a.bonds: 
374                    j = b.atom1 
375                    if j==a: j = b.atom2 
376                    bond_length = distance(a.coords, j.coords) 
377                    if bond_length > V2_C3_C_CUTOFF and j.babel_type[0] == 'C': 
378                        a.babel_type = "C3" 
379                    elif bond_length > V2_C3_N_CUTOFF and j.babel_type[0] == 'N': 
380                        a.babel_type = "C3" 
381                    elif bond_length > V2_C3_O_CUTOFF and j.babel_type[0] == 'O': 
382                        a.babel_type = "C3" 
383   
384            elif a._babel_redo==2: 
385                for b in a.bonds: 
386                    j = b.atom1 
387                    if j==a: j = b.atom2 
388                    bond_length = distance(a.coords, j.coords) 
389                    if bond_length <= V2_N2_C_CUTOFF and j.babel_type[0] == 'C': 
390                        a.babel_type = "Npl" 
391                    elif bond_length <= V2_N2_N_CUTOFF and j.babel_type[0] =='N': 
392                        a.babel_type = "Npl" 
393   
394            elif a._babel_redo==3: 
395                flag = 0 
396                for b in a.bonds: 
397                    j = b.atom1 
398                    if j==a: j = b.atom2 
399                    bond_length = distance(a.coords, j.coords) 
400   
401                    if bond_length <= V2_C2_C_CUTOFF and j.babel_type[0] == 'C': 
402                        a.babel_type = "C2" 
403                        flag = 1 
404                    elif bond_length <= V2_C2_N_CUTOFF and j.babel_type[0] =='N': 
405                        a.babel_type = "C2" 
406                        flag = 1 
407   
408                if flag == 0: 
409                    for b in a.bonds: 
410                        j = b.atom1 
411                        if j==a: j = b.atom2 
412                        bond_length = distance(a.coords, j.coords) 
413                        if bond_length > V2_C3_C_CUTOFF and j.babel_type[0]=='C': 
414                            a.babel_type = "C3" 
415                            flag = 1 
416                        elif bond_length>V2_C3_N_CUTOFF and j.babel_type[0]=='N': 
417                            a.babel_type = "C3" 
418                            flag = 1 
419                        elif bond_length>V2_C3_O_CUTOFF and j.babel_type[0]=='O': 
420                            a.babel_type = "C3" 
421                            flag = 1 
422                        elif flag == 0: 
423                            if bond_length > GEN_C3_C_CUTOFF and \ 
424                               j.babel_type[0] == 'C': 
425                                a.babel_type = "C3" 
426                                flag = 1 
427   
428 -    def phase5(self): 
429   
430          for a in self.atoms: 
431              if a.babel_type == "C2": 
432                  flag = 0; 
433                  for b in a.bonds: 
434                      j = b.atom1 
435                      if j==a: j = b.atom2 
436                      if not j.babel_type in ["C3", "DC", "HC", "N3", "N3+", "O3"] and\ 
437                         not j.babel_type in ["Pac", "Sac", "Sox", "C1", "S3", "Cac" ]: 
438                          flag = 1 
439                  if flag == 0: 
440                      a.babel_type = "C3" 
441   
442   
443 -    def phase6(self): 
444   
445          for a in self.atoms: 
446              no_plus = 1 
447              protonated = 1 
448              if a.babel_type=="N3": 
449                  for b in a.bonds: 
450                      conn = b.atom1 
451                      if conn==a: conn = b.atom2 
452                      # If an unsaturated atom is attached to this nitrogen then 
453                      # it should be Npl 
454                      if len(a.bonds) == 2 and \ 
455                         conn.babel_type in ["Car","C2","Sox","Sac","Pac","So2"]: 
456                          protonated = 0 
457                          a.babel_type = "Npl" 
458   
459                      # If the attached atom is something other that C3, 
460                      # H or D the nitrogen is not positively charged 
461                      if conn.babel_type != "C3" and conn.babel_atomic_number != 1: 
462                          protonated = 0 
463                  if protonated: a.babel_type = "N3+" 
464   
465              # look for guanadinium nitrogens 
466              elif a.babel_type == "C2": 
467                  # First see if we have an sp2 carbon surrounded by 3 sp2 
468                  # nitrogens 
469   
470                  m = 0; 
471                  for b in a.bonds: 
472                      k = b.atom1 
473                      if k==a: k = b.atom2 
474                      if k.babel_type=="Npl" or k.babel_type=="N2" or k.babel_type=="Ng+": m=m+1 
475   
476                  if m == 3: 
477                      a.babel_type = "C+" 
478                      for b in a.bonds: 
479                          k = b.atom1 
480                          if k==a: k = b.atom2 
481                          k.babel_type = "Ng+" 
482   
483              elif a.babel_type == "Cac": 
484                  for b in a.bonds: 
485                      k = b.atom1 
486                      if k==a: k = b.atom2 
487                      if k.babel_type[0]=="O" and self.count_heavy_atoms(k) == 1: 
488                          k.babel_type = "O-" 
489   
490   
491 -    def check_for_carbonyl(self, atom): 
492          for b in atom.bonds: 
493              bonded_atom = b.atom1 
494              if bonded_atom==atom: bonded_atom = b.atom2 
495              if bonded_atom.babel_type=="O2" or bonded_atom.babel_type=="S2": 
496                  return 3 
497          return 2 
498   
499   
500 -    def check_for_amides(self): 
501     
502          for a in self.atoms: 
503              if a.babel_type=="Npl": 
504                  for b in a.bonds: 
505                      conn = b.atom1 
506                      if conn==a: conn = b.atom2 
507                      if conn.babel_type=="Cac" or conn.babel_type=="Sox" or \ 
508                         conn.babel_type=="So2": 
509                          a.babel_type = "Nam" 
510                          break 
511   
512                      if conn.babel_type=="C2": 
513                          if self.check_for_carbonyl(conn) == 3: 
514                              a.babel_type = "Nam" 
515                              break 
516   
517   
518 -    def getProperty(self, property, elements): 
519          """list <- getProperty(property, elements) 
520   
521          property has to be 'num', 'cov_rad', 'bond_ord_rad', 'vdw_rad', 
522                             'bs_rad', 'max_bonds' or 'rgb' 
523          elements is a list of 1 or 2 character(s) strings 
524          """ 
525          if property not in babel_elements["C"].keys(): 
526              raise RuntimeError("Invalid property %s, has to be in %s\n" % \ 
527                                 (property, babel_elements["C"].keys())) 
528          prop = [] 
529          for el in elements: 
530              prop.append(babel_elements[el][property]) 
531          return prop 
532   
533       
534                                  
535 -class TypeConverter: 
536   
537 -    def __init__(self, outputType): 
538          if outputType not in babel_types.keys(): 
539              raise RuntimeError("Invalid format %s\n"%outputType) 
540          self.outputType = outputType 
541   
542   
543 -    def convert(self, input, mode='all_caps'): 
544   
545          try: 
546              i = babel_types['INT'].index(input) 
547              return babel_types[self.outputType][i] 
548          except: 
549              print "Unable to assign %s type to atom %s"%(self.outputType,input) 
550              if mode=='zero': 
551                  return 0 
552              elif mode=='dummy': 
553                  i = babel_types['INT'].index("X") 
554                  return babel_types[self.outputType][i] 
555              elif mode=='all_caps': 
556                  return string.upper(input) 
557              else: return input 
558   
559 -    def clean_atom_type(self, type_name): 
560         name = string.upper(type_name[0]) 
561         if len(type_name) > 1: 
562             name = name + string.lower(type_name[1]) 
563             if name[1] not in string.letters: 
564                 return name[0] 
565         return name 
566   
567   
568  if __name__ == '__main__': 
569      import pdb, sys 
570      from cycle import RingFinder 
571      from bo import BondOrder 
572      from aromatic import Aromatic 
573       
574      from MolKit.pdbParser import NewPdbParser 
575      parser = NewPdbParser("/tsri/pdb/struct/%s.pdb"%sys.argv[1]) 
576      mols = parser.parse() 
577      mol = mols[0] 
578      mol.buildBondsByDistance() 
579      allAtoms = mol.chains.residues.atoms 
580      bonds = allAtoms.bonds[0] 
581   
582      print "assigning atom types" 
583      babel = AtomHybridization() 
584      babel.assignHybridization(allAtoms) 
585