Openbabel error in ReadFile

Next Topic
 
classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Openbabel error in ReadFile

Kohsei Tsuji
This post has NOT been accepted by the mailing list yet.
Hi dear,

I would like to perform covalent docking using the flexible side chain method as recently described in the paper "Covalent docking using autodock: Two-point attractor and flexible side chain methods".

I have used the script provided in http://autodock.scripps.edu/resources/covalentdocking, and

I have followed the procedure in the README text (attached).

But, I got the following openbabel error.

~ /Desktop/adCovalentDockResidue/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
>                   --ligindices 1,2\
>                   --receptor 3upo_protein.pdb\
>                   --residue B:SER222\
>                   --outputfile ligcovalent.pdb
==============================
*** Open Babel Error  in ReadFile
  Cannot read from 3upo_protein.pdb
[ ERROR ] Fail to read molecule "3upo_protein.pdb"

I guess openbabel (ver 2.4.1) (pybel)  has been properly installed in Mac OS X 10.10.5, and imported in preparecovalent.py (attached).

What is the reason of this openbabel error, and how can I fix this error? Any suggestions or comments are highly appreciated.

Thank you.




----------------------------------------------------------------------------------

README.txt

This tutorial describes how to apply the 'flexible side chain' covalent docking as described in the paper:

   "Covalent docking using autodock: Two-point attractor and flexible side chain methods"
   http://www.ncbi.nlm.nih.gov/pubmed/26103917

If you use this method, please cite the paper.
Happy dockings!

==============
REQUIREMENTS
==============
- Python  >=2.6
- OpenBabel Python bindings
   (On Debian/Ubuntu derivatives, it can be installed with "apt-get install python-openbabel")
- Basic experience with docking with AutoDock

==============
INSTALLATION
==============
- Copy the directory 'adcovalent' in a location that is accessible from the user (i.e., '~/local/adcovalent').

==============
USAGE
==============

The preparation code is used to generate the geometry of the ligand covalently bound to the
residue. Once the ligand initial coordinate is generated, the PDBQT files of the receptor and
flexible ligand need to be generated using prepare_flexreceptor4.py script in AutoDockTools.

The directory '3upo_test' contains an example case of a ligand ready to processed.
The following examples assume a terminal opened in that location.


1. Generate the alignment
--------------------------
The input ligand structure must be modeled with a portion of the alkylated residue already present,
i.e. for a ligand bound to a serine:

         LIGAND-OG-CB-CA-RESIDUE  [MOL1]

the ligand structure to process should be:
 
         LIGAND-O-C [MOL2]

where C and O are the two atoms that are going to be used for the alignment. More in general:

        0. Initial configuration

                    B
                   /
           [L]----A      I---J--[R]
                        /
                       X

        1. Translate [L] to overlap A->I

                    B
                   /
           [L]----A  ->  I---J--[R]
                        /
                       X

                          B
                         /
                 [L]----AI---J--[R]
                        /
                       X

        2. Rotate [L] around the normal of plane (B,A,J) to overlap B->J
             
                [L]
                  \
                   AI--BJ--[R]
                   /
                  X

The alignment could be skipped if the ligand is already attached to the receptor structure
(i.e., a re-docking experiment from a PDB is going to be perfoemd).

To perform the alignment, a receptor file and a residue name (Chain:ResNum) must
be specified, then a ligand file and an alignment criterion.

The ligand alignment can be defined in two ways. The first is by specifying
the atom indices in the ligand file (i.e., first and second atoms):

    python  ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
                  --ligindices 1,2\
                  --receptor 3upo_protein.pdb\
                  --residue B:SER222\
                  --outputfile ligcovalent.pdb

or, in a more general approach, by specifying a SMARTS pattern and the indices
of the alignment atoms in the pattern itself:

    python  ~/local/adcovalent/prepareCovalent.py --ligand ligand.mol2 \
                  --ligindices 4,3\
                  --ligsmart "C(=O)-O-C"\
                  --receptor 3upo_protein.pdb\
                  --residue B:SER222\
                  --outputfile ligcovalent.pdb

2. Generate PDBQT files
--------------------------
The standard PDBQT files for AutoDock need to be generated for both the receptor structure. The standard ADT scripts
are going to be used, assuming that $MGLROOT is the path where ADT has been installed.

*** NOTE: if a covalent ligand is already bound to the residue in the protein structure, it must be removed! ***

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r 3upo_protein.pdb -A hydrogens
       [ this generates the file "3upo_protein.pdbqt" ]

and the covalent ligand that has been aligned:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_receptor4.py -r ligcovalent.pdb
       [ this generates the file "ligcovalent.pdbqt" ]


3. Generate flexible/rigid PDBQT files
---------------------------------------
The PDBQT files are going to be used to generate the rigid and flexible components that are going to be used
for the docking.
First, the receptor is processed to extract the rigid part by specifying which residue is going to be made flexible:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r 3upo_protein.pdbqt -s 3upo_protein:B:SER222
       [ this generates the files "3upo_protein_rigid.pdbqt" and  "3upo_protein_flex.pdbqt" ]

The same thing is going to be done for the processed ligand:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_flexreceptor4.py -r ligcovalent.pdbqt -s ligcovalent:B:SER222
       [ this generates the files "ligcovalent_rigid.pdbqt" and  "ligcovalent_flex.pdbqt" ]

NOTE: if the ligand has not been prepared in step 1, make sure that all atoms in the ligand have the poper residue id
corresponding to the residue modified (i.e., Ser222A).


4. Generate GPF and DPF files
------------------------------
The parameter files for the actual calculation are going to be generated. The GPF for AutoGrid is generated with:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_gpf4.py -r 3upo_protein_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -y -I 20\
                -o 3upo_priotein.gpf

The DPF for AutoDock is generated with:

    $MGLROOT/bin/pythonsh $MGLROOT/MGLToolsPckgs/AutoDockTools/Utilities24/prepare_dpf4.py -r 3upo_protein_rigid.pdbqt\
                -x ligcovalent_flex.pdbqt\
                -l ligcovalent_flex.pdbqt\
                -o ligcovalent_3upo_protein.dpf\
                -p move='empty'

This command instructs the script to dock the covalent ligand as a flexible residue and ignore any 'true' ligand ("move='empty'"). To do this,
an empty file must be created, and it can be done with the following Unix command:
   
    touch empty

Also, in order to include the ligand atoms in the calculation of RMSD and clustering, the following parameters need to be added:

    rmsdatoms all

which will include flexible residue (and ligand) atoms in the RMSD calculation.
 Finally, the DPF file must be manually edited to set the appropriate energy model so that the docking score corresponds to the interaction between
the flexible residue (the ligand) and the rigid receptor. For this, the entry:
   
    unbound_model bound

