#-------------------------------------------------------------------------------
# Name:        AlterOriginSymmate
# Purpose:
#
# Author:      Robert Oeffner
#
# Created:     26/12/2010
# Copyright:   (c) oeffner 2010
# Licence:     <your licence>
#-------------------------------------------------------------------------------
#!/usr/bin/env python



from copy import deepcopy
from fractions import Fraction
import itertools
from boost import rational
from types import *
import sys, os, math, re
import mmtbx.alignment
from iotbx import pdb
from iotbx import mtz
from cctbx import sgtbx
from scitbx.array_family import flex
import ccp4io_adaptbx, cctbx
import phaser.test_phaser_parsers
from phaser import ensembler
from phaser import tbx_utils
from iotbx.pdb import amino_acid_codes
import scitbx.math.euler_angles
import logging
import pickle
from phaser import phenix_interface
from libtbx.phil import parse
import getopt
from scitbx import golden_section_search
import cStringIO
from libtbx.utils import Sorry
from phaser.utils2 import *
from mmtbx.command_line import sort_hetatms
import random, string

chainid = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
debugfolder = "AltOrigSymMatesFiles"




helpstr = """
For different molecular replacement solutions from the same dataset %s
finds for each chain in the moving solution the copy closest to chains in the fixed solution with
respect to all symmetry operations and alternative origin shifts permitted by the
spacegroup of the crystal.

Usage:
    %s moving_pdb=pdbfile1 fixed_pdb=pdbfile2

or more generally

    %s input.txt

    where the content of input.txt file is in PHIL format:

%s


The PHIL input specifies the MR solution files either as \"moving_pdb\" and
\"fixed_pdb\" or as the scopes, \"moving\" and \"fixed\". Both of the scopes
\"moving\" and \"fixed\" must be valid. These two scopes hold the parameter
\"xyzfname\" and the scopes \"mrsolution\" and \"pickle_solution\". An MR solution
is specified by following only one of the procedures a, b and c:

    (a) Assign the "moving_pdb" or "fixed_pdb" parameter to the pdb file from one
    of the MR solutions in question. This is useful for simple cases such as when
    the solution comes from the PHENIX AutoMR GUI is or when it is from a different
    MR program other than Phaser.

    (b) Specify a Phaser MR solution by assigning the file name of the solution
    file to the "moving.mrsolution.solfname" or "fixed.mrsolution.solfname"
    parameter as well as assigning the IDs of the ensembles and the corresponding
    MR models to the parameters "moving.mrsolution.ensembles.name" or
    "fixed.mrsolution.ensembles.name" and "moving.mrsolution.ensembles.xyzfname"
    or "fixed.mrsolution.ensembles.xyzfname" respectively. Multiple components
    are specified as multiple ensembles. This way of specifying solutions is useful
    for solution files produced when Phaser is run from the command line.

    (c) Use a solution from the PhaserMR GUI in PHENIX. In that case the parameter
    "moving.pickle_solution.pklfname" or "fixed.pickle_solution.pklfname" is assigned
    to the solution file that is produced by PHENIX after an MR calculation. The
    "moving.pickle_solution.philfname" or "fixed.pickle_solution.philfname" is then
    assigned to the input file for the MR calculation. This way of specifying
    solutions is useful when MR solutions are available from having run Phaser
    through the PHENIX interface.

There is no restriction on what method to use for the fixed scope depending on
the method used on the moving scope, and vice versa. If space group and the unit
cell dimensions is not available when for instance (b) is used for both the moving
and the fixed scopes then these need to be specified by assigning the parameter,
\"spacegroupfname\", to a PDB file with a CRYST1 record or to an MTZ file with that
information. Typically this would be the data file used for the molecular replacement
calculation.

Output:
    The closest match between chains in the moving section to chains in the fixed
    section is saved with the name of the reference_pdb,
    concatenated with the name of the moving_pdb but prepended with \"MinMLAD_\".
    A log file with the name of the reference_pdb, concatenated with the name of the
    moving_pdb but prepended with \"AltOrigSymMLAD_\" contains all screen output.
    All files are saved in the folder containing reference_pdb. A good match
    between two chains usually have a MLAD value less than 1.5
    whereas a bad match usually have a value larger than 2. This is a rule of thumb
    and exceptions do occur.

If a Phaser MR solution or a pickle solution is supplied, %s
iterates through each solution and does the respective translation and rotation of
model_pdb as to realize each MR solution as a temporary pdb file. These are in turn
used as moving or as fixed pdb files.

For poor homologues or MR solutions try using \"use_all_SSM=True\".

Algorithm:
    %s computes configurations by looping over all
    symmetry operations and alternative origin shifts. These are computed using
    procedures outlined in R. W. Grosse-Kunstleve, Acta Cryst. 1999, A55, 383-395.
    An alignment between C-alpha atoms from moving_pdb and reference_pdb is
    computed using SSM (E. Krissinel* and K. Henrick, Acta Cryst. 2004. D60, 2256-2268).
    If SSM fails or if the MLAD score achieved with the SSM alignment  is larger
    than 2.0 then a sequence alignment algorithm in the CCTBX is used instead.
    To estimate the best match a distance measure between the aligned atoms is
    computed for each configuration. The mean log absolute deviation (MLAD) is
    defined as:

      MLAD(dR) = Sum( log(dr*dr/(|dr| + 0.9) + max(0.9, min(dr*dr,1)))),

    where dr is the difference vector between a pair of aligned C-alpha atoms and the sum
    is taken over all atom pairs in the alignment. MLAD is used as the distance measure.
    MLAD will downplay the contributions of spatially distant C-alpha atom pairs which
    may otherwise lead to incorrect matches when the space group has a floating origin
    and the chains tested against one another consists of large unaligned domains.

    %s will for each chain in reference_pdb find the
    smallest MLAD with a copy of each chain in moving_pdb for a given symmetry
    operation and alternative origin. When all chains in reference_pdb have been
    tested these copies will be saved to a file. A warning is issued if not all
    saved copies of chains in moving_pdb have the same origin.

    For spacegroups with floating origin the minimum MLAD is found
    by doing a Golden Sectioning minimization along the polar axis for each copy
    of chains in moving_pdb.

Debug mode:
    Invoking the debug flag produces individual pdb files of the C-alpha atoms used for
    each SSM alignment of moving_pdb for each permitted symmetry operation and alternative
    origin. A gold atom is placed at the centroid of the C-alpha atoms. Similar files are
    produced for the reference_pdb. These files are stored in a subfolder named
    "AltOrigSymMatesFiles".

    If a floating origin is present in the space group a table of MLAD values is produced by
    sliding a copy of moving_pdb along the polar axis in the fractional interval [-0.5, 0.5]
    for each permitted symmetry operation and alternative origin.

    """

def usage(progname, mphilstr):
    print helpstr %(progname, progname, progname, mphilstr, progname, progname, progname)


def RunFromCommandline (parser, args=sys.argv):
    philargs = parser.parse_args(args).phils
    accum = tbx_utils.PhilAccumulator( master_phil = parser.master_phil )
    output_file_names = []
    logstr =""

    if len(philargs) > 0:
        for argument in philargs:
            argument.process( accumulator = accum )

        merged_phil = accum.merge()
        params = merged_phil.extract().AltOrigSymMates

        logstr = merged_phil.as_str()
        verbose = params.verbose
        if verbose:
            print logstr

        usebestSSM = True
        if params.use_all_SSM:
            usebestSSM = False
        chosenorigin = params.chosenorigin
        usesymmop = True
        if params.no_symmetry_operations:
            usesymmop = False
        movehetatms = params.movehetatms
        debug = params.debug
        mediocre_mlad = 2.0

        if debug and os.path.exists(debugfolder) == False:
            os.mkdir(debugfolder)

# get crystal dimensions either from a pdb file or from an mtz file if it is present
        xtaldimensionsources = [params.fixed_pdb, params.moving_pdb,
                                   params.spacegroupfname]
        pxtalsym = [None]
        xtalsym = None
