#!/usr/bin/env python
# -*- coding: utf-8 -*-
#

import sys
import os
import glob
import time
import itertools
import ctypes
from numpy.ctypeslib import ndpointer
from ctypes import c_double

from heapq import heappush, heappop

# Compatibility with Python < 2.6
#
try:
  from heapq import heappushpop
except ImportError:
  def heappushpop(heap, item):
    heappush(heap, item)
    return heappop(heap)

import numpy
import math
import pickle

from .protein import Pdb, ResidueList, Protein
from .pdb3d import dist
from .geohash import GeometricHash

from . import ccd
from .loopmodel import ANCHOR_LENGTH, get_max_loop_length, describe_anchors, iterate_database, get_loop_structure, add_oxygens, is_loop_length_in_database, get_min_loop_length_in_database, relabel_loop, is_clash, calculate_rmsd, is_consecutive, make_contact_gh, find_contacts, get_contact_class, find_contacts_simple, get_native_contacts
from .esst import load_esst

# Define consecutive residues very loosely, for purposes of checking the framework
def is_framework_consecutive(a, b):
  # Set max distance for consecutive C and N atoms to 2.8 Angstroms
  # ... which is very loose
  #
  return is_consecutive(a, b, maxlen_sq=7.84)



def _heapadd(heap, item, maxitems):
  if len(heap) < maxitems:
    heappush(heap, item)
  else:
    heappushpop(heap, item)



class FreadError(RuntimeError):
  pass

class BadInputError(FreadError):
  pass

class UnsupportedLoopLengthError(BadInputError):
  pass

class LoopStretchError(BadInputError):
  pass

class NonConsecutiveAnchorError(BadInputError):
  pass
  
class Data:
  pass