must be replaced with:

    unbound_energy 0.0


5. Run AutoGrid and AutoDock
------------------------------
Both programs can be executed using the standard procedure:

    autogrid4 -p 3upo_priotein.gpf -l 3upo_priotein.glg
    autodock4 -p ligcovalent_3upo_protein.dpf -l  ligcovalent_3upo_protein.dlg



The output generated at each step of the process is available in the directory '3upo_test/output'.










------------------------------------------------------------------------------------------------------------
prepareCovalent.py

#!/usr/bin/env python
#
# Covalent Docking preparation for AutoDock
#
# v.1.0c  Stefano Forli
#
# Copyright 2014, Molecular Graphics Lab
#     The Scripps Research Institute
#        _  
#       (,)  T  h e
#      _/
#     (.)    S  c r i p p s
#      '\_
#       (,)  R  e s e a r c h
#      ./'
#     ( )    I  n s t i t u t e
#      "
#
#################################################################################
#
#     This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>
#

import sys, os, math
import pybel
ob = pybel.ob
import argparse
import numpy as np
from math import sqrt, sin, cos, acos, degrees
import ResidueProfiler
from CopyMol import MoleculeDuplicator


# TODO improve the method to build the bond directly, calculating
#      bond lenght and its placement

#### Debug functions
def makePdb(coord, keyw = "ATOM  ", at_index = 1, res_index = 1, atype = 'X', elem = None,
            res = "CNT", chain  ="Z", bfactor = 10,pcharge = 0.0):
    if not elem: elem = atype
    # padding bfactor
    bfactor = "%2.2f" % bfactor
    if len(bfactor.split(".")[0]) == 1:
        bfactor = " "+bfactor
    if len(atype) == 1:
        atype = atype + " "
    #atom = "%s%5d  %2s  %3s %1s%4d    %8.3f%8.3f%8.3f  1.00 %02.2f  %8.3f %1s" % (keyw,
    # SOURCE: http://cupnet.net/pdb-format/
    "%-6s%5d %4s%1s%3s %1s%4d%1s   %8.3f%8.3f%8.3f%6.2f%6.2f          %2s%2s"
    atom = "%s%5d %2s   %3s %1s%4d    %8.3f%8.3f%8.3f  1.00 %s  %8.3f %1s" % (keyw,
            at_index, elem, res, chain, res_index,
            coord[0], coord[1], coord[2], bfactor, pcharge, atype)
    return atom

def writeList(filename, inlist, mode = 'w', addNewLine = False):
    if addNewLine: nl = "\n"
    else: nl = ""
    fp = open(filename, mode)
    for i in inlist:
        fp.write(str(i)+nl)
    fp.close()


def vector(p1 , p2 = None, norm = 0): # TODO use Numpy?
    if not p2 == None:
        vec = np.array([p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]],'f')
    else:
        vec = np.array([p1[0], p1[1], p1[2] ], 'f' )

    if norm:
        return normalize(vec)
    else:
        return vec

def dot(vector1, vector2):  # TODO remove and use Numpy
    dot_product = 0.
    for i in range(0, len(vector1)):
        dot_product += (vector1[i] * vector2[i])
    return dot_product

def vecAngle(v1, v2, rad=1): # TODO remove and use Numpy?
    angle = dot(normalize(v1), normalize(v2))
    try:
        if rad:
            return acos(angle)
        else:
            return degrees(acos(angle))
    except:
        print "#vecAngle> CHECK TrottNormalization"
        return 0

def norm(A): # TODO use Numpy
        "Return vector norm"
        return math.sqrt(sum(A*A))

def normalize(A): # TODO use Numpy
        "Normalize the Vector"
        return A/norm(A)



def calcPlane(p1, p2, p3):
    # returns the plane containing the 3 input points
    v12 = vector(p1,p2)
    v13 = vector(p3,p2)
    return normalize(np.cross(v12, v13))



def rotatePoint(pt,m,ax):
    """
    Rotate a point applied in m around a pivot ax ?


    pt = point that is rotated
    ax = vector around wich rotation is performed
    """
    # point
    x=pt[0]
    y=pt[1]
    z=pt[2]
    # rotation pivot
    u=ax[0]
    v=ax[1]
    w=ax[2]
    ux=u*x
    uy=u*y
    uz=u*z
    vx=v*x
    vy=v*y
    vz=v*z
    wx=w*x
    wy=w*y
    wz=w*z
    sa=sin(ax[3])
    ca=cos(ax[3])
    p0 =(u*(ux+vy+wz)+(x*(v*v+w*w)-u*(vy+wz))*ca+(-wy+vz)*sa)+ m[0]
    p1=(v*(ux+vy+wz)+(y*(u*u+w*w)-v*(ux+wz))*ca+(wx-uz)*sa)+ m[1]
    p2=(w*(ux+vy+wz)+(z*(u*u+v*v)-w*(ux+vy))*ca+(-vx+uy)*sa)+ m[2]
    return np.array([ p0, p1, p2])