# first see if fixed and moving scope explicitly provide pdb files with spg info
        for i,e in enumerate(xtaldimensionsources):
            if e != None and i < 2:
                xtal = pdb.input(file_name = e)
                pxtalsym[0] = xtal.crystal_symmetry()
                if pxtalsym[0] != None:
                    break
# otherwise see if spacegroupfname has the spg info
            if e!= None and i==2:
                try:
                    mtzobj = mtz.object(e)
                    xtalsym = mtzobj.crystals()[0].crystal_symmetry()
                except Exception, m:
                    try:
                        xtal = pdb.input(file_name = e)
                        xtalsym = xtal.crystal_symmetry()
                        assert(type(xtalsym) != NoneType)
                    except Exception, m:
                        outstr = "%s is not an MTZ file or a PDB file with a CRYST1 record" %e
                        logstr += outstr + "\n"
                        raise Sorry, m

                pxtalsym[0] = xtalsym

        (fixsoltxts, pxtalsym, fixensmblspdbs, fixlogfname) = GetSolutionsAsPdbs(params.fixed.pickle_solution.philfname,
          params.fixed.pickle_solution.pklfname, params.fixed.mrsolution.solfname,
          params.fixed.mrsolution.ensembles, params.fixed_pdb, pxtalsym)

        (movsoltxts, pxtalsym, movensmblspdbs, movlogfname) = GetSolutionsAsPdbs(params.moving.pickle_solution.philfname,
          params.moving.pickle_solution.pklfname, params.moving.mrsolution.solfname,
          params.moving.mrsolution.ensembles, params.moving_pdb, pxtalsym)

        fixpdbs = [params.fixed_pdb]
        if len(fixsoltxts) > 0:
            fixpdbs = RealizePhaserSolutions(fixsoltxts, pxtalsym[0], fixensmblspdbs, "fix_")

        movpdbs = [params.moving_pdb]
        if len(movsoltxts) > 0:
            movpdbs = RealizePhaserSolutions(movsoltxts, pxtalsym[0], movensmblspdbs, "mov_")

        enumfixpdbs = list(enumerate(fixpdbs))
        enummovpdbs = list(enumerate(movpdbs))
        minAllAltOrigPdbFnames = None
        minoutput_file = None
        minsolmlad = 9e99
        minbestsols = [0,0]
        mini = 0
        minj = 0
        minmovchainid = ""
        minfixchainid = ""
        logstr += "\n"

        for (i,fixpdb) in enumfixpdbs: # iterating over solutions
            for (j,movpdb) in enummovpdbs: # iterating over solutions
                outstr = "\nTesting %d. fixed solution against %d. moving solution" %(i, j)
                outstr += "\nfixed_pdb = " + str(fixpdb) + "\nmoving_pdb = " + str(movpdb)
                logstr += outstr + "\n"

                AllAltOrigPdbFnames, output_file, outstr = GetAltOrigSymmates(movpdb, fixpdb,
                  pxtalsym, chosenorigin, usesymmop, debug, mediocre_mlad = mediocre_mlad,
                  usebestSSM = usebestSSM, verbose = verbose, movehetatms = movehetatms)
                logstr += outstr
                output_file_names.append(output_file)
                for (mlad, movchainid, fixchainid, dummy1, dummy2, dummy3) in AllAltOrigPdbFnames:
                    #outstr =  "\nMLAD of %f for chain %s in %s with chain %s in %s is achieved with solution %d" \
                    #  "\nsaved as %s" %(mlad, movchainid, movpdb, fixchainid, fixpdb, j, output_file)
                    #logstr += outstr
                    if minsolmlad > mlad:
                        minsolmlad = mlad
                        mini = i
                        minj = j
                        minoutput_file = output_file
                        minmovchainid = movchainid
                        minfixchainid = fixchainid
                        minAllAltOrigPdbFnames = AllAltOrigPdbFnames

        logstr += "\n"

        if len(enumfixpdbs) > 1 or len(enummovpdbs) > 1 :
            print " "
            for (mlad, movchainid, fixchainid, dummy1, dummy2, dummy3) in minAllAltOrigPdbFnames:
                outstr = "\nFound best MLAD = %f for chain %s in %s of solution %d with chain %s in %s of solution %d" \
                  %(mlad, movchainid, minoutput_file, minj, fixchainid, fixpdb, mini)
                logstr += outstr

            logstr += "\n"

        logstr += \
"""\nA good match between two chains usually have a MLAD << 1.5
whereas a poor match usually have a MLAD > 2.
The alternate origin for superpositions of solutions with more chains
should be the same for all chain pairs with a good match.\n"""
        fname = os.path.join(os.getcwd(), "AltOrigSymMLAD_%s_%s.log" %(fixlogfname, movlogfname))
        logfile= open(fname,"wb")
        logfile.write(logstr.encode('utf-8'))
        logfile.close()
        print logstr
        print "Logfile written to %s" %fname


        if movpdbs == [None] or fixpdbs == [None] :
            usage(parser.prog, accum.master_phil.as_str())
    else:
        usage(parser.prog, accum.master_phil.as_str())

    return output_file_names




def GetSolutionsAsPdbs(philfname, picklefname, mrsolfname, ensembles, xyzfname, pxtalsym):
# get list of pdb files either from pickle and phil, a phaser Mr solution file or just a pdb file supplied
    philfile = philfname
    pklfile = picklefname
    soltxts = []
    solfname = mrsolfname
    pdbfname = xyzfname
    ensmblspdbs = []
    if xyzfname is not None:
        halflogfname = os.path.basename(xyzfname).split(".")[0]

    if solfname != None and pdbfname == None and (pklfile == None or philfile == None):
        soltxts = [open(solfname,"r").read()]

    if pklfile != None and philfile != None and solfname == None and pdbfname == None:
        mstr = open(philfile, "r").read()
        user_phil = parse(mstr)
#        working_phil = phenix_interface.master_params.fetch(sources = [user_phil])
        working_phil = phenix_interface.master_phil().fetch(sources = [user_phil])
        philsxtrct = working_phil.extract()
# if using phil file get the crystal symmetry from the mtz file specified in the phil file
        if pxtalsym[0] == None:
            philmtz = philsxtrct.phaser.hklin
            mtzobj = mtz.object(philmtz)
            pxtalsym[0] = mtzobj.crystals()[0].crystal_symmetry()

        for ens in philsxtrct.phaser.ensemble:
            pdbs = []
            for coord in ens.coordinates:
                pdbs.append(coord.pdb)

            ensmblspdbs.append((ens.model_id, pdbs))
# get the MR solutions part of the pickle file
        pkl = pickle.load(open(pklfile,"rb"))
        halflogfname = os.path.basename(pklfile).split(".")[0]
        result = pkl.get_phaser_object()
        soltxts = [s.unparse() for s in result.getDotSol()]

    if solfname != None and pdbfname == None and (pklfile == None or philfile == None):
        mrparser = phaser.test_phaser_parsers.SolutionFileParser(file(solfname))
        halflogfname = os.path.basename(solfname).split(".")[0]
# make list of mrsolution strings
        soltxts = [sol.unparse() for sol in mrparser.get_dot_sol()]

        for ens in ensembles:
            ensmblspdbs.append((ens.name, ens.xyzfname))

    return (soltxts, pxtalsym, ensmblspdbs, halflogfname)




def RealizePhaserSolutions(soltxts, sym, ensmblspdbs, prefix, wpath = os.curdir):
    fnames = []
    if type(soltxts) != ListType or type(ensmblspdbs) != ListType \
     or isinstance(sym, cctbx.crystal.symmetry) == False:
        raise Sorry

    for i,soltxt in enumerate(soltxts):
        mrsolhierarchies = []
        for ensmblpdbs in ensmblspdbs:
            ensemblename = ensmblpdbs[0]

            for pdbmodelfname in ensmblpdbs[1]:# in case of more models constituting this ensemble
