How to use pdb2pqr - 10 common examples

def protonateBackBone(self, residue, C, O):
        """
        Protonates an atom, X1, given a direction (X2 -> X3) [X1, X2, X3]
        """
        N = residue.getAtom(name='N')

        if C == None and O == None:
            """ do nothing, first residue """
        elif N == None:
            pka_print("could not find N atom in '%s' (protonateBackBone())" % (residue.label))
            exit(9)
        elif residue.resName == "PRO":
            """ do nothing, proline doesn't have a proton """
        else:
            H = self.protonateDirection(atoms=[N, O, C])
            residue.atoms.append(H)

        return residue.getAtom(name='C'), residue.getAtom(name='O')

def __init__(self, atom1=0, atom2=0):
        self.vectors = []
        self.keys = []

        # store vectors for all configurations of atoms
        if atom1 != 0:
            self.keys = lib.get_sorted_configurations(atom1.configurations.keys())
            if atom2 != 0:
                keys2 = lib.get_sorted_configurations(atom2.configurations.keys())
                if self.keys != keys2:
                    pka_print("ERROR: Inequivalent configurations for atoms, please correct your pdbfile to single configuration")
                    pka_print("%s\n%s" % (atom1, atom2))
                    exit(8)
                    #raise 'Cannot make multi vector: Atomic configurations mismatch for\n   %s\n   %s\n'%(atom1,atom2)
            for key in self.keys:
                atom1.setConfiguration(key)
                if atom2 != 0:
                    atom2.setConfiguration(key)
                v = vector(atom1=atom1, atom2=atom2)
                self.vectors.append(v)
                # pka_print(key,v)
        return

for residue in self.residues:
            if not isinstance(residue, aa.Amino):
                continue
            residue.dihedrals = []

            refangles = residue.reference.dihedrals
            for dihed in refangles:
                coords = []
                atoms = dihed.split()
                for i in range(4):
                    atomname = atoms[i]
                    if residue.has_atom(atomname):
                        coords.append(residue.get_atom(atomname).coords)

                if len(coords) == 4:
                    angle = util.dihedral(coords[0], coords[1], coords[2], coords[3])
                else:
                    angle = None

                residue.add_dihedral_angle(angle)

closeresidue = closeatom.residue
            if closeresidue == residue:
                continue
            if not isinstance(closeresidue, (aa.Amino, aa.WAT)):
                continue
            if isinstance(residue, aa.CYS):
                if residue.ss_bonded_partner == closeatom:
                    continue

            # Also ignore if this is a donor/acceptor pair
            if atom.is_hydrogen and atom.bonds[0].hdonor and closeatom.hacceptor:
                continue
            if closeatom.is_hydrogen and closeatom.bonds[0].hdonor and atom.hacceptor:
                continue

            dist = util.distance(atom.coords, closeatom.coords)

            if isinstance(closeresidue, aa.WAT):
                if dist < bestwatdist:
                    bestwatdist = dist
                    bestwatatom = closeatom
            else:
                if dist < bestdist:
                    bestdist = dist
                    bestatom = closeatom

        if bestdist > bestwatdist:
            txt = ("Skipped atom during water optimization: %s in %s skipped "
                   "when optimizing %s in %s") % (bestwatatom.name,
                                                  bestwatatom.residue,
                                                  atom.name, residue)
            _LOGGER.warning(txt)

else:
                hname1 = name

        bondatom1 = residue.get_atom(optinstance.map[hname1].bond)
        bondatom2 = residue.get_atom(optinstance.map[hname2].bond)
        longflag = 0

        # If one bond in the group is significantly (0.05 A)
        # longer than the other, use that group only
        for pivotatom in bondatom1.bonds:
            if not pivotatom.is_hydrogen:
                the_pivatom = pivotatom
                break

        dist1 = distance(the_pivatom.coords, bondatom1.coords)
        dist2 = distance(the_pivatom.coords, bondatom2.coords)

        order = [hname1, hname2]

        if dist2 > dist1 and abs(dist1 - dist2) > 0.05:
            longflag = 1
            order = [hname2, hname1]

        elif dist1 > dist2 and abs(dist1 - dist2) > 0.05:
            longflag = 1
            order = [hname1, hname2]

        for hname in order:
            bondatom = residue.get_atom(optinstance.map[hname].bond)

            # First mirror the hydrogen about the same donor
            for dihedral in residue.reference.dihedrals:

if self.has_atom(bond) and bond.startswith("H"):
                    hatoms.append(self.get_atom(bond))

            # If this is more than two something is wrong
            if len(hatoms) != 2:
                return False

            # Use the existing hydrogen and rotate about the bond
            self.rotate_tetrahedral(nextatom, bondatom, 120)
            newcoords1 = hatoms[0].coords
            self.rotate_tetrahedral(nextatom, bondatom, 120)
            newcoords2 = hatoms[0].coords
            self.rotate_tetrahedral(nextatom, bondatom, 120)

            # Determine which one hatoms[1] is not in
            if util.distance(hatoms[1].coords, newcoords1) > 0.1:
                self.create_atom(atomname, newcoords1)
            else:
                self.create_atom(atomname, newcoords2)

            return True
        return False

continue
            if isinstance(residue, aa.CYS):
                if residue.ss_bonded_partner == closeatom:
                    continue

            # Also ignore if this is a donor/acceptor pair
            pair_ignored = False
            if atom.is_hydrogen and len(atom.bonds) != 0 and atom.bonds[0].hdonor \
               and closeatom.hacceptor:
                continue

            if closeatom.is_hydrogen and len(closeatom.bonds) != 0 and closeatom.bonds[0].hdonor \
                   and atom.hacceptor:
                continue

            dist = util.distance(atom.coords, closeatom.coords)
            other_size = BUMP_HYDROGEN_SIZE if closeatom.is_hydrogen else BUMP_HEAVY_SIZE
            cutoff = atom_size + other_size
            if dist < cutoff:
                bumpscore = bumpscore + 1000.0
                if pair_ignored:
                    _LOGGER.debug('This bump is a donor/acceptor pair.')
        _LOGGER.debug('BUMPSCORE %s', str(bumpscore))
        return bumpscore

# We want to ignore the Hs on the acceptor
        if self.is_carboxylic_hbond(donor, acc):

            # Eliminate the closer hydrogen
            hyds = []
            dist = None
            donorhatom = None
            for hatom in self.hlist:
                if hatom.is_hydrogen:
                    hyds.append(hatom)

            if len(hyds) < 2:
                return 1

            dist = distance(hyds[0].coords, donor.coords)
            dist2 = distance(hyds[1].coords, donor.coords)
            # Eliminate hyds[0]
            if dist < dist2:
                self.hlist.remove(hyds[0])
                self.routines.cells.remove_cell(hyds[0])
                residue.remove_atom(hyds[0].name)
                donorhatom = residue.get_atom(hyds[1].name)
            elif hyds[1] in self.hlist:
                self.hlist.remove(hyds[1])
                self.routines.cells.remove_cell(hyds[1])
                residue.remove_atom(hyds[1].name)
                if residue.has_atom(hyds[0].name):
                    donorhatom = residue.get_atom(hyds[0].name)
                elif len(self.hlist) != 0 and residue.has_atom(self.hlist[0].name):
                    donorhatom = residue.get_atom(self.hlist[0].name)

            # If only one H is left, we're done

addname = "H1"
        else:
            addname = "H2"

        _LOGGER.debug("Finalizing %s by adding %s (%i current O bonds)",
                      residue, addname, len(atom.bonds))

        if len(atom.bonds) == 0:

            newcoords = []

            # Build hydrogen away from closest atom
            closeatom = self.routines.get_closest_atom(atom)
            if closeatom != None:
                vec = util.subtract(atom.coords, closeatom.coords)
                dist = util.distance(atom.coords, closeatom.coords)

                for i in range(3):
                    newcoords.append(vec[i]/dist + atom.coords[i])

            else:
                newcoords = util.add(atom.coords, [1.0, 0.0, 0.0])

            residue.create_atom(addname, newcoords)
            self.routines.cells.add_cell(residue.get_atom(addname))

            self.finalize()

        elif len(atom.bonds) == 1:

            # Initialize variables
            pivot = atom.bonds[0]

continue

            # Check the H(D)-A distance
            dist = distance(donorhatom.coords, acc.coords)
            if dist > DIST_CUTOFF:
                continue

            # Ensure no conflicts if H(A)s if present
            flag = 1
            for acchatom in acc.bonds:
                if not acchatom.is_hydrogen:
                    continue
                flag = 0

                # Check the H(D)-H(A) distance
                hdist = distance(donorhatom.coords, acchatom.coords)
                if hdist < 1.5:
                    continue

                # Check the H(D)-H(A)-A angle
                angle = self.get_hbond_angle(donorhatom, acchatom, acc)
                if angle < 110.0:
                    flag = 1

            if flag == 0:
                continue

            # Check the A-D-H(D) angle
            angle = self.get_hbond_angle(acc, donor, donorhatom)
            if angle <= ANGLE_CUTOFF:
                _LOGGER.debug("Found HBOND! %.4f %.4f", dist, angle)
                return 1