class Fread(object):
  def __init__(self, db=None, subst_table_path=None, score_cutoff=25, open_rmsd_cutoff=1.0, verbose=False, errstream=sys.stderr):
    self.score_cutoff = score_cutoff
    self.open_rmsd_cutoff = open_rmsd_cutoff
    self.verbose = verbose
    
    self.warnings = []
    
    if db:
      self.set_db(db)
    
    self.set_subst_tables(subst_table_path)
    
    self.results = []
    self.counter = itertools.count(1)
    self.idecoy = 0
  
  def set_options(self, opts):
    for o in opts:
      if o in self.__dict__ and not o.startswith('_'):
        self.__dict__[o] = opts[o]
  
  def set_silent(self):
    #self.verbose = False
    self.errstream = open(os.devnull, 'w')
  
  def note(self, message, force=False, prefix=""):
    if force or self.verbose:
      self.errstream.write(prefix+str(message)+"\n")

  def warn(self, exception, force=False, prefix="WARNING: "):
    if isinstance(exception, str):
      exception = FreadError(exception)
    self.note(str(exception), force, prefix)
    self.warnings.append(exception)
    return exception

  def die(self, exception, force=True, prefix="ERROR: "):
    exc = self.warn(exception, force, prefix)
    exit()
  
  
  def set_subst_tables(self, subst_table_path):
    self.subst_tables = load_esst(subst_table_path)
  
  def set_db(self, db):
    self.min_database_loop_length = get_min_loop_length_in_database(db)
    self.db = db
  
  def set_structure(self, structure):
    self.clear_results()
    p = Pdb(structure)
    self.ligand_atoms = p.ligands
    self.residues = ResidueList(p)
    add_oxygens(self.residues)
  
  
  def get_loop_index(self, residue_number, inscode="", chain=None):
    if inscode:
      inscode = inscode.strip()
    if chain:
      chain = chain.strip()
    
    if isinstance(residue_number, str) and residue_number[-1].isalpha():
        ires = int(residue_number[:-1])
        inscode = residue_number[-1]
    else:
        ires = int(residue_number)
    
    start_of_loop=None
    for i,r in enumerate(self.residues):
      if r.ires == ires and r.inscode == inscode and (not chain or chain == r.chain):
          start_of_loop = i+1
          break
    
    return start_of_loop
  
  
  def get_structure_sequence(self, startindex=0, endindex=None):
    return self.residues[startindex:endindex].get_seq()
  
  
  def set_loop(self, start_of_loop, loop_sequence, resnum=False, chain=None):
    self.clear_results()
    if chain is not None:
      chain = chain.strip()
    
    # We've been given a residue number instead of a start index
    if resnum:
      if isinstance(start_of_loop, str):
        if start_of_loop[-1].isalpha():
          ires = int(start_of_loop[:-1])
          inscode = start_of_loop[-1]
        else:
          ires = int(start_of_loop)
          inscode = ""
      else:
        ires = int(start_of_loop)
        inscode = ""
      
      found=False
      for i,r in enumerate(self.residues):
        if r.ires == ires and r.inscode == inscode and (not chain or chain == r.chain):
            start_of_loop = i+1
            found = True
            self.note("Located residue left of loop to be modelled: %d%s. Loop starting index: %d. Sequence surrounding start of loop: %s|%s" % (r.ires, r.inscode, start_of_loop, self.residues[max(0,start_of_loop-8):start_of_loop].get_seq(), self.residues[start_of_loop:start_of_loop+8].get_seq()))
            break
      
      if not found:
        msg = "Residue before start of loop (residue number %d%s" % (ires, inscode)
        if chain:
          msg += ", chain '%s'" % chain
        msg += ") not found in query structure"
        self.die(BadInputError(msg))
    
    if start_of_loop < ANCHOR_LENGTH or start_of_loop > len(self.residues) - ANCHOR_LENGTH:
      self.die(BadInputError("Cannot model loop closer than %d residues to the terminus (need a complete anchor on both sides)." % (ANCHOR_LENGTH)))
    
    # Assuming no coordinates are present in file, start and end of the loop are the same
    end_of_loop = start_of_loop
    
    # Are the coordinates of the loop in the file?
    self.loop_structure_is_known = is_framework_consecutive(self.residues[start_of_loop-1],
                                                    self.residues[start_of_loop])
    self.note("Is loop structure present: "+str(self.loop_structure_is_known))
    
    # If we have the native loop structure, adjust end_of_loop.
    # Also do some additional sanity checks
    #
    if self.loop_structure_is_known:
      end_of_loop += len(loop_sequence)
      strucseq = ""
      for i,r in enumerate(self.residues[start_of_loop:end_of_loop]):
          strucseq += r.get_seq()
      if (loop_sequence != strucseq):
        self.die(BadInputError("Residues differ between sequence and structure input: %s, %s" % (loop_sequence, strucseq)))
      self.note("Loop sequence of given ATOM co-ordinates: %s" % (strucseq))
      del strucseq
    
    if end_of_loop != start_of_loop:
      if len(loop_sequence) != end_of_loop - start_of_loop:
        self.die(BadInputError("Loop co-ordinates present in input, but number of residues (%d) does not match length of input loop sequence (%d)." % (end_of_loop - start_of_loop, len(loop_sequence))))
      if end_of_loop > len(self.residues) - ANCHOR_LENGTH:
        self.die(BadInputError("Cannot model loop closer than %d residues to the terminus (need a complete anchor on both sides)." % (ANCHOR_LENGTH)))
    
    self.seq = self.loop_sequence = loop_sequence
    self.s = self.start_of_loop = start_of_loop
    self.e = self.end_of_loop = end_of_loop
    
    self.verify_anchors()
  
  
  def verify_anchors(self):
    # Ensure anchors are consecutive stretches of amino acids
    for x in range(self.s-ANCHOR_LENGTH+1, self.s):
      if not is_framework_consecutive(self.residues[x-1], self.residues[x]):
        self.die(NonConsecutiveAnchorError("Anchor residues not consecutive in framework structure: residue index %s, %s"%(self.s, self.residues.code)))
    for x in range(self.e, self.e+ANCHOR_LENGTH-1):
      if not is_framework_consecutive(self.residues[x], self.residues[x+1]):
        self.die(NonConsecutiveAnchorError("Anchor residues not consecutive in framework structure: residue index %s, %s"%(self.s, self.residues.code)))
  
  def verify_stretch(self):
    # Ensure anchors are not too far apart for loop to stretch gap
    anchor_distance = dist(self.residues[self.s-1].C, self.residues[self.e].N)
    if anchor_distance > get_max_loop_length(len(self.seq)) * 1.05:
      self.die(LoopStretchError("Loop span (%.2f Angstrom) too large to be closed using %d residues. Trying to extend gap." % (anchor_distance, len(self.seq))))
  
  def verify_length(self):
    # Ensure loop is not too long for the database    
    if not is_loop_length_in_database(self.db, len(self.seq)):
      self.die(UnsupportedLoopLengthError("Cannot model loop of length %d, due to database limitations"%(len(self.seq))))
  
  def verify(self):
    self.verify_anchors()
    self.verify_stretch()
    self.verify_length()  
  
  def model(self, pdb_file, MECHANOstart, MECHANOend, results_file, decoy_directory, chain, scriptpath, no_hits, ignore_list = [], sphinx_version="general", ranking_method="soaploop", sidechain_method="pears"):
              
    residues = self.residues
                
    log_file = open(results_file, "w")
    
    log_file.write("Sphinx log file\n\n"
                   "    Structure        : %s\n"
                   #"    PDB file         : %s\n"
                   #"    Database         : %s\n"
                   "    Start residue    : %s\n"
                   "    End residue      : %s\n"
                   "    Chain            : %s\n"
                   "    Sequence         : %s\n" 
                   "    Start Time       : %s\n"
                   "    Sphinx Version   : %s\n"
                   "    Ranking Method   : %s\n"
                   "    Sidechain Method : %s\n"
                #%(self.residues.code, pdb_file, self.db, MECHANOstart, MECHANOend, chain, 
                %(self.residues.code, MECHANOstart, MECHANOend, chain,
                self.seq, time.strftime("%c"), sphinx_version, ranking_method, sidechain_method))
                   
    log_file.close() 
    
    MECHANO_PATH = os.path.join( scriptpath, "mechano" )
    # Hacky fix!
    so_file = glob.glob(os.path.join(MECHANO_PATH, "MECHANO*.so"))[0]
    MECHANO = ctypes.CDLL( so_file )
      
    def calculate_dihedral(A,B,C,D):
      b1 = B-A
      b2 = C-B
      b3 = D-C
                    
      n1 = numpy.cross(b1,b2)
      n2 = numpy.cross(b2,b3)
                  
      n1_unit = n1/numpy.linalg.norm(n1)
      n2_unit = n2/numpy.linalg.norm(n2)
                    
      b2_unit = b2/numpy.linalg.norm(b2)
      m1 = numpy.cross(n1,b2_unit)
                    
      x = numpy.dot(n1,n2)
      y = numpy.dot(m1,n2)
                    
      angle = -math.atan2(y,x) * 180 / math.pi
    
      return angle      
      
    """Model loop using the FREAD algorithm. This is the method handling most of the work. Raises UnsupportedLoopLengthError if loop length is not supported.
    """
    
    self.verify()
    
    verbose = self.verbose
    start_of_loop = self.s
    end_of_loop = self.e
    loop_sequence = self.seq
    
    start = start_of_loop - ANCHOR_LENGTH
    end = end_of_loop + ANCHOR_LENGTH
    loop_length = len(loop_sequence)
    total_length = loop_length + 2*ANCHOR_LENGTH
    
    self.note("Loop region: N(%4d, %4d)  C(%4d,%4d)  Loop(%3d,'%s')" % (start, start_of_loop, end_of_loop, end, loop_length, loop_sequence))
    
    
    ############################################################################
    
                                 # Get query anchors #
    
    
    # Get anchor coordinates
    #
    anchor_N = residues[start:start_of_loop]
    anchor_C = residues[end_of_loop:end]
    
    # Describe anchors in query structure
    #
    anchor_description, query_transform = describe_anchors(anchor_N, anchor_C, loop_length)    
    
    ############################################################################
    
    #Create coordinate arrays for MECHANO
    native_coordinates = []
    C_anchor_coordinates = []
    N_anchor_coordinates = []
    
    structure_coordinates = Data()
    structure_coordinates.N = []
    structure_coordinates.CA = []
    structure_coordinates.C = []
    structure_coordinates.O = []
    
    if self.loop_structure_is_known:
        self.actual_loop_start = str(residues[start_of_loop].ires) + residues[start_of_loop].inscode
    else:
        self.actual_loop_start = residues[start_of_loop-1].ires + 1
    
    for i in residues[start_of_loop:end_of_loop]:
        native_coordinates.extend(list(i.N.xyz))
        native_coordinates.extend(list(i.CA.xyz))
        native_coordinates.extend(list(i.C.xyz))
        native_coordinates.extend(list(i.O.xyz))        
        
    for i in residues[start_of_loop-3:start_of_loop]:
        N_anchor_coordinates.extend(list(i.N.xyz))
        N_anchor_coordinates.extend(list(i.CA.xyz))
        N_anchor_coordinates.extend(list(i.C.xyz))
        N_anchor_coordinates.extend(list(i.O.xyz))
        
    for i in residues[end_of_loop:end_of_loop+3]:
        C_anchor_coordinates.extend(list(i.N.xyz))
        C_anchor_coordinates.extend(list(i.CA.xyz))
        C_anchor_coordinates.extend(list(i.C.xyz))
        C_anchor_coordinates.extend(list(i.O.xyz))
            
    for i in residues[:start_of_loop-1]:
        for atom in i.iter_backbone():
            if atom.atom == "N":
                structure_coordinates.N.extend(list(i.N.xyz))
            if atom.atom == "CA":
                structure_coordinates.CA.extend(list(i.CA.xyz))
            if atom.atom == "C":
                structure_coordinates.C.extend(list(i.C.xyz))
            if atom.atom == "O":
                structure_coordinates.O.extend(list(i.O.xyz))
        
    for i in residues[end_of_loop+1:]:
        for atom in i.iter_backbone():
            if atom.atom == "N":
                structure_coordinates.N.extend(list(i.N.xyz))
            if atom.atom == "CA":
                structure_coordinates.CA.extend(list(i.CA.xyz))
            if atom.atom == "C":
                structure_coordinates.C.extend(list(i.C.xyz))
            if atom.atom == "O":
                structure_coordinates.O.extend(list(i.O.xyz))
        
    atom_frequencies = [0,0,0,0]
    
    atom_frequencies[0] = len(structure_coordinates.N)//3
    atom_frequencies[1] = len(structure_coordinates.CA)//3
    atom_frequencies[2] = len(structure_coordinates.C)//3
    atom_frequencies[3] = len(structure_coordinates.O)//3
    
    info_array = [0,0]
    
    if self.loop_structure_is_known:
        info_array[0] = 1
        
    blah = residues[start_of_loop-1]
    
    if blah.rest != []:
        info_array[1] = blah.rest[-1].iatom
    elif blah.res == "GLY":
        info_array[1] = blah.O.iatom
    else:
        info_array[1] = blah.CB.iatom

    C_native_coords  = (ctypes.c_double * len(native_coordinates))(*native_coordinates)
    C_Nanchor_coords = (ctypes.c_double * len(N_anchor_coordinates))(*N_anchor_coordinates)
    C_Canchor_coords = (ctypes.c_double * len(C_anchor_coordinates))(*C_anchor_coordinates)
    
    C_Ncoords  = (ctypes.c_double * len(structure_coordinates.N))(*structure_coordinates.N)
    C_CAcoords  = (ctypes.c_double * len(structure_coordinates.CA))(*structure_coordinates.CA)
    C_Ccoords  = (ctypes.c_double * len(structure_coordinates.C))(*structure_coordinates.C)
    C_Ocoords  = (ctypes.c_double * len(structure_coordinates.O))(*structure_coordinates.O)
    
    C_atom_freqs  = (ctypes.c_int * len(atom_frequencies))(*atom_frequencies)
    
    C_info_array = (ctypes.c_int * len(info_array))(*info_array)
    
    N_seq = residues[start_of_loop-3:start_of_loop].get_seq()
    C_seq = residues[end_of_loop:end_of_loop+3].get_seq() 

    ################################# Search the database #################################
    
    tables = self.subst_tables.tables
    ascii2index = self.subst_tables.ascii2index

    #######################################################################################
    
    no_decoys = 100
    
    if len(self.seq) <= 12:
        length_limit = 3
    else:
        length_limit = 4
        
