#-------------------------------------------------------------------------------
# Name:        module1
# Purpose:
#
# Author:      oeffner
#
# Created:     05/04/2011
# Copyright:   (c) oeffner 2011
# Licence:     <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python

import os, re, sys, shutil, glob, traceback
import time
import cStringIO
import random, string
from iotbx import pdb, mtz
import scitbx.math.euler_angles
from phaser import AlterOriginSymmate
from phaser import *
import phaser.test_phaser_parsers
from phenix.command_line import get_cc_mtz_mtz
from phenix.command_line import get_cc_mtz_pdb
from phaser.utils2 import *



def ShiftMRsols(solfname, outfname):
    def mymatch(m):
        ( xs, ys, zs ) = m.group( 1, 2, 3)
        return "FRAC %.5f %.5f %.5f" % (float(xs)%1, float(ys)%1, float(zs)%1)

    solsstr = open(solfname,"r").read()
    newsolstr = re.sub(r"FRAC\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)",
     mymatch, solsstr)
    mfile = open(outfname,"w")
    mfile.write(newsolstr)
    mfile.close()



def RealizeSoluPDB(xtalsym, pdbfname, euler, frac, outfname):
    mstr = "\nRealizing MR solution with model %s at Euler: %s and Frac: %s\n" \
      %(pdbfname, str(euler), str(frac))
    print mstr
    pdbmrsol = pdb.input(file_name= pdbfname)
    rot = scitbx.math.euler_angles.zyz_matrix(float(euler[0]), float(euler[1]), float(euler[2]))
    trans = (float(frac[0]), float(frac[1]), float(frac[2]))
    mrsolhiearchy = pdbmrsol.construct_hierarchy().deep_copy()
    modelatoms = mrsolhiearchy.atoms()
    # do the rotation in cartesian space
    rotatomsites = rot * modelatoms.extract_xyz()
    ucell = xtalsym.unit_cell()
    # get into fractional space
    fracatomsites = ucell.fractionalize(sites_cart= rotatomsites )
    # position and rotate pdb to the mr solution
    mrsolsites = fracatomsites + trans
    # cast back to cartesian coordinates
    modelatoms.set_xyz(ucell.orthogonalize(sites_frac= mrsolsites))
    lastpdbfname = outfname
    mrsolpdbfile = open(lastpdbfname, "w")
    mrsolpdbfile.write(mrsolhiearchy.as_pdb_string())
    mrsolpdbfile.close()
    return mstr




def getCC(mtz1, fname2, temp):
    outstr = "Calculating CC between %s and %s\n" %(mtz1, fname2)
    print outstr
    dmpstr = cStringIO.StringIO()
    myrand = random.SystemRandom()
    randstr = "".join(myrand.sample(string.letters+string.digits, 5))
    randtemp = temp + "_" + randstr # ensure uniqueness of tempfolder to avoid race conditions
    ret = (None, None)

    try:
        if fname2[-4:] == ".mtz":
            cmdargs = ["mtz_1="+ mtz1, "mtz_2=" + fname2,
             "temp_dir=" + randtemp, "output_dir=" + randtemp]
            myobj = get_cc_mtz_mtz.get_cc_mtz_mtz(cmdargs, quiet=True, out=dmpstr)
            outstr += dmpstr.getvalue()
            ret = (Roundoff(float(myobj.final_cc),4), 0.0)
        else:
            cmdargs = ["mtz_in="+ mtz1, "pdb_in=" + fname2,
              "temp_dir=" + randtemp, "output_dir=" + randtemp]
            myobj = get_cc_mtz_pdb.get_cc_mtz_pdb(cmdargs, quiet=False, out=dmpstr)
            outstr += dmpstr.getvalue()
            ret = (Roundoff(float(myobj.found_overall),4), Roundoff(float(myobj.found_region),4))

        shutil.rmtree(randtemp)
    except Exception:
        outstr += traceback.format_exc()

    return outstr, ret