class CovalentDockingMaker:
    """ This class generate the input for covalent dockings.
        Required input are :

            - ligand structure (any OB-supported 3D format)
            - receptor structure (any OB-supported 3D format)
       
        Ligand and receptor structures should share at least 2 atoms
        to be overlapped, e.g.:

                Cys-Ca-Cb-S     Cb-S-Lig

        Alignment can be performed using one of the following criteria:

            - indices: two indices pairs for ligand and receptor

            - SMARTS: a pattern and two atom indices must be defied for the ligand;
                      if nothing is defined for the receptor (default), the residue
                      type will be used to pick its strongst nucleophyle atom,
                      i.e. Cys-SG, Ser/Thr-OG, Lys-NZ.

            - bond: a bond for each structure is provided and indices
                    calculated from them; bonds must be terminal bonds
                    (not counting H's)
                    (not working yet)
    """
    def __init__(self,
            lig,    # ligand openbabel molecule
            rec,    # receptor openbabel molecule
            resName,# residue id to be modifled (i.e. A:THR276)

            smartsLig=None, smartsRec=None,   # use SMARTS patterns
                                              # patterns should be aligned
            smartsIdxLig=(), smartsIdxRec=(), # align using atomIdx from SMARTS patterns

            bondIdxLig=None, bondIdxRec=None, # use bonds

            indicesLig=None, indicesRec=None, # use atom indices in both

            genFullRec = False,               # add the covalent ligand to the receptor
                                              # as modified residue; this is *NOT*
                                              # recommended, because it can lead to
                                              # spurious bonds between the mod res
                                              # and the rest of the receptor
            genLigand = False,                  


            auto=True,                        # start automatically
            debug = False,
            verbose=False):

        self.verbose = verbose
        self.debug = debug
        # LIGAND
        dup = MoleculeDuplicator(lig)
        self.lig = dup.getCopy()
        # parameters(REQUIRED)
        self.smartsLig = smartsLig       # smarts
        self.smartsIdxLig = smartsIdxLig # smarts indices
        self.bondIdxLig = bondIdxLig     # bonds
        self.indicesLig = indicesLig     # indices (direct)

        # RECEPTOR
        self.fullRec = rec
        self.resName = resName.upper()
        # parameters (OPTIONAL)
        self.smartsRec = smartsRec        # smarts          
        self.smartsIdxRec = smartsIdxRec  # smarts indices
        self.bondIdxRec = bondIdxRec      # bonds
        self.indicesRec = indicesRec      # indices (direct)
        # collect receptor information
        self.genFullRec = genFullRec
        self.setMode()
        self.initResidue()
        if auto:
            self.process()


    def vprint(self, msg):
        """ """
        if not self.verbose:
            return
        print("VERBOSE " + msg)
   

    def setMode(self):  
        """ """
        # ligand mode
        if self.indicesLig:
            self._ligmode = 'idx'
        elif self.smartsLig:
            self._ligmode = 'smarts'
        elif self.bondIdxLig:
            self._ligmode = 'bond'
        self.vprint("[setMode] ligand alignment mode: %s" % self._ligmode)
        # receptor mode
        if self.indicesRec:
            self._recmode = 'idx'
        elif self.smartsRec:
            self._recmode = 'smarts'
        elif self.bondIdxRec:
            self._recmode = 'bond'
        else:
            self._recmode = 'smarts'
            self.vprint("[setMode] no mode specified-> switching to default")
        self.vprint("[setMode] receptor alignment mode: %s" % self._recmode)
   
    def initResidue(self):
        """ initialize residue information"""
        string = self.resName
        chain, res = string.split(":")
        self.resInfo = { 'chain' : chain,
                         'name'  : res[0:3],
                         'num': int(res[3:]),
                       }
        if self.debug: print "PROCESSING RESIDUE", self.resName
        self.residueProfile()

    def process(self):
        """ do the stuff"""
        if self.verbose: print "Aligning ligand on residue...", self.resName
        self.align()
        self.cleanup()

    def residueProfile(self):
        """ tries to recognize the receptor residue
            and rename atoms to facilitate Calpha-Cbeta identification
        """
        self.receptorProfiler = ResidueProfiler.AminoAcidProfiler(self.fullRec,
                resId=self.resName, setLabels=True,
                auto=True, debug=self.debug)
        self.resType = self.receptorProfiler.residues[self.resName]['type']
        if self.verbose:
            print "Selected residue recognized as:", self.resType.upper()
        dup = MoleculeDuplicator(self.fullRec, resList=[self.resName])
        self.modifiedRes = dup.getCopy()
        if self._recmode == 'idx':
            org = self.indicesRec[:]
            self.indicesRec = [ self.modifiedRes._numbering[x] for x in self.indicesRec]
            self.vprint( ('converting full receptor indices [%s] to '
                 'residue copy indices [%s]' % (org, self.indicesRec)))
        self.setResidueParms()

    def setResidueParms(self):
        """ populate automatically SMARTS and indices for
            the residue if nothing was specified already
        """
        if not self._recmode == 'smarts':
            msg = "[setResidueParm] receptor alignment mode is (%s), skipping"
            self.vprint (msg % self._recmode)
            return
        if self.smartsRec:
            msg = ("[setResidueParms] receptor SMARTS already defined [%s], skipping")
            self.vprint(msg % self.smartsRec)
            return
        self._residueParms = {'cys' : { 'patternIdx' : [3,2],     # when using SMARTS patterns from ResProfiler
                                 'atNames': ['SG', 'CB'],   # when looking for [setMode] atom names in the structure
                                },
                        'ser' : { 'patternIdx' : [3,2],
                                 'atNames' : ['OG', 'CB'],
                                },
                        'thr' : { 'patternIdx' : [4,2],
                                 'atNames' :  ['OG', 'CB'],
                                },
                        'lys' : { 'patternIdx': [6,5],
                                 'atNames': ['NZ','CE'],
                                },
                        }
        resInfo = self._residueParms[self.resType]
        self.smartsRec = self.receptorProfiler.getResPattern(self.resType)
        self.smartsIdxRec = resInfo['patternIdx']
        if self.verbose:
            print "Auto-parameters for residue (%s)" % (self.resType)
            print "    SMARTS pattern : %s" % self.smartsRec
            print "    atom indices   : %s" % self.smartsIdxRec
            print "    atom names     : %s" % resInfo['atNames']

    def align(self):
        """ where the magic happens... perform the alignment"""
        self.populateIndices()
        if not self.ready:
            print "Molecule not ready... Aborting."
            return
        aligner = VectorMolAligner(lig=self.lig, rec=self.modifiedRes,
            ligIndices=self.indicesLig, recIndices=self.indicesRec,
            verbose=self.verbose, debug=self.debug)

    def populateIndices(self):
        """ populate indices pairs accordingly to the modes"""
        # XXX add a check for missing matches
        # ligand indices:
        self.ready = True
        if self._ligmode == 'smarts':
            smarter = IndicesFromSMARTS(mol = self.lig,
                pattern = self.smartsLig, indices = self.smartsIdxLig,
                firstOnly = False, debug = self.debug)
            if not smarter.matches:
                print "*** WARNING *** Pattern not found in the ligand"
                self.ready = False
                return
            self.indicesLig = smarter.matches[0]
        elif self._ligmode == 'bond':
            pass
        else:
            self.vprint("[align] Using provided indices for ligand")
        # receptor indices
        if self._recmode == 'smarts':
            smarter = IndicesFromSMARTS(mol = self.modifiedRes,
                pattern = self.smartsRec, indices = self.smartsIdxRec,
                debug = self.debug)
            if not smarter.matches:
                print "*** WARNING *** Pattern not found in the receptor"
                self.ready = False
                return
            self.indicesRec = smarter.matches[0]
        elif self._recmode == 'bond':
            pass
        else:
            self.vprint("[align] Using provided indices for receptor")
       
    def cleanup(self):
        """ Lig: remove from A on...
            Rec: remove from I on...
            merge from CA on to the residue
        """
        if not self.ready:
            return
        self.cleanLigand()
        self.cleanReceptor()
        self.cleanModRes()
        if self.genFullRec:
            self.generateFullRec()

    def cleanLigand(self):
        """ remove atoms A and B from the ligand"""
        # delete ligand atoms used to overlap
        for a in [ self.lig.GetAtom(i) for i in self.indicesLig ]:
            self.lig.DeleteAtom(a)
        # delete hydrogens
        self.lig.DeleteHydrogens()
        # create new molecule (XXX maybe it should be moved at the beginning?
        # before deleting atoms
        dup = MoleculeDuplicator(self.lig, name='Covalent docking of...')
        self.covalentLigand = dup.getCopy()
        # XXX BUG this triggers the creation of a residue "LIG" for ligands
        # with no residues
        self.covalentLigand.GetAtom(1).GetResidue()
        covalentRes = self.covalentLigand.GetResidue(0)
        #if covalentRes:
        covalentRes.SetName(self.resType.upper())
        # renaming ligand atoms 'C -> CX' to avoid issues later (ADT...)
        for atom in ob.OBResidueAtomIter(covalentRes):
            oldIdFull = covalentRes.GetAtomID(atom)
            # keep the element letter, throw the rest...
            oldId = oldIdFull[0]
            covalentRes.SetAtomID(atom, oldId+'X')

        if self.debug:
            pybel.Molecule(self.lig).write('pdb', 'DEBUG_cleanLig.pdb', overwrite=1)

    def cleanReceptor(self):
        """ use PDB atom labels to decide atoms to be kept
        """
        # XXX CORRECT XXX
        #accepted = ['CA', 'N'] + self.receptorProfiler.getResLabels(self.resType,
        # backbone=False)
        accepted = ['CA'] + self.receptorProfiler.getResLabels(self.resType,
            backbone=True)
        if self.debug:
            print "Accepted residue atom names:" , accepted
        covalentRes = self.covalentLigand.GetResidue(0)
        for atom in ob.OBMolAtomIter(self.modifiedRes):
            aName = self.getAtomName(atom, True)
            if self.debug: print "Scanning atom [%s]" % aName
            if aName in accepted:
                if self.debug: print "[ accepted ]"
                newName = ' {0:3}'.format(aName)
                newAtom = self.covalentLigand.NewAtom()
                newAtom.Duplicate(atom)
                covalentRes.AddAtom(newAtom)
                # residue atoms will be ATOM
                covalentRes.SetAtomID(newAtom, newName)
        if self.debug:
            mol = pybel.Molecule(self.covalentLigand)
            mol.write('pdb', 'final.pdb', overwrite=1)
   
    def cleanModRes(self, setHet=True, connect=True, addH=True):
        """ perform last clean-up operations:
                conect : autodetect bonds
                addH   : add hydrogens
                setHet : set all atom records as HETATM
        """
        # autodetect bonds
        if connect:
            self.covalentLigand.ConnectTheDots()
        # add hydrogens
        if addH:
            self.addHydrogens()
        # set all atoms to be HETATM records
        covalentRes = [ x for x in ob.OBResidueIter(self.covalentLigand)][0]
        if setHet:
            for a in ob.OBMolAtomIter(self.covalentLigand):
                covalentRes.SetHetAtom(a, True)

        covalentRes.SetName(self.resType.upper())
        chain, res = self.resName.split(':',1)
        print "CHAIN", chain
        rnum = res[3:]
        covalentRes.SetNum(int(rnum))
        covalentRes.SetChain(chain)

    def generateFullRec(self):
        """ replace the original residue with the
            modified one
        """
        rName = self.resInfo['name']
        rChain = self.resInfo['chain']
        rNum = self.resInfo['num']
        for res in ob.OBResidueIter(self.rec):
            if (res.GetChain() == rChain):
                if (res.GetNum() == rNum):
                        self.rec.DeleteResidue(res)
                        if self.verbose:
                            print "deleted residue [%s:%s%s]" % (rChain, rName, rNum)
                        break
        # add modified residue
        #covalentRes = [ x for x in ob.OBResidueIter(self.covalentLigand)][0]
        #self.rec.AddResidue(covalentRes)


    def addHydrogens(self):
        """ manages all the gymnastics for adding hydrogens at pH 7.4 (DEFAULT)"""
        #self._temporaryPHFix()
        self.covalentLigand.DeleteHydrogens()
        self.covalentLigand.UnsetFlag(ob.OB_PH_CORRECTED_MOL)
        for a in ob.OBMolAtomIter(self.covalentLigand):
            a.SetFormalCharge(0)
        self.covalentLigand.SetAutomaticFormalCharge(True)
        #self.covalentLigand.PerceiveBondOrders()
        self.covalentLigand.AddHydrogens(False, True )#, pH)

        #self.covalentLigand.AddHydrogens(False,True)

    def _temporaryPHFix(self):
        """ this is a nasty workaround for an even nastier bug in OB:
            http://sourceforge.net/p/openbabel/bugs/710/
        """
        t = 'mol2'
        n = '.tmpmol'
        m = pybel.Molecule(self.covalentLigand)
        m.write(t, n, overwrite=1)
        self.covalentLigand = pybel.readfile(t, n).next().OBMol
        os.remove(n)
       
    def getAtomName(self, atom, strip=False):
        """ return the PDB atom name """
        res = atom.GetResidue()
        if strip:
            return res.GetAtomID(atom).strip()
        return res.GetAtomID(atom)


    def writeCovalent(self, outFname, _format='pdb'):
        """ _format could be any OB format
       
            NOTE: test mol2 to avoid spurious bonds?
        """
        if self.genFullRec:
            mol = self.rec
        else:
            mol = self.covalentLigand

        if _format == 'pdbqt':
            self._writePdbqt(outFname+'.pdbqt')
        else:
            #fname = "%s.%s" % ( outFname, _format.lower() )
            fname = outFname
            conv = ob.OBConversion()
            conv.SetOutFormat(_format)
            conv.WriteFile(mol, fname)
            if self.verbose: print "Written: %s" % (fname)
           
    def _writePdbqt(self, outFname):
        """ not necessary, for now, we will use prepare_flexres code later"""
        if self.genFullRec:
            mol = self.rec
        else:
            mol = self.covalentLigand