# if processing solutionfilepath generated on windows but now running on unix
                pdbmodelfname = pdbmodelfname.replace("\\", os.path.sep)
                pdbmrsol = pdb.input(file_name= pdbmodelfname)
                solstrs = re.findall(r"(SOLU SET.+\n(SOLU SPAC.+\n)*SOLU 6DIM.+)",soltxt, re.MULTILINE)
                ensemblesolutions = re.findall(r"SOLU 6DIM\s+ENSE\s+%s (.+)" %ensemblename, soltxt, re.MULTILINE)
                for ensemblesolution in ensemblesolutions:
                    eulers = re.findall(r"EULER\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+FRAC",
                      ensemblesolution, re.MULTILINE)[0]
                    fracs = re.findall(r"EULER\s+.+\s+FRAC\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+([-+]?[0-9]*\.?[0-9]+)\s+BFAC",
                      ensemblesolution, re.MULTILINE)[0]

                    rot = scitbx.math.euler_angles.zyz_matrix(float(eulers[0]), float(eulers[1]), float(eulers[2]))
                    trans = (float(fracs[0]), float(fracs[1]), float(fracs[2]))
                    mrsolhiearchy = pdbmrsol.construct_hierarchy().deep_copy()
                    modelatoms = mrsolhiearchy.atoms()
            # do the rotation in cartesian space
                    rotatomsites = rot * modelatoms.extract_xyz()
                    ucell = sym.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))
                    mrsolhierarchies.append( mrsolhiearchy.deep_copy() )
                    print "\nRealizing MR solution for %s:\n%s" %(ensemblename, ensemblesolution)

        fnames.append(os.path.abspath(os.path.join(wpath, prefix + "sol_" + str(i) + ".pdb")))
        jointhierarchy = pdb.hierarchy.join_roots( mrsolhierarchies, None )

        for ci, hierarchy in enumerate(jointhierarchy.models()[0].chains()):
            hierarchy.id = chainid[ci]

        mrsolpdbfile = open(fnames[-1], "w")
        mrsolpdbfile.write(jointhierarchy.as_pdb_string(sym))
        mrsolpdbfile.close()

    return fnames




def GetAltOrigSymmates(
        movingname,
        pdbfixedrefname,
        pxtalsym,
        chosenorigin = "Any",
        with_symm_operations=True,
        Debug=False,
        usebestSSM=True,
        mediocre_mlad=2.0,
        verbose=False,
        pdb_out=None,
        log_out=None,
        movehetatms=True,
        out=sys.stdout):

    if movehetatms:
        myrand = random.SystemRandom()
# label with unique random string to ensure another process running GetAltOrigSymmates
# on the same file doesn't remove the file pdbxtalmovingname
        randstr = "".join(myrand.sample(string.letters+string.digits, 5))
        pdbxtalmovingname = os.path.split(movingname)[1].split(".pdb")[0] + "_" \
         + randstr + "_movhetatms.pdb"
        sort_hetatms.run(args=["file_name=" + movingname, "output_file=" \
         + pdbxtalmovingname])
    else:
        pdbxtalmovingname = movingname
# use movingname in the log file rather than pdbxtalmovingname which may contain random characters
    xtalmovingname = os.path.split(movingname)[1].split(".pdb")[0]
    fixedrefname = os.path.split(pdbfixedrefname)[1].split(".pdb")[0]
    xtalsym = pxtalsym[0]
    pdbxtal = pdb.input(file_name= pdbxtalmovingname)
    pdbtest = pdb.input(file_name= pdbfixedrefname)
    h = tbx_utils.NullHandler()
    logger = logging.getLogger()
    logger.addHandler(h)

    sg = xtalsym.space_group()
    sgnumber = xtalsym.space_group_info().type().number()
    uc = xtalsym.unit_cell()
    movinghierarchy = pdbxtal.construct_hierarchy()
    fixedhierarchy = pdbtest.construct_hierarchy()

    nmovchains = 0
    for chain in movinghierarchy.models()[0].chains():
        if chain.is_protein()==False:
            continue
        nmovchains +=1

    nfixchains = 0
    for chain in fixedhierarchy.models()[0].chains():
        if chain.is_protein()==False:
            continue
        nfixchains +=1

    logstr =""
    stdoutstr = "%s has %d chain(s)" %(xtalmovingname, nmovchains)
    stdoutstr += "\n%s has %d chain(s)" %(fixedrefname, nfixchains)
    stdoutstr += "\nSpacegroup in %s is %s\n\n" %(xtalmovingname, xtalsym.space_group_info().symbol_and_number())
    logstr += stdoutstr + "\n"
    print >> out, stdoutstr
# bail out if unequal number of chains in moving and fixed PDB
    if nfixchains != nmovchains:
        stdoutstr = "Warning: Ambiguous input. Not the same number of chains in the two PDB files.\n"
        stdoutstr += "Results may be unexpected"
        logstr += stdoutstr + "\n"
        print >> out, stdoutstr

    AllAltOrigPdbFnames = []
    AltOrigSymMLADs = [[]]
    pdbminmladmovhierarchies = []
    pdbminmladmovchainids = []

# ------------------ looping over fixed chains -----------------------------------
    for (ifixch,thischain) in enumerate(fixedhierarchy.models()[0].chains()):
        if thischain.is_protein() == False: # ignore HETATM chains such as waters, ligands or heme groups
            continue

        testhierarchy = pdb.hierarchy.new_hierarchy_from_chain( thischain )
        fixed_chainid = testhierarchy.models()[0].chains()[0].id
        testatoms = testhierarchy.atoms()
        testcentre = GetCentreOfMolecule(testatoms)
# compute centre of coordinates in fractional space
        fractesthierarchy = testhierarchy.deep_copy()
        fractestatoms = fractesthierarchy.atoms()
        testsites_frac = uc.fractionalize(sites_cart=fractestatoms.extract_xyz())
        fractestatoms.set_xyz(testsites_frac)
# get mean fractional position of the test atoms
        testsites_frac = uc.fractionalize(sites_cart= testatoms.extract_xyz())
        meantestsitesfrac = testsites_frac.mean()
        symops = sg.all_ops()
        nsymop = len(symops)

# if there is a polar axis it will be present in any of the alternative origins
        firstorig = list_origins(sgnumber)[0][0]
# we inspect firstorig below for that
        altorigs = [([0.0, 0.0, 0.0],["0","0","0"])]
        if chosenorigin == "Any":
            altorigs = GetAllOriginsAsFloatsnStrings(sgnumber)
        else:
# parse user's explicit alternative origin string
# change a string like "1/2, 2 /3 , 5/ 6" to ["1/2", "2/3", "5/6"]
            strvec = [ s.strip() for s in chosenorigin.split(",") ]
            def zero_or_rational(v):
                try:
                    return rational.from_string(v)
                except ValueError:
                    return 0.0

            frcvc = [zero_or_rational(e) for e in strvec]
            floatvec = [float(e) for e in frcvc]
# replace components with x,y,z if there is a polar axis along those directions
            fvecpolar = [
              firstorig[i] if type(firstorig[i]) is StringType else floatvec[i] for i in range(0,3)
            ]
            strvec1 = [
              firstorig[i] if type(firstorig[i]) is StringType else strvec[i] for i in range(0,3)
            ]
# form tuple of floatvector and string vector of user's explicit alternative origin
            altorigs = [(fvecpolar, strvec1)]

        norig = len(altorigs)

        minmladstr =""
        altorigstr =""
        movhierarchyminmlad = 9e99
        movminmladstr = ""
        movminaltorigstr = ""

        pdbminmladmovhierarchy = None
