#!/usr/bin/env python
#
# pKa calculations with APBS
#
# Copyright University College Dublin & Washington University St. Louis 2004-2007
# All rights reserved
#

#
# Get paths
#
debug=False
import optparse
import sys, os

from src.definitions import Definition
from src.forcefield import Forcefield
from src.routines import Routines
from src.protein import getPDBFile, readPDB, Protein, Amino, Nucleic
from src.aa import LIG
from src.errors import PDB2PKAError
from pdb2pka import inputgen_pKa

from pdb2pka.pka_routines import smooth
from pdb2pka import pka_help

def usage(x):
    #
    # Print usage guidelines
    #
    print 'Usage: pka.py --ff <forcefield> --lig <ligand in MOL2> --pdie <protein diel cons> --maps <1 for using provided 3D maps; 2 for genereting new maps>'
    print '--xdiel <xdiel maps> --ydiel <ydiel maps> --zdiel <zdiel maps> --kappa <ion-accessibility map> '
    print '--smooth <st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev> <pdbfile>'
    print 'Force field can be amber, charmm and parse'
    print
    return

#
# --------------------------------------------------
#

def startpKa():
    """
        Function for starting pKa script from the command line.

        Returns
            protein:    The protein object as generated by PDB2PQR
            routines:   The routines object as generated by PDB2PQR
            forcefield: The forcefield object as generated by PDB2PQR
    """
    print
    print 'PDB2PQR pKa calculations'
    print

    parser = optparse.OptionParser()

    ##
    ## set optparse options
    ##
    parser.add_option(
        '-v','--verbose',
        dest='verbose',
        action="store_true",
        default=False,
        )
    parser.add_option(
        '--pdie',
        dest='pdie',
        default=8,
        type='int',
        help='<protein dielectric constant>',
        )
    parser.add_option(
        '--sdie',
        dest='sdie',
        default=80,
        type='int',
        help='<solvent dielectric constant>',
        )
    parser.add_option(
        '--ff',
        dest='ff',
        type='choice',
        default='parse',
        choices=("amber","AMBER","charmm","CHARMM","parse","PARSE",),
        help='<force field (amber, charmm, parse)>',
        )
    parser.add_option(
        '--resume',
        dest='resume',
        action="store_true",
        default=False,
        help='resume run from saved state.',
        )
    parser.add_option(
        '--ligand',
        dest='ligand',
        type='str',
        help='<ligand in MOL2 format>',
        )
    parser.add_option(
        '--maps',
        dest='maps',
        default=None,
        type='int',
        help='<1 for using provided 3D maps; 2 for genereting new maps>',
        )
    parser.add_option(
        '--xdiel',
        dest='xdiel',
        default=None,
        type='str',
        help='<xdiel maps>',
        )
    parser.add_option(
        '--ydiel',
        dest='ydiel',
        default=None,
        type='str',
        help='<ydiel maps>',
        )
    parser.add_option(
        '--zdiel',
        dest='zdiel',
        default=None,
        type='str',
        help='<zdiel maps>',
        )
    parser.add_option(
        '--kappa',
        dest='kappa',
        default=None,
        type='str',
        help='<ion-accessibility map>',
        )
    parser.add_option(
        '--smooth',
        dest='sd',
        default=None,
        type='float',
        help='<st.dev [A] of Gaussian smooting of 3D maps at the boundary, bandthwith=3 st.dev>',
        )
    #
    # Cut off energy for calculating non-charged-charged interaction energies
    #
    parser.add_option('--pairene',dest='pairene',type='float',default=1.0,
                      help='Cutoff energy in kT for calculating non charged-charged interaction energies. Default: %default')
    #
    # Options for doing partial calculations
    #
    parser.add_option('--res_energy',
                      dest='desolvation_res',
                      default=[],
                      action='append',
                      type='string',
                      help='Calculate desolvation energy and interaction energy for this residue in its default protonation state. Protonation states can be specified with the --protonation_state argument')
    parser.add_option('--PS_file',dest='PS_file',default='',type='string',action='store',help='Set protonation states according to the pdb2pka protonation state file (option --PS_file)')
    (options,args,) = parser.parse_args()

    ##
    ## parse optparse options
    ##
    ff = options.ff.lower()
    pdie = options.pdie
    verbose = options.verbose
    sdie = options.sdie
    maps = options.maps
    xdiel = options.xdiel
    ydiel = options.ydiel
    zdiel = options.zdiel
    kappa = options.kappa
    sd = options.sd

    #
    # Find the PDB file
    #
    if len(args) != 2:
        parser.error("Usage: pka.py [options] <pdbfile> <output directory>\n")
    input_path = args[0]
    output_path = args[1]

    ligand = None
    if options.ligand is not None:
        try:
            ligand = open(options.ligand, 'rU')
        except IOError:
            print 'Unable to find ligand file %s! Skipping...' % options.ligand

    #Set up the protien object
    #In the standalone version of pdb2pka this is redundent but needed so we emulate the
    #interface needed by pdb2pqr

    pdbfile = getPDBFile(input_path)
    pdblist, errlist = readPDB(pdbfile)
    if len(errlist) != 0 and verbose:
        print "Warning: %s is a non-standard PDB file.\n" %input_path
        print errlist
    #
    # Read the definition file
    #
    myDefinition = Definition()
    #
    #
    # Choose whether to include the ligand or not
    #
    # Add the ligand to the pdb2pqr arrays
    #
    if ligand is None:
        myProtein = Protein(pdblist, myDefinition)
    else:
        from pdb2pka.ligandclean import ligff
        myProtein, _, _ = ligff.initialize(myDefinition, ligand, pdblist, verbose)

    #
    # Call the pre_init function
    #
    return pre_init(protein=myProtein,
                    output_dir=output_path,
                    ff=ff,
                    verbose=verbose,
                    pdie=pdie,
                    sdie=sdie,
                    maps=maps,
                    xdiel=xdiel,
                    ydiel=ydiel,
                    zdiel=zdiel,
                    kappa=kappa,
                    sd=sd,
                    ligand=ligand),options