class IndicesFromSMARTS:
    """ Given a SMARTS pattern, return the atom indices matching
        the specified atoms in the pattern
    """
    def __init__(self, mol, pattern, indices=(), firstOnly = True, debug=False):
        """ mol : OBMol()
            pattern = SMARTS string
            [ indices: indices of the SMARTS atoms that need to be identified ]
            [ firstOnly: by default only first match is returned; multiple are possible]
        """
        self.debug = debug
        self.mol = mol
        if self.debug:
            print "IndicesFromSMARTS> MOLECULE NAME", self.mol.GetTitle()
            print "IndicesFromSMARTS> MOLECULE ", self.mol
            print "PATTERN", pattern
        self.pattern = pattern
        self.indices = indices
        if self.indices == ():
            self.indices = (0,1)
            print "pointMatch> WARNING! missing indices, using default (0,1)"
        self.matches = []
        self.pointMatch()

    def pointMatch(self):
        """ correct the indexing from the SMARTS pattern indexing to
            the 0-based used to retrieve atoms
        """
        result = self.findSmarts()
        if not len(result):
            print "IndexFromSMARTS: NO MATCH FOUND! troubles ahead..."
            return
        if self.debug:
            print "[IndicesFromSMARTS] >>> raw out", result
        if len(result) > 1:
            print "TMP: MULTIPLE MATCHES FOUND!", len(result)
        i,j = self.indices
        if self.debug:
            print "[IndicesFromSMARTS]> extracting indices [%d,%d] from smarts[%s]" % (i,j,result)
        self.matches = [ (r[i],r[j]) for r in result ]

        if self.debug:
            print "FOUND!", self.matches
            for a in ob.OBMolAtomIter(self.mol):
                if a.GetIdx() in self.matches[0]:
                    res = a.GetResidue()
                    print "RESNAME", res.GetName()
        return self.matches

    def findSmarts(self):
        """ OB SMARTS matcher """
        obpat = ob.OBSmartsPattern()
        obpat.Init(self.pattern)
        obpat.Match(self.mol)
        if self.debug:
            print "PROCESSING MOLECULE", self.getSmi()
        return [ x for x in obpat.GetUMapList() ]
   
    def getSmi(self):
        """ """
        m = pybel.Molecule(self.mol)
        return m.__str__()
       