# ------------------ looping over moving chains -----------------------------------
        for thischain in movinghierarchy.models()[0].chains():
            if thischain.is_protein() == False: # ignore HETATM records such as waters or heme groups
                continue

            minmlad = 9e99
            use_mmtbx_alignment = False
            ssm_crashed = False
            pdbminmlad = None

            if movehetatms:
                pdbhetatmhierarchy = GetHETATM(thischain.id, movinghierarchy)

            pdb_hierarchy = pdb.hierarchy.new_hierarchy_from_chain( thischain )
            pdb_chain_id = pdb_hierarchy.models()[0].chains()[0].id
            original_sites = pdb_hierarchy.atoms().extract_xyz()
# compute ssm alignments once, before using them in subsequent for-loops
            mov_ref_chains = [pdb_hierarchy.models()[0].chains()[0], testhierarchy.models()[0].chains()[0]]
            alignments = []
            sortQvals = []

            try:
                ssm = ccp4io_adaptbx.SecondaryStructureMatching( moving = mov_ref_chains[0],
                  reference = mov_ref_chains[1])

                Qvalues = ssm.GetQvalues()
                Qvalsindex = [(i,e,"SSM") for (i,e) in enumerate(Qvalues)]
# make list with the best SSM Q-score first
                sortQvals = sorted(Qvalsindex, key= lambda Qvindx: Qvindx[1], reverse=True)

                if usebestSSM: # then use the SSM alignment with the highest Q-score
                    sortQvals = [sortQvals[0]]
# make list of alignments and Q-scores
                for qi in sortQvals:
                    ssm.AlignSelectedMatch(qi[0])
                    alignments.append( (ccp4io_adaptbx.SSMAlignment.residue_groups(match = ssm).pairs, qi) )

            except Exception, m:
                stdoutstr = "SSM alignment of chains %s and %s failed. Trying mmtbx alignment instead." \
                  %(fixed_chainid, pdb_chain_id)
                logstr += stdoutstr + "\n"
                print >> out, stdoutstr
                use_mmtbx_alignment = True
                ssm_crashed = True

# always append an mmtbx alignment just in case a gene duplication in the SSM alignment yields a poor MLAD value
            finally:
                resmap2 = mmtbxalign(mov_ref_chains[0], mov_ref_chains[1])
                alignments.append((resmap2, (0, 1.0,"MMTBX")))

            stdoutstr = "\nIdentifying smallest MLAD between chain %s in %s and chain %s in %s" \
               %(pdb_chain_id, xtalmovingname, fixed_chainid, fixedrefname)
            stdoutstr += "\nfor the %d allowed symmetry operations and %d allowed alternative origin(s) using %d alignment(s)." \
               %(nsymop, norig, len(sortQvals))
            logstr += stdoutstr + "\n"
            print >> out, stdoutstr
# ------------------ looping over ssm alignments -----------------------------------
            for j,ali in enumerate(alignments):
# try the mmtbx alignment if the SSM alignments produces a bad MLAD value
                if not use_mmtbx_alignment and j == (len(alignments)-1):
                    break
                else:
                    if use_mmtbx_alignment and j == (len(alignments)-1) and ssm_crashed == False:
                        stdoutstr = "\nCurrent best MLAD value is poor. Trying an mmtbx alignment."
                        logstr += stdoutstr + "\n"
                        print >> out, stdoutstr

                n = 0
                ci = 0
                qi = ali[1]
                nssm = (qi[0]+1) # number of this alignment
                alimethod = qi[2]

                allmlads = []
                minshortoutstr = ""
                minaltorigsymopmlad = 9e99
                stdoutstr = "#symm_op, symm_op, #alt_origin, alt_origin, UC_shift, #%s, %s(Q), label, mlad" \
                  %(alimethod, alimethod)
                logstr += stdoutstr + "\n"
                print >> out, stdoutstr

                if Debug:
                    fname = "%s_UC%d_%s_align%d.pdb" %(xtalmovingname, n, fixedrefname, nssm)
                    allsitesfname = os.path.join(debugfolder, fname)
                    stdoutstr = "\nWriting sites of the %d. alignment to %s" %(nssm, allsitesfname)
                    logstr += stdoutstr + "\n"
                    print >> out, stdoutstr
                    allsites = open(allsitesfname, "w")
# ------------------ looping over symmetry operations -----------------------------------
                for (symcount, symop) in enumerate(symops):
# if we only want to study alternative origins then we only use the identity operation
                    if with_symm_operations == False and symcount > 0:
                        break

                    rot = symop.r()
                    trans = symop.t()

                    altorig = [0.0, 0.0, 0.0]
                    isFloatingXorig = False
                    isFloatingYorig = False
                    isFloatingZorig = False

# ------------------ looping over alternate origins -----------------------------------
                    for (oc, genuninealtorig) in enumerate(altorigs):
                        if type(genuninealtorig[0][0]) is not StringType:
                            altorig[0] = genuninealtorig[0][0]
                        else:
                            altorig[0] = 0.0
                            isFloatingXorig = True

                        if type(genuninealtorig[0][1]) is not StringType:
                            altorig[1] = genuninealtorig[0][1]
                        else:
                            altorig[1] = 0.0
                            isFloatingYorig = True

                        if type(genuninealtorig[0][2]) is not StringType:
                            altorig[2] = genuninealtorig[0][2]
                        else:
                            altorig[2] = 0.0
                            isFloatingZorig = True

                        atoms = pdb_hierarchy.atoms()
# put all xtal atoms back to original position
                        atoms.set_xyz(original_sites)
# get us into fractional coordinates
                        sites_frac = uc.fractionalize(sites_cart=atoms.extract_xyz())
# get the symmetry mate w.r.t. this space group
                        new_sites = rot.as_double() * sites_frac + trans.as_double()
# translate atom sites to alternative origin
                        altorig_new_sites = new_sites + altorig
# cast new origin back to unit cell around the test atoms
                        na = 0; nb = 0; nc = 0
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[0] > 0.5:
                            altorig_new_sites = altorig_new_sites - (1, 0, 0)
                            na -=1
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[0] < -0.5:
                            altorig_new_sites = altorig_new_sites + (1, 0, 0)
                            na +=1
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[1] > 0.5:
                            altorig_new_sites = altorig_new_sites - (0, 1, 0)
                            nb -=1
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[1] < -0.5:
                            altorig_new_sites = altorig_new_sites + (0, 1, 0)
                            nb +=1
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[2] > 0.5:
                            altorig_new_sites = altorig_new_sites - (0, 0, 1)
                            nc -=1
                        while DiffVec(altorig_new_sites.mean(), meantestsitesfrac)[2] < -0.5:
                            altorig_new_sites = altorig_new_sites + (0, 0, 1)
                            nc +=1

                        unitcellshifts = (na,nb,nc) # count how many unit cell shifts we do
# translate xtalpdb to alternative origin as well as integral numbers of unitcells vectors
                        atoms.set_xyz(altorig_new_sites)
# cast back to cartesian coordinates
                        atoms.set_xyz(uc.orthogonalize(sites_frac= altorig_new_sites))
# compute vector between alignments in fractional space for shifting atoms along floating axis
# compute distance between the two structures using ssm
                        resmap = ali[0]
# map C-alpha atoms with one another, some possibly to none
                        atommap =ensembler.to_atom_mapping(resmap, lambda atom:atom.name==' CA ', logger)
# don't use gap atoms
                        CAchains = ensembler.intersect_selection(atommap, logger)
# identify the same atom in the moved xtal
                        CAchainzip = zip(*CAchains)

                        #import code; code.interact(local=dict(globals(), **locals()))
                        xtalalignmentcentre = GetCentreOfMolecule(CAchainzip[0]) # only used if we have a floating origin
                        testalignmentcentre = GetCentreOfMolecule(CAchainzip[1]) # only used if we have a floating origin
                        xtalfraccentre = uc.fractionalize(xtalalignmentcentre)
                        fractestcentre = uc.fractionalize(testalignmentcentre)

                        mlad = GetMLADfromAtoms(CAchains)
                        fracdiff = DiffVec(xtalfraccentre, fractestcentre)
