hdiff output

r29528/cistrans.py 2015-11-26 13:30:05.898523122 +0000 r29527/cistrans.py 2015-11-26 13:30:06.094525708 +0000
  1: #!/usr/bin/python  1: #!/usr/bin/python
  2:   2: 
  3: import read_amber as ra  3: import read_amber as ra
  4: import networkx as nx  4: import networkx as nx
  5: import sys  5: import sys
  6: from collections import defaultdict 
  7:   6: 
  8: class GhostAtom(ra.Atom):  7: class GhostAtom(ra.Atom):
  9:   pass  8:   pass
 10:   9: 
 11: # Set the default valence to 0, if the element isn't present in the defined_valences dictionary. 10: valences = { "H": 1,
 12: valences = defaultdict(int) 
 13: defined_valences = { "H": 1, 
 14:              "C": 4, 11:              "C": 4,
 15:              "N": 3, 12:              "N": 3,
 16:              "O": 2, 13:              "O": 2,
 17:              "S": 2, 14:              "S": 2,
 18:              "Se": 2, 15:              "Se": 2,
 19:              "P": 5 } 16:              "P": 5 }
 20:  17: 
 21: for k, v in defined_valences.items(): 
 22:     valences[k] = v 
 23:  
 24: def peptide_bonds(atoms): 18: def peptide_bonds(atoms):
 25:   peptide_bond_list = [] 19:   peptide_bond_list = []
 26:   amide_carbons = [] 20:   amide_carbons = []
 27: # Identify amide carbons. 21: # Identify amide carbons.
 28:   for carbon in [atom for atom in atoms if atom.element is "C"]: 22:   for carbon in [atom for atom in atoms if atom.element is "C"]:
 29:     if isinstance(carbon, GhostAtom): 23:     if isinstance(carbon, GhostAtom):
 30:       continue 24:       continue
 31:     neighbors = atoms.neighbors(carbon) 25:     neighbors = atoms.neighbors(carbon)
 32:     nitrogens = [neighbor for neighbor in neighbors if neighbor.element is "N"] 26:     nitrogens = [neighbor for neighbor in neighbors if neighbor.element is "N"]
 33:     oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and  27:     oxygens = [neighbor for neighbor in neighbors if neighbor.element is "O" and 
 46:         continue 40:         continue
 47:     n_neighbors = atoms.neighbors(nitrogen) 41:     n_neighbors = atoms.neighbors(nitrogen)
 48:     hydrogens = [neighbor for neighbor in n_neighbors if neighbor.element is "H"] 42:     hydrogens = [neighbor for neighbor in n_neighbors if neighbor.element is "H"]
 49:     # Skip amide groups (where the N has two H neighbours) 43:     # Skip amide groups (where the N has two H neighbours)
 50:     if len(hydrogens) > 1: 44:     if len(hydrogens) > 1:
 51:         continue 45:         continue
 52:     hydrogen = hydrogens[0] 46:     hydrogen = hydrogens[0]
 53:     peptide_bond_list.append((oxygen, carbon, nitrogen, hydrogen)) 47:     peptide_bond_list.append((oxygen, carbon, nitrogen, hydrogen))
 54:   return peptide_bond_list 48:   return peptide_bond_list
 55:  49: 
 56: def not_in_dict(atoms): 
 57:   """ 
 58:   Writes an error message to `cistrans_warning` and stderr, if an atom in the topology file is not present in the 
 59:   `defined_valences` dictionary. 
 60:  
 61:   :param atoms: atom graph 
 62:   """ 
 63:   missing = [] 
 64:   for atom in atoms.nodes(): 
 65:     try: 
 66:       _ = defined_valences[atom.element] 
 67:     except KeyError: 
 68:       missing.append(str(atom)) 
 69:   if missing: 
 70:     with open("cistrans_warning", "w") as error_file: 
 71:       error_file.write("The following atoms were not present in the defined_valences dictionary\n") 
 72:       error_file.write("of the cistrans script, and given a valence of 0 (not bonded):\n\n") 
 73:       sys.stderr.write("The following atoms were not present in the defined_valences dictionary\n") 
 74:       sys.stderr.write("of the cistrans script, and given a valence of 0 (not bonded):\n\n") 
 75:       for atom_str in missing: 
 76:         error_file.write(atom_str + "\n") 
 77:       for atom_str in missing: 
 78:         sys.stderr.write(atom_str + "\n") 
 79:  
 80: def multi_bonds(atoms): 50: def multi_bonds(atoms):
 81:   not_in_dict(atoms) 
 82:   multibonded = [atom for atom in atoms.nodes()  51:   multibonded = [atom for atom in atoms.nodes() 
 83:                  if len(nx.edges(atoms, atom)) < valences[atom.element]] 52:                  if len(nx.edges(atoms, atom)) < valences[atom.element]]
 84:   for i, atom in enumerate(multibonded): 53:   for i, atom in enumerate(multibonded):
 85:     paired = False 54:     paired = False
 86:     for other in atoms.neighbors(atom): 55:     for other in atoms.neighbors(atom):
 87:       if isinstance(other, GhostAtom): 56:       if isinstance(other, GhostAtom):
 88:         paired = True 57:         paired = True
 89:         continue 58:         continue
 90:       if len(nx.edges(atoms, other)) < valences[other.element]: 59:       if len(nx.edges(atoms, other)) < valences[other.element]:
 91:         ghost_atom = GhostAtom(**(atom.__dict__)) 60:         ghost_atom = GhostAtom(**(atom.__dict__))


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0