How to use the parmed.amber.titratable_residues function in ParmEd

from argparse import ArgumentParser
    from parmed.amber import AmberParm
    from parmed.amber import titratable_residues as residues
    from parmed.exceptions import ParmedError, AmberError
    from parmed.tools.actions import changeRadii, change
except ImportError:
    if sys.version_info < (2, 7):
        raise ImportError('%s requires Python 2.7 or later' %
                          (os.path.split(sys.argv[0])[1]))
    amberhome = os.getenv('AMBERHOME') or '$AMBERHOME'
    raise ImportError('Could not import Amber Python modules. Please make sure '
                      'you have sourced %s/amber.sh (if you are using sh/ksh/'
                      'bash/zsh) or %s/amber.csh (if you are using csh/tcsh)' %
                      (amberhome, amberhome))
else:
    TitratableResidueList = residues.TitratableResidueList
    LineBuffer = residues._LineBuffer

parser = ArgumentParser(epilog='''This program will read a topology file and
                        generate a cein file for constant Redox Potential simulations with
                        sander or pmemd''', usage='%(prog)s [Options]')
parser.add_argument('-v', '--version', action='version', version='%s: %s' %
                    (parser.prog, __version__))
parser.add_argument('-d', '--debug', dest='debug', action='store_const',
                    help='Enable verbose tracebacks to debug this program',
                    const=True, default=False)
group = parser.add_argument_group('Output files')
group.add_argument('-o', '--output', dest='output', metavar='FILE',
                   help='Output file. Defaults to standard output')
group = parser.add_argument_group('Required Arguments')
group.add_argument('-p', dest='prmtop', metavar='FILE', required=False,
                   help='Topology file to be used in constant Redox Potential simulation',

def print_residues(resnames,mode):
    for resname in resnames:
        if not hasattr(residues, resname):
            print ('%s is not titratable\n' % resname)
            sys.exit(0)
        if not getattr(residues, resname).typ == "redox" and mode == 1:
            print ('%s is not a redox-active titratable residue\n' % resname)
            sys.exit(0)
        if getattr(residues, resname).typ == "redox":
            print (str(getattr(residues, resname)) + '\n')

solvated = False
    first_solvent = 0
    if 'WAT' in parm.parm_data['RESIDUE_LABEL']:
        solvated = True
        for i, res in enumerate(parm.parm_data['RESIDUE_LABEL']):
            if res in solvent_ions:
                first_solvent = parm.parm_data['RESIDUE_POINTER'][i]
                break
    main_reslist = TitratableResidueList(system_name=opt.system,
                        solvated=solvated, first_solvent=first_solvent)
    trescnt = 0
    for resnum in resnums:
        resname = parm.parm_data['RESIDUE_LABEL'][resnum-1]
        if not resname in titratable_residues: continue
        res = getattr(residues, resname)
        # Filter out termini (make sure the residue in the prmtop has as many
        # atoms as the titratable residue defined in residues.py)
        if resnum == parm.ptr('nres'):
            natoms = (parm.ptr('natom') + 1 -
                      parm.parm_data['RESIDUE_POINTER'][resnum-1])
        else:
            natoms = (parm.parm_data['RESIDUE_POINTER'][resnum] -
                      parm.parm_data['RESIDUE_POINTER'][resnum-1])
        if natoms != len(res.atom_list): continue
        # If we have gotten this far, add it to the list.
        main_reslist.add_residue(res, resnum,
                                 parm.parm_data['RESIDUE_POINTER'][resnum-1])
        trescnt += 1

    # Prints a warning if the number of titratable residues is larger than 50
    if trescnt > 50: sys.stderr.write('Warning: Your CPIN file has more than 50 titratable residues! pmemd and sander have a\n'

# Open the output file
    if opt.output is None:
        output = sys.stdout
    else:
        output = open(opt.output, 'w')

    main_reslist.write_cpin(output, opt.igb, opt.intdiel, opt.oldfmt, "ph")

    if opt.output is not None:
        output.close()

    if solvated:
        if opt.outparm is None:
            has_carboxylate = False
            for res in main_reslist:
                if res is residues.AS4 or res is residues.GL4 or res is residues.PRN:
                    has_carboxylate = True
                    break
            if has_carboxylate:
                sys.stderr.write(
                        'Warning: Carboxylate residues in explicit solvent '
                        'simulations require a modified topology file!\n'
                        '         Use the -op flag to print one.\n'
                )
        else:
            changeRadii(parm, 'mbondi2').execute()
            change(parm, 'RADII', ':AS4,GL4,PRN@OD=,OE=,O1=,O2=', 1.3).execute()
            parm.overwrite = True
            parm.write_parm(opt.outparm)
    else:
        if opt.outparm is not None:
            sys.stderr.write(

solvated = False
    first_solvent = 0
    if 'WAT' in parm.parm_data['RESIDUE_LABEL']:
        solvated = True
        for i, res in enumerate(parm.parm_data['RESIDUE_LABEL']):
            if res in solvent_ions:
                first_solvent = parm.parm_data['RESIDUE_POINTER'][i]
                break
    main_reslist = TitratableResidueList(system_name=opt.system,
                        solvated=solvated, first_solvent=first_solvent)
    trescnt = 0
    for resnum in resnums:
        resname = parm.parm_data['RESIDUE_LABEL'][resnum-1]
        if not resname in titratable_residues: continue
        res = getattr(residues, resname)
        # Filter out termini (make sure the residue in the prmtop has as many
        # atoms as the titratable residue defined in residues.py)
        if resnum == parm.ptr('nres'):
            natoms = (parm.ptr('natom') + 1 -
                      parm.parm_data['RESIDUE_POINTER'][resnum-1])
        else:
            natoms = (parm.parm_data['RESIDUE_POINTER'][resnum] -
                      parm.parm_data['RESIDUE_POINTER'][resnum-1])
        if natoms != len(res.atom_list): continue
        # If we have gotten this far, add it to the list.
        main_reslist.add_residue(res, resnum,
                                 parm.parm_data['RESIDUE_POINTER'][resnum-1])
        trescnt += 1

    # Prints a warning if the number of titratable residues is larger than 50
    if trescnt > 50: sys.stderr.write('Warning: Your CPEIN file has more than 50 titratable residues! pmemd and sander have a\n'

raise AmberError('Cannot specify -resnames and -notresnames together')

    if opt.intdiel != 1 and opt.intdiel != 2:
        raise AmberError('-intdiel must be either 1 or 2 currently')

    # Set the list of residue names we will be willing to titrate
    titratable_residues = []
    if notresnames is not None:
        for resname in residues.titratable_residues:
            if resname in notresnames: continue
            titratable_residues.append(resname)
    elif resnames is not None:
        for resname in resnames:
            if not resname in residues.titratable_residues:
                raise AmberError('%s is not a titratable residue!' % resname)
            elif not getattr(residues, resname).typ == "redox":
                raise AmberError('%s is not a redox-active titratable residue!' % resname)
            titratable_residues.append(resname)
    else:
        for resname in residues.titratable_residues:
            if getattr(residues, resname).typ == "redox":
                titratable_residues.append(resname)

    solvent_ions = ['WAT', 'Na+', 'Br-', 'Cl-', 'Cs+', 'F-', 'I-', 'K+', 'Li+',
                    'Mg+', 'Rb+', 'CIO', 'IB', 'MG2']

    # Filter titratable residues based on min and max pKa
    new_reslist = []
    for res in titratable_residues:
        if getattr(residues, res).Eo < mineo: continue
        if getattr(residues, res).Eo > maxeo: continue
        new_reslist.append(res)