# if alternative origin has a floating axis then shift it back along this axis
                        floatingoffset = [0.0, 0.0, 0.0]
                        unitfloataxis = [0.0, 0.0, 0.0]

                        if isFloatingXorig:
                            floatingoffset[0] = fracdiff[0]
                            unitfloataxis[0] = 1.0
                        if isFloatingYorig:
                            floatingoffset[1] = fracdiff[1]
                            unitfloataxis[1] = 1.0
                        if isFloatingZorig:
                            floatingoffset[2] = fracdiff[2]
                            unitfloataxis[2] = 1.0
# place moving pdb near the reference pdb
                        cartfloatingoffset = uc.orthogonalize(floatingoffset)
                        atoms.set_xyz(atoms.extract_xyz() - tuple(cartfloatingoffset))

                        rpolarmin = 0.0
                        if isFloatingXorig or isFloatingYorig or isFloatingZorig:
                            tmpsites = atoms.extract_xyz()
# polaraxisoverlap computes GetMLADfromAtoms for this symmetry copy positioned at r along the polar axis
                            def polaraxisoverlap(r):
                                floatvec = [unitfloataxis[0]*r, unitfloataxis[1]*r,
                                  unitfloataxis[2]*r]
                                cartfloatingoffset = uc.orthogonalize(floatvec)
                                atoms.set_xyz(tmpsites - tuple(cartfloatingoffset))
                                return GetMLADfromAtoms(CAchains)

                            if Debug:
# make table of MLAD values along floating axis for each symops and origins
                                r_ini = -0.5
                                r1 = r_ini
                                mlads = [ "sym%dorig%d" %(symcount+1, oc+1)] # column label for table
                                incr = 200
                                invincr = 1.0/incr
                                for i in range(1,incr):
                                    r1 = r1 + invincr
                                    mlad = polaraxisoverlap(r1)
                                    mlads.extend([mlad])

                                allmlads.append(mlads)
# find the best overlap along the polar axis with a golden section minimisation which is fast and robust
                            rpolarmin = -9e99
                            if not Debug:
                                rpolarmin = golden_section_search.gss(polaraxisoverlap, -0.5, 0.5)
                            else:
                                dmpstr = cStringIO.StringIO()
                                rpolarmin = golden_section_search.gss(polaraxisoverlap,
                                   -0.5, 0.5, out = dmpstr, monitor_progress=True)
                                tmpstr = dmpstr.getvalue()
                                dmpstr.close()
                                # get rid of double line breaks from the gss output
                                stdoutstr = "\n".join( tmpstr.split("\n\n"))
                                logstr += stdoutstr
                                print >> out, "\n" + stdoutstr

                            mlad = polaraxisoverlap(rpolarmin)
# length of floatingoffset vector is the length of the one component along which the polar axis is aligned
                            floatoffsetlen= math.sqrt(floatingoffset[0]*floatingoffset[0] +
                             floatingoffset[1]*floatingoffset[1] + floatingoffset[2]*floatingoffset[2])
                            polaraxisminstr = "MLAD minimum= %f at fractional polar axis r= %f" %(mlad, rpolarmin + floatoffsetlen)
                            if Debug:
                                logstr += polaraxisminstr + "\n"
                                print >> out, polaraxisminstr

                        floatvec = [unitfloataxis[0]*rpolarmin + floatingoffset[0],
                          unitfloataxis[1]*rpolarmin + floatingoffset[1],
                          unitfloataxis[2]*rpolarmin + floatingoffset[2]]

                        v = flex.vec3_double([altorig]) + flex.vec3_double([floatvec])
                        altorigsymtrans = (symop, v, flex.vec3_double([unitcellshifts]))

# compute mlad between the two structures
                        mlad = GetMLADfromAtoms(CAchains)
# make a deep copy as we are about to change the chain id which otherwise would crash SSM that assumes chain ids remain constant
                        filepdbhiearchy = pdb_hierarchy.deep_copy()
                        if Debug: # because we print out all transformations
                            filepdbhiearchy.models()[0].chains()[0].id += chainid[ci]

                        symopname = symop.as_xyz()
                        altorigstr = "(%s, %s, %s)" %(genuninealtorig[1][0], genuninealtorig[1][1],
                          genuninealtorig[1][2])
                        outstr = "%d. symm.op: (%s), %d. origin: %s, UC shift: %s, %d. %s Q: %2.3f, %s" \
                          %(symcount+1, symopname, oc+1, altorigstr, str(unitcellshifts), nssm, qi[2], qi[1],
                           filepdbhiearchy.models()[0].chains()[0].id)
                        shortoutstr = "%d. (%s), %d. %s, %s, %d. %2.3f, %s, %3.3f" \
                          %(symcount+1, symopname, oc+1, altorigstr, str(unitcellshifts), nssm, qi[1],
                           filepdbhiearchy.models()[0].chains()[0].id, mlad)

                        alignstr = "\n" + str(nssm) + ". alignment. " + MakeAlignstr(atommap)

                        if minaltorigsymopmlad > mlad:
                            minaltorigsymopmlad = mlad
                            minshortoutstr = shortoutstr

                        if minmlad > mlad:
                            minmlad = mlad
                            minaltorigstr = altorigstr
                            minmladstr = outstr
                            if isFloatingXorig or isFloatingYorig or isFloatingZorig:
                                minmladstr = outstr + "\n" + polaraxisminstr + "\n"

                            pdbminmlad = filepdbhiearchy.deep_copy()
                            minmladaltorigsymtrans = deepcopy(altorigsymtrans)

                        stdoutstr = shortoutstr
                        if Debug or verbose:
                            logstr += stdoutstr + "\n"
                        print >> out, "\r" + stdoutstr + "                   ",

                        if Debug:
# place a gold atom on the centroid for the aligned moving C-alpha atoms and dump the lot as a pdb file
                            newxtalalignmentcentre = GetCentreOfMolecule(CAchainzip[0])
                            matom = pdb.hierarchy.atom()
                            matom.set_xyz((newxtalalignmentcentre[0], newxtalalignmentcentre[1],
                              newxtalalignmentcentre[2]))
                            matom.set_name("AU")
                            matom.segid = "  AU"
                            matom.hetero = True
                            CAxtalhierarchy = atoms2pdb(list(CAchainzip[0]) + [matom])

                            fname = "%s_Ca_ssm%d_sym%d_orig%d.pdb" %(xtalmovingname, nssm,
                              symcount+1, oc+1)
                            cntrfname = os.path.join(debugfolder, fname)
                            cntrout = open(cntrfname,"w")
                            cntrout.write(CAxtalhierarchy.as_pdb_string(xtalsym))
                            cntrout.close()

                            allsites.write(filepdbhiearchy.as_pdb_string(xtalsym))
                            ci += 1
                            if ci > 25: # can't exceed more than 26 chains
                                allsites.close() # close previous file
                                n += 1 # prepend filename with n++
                                ci=0 # so we can use with ABCD... again
                                fname = "%s_UC%d_%s_align%d.pdb" %(xtalmovingname, n, fixedrefname, nssm)
                                allsitesfname = os.path.join(debugfolder, fname)
                                stdoutstr = "\nWriting remaining sites of the %d. alignment to %s" \
                                  %(nssm, allsitesfname)
                                logstr += stdoutstr + "\n"
                                print >> out, stdoutstr
                                allsites = open(allsitesfname, "w")

                if Debug == False: # print intermediate symop and altorig with best mlad score to logfile
                    logstr += minshortoutstr + "\n"
                print >> out, " "

                if Debug:
                    allsites.close()
# place a gold atom on the centroid for the aligned C-alpha reference atoms and dump the lot as a pdb file
                    matom = pdb.hierarchy.atom()
                    matom.set_xyz((testalignmentcentre[0], testalignmentcentre[1],
                      testalignmentcentre[2]))
                    matom.set_name("AU")
                    matom.segid = "  AU"
                    matom.hetero = True
                    CAtesthierarchy = atoms2pdb(list(CAchainzip[1]) + [matom])
                    fname = fixedrefname + "_Ca_ssm" + str(nssm) + ".pdb"
                    refpdboutfname = os.path.join(debugfolder, fname)
                    refpdbout = open(refpdboutfname,"w")
                    refpdbout.write(CAtesthierarchy.as_pdb_string(xtalsym))
                    refpdbout.close()

                    if isFloatingXorig or isFloatingYorig or isFloatingZorig:
# print table of MLAD values along floating axis for each symops and origins
                        r1 = r_ini
                        stdoutstr = "\nMLAD values for all symmetry operations and alternative origins when sliding moving_pdb along polar axis\n\n"
                        for (i,mlad) in enumerate(allmlads[0]):
                            if i==0:
                                stdoutstr += "polar-axis "
                                for one in allmlads:
                                    stdoutstr += "%s " %one[i]
                            else:
                                stdoutstr += "%f " %r1
                                for one in allmlads:
                                    stdoutstr +=  "%2.4f " %one[i]
                            stdoutstr +=  "\n"
                            r1 = r1 + invincr

                        logstr += stdoutstr + "\n"
                        print >> out, stdoutstr

                if verbose:
                    logstr += alignstr
                    print >> out, "\n" + alignstr,

                if minmlad > mediocre_mlad:
                    use_mmtbx_alignment = True
                else:
                    use_mmtbx_alignment = False

            movhierarchyminmlad = minmlad
            tmphierarchy = pdbminmlad.deep_copy()
            minmladmovchainid = pdb_chain_id
            movminmladstr = minmladstr
            movminaltorigstr = minaltorigstr
            movminmladaltorigsymtrans = deepcopy(minmladaltorigsymtrans)
            if movehetatms and len(pdbhetatmhierarchy.models()) > 0:
                (symm, altorg, ucshift) = movminmladaltorigsymtrans
                hetatms = pdbhetatmhierarchy.atoms()
# get us into fractional coordinates
                sites_frac = uc.fractionalize(sites_cart=hetatms.extract_xyz())
# get the symmetry mate w.r.t. this space group
                new_sites = symm.r().as_double() * sites_frac + symm.t().as_double() + altorg[0] + ucshift[0]
                hetatms.set_xyz(new_sites)
# cast back into cartesian coordinates and add floatvec sliding hetatoms along polar axis if present
                hetatms.set_xyz(uc.orthogonalize(sites_frac= (hetatms.extract_xyz())))
                pdbminmladmovhierarchylst = [tmphierarchy, pdbhetatmhierarchy.deep_copy()]
            else:
                pdbminmladmovhierarchylst = [tmphierarchy]

# store MLAD scores indexed by fixed and moving chain ids in a matrix
            pdb_chain_id = pdbminmladmovhierarchylst[0].models()[0].chains()[0].id
            if movehetatms and len(pdbhetatmhierarchy.models()) > 0:
                pdbminmladmovhierarchylst[1].models()[0].chains()[0].id = pdb_chain_id
            pdbminmladmovhierarchylst[0].models()[0].chains()[0].id = minmladmovchainid

            AltOrigSymMLADs[ifixch].append([movhierarchyminmlad, minmladmovchainid,
              fixed_chainid, movminmladstr, movminaltorigstr, pdbminmladmovhierarchylst,
              movminmladaltorigsymtrans])

        if ifixch < (nfixchains-1):
            AltOrigSymMLADs.append([])

    stdoutstr = "Considering matrix of MLAD values between fixed and moving chains as to\n"\
      "ensure chain matching from fixed to moving pdb is an injective function:"
    stdoutstr += "\nUnsorted matrix:\n"

    def PrintColMatrix(mtrx):
        mstr = ""
        for f in mtrx:
            mstr += f[0][2] + ", "
            for g in f:
                mstr += "%2.3f %s%s, " %(g[0], g[1], g[2])
            mstr += "\n"
        return mstr

    stdoutstr += PrintColMatrix(AltOrigSymMLADs)
# now identify best matches between chains from fixed and moving pdbs
# first sort matrix row-wise according to MLAD value
    sortmatrix = []
    for f in AltOrigSymMLADs:
        sortrow = sorted(f, key = lambda row: row[0])
        sortmatrix.append(sortrow)

    stdoutstr += "\nSorted matrix:\n"
    stdoutstr += PrintColMatrix(sortmatrix)
# now sort according to first column MLAD value
    colsortmatrix = sorted(sortmatrix, key = lambda row: row[0][0])
    stdoutstr += "\nColumn sorted matrix:\n"
    stdoutstr += PrintColMatrix(colsortmatrix)
# ensure chain matching from fixed to moving pdb is an injective function
    injective = True
    for (i,f) in enumerate(colsortmatrix):
        for (j,g) in enumerate(colsortmatrix):
            if j > (len(f)-1): # skip moving chains that are in excess of fixed chains
                break
            if f[0][1] == g[0][1] and j > i:
                injective = False
# keep selecting next best matches as long as one chain accidentally is matching more than one
    while injective == False:
        sortmatrix = []
        for f in colsortmatrix:
            sortrow = sorted(f, key = lambda row: row[0])
            sortmatrix.append(sortrow)
        colsortmatrix = sortmatrix
# now sort according to first column MLAD value
        colsortmatrix = sorted(sortmatrix, key = lambda row: row[0][0])
# demote duplicate chainid that have a lower MLAD score from the first column
        for (i,f) in enumerate(colsortmatrix):
            for (j,g) in enumerate(colsortmatrix):
                if f[0][1] == g[0][1] and j > i:
                    g[0][0] = 100.0 # set to a high value as to demote it during next sort

        injective = True
        for (i,f) in enumerate(colsortmatrix):
            for (j,g) in enumerate(colsortmatrix):
                if j > (len(f)-1): # skip moving chains that are in excess of fixed chains
                    break
                if f[0][1] == g[0][1] and j > i:
                    injective = False

    stdoutstr += "\nDemoted column sorted matrix:\n"
    stdoutstr += PrintColMatrix(colsortmatrix)

    if Debug or verbose:
        logstr += stdoutstr + "\n"
        print >> out, stdoutstr

    for (i,f) in enumerate(colsortmatrix):
        if i > (len(f)-1): # skip moving chains that are in excess of fixed chains
            break
        AllAltOrigPdbFnames.append((f[0][0], f[0][1], f[0][2], f[0][3], f[0][4], f[0][6]))
        pdbminmladmovhierarchies.extend( f[0][5] )

    logstr += "\n"
    print >> out, "\n"
    notprinted = True

    movingnamebase = os.path.split(movingname)[1].split(".pdb")[0]
    jointpdbminmladmov = pdb.hierarchy.join_roots( pdbminmladmovhierarchies, None)
    altorigpdbfname = pdb_out
    if (pdb_out is None) :
      altorigpdbfname = "MinMLAD_%s_%s.pdb" %(fixedrefname, movingnamebase)
    altorigpdbfname = os.path.abspath(altorigpdbfname)

    resultsmtrxstr = "\nMLAD  \t%s\t%s\t  Symmetry operation  \t  Alternate origin  \t Unit cell shift\n" %(xtalmovingname, fixedrefname)
    for (mlad, movchainid, fixed_chainid, movstr, altorig, altorigsymtrans) in AllAltOrigPdbFnames:
        stdoutstr = "\nMLAD of %f for chain %s in %s with chain %s in %s using\n%s" \
           %(mlad, movchainid, movingname, fixed_chainid, pdbfixedrefname, movstr)
        logstr += stdoutstr
        print >> out, stdoutstr
        (symm, altorg, ucshift) = altorigsymtrans
        resultsmtrxstr += "%s\t   %s\t   %s\t %s\t %s\t %s\n" %(str(Roundoff(mlad,4)),
         movchainid, fixed_chainid, symm.as_xyz(), Roundoff(list(altorg[0]),4),
         Roundoff(list(ucshift[0]),4))

        if mlad > mediocre_mlad:
            stdoutstr = "These two chains are unlikely to represent the same MR solution.\n\n"
            logstr += stdoutstr
            print >> out, stdoutstr
        for f in AllAltOrigPdbFnames:
            if altorig != f[4] and notprinted:
                stdoutstr = "\nWARNING! Not all chains are matched on the same origin. " \
                 "At least one of the matches between chains is spurious!"
                logstr += stdoutstr
                print >> out, stdoutstr
                notprinted = False
        print " "
        logstr += "\n"

    print >> out, resultsmtrxstr
    logstr += "\nSuperposed MR solution saved in %s\n" %altorigpdbfname
    logstr += resultsmtrxstr

    pdbminRMSDfile = open(altorigpdbfname, "w")
    if (len(pdbxtal.remark_section()) > 0) :
      pdbminRMSDfile.write("%s\n" % "\n".join(pdbxtal.remark_section()))
    pdbminRMSDfile.write(jointpdbminmladmov.as_pdb_string(xtalsym))
    pdbminRMSDfile.close()

    if movehetatms: # delete temporary file with hetatoms
        os.remove(pdbxtalmovingname)

    return AllAltOrigPdbFnames, altorigpdbfname, logstr