def CreateZeroRNPcycleMap(datrun, rano, rootname, pdbfname, targseqnames,
  modelrmsd, overwrite=False, llgref = None, debug = False):
    mtzmapname = rootname + ".1.mtz"
    outstr = "Making map %s from 0 cycle RNP of %s\n" %(mtzmapname, pdbfname)
    if os.path.exists(mtzmapname) and not overwrite:
        outstr += mtzmapname + " exists already.\n"
        return outstr

    rnptargetinput = InputMR_RNP()
    rnptargetinput.setSPAC_HALL(datrun.getSpaceGroupHall())
    rnptargetinput.setCELL6(datrun.getUnitCell())
    #rnptargetinput.setREFL(rano.getMiller(),rano.getCorrectedF(),rano.getCorrectedSIGF())
    rnptargetinput.setREFL_DATA(datrun.getDATA())
    rnptargetinput.setNORM(rano.getBINS(), rano.getANISO(), rano.getSOLK(), rano.getSOLB())
    rnptargetinput.setROOT(rootname)
    rnptargetinput.setJOBS(1)
    rnptargetinput.setXYZO(False) # only want the map
    rnptargetinput.addENSE_PDB_RMS("RNP0cycle", pdbfname, modelrmsd)
    rnptargetinput.addSOLU_ORIG_ENSE("RNP0cycle")
    for i in range(len(targseqnames)):
        rnptargetinput.addCOMP_PROT_SEQ_NUM(targseqnames[i], 1)

    if debug:
        rnptargetinput.setMUTE(False)
    else:
        rnptargetinput.setMUTE(True)
# 0 macrocycles, no refinement
    rnptargetinput.setMACM_PROT("OFF")
    rnptarget = runMR_RNP(rnptargetinput)
    if not rnptarget.Success():
        outstr += rnptarget.ErrorName() + " TargetPDB RNP ERROR: " + rnptarget.ErrorMessage() + "\n"
        #raise Exception, m
    else:
        if llgref != None:
            llgref.append(rnptarget.getValues()[0] )
    open(rootname + ".log","w").write(rnptarget.logfile())
    return outstr



def Sol2EulerFrac(solstr):
    enseulerfrac = []
    MRsolens = re.findall("SOLU 6DIM ENSE\s+(.+)\s+(EULER.+)\s+(FRAC.+)\s+BFAC.+", solstr, re.MULTILINE)
    for mrsolens in MRsolens:
        euler = re.findall(r"EULER\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)",
          mrsolens[1], re.MULTILINE)
        frac = re.findall(r"FRAC\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)",
          mrsolens[2], re.MULTILINE)
        enseulerfrac.append((mrsolens[0], euler[0], frac[0]))
    return enseulerfrac



def CalcCCFromMRsolution(targetpdbid, MRfolder, calclabel, workdir, partmodels,
  bCalcCC = False, mylock=None, nthreads=1, debug = False):
    stroutput = time.asctime( time.localtime(time.time())) + "\n\n"
    calcdir = os.getcwd()
    MRfolder = PosixSlashFilePath(MRfolder)
    MRfolder = RecaseFilename2(MRfolder)
    mrinput = open(os.path.join(MRfolder, calclabel + "MRinput.txt"), "r").read()
    mtzfname = re.findall(r"^HKLIN\s+(\S+)", mrinput, re.MULTILINE)[0] # to get rid of any \r character
    # if filepath generated on windows but now running on unix then change \ to /
    mtzfname = mtzfname.replace("\\", os.path.sep)

    mtzfcol = re.findall(r"^LABIN\s*F\s*=\s*(\S*)\s*SIGF\s*=\s*(\S+)\s*", mrinput, re.MULTILINE)
    mtzicol = re.findall(r"^LABIN.*I\s*=\s*(\S*)\s*SIGI\s*=\s*(\S*)\s*", mrinput, re.MULTILINE)
    myinput = InputMR_DAT()
    mtzfname = PosixSlashFilePath(mtzfname)
    mtzfname= RecaseFilename2(mtzfname)
    mtzobj = mtz.object(mtzfname)
    xtalsym = mtzobj.crystals()[0].crystal_symmetry()

    myinput.setHKLI(mtzfname)
    if len(mtzicol) > 0:
# use intensity columns if present
        stroutput += "Using columns: " + str(mtzicol[0])
        myinput.setLABI_I(mtzicol[0][0])
        myinput.setLABI_SIGI(mtzicol[0][1])
    else:
        myinput.setLABI_F(mtzfcol[0][0])
        myinput.setLABI_SIGF(mtzfcol[0][1])
        stroutput += "Using columns: " + str(mtzfcol[0])

    if debug:
        myinput.setMUTE(False)
    else:
        myinput.setMUTE(True)
    datrun = runMR_DAT(myinput)

    if datrun.Success():
        solfname = calclabel + "{" + workdir + "}" + partmodels + ".sol"
        fname = os.path.join(MRfolder, solfname)
        mstr = "Calculating LLGs from %s\n" %solfname
        stroutput += mstr
        print mstr
        fname = PosixSlashFilePath(fname)
        fname= RecaseFilename2(fname)
        if os.path.exists(fname) == False:
            mstr = "Error: %s doesn't exist.\n" %fname
            stroutput += mstr
            print mstr
            return stroutput
