How to use the pdb2pqr.pdb2pqr.utilities function in pdb2pqr

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)

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

for i in range(len(chain.residues) - 1):
                res1 = chain.residues[i]
                res2 = chain.residues[i + 1]
                if not isinstance(res1, aa.Amino) or not isinstance(res2, aa.Amino):
                    continue
                atom1 = res1.get_atom("C")
                atom2 = res2.get_atom("N")

                if atom1 is not None:
                    res2.peptide_c = atom1
                if atom2 is not None:
                    res1.peptide_n = atom2
                if atom1 is None or atom2 is None:
                    continue

                if util.distance(atom1.coords, atom2.coords) > PEPTIDE_DIST:
                    text = "Gap in backbone detected between %s and %s!" % \
                           (res1, res2)
                    _LOGGER.warning(text)
                    res2.peptide_c = None
                    res1.peptide_n = None

options:      The name of the forcefield (string)

        Returns
            header:  The PQR file header (string)
            lines:   The PQR file atoms (list)
            missedligandresidues:  A list of ligand residue names whose charges could
                     not be assigned (ligand)
            protein: The protein object
    """

    pkaname = ""
    lines = []
    Lig = None
    atomcount = 0   # Count the number of ATOM records in pdb

    outroot = utilities.getPQRBaseFileName(options.output_pqr)

    if options.pka_method == 'propka':
        pkaname = outroot + ".propka"
        #TODO: What? Shouldn't it be up to propka on how to handle this?
        if os.path.isfile(pkaname):
            os.remove(pkaname)

    start = time.time()
    _LOGGER.info("Beginning PDB2PQR...")

    myDefinition = definitions.Definition()
    _LOGGER.info("Parsed Amino Acid definition file.")

    if options.drop_water:
        # Remove the waters
        pdblist_new = []

def rotate_tetrahedral(cls, atom1, atom2, angle):
        """Rotate about the atom1-atom2 bond by a given angle All atoms connected
        to atom2 will rotate.

        Args:
            atom1:  The first atom of the bond to rotate about (atom)
            atom2:  The second atom of the bond to rotate about (atom)
            angle:  The number of degrees to rotate (float)
        """
        moveatoms = []
        movecoords = []
        initcoords = util.subtract(atom2.coords, atom1.coords)

        # Determine which atoms to rotate
        for atom in atom2.bonds:
            if atom == atom1:
                continue
            moveatoms.append(atom)
            movecoords.append(util.subtract(atom.coords, atom1.coords))

        newcoords = quat.qchichange(initcoords, movecoords, angle)
        for iatom, atom in enumerate(moveatoms):
            x = newcoords[iatom][0] + atom1.x
            y = newcoords[iatom][1] + atom1.y
            z = newcoords[iatom][2] + atom1.z
            atom.x = x
            atom.y = y
            atom.z = z

y = (newcoords[iatom][1] + coordlist[1][1])
            z = (newcoords[iatom][2] + coordlist[1][2])
            atom.x = x
            atom.y = y
            atom.z = z
            self.cells.add_cell(atom)

        # Set the new angle
        coordlist = []
        for atomname in atomnames:
            if residue.has_atom(atomname):
                coordlist.append(residue.get_atom(atomname).coords)
            else:
                raise ValueError("Error occurred while trying to debump!")

        dihed = util.dihedral(coordlist[0], coordlist[1], coordlist[2], coordlist[3])
        residue.dihedrals[anglenum] = dihed

closeatoms = self.debumper.cells.get_near_cells(atom)
                for closeatom in closeatoms:

                    # Conditions for continuing
                    if atom.residue == closeatom.residue:
                        continue
                    if not (closeatom.hacceptor or closeatom.hdonor):
                        continue
                    if atom.hdonor and not atom.hacceptor:
                        if not closeatom.hacceptor:
                            continue
                    if atom.hacceptor:
                        if not atom.hdonor and not closeatom.hdonor:
                            continue

                    dist = util.distance(atom.coords, closeatom.coords)
                    if dist < 4.3:
                        residue = atom.residue
                        hbond = structures.PotentialBond(atom, closeatom, dist)

                        # Store the potential bond
                        obj.hbonds.append(hbond)

                        # Keep track of connectivity
                        if closeatom in self.atomlist:
                            closeobj = self.resmap[closeatom.residue]
                            if closeobj not in connectivity[obj]:
                                connectivity[obj].append(closeobj)

            progress += increment
            while progress >= 0.0499:
                progress -= 0.05

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

            # Also ignore if this is a donor/acceptor pair
            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:
                nearatoms[closeatom] = cutoff - dist

        return nearatoms

return 0
        elif residue.has_atom("LP1"):
            newname = "LP2"
        else:
            newname = "LP1"

        _LOGGER.debug("Working on %s %s (acceptor) to %s %s (donor)",
                      acc.residue, acc.name, donor.residue, donor.name)

        # Act depending on the number of bonds
        if len(acc.bonds) == 0:

            if self.is_hbond(donor, acc):

                # Find the best donor hydrogen and use that
                bestdist = util.distance(acc.coords, donor.coords)
                for donorh in donor.bonds:
                    dist = util.distance(acc.coords, donorh.coords)
                    if dist < bestdist:
                        bestdist = dist

                # Point the LP to the best H
                # TODO - this looks like a bug...
                # donorh ends up being the last item in donor.bonds
                # This may be fixed by setting a best_donorh to go with bestdist
                # and using best_donorh in the function below
                self.make_atom_with_no_bonds(acc, donorh, newname)
                _LOGGER.warning("The best donorH was not picked (BUG?).")
                _LOGGER.debug("Added %s to %s", newname, acc.residue)
                return 1
            return 0

if not residue.has_atom("H1"):
            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