def GetHETATM(id, mhierarchy): # make a pdbhierachy from hetatoms having the same chain.id as id
    hetatmhierarchies = []
    for chn in mhierarchy.models()[0].chains():
        if chn.is_protein() == False and chn.id == id:
            hetatmhierarchies.append( pdb.hierarchy.new_hierarchy_from_chain( chn ))
    hetatmhierarchy = pdb.hierarchy.join_roots( hetatmhierarchies, None)
    return hetatmhierarchy




def mmtbxalign(movingchain, fixedchain):
# use mmtbx alignment as a fall back if SSM craps out
    movres = []
    fixres = []

    for rg in movingchain.residue_groups():
        movres.append(amino_acid_codes.one_letter_given_three_letter.get(rg.unique_resnames()[0],"X"))
    movres = "".join(movres)

    for rg in fixedchain.residue_groups():
        fixres.append(amino_acid_codes.one_letter_given_three_letter.get(rg.unique_resnames()[0],"X"))
    fixres = "".join(fixres)

    myalign = mmtbx.alignment.align(movres, fixres)
    mmtbxali =  myalign.extract_alignment()
    resmap2 = []
    for i, i_seq1 in enumerate(mmtbxali.i_seqs_a):
        i_seq2 = mmtbxali.i_seqs_b[i]
        res1 = movingchain.residue_groups()[i_seq1] if i_seq1 != None else None
        res2 = fixedchain.residue_groups()[i_seq2] if i_seq2 != None else None
        resmap2.append((res1, res2))

    return resmap2




def atoms2pdb(atoms):
    r = pdb.hierarchy.root()
    r.append_model(pdb.hierarchy.model())
    r.models()[0].append_chain(pdb.hierarchy.chain(id="A"))
    r.models()[0].chains()[0].append_residue_group(pdb.hierarchy.residue_group())
    r.models()[0].chains()[0].residue_groups()[0].append_atom_group(pdb.hierarchy.atom_group())

    for (i,atm) in enumerate(atoms):
        r.models()[0].chains()[0].residue_groups()[0].append_atom_group(pdb.hierarchy.atom_group())
        r.models()[0].chains()[0].residue_groups()[0].atom_groups()[i].append_atom_with_other_parent(atm)
        if atm.hetero == False:
            r.models()[0].chains()[0].residue_groups()[0].atom_groups()[i].resname = atm.fetch_labels().resname
        else:
            r.models()[0].chains()[0].residue_groups()[0].atom_groups()[i].resname = "HEM"

        r.models()[0].chains()[0].residue_groups()[0].resseq = i

    return r




def MakeAlignstr(CAchains):
# make alignment strings of residue codes
    Ca1rescodes = ""
    Ca2rescodes = ""
    errstr = ""
    for Capair in CAchains:
        try:
            if type(Capair[0]) is not NoneType:
                Ca1rescodes += amino_acid_codes.one_letter_given_three_letter[Capair[0].parent().resname]
            else:
                Ca1rescodes += "-"
        except Exception, m:
            errstr += "Unknown residue: %s\n" %m
            Ca1rescodes += "?"

        try:
            if type(Capair[1]) is not NoneType:
                Ca2rescodes += amino_acid_codes.one_letter_given_three_letter[Capair[1].parent().resname]
            else:
                Ca2rescodes += "-"
        except Exception, m:
            errstr += "Unknown residue: %s\n" %m
            Ca2rescodes += "?"

    niceprint = ""
    n = 0
    nchar = 80
    while len(CAchains) > n:
        niceprint += Ca1rescodes[n:(n+nchar)] + "\n"
        niceprint += Ca2rescodes[n:(n+nchar)] + "\n"
        niceprint += "\n"
        n +=nchar

    alignstr = "Aligned residues:\n" + niceprint + errstr
    return alignstr




def DiffVec(vec1,vec2):
    d=[0,0,0]
    d[0] = vec1[0] - vec2[0]
    d[1] = vec1[1] - vec2[1]
    d[2] = vec1[2] - vec2[2]
    return d



def GetRMSDfromAtoms(atompairs):
    ln = len(atompairs)

    D = 0.0
    for atompair in atompairs:
        dx = atompair[0].xyz[0] - atompair[1].xyz[0]
        dy = atompair[0].xyz[1] - atompair[1].xyz[1]
        dz = atompair[0].xyz[2] - atompair[1].xyz[2]
        D += dx*dx + dy*dy + dz*dz

    return math.sqrt(D/ln)


# Mean Log Absolute Deviation
# is better than RMSD for matching overlapping chains as it downplays atom pairs which are very far apart
def GetMLADfromAtoms(atompairs):
    ln = len(atompairs)

    D = 0.0
    p = 0.9
    for atompair in atompairs:
        dx = atompair[0].xyz[0] - atompair[1].xyz[0]
        dy = atompair[0].xyz[1] - atompair[1].xyz[1]
        dz = atompair[0].xyz[2] - atompair[1].xyz[2]
        sqrD = dx*dx + dy*dy + dz*dz
        D += math.log(sqrD/(math.sqrt(sqrD) + p) + max(p,min(sqrD,1)))

    return D/ln - math.log(p) # offset to avoid values less than zero



def GetCentreOfMolecule(atoms):
    c = [0.0,0.0,0.0]
    for atom in atoms:
        c[0] += float(atom.xyz[0])
        c[1] += float(atom.xyz[1])
        c[2] += float(atom.xyz[2])

    c[0] /= len(atoms)
    c[1] /= len(atoms)
    c[2] /= len(atoms)

    return c





def make_nice( t,d):
    if t < 0:
        t=t+d
        make_nice(t,d)
    if t >= d:
        t = t-d
        make_nice(t,d)
    return t,d



def translate_v_m(v,m):
    res = [None, None ,None]
#    result = "( "
    xyz=[ 'x','y','z' ]
    isfloatorig = False
    if m==0:
    #find the first no-zero element of the vector
        this_one = 0
        for ii in range(3):
            if v[ii]!=0:
                this_one = ii
                break

        for ii in range(3):
            if v[ii]==0:
#                result += "0 "
                res[ii] = 0
            elif v[ii]==1:
#                result += str(xyz[this_one]) + " "
                res[ii] = str(xyz[this_one])
                isfloatorig = True
            else:
#                result += str(xyz[this_one])+str(v[ii]) + " "
                res[ii] = str(xyz[this_one])+str(v[ii])
                isfloatorig = True

#            if ii<2:
#                result += ", "

#        result += ")"
#        print result
        return (res, isfloatorig)

    else:
        for ii in range(3):
            if v[ii]==0:
#                result += "0 "
                res[ii] = 0
            else:
                t,d = make_nice(v[ii],m)