# filename for mrsol file where fractional coordinates have been cast into unitcell near origo
        newsolfname = os.path.join(MRfolder, calclabel + "_shift{" + workdir + "}" + partmodels + ".sol")
        ShiftMRsols(fname, newsolfname)

        #mrparser = phaser.test_phaser_parsers.SolutionFileParser(file(fname))
        mrparser = phaser.test_phaser_parsers.SolutionFileParser(file(newsolfname))
        mrsols = mrparser.get_dot_sol()
# get ensemble name and sequence file name from some phaser input file
        MRpdbens = re.findall(r"^ENSEMBLE\s+(\S+)\s+PDB\s+(\S+)\s+IDENTITY\s+(\S+)\s+",
          mrinput, re.MULTILINE)
        seqfnames = re.findall(r"^COMPOSITION PROTEIN SEQUENCE\s+(\S+)\s+NUM.+",
          mrinput, re.MULTILINE)
# because strings are immutable we iterate with int index
# rather than list elements in order to change the strings
        for i in range(len(seqfnames)):
            seqfnames[i] = PosixSlashFilePath(seqfnames[i])
            seqfnames[i] = RecaseFilename2(seqfnames[i])

        myano = InputANO()
        myano.setSPAC_HALL(datrun.getSpaceGroupHall())
        myano.setCELL6(datrun.getUnitCell())
        #myano.setREFL(datrun.getMiller(),datrun.getFobs(),datrun.getSigFobs())
        myano.setREFL_DATA(datrun.getDATA())
        if debug:
            myano.setMUTE(False)
        else:
            myano.setMUTE(True)
        rano = phaser.runANO( myano )
        if not rano.Success():
            mstr = rano.ErrorName() + "ANO ERROR: " + rano.ErrorMessage() + "\n"
            stroutput += mstr
            return stroutput

# make map of targetpdb with its own mtz file by doing RNP with 0 cycles with an origin solution
        rnptargetrootname = targetpdbid + "_map"
        tmpname = targetpdbid + ".pdb"
        tmpname2 = PosixSlashFilePath(tmpname)
        targetpdb = RecaseFilename2(tmpname2)
        if os.path.exists(rnptargetrootname + ".1.mtz") == False:
            stroutput += CreateZeroRNPcycleMap(datrun, rano, rnptargetrootname,
              targetpdb, seqfnames, 0.1, debug = debug)

# setup RNP refinement of the MR model
        myinput = InputMR_RNP()
        myinput.setSPAC_HALL(datrun.getSpaceGroupHall())
        myinput.setCELL6(datrun.getUnitCell())
        #myinput.setREFL(rano.getMiller(),rano.getCorrectedF(),rano.getCorrectedSIGF())
        myinput.setREFL_DATA(datrun.getDATA())
        myinput.setNORM(rano.getBINS(), rano.getANISO(), rano.getSOLK(), rano.getSOLB())

        rootname = os.path.join(MRfolder, calclabel + "_VRMS{" + workdir + "}" + partmodels)
        myinput.setROOT(rootname)
        myinput.setHKLO(False)
        myinput.setXYZO(False)
# refine all solutions in one go
        myinput.setSOLU(mrsols)
        #if bfullres :
        #    myinput.setRESO_HIGH(0.0) # use all reflections
        myinput.setJOBS(nthreads)
        for ensemble in MRpdbens:
            modelfname = ensemble[1].replace("\\", os.path.sep).strip() # to get rid of any \r character
            modelfname = PosixSlashFilePath(modelfname)
            modelfname = RecaseFilename2(modelfname)
            myinput.addENSE_PDB_ID(ensemble[0], modelfname, float(ensemble[2]))

        for seqfname in seqfnames:
            myinput.addCOMP_PROT_SEQ_NUM(seqfname, 1)

        if debug:
            myinput.setMUTE(False)
        else:
            myinput.setMUTE(True)
# RNP vrms refinement, no rotation or translation refinement
        #myinput.setMACM_PROT("OFF")
        myinput.setMACM_PROT("CUSTOM")
        myinput.addMACM(False, False, False, True )
        #myinput.setVERB(True)
        #myinput.setDEBU(True)
        #del(run)

        run = runMR_RNP(myinput)
        if run.Success():