#     if len(self.seq)%2 == 0:
#         no_hits = (len(self.seq)-2) * 25 / 2
#     else:
#         no_hits = (len(self.seq)-1) * 25 / 2
    
    results_count = {}    
    dihedrals_dictionary = {}
    dihedrals_set = {}

    lower_limit = max(3, len(self.seq)-length_limit)
    upper_limit = len(self.seq)
    
#    already_done = []
#    previous = open("/data/cockatrice/not-backed-up/marks/Sphinx-General/Training/LoopsOnlyDatabase/%s_%s_%s_log.txt" %(residues.code, str(self.actual_loop_start), MECHANOend)).readlines()
#    for line in previous:
#        if line.startswith(">"):
#            already_done.append(line.rstrip().split()[-1])
     
    all_decoys = []
    for loop_len in range(lower_limit,upper_limit):        
        total_length = loop_len + 2*ANCHOR_LENGTH
    
        for decoy in iterate_database(self.db, loop_len, self.subst_tables, anchor_description, loop_sequence, self.open_rmsd_cutoff, self.score_cutoff):
            if str(decoy.struc)[:4] in ignore_list:
                continue
            else:
                all_decoys.append(decoy)
    
    sorted_decoys = sorted(all_decoys, key=lambda x: x.score, reverse=True)

    #Get score of the nth decoy (where n = number of hits required, ranked by ESSS)
    if len(sorted_decoys) < no_hits:
        chosen_decoys = [decoy for decoy in all_decoys if decoy.score >= 25]
    else:
        min_score = sorted_decoys[no_hits-1].score
        min_score = max(25, min_score)
        chosen_decoys = [decoy for decoy in all_decoys if decoy.score >= min_score]
    
    decoy_file = str(os.path.join(decoy_directory, "all_decoy_structures.pdb"))
    self.decoy_file = decoy_file
    f = open(decoy_file, "w")
    f.close()
    
    count = 0
     
    #required_decoys = [str("%s_%d_%d" %(decoy.struc, decoy.start, len(decoy.seq))) for decoy in chosen_decoys]
    #write_list_to_file("/data/cockatrice/not-backed-up/marks/Sphinx-General/Training/LoopsOnlyDatabase/DoubleFragments/%s_%s_%s_required_decoys.txt" %(residues.code, MECHANOstart, MECHANOend), required_decoys)
    
    for decoy in chosen_decoys:
        count += 1
        
        #decoy_id = str("%s_%d_%d" %(decoy.struc, decoy.start, len(decoy.seq)))
        #if decoy_id in already_done:
        #    continue
     
        # Retrieve loop structure from database
        # Decoy residues + 2 anchors either side
        
        try:
          decoy_residues = get_loop_structure(self.db, decoy.struc, decoy.start, len(decoy.seq)+ANCHOR_LENGTH*2)
        except IOError:
          self.note("%s_%d_%d : database inconsistency : protein structure not found in db '%s'"%(decoy.struc, decoy.start, loop_length, self.db))
          continue
        
        assert len(decoy_residues) == len(decoy.seq)+ANCHOR_LENGTH*2
        
        # Superimpose anchors and check anchor RMSD before starting
        anchor_rmsd_open = ccd.superimpose(decoy_residues[:ANCHOR_LENGTH]+decoy_residues[-ANCHOR_LENGTH:], anchor_N+anchor_C, decoy_residues)
          
        #Find which residues we are going to copy
        alignment = {}
        target_index = 0
        template_index = 2
        for i in range(len(decoy.aligned1)):
          if decoy.aligned1[i] == "-":
            template_index += 1
            continue
          elif decoy.aligned2[i] == "-":
            alignment[target_index] = None
            target_index += 1
          elif decoy.aligned1[i] != "-" and decoy.aligned2[i] != "-":
            alignment[target_index] = template_index
            target_index += 1
            template_index += 1
            
        bitstring = ""
        for i in range(len(loop_sequence)):
          if alignment[i] == None:
            bitstring += "0"
          else:
            bitstring += "1"
                          
        #--------------------------------------------------------------------------------------------------       
        #--------------------------------------------------------------------------------------------------
        #Calculate bond lengths, angles and dihedrals for aligned residues
        #building_data = dictionary containing bond lengths, angles and dihedrals for all loop residues and
        #one anchor residue either side
        building_data = {}
        for i in range(1,len(decoy_residues)-1):
          building_data[i] = Data()
          
          prev = decoy_residues[i-1]
          curr = decoy_residues[i]
          next = decoy_residues[i+1]
              
          cac = prev.CA.xyz - prev.C.xyz
          cN  = prev.C.xyz - curr.N.xyz
          NCA = curr.N.xyz - curr.CA.xyz
          CAC = curr.CA.xyz - curr.C.xyz
          CO  = curr.C.xyz - curr.O.xyz
          Cn  = curr.C.xyz - next.N.xyz 
          
          #Calculate bond lengths
          building_data[i].length1 = numpy.linalg.norm(cN)
          building_data[i].length2 = numpy.linalg.norm(NCA)
          building_data[i].length3 = numpy.linalg.norm(CAC)
          building_data[i].length4 = numpy.linalg.norm(CO)
              
          #Calculate bond angles
          building_data[i].angle1 = numpy.arccos(numpy.dot(-cac,cN)  / (numpy.linalg.norm(cac) * numpy.linalg.norm(cN)))
          building_data[i].angle2 = numpy.arccos(numpy.dot(-cN,NCA)  / (numpy.linalg.norm(cN)  * numpy.linalg.norm(NCA)))
          building_data[i].angle3 = numpy.arccos(numpy.dot(-NCA,CAC) / (numpy.linalg.norm(NCA) * numpy.linalg.norm(CAC)))
          building_data[i].angle4 = numpy.arccos(numpy.dot(-CAC,CO)  / (numpy.linalg.norm(CAC) * numpy.linalg.norm(CO)))
          building_data[i].angle5 = numpy.arccos(numpy.dot(Cn,CO)    / (numpy.linalg.norm(Cn)  * numpy.linalg.norm(CO)))
          
          #Calculate dihedral angles
          building_data[i].phi   = calculate_dihedral(prev.C.xyz, curr.N.xyz, curr.CA.xyz, curr.C.xyz)
          building_data[i].psi   = calculate_dihedral(curr.N.xyz, curr.CA.xyz, curr.C.xyz, next.N.xyz)
          building_data[i].omega = calculate_dihedral(curr.CA.xyz, curr.C.xyz, next.N.xyz, next.CA.xyz)
          building_data[i].oxy   = calculate_dihedral(curr.N.xyz, curr.CA.xyz, curr.C.xyz, curr.O.xyz)
          building_data[i].oxy2  = calculate_dihedral(curr.O.xyz, curr.C.xyz, next.N.xyz, next.CA.xyz)

        #N anchor data (from ACTUAL coordinates)
        prev = anchor_N[0]
        curr = anchor_N[1]
        
        cac = prev.CA.xyz - prev.C.xyz
        cN  = prev.C.xyz - curr.N.xyz
        NCA = curr.N.xyz - curr.CA.xyz
        CAC = curr.CA.xyz - curr.C.xyz
        CO  = curr.C.xyz - curr.O.xyz
        
        building_data["realN"] = Data()
        
        building_data["realN"].length1 = numpy.linalg.norm(cN)
        building_data["realN"].length2 = numpy.linalg.norm(NCA)
        building_data["realN"].length3 = numpy.linalg.norm(CAC)
        building_data["realN"].length4 = numpy.linalg.norm(CO)
        
        building_data["realN"].angle1 = numpy.arccos(numpy.dot(-cac,cN)  / (numpy.linalg.norm(cac)*numpy.linalg.norm(cN)))
        building_data["realN"].angle2 = numpy.arccos(numpy.dot(-cN,NCA)  / (numpy.linalg.norm(cN)*numpy.linalg.norm(NCA)))
        building_data["realN"].angle3 = numpy.arccos(numpy.dot(-NCA,CAC) / (numpy.linalg.norm(NCA)*numpy.linalg.norm(CAC)))
        building_data["realN"].angle4 = numpy.arccos(numpy.dot(-CAC,CO)  / (numpy.linalg.norm(CAC)*numpy.linalg.norm(CO)))
        
        building_data["realN"].phi   = calculate_dihedral(prev.C.xyz, curr.N.xyz, curr.CA.xyz, curr.C.xyz)
        building_data["realN"].oxy   = calculate_dihedral(curr.N.xyz, curr.CA.xyz, curr.C.xyz, curr.O.xyz)
        
        
        #C anchor data (from ACTUAL coordinates)
        curr = anchor_C[0]
        next = anchor_C[1]
   
        cN  = prev.C.xyz - curr.N.xyz
        NCA = curr.N.xyz - curr.CA.xyz
        CAC = curr.CA.xyz - curr.C.xyz
        CO  = curr.C.xyz - curr.O.xyz
        Cn  = curr.C.xyz - next.N.xyz 
        
        building_data["realC"] = Data()
      
        #Calculate bond lengths
        building_data["realC"].length1 = numpy.linalg.norm(cN)
        building_data["realC"].length2 = numpy.linalg.norm(NCA)
        building_data["realC"].length3 = numpy.linalg.norm(CAC)
        building_data["realC"].length4 = numpy.linalg.norm(CO)
          
        #Calculate bond angles
        building_data["realC"].angle3 = numpy.arccos(numpy.dot(-NCA,CAC) / (numpy.linalg.norm(NCA)*numpy.linalg.norm(CAC)))
        building_data["realC"].angle4 = numpy.arccos(numpy.dot(-CAC,CO)  / (numpy.linalg.norm(CAC)*numpy.linalg.norm(CO)))
        building_data["realC"].angle5 = numpy.arccos(numpy.dot(Cn,CO)   / (numpy.linalg.norm(Cn)*numpy.linalg.norm(CO)))
      
        #Calculate dihedral angles
        building_data["realC"].psi   = calculate_dihedral(curr.N.xyz, curr.CA.xyz, curr.C.xyz, next.N.xyz)
        building_data["realC"].omega = calculate_dihedral(curr.CA.xyz, curr.C.xyz, next.N.xyz, next.CA.xyz)
        building_data["realC"].oxy   = calculate_dihedral(curr.N.xyz, curr.CA.xyz, curr.C.xyz, curr.O.xyz)
        building_data["realC"].oxy2  = calculate_dihedral(curr.O.xyz, curr.C.xyz, next.N.xyz, next.CA.xyz)
                    
        decoy_id = str("%s_%d_%d" %(decoy.struc, decoy.start, len(decoy.seq)))
            
        #--------------------------------------------------------------------------------------------------
        #--------------------------------------------------------------------------------------------------
        #Initialise MECHANO arrays 
          
        MECHANO_arrays = {}
        
        no_atoms = (len(loop_sequence)+1) * 4 
        
        MECHANO_arrays["lengthsN"] = [0] * no_atoms
        MECHANO_arrays["anglesN"] = [0] * no_atoms
        MECHANO_arrays["dihedralsN"] = [0] * no_atoms
        
        MECHANO_arrays["lengthsC"] = [0] * no_atoms
        MECHANO_arrays["anglesC"] = [0] * no_atoms
        MECHANO_arrays["dihedralsC"] = [0] * no_atoms
                 
        #--------------------------------------------------------------------------------------------------
        #--------------------------------------------------------------------------------------------------
        #Insert values into arrays - N direction
        
        MECHANO_arrays["dihedralsN"][0] = (building_data["realN"].oxy + 180) * numpy.pi/180
        MECHANO_arrays["dihedralsN"][1] = building_data[1].omega * numpy.pi/180
        
        keys = list(alignment.keys())
        for target_index in keys:
            if alignment[target_index] != None:
                match = alignment[target_index]
                arr_index = target_index*4
                
                MECHANO_arrays["lengthsN"][arr_index]   = building_data[match].length1
                MECHANO_arrays["lengthsN"][arr_index+1] = building_data[match].length2
                MECHANO_arrays["lengthsN"][arr_index+2] = building_data[match].length3
                MECHANO_arrays["lengthsN"][arr_index+3] = building_data[match].length4
                
                MECHANO_arrays["anglesN"][arr_index]   = building_data[match].angle1
                MECHANO_arrays["anglesN"][arr_index+1] = building_data[match].angle2
                MECHANO_arrays["anglesN"][arr_index+2] = building_data[match].angle3
                MECHANO_arrays["anglesN"][arr_index+3] = building_data[match].angle4
                
                MECHANO_arrays["dihedralsN"][arr_index+2] = building_data[match].phi * numpy.pi/180
                MECHANO_arrays["dihedralsN"][arr_index+3] = building_data[match].oxy * numpy.pi/180
                MECHANO_arrays["dihedralsN"][arr_index+4] = building_data[match].psi * numpy.pi/180
                MECHANO_arrays["dihedralsN"][arr_index+5] = building_data[match].omega * numpy.pi/180
                   
        MECHANO_arrays["lengthsN"][no_atoms-4] = building_data[len(decoy_residues)-2].length1
        MECHANO_arrays["lengthsN"][no_atoms-3] = building_data["realC"].length2
        MECHANO_arrays["lengthsN"][no_atoms-2] = building_data["realC"].length3
        MECHANO_arrays["lengthsN"][no_atoms-1] = building_data["realC"].length4
        
        MECHANO_arrays["anglesN"][no_atoms-4] = building_data[len(decoy_residues)-2].angle1
        MECHANO_arrays["anglesN"][no_atoms-3] = building_data[len(decoy_residues)-2].angle2
        MECHANO_arrays["anglesN"][no_atoms-2] = building_data["realC"].angle3
        MECHANO_arrays["anglesN"][no_atoms-1] = building_data["realC"].angle4
        
        MECHANO_arrays["dihedralsN"][no_atoms-2] = building_data[len(decoy_residues)-2].phi * numpy.pi/180
        MECHANO_arrays["dihedralsN"][no_atoms-1] = building_data["realC"].oxy * numpy.pi/180
        
        
        
        C_lengthsN = (ctypes.c_double * len(MECHANO_arrays["lengthsN"]))(*MECHANO_arrays["lengthsN"])
        C_anglesN = (ctypes.c_double * len(MECHANO_arrays["anglesN"]))(*MECHANO_arrays["anglesN"])
        C_dihedralsN = (ctypes.c_double * len(MECHANO_arrays["dihedralsN"]))(*MECHANO_arrays["dihedralsN"])
        
        #--------------------------------------------------------------------------------------------------
        #--------------------------------------------------------------------------------------------------
        #Insert values into arrays - C direction
        keys.reverse()
        MECHANO_arrays["lengthsC"][0]   = building_data[len(decoy_residues)-2].length1
        
        MECHANO_arrays["anglesC"][0]    = building_data[len(decoy_residues)-2].angle2
        MECHANO_arrays["anglesC"][2]    = building_data[len(decoy_residues)-2].angle1
        
        MECHANO_arrays["dihedralsC"][0] = building_data[len(decoy_residues)-2].phi * numpy.pi/180
        
        
        for target_index in keys:
            if alignment[target_index] != None:
                match = alignment[target_index]
                arr_index = no_atoms - (4*target_index) - 7
                
                MECHANO_arrays["lengthsC"][arr_index]   = building_data[match].length4
                MECHANO_arrays["lengthsC"][arr_index+1] = building_data[match].length3
                MECHANO_arrays["lengthsC"][arr_index+2] = building_data[match].length2
                MECHANO_arrays["lengthsC"][arr_index+3] = building_data[match].length1
                
                MECHANO_arrays["anglesC"][arr_index]   = building_data[match].angle5
                MECHANO_arrays["anglesC"][arr_index+2] = building_data[match].angle3
                MECHANO_arrays["anglesC"][arr_index+3] = building_data[match].angle2
                MECHANO_arrays["anglesC"][arr_index+5] = building_data[match].angle1
                
                MECHANO_arrays["dihedralsC"][arr_index]   = building_data[match].oxy2 * numpy.pi/180
                MECHANO_arrays["dihedralsC"][arr_index+1] = building_data[match].omega * numpy.pi/180
                MECHANO_arrays["dihedralsC"][arr_index+2] = building_data[match].psi * numpy.pi/180
                MECHANO_arrays["dihedralsC"][arr_index+3] = building_data[match].phi * numpy.pi/180
                
        
        MECHANO_arrays["lengthsC"][no_atoms-3] = building_data["realN"].length4
        MECHANO_arrays["lengthsC"][no_atoms-2] = building_data["realN"].length3
        MECHANO_arrays["lengthsC"][no_atoms-1] = building_data["realN"].length2
        
        MECHANO_arrays["anglesC"][no_atoms-3] = building_data[1].angle5
        MECHANO_arrays["anglesC"][no_atoms-1] = building_data["realN"].angle3
        
        MECHANO_arrays["dihedralsC"][no_atoms-3] = building_data[1].oxy2 * numpy.pi/180
        MECHANO_arrays["dihedralsC"][no_atoms-2] = building_data[1].omega * numpy.pi/180
        MECHANO_arrays["dihedralsC"][no_atoms-1] = building_data[1].psi * numpy.pi/180
                   
        C_lengthsC = (ctypes.c_double * len(MECHANO_arrays["lengthsC"]))(*MECHANO_arrays["lengthsC"])
        C_anglesC = (ctypes.c_double * len(MECHANO_arrays["anglesC"]))(*MECHANO_arrays["anglesC"])
        C_dihedralsC = (ctypes.c_double * len(MECHANO_arrays["dihedralsC"]))(*MECHANO_arrays["dihedralsC"])
            
        #--------------------------------------------------------------------------------------------------
        #--------------------------------------------------------------------------------------------------
        #Run Stuff!!!
        
        if "0" in bitstring:
            pass
        else:
            no_decoys = 1                                                                    
                              
        #Write some stuff to the log file here!
        log_file = open(results_file, "a")
        
        log_file.write("\n>Fragment ID: %s_%d_%d\n"
                       "  Alignment:\n"
                       "      %s\n"
                       "      %s\n"
                       "  ESSS = %d\n"  %(decoy.struc, decoy.start, len(decoy.seq), decoy.aligned1, decoy.aligned2, decoy.score))
        
        log_file.close()
        
        decoy_code = str("%s_%d_%d" %(decoy.struc, decoy.start, len(decoy.seq)))
  
        #CALL MECHANO!! 
        print("Fragment %d/%d: %s" %(count, len(chosen_decoys), decoy_code))
        print("    Target  : %s" %decoy.aligned1)
        print("    Fragment: %s" %decoy.aligned2)
        MECHANO.main(pdb_file.encode('utf-8'), chain, loop_sequence.encode('utf-8'), len(loop_sequence), int(MECHANOstart)+1, 
                     int(MECHANOend), no_decoys, C_lengthsN, C_anglesN, C_dihedralsN, C_lengthsC, 
                     C_anglesC, C_dihedralsC, bitstring, results_file.encode('utf-8'), decoy_file.encode('utf-8'), C_native_coords, 
                     C_Nanchor_coords, C_Canchor_coords, C_Ncoords, C_CAcoords, C_Ccoords, C_Ocoords, 
                     C_atom_freqs, C_info_array, N_seq.encode('utf-8'), C_seq.encode('utf-8'), decoy_code.encode('utf-8'), sphinx_version.encode('utf-8'))
     
  
  
  def clear_results(self):
    self.results = []
    self.warnings = []
    self.counter = itertools.count(1)
    self.idecoy = 0


  def open_errstream(self, filename):
      if filename == "-":
        self.errstream = sys.stderr
      else:
        self.errstream = open(filename, "w")


  def close_errstream(self):
      if self.errstream != sys.stderr:
        self.errstream.close()
        self.set_silent()
  
  
  def automodel_loop(self, start_of_loop, end_of_loop, loop_sequence, pdb_file, results_file, 
                    decoy_directory,
                    scriptpath,
                    no_fragments,
                    dbdir = '.',
                    resnum = False,
                    chain = None,
                    ignore_list = [],
                    sphinx_version = "general",
                    ranking_method = "soaploop",
                    sidechain_method = "pears",
                    **kwargs):
      """Search the given database and return a list of decoys for the given loop.
      """
      
      # Set options defined in the Fread constructor
      self.set_options(kwargs)
      
      for db_group in dbdir.split(":"): # try these databases until we find something
        databases = db_group.split("|") # combine results from all these databases
        self.set_db(databases[0])
        
        try:
          self.set_loop(start_of_loop, loop_sequence.upper(), resnum=resnum, chain=chain)
        except FreadError:
          break
        
        # Try to model, and extend loop if nothing was found
        anyresults = False
        for db in databases:
          self.set_db(db)
          try:
            anyresults |= bool(self.model(pdb_file, start_of_loop, end_of_loop, results_file, decoy_directory, chain, scriptpath, no_fragments, ignore_list, sphinx_version, ranking_method=ranking_method, sidechain_method=sidechain_method))
                    
          except FreadError:
            break
              
      return self.results