#                result += str(t)+"/"+str(d)
                #res[ii] = float(t)/float(d)
                #print "v[ii],t,d= %s, %s, %s" %(v[ii], t, d)
                #print "type v[ii],t,d= %s, %s, %s" %(type(v[ii]),type(t), type(d))
                if type(t) == int:
                    t= rational.int(t,1)
                if type(d) == int:
                    d= rational.int(d,1)

                obj = t/d
                res[ii] = Fraction(obj.numerator(), obj.denominator())

#            if ii<2:
#                result += ", "

#        result += ")"
#        print result
        return (res, isfloatorig)





#get alternative origins for spage group
def list_origins(sgnumbersymbol):
    sginfo = sgtbx.space_group_info(symbol=sgnumbersymbol, space_group_t_den=144)
    #print
    #print
    #sginfo.show_summary()
    #print "Allowed additional origins:"
    vecs = sginfo.structure_seminvariants().vectors_and_moduli()

    sg_type = sginfo.type()
    k2l = sg_type.addl_generators_of_euclidean_normalizer(True, False)

    results = []

    if len(vecs)>0:
        for ii in vecs:
            r = translate_v_m(ii.v, ii.m)
            if r not in results:
                results.append(r)
            """
            if ii.m is not 0:
                if (len(k2l)):
                    inv = sgtbx.translation_part_info(k2l[0]).origin_shift().as_rational()
                    nm = []
                    for a,b in zip(ii.v,inv):
                        nm.append( -(a-b) )

                    #print "k2l[0], inv, nm = %s, %s %s" %(k2l[0], inv, nm)
                    r = translate_v_m(nm, ii.m)
                    if r not in results:
                        results.append(r)
            """
    hasfloatingorigin = False
    resultsoriginf = [[0.0, 0.0, 0.0]] # when only centering vectors provide allowed origin shifts
    for res in results:
        if res[1]:
            hasfloatingorigin = True
        resultsoriginf.append(res[0])

    return (resultsoriginf, hasfloatingorigin)



# cast x,y,z components within the interval [0;1]
def Remainder(v):
    t = list(v)
    if type(t[0]) is not StringType:
        t[0] = divmod(t[0], 1)[1]
    if type(t[1]) is not StringType:
        t[1] = divmod(t[1], 1)[1]
    if type(t[2]) is not StringType:
        t[2] = divmod(t[2], 1)[1]
    return t



def GenerateEquivalentOriginsForVector(origvlst):
    allorigins = set([])
    for origv in origvlst:
# equivalent origins are those that are translated along each origv
        modnewv = tuple(Remainder(origv))
        n=1
        while modnewv not in allorigins:
            allorigins.add(modnewv)
            newv = [0,0,0]
            n += 1

            if type(origv[0]) is StringType:
                newv[0] = origv[0]
            else:
                newv[0] = n * origv[0]

            if type(origv[1]) is StringType:
                newv[1] = origv[1]
            else:
                newv[1] = n * origv[1]

            if type(origv[2]) is StringType:
                newv[2] = origv[2]
            else:
                newv[2] = n * origv[2]

            modnewv = tuple(Remainder(newv))

    return allorigins



def GetOriginVectorCombinations(originvlst):
    origvClst = []
    for n in range(1,len(originvlst)+1):
# make a K combination for each n
        Klst = list( itertools.combinations(originvlst,n))

# want origins translated along the combination of these individual vectors, i.e. the sum of each vector in each combination
        for K in Klst:
            combinedorigv = [0,0,0]
            for v in K:
                if type(v[0]) is StringType:
# if this is either x, y or z then it is a floating origin along that axis
                    combinedorigv[0] = v[0]
                else:
                    if type(combinedorigv[0]) is not StringType:
                        combinedorigv[0] += v[0]

                if type(v[1]) is StringType:
# if this is either x, y or z then it is a floating origin along that axis
                    combinedorigv[1] = v[1]
                else:
                    if type(combinedorigv[1]) is not StringType:
                        combinedorigv[1] += v[1]

                if type(v[2]) is StringType:
# if this is either x, y or z then it is a floating origin along that axis
                    combinedorigv[2] = v[2]
                else:
                    if type(combinedorigv[2]) is not StringType:
                        combinedorigv[2] += v[2]

            origvClst.append(combinedorigv)

    return origvClst




def GetCentringVectors(sgnumber):
    cvecs = []
    sginfo = sgtbx.space_group_info(symbol=sgnumber, space_group_t_den=144)

    if sginfo.type().lookup_symbol()[0] == "A":
        cvecs.append([0.0, Fraction(1,2), Fraction(1,2)])

    if sginfo.type().lookup_symbol()[0] == "C":
        cvecs.append([Fraction(1,2), 0.0, Fraction(1,2)])

    if sginfo.type().lookup_symbol()[0] == "C":
        cvecs.append([Fraction(1,2), Fraction(1,2), 0.0])

    if sginfo.type().lookup_symbol()[0] == "I":
        cvecs.append([Fraction(1,2), Fraction(1,2), Fraction(1,2)])

    if sginfo.type().lookup_symbol()[0] == "F":
        cvecs.append([Fraction(1,2), Fraction(1,2), 0.0])
        cvecs.append([0.0, Fraction(1,2), Fraction(1,2)])
        cvecs.append([Fraction(1,2), 0.0, Fraction(1,2)])

    if sginfo.type().lookup_symbol()[0] == "R":
        # centred rhombohedral: obverse setting
        cvecs.append([Fraction(2,3), Fraction(1,3), Fraction(1,3)])
        cvecs.append([Fraction(1,3), Fraction(2,3), Fraction(2,3)])
        # centred rhombohedral: reverse setting is unconventional so bin it for now
        #cvecs.append([Fraction(1,3), Fraction(2,3), Fraction(1,3)])
        #cvecs.append([Fraction(2,3), Fraction(1,3), Fraction(2,3)])

    if sginfo.type().lookup_symbol()[0] == "H":
        # centred hexagonal
        cvecs.append([Fraction(2,3), Fraction(1,3), 0.0])
        cvecs.append([Fraction(1,3), Fraction(2,3), 0.0])

    return cvecs




def GetAllOrigins(sgnumber):
    origlst = list_origins(sgnumber)[0]
    origlst.extend(GetCentringVectors(sgnumber))
    Coriglst = GetOriginVectorCombinations(origlst)
# impose the floating origin on all other vectors if present
    origfloatX = False; origfloatY = False; origfloatZ = False
    for elm in Coriglst:
        if elm[0] == "x":
            origfloatX = True
        if elm[1] == "y":
            origfloatY = True
        if elm[2] == "z":
            origfloatZ = True

    for elm in Coriglst:
        if origfloatX:
            elm[0] = "x"

    for elm in Coriglst:
        if origfloatY:
            elm[1] = "y"

    for elm in Coriglst:
        if origfloatZ:
            elm[2] = "z"

    eqoriglst = GenerateEquivalentOriginsForVector(Coriglst)

    return list(eqoriglst)




def ListAsFloats(vs):
    lfv = []
    for v in vs:
        fv = [0,0,0]
        if type(v[0]) is not StringType:
            fv[0] = float(v[0])
        else:
            fv[0] = v[0]

        if type(v[1]) is not StringType:
            fv[1] = float(v[1])
        else:
            fv[1] = v[1]

        if type(v[2]) is not StringType:
            fv[2] = float(v[2])
        else:
            fv[2] = v[2]

        lfv.append(fv)

    return lfv





def GetAllOriginsAsFloatsnStrings(sgnumber):
    lst = GetAllOrigins(sgnumber)
    floatlst = ListAsFloats(lst)
    strlst = [[str(e[0]), str(e[1]), str(e[2])] for e in lst]
    fslst = zip(floatlst, strlst)
    return fslst




def PrettyPrintAllOrigins(sgnumber):
    lst = GetAllOrigins(sgnumber)
    strlst = [("(" + str(e[0]) + ", " + str(e[1]) + ", " + str(e[2]) + ")") for e in lst]
    for l in strlst:
        print l




if __name__ == "__main__" :
    RunFromCommandline(sys.argv)