#save result to log and sol files
            logfname = rootname + ".log"
            logfile = open(logfname,"w")
            logfile.write(run.logfile())
            logfile.close()
            rnpsolfname = rootname + ".sol"
            solfile = open(rnpsolfname,"w")
            solfile.write(run.getDotSol().unparse())
            solfile.close()
            fname = os.path.join(MRfolder, calclabel +"LLG{" + workdir + "}" + partmodels + ".txt")
            llgfile = open(fname,"w")
            llgs =[]
            for (i,llg) in enumerate(run.getValues()):
                llgs.append(llg)
                mstr = "llg= %s\n" %llg
                #print stroutput
                stroutput += mstr
                print mstr

            llgfile.write(str(llgs))
            llgfile.close()
        else:
            msg = run.ErrorName() + "Model RNP ERROR: " + run.ErrorMessage() + "\n"
            print msg
            return msg
# Set up Correlation Coefficient calculation
        rnpsols = run.getDotSol()
        myCCs = []
        if not bCalcCC:
            return

        mstr += "Calculating correlation coefficients and MLADs of %s with %d solutions\n" \
          %(rnpsolfname, len(rnpsols))
        stroutput += mstr
        print mstr
# set up a temp folder for getCC
        CCtemp = ""
        tmpdir = targetpdbid + "_" + workdir + "_"+ partmodels
        if "TMP" in os.environ:
            CCtemp = os.path.join(os.environ["TMP"], tmpdir)
        if "TEMP" in os.environ:
            CCtemp = os.path.join(os.environ["TEMP"], tmpdir)
        #if os.path.exists("/scratch"):
        #    CCtemp = os.path.join("/scratch", tmpdir) # telemachus cluster solution
        if CCtemp == "":
            CCtemp = os.path.join(os.environ["HOME"], "CCTEMP", tmpdir) #last resort

        nsols = len(rnpsols)
        for (j, rnpsol) in enumerate(rnpsols):
            i = j+1
            mstr = "\nRNP solution %d out of %d\n" %(i, nsols)
            stroutput += mstr
            print mstr
            root = os.path.join(MRfolder, calclabel + "CalcCC_%d" %i)
            mtzfname = root + ".1.mtz"

            myinput = InputMR_RNP()
# RNP vrms refinement, no rotation, no translation, no Bfactor refinement
            #myinput.setMACM_PROT("OFF")
            myinput.setMACM_PROT("CUSTOM")
            myinput.addMACM(False, False, False, True )
            myinput.setSPAC_HALL(datrun.getSpaceGroupHall())
            myinput.setCELL6(datrun.getUnitCell())
            #myinput.setREFL(rano.getMiller(),rano.getCorrectedF(),rano.getCorrectedSIGF())
            myinput.setREFL_DATA(datrun.getDATA())
            myinput.setNORM(rano.getBINS(), rano.getANISO(), rano.getSOLK(), rano.getSOLB())

            myinput.setROOT(root)
            myinput.setXYZO(False)
            myinput.setHKLO(True)
            myinput.setJOBS(1)
            if debug:
                myinput.setMUTE(False)
            else:
                myinput.setMUTE(True)
            solu = phaser.mr_solution([rnpsol]) # convert this rnpsol into a string
            myinput.setSOLU(solu)
            for seqfname in seqfnames:
                myinput.addCOMP_PROT_SEQ_NUM(seqfname, 1)

            lastmodelname = ""
            lastsolcomponent = ""
            lastvrms = 0.0

            ensvrmsdic = {} # make a dictionary of vrms values with ensemblename as keys
            for vitem in rnpsol.NEWVRMS.items():
                ensvrmsdic[vitem[0]] = vitem[1][0]

            ensset = set([])
            CCcur = []
            solstr = rnpsol.unparse()
            EnseulerFracs = Sol2EulerFrac(solstr)
            for enseulerfrac in EnseulerFracs:
                print "Ensemble: %s" %enseulerfrac[0]
                vrms = ensvrmsdic[enseulerfrac[0]]
# find filename of the pdb file with this ensemble name from the MR input file
                for ens in MRpdbens:
                    if ens[0] == enseulerfrac[0]:
                        modelfname = ens[1].replace("\\", os.path.sep).strip() # to get rid of any \r character
                        modelfname = PosixSlashFilePath(modelfname)
                        modelfname= RecaseFilename2(modelfname)
# prepare for refinement by adding all components
                        if enseulerfrac[0] not in ensset: # don't over-define the same ensemble
                            myinput.addENSE_PDB_RMS(ens[0], modelfname, vrms)
                        ensset.add(enseulerfrac[0])
                        #break