# XXX INACCURATE! PROBABLY TO BE DROPPED?
class IndicesFromBond:
    """ identify atom indices involved in a bond"""
    def __init__(self, mol, bond):
        """ The requirement
        XXX the requirement is that there must be a unique atom type that's
       
        common?  -X-S  S-Y-

        Some smartness could be used to guess what is the 'terminal'
        direction to infer the proper orientation of the molecules

        """
        pass


class VectorMolAligner:
    """ """
    def __init__(self, lig, rec, ligIndices=(), recIndices=(), auto=True,
        verbose=False, debug=False):
        """
        Given two OB Molecules and two pairs of atom indices
        (A,B) and (I,J) perform the covalent alignment

        0. Initial configuration

                    B
                   /
           [L]----A      I---J--[R]
                        /
                       X


        1. Translate [L] to overlap A->I


                    B
                   /
           [L]----A  ->  I---J--[R]
                        /
                       X

                          B
                         /
                 [L]----AI---J--[R]
                        /
                       X


        2. Rotate [L] around the normal of plane (B,A,J) to overlap B->J
             
                [L]
                  \
                   AI--BJ--[R]
                   /
                  X


        Indexing provided must be as following:
                ligIndices = (B, A)
                recIndices = (I, J)

        """
        self.debug = debug
        self.verbose = verbose
        self.lig = lig
        self.rec = rec
        self.ligIndices = ligIndices
        self.recIndices = recIndices
        if auto:
            self.align()

    def align(self):
        """ perform the alignment based on vectors
        """
        self.getVectors()
        self.translate()
        self.rotate()


    def getVectors(self):
        """ convert all lig and rec indices to vectors

        """
        # XXX ADAPT THIS TO WORK WITH MULTIPLE MATCHES
        if self.debug: print "GETTING LIGAND VECTORS..."
        self._ligB = self._idxToVec(self.lig, self.ligIndices[0])
        self._ligA = self._idxToVec(self.lig, self.ligIndices[1])


        if self.debug: print "GETTING RECEPTOR VECTORS..."
        self._recI = self._idxToVec(self.rec, self.recIndices[0], translate=0) #  , reverse=)
        self._recJ = self._idxToVec(self.rec, self.recIndices[1], translate=0) # , reverse=True)
        #self._recI = self._idxToVec(self.rec, self.recIndices[0], translate=True, reverse=True)
        #self._recJ = self._idxToVec(self.rec, self.recIndices[1], translate=True, reverse=True)


        #self._recI = self._idxToVec(self.rec, self.recIndices[0], translate=True) #  , reverse=)
        #self._recJ = self._idxToVec(self.rec, self.recIndices[1], translate=True) # , reverse=True)

        if self.debug:
            c = 1
            buff = ['MODEL']
            buff.append( makePdb(coord=self._ligB, at_index=c, atype='P'))
            c+=1
            buff.append( makePdb(coord=self._ligA, at_index=c, atype='P'))
            c+=1
            buff.append( makePdb(coord=self._recI, at_index=c, atype='S'))
            c+=1
            buff.append( makePdb(coord=self._recJ, at_index=c, atype='S'))
            buff.append('ENDMDL')
            writeList('DEBUG_getVectors.pdb', buff, addNewLine=1, mode='a')

    def _idxToVec(self, mol, atIdx, convert=True,translate=False):
        """ return the vector of the atom at index atIdx
            if convert requested, return Numpy array instead of OB.Vector3
        """
        #print "CALLED BY", mol, atIdx, convert, translate
        atIdx = int(atIdx)
        if translate:
            # DuplicateMolecule obj
            atom = mol.GetAtom(atIdx, translate=True)
        else:
            # OBMol object
            atom = mol.GetAtom(atIdx)
        try:
            vec = atom.GetVector()
        except:
            print "*** ERROR:", sys.exc_info()[1]
            print "*** REQUESTED", atIdx
            print "*** FOUND",
            for a in ob.OBMolAtomIter(mol):
                print a.GetIdx(),
            print " "
            return
        if not convert:
            return vec
        return self.obVecToNumpy( vec )


    def translate(self):
        """ perform the 1. translation step """
        if self.verbose: print "Translation...",
        translation = vector( self._ligA, self._recI)
        self._translateMol(self.lig, translation)
        # update vector coords after translation
        self.getVectors()
        if self.debug:
            pybel.Molecule(self.lig).write('pdb',
                'DEBUG_translated.pdb', overwrite=True)
        if self.verbose: print "[ DONE ]"

    def rotate(self):
        """ perform the 2. rotation step"""
        import math
        if self.verbose: print "Rotation...",
        ## rotation
        p1 = self._ligB
        p2 = self._ligA
        p3 = self._recJ
        p4 = self._recI
        v1 = vector(p1,p2)
        v2 = vector(p3,p2)
        angle = vecAngle(v1, v2) # rad
        center = vector(p2)
        normVec = calcPlane( p1, p2, p3)
        if angle < 0.09:
            if self.verbose:
                print ( '[ skipping rotation, angle near-zero (%1.3f ), '
                        'ligand possibly already in place' % (angle) )
            #return
            pass
        if self.debug:
            print "Rotation angle:", math.degrees(angle)
            pdb1 = makePdb(coord=p1, at_index = 1, atype ='S')
            pdb2 = makePdb(coord=p2, at_index = 2, atype ='S')
            pdb3 = makePdb(coord=p3, at_index = 3, atype ='S')
            pdb4 = makePdb(coord=normVec+p2, at_index = 4, atype ='N')
            self._testrot(p1, p2, normVec) # writes rotor.pdb
            writeList('DEBUG_plane.pdb', [pdb1, pdb2, pdb3, pdb4], addNewLine=1)
        rotor = ( normVec[0], normVec[1], normVec[2], angle )
        self._rotateMol( self.lig, p2, rotor)
        if self.debug:
            pybel.Molecule(self.lig).write('pdb', 'DEBUG_rotated.pdb', overwrite=True)
        if self.verbose: print "[ DONE ]"
        return


    def _translateMol(self, mol, v):
        """ translate the molecule using vector3"""
        for i in range(mol.NumAtoms()):
            atom = mol.GetAtom(i+1)
            pre = (atom.GetX(), atom.GetY(), atom.GetZ())
            post = (pre[0]+v[0], pre[1]+v[1], pre[2]+v[2])
            obvec = ob.vector3()
            obvec.Set(post[0], post[1], post[2])
            atom.SetVector(obvec)
        return

    def _rotateMol(self, mol, center, rotor):
        """ rotate molecule atoms applying vector around center"""
        for a in ob.OBMolAtomIter(mol):
            coord = self.obVecToNumpy( a.GetVector() )
            vec = vector(center, coord)
            newCoord = rotatePoint( vec, center, rotor)
            x,y,z = self.vecToDouble(newCoord)
            vec = ob.vector3()
            vec.Set(x,y,z)
            a.SetVector(vec)

    def _testrot(self, point, center, norm):
        """ debug function to depict the rotor applied in the mol rotation"""
        buff = []
        buff = [  makePdb((0.,0.,0.), at_index =1, atype='S') ]
        buff.append(  makePdb((1.8,0.,0.), at_index =2, atype='S') )
        buff.append(  makePdb((0.,1.8,0.), at_index =3, atype='S') )
        buff.append(  makePdb((0.,0.,1.8), at_index =4, atype='S') )
        buff.append(  makePdb(point, at_index =5, atype='S') )
        c = 6
        for a in range(0, 360, 30):
            a = math.radians(a)
            rotor = ( norm[0], norm[1], norm[2], a)
            coord = vector(point,center)
            new = rotatePoint( coord, center, rotor)
            atom = makePdb(new, at_index =c)
            buff.append(atom)
        writeList('DEBUG_rotor.pdb', buff, addNewLine=1)
       
       
    def obVecToNumpy(self, ObVec):
        """ """
        return np.array( ( ObVec.GetX(), ObVec.GetY(), ObVec.GetZ() ), 'f')

    def vecToDouble(self, vector):
        """ """
        return float(vector[0]), float(vector[1]), float(vector[2])
       