def pre_init(original_pdb_list=None,
             output_dir=None,
             ff=None,
             verbose=False,
             pdie=8.0,
             sdie=80,
             maps=None,
             xdiel=None,
             ydiel=None,
             zdiel=None,
             kappa=None,
             sd=None,
             ligand=None):
    """This function cleans the PDB and prepares the APBS input file

    Prepares the output folder."""

    #prepare the output directory

    output_dir = os.path.abspath(output_dir)

    try:
        os.makedirs(output_dir)
    except OSError:
        if not os.path.isdir(output_dir):
            raise ValueError('Target directory is a file! Aborting.')

    workspace_dir = os.path.join(output_dir,'workspace')

    try:
        os.makedirs(workspace_dir)
    except OSError:
        if not os.path.isdir(output_dir):
            raise ValueError('Target directory is a file! Aborting.')

    #
    # remove hydrogen atoms
    #

    working_pdb_filename = os.path.join(workspace_dir,'working.pdb')

    pka_help.dump_protein_no_hydrogens(original_pdb_list, working_pdb_filename)
    #
    # Get the PDBfile
    #
    pdbfile = getPDBFile(working_pdb_filename)
    pdblist, errlist = readPDB(pdbfile)

    if verbose:
        print "Beginning PDB2PKA...\n"
    #
    # Read the definition file
    #
    myDefinition = Definition()
    ligand_titratable_groups=None
    #
    #
    # Choose whether to include the ligand or not
    #
    # Add the ligand to the pdb2pqr arrays
    #
    Lig=None
    if ligand is None:
        myProtein = Protein(pdblist, myDefinition)
    else:
        from pdb2pka.ligandclean import ligff
        myProtein, myDefinition, Lig = ligff.initialize(myDefinition, ligand, pdblist, verbose)
    #
    # =======================================================================
    #
    # We have identified the structural elements, now contiue with the setup
    #
    # Print something for some reason?
    #
    if verbose:
        print "Created protein object -"
        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
    #
    # Set up all other routines
    #
    myRoutines = Routines(myProtein, verbose) #myDefinition)
    myRoutines.updateResidueTypes()
    myRoutines.updateSSbridges()
    myRoutines.updateBonds()
    myRoutines.setTermini()
    myRoutines.updateInternalBonds()

    myRoutines.applyNameScheme(Forcefield(ff, myDefinition, None))
    myRoutines.findMissingHeavy()
    myRoutines.addHydrogens()
    myRoutines.debumpProtein()

    #myRoutines.randomizeWaters()
    myProtein.reSerialize()
    #
    # Inject the information on hydrogen conformations in the HYDROGENS.DAT arrays
    # We get this information from ligand_titratable_groups
    #
    from src.hydrogens import hydrogenRoutines
    myRoutines.updateInternalBonds()
    myRoutines.calculateDihedralAngles()
    myhydRoutines = hydrogenRoutines(myRoutines)
    #
    # Here we should inject the info!!
    #
    myhydRoutines.setOptimizeableHydrogens()
    myhydRoutines.initializeFullOptimization()
    myhydRoutines.optimizeHydrogens()
    myhydRoutines.cleanup()
    myRoutines.setStates()

    #
    # Choose the correct forcefield
    #
    myForcefield = Forcefield(ff, myDefinition, None)
    if Lig:
        hitlist, misslist = myRoutines.applyForcefield(myForcefield)
        #
        # Can we get charges for the ligand?
        #
        templist=[]
        ligsuccess=False
        for residue in myProtein.getResidues():
            if isinstance(residue, LIG):
                templist = []
                Lig.make_up2date(residue)
                net_charge=0.0
                print 'Ligand',residue
                print 'Atom\tCharge\tRadius'
                for atom in residue.getAtoms():
                    if atom.mol2charge:
                        atom.ffcharge=atom.mol2charge
                    else:
                        atom.ffcharge = Lig.ligand_props[atom.name]["charge"]
                    #
                    # Find the net charge
                    #
                    net_charge=net_charge+atom.ffcharge
                    #
                    # Assign radius
                    #
                    atom.radius = Lig.ligand_props[atom.name]["radius"]
                    print '%s\t%6.4f\t%6.4f' %(atom.name,atom.ffcharge,atom.radius)
                    if atom in misslist:
                        misslist.pop(misslist.index(atom))
                        templist.append(atom)
                    #
                    # Store the charge and radius in the atom instance for later use
                    # This really should be done in a nicer way, but this will do for now
                    #
                    atom.secret_radius=atom.radius
                    atom.secret_charge=atom.ffcharge
                    #
                    #

                charge = residue.getCharge()
                if abs(charge - round(charge)) > 0.01:
                    # Ligand parameterization failed
                    myProtein.residues.remove(residue)
                    raise Exception('Non-integer charge on ligand: %8.5f' %charge)
                else:
                    ligsuccess = 1
                    # Mark these atoms as hits
                    hitlist = hitlist + templist
                #
                # Print the net charge
                #
                print 'Net charge for ligand %s is: %5.3f' %(residue.name,net_charge)
        #
        # Temporary fix; if ligand was successful, pull all ligands from misslist
        # Not sure if this is needed at all here ...? (Jens wrote this)
        #
        if ligsuccess:
            templist = misslist[:]
            for atom in templist:
                if isinstance(atom.residue, Amino) or isinstance(atom.residue, Nucleic):
                    continue
                misslist.remove(atom)

    if verbose:
        print "Created protein object (after processing myRoutines) -"
        print "\tNumber of residues in protein: %s" % myProtein.numResidues()
        print "\tNumber of atoms in protein   : %s" % myProtein.numAtoms()
    #
    # Create the APBS input file
    #
    import src.psize
    size=src.psize.Psize()

    method=""
    async=0
    split=0

    igen = inputgen_pKa.inputGen(working_pdb_filename)
    #
    # For convenience
    #
    igen.pdie = pdie
    print 'Setting protein dielectric constant to ',igen.pdie
    igen.sdie=sdie
    igen.maps=maps
    if maps==1:
        print "Using dielectric and mobile ion-accessibility function maps in PBE"
        if xdiel:
            igen.xdiel = xdiel
        else:
            raise PDB2PKAError('X dielectric map is missing')
        if ydiel:
            igen.ydiel = ydiel
        else:
            raise PDB2PKAError("Y dielectric map is missing\n")
        if zdiel:
            igen.zdiel = zdiel
        else:
            raise PDB2PKAError("Z dielectric map is missing\n")

        print 'Setting dielectric function maps: %s, %s, %s'%(igen.xdiel,igen.ydiel,igen.zdiel)

        if kappa:
            igen.kappa = kappa
        else:
            raise PDB2PKAError("Mobile ion-accessibility map is missing\n")

        print 'Setting mobile ion-accessibility function map to: ',igen.kappa

        if sd:
            xdiel_smooth, ydiel_smooth, zdiel_smooth = smooth(xdiel,ydiel,zdiel)
            igen.xdiel = xdiel_smooth
            igen.ydiel = ydiel_smooth
            igen.zdiel = zdiel_smooth
    #
    # Return all we need
    #
    return output_dir, myProtein, myRoutines, myForcefield,igen, ligand_titratable_groups, maps, sd

#
# --------------
#

if __name__ == "__main__":

    (output_dir, protein, routines, forcefield,apbs_setup, ligand_titratable_groups, maps, sd), options = startpKa()
    from pdb2pka import pka_routines
    mypkaRoutines = pka_routines.pKaRoutines(protein, routines, forcefield, apbs_setup, output_dir, maps, sd,
                                             restart=not options.resume, pairene=options.pairene)

    print 'Doing full pKa calculation'
    mypkaRoutines.runpKa()