# create map of this component from RNP 0 cycles with phaser
                        lastroot = os.path.join(MRfolder, ens[0] + "." + str(i))
                        curpdbfname = lastroot +".pdb"
                        stroutput += RealizeSoluPDB(xtalsym, modelfname, enseulerfrac[1],
                         enseulerfrac[2], curpdbfname)
                        llgs = []
                        stroutput += CreateZeroRNPcycleMap(datrun, rano, lastroot,
                          curpdbfname, seqfnames, vrms, overwrite=True, llgref=llgs,
                          debug = debug)
# get CC with that map and target map
                        lastrootmtz = lastroot+ ".1.mtz"
                        outstr, cc = getCC(rnptargetrootname + ".1.mtz", lastrootmtz, CCtemp)
                        stroutput += outstr
                        stroutput += "CCmap of %s component:  %s\n" %(ens[0], str(cc[0]))
# calculate CC between target and pdbfile of last component
                        outstr, ccpdb = getCC(rnptargetrootname + ".1.mtz", curpdbfname, CCtemp)
                        stroutput += outstr
                        stroutput += "CCpdb of %s component: %s\n" %(ens[0], str(ccpdb[0]))
                        AltOrigres, output_file, outstr = AlterOriginSymmate.GetAltOrigSymmates(
                          pdbfixedrefname = curpdbfname, movingname = targetpdb, pxtalsym = [xtalsym])
                        CCcur.append((cc[0],ccpdb[0], Roundoff(AltOrigres[0][0],4), Roundoff(llgs[0],3)))
                        stroutput += outstr
                        stroutput += "\n"
                        if not debug:
                            files = glob.glob(lastroot + "*")
                            for f in files:
                                try:
                                    os.remove(f)
                                except Exception, m:
                                    stroutput += traceback.format_exc()

            llgrun = runMR_RNP(myinput)
# get CC of map of all components against target map
            if llgrun.Success():
                open(root + ".log","w").write(llgrun.logfile())
                outstr, CCtotal = getCC(rnptargetrootname + ".1.mtz", mtzfname, CCtemp)
                stroutput += outstr
                stroutput += "CC of all components: %s\n\n" %CCtotal[0]
# append all CCs and MLAD to list
                llgall = Roundoff(llgrun.getValues()[0], 3)
                CCall =[CCtotal[0], CCtotal[1], llgall]
                CCall.extend(CCcur)
                myCCs.append(CCall)
# Format of correlation coefficient text file is like
# [[0.3576, 0.0, 89.639, (0.3105, 0.203, 0.8105, 45.367), (0.2222, 0.136, 1.0067, 23.764)],
# [0.283, 0.0, 37.207, (0.3056, 0.203, 0.8174, 12.976), (0.0381, 0.011, 3.4518, 62.891)]]
# Each list element are the CCs and LLG values corrresponding to a solution
# First 2 numbers in a list are the overall and local CCs and LLG value for the solution
# including all ensembles
# Subsequent tuples are CCmap, CCpdb, MLAD and llg value of each individual ensemble in that solution
            else:
                mstr = llgrun.ErrorName() + " MR_LLG ERROR: " + llgrun.ErrorMessage() + "\n"
                #raise Exception, m
                stroutput += mstr
                print mstr
                #return stroutput
            ccroot = os.path.join(MRfolder, calclabel + "CalcCC_")
            fname = ccroot + "{" + workdir + "}" + partmodels + ".txt"
            CCfile = open(fname,"w")
            CCfile.write(str(myCCs))
            CCfile.close()

            if not debug:
                files = glob.glob(root + "*.mtz")
                for f in files:
                    try:
                        os.remove(f)
                    except Exception, m:
                        stroutput += traceback.format_exc()

            if os.path.exists(CCtemp): # clean up just in case a previous calculation crashed
                try:
                    shutil.rmtree(CCtemp)
                except Exception:
                    stroutput += traceback.format_exc()
    else:
        mstr = datrun.ErrorName() + " MR_DAT ERROR: " + datrun.ErrorMessage() + "\n"
        #raise Exception, m
        stroutput += mstr

    if not debug:
        files = glob.glob("MinMLAD*")
        for f in files:
            try:
                os.remove(f)
            except Exception:
                stroutput += traceback.format_exc()

    return stroutput + "\n"