class CovalentReaction:
    def __init__(self):
        """ simmeglio """
        """ IT SHOULD be used to attach the next group to the reaction
        to be used for the alignment """
        pass

           


# ########################


class CovalentDockingMaster:
    """ manage settings for the CovalentDocking object"""

    def __init__(self, debug=False, verbose=False):
        """ """
        self.debug = debug
        self.verbose = verbose
        self.initOpts()
        self.parseOpts()
        self.start()

    def getFnameExt(self,string):
        """ """
        name, ext = os.path.splitext(string)
        ext = ext[1:]
        return name, ext

    def initOpts(self):
        """ """
        self.initDefaults()
        self.initDocs()
        self.opts = {

        '--ligand' : {'help': 'ligand file (OPTIONAL)',
            'action':'store', #'required':True,
            'type':str},

        '--receptor' : {'help':'receptor file (REQUIRED)',
            'action':'store', 'required':True, 'type':str},

        '--residue' : {'help': ('residue to modify; the format is "chain:resNUM (es. "B:THR276")'
                               '; multiple residues can be specified by repeating --residue (note: '
                                'case insensitive)'), 'action':'append', 'default':[], 'type':str},

###        '--residuelist' : {'help' : 'modify residues read from file (one per line)',
###                            'action':'store', 'metavar': 'FILE', 'default':None},
###
###        '--covresidue' : {'help':'prepare modified residue for the covalent residue mode (default)',
###            'default': True, 'action': 'store_true'},
###
###        '--covmap' : {'help': 'calculate coordinates for Z-map covalent method',
###            'default': False, 'action': 'store_true'},
###
###
        '--ligsmarts' : {'help':('ligand SMARTS pattern; by default, the first two atoms of '
                         'the pattern are used; use --ligindices to specify different atoms'),
            'action':'store', 'type':str, 'metavar':'SMARTS', 'default':None},

        '--ligindices' : {'help':('indices of ligand atoms to use for alignment (i.e X-Y-lig)'
            'starting from 1; by default, indices refer to atoms in the ligand file; if --ligsmarts '
            'is used, indices specify atoms in the SMARTS pattern'),
            'action':'store', 'type':str, 'metavar':'X,Y'},

        '--recsmarts' : {'help':('receptor SMARTS pattern; by default, the first two atoms of '
                         'the pattern are used; use --recindices to specify different atoms'),
            'action':'store', 'type':str, 'metavar':'SMARTS', 'default':None},

        '--recindices' : {'help':('indices of receptor atoms to use for alignment (i.e. I-J-rec) '
            'starting from 1; by default, indices refer to atoms in the ligand file; if --recsmarts '
            'is used, indices specify atoms in the SMARTS pattern'),
            'action':'store', 'type':str, 'metavar':'I,J'},

        '--genfullrec' :  {'help':('generate combined receptor-ligand structure; '
                       'by default, only the modified residue (including the ligand) is written '
                        '(default : %s)' % self.default_genfullrec ) ,
                'action':'store_true' },

        '--outputfile' :  {'help':('save the output in the file specified'
                       '(default : use input file'),
            'action':'store', 'type':str, 'metavar':'FILENAME'},

###        '--generatelig' : {'help' : ('generate structure of flexible covalent ligand for AutoDock'),
###                       'action':'store_true', 'default':self.default_flex},
###
        '--generaterec' : {'help' : ('generate structure')},

        '--verbose':{'help': ('enable verbose mode'), 'action':'store_true',
                    'default':False},

        '--log' : {'help' : ('write log file'), 'action':'store', 'metavar' : 'LOGFILE'},

        '--help_advanced':{'help': ('show advanced help and documentation'),
            'action': 'store_true', 'default':False}
            }

        self.groups_order = ['INPUT STRUCTURES', 'RESIDUE SPECIFICATION',# 'MODE',
            'LIGAND ALIGNMENT DATA', 'RECEPTOR ALIGNMENT DATA (OPTIONAL)', 'OUTPUT', 'LOGGING/HELP']

        self.groups = { 'INPUT STRUCTURES' : ['--receptor', '--ligand'],
           
                        'RESIDUE SPECIFICATION' : ['--residue'  ],

                        #'MODE' : ['--covresidue', '--covmap'],

                        'LIGAND ALIGNMENT DATA': ['--ligindices', '--ligsmarts'],


                        'RECEPTOR ALIGNMENT DATA (OPTIONAL)': ['--recindices', '--recsmarts'],

                        'OUTPUT' : [ '--generaterec', '--genfullrec', '--outputfile'],
                        'LOGGING/HELP':['--log', '--verbose', '--help_advanced'],
                        }
        self.usage = '%s --ligand xxxx.xxx --receptor xxxx.xxx [ alignment mode ]'
        self.usage = None
       
    def initDocs(self):
        """ """
        self.description = (
        """Raccoon PrepareCovalent: typical usages

 Create input ligand/receptor structures for covalent docking:
    %(name)s --receptor receptor.pdb --ligand ligand.pdb --residue A:CYS113 --ligsmarts CSCC --
       
 Create input ligand/receptor structures for covalent docking:
    %(name)s --receptor receptor.pdb --ligand ligand.pdb --residue A:CYS113 --ligsmarts CSCC
        """% { 'name': sys.argv[0] } )
        self.epilog = "(C) 2014 Stefano Forli, MGL, TSRI\n Please cite..."

    def initDefaults(self):
        """ """
        self.default_ligPatternIdx = (0,1)
        self.default_recPatternIdx = None #(0,1)
        self.default_genfullrec = False
        self.default_ligpatternidx = (0,1)
        self.default_recpatternidx = None #(0,1)
        self.default_outputfile = None
        self.default_flex = False

    def showAdvancedHelp(self):
        """ """
        advHelp = ('[ADVANCED HELP GOES HERE')
        print advHelp
        sys.exit(0)

    def parseForeignOpts(self):
        """ remove foreign options from the command line"""
        foreign = {'prepare': []}
        # debug
        if '--debug' in sys.argv[1:]:
            print "\n**** Debug mode activated ***\n\n"
            self.debug = True
            sys.argv.remove('--debug')
        # prepare
        for k in foreign.keys():
            remove = []
            for opt in sys.argv[1:]:
                if opt[2:].startswith(k):
                    foreign[k] = opt
                    remove.append(opt)
            for opt in remove:
                sys.argv.remove(opt)

    def parseOpts(self):
        """ parse command-line options"""
        self.parseForeignOpts()

        self.parser = argparse.ArgumentParser(description = self.description,
            usage=self.usage, epilog=self.epilog,
            formatter_class = argparse.RawDescriptionHelpFormatter)
        for g in self.groups_order:
            group = self.parser.add_argument_group(g)
            for i in self.groups[g]:
                    opts = self.opts[i]
                    group.add_argument(i, **opts)
        self.args = self.parser.parse_args()
        if self.args.help_advanced:
            self.showAdvancedHelp()
        self.alignerArgs = {}
        self.verbose =  self.alignerArgs['verbose'] = self.args.verbose
        self.validateOpts()
        if hasattr(self.args, 'outputfile'):
            self._output = self.args.outputfile
        else:
            self._output = None

    def validateOpts(self):
        """ validate all options and set the modes for the actual calculation"""
        self.validateMode()
        self.setLigandMode()
        self.setReceptorMode()

    def setLigandMode(self):
        """ guess the ligand mode basing on the spefication of indices or SMARTS
            --ligindices
            --ligsmarts [ --ligindices ]
        """
        args = self.alignerArgs #= {}
        self._ligmode = 'idx'
        idxError = ('Ligand indices (--ligindices) must be specified as "number,number" (starting from 1)')
        # ligand mode
        if self.args.ligsmarts:
            self._ligmode = 'smarts'
            args['smartsLig'] = self.args.ligsmarts
            if self.args.ligindices:
                try:
                    # smarts indices are 0-based internally
                    args['smartsIdxLig'] = [int(x)-1 for x in  self.args.ligindices.split(',')]
                except:
                    self.showError(idxError)
            else:
                args['smartsIdxLig'] = self.default_ligPatternIdx
        else:
            if not self.args.ligindices:
                msg = ('Ligand indices (--ligindices) are required, or use SMARTS (--ligsmarts)')
                self.showError(msg)
            else:
                try:
                    args['indicesLig'] = [int(x) for x in  self.args.ligindices.split(',')]
                except:
                    self.showError(idxError)
        self.vprint("LIGAND MODE:"+self._ligmode)

    def setReceptorMode(self):
        """ set receptor alignment mode (OPTIONAL)"""
        args = self.alignerArgs #= {}
        self._recmode = None
        idxError = ('Receptor indices (--recindices) must be specified as "number,number" (starting from 1)')
        # receptor mode
        if self.args.recsmarts:
            self._recmode = 'smarts'
            args['smartsRec'] = self.args.recsmarts
            if self.args.recindices:
                try:
                    # smarts indices are 0-based internally
                    args['smartsIdxRec'] = [int(x)-1 for x in  self.args.recindices.split(',')]
                except:
                    self.showError(idxError)
            else:
                args['smartsIdxRec'] = self.default_recPatternIdx
        else:
            if self.args.recindices:
                self._recmode = 'idx'
                try:
                    args['indicesRec'] = [int(x) for x in  self.args.recindices.split(',')]
                except:
                    self.showError(idxError)
        # check that specified residue is managed, then no SMARTS/indices is required
        if (self._recmode == None):
            unsupported = self.checkManagedResidue()
            if len(unsupported):
                msg = ('The specified residue types are not supported automatically: %s. Use either '
                        '--recsmarts or --recindices.' % (','.join(unsupported)))
                self.showError(msg)
        args['genFullRec'] = self.args.genfullrec
         
    def validateMode(self):
        """ set session mode to be:
                - normal (lig + receptor + residue|residuetype)
        """
        self.mode = 'normal'
        if not self.args.residue:
            msg = ('No residue specification. Use at least --residue to specify which residue to modify')
            self.showError(msg)
        if self.args.residue and not self.args.ligand:
            msg = ('option --residue requires --ligand to be specified')
            self.showError(msg)
   
    def checkManagedResidue(self):
        """ check that requested residue type is
            among managed residues: CYS, SER, THR, LYS
        """
        supported = ['cys', 'ser', 'thr', 'lys']
        unsupported = []
        for r in self.args.residue:
            if ':' in r:
                r = r.split(":")[1][:3].lower()
                if not r in supported:
                    unsupported.append(r)
        return set([ x.upper() for x in unsupported])

       



    def parseLigAlignmentOpts(self):
        """ validate parameters for ligand alignment """
        self._ligmode = None
        # SMART INDICES
        self.ligsmarts = self.args.ligsmarts
        if self.ligsmarts:
            self._ligmode = 'smarts'
        self.ligsmartsidx = self.args.ligsmartsidx
        try:
            if not self.ligsmartsidx == None:
                self.ligsmartsidx = [ int(x) for x in self.ligsmartsidx.split(',') ]
        except:
            print "ERROR (--ligsmartsidx): indices must be specified with number pairs, e.g. 1,3"
            sys.exit(1)
        self._ligmode = 'smarts'

        # ATOM INDICES
        self.ligatomidx = self.args.ligatomidx
        if self.ligatomidx and self._ligmode == 'smarts':
            msg = ('Multiple alignment modes selected: use either SMARTS (\'--ligsmarts\')'
                    'or atom indices(\'--ligatomidx\')')
            self.showError(msg)
        # INDICES
        try:    
            if not self.ligatomidx == None:
                self.ligatomidx = [ int(x) for x in self.ligatomidx.split(',') ]
        except:
            print "ERROR (--ligatomidx): indices must be specified with number pairs, e.g. 1,3"
            sys.exit(1)


    def parseRecAlignmentOpts(self):
        """ """
        # validate receptor SMARTS indices
        try:
            self.recsmartsidx = self.args.recsmartsidx
            if not  self.recsmartsidx == None:
                self.recsmartsidx = [ int(x) for x in self.recsmartsidx.split(',') ]
        except:
            print "ERROR (--recsmartsidx): indices must be specified with number pairs, e.g. 1,3"
            sys.exit(1)
        # validate receptor atom indices  
        try:
            self.recatomidx = self.args.recatomidx
            if not self.recatomidx == None:
                self.recatomidx = [ int(x) for x in self.recatomidx.split(',') ]
        except:
            print "ERROR (--recatomidx): indices must be specified with number pairs, e.g. 1,3"
            sys.exit(1)

    def showError(self, msg, showhelp=True):
        """ error reporting facility (modify to STDERR)"""
        tag = '[ ERROR ] '
        msg = msg.replace('\n', ('\n%s' % (" " * len(tag) )  ) )
        print("%s%s" % (tag, msg))
        if showhelp:
            print ""
            self.parser.print_help()
        sys.exit(1)

    def vprint(self, msg):
        """ verbose printer"""
        if not self.verbose:
            return
        print("VERBOSE: "+msg)

    def getMolecule(self, fname, perceiveChains=False):
        """ molecule parser"""
        name, ext = self.getFnameExt(fname)
        self.molParser.SetInFormat(ext)
        mol = ob.OBMol()
        result = False
        try:
            result = self.molParser.ReadFile(mol, fname)
        except:
            msg = ('Fail to read molecule "%s" : %s' % (fname, sys.exc_info()[1]))
            self.showError(msg, showhelp=False)
        if not result:
            msg = ('Fail to read molecule "%s"' % (fname))
            self.showError(msg, showhelp=False)
        if perceiveChains:
            self.chainParser.PerceiveChains(mol)
        return mol, name

    def loadMolecules(self):
        """ import ligand and receptor structures"""
        args = self.alignerArgs
        args['rec'], self.recName = self.getMolecule(self.args.receptor)
        args['lig'], self.ligName = self.getMolecule(self.args.ligand, True)


    def getResidues(self, log=None):
        """ find which residue(s) will be modified:
            - specific residues
            - residue types (surface-accessible residues)
        """
        residue = self.args.residue
        if not residue:
            msg = ('No residues specified. Use either to process!')
            self.showError(msg)
        self.residuePool = residue
        chain = None
        if len(self.residuePool) == 0:
            msg = ('No residues to process!')
            self.showError(msg)
        self.residuePool.sort()
        if log:
            fp = open(log, 'w')
            for r in self.residuePool:
                #fp.write( '%s:%s\n' % (self.recName,r))
                fp.write( '%s\n' % (r))
            fp.close()
        #print "========================================================"
        #print "Total number of residues to process [%d]" % len(self.residuePool)
               
    def initParsers(self):
        """ initialize molecule parser and chain perception"""
        self.molParser = ob.OBConversion()
        self.chainParser = ob.OBChainsParser()

    def start(self):
        """ """
        #  XXX validate here that the indices are within the lenght of matches?
        self.initParsers()
        self.loadMolecules()
        self.getResidues('residues.log')
        self.processResidues()

    def processResidues(self):
        """ perform full covalent protocol"""
        outType = 'pdb'
        outType = 'mol2'

        for count, residue in enumerate(self.residuePool):
            #print "Processing residue %s\t (%d/%d)" % (residue,
            #    count+1, len(self.residuePool)),
            print "Processing residue %s" % residue
            if not self._output == None:
                outname = self._output
                outType = os.path.splitext(self._output)[1][1:]
            else:
                outname = '%s_%s_cov_%s' % (os.path.basename(self.ligName),
                    os.path.basename(self.recName),
                        residue.replace(':','').upper())
                outname = outname + "." + outType
            self.vprint("[start] output filename is: "+ outname)
            print ("[start] output filename is: "+ outname)

            self.alignerArgs['resName'] = residue
            if self.verbose:
                print "==============="
                print "[start] Calling CovalentDockingMaker with settings:"
                for k,v in self.alignerArgs.items():
                    if isinstance(v, ob.OBMol):
                        print "KW[%s]: %d atoms" % (k, v.NumAtoms())
                    else:
                        print "KW[%s]:" % k, v

                print "==============="
            #print "\n==============\n", self.alignerArgs
            if self.debug:
                self.alignerArgs['debug'] = True
            aligner = self.x = CovalentDockingMaker(**self.alignerArgs)
            if not aligner.ready:
                print "\t\t\tSKIPPING?"
                continue

            #print "\t=>", outname + '.pdb'
            # XXX LOCAL FILE HERE
            fname = os.path.basename(outname) # + '.pdb'
            #aligner.writeCovalent(outname, 'pdb')
            aligner.writeCovalent(fname, outType)
            print "Writing output filename: %s" % fname

            #aligner.writeCovalent(outname, 'pdbqt')
            # WORKING PATTERNS
            #           LIG     REC
            #    -------------------------------
            #           CSCC    '[SH]C'
            #                   '[$(SH1);$(SC)]'
            #           SCC     '[$([SH1]);$(SC)]'
   
if __name__ == '__main__':
    x = CovalentDockingMaster()