Source code for pdb2pqr.definitions

"""XML handling for biomolecular residue topology definitions.

.. codeauthor::  Jens Erik Nielsen
.. codeauthor::  Todd Dolinsky
.. codeauthor::  Yong Huang
import copy
import re
import logging
from xml import sax
from . import structures
from . import residue

_LOGGER = logging.getLogger(__name__)

[docs]class DefinitionHandler(sax.ContentHandler): """Handle definition XML file content."""
[docs] def __init__(self): self.curelement = "" self.curatom = None self.curholder = None self.curobj = None = {} self.patches = [] self.content = ""
[docs] def startElement(self, name, _): """Start XML element parsing. :param name: element name :type name: str """ if name == "residue": obj = DefinitionResidue() self.curholder = obj self.curobj = obj elif name == "patch": obj = Patch() self.curholder = obj self.curobj = obj elif name == "atom": obj = DefinitionAtom() self.curatom = obj self.curobj = obj else: self.curelement = name
[docs] def endElement(self, name): """End XML element parsing. :param name: element name :type name: str :raises KeyError: for invalid atom or residue names :raises RuntimeError: for unexpected XML features or format """ if name == "residue": # Complete Residue object residue_ = self.curholder if not isinstance(residue_, DefinitionResidue): raise RuntimeError("Internal error parsing XML!") resname = if resname == "": raise KeyError("Residue name not set in XML!") else:[resname] = residue_ self.curholder = None self.curobj = None elif name == "patch": # Complete patch object patch = self.curholder if not isinstance(patch, Patch): raise RuntimeError("Internal error parsing XML!") patchname = if patchname == "": raise KeyError("Residue name not set in XML!") else: self.patches.append(patch) self.curholder = None self.curobj = None elif name == "atom": # Complete atom object atom = self.curatom if not isinstance(atom, DefinitionAtom): raise RuntimeError("Internal error parsing XML!") atomname = if atomname == "": raise KeyError("Atom name not set in XML!") else:[atomname] = atom self.curatom = None self.curobj = self.curholder elif (self.curobj is not None) and (self.content.strip() != ""): self.content = self.content.strip() _LOGGER.debug( f"Got text for {} <{name}>: {self.content}" ) if name == "bond": self.curobj.bonds.append(self.content) elif name == "dihedral": self.curobj.dihedrals.append(self.content) elif name == "altname": self.curholder.altnames[self.content] = elif name == "remove": self.curobj.remove.append(self.content) elif name == "name": = self.content elif self.curelement != "": try: self.content = float(self.content) except ValueError: pass setattr(self.curobj, self.curelement, self.content) self.content = "" self.curelement = "" return
[docs] def characters(self, text): """Parse text data in XML. :param text: text data to parse :type text: str """ self.content += text
[docs]class Definition: """Force field topology definitions. The Definition class contains the structured definitions found in the files and several mappings for easy access to the information. """
[docs] def __init__(self, aa_file, na_file, patch_file): """Initialize object. :param aa_file: file-like object with amino acid definitions :type aa_file: file :param na_file: file-like object with nucleic acid definitions :type na_file: file :param patch_file: file-like object with patch definitions :type patch_file: file """ = {} self.patches = {} handler = DefinitionHandler() sax.make_parser() for def_file in [aa_file, na_file]: sax.parseString(, handler) = {} sax.parseString(, handler) # Apply specific patches to the reference object, allowing users # to specify protonation states in the PDB file for patch in handler.patches: if patch.newname != "": # Find all residues matching applyto resnames = list( for name in resnames: regexp = re.compile(patch.applyto).match(name) if not regexp: continue newname = patch.newname.replace("*", name) self.add_patch(patch, name, newname) # Either way, make sure the main patch name is available self.add_patch(patch, patch.applyto,
[docs] def add_patch(self, patch, refname, newname): """Add a patch to a topology definition residue. :param patch: the patch object to add :type patch: Patch :param refname: the name of the object to add the patch to :type refname: str :param newname: the name of the new (patched) object :type newname: str """ try: aadef =[refname] # The reference patch_residue = copy.deepcopy(aadef) # Add atoms from patch for atomname in[atomname] =[atomname] for bond in[atomname].bonds: if bond not in continue if atomname not in[bond].bonds:[bond].bonds.append(atomname) # Rename atoms as directed for key in patch.altnames: patch_residue.altnames[key] = patch.altnames[key] # Remove atoms as directed for remove in patch.remove: if not patch_residue.has_atom(remove): continue removebonds =[remove].bonds del[remove] for bond in removebonds: if remove in[bond].bonds:[bond].bonds.remove(remove) # Add the new dihedrals for dihedral in patch.dihedrals: patch_residue.dihedrals.append(dihedral) # Point at the new reference[newname] = patch_residue # Store the patch self.patches[newname] = patch except KeyError: # Just store the patch self.patches[newname] = patch
[docs]class Patch: """Residue patches for structure topologies."""
[docs] def __init__(self): = "" self.applyto = "" = {} self.remove = [] self.altnames = {} self.dihedrals = [] self.newname = ""
def __str__(self): text = f"{}\n" text += f"Apply to: {self.applyto}\n" text += "Atoms to add: \n" for atom in text += f"\t{str([atom])}\n" text += "Atoms to remove: \n" for remove in self.remove: text += f"\t{remove}\n" text += "Alternate naming map: \n" text += f"\t{self.altnames}\n" return text
[docs]class DefinitionResidue(residue.Residue): """Force field toplogy representation for a residue."""
[docs] def __init__(self): = "" self.dihedrals = [] = {} self.altnames = {}
def __str__(self): text = f"{}\n" text += "Atoms: \n" for atom in text += f"\t{str([atom])}\n" text += "Dihedrals: \n" for dihedral in self.dihedrals: text += f"\t{dihedral}\n" text += "Alternate naming map: \n" text += f"\t{self.altnames}\n" return text
[docs] def get_nearest_bonds(self, atomname): """Get bonded atoms near a given atom. :param atomname: name of specific atom :type atomname: str :return: list of nearby bonded atom names :rtype: [str] """ bonds = [] lev2bonds = [] atom =[atomname] # Get directly bonded (length = 1) atoms for bondedatom in atom.bonds: if bondedatom not in bonds: bonds.append(bondedatom) # Get bonded atoms 2 bond lengths away for bondedatom in atom.bonds: for bond2 in[bondedatom].bonds: if bond2 not in bonds and bond2 != atomname: bonds.append(bond2) lev2bonds.append(bond2) # Get bonded atoms 3 bond lengths away for lev2atom in lev2bonds: for bond3 in[lev2atom].bonds: if bond3 not in bonds: bonds.append(bond3) return bonds
[docs]class DefinitionAtom(structures.Atom): """Store force field atom topology definitions."""
[docs] def __init__(self, name=None, x=None, y=None, z=None): """Initialize class. :param name: atom name :type name: str :param x: x-coordinate :type x: float :param y: y-coordinate :type y: float :param z: z-coordinate :type z: float """ = name self.x = x self.y = y self.z = z if name is None: = "" if x is None: self.x = 0.0 if y is None: self.y = 0.0 if z is None: self.z = 0.0 self.bonds = []
def __str__(self): text = f"{}: {self.x:.3f} {self.y:.3f} {self.z:.3f}" for bond in self.bonds: text += f" {bond}" return text @property def is_backbone(self): """Identify whether atom is in backbone. :return: true if atom name is in backbone, otherwise false :rtype: bool """ return in structures.BACKBONE