def CalcCCFromSCOP_MRcalcs(selectedpdbfname, scoppdbsfname, calcdir):
    mfile = open(selectedpdbfname,"r")
    selectedpdbs = eval(mfile.read())
    mfile.close()

    mfile = open(scoppdbsfname,"r")
    scoppdbs = eval(mfile.read())
    mfile.close()

    logfname = os.path.join(calcdir, os.path.basename(selectedpdbfname) + ".log")

    n = 0
    l = len(selectedpdbs)
    for (n,pdbs) in enumerate(selectedpdbs):
        mstr = ""
        for scoppdb in scoppdbs:
            if scoppdb[0].lower().find(pdbs[0].lower()) > -1:
                scopid = scoppdb[1]
# only want directories named with SCOP fold id
                if scopid == "UnSCOPed":
                    foldID = "UnSCOPed"
                else:
                    foldID = scopid.split(".")[0] + "." + scopid.split(".")[1]

                #foldID = scopid.split(".")[0] + "." + scopid.split(".")[1]
                fold_targetpdb = os.path.join(calcdir,foldID,pdbs[0])
                if not os.path.exists(fold_targetpdb):
                    fold_targetpdb = os.path.join(calcdir,foldID,pdbs[0].lower())

                BLASTtxtchainfname =  os.path.join(fold_targetpdb, "Strider" + pdbs[0] + "_BLASTchain_ids.txt")
                psearches = eval(open(BLASTtxtchainfname,"r").read())

                mstr += "CalcCCFromMRcalcs %s, %s, %d. out of %d\n" %(fold_targetpdb, pdbs[0], n, l)

                for psearch in psearches:
                    workdir = psearch[1][-1] # the last element in ['2WBJ_C0', '2WBJ_C0,1NFD_B0']
                    for partmodels in psearch[1]:

                        MRfolder = os.path.join(fold_targetpdb,"MR", workdir)

                        try:
                            os.chdir(fold_targetpdb)
                            mstr += "Calculating CCs from %s\n" %solfname
                            mstr += CalcCCFromMRsolution(pdbs[0], MRfolder, "Strider", workdir,
                              partmodels, bCalcCC = True)

                        except Exception, m:
                            mstr += "Error in %s, %s, %s\n" %(selectedpdbfname, scoppdbsfname, m)
                            #raise

                        finally:
                            os.chdir(calcdir)

        logfile = open(logfname,"a")
        logfile.write(mstr)
        logfile.close()

    mstr = "Done with " + selectedpdbfname + "\n"
    logfile = open(logfname,"a")
    logfile.write(mstr)
    logfile.close()




def CalcRNP_FromMRcalcs(pdbid, calclabel):
    mstr = ""

    BLASTtxtchainfname = os.path.join(pdbid, calclabel + pdbid + "_BLASTchain_ids.txt")
    if os.path.exists(BLASTtxtchainfname):
        psearches = eval(open(BLASTtxtchainfname,"r").read())

        mstr += "CalcCCFromMRcalcs %s\n" %pdbid
        os.chdir(pdbid)
        for psearch in psearches:
            workdir = psearch[1][-1] # the last element in ['2WBJ_C0', '2WBJ_C0,1NFD_B0']
            partmodels = workdir # on this occassion we're only interested in multicomponent refinements

            MRfolder = os.path.join("MR", workdir)
            templdir = os.path.join("TemplateSolutions", workdir)

            try:
                #print "MRfolder= %s, partmodels = %s" %(MRfolder, partmodels)
                mstr += CalcCCFromMRsolution(pdbid, MRfolder, calclabel, workdir,
                  partmodels, bCalcCC = True, nthreads=6)

                mstr += CalcCCFromMRsolution(pdbid,
                    templdir, calclabel + "RNP", workdir, ".superposed",
                    bCalcCC = True)

            except Exception, m:
                mstr += "Error in %s, %s\n" %(pdbid, m)
                print "Error in %s, %s\n" %(pdbid, m)
                #raise

    else:
        mstr += "Error: %s doesn't exist" %BLASTtxtchainfname
        print "Error: %s doesn't exist" %BLASTtxtchainfname

    logfname = os.path.basename(pdbid) + "CalcCC.log"
    logfile = open(logfname,"a")
    logfile.write(mstr)
    mstr = "Done with " + pdbid + "\n"
    logfile.write(mstr)
    logfile.close()




if __name__=="__main__":
    if len(sys.argv) > 3:
        CalcCCFromSCOP_MRcalcs(sys.argv[1], sys.argv[2], sys.argv[3])
    else:
        CalcRNP_FromMRcalcs(sys.argv[1], sys.argv[2])

    print "Done for today"