scoring.py

Go to the documentation of this file.

import os
from ost import mol
from ost import seq
from ost import io
from ost import conop
from ost import settings
from ost import geom
from ost import LogScript, LogInfo, LogDebug
from ost.mol.alg import lddt
from ost.mol.alg import qsscore
from ost.mol.alg import chain_mapping
from ost.mol.alg import stereochemistry
from ost.mol.alg import dockq
from ost.mol.alg.lddt import lDDTScorer
from ost.mol.alg.qsscore import QSScorer
from ost.mol.alg.contact_score import ContactScorer
from ost.mol.alg.contact_score import ContactEntity
from ost.mol.alg import GDT
from ost.mol.alg import Molck, MolckSettings
from ost import bindings
from ost.bindings import cadscore
from ost.bindings import tmtools
import numpy as np
 
class lDDTBSScorer:
    """Scorer specific for a reference/model pair
 
    Finds best possible binding site representation of reference in model given
    lDDT score. Uses :class:`ost.mol.alg.chain_mapping.ChainMapper` to deal with
    chain mapping.
 
    :param reference: Reference structure
    :type reference: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
    :param model: Model structure
    :type model: :class:`ost.mol.EntityView`/:class:`ost.mol.EntityHandle`
    :param residue_number_alignment: Passed to ChainMapper constructor
    :type residue_number_alignment: :class:`bool`
    """
    def __init__(self, reference, model,
                 residue_number_alignment=False):
        self.chain_mapperchain_mapper = chain_mapping.ChainMapper(reference,
            resnum_alignments=residue_number_alignment)
        self.refref = self.chain_mapperchain_mapper.target
        self.mdlmdl = model
 
    def ScoreBS(self, ligand, radius = 4.0, lddt_radius=10.0):
        """Computes binding site lDDT score given *ligand*. Best possible
        binding site representation is selected by lDDT but other scores such as
        CA based RMSD and GDT are computed too and returned.
 
        :param ligand: Defines the scored binding site, i.e. provides positions
                       to perform proximity search
        :type ligand: r'((Residue)|(Chain)|(Entity))((View)|(Handle))'
        :param radius: Reference residues with any atom position within *radius*
                       of *ligand* consitute the scored binding site
        :type radius: :class:`float`
        :param lddt_radius: Passed as *inclusion_radius* to
                            :class:`ost.mol.alg.lddt.lDDTScorer`
        :type lddt_radius: :class:`float`
        :returns: Object of type :class:`ost.mol.alg.chain_mapping.ReprResult`
                  containing all atom lDDT score and mapping information.
                  None if no representation could be found.
        """
 
        # create view of reference binding site
        ref_residues_hashes = set() # helper to keep track of added residues
        for ligand_at in ligand.atoms:
            close_atoms = self.refref.FindWithin(ligand_at.GetPos(), radius)
            for close_at in close_atoms:
                ref_res = close_at.GetResidue()
                h = ref_res.handle.GetHashCode()
                if h not in ref_residues_hashes:
                    ref_residues_hashes.add(h)
 
        # reason for doing that separately is to guarantee same ordering of
        # residues as in underlying entity. (Reorder by ResNum seems only
        # available on ChainHandles)
        ref_bs = self.refref.CreateEmptyView()
        for ch in self.refref.chains:
            for r in ch.residues:
                if r.handle.GetHashCode() in ref_residues_hashes:
                    ref_bs.AddResidue(r, mol.ViewAddFlag.INCLUDE_ALL)
 
        # gogogo
        bs_repr = self.chain_mapperchain_mapper.GetRepr(ref_bs, self.mdlmdl,
                                            inclusion_radius = lddt_radius)
        if len(bs_repr) >= 1:
            return bs_repr[0]
        else:
            return None
 
 
class Scorer:
    """ Helper class to access the various scores available from ost.mol.alg
 
    Deals with structure cleanup, chain mapping, interface identification etc.
    Intermediate results are available as attributes.
 
    :param model: Model structure - a deep copy is available as :attr:`model`.
                  Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
                  is applied.
    :type model: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
    :param target: Target structure - a deep copy is available as :attr:`target`.
                  Additionally, :func:`ost.mol.alg.Molck` using *molck_settings*
                  is applied.
    :type target: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
    :param resnum_alignments: Whether alignments between chemically equivalent
                              chains in *model* and *target* can be computed
                              based on residue numbers. This can be assumed in
                              benchmarking setups such as CAMEO/CASP.
    :type resnum_alignments: :class:`bool`
    :param molck_settings: Settings used for Molck on *model* and *target*, if
                           set to None, a default object is constructed by
                           setting everything except rm_zero_occ_atoms and
                           colored to True in
                           :class:`ost.mol.alg.MolckSettings` constructor.
    :type molck_settings: :class:`ost.mol.alg.MolckSettings`
    :param cad_score_exec: Explicit path to voronota-cadscore executable from
                           voronota installation from 
                           https://github.com/kliment-olechnovic/voronota. If
                           not given, voronota-cadscore must be in PATH if any
                           of the CAD score related attributes is requested.
    :type cad_score_exec: :class:`str`
    :param custom_mapping: Provide custom chain mapping between *model* and
                           *target*. Dictionary with target chain names as key
                           and model chain names as value.
                           :attr:`~mapping` is constructed from this.
    :type custom_mapping: :class:`dict`
    :param custom_rigid_mapping: Provide custom chain mapping between *model*
                                 and *target*. Dictionary with target chain
                                 names as key and model chain names as value.
                                 :attr:`~rigid_mapping` is constructed from
                                 this.
    :type custom_rigid_mapping: :class:`dict`
    :param usalign_exec: Explicit path to USalign executable used to compute
                         TM-score. If not given, TM-score will be computed
                         with OpenStructure internal copy of USalign code.
    :type usalign_exec: :class:`str`
    :param lddt_no_stereochecks: Whether to compute lDDT without stereochemistry
                                checks
    :type lddt_no_stereochecks: :class:`bool`
    :param n_max_naive: Parameter for chain mapping. If the number of possible
                        mappings is <= *n_max_naive*, the full
                        mapping solution space is enumerated to find the
                        the optimum. A heuristic is used otherwise. The default
                        of 40320 corresponds to an octamer (8! = 40320).
                        A structure with stoichiometry A6B2 would be
                        6!*2! = 1440 etc.
    :type n_max_naive: :class:`int`
    :param oum: Override USalign Mapping. Inject rigid_mapping of
                :class:`Scorer` object into USalign to compute TM-score.
                Experimental feature with limitations.
    :type oum: :class:`bool`
    :param min_pep_length: Relevant parameter if short peptides are involved in
                           scoring. Minimum peptide length for a chain in the
                           target structure to be considered in chain mapping.
                           The chain mapping algorithm first performs an all vs.
                           all pairwise sequence alignment to identify \"equal\"
                           chains within the target structure. We go for simple
                           sequence identity there. Short sequences can be
                           problematic as they may produce high sequence identity
                           alignments by pure chance.
    :type min_pep_length: :class:`int`
    :param min_nuc_length: Relevant parameter if short nucleotides are involved
                           in scoring. Minimum nucleotide length for a chain in
                           the target structure to be considered in chain
                           mapping. The chain mapping algorithm first performs
                           an all vs. all pairwise sequence alignment to
                           identify \"equal\" chains within the target
                           structure. We go for simple sequence identity there.
                           Short sequences can be problematic as they may
                           produce high sequence identity alignments by pure
                           chance.
    :type min_nuc_length: :class:`int`
    :param lddt_add_mdl_contacts: lDDT specific flag. Only using contacts in
                                  lDDT that are within a certain distance 
                                  threshold in the target does not penalize
                                  for added model contacts. If set to True, this
                                  flag will also consider target contacts
                                  that are within the specified distance
                                  threshold in the model but not necessarily in
                                  the target. No contact will be added if the
                                  respective atom pair is not resolved in the
                                  target.
    :type lddt_add_mdl_contacts: :class:`bool`
    :param dockq_capri_peptide: Flag that changes two things in the way DockQ
                                and its underlying scores are computed which is
                                proposed by the CAPRI community when scoring
                                peptides (PMID: 31886916).
                                ONE: Two residues are considered in contact if 
                                any of their atoms is within 5A. This is
                                relevant for fnat and fnonat scores. CAPRI
                                suggests to lower this threshold to 4A for
                                protein-peptide interactions.
                                TWO: irmsd is computed on interface residues. A
                                residue is defined as interface residue if any
                                of its atoms is within 10A of another chain.
                                CAPRI suggests to lower the default of 10A to
                                8A in combination with only considering CB atoms
                                for protein-peptide interactions.
                                This flag has no influence on patch_dockq
                                scores.
    :type dockq_capri_peptide: :class:`bool`
    """
    def __init__(self, model, target, resnum_alignments=False,
                 molck_settings = None, cad_score_exec = None,
                 custom_mapping=None, custom_rigid_mapping=None,
                 usalign_exec = None, lddt_no_stereochecks=False,
                 n_max_naive=40320, oum=False, min_pep_length = 6,
                 min_nuc_length = 4, lddt_add_mdl_contacts=False,
                 dockq_capri_peptide=False):
 
        self._target_orig_target_orig = target
        self._model_orig_model_orig = model
 
        if isinstance(self._model_orig_model_orig, mol.EntityView):
            self._model_model = mol.CreateEntityFromView(self._model_orig_model_orig, False)
        else:
            self._model_model = self._model_orig_model_orig.Copy()
 
        if isinstance(self._target_orig_target_orig, mol.EntityView):
            self._target_target = mol.CreateEntityFromView(self._target_orig_target_orig, False)
        else:
            self._target_target = self._target_orig_target_orig.Copy()
 
        if molck_settings is None:
            molck_settings = MolckSettings(rm_unk_atoms=True,
                                           rm_non_std=False,
                                           rm_hyd_atoms=True,
                                           rm_oxt_atoms=True,
                                           rm_zero_occ_atoms=False,
                                           colored=False,
                                           map_nonstd_res=True,
                                           assign_elem=True)
        LogScript("Cleaning up input structures")
        Molck(self._model_model, conop.GetDefaultLib(), molck_settings)
        Molck(self._target_target, conop.GetDefaultLib(), molck_settings)
        self._model_model = self._model_model.Select("peptide=True or nucleotide=True")
        self._target_target = self._target_target.Select("peptide=True or nucleotide=True")
 
        # catch models which have empty chain names
        for ch in self._model_model.chains:
            if ch.GetName().strip() == "":
                raise RuntimeError("Model chains must have valid chain names")
            if ch.GetName().strip() == "'" or ch.GetName().strip() == '"':
                raise RuntimeError("Model chains cannot be named with quote signs (' or \"\")")
        
        # catch targets which have empty chain names
        for ch in self._target_target.chains:
            if ch.GetName().strip() == "":
                raise RuntimeError("Target chains must have valid chain names")
            if ch.GetName().strip() == "'" or ch.GetName().strip() == '"':
                raise RuntimeError("Target chains cannot be named with quote signs (' or \"\")")
 
        if resnum_alignments:
            # In case of resnum_alignments, we have some requirements on 
            # residue numbers in the chain mapping: 1) no ins codes 2) strictly
            # increasing residue numbers.
            for ch in self._model_model.chains:
                ins_codes = [r.GetNumber().GetInsCode() for r in ch.residues]
                if len(set(ins_codes)) != 1 or ins_codes[0] != '\0':
                    raise RuntimeError("Residue numbers in each model chain "
                                       "must not contain insertion codes if "
                                       "resnum_alignments are enabled")
                nums = [r.GetNumber().GetNum() for r in ch.residues]
                if not all(i < j for i, j in zip(nums, nums[1:])):
                    raise RuntimeError("Residue numbers in each model chain "
                                       "must be strictly increasing if "
                                       "resnum_alignments are enabled")
 
            for ch in self._target_target.chains:
                ins_codes = [r.GetNumber().GetInsCode() for r in ch.residues]
                if len(set(ins_codes)) != 1 or ins_codes[0] != '\0':
                    raise RuntimeError("Residue numbers in each target chain "
                                       "must not contain insertion codes if "
                                       "resnum_alignments are enabled")
                nums = [r.GetNumber().GetNum() for r in ch.residues]
                if not all(i < j for i, j in zip(nums, nums[1:])):
                    raise RuntimeError("Residue numbers in each target chain "
                                       "must be strictly increasing if "
                                       "resnum_alignments are enabled")
 
        if usalign_exec is not None:
            if not os.path.exists(usalign_exec):
                raise RuntimeError(f"USalign exec ({usalign_exec}) "
                                   f"not found")
            if not os.access(usalign_exec, os.X_OK):
                raise RuntimeError(f"USalign exec ({usalign_exec}) "
                                   f"is not executable")
 
        self.resnum_alignmentsresnum_alignments = resnum_alignments
        self.cad_score_execcad_score_exec = cad_score_exec
        self.usalign_execusalign_exec = usalign_exec
        self.lddt_no_stereocheckslddt_no_stereochecks = lddt_no_stereochecks
        self.n_max_naiven_max_naive = n_max_naive
        self.oumoum = oum
        self.min_pep_lengthmin_pep_length = min_pep_length
        self.min_nuc_lengthmin_nuc_length = min_nuc_length
        self.lddt_add_mdl_contactslddt_add_mdl_contacts = lddt_add_mdl_contacts
        self.dockq_capri_peptidedockq_capri_peptide = dockq_capri_peptide
 
        # lazily evaluated attributes
        self._stereochecked_model_stereochecked_model = None
        self._stereochecked_target_stereochecked_target = None
        self._model_clashes_model_clashes = None
        self._model_bad_bonds_model_bad_bonds = None
        self._model_bad_angles_model_bad_angles = None
        self._target_clashes_target_clashes = None
        self._target_bad_bonds_target_bad_bonds = None
        self._target_bad_angles_target_bad_angles = None
        self._chain_mapper_chain_mapper = None
        self._mapping_mapping = None
        self._rigid_mapping_rigid_mapping = None
        self._model_interface_residues_model_interface_residues = None
        self._target_interface_residues_target_interface_residues = None
        self._aln_aln = None
        self._stereochecked_aln_stereochecked_aln = None
        self._pepnuc_aln_pepnuc_aln = None
 
        # lazily constructed scorer objects
        self._lddt_scorer_lddt_scorer = None
        self._bb_lddt_scorer_bb_lddt_scorer = None
        self._qs_scorer_qs_scorer = None
        self._contact_scorer_contact_scorer = None
 
        # lazily computed scores
        self._lddt_lddt = None
        self._local_lddt_local_lddt = None
        self._bb_lddt_bb_lddt = None
        self._bb_local_lddt_bb_local_lddt = None
        self._ilddt_ilddt = None
 
        self._qs_global_qs_global = None
        self._qs_best_qs_best = None
        self._qs_target_interfaces_qs_target_interfaces = None
        self._qs_model_interfaces_qs_model_interfaces = None
        self._qs_interfaces_qs_interfaces = None
        self._per_interface_qs_global_per_interface_qs_global = None
        self._per_interface_qs_best_per_interface_qs_best = None
 
        self._contact_target_interfaces_contact_target_interfaces = None
        self._contact_model_interfaces_contact_model_interfaces = None
        self._native_contacts_native_contacts = None
        self._model_contacts_model_contacts = None
        self._ics_precision_ics_precision = None
        self._ics_recall_ics_recall = None
        self._ics_ics = None
        self._per_interface_ics_precision_per_interface_ics_precision = None
        self._per_interface_ics_recall_per_interface_ics_recall = None
        self._per_interface_ics_per_interface_ics = None
        self._ips_precision_ips_precision = None
        self._ips_recall_ips_recall = None
        self._ips_ips = None
        self._per_interface_ics_precision_per_interface_ics_precision = None
        self._per_interface_ics_recall_per_interface_ics_recall = None
        self._per_interface_ics_per_interface_ics = None
 
        self._dockq_target_interfaces_dockq_target_interfaces = None
        self._dockq_interfaces_dockq_interfaces = None
        self._fnat_fnat = None
        self._fnonnat_fnonnat = None
        self._irmsd_irmsd = None
        self._lrmsd_lrmsd = None
        self._nnat_nnat = None
        self._nmdl_nmdl = None
        self._dockq_scores_dockq_scores = None
        self._dockq_ave_dockq_ave = None
        self._dockq_wave_dockq_wave = None
        self._dockq_ave_full_dockq_ave_full = None
        self._dockq_wave_full_dockq_wave_full = None
 
        self._mapped_target_pos_mapped_target_pos = None
        self._mapped_model_pos_mapped_model_pos = None
        self._transformed_mapped_model_pos_transformed_mapped_model_pos = None
        self._n_target_not_mapped_n_target_not_mapped = None
        self._transform_transform = None
 
        self._rigid_mapped_target_pos_rigid_mapped_target_pos = None
        self._rigid_mapped_model_pos_rigid_mapped_model_pos = None
        self._rigid_transformed_mapped_model_pos_rigid_transformed_mapped_model_pos = None
        self._rigid_n_target_not_mapped_rigid_n_target_not_mapped = None
        self._rigid_transform_rigid_transform = None
 
        self._gdt_05_gdt_05 = None
        self._gdt_1_gdt_1 = None
        self._gdt_2_gdt_2 = None
        self._gdt_4_gdt_4 = None
        self._gdt_8_gdt_8 = None
        self._gdtts_gdtts = None
        self._gdtha_gdtha = None
        self._rmsd_rmsd = None
 
        self._cad_score_cad_score = None
        self._local_cad_score_local_cad_score = None
 
        self._patch_qs_patch_qs = None
        self._patch_dockq_patch_dockq = None
 
        self._tm_score_tm_score = None
        self._usalign_mapping_usalign_mapping = None
 
        if custom_mapping is not None:
            self._mapping_mapping = self._construct_custom_mapping_construct_custom_mapping(custom_mapping)
 
        if custom_rigid_mapping is not None:
            self._rigid_mapping_rigid_mapping = \
            self._construct_custom_mapping_construct_custom_mapping(custom_rigid_mapping)
        LogDebug("Scorer sucessfully initialized")
 
    @property
    def model(self):
        """ Model with Molck cleanup
 
        :type: :class:`ost.mol.EntityHandle`
        """
        return self._model_model
 
    @property
    def model_orig(self):
        """ The original model passed at object construction
 
        :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
        """
        return self._model_orig_model_orig
 
    @property
    def target(self):
        """ Target with Molck cleanup
 
        :type: :class:`ost.mol.EntityHandle`
        """
        return self._target_target
 
    @property
    def target_orig(self):
        """ The original target passed at object construction
 
        :type: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
        """
        return self._target_orig_target_orig
 
    @property
    def aln(self):
        """ Alignments of :attr:`model`/:attr:`target` chains
 
        Alignments for each pair of chains mapped in :attr:`mapping`.
        First sequence is target sequence, second sequence the model sequence.
 
        :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
        """
        if self._aln_aln is None:
            self._compute_aln_compute_aln()
        return self._aln_aln
 
    @property
    def stereochecked_aln(self):
        """ Stereochecked equivalent of :attr:`aln`
 
        The alignments may differ, as stereochecks potentially remove residues
 
        :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
        """
        if self._stereochecked_aln_stereochecked_aln is None:
            self._compute_stereochecked_aln_compute_stereochecked_aln()
        return self._stereochecked_aln_stereochecked_aln
 
    @property
    def pepnuc_aln(self):
        """ Alignments of :attr:`model_orig`/:attr:`target_orig` chains
 
        Selects for peptide and nucleotide residues before sequence
        extraction. Includes residues that would be removed by molck in
        structure preprocessing.
 
        :type: :class:`list` of :class:`ost.seq.AlignmentHandle`
        """
        if self._pepnuc_aln_pepnuc_aln is None:
            self._compute_pepnuc_aln_compute_pepnuc_aln()
        return self._pepnuc_aln_pepnuc_aln
 
    @property
    def stereochecked_model(self):
        """ View of :attr:`~model` that has stereochemistry checks applied
 
        First, a selection for peptide/nucleotide residues is performed,
        secondly peptide sidechains with stereochemical irregularities are
        removed (full residue if backbone atoms are involved). Irregularities
        are clashes or bond lengths/angles more than 12 standard deviations
        from expected values.
 
        :type: :class:`ost.mol.EntityView`
        """
        if self._stereochecked_model_stereochecked_model is None:
            self._do_stereochecks_do_stereochecks()
        return self._stereochecked_model_stereochecked_model
 
    @property
    def model_clashes(self):
        """ Clashing model atoms
 
        :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
        """
        if self._model_clashes_model_clashes is None:
            self._do_stereochecks_do_stereochecks()
        return self._model_clashes_model_clashes
 
    @property
    def model_bad_bonds(self):
        """ Model bonds with unexpected stereochemistry
 
        :type: :class:`list` of
               :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
        """
        if self._model_bad_bonds_model_bad_bonds is None:
            self._do_stereochecks_do_stereochecks()
        return self._model_bad_bonds_model_bad_bonds
 
    @property
    def model_bad_angles(self):
        """ Model angles with unexpected stereochemistry
 
        :type: :class:`list` of
               :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
        """
        if self._model_bad_angles_model_bad_angles is None:
            self._do_stereochecks_do_stereochecks()
        return self._model_bad_angles_model_bad_angles
 
    @property
    def stereochecked_target(self):
        """ Same as :attr:`~stereochecked_model` for :attr:`~target`
 
        :type: :class:`ost.mol.EntityView`
        """
        if self._stereochecked_target_stereochecked_target is None:
            self._do_stereochecks_do_stereochecks()
        return self._stereochecked_target_stereochecked_target
 
    @property
    def target_clashes(self):
        """ Clashing target atoms
 
        :type: :class:`list` of :class:`ost.mol.alg.stereochemistry.ClashInfo`
        """
        if self._target_clashes_target_clashes is None:
            self._do_stereochecks_do_stereochecks()
        return self._target_clashes_target_clashes
 
    @property
    def target_bad_bonds(self):
        """ Target bonds with unexpected stereochemistry
 
        :type: :class:`list` of
               :class:`ost.mol.alg.stereochemistry.BondViolationInfo`
        """
        if self._target_bad_bonds_target_bad_bonds is None:
            self._do_stereochecks_do_stereochecks()
        return self._target_bad_bonds_target_bad_bonds
 
    @property
    def target_bad_angles(self):
        """ Target angles with unexpected stereochemistry
 
        :type: :class:`list` of
               :class:`ost.mol.alg.stereochemistry.AngleViolationInfo`
        """
        if self._target_bad_angles_target_bad_angles is None:
            self._do_stereochecks_do_stereochecks()
        return self._target_bad_angles_target_bad_angles
 
    @property
    def chain_mapper(self):
        """ Chain mapper object for given :attr:`target`
 
        :type: :class:`ost.mol.alg.chain_mapping.ChainMapper`
        """
        if self._chain_mapper_chain_mapper is None:
            self._chain_mapper_chain_mapper = chain_mapping.ChainMapper(self.targettarget,
                                                           n_max_naive=1e9,
                                                           resnum_alignments=self.resnum_alignmentsresnum_alignments,
                                                           min_pep_length=self.min_pep_lengthmin_pep_length,
                                                           min_nuc_length=self.min_nuc_lengthmin_nuc_length)
        return self._chain_mapper_chain_mapper
 
    @property
    def mapping(self):
        """ Full chain mapping result for :attr:`target`/:attr:`model`
 
        Computed with :func:`ost.mol.alg.ChainMapper.GetMapping`
 
        :type: :class:`ost.mol.alg.chain_mapping.MappingResult` 
        """
        if self._mapping_mapping is None:
            LogScript("Computing chain mapping")
            self._mapping_mapping = \
            self.chain_mapperchain_mapper.GetMapping(self.modelmodel,
                                         n_max_naive = self.n_max_naiven_max_naive)
        return self._mapping_mapping
 
    @property
    def rigid_mapping(self):
        """ Full chain mapping result for :attr:`target`/:attr:`model`
 
        Computed with :func:`ost.mol.alg.ChainMapper.GetRMSDMapping`
 
        :type: :class:`ost.mol.alg.chain_mapping.MappingResult` 
        """
        if self._rigid_mapping_rigid_mapping is None:
            LogScript("Computing rigid chain mapping")
            self._rigid_mapping_rigid_mapping = \
            self.chain_mapperchain_mapper.GetRMSDMapping(self.modelmodel)
        return self._rigid_mapping_rigid_mapping
 
    @property
    def model_interface_residues(self):
        """ Interface residues in :attr:`~model`
 
        Thats all residues having a contact with at least one residue from
        another chain (CB-CB distance <= 8A, CA in case of Glycine)
 
        :type: :class:`dict` with chain names as key and and :class:`list`
                with residue numbers of the respective interface residues.
        """
        if self._model_interface_residues_model_interface_residues is None:
            self._model_interface_residues_model_interface_residues = \
            self._get_interface_residues_get_interface_residues(self.modelmodel)
        return self._model_interface_residues_model_interface_residues
 
    @property
    def target_interface_residues(self):
        """ Same as :attr:`~model_interface_residues` for :attr:`~target`
 
        :type: :class:`dict` with chain names as key and and :class:`list`
                with residue numbers of the respective interface residues.
        """
        if self._target_interface_residues_target_interface_residues is None:
            self._target_interface_residues_target_interface_residues = \
            self._get_interface_residues_get_interface_residues(self.targettarget)
        return self._target_interface_residues_target_interface_residues
 
    @property
    def lddt_scorer(self):
        """ lDDT scorer for :attr:`~stereochecked_target` (default parameters)
 
        :type: :class:`ost.mol.alg.lddt.lDDTScorer`
        """
        if self._lddt_scorer_lddt_scorer is None:
            if self.lddt_no_stereocheckslddt_no_stereochecks:
                self._lddt_scorer_lddt_scorer = lDDTScorer(self.targettarget)
            else:
                self._lddt_scorer_lddt_scorer = lDDTScorer(self.stereochecked_targetstereochecked_target)
        return self._lddt_scorer_lddt_scorer
 
    @property
    def bb_lddt_scorer(self):
        """ Backbone only lDDT scorer for :attr:`~target`
 
        No stereochecks applied for bb only lDDT which considers CA atoms
        for peptides and C3' atoms for nucleotides.
 
        :type: :class:`ost.mol.alg.lddt.lDDTScorer`
        """
        if self._bb_lddt_scorer_bb_lddt_scorer is None:
            self._bb_lddt_scorer_bb_lddt_scorer = lDDTScorer(self.targettarget, bb_only=True)
        return self._bb_lddt_scorer_bb_lddt_scorer
 
    @property
    def qs_scorer(self):
        """ QS scorer constructed from :attr:`~mapping`
 
        The scorer object is constructed with default parameters and relates to
        :attr:`~model` and :attr:`~target` (no stereochecks).
 
        :type: :class:`ost.mol.alg.qsscore.QSScorer`
        """
        if self._qs_scorer_qs_scorer is None:
            self._qs_scorer_qs_scorer = QSScorer.FromMappingResult(self.mappingmapping)
        return self._qs_scorer_qs_scorer
 
    @property
    def contact_scorer(self):
        if self._contact_scorer_contact_scorer is None:
            self._contact_scorer_contact_scorer = ContactScorer.FromMappingResult(self.mappingmapping)
        return self._contact_scorer_contact_scorer
    
 
    @property
    def lddt(self):
        """ Global lDDT score in range [0.0, 1.0]
 
        Computed based on :attr:`~stereochecked_model`. In case of oligomers,
        :attr:`~mapping` is used.
 
        :type: :class:`float`
        """
        if self._lddt_lddt is None:
            self._compute_lddt_compute_lddt()
        return self._lddt_lddt
    
    @property
    def local_lddt(self):
        """ Per residue lDDT scores in range [0.0, 1.0]
 
        Computed based on :attr:`~stereochecked_model` but scores for all 
        residues in :attr:`~model` are reported. If a residue has been removed
        by stereochemistry checks, the respective score is set to 0.0. If a
        residue is not covered by the target or is in a chain skipped by the
        chain mapping procedure (happens for super short chains), the respective
        score is set to None. In case of oligomers, :attr:`~mapping` is used.
 
        :type: :class:`dict`
        """
        if self._local_lddt_local_lddt is None:
            self._compute_lddt_compute_lddt()
        return self._local_lddt_local_lddt
 
    @property
    def bb_lddt(self):
        """ Backbone only global lDDT score in range [0.0, 1.0]
 
        Computed based on :attr:`~model` on backbone atoms only. This is CA for
        peptides and C3' for nucleotides. No stereochecks are performed. In case
        of oligomers, :attr:`~mapping` is used.
 
        :type: :class:`float`
        """
        if self._bb_lddt_bb_lddt is None:
            self._compute_bb_lddt_compute_bb_lddt()
        return self._bb_lddt_bb_lddt
    
    @property
    def bb_local_lddt(self):
        """ Backbone only per residue lDDT scores in range [0.0, 1.0]
 
        Computed based on :attr:`~model` on backbone atoms only. This is CA for
        peptides and C3' for nucleotides. No stereochecks are performed. If a
        residue is not covered by the target or is in a chain skipped by the
        chain mapping procedure (happens for super short chains), the respective
        score is set to None. In case of oligomers, :attr:`~mapping` is used.
 
        :type: :class:`dict`
        """
        if self._bb_local_lddt_bb_local_lddt is None:
            self._compute_bb_lddt_compute_bb_lddt()
        return self._bb_local_lddt_bb_local_lddt
 
    @property
    def ilddt(self):
        """ Global interface lDDT score in range [0.0, 1.0]
 
        This is lDDT only based on inter-chain contacts. Value is None if no
        such contacts are present. For example if we're dealing with a monomer.
        Computed based on :attr:`~stereochecked_model` and :attr:`~mapping` for
        chain mapping.
 
        :type: :class:`float`
        """
        if self._ilddt_ilddt is None:
            # the whole None business kind of invalidates the idea of lazy
            # evaluation. The assumption is that this is called only once...
            self._compute_ilddt_compute_ilddt()
        return self._ilddt_ilddt
    
 
    @property
    def qs_global(self):
        """  Global QS-score
 
        Computed based on :attr:`model` using :attr:`mapping`
 
        :type: :class:`float`
        """
        if self._qs_global_qs_global is None:
            self._compute_qs_compute_qs()
        return self._qs_global_qs_global
 
    @property
    def qs_best(self):
        """  Global QS-score - only computed on aligned residues
 
        Computed based on :attr:`model` using :attr:`mapping`. The QS-score
        computation only considers contacts between residues with a mapping
        between target and model. As a result, the score won't be lowered in
        case of additional chains/residues in any of the structures.
 
        :type: :class:`float`
        """
        if self._qs_best_qs_best is None:
            self._compute_qs_compute_qs()
        return self._qs_best_qs_best
 
    @property
    def qs_target_interfaces(self):
        """ Interfaces in :attr:`~target` with non-zero contribution to
        :attr:`~qs_global`/:attr:`~qs_best`
 
        Chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 2 elements each:
               (trg_ch1, trg_ch2)
        """
        if self._qs_target_interfaces_qs_target_interfaces is None:
            self._qs_target_interfaces_qs_target_interfaces = self.qs_scorerqs_scorer.qsent1.interacting_chains
            self._qs_target_interfaces_qs_target_interfaces = \
            [(min(x[0],x[1]), max(x[0],x[1])) for x in self._qs_target_interfaces_qs_target_interfaces]
        return self._qs_target_interfaces_qs_target_interfaces
 
    @property
    def qs_model_interfaces(self):
        """ Interfaces in :attr:`~model` with non-zero contribution to
        :attr:`~qs_global`/:attr:`~qs_best`
 
        Chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 2 elements each:
               (mdl_ch1, mdl_ch2)
        """
        if self._qs_model_interfaces_qs_model_interfaces is None:
            self._qs_model_interfaces_qs_model_interfaces = self.qs_scorerqs_scorer.qsent2.interacting_chains
            self._qs_model_interfaces_qs_model_interfaces = \
            [(min(x[0],x[1]), max(x[0],x[1])) for x in self._qs_model_interfaces_qs_model_interfaces]
 
        return self._qs_model_interfaces_qs_model_interfaces
 
    @property
    def qs_interfaces(self):
        """ Interfaces in :attr:`~qs_target_interfaces` that can be mapped
        to :attr:`~model`.
 
        Target chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 4 elements each:
               (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
        """
        if self._qs_interfaces_qs_interfaces is None:
            self._qs_interfaces_qs_interfaces = list()
            flat_mapping = self.mappingmapping.GetFlatMapping()
            for i in self.qs_target_interfacesqs_target_interfaces:
                if i[0] in flat_mapping and i[1] in flat_mapping:
                    self._qs_interfaces_qs_interfaces.append((i[0], i[1],
                                                flat_mapping[i[0]],
                                                flat_mapping[i[1]]))
        return self._qs_interfaces_qs_interfaces
    
    @property
    def per_interface_qs_global(self):
        """ QS-score for each interface in :attr:`qs_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_qs_global_per_interface_qs_global is None:
            self._compute_per_interface_qs_scores_compute_per_interface_qs_scores()
        return self._per_interface_qs_global_per_interface_qs_global
    
    @property
    def per_interface_qs_best(self):
        """ QS-score for each interface in :attr:`qs_interfaces`
 
        Only computed on aligned residues
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_qs_best_per_interface_qs_best is None:
            self._compute_per_interface_qs_scores_compute_per_interface_qs_scores()
        return self._per_interface_qs_best_per_interface_qs_best
    
    @property
    def native_contacts(self):
        """ Native contacts
 
        A contact is a pair or residues from distinct chains that have
        a minimal heavy atom distance < 5A. Contacts are specified as
        :class:`tuple` with two strings in format:
        <cname>.<rnum>.<ins_code>
 
        :type: :class:`list` of :class:`tuple`
        """
        if self._native_contacts_native_contacts is None:
            self._native_contacts_native_contacts = self.contact_scorercontact_scorer.cent1.hr_contacts
        return self._native_contacts_native_contacts
 
    @property
    def model_contacts(self):
        """ Same for model
        """
        if self._model_contacts_model_contacts is None:
            self._model_contacts_model_contacts = self.contact_scorercontact_scorer.cent2.hr_contacts
        return self._model_contacts_model_contacts
 
    @property
    def contact_target_interfaces(self):
        """ Interfaces in :class:`target` which have at least one contact
 
        Contact as defined in :attr:`~native_contacts`,
        chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 2 elements each
               (trg_ch1, trg_ch2)
        """
        if self._contact_target_interfaces_contact_target_interfaces is None:
            tmp = self.contact_scorercontact_scorer.cent1.interacting_chains
            tmp = [(min(x[0],x[1]), max(x[0],x[1])) for x in tmp]
            self._contact_target_interfaces_contact_target_interfaces = tmp
        return self._contact_target_interfaces_contact_target_interfaces
 
    @property
    def contact_model_interfaces(self):
        """ Interfaces in :class:`model` which have at least one contact
 
        Contact as defined in :attr:`~native_contacts`,
        chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 2 elements each
               (mdl_ch1, mdl_ch2)
        """
        if self._contact_model_interfaces_contact_model_interfaces is None:
            tmp = self.contact_scorercontact_scorer.cent2.interacting_chains
            tmp = [(min(x[0],x[1]), max(x[0],x[1])) for x in tmp]
            self._contact_model_interfaces_contact_model_interfaces = tmp
        return self._contact_model_interfaces_contact_model_interfaces
 
    @property
    def ics_precision(self):
        """ Fraction of model contacts that are also present in target
 
        :type: :class:`float`
        """
        if self._ics_precision_ics_precision is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._ics_precision_ics_precision
    
    @property
    def ics_recall(self):
        """ Fraction of target contacts that are correctly reproduced in model
 
        :type: :class:`float`
        """
        if self._ics_recall_ics_recall is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._ics_recall_ics_recall
 
    @property
    def ics(self):
        """ ICS (Interface Contact Similarity) score
 
        Combination of :attr:`~ics_precision` and :attr:`~ics_recall`
        using the F1-measure
 
        :type: :class:`float`
        """
        if self._ics_ics is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._ics_ics
 
    @property
    def per_interface_ics_precision(self):
        """ Per-interface ICS precision
 
        :attr:`~ics_precision` for each interface in
        :attr:`~contact_target_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_ics_precision_per_interface_ics_precision is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._per_interface_ics_precision_per_interface_ics_precision
 
 
    @property
    def per_interface_ics_recall(self):
        """ Per-interface ICS recall
 
        :attr:`~ics_recall` for each interface in
        :attr:`~contact_target_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_ics_recall_per_interface_ics_recall is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._per_interface_ics_recall_per_interface_ics_recall
 
    @property
    def per_interface_ics(self):
        """ Per-interface ICS (Interface Contact Similarity) score
 
        :attr:`~ics` for each interface in 
        :attr:`~contact_target_interfaces`
 
        :type: :class:`float`
        """
 
        if self._per_interface_ics_per_interface_ics is None:
            self._compute_ics_scores_compute_ics_scores()
        return self._per_interface_ics_per_interface_ics
    
 
    @property
    def ips_precision(self):
        """ Fraction of model interface residues that are also interface
        residues in target
 
        :type: :class:`float`
        """
        if self._ips_precision_ips_precision is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._ips_precision_ips_precision
    
    @property
    def ips_recall(self):
        """ Fraction of target interface residues that are also interface
        residues in model
 
        :type: :class:`float`
        """
        if self._ips_recall_ips_recall is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._ips_recall_ips_recall
 
    @property
    def ips(self):
        """ IPS (Interface Patch Similarity) score
 
        Jaccard coefficient of interface residues in target and their mapped
        counterparts in model
 
        :type: :class:`float`
        """
        if self._ips_ips is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._ips_ips
 
    @property
    def per_interface_ips_precision(self):
        """ Per-interface IPS precision
 
        :attr:`~ips_precision` for each interface in
        :attr:`~contact_target_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_ips_precision_per_interface_ips_precision is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._per_interface_ips_precision_per_interface_ips_precision
 
 
    @property
    def per_interface_ips_recall(self):
        """ Per-interface IPS recall
 
        :attr:`~ips_recall` for each interface in
        :attr:`~contact_target_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
        if self._per_interface_ics_recall_per_interface_ics_recall is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._per_interface_ips_recall_per_interface_ips_recall
 
    @property
    def per_interface_ips(self):
        """ Per-interface IPS (Interface Patch Similarity) score
 
        :attr:`~ips` for each interface in 
        :attr:`~contact_target_interfaces`
 
        :type: :class:`list` of :class:`float`
        """
 
        if self._per_interface_ips_per_interface_ips is None:
            self._compute_ips_scores_compute_ips_scores()
        return self._per_interface_ips_per_interface_ips
 
    @property
    def dockq_target_interfaces(self):
        """ Interfaces in :attr:`target` that are relevant for DockQ
 
        All interfaces in :attr:`~target` with non-zero contacts that are
        relevant for DockQ. Chain names are lexicographically sorted.
 
        :type: :class:`list` of :class:`tuple` with 2 elements each:
               (trg_ch1, trg_ch2)
        """
        if self._dockq_target_interfaces_dockq_target_interfaces is None:
            
            # interacting chains are identified with ContactEntity
            contact_d = 5.0
            if self.dockq_capri_peptidedockq_capri_peptide:
                contact_d = 4.0
            cent = ContactEntity(self.targettarget, contact_mode = "aa",
                                 contact_d = contact_d)
 
            # fetch lexicographically sorted interfaces
            interfaces = cent.interacting_chains
            interfaces = [(min(x[0],x[1]), max(x[0],x[1])) for x in interfaces]
 
            # select the ones with only peptides involved
            pep_seqs = set([s.GetName() for s in self.chain_mapperchain_mapper.polypep_seqs])
            self._dockq_target_interfaces_dockq_target_interfaces = list()
            for interface in interfaces:
                if interface[0] in pep_seqs and interface[1] in pep_seqs:
                    self._dockq_target_interfaces_dockq_target_interfaces.append(interface)
 
        return self._dockq_target_interfaces_dockq_target_interfaces
 
    @property
    def dockq_interfaces(self):
        """ Interfaces in :attr:`dockq_target_interfaces` that can be mapped
        to model
 
        Target chain names are lexicographically sorted
 
        :type: :class:`list` of :class:`tuple` with 4 elements each:
               (trg_ch1, trg_ch2, mdl_ch1, mdl_ch2)
        """
        if self._dockq_interfaces_dockq_interfaces is None:
            self._dockq_interfaces_dockq_interfaces = list()
            flat_mapping = self.mappingmapping.GetFlatMapping()
            for i in self.dockq_target_interfacesdockq_target_interfaces:
                if i[0] in flat_mapping and i[1] in flat_mapping:
                    self._dockq_interfaces_dockq_interfaces.append((i[0], i[1],
                                                   flat_mapping[i[0]],
                                                   flat_mapping[i[1]]))
        return self._dockq_interfaces_dockq_interfaces
    
    @property
    def dockq_scores(self):
        """ DockQ scores for interfaces in :attr:`~dockq_interfaces` 
 
        :class:`list` of :class:`float`
        """
        if self._dockq_scores_dockq_scores is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._dockq_scores_dockq_scores
 
    @property
    def fnat(self):
        """ fnat scores for interfaces in :attr:`~dockq_interfaces` 
 
        fnat: Fraction of native contacts that are also present in model
 
        :class:`list` of :class:`float`
        """
        if self._fnat_fnat is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._fnat_fnat
 
    @property
    def nnat(self):
        """ N native contacts for interfaces in :attr:`~dockq_interfaces` 
 
        :class:`list` of :class:`int`
        """
        if self._nnat_nnat is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._nnat_nnat
 
    @property
    def nmdl(self):
        """ N model contacts for interfaces in :attr:`~dockq_interfaces` 
 
        :class:`list` of :class:`int`
        """
        if self._nmdl_nmdl is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._nmdl_nmdl
 
    @property
    def fnonnat(self):
        """ fnonnat scores for interfaces in :attr:`~dockq_interfaces` 
 
        fnat: Fraction of model contacts that are not present in target
 
        :class:`list` of :class:`float`
        """
        if self._fnonnat_fnonnat is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._fnonnat_fnonnat
 
    @property
    def irmsd(self):
        """ irmsd scores for interfaces in :attr:`~dockq_interfaces` 
 
        irmsd: RMSD of interface (RMSD computed on N, CA, C, O atoms) which
        consists of each residue that has at least one heavy atom within 10A of
        other chain.
 
        :class:`list` of :class:`float`
        """
        if self._irmsd_irmsd is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._irmsd_irmsd
 
    @property
    def lrmsd(self):
        """ lrmsd scores for interfaces in :attr:`~dockq_interfaces` 
 
        lrmsd: The interfaces are superposed based on the receptor (rigid
        min RMSD superposition) and RMSD for the ligand is reported.
        Superposition and RMSD are based on N, CA, C and O positions,
        receptor is the chain contributing to the interface with more
        residues in total.
 
        :class:`list` of :class:`float`
        """
        if self._lrmsd_lrmsd is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._lrmsd_lrmsd
        
    @property
    def dockq_ave(self):
        """ Average of DockQ scores in :attr:`dockq_scores`
 
        In its original implementation, DockQ only operates on single
        interfaces. Thus the requirement to combine scores for higher order
        oligomers.
 
        :type: :class:`float`
        """
        if self._dockq_ave_dockq_ave is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._dockq_ave_dockq_ave
    
    @property
    def dockq_wave(self):
        """ Same as :attr:`dockq_ave`, weighted by native contacts
 
        :type: :class:`float`
        """
        if self._dockq_wave_dockq_wave is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._dockq_wave_dockq_wave
        
    @property
    def dockq_ave_full(self):
        """ Same as :attr:`~dockq_ave` but penalizing for missing interfaces
 
        Interfaces that are not covered in model are added as 0.0
        in average computation.
 
        :type: :class:`float`
        """
        if self._dockq_ave_full_dockq_ave_full is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._dockq_ave_full_dockq_ave_full
    
    @property
    def dockq_wave_full(self):
        """ Same as :attr:`~dockq_ave_full`, but weighted
 
        Interfaces that are not covered in model are added as 0.0 in
        average computations and the respective weights are derived from
        number of contacts in respective target interface. 
        """
        if self._dockq_wave_full_dockq_wave_full is None:
            self._compute_dockq_scores_compute_dockq_scores()
        return self._dockq_wave_full_dockq_wave_full
 
    @property
    def mapped_target_pos(self):
        """ Mapped representative positions in target
 
        Thats CA positions for peptide residues and C3' positions for
        nucleotides. Has same length as :attr:`~mapped_model_pos` and mapping
        is based on :attr:`~mapping`.
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._mapped_target_pos_mapped_target_pos is None:
            self._extract_mapped_pos_extract_mapped_pos()
        return self._mapped_target_pos_mapped_target_pos
 
    @property
    def mapped_model_pos(self):
        """ Mapped representative positions in model
 
        Thats CA positions for peptide residues and C3' positions for
        nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
        is based on :attr:`~mapping`.
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._mapped_model_pos_mapped_model_pos is None:
            self._extract_mapped_pos_extract_mapped_pos()
        return self._mapped_model_pos_mapped_model_pos
 
    @property
    def transformed_mapped_model_pos(self):
        """ :attr:`~mapped_model_pos` with :attr:`~transform` applied
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._transformed_mapped_model_pos_transformed_mapped_model_pos is None:
            self._transformed_mapped_model_pos_transformed_mapped_model_pos = \
            geom.Vec3List(self.mapped_model_posmapped_model_pos)
            self._transformed_mapped_model_pos_transformed_mapped_model_pos.ApplyTransform(self.transformtransform)
        return self._transformed_mapped_model_pos_transformed_mapped_model_pos
 
    @property
    def n_target_not_mapped(self):
        """ Number of target residues which have no mapping to model
 
        :type: :class:`int`
        """
        if self._n_target_not_mapped_n_target_not_mapped is None:
            self._extract_mapped_pos_extract_mapped_pos()
        return self._n_target_not_mapped_n_target_not_mapped
 
    @property
    def transform(self):
        """ Transform: :attr:`~mapped_model_pos` onto :attr:`~mapped_target_pos`
 
        Computed using Kabsch minimal rmsd algorithm
 
        :type: :class:`ost.geom.Mat4`
        """
        if self._transform_transform is None:
            try:
                res = mol.alg.SuperposeSVD(self.mapped_model_posmapped_model_pos,
                                           self.mapped_target_posmapped_target_pos)
                self._transform_transform = res.transformation
            except:
                self._transform_transform = geom.Mat4()
        return self._transform_transform
 
    @property
    def rigid_mapped_target_pos(self):
        """ Mapped representative positions in target
 
        Thats CA positions for peptide residues and C3' positions for
        nucleotides. Has same length as :attr:`~rigid_mapped_model_pos` and mapping
        is based on :attr:`~rigid_mapping`.
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._rigid_mapped_target_pos_rigid_mapped_target_pos is None:
            self._extract_rigid_mapped_pos_extract_rigid_mapped_pos()
        return self._rigid_mapped_target_pos_rigid_mapped_target_pos
 
    @property
    def rigid_mapped_model_pos(self):
        """ Mapped representative positions in model
 
        Thats CA positions for peptide residues and C3' positions for
        nucleotides. Has same length as :attr:`~mapped_target_pos` and mapping
        is based on :attr:`~rigid_mapping`.
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._rigid_mapped_model_pos_rigid_mapped_model_pos is None:
            self._extract_rigid_mapped_pos_extract_rigid_mapped_pos()
        return self._rigid_mapped_model_pos_rigid_mapped_model_pos
 
    @property
    def rigid_transformed_mapped_model_pos(self):
        """ :attr:`~rigid_mapped_model_pos` with :attr:`~rigid_transform` applied
 
        :type: :class:`ost.geom.Vec3List`
        """
        if self._rigid_transformed_mapped_model_pos_rigid_transformed_mapped_model_pos is None:
            self._rigid_transformed_mapped_model_pos_rigid_transformed_mapped_model_pos = \
            geom.Vec3List(self.rigid_mapped_model_posrigid_mapped_model_pos)
            self._rigid_transformed_mapped_model_pos_rigid_transformed_mapped_model_pos.ApplyTransform(self.rigid_transformrigid_transform)
        return self._rigid_transformed_mapped_model_pos_rigid_transformed_mapped_model_pos
 
    @property
    def rigid_n_target_not_mapped(self):
        """ Number of target residues which have no rigid mapping to model
 
        :type: :class:`int`
        """
        if self._rigid_n_target_not_mapped_rigid_n_target_not_mapped is None:
            self._extract_rigid_mapped_pos_extract_rigid_mapped_pos()
        return self._rigid_n_target_not_mapped_rigid_n_target_not_mapped
 
    @property
    def rigid_transform(self):
        """ Transform: :attr:`~rigid_mapped_model_pos` onto :attr:`~rigid_mapped_target_pos`
 
        Computed using Kabsch minimal rmsd algorithm
 
        :type: :class:`ost.geom.Mat4`
        """
        if self._rigid_transform_rigid_transform is None:
            try:
                res = mol.alg.SuperposeSVD(self.rigid_mapped_model_posrigid_mapped_model_pos,
                                           self.rigid_mapped_target_posrigid_mapped_target_pos)
                self._rigid_transform_rigid_transform = res.transformation
            except:
                self._rigid_transform_rigid_transform = geom.Mat4()
        return self._rigid_transform_rigid_transform
 
    @property
    def gdt_05(self):
        """ Fraction CA (C3' for nucleotides) that can be superposed within 0.5A
 
        Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
        Similar iterative algorithm as LGA tool
 
        :type: :class:`float` 
        """
        if self._gdt_05_gdt_05 is None:
            n, m = GDT(self.rigid_mapped_model_posrigid_mapped_model_pos, self.rigid_mapped_target_posrigid_mapped_target_pos, 7, 1000, 0.5)
            n_full = len(self.rigid_mapped_target_posrigid_mapped_target_pos) + self.rigid_n_target_not_mappedrigid_n_target_not_mapped
            if n_full > 0:
                self._gdt_05_gdt_05 = float(n) / n_full
            else:
                self._gdt_05_gdt_05 = 0.0
        return self._gdt_05_gdt_05
 
    @property
    def gdt_1(self):
        """ Fraction CA (C3' for nucleotides) that can be superposed within 1.0A
 
        Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
        Similar iterative algorithm as LGA tool
 
        :type: :class:`float` 
        """
        if self._gdt_1_gdt_1 is None:
            n, m = GDT(self.rigid_mapped_model_posrigid_mapped_model_pos, self.rigid_mapped_target_posrigid_mapped_target_pos, 7, 1000, 1.0)
            n_full = len(self.rigid_mapped_target_posrigid_mapped_target_pos) + self.rigid_n_target_not_mappedrigid_n_target_not_mapped
            if n_full > 0:
                self._gdt_1_gdt_1 = float(n) / n_full
            else:
                self._gdt_1_gdt_1 = 0.0
        return self._gdt_1_gdt_1
 
    @property
    def gdt_2(self):
        """ Fraction CA (C3' for nucleotides) that can be superposed within 2.0A
 
        Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
        Similar iterative algorithm as LGA tool
 
 
        :type: :class:`float` 
        """
        if self._gdt_2_gdt_2 is None:
            n, m = GDT(self.rigid_mapped_model_posrigid_mapped_model_pos, self.rigid_mapped_target_posrigid_mapped_target_pos, 7, 1000, 2.0)
            n_full = len(self.rigid_mapped_target_posrigid_mapped_target_pos) + self.rigid_n_target_not_mappedrigid_n_target_not_mapped
            if n_full > 0:
                self._gdt_2_gdt_2 = float(n) / n_full
            else:
                self._gdt_2_gdt_2 = 0.0
        return self._gdt_2_gdt_2
 
    @property
    def gdt_4(self):
        """ Fraction CA (C3' for nucleotides) that can be superposed within 4.0A
 
        Uses :attr:`~rigid_mapped_model_pos` and :attr:`~rigid_mapped_target_pos`.
        Similar iterative algorithm as LGA tool
 
        :type: :class:`float` 
        """
        if self._gdt_4_gdt_4 is None:
            n, m = GDT(self.rigid_mapped_model_posrigid_mapped_model_pos, self.rigid_mapped_target_posrigid_mapped_target_pos, 7, 1000, 4.0)
            n_full = len(self.rigid_mapped_target_posrigid_mapped_target_pos) + self.rigid_n_target_not_mappedrigid_n_target_not_mapped
            if n_full > 0:
                self._gdt_4_gdt_4 = float(n) / n_full
            else:
                self._gdt_4_gdt_4 = 0.0
        return self._gdt_4_gdt_4
 
    @property
    def gdt_8(self):
        """ Fraction CA (C3' for nucleotides) that can be superposed within 8.0A
 
        Similar iterative algorithm as LGA tool
 
        :type: :class:`float` 
        """
        if self._gdt_8_gdt_8 is None:
            n, m = GDT(self.rigid_mapped_model_posrigid_mapped_model_pos, self.rigid_mapped_target_posrigid_mapped_target_pos, 7, 1000, 8.0)
            n_full = len(self.rigid_mapped_target_posrigid_mapped_target_pos) + self.rigid_n_target_not_mappedrigid_n_target_not_mapped
            if n_full > 0:
                self._gdt_8_gdt_8 = float(n) / n_full
            else:
                self._gdt_8_gdt_8 = 0.0
        return self._gdt_8_gdt_8
    
 
    @property
    def gdtts(self):
        """ avg GDT with thresholds: 8.0A, 4.0A, 2.0A and 1.0A
 
        :type: :class:`float`
        """
        if self._gdtts_gdtts is None:
            LogScript("Computing GDT-TS score")
            self._gdtts_gdtts = (self.gdt_1gdt_1 + self.gdt_2gdt_2 + self.gdt_4gdt_4 + self.gdt_8gdt_8) / 4
        return self._gdtts_gdtts
 
    @property
    def gdtha(self):
        """ avg GDT with thresholds: 4.0A, 2.0A, 1.0A and 0.5A
 
        :type: :class:`float`
        """
        if self._gdtha_gdtha is None:
            LogScript("Computing GDT-HA score")
            self._gdtha_gdtha = (self.gdt_05gdt_05 + self.gdt_1gdt_1 + self.gdt_2gdt_2 + self.gdt_4gdt_4) / 4
        return self._gdtha_gdtha
 
    @property
    def rmsd(self):
        """ RMSD
 
        Computed on :attr:`~rigid_transformed_mapped_model_pos` and
        :attr:`rigid_mapped_target_pos`
 
        :type: :class:`float`
        """
        if self._rmsd_rmsd is None:
            LogScript("Computing RMSD")
            self._rmsd_rmsd = \
            self.rigid_mapped_target_posrigid_mapped_target_pos.GetRMSD(self.rigid_transformed_mapped_model_posrigid_transformed_mapped_model_pos)
        return self._rmsd_rmsd
 
    @property
    def cad_score(self):
        """ The global CAD atom-atom (AA) score
 
        Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
        is used.
 
        :type: :class:`float`
        """
        if self._cad_score_cad_score is None:
            self._compute_cad_score_compute_cad_score()
        return self._cad_score_cad_score
 
    @property
    def local_cad_score(self):
        """ The per-residue CAD atom-atom (AA) scores
 
        Computed based on :attr:`~model`. In case of oligomers, :attr:`~mapping`
        is used.
 
        :type: :class:`dict`
        """
        if self._local_cad_score_local_cad_score is None:
            self._compute_cad_score_compute_cad_score()
        return self._local_cad_score_local_cad_score
 
    @property
    def patch_qs(self):
        """ Patch QS-scores for each residue in :attr:`model_interface_residues`
 
        Representative patches for each residue r in chain c are computed as
        follows:
    
        * mdl_patch_one: All residues in c with CB (CA for GLY) positions within
          8A of r and within 12A of residues from any other chain.
        * mdl_patch_two: Closest residue x to r in any other chain gets
          identified. Patch is then constructed by selecting all residues from
          any other chain within 8A of x and within 12A from any residue in c.
        * trg_patch_one: Chain name and residue number based mapping from
          mdl_patch_one
        * trg_patch_two: Chain name and residue number based mapping from
          mdl_patch_two
 
        Results are stored in the same manner as
        :attr:`model_interface_residues`, with corresponding scores instead of
        residue numbers. Scores for residues which are not
        :class:`mol.ChemType.AMINOACIDS` are set to None. Additionally,
        interface patches are derived from :attr:`model`. If they contain
        residues which are not covered by :attr:`target`, the score is set to
        None too.
 
        :type: :class:`dict` with chain names as key and and :class:`list`
                with scores of the respective interface residues.
        """
        if self._patch_qs_patch_qs is None:
            self._compute_patchqs_scores_compute_patchqs_scores()
        return self._patch_qs_patch_qs
 
    @property
    def patch_dockq(self):
        """ Same as :attr:`patch_qs` but for DockQ scores
        """
        if self._patch_dockq_patch_dockq is None:
            self._compute_patchdockq_scores_compute_patchdockq_scores()
        return self._patch_dockq_patch_dockq
 
    @property
    def tm_score(self):
        """ TM-score computed with USalign
 
        USalign executable can be specified with usalign_exec kwarg at Scorer
        construction, an OpenStructure internal copy of the USalign code is
        used otherwise.
 
        :type: :class:`float`
        """
        if self._tm_score_tm_score is None:
            self._compute_tmscore_compute_tmscore()
        return self._tm_score_tm_score
 
    @property
    def usalign_mapping(self):
        """ Mapping computed with USalign
 
        Dictionary with target chain names as key and model chain names as
        values. No guarantee that all chains are mapped. USalign executable
        can be specified with usalign_exec kwarg at Scorer construction, an
        OpenStructure internal copy of the USalign code is used otherwise.
 
        :type: :class:`dict`
        """
        if self._usalign_mapping_usalign_mapping is None:
            self._compute_tmscore_compute_tmscore()
        return self._usalign_mapping_usalign_mapping
 
    def _aln_helper(self, target, model):
        # perform required alignments - cannot take the alignments from the
        # mapping results as we potentially remove stuff there as compared
        # to self.model and self.target
        trg_seqs = dict()
        for ch in target.chains:
            cname = ch.GetName()
            s = ''.join([r.one_letter_code for r in ch.residues])
            s = seq.CreateSequence(ch.GetName(), s)
            s.AttachView(target.Select(f"cname={mol.QueryQuoteName(cname)}"))
            trg_seqs[ch.GetName()] = s
        mdl_seqs = dict()
        for ch in model.chains:
            cname = ch.GetName()
            s = ''.join([r.one_letter_code for r in ch.residues])
            s = seq.CreateSequence(cname, s)
            s.AttachView(model.Select(f"cname={mol.QueryQuoteName(cname)}"))
            mdl_seqs[ch.GetName()] = s
 
        alns = list()
        trg_pep_chains = [s.GetName() for s in self.chain_mapperchain_mapper.polypep_seqs]
        trg_nuc_chains = [s.GetName() for s in self.chain_mapperchain_mapper.polynuc_seqs]
        trg_pep_chains = set(trg_pep_chains)
        trg_nuc_chains = set(trg_nuc_chains)
        for trg_ch, mdl_ch in self.mappingmapping.GetFlatMapping().items():
            if mdl_ch in mdl_seqs and trg_ch in trg_seqs:
                if trg_ch in trg_pep_chains:
                    stype = mol.ChemType.AMINOACIDS
                elif trg_ch in trg_nuc_chains:
                    stype = mol.ChemType.NUCLEOTIDES
                else:
                    raise RuntimeError("Chain name inconsistency... ask "
                                       "Gabriel")
                alns.append(self.chain_mapperchain_mapper.Align(trg_seqs[trg_ch],
                                                    mdl_seqs[mdl_ch],
                                                    stype))
                alns[-1].AttachView(0, trg_seqs[trg_ch].GetAttachedView())
                alns[-1].AttachView(1, mdl_seqs[mdl_ch].GetAttachedView())
        return alns
 
    def _compute_pepnuc_aln(self):
        query = "peptide=true or nucleotide=true"
        pep_nuc_target = self.target_origtarget_orig.Select(query)
        pep_nuc_model = self.model_origmodel_orig.Select(query)
        self._pepnuc_aln_pepnuc_aln = self._aln_helper_aln_helper(pep_nuc_target, pep_nuc_model)
 
    def _compute_aln(self):
        self._aln_aln = self._aln_helper_aln_helper(self.targettarget, self.modelmodel)
 
    def _compute_stereochecked_aln(self):
        self._stereochecked_aln_stereochecked_aln = self._aln_helper_aln_helper(self.stereochecked_targetstereochecked_target,
                                                   self.stereochecked_modelstereochecked_model)
 
    def _compute_lddt(self):
        LogScript("Computing all-atom lDDT")
        # lDDT requires a flat mapping with mdl_ch as key and trg_ch as value
        flat_mapping = self.mappingmapping.GetFlatMapping(mdl_as_key=True)
 
        # make alignments accessible by mdl seq name
        stereochecked_alns = dict()
        for aln in self.stereochecked_alnstereochecked_aln:
            mdl_seq = aln.GetSequence(1)
            stereochecked_alns[mdl_seq.name] = aln
        alns = dict()
        for aln in self.alnaln:
            mdl_seq = aln.GetSequence(1)
            alns[mdl_seq.name] = aln
 
        # score variables to be set
        lddt_score = None
        local_lddt = None
 
        if self.lddt_no_stereocheckslddt_no_stereochecks:
            lddt_chain_mapping = dict()
            for mdl_ch, trg_ch in flat_mapping.items():
                if mdl_ch in alns:
                    lddt_chain_mapping[mdl_ch] = trg_ch
            lddt_score = self.lddt_scorerlddt_scorer.lDDT(self.modelmodel,
                                               chain_mapping = lddt_chain_mapping,
                                               residue_mapping = alns,
                                               check_resnames=False,
                                               local_lddt_prop="lddt",
                                               add_mdl_contacts = self.lddt_add_mdl_contactslddt_add_mdl_contacts)[0]
            local_lddt = dict()
            for r in self.modelmodel.residues:
                cname = r.GetChain().GetName()
                if cname not in local_lddt:
                    local_lddt[cname] = dict()
                if r.HasProp("lddt"):
                    score = round(r.GetFloatProp("lddt"), 3)
                    local_lddt[cname][r.GetNumber()] = score
                else:
                    # not covered by trg or skipped in chain mapping procedure
                    # the latter happens if its part of a super short chain
                    local_lddt[cname][r.GetNumber()] = None
 
        else:
            lddt_chain_mapping = dict()
            for mdl_ch, trg_ch in flat_mapping.items():
                if mdl_ch in stereochecked_alns:
                    lddt_chain_mapping[mdl_ch] = trg_ch
            lddt_score = self.lddt_scorerlddt_scorer.lDDT(self.stereochecked_modelstereochecked_model,
                                               chain_mapping = lddt_chain_mapping,
                                               residue_mapping = stereochecked_alns,
                                               check_resnames=False,
                                               local_lddt_prop="lddt",
                                               add_mdl_contacts = self.lddt_add_mdl_contactslddt_add_mdl_contacts)[0]
            local_lddt = dict()
            for r in self.modelmodel.residues:
                cname = r.GetChain().GetName()
                if cname not in local_lddt:
                    local_lddt[cname] = dict()
                if r.HasProp("lddt"):
                    score = round(r.GetFloatProp("lddt"), 3)
                    local_lddt[cname][r.GetNumber()] = score
                else:
                    # rsc => residue stereo checked...
                    mdl_res = self.stereochecked_modelstereochecked_model.FindResidue(cname, r.GetNumber())
                    if mdl_res.IsValid():
                        # not covered by trg or skipped in chain mapping procedure
                        # the latter happens if its part of a super short chain
                        local_lddt[cname][r.GetNumber()] = None
                    else:
                        # opt 1: removed by stereochecks => assign 0.0
                        # opt 2: removed by stereochecks AND not covered by ref
                        #        => assign None
                        # fetch trg residue from non-stereochecked aln
                        trg_r = None
                        if cname in flat_mapping:
                            for col in alns[cname]:
                                if col[0] != '-' and col[1] != '-':
                                    if col.GetResidue(1).number == r.number:
                                        trg_r = col.GetResidue(0)
                                        break
                        if trg_r is None:
                            local_lddt[cname][r.GetNumber()] = None
                        else:
                            local_lddt[cname][r.GetNumber()] = 0.0
 
        self._lddt_lddt = lddt_score
        self._local_lddt_local_lddt = local_lddt
 
    def _compute_bb_lddt(self):
        LogScript("Computing backbone lDDT")
        # make alignments accessible by mdl seq name
        alns = dict()
        for aln in self.alnaln:
            mdl_seq = aln.GetSequence(1)
            alns[mdl_seq.name] = aln
 
        # lDDT requires a flat mapping with mdl_ch as key and trg_ch as value
        flat_mapping = self.mappingmapping.GetFlatMapping(mdl_as_key=True)
        lddt_chain_mapping = dict()
        for mdl_ch, trg_ch in flat_mapping.items():
            if mdl_ch in alns:
                lddt_chain_mapping[mdl_ch] = trg_ch
 
        lddt_score = self.bb_lddt_scorerbb_lddt_scorer.lDDT(self.modelmodel,
                                              chain_mapping = lddt_chain_mapping,
                                              residue_mapping = alns,
                                              check_resnames=False,
                                              local_lddt_prop="bb_lddt",
                                              add_mdl_contacts = self.lddt_add_mdl_contactslddt_add_mdl_contacts)[0]
        local_lddt = dict()
        for r in self.modelmodel.residues:
            cname = r.GetChain().GetName()
            if cname not in local_lddt:
                local_lddt[cname] = dict()
            if r.HasProp("bb_lddt"):
                score = round(r.GetFloatProp("bb_lddt"), 3)
                local_lddt[cname][r.GetNumber()] = score
            else:
                # not covered by trg or skipped in chain mapping procedure
                # the latter happens if its part of a super short chain
                local_lddt[cname][r.GetNumber()] = None
 
        self._bb_lddt_bb_lddt = lddt_score
        self._bb_local_lddt_bb_local_lddt = local_lddt
 
    def _compute_ilddt(self):
        LogScript("Computing all-atom ilDDT")
        # lDDT requires a flat mapping with mdl_ch as key and trg_ch as value
        flat_mapping = self.mappingmapping.GetFlatMapping(mdl_as_key=True)
 
        if self.lddt_no_stereocheckslddt_no_stereochecks:
            alns = dict()
            for aln in self.alnaln:
                mdl_seq = aln.GetSequence(1)
                alns[mdl_seq.name] = aln
            lddt_chain_mapping = dict()
            for mdl_ch, trg_ch in flat_mapping.items():
                if mdl_ch in alns:
                    lddt_chain_mapping[mdl_ch] = trg_ch
            self._ilddt_ilddt = self.lddt_scorerlddt_scorer.lDDT(self.modelmodel,
                                                chain_mapping = lddt_chain_mapping,
                                                residue_mapping = alns,
                                                check_resnames=False,
                                                local_lddt_prop="lddt",
                                                add_mdl_contacts = self.lddt_add_mdl_contactslddt_add_mdl_contacts,
                                                no_intrachain=True)[0]
        else:
            alns = dict()
            for aln in self.stereochecked_alnstereochecked_aln:
                mdl_seq = aln.GetSequence(1)
                alns[mdl_seq.name] = aln
            lddt_chain_mapping = dict()
            for mdl_ch, trg_ch in flat_mapping.items():
                if mdl_ch in alns:
                    lddt_chain_mapping[mdl_ch] = trg_ch
            self._ilddt_ilddt = self.lddt_scorerlddt_scorer.lDDT(self.stereochecked_modelstereochecked_model,
                                                chain_mapping = lddt_chain_mapping,
                                                residue_mapping = alns,
                                                check_resnames=False,
                                                local_lddt_prop="lddt",
                                                add_mdl_contacts = self.lddt_add_mdl_contactslddt_add_mdl_contacts,
                                                no_intrachain=True)[0]
 
 
    def _compute_qs(self):
        LogScript("Computing global QS-score")
        qs_score_result = self.qs_scorerqs_scorer.Score(self.mappingmapping.mapping)
        self._qs_global_qs_global = qs_score_result.QS_global
        self._qs_best_qs_best = qs_score_result.QS_best
 
    def _compute_per_interface_qs_scores(self):
        LogScript("Computing per-interface QS-score")
        self._per_interface_qs_global_per_interface_qs_global = list()
        self._per_interface_qs_best_per_interface_qs_best = list()
 
        for interface in self.qs_interfacesqs_interfaces:
            trg_ch1 = interface[0]
            trg_ch2 = interface[1]
            mdl_ch1 = interface[2]
            mdl_ch2 = interface[3]
            qs_res = self.qs_scorerqs_scorer.ScoreInterface(trg_ch1, trg_ch2,
                                                   mdl_ch1, mdl_ch2)
            self._per_interface_qs_best_per_interface_qs_best.append(qs_res.QS_best)
            self._per_interface_qs_global_per_interface_qs_global.append(qs_res.QS_global)
 
    def _compute_ics_scores(self):
        LogScript("Computing ICS scores")
        contact_scorer_res = self.contact_scorercontact_scorer.ScoreICS(self.mappingmapping.mapping)
        self._ics_precision_ics_precision = contact_scorer_res.precision
        self._ics_recall_ics_recall = contact_scorer_res.recall
        self._ics_ics = contact_scorer_res.ics
        self._per_interface_ics_precision_per_interface_ics_precision = list()
        self._per_interface_ics_recall_per_interface_ics_recall = list()
        self._per_interface_ics_per_interface_ics = list()
        flat_mapping = self.mappingmapping.GetFlatMapping()
        for trg_int in self.contact_target_interfacescontact_target_interfaces:
            trg_ch1 = trg_int[0]
            trg_ch2 = trg_int[1]
            if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
                mdl_ch1 = flat_mapping[trg_ch1]
                mdl_ch2 = flat_mapping[trg_ch2]
                res = self.contact_scorercontact_scorer.ScoreICSInterface(trg_ch1, trg_ch2,
                                                            mdl_ch1, mdl_ch2)
                self._per_interface_ics_precision_per_interface_ics_precision.append(res.precision)
                self._per_interface_ics_recall_per_interface_ics_recall.append(res.recall)
                self._per_interface_ics_per_interface_ics.append(res.ics)
            else:
                self._per_interface_ics_precision_per_interface_ics_precision.append(None)
                self._per_interface_ics_recall_per_interface_ics_recall.append(None)
                self._per_interface_ics_per_interface_ics.append(None)
 
    def _compute_ips_scores(self):
        LogScript("Computing IPS scores")
        contact_scorer_res = self.contact_scorercontact_scorer.ScoreIPS(self.mappingmapping.mapping)
        self._ips_precision_ips_precision = contact_scorer_res.precision
        self._ips_recall_ips_recall = contact_scorer_res.recall
        self._ips_ips = contact_scorer_res.ips
 
        self._per_interface_ips_precision_per_interface_ips_precision = list()
        self._per_interface_ips_recall_per_interface_ips_recall = list()
        self._per_interface_ips_per_interface_ips = list()
        flat_mapping = self.mappingmapping.GetFlatMapping()
        for trg_int in self.contact_target_interfacescontact_target_interfaces:
            trg_ch1 = trg_int[0]
            trg_ch2 = trg_int[1]
            if trg_ch1 in flat_mapping and trg_ch2 in flat_mapping:
                mdl_ch1 = flat_mapping[trg_ch1]
                mdl_ch2 = flat_mapping[trg_ch2]
                res = self.contact_scorercontact_scorer.ScoreIPSInterface(trg_ch1, trg_ch2,
                                                            mdl_ch1, mdl_ch2)
                self._per_interface_ips_precision_per_interface_ips_precision.append(res.precision)
                self._per_interface_ips_recall_per_interface_ips_recall.append(res.recall)
                self._per_interface_ips_per_interface_ips.append(res.ips)
            else:
                self._per_interface_ips_precision_per_interface_ips_precision.append(None)
                self._per_interface_ips_recall_per_interface_ips_recall.append(None)
                self._per_interface_ips_per_interface_ips.append(None)
 
    def _compute_dockq_scores(self):
        LogScript("Computing DockQ")
        # lists with values in contact_target_interfaces
        self._dockq_scores_dockq_scores = list()
        self._fnat_fnat = list()
        self._fnonnat_fnonnat = list()
        self._irmsd_irmsd = list()
        self._lrmsd_lrmsd = list()
        self._nnat_nnat = list()
        self._nmdl_nmdl = list()
 
        dockq_alns = dict()
        for aln in self.alnaln:
            trg_s = aln.GetSequence(0)
            mdl_s = aln.GetSequence(1)
            dockq_alns[(trg_s.GetName(), mdl_s.GetName())] = aln
 
        for interface in self.dockq_interfacesdockq_interfaces:
            trg_ch1 = interface[0]
            trg_ch2 = interface[1]
            mdl_ch1 = interface[2]
            mdl_ch2 = interface[3]
            aln1 = dockq_alns[(trg_ch1, mdl_ch1)]
            aln2 = dockq_alns[(trg_ch2, mdl_ch2)]
            if self.dockq_capri_peptidedockq_capri_peptide:
                res = dockq.DockQ(self.modelmodel, self.targettarget, mdl_ch1, mdl_ch2,
                                  trg_ch1, trg_ch2, ch1_aln=aln1,
                                  ch2_aln=aln2, contact_dist_thresh = 4.0,
                                  interface_dist_thresh=8.0,
                                  interface_cb = True)
            else:
                res = dockq.DockQ(self.modelmodel, self.targettarget, mdl_ch1, mdl_ch2,
                                  trg_ch1, trg_ch2, ch1_aln=aln1,
                                  ch2_aln=aln2)
 
            self._fnat_fnat.append(res["fnat"])
            self._fnonnat_fnonnat.append(res["fnonnat"])
            self._irmsd_irmsd.append(res["irmsd"])
            self._lrmsd_lrmsd.append(res["lrmsd"])
            self._dockq_scores_dockq_scores.append(res["DockQ"])
            self._nnat_nnat.append(res["nnat"])
            self._nmdl_nmdl.append(res["nmdl"])
 
        # keep track of native counts in target interfaces which are
        # not covered in model in order to compute
        # dockq_ave_full/dockq_wave_full in the end
        not_covered_counts = list()
        proc_trg_interfaces = set([(x[0], x[1]) for x in self.dockq_interfacesdockq_interfaces])
        for interface in self.dockq_target_interfacesdockq_target_interfaces:
            if interface not in proc_trg_interfaces:
                # let's run DockQ with trg as trg/mdl in order to get the native
                # contacts out - no need to pass alns as the residue numbers
                # match for sure
                trg_ch1 = interface[0]
                trg_ch2 = interface[1]
 
                if self.dockq_capri_peptidedockq_capri_peptide:
                    res = dockq.DockQ(self.targettarget, self.targettarget,
                                      trg_ch1, trg_ch2, trg_ch1, trg_ch2,
                                      contact_dist_thresh = 4.0,
                                      interface_dist_thresh=8.0,
                                      interface_cb = True)
                else:
                    res = dockq.DockQ(self.targettarget, self.targettarget,
                                      trg_ch1, trg_ch2, trg_ch1, trg_ch2)
 
                not_covered_counts.append(res["nnat"])
  
        # there are 4 types of combined scores
        # - simple average
        # - average weighted by native_contacts
        # - the two above including nonmapped_contact_interfaces => set DockQ to 0.0
        scores = np.array(self._dockq_scores_dockq_scores)
        weights = np.array(self._nnat_nnat)
        if len(scores) > 0:
            self._dockq_ave_dockq_ave = np.mean(scores)
        else:
            self._dockq_ave_dockq_ave = 0.0
        self._dockq_wave_dockq_wave = np.sum(np.multiply(weights/np.sum(weights), scores))
        scores = np.append(scores, [0.0]*len(not_covered_counts))
        weights = np.append(weights, not_covered_counts)
        if len(scores) > 0:
            self._dockq_ave_full_dockq_ave_full = np.mean(scores)
        else:
            self._dockq_ave_full_dockq_ave_full = 0.0
        self._dockq_wave_full_dockq_wave_full = np.sum(np.multiply(weights/np.sum(weights),
                                                   scores))
 
    def _extract_mapped_pos(self):
        self._mapped_target_pos_mapped_target_pos = geom.Vec3List()
        self._mapped_model_pos_mapped_model_pos = geom.Vec3List()
        self._n_target_not_mapped_n_target_not_mapped = 0
        processed_trg_chains = set()
        for trg_ch, mdl_ch in self.mappingmapping.GetFlatMapping().items():
            processed_trg_chains.add(trg_ch)
            aln = self.mappingmapping.alns[(trg_ch, mdl_ch)]
            for col in aln:
                if col[0] != '-' and col[1] != '-':
                    trg_res = col.GetResidue(0)
                    mdl_res = col.GetResidue(1)
                    trg_at = trg_res.FindAtom("CA")
                    mdl_at = mdl_res.FindAtom("CA")
                    if not trg_at.IsValid():
                        trg_at = trg_res.FindAtom("C3'")
                    if not mdl_at.IsValid():
                        mdl_at = mdl_res.FindAtom("C3'")
                    self._mapped_target_pos_mapped_target_pos.append(trg_at.GetPos())
                    self._mapped_model_pos_mapped_model_pos.append(mdl_at.GetPos())
                elif col[0] != '-':
                    self._n_target_not_mapped_n_target_not_mapped += 1
        # count number of trg residues from non-mapped chains
        for ch in self.mappingmapping.target.chains:
            if ch.GetName() not in processed_trg_chains:
                self._n_target_not_mapped_n_target_not_mapped += len(ch.residues)
 
 
    def _extract_rigid_mapped_pos(self):
        self._rigid_mapped_target_pos_rigid_mapped_target_pos = geom.Vec3List()
        self._rigid_mapped_model_pos_rigid_mapped_model_pos = geom.Vec3List()
        self._rigid_n_target_not_mapped_rigid_n_target_not_mapped = 0
        processed_trg_chains = set()
        for trg_ch, mdl_ch in self.rigid_mappingrigid_mapping.GetFlatMapping().items():
            processed_trg_chains.add(trg_ch)
            aln = self.rigid_mappingrigid_mapping.alns[(trg_ch, mdl_ch)]
            for col in aln:
                if col[0] != '-' and col[1] != '-':
                    trg_res = col.GetResidue(0)
                    mdl_res = col.GetResidue(1)
                    trg_at = trg_res.FindAtom("CA")
                    mdl_at = mdl_res.FindAtom("CA")
                    if not trg_at.IsValid():
                        trg_at = trg_res.FindAtom("C3'")
                    if not mdl_at.IsValid():
                        mdl_at = mdl_res.FindAtom("C3'")
                    self._rigid_mapped_target_pos_rigid_mapped_target_pos.append(trg_at.GetPos())
                    self._rigid_mapped_model_pos_rigid_mapped_model_pos.append(mdl_at.GetPos())
                elif col[0] != '-':
                    self._rigid_n_target_not_mapped_rigid_n_target_not_mapped += 1
        # count number of trg residues from non-mapped chains
        for ch in self.rigid_mappingrigid_mapping.target.chains:
            if ch.GetName() not in processed_trg_chains:
                self._rigid_n_target_not_mapped_rigid_n_target_not_mapped += len(ch.residues)
 
 
    def _compute_cad_score(self):
        if not self.resnum_alignmentsresnum_alignments:
            raise RuntimeError("CAD score computations rely on residue numbers "
                               "that are consistent between target and model "
                               "chains, i.e. only work if resnum_alignments "
                               "is True at Scorer construction.")
        try:
            LogScript("Computing CAD score")
            cad_score_exec = \
            settings.Locate("voronota-cadscore",
                            explicit_file_name=self.cad_score_execcad_score_exec)
        except Exception as e:
            raise RuntimeError("voronota-cadscore must be in PATH for CAD "
                               "score scoring") from e
        cad_bin_dir = os.path.dirname(cad_score_exec)
        m = self.mappingmapping.GetFlatMapping(mdl_as_key=True)
        cad_result = cadscore.CADScore(self.modelmodel, self.targettarget,
                                       mode = "voronota",
                                       label="localcad",
                                       old_regime=False,
                                       cad_bin_path=cad_bin_dir,
                                       chain_mapping=m)
 
        local_cad = dict()
        for r in self.modelmodel.residues:
            cname = r.GetChain().GetName()
            if cname not in local_cad:
                local_cad[cname] = dict()
            if r.HasProp("localcad"):
                score = round(r.GetFloatProp("localcad"), 3)
                local_cad[cname][r.GetNumber()] = score
            else:
                local_cad[cname][r.GetNumber()] = None
 
        self._cad_score_cad_score = cad_result.globalAA
        self._local_cad_score_local_cad_score = local_cad
 
    def _get_repr_view(self, ent):
        """ Returns view with representative peptide atoms => CB, CA for GLY
    
        Ensures that each residue has exactly one atom with assertions
    
        :param ent: Entity for which you want the representative view
        :param ent: :class:`ost.mol.EntityHandle`/:class:`ost.mol.EntityView`
        :returns: :class:`ost.mol.EntityView` derived from ent
        """
        repr_ent = ent.Select("(aname=\"CB\" or (rname=\"GLY\" and aname=\"CA\"))")
        for r in repr_ent.residues:
            assert(len(r.atoms) == 1)
        return repr_ent
 
    def _get_interface_residues(self, ent):
        """ Get interface residues
    
        Thats all residues having a contact with at least one residue from another
        chain (CB-CB distance <= 8A, CA in case of Glycine)
    
        :param ent: Model for which to extract interface residues
        :type ent: :class:`ost.mol.EntityView`
        :returns: :class:`dict` with chain names as key and and :class:`list`
                  with residue numbers of the respective interface residues.
        """
        # select for representative positions => CB, CA for GLY 
        repr_ent = self._get_repr_view_get_repr_view(ent)
        result = {ch.GetName(): list() for ch in ent.chains}
        for ch in ent.chains:
            cname = ch.GetName()
            sel = repr_ent.Select(f"(cname={mol.QueryQuoteName(cname)} and 8 <> [cname!={mol.QueryQuoteName(cname)}])")
            result[cname] = [r.GetNumber() for r in sel.residues]
        return result
 
    def _do_stereochecks(self):
        """ Perform stereochemistry checks on model and target
        """
        LogInfo("Performing stereochemistry checks on model and target")
        data = stereochemistry.GetDefaultStereoData()
        l_data = stereochemistry.GetDefaultStereoLinkData()
 
        a, b, c, d = stereochemistry.StereoCheck(self.modelmodel, stereo_data = data,
                                                 stereo_link_data = l_data)
        self._stereochecked_model_stereochecked_model = a
        self._model_clashes_model_clashes = b
        self._model_bad_bonds_model_bad_bonds = c
        self._model_bad_angles_model_bad_angles = d
 
        a, b, c, d = stereochemistry.StereoCheck(self.targettarget, stereo_data = data,
                                                 stereo_link_data = l_data)
        self._stereochecked_target_stereochecked_target = a
        self._target_clashes_target_clashes = b
        self._target_bad_bonds_target_bad_bonds = c
        self._target_bad_angles_target_bad_angles = d
 
 
    def _get_interface_patches(self, mdl_ch, mdl_rnum):
        """ Select interface patches representative for specified residue
    
        The patches for specified residue r in chain c are selected as follows:
    
        * mdl_patch_one: All residues in c with CB (CA for GLY) positions within 8A
          of r and within 12A of residues from any other chain.
        * mdl_patch_two: Closest residue x to r in any other chain gets identified.
          Patch is then constructed by selecting all residues from any other chain
          within 8A of x and within 12A from any residue in c.
        * trg_patch_one: Chain name and residue number based mapping from
          mdl_patch_one
        * trg_patch_two: Chain name and residue number based mapping from
          mdl_patch_two
    
        :param mdl_ch: Name of chain in *self.model* of residue of interest
        :type mdl_ch: :class:`str`
        :param mdl_rnum: Residue number of residue of interest
        :type mdl_rnum: :class:`ost.mol.ResNum`
        :returns: Tuple with 5 elements: 1) :class:`bool` flag whether all residues
                  in *mdl* patches are covered in *trg* 2) mtl_patch_one
                  3) mdl_patch_two 4) trg_patch_one 5) trg_patch_two
        """
        # select for representative positions => CB, CA for GLY 
        repr_mdl = self._get_repr_view_get_repr_view(self.modelmodel.Select("peptide=true"))
    
        # get position for specified residue
        r = self.modelmodel.FindResidue(mdl_ch, mdl_rnum)
        if not r.IsValid():
            raise RuntimeError(f"Cannot find residue {mdl_rnum} in chain {mdl_ch}")
        if r.GetName() == "GLY":
            at = r.FindAtom("CA")
        else:
            at = r.FindAtom("CB")
        if not at.IsValid():
            raise RuntimeError("Cannot find interface views for res without CB/CA")
        r_pos = at.GetPos()
    
        # mdl_patch_one contains residues from the same chain as r
        # => all residues within 8A of r and within 12A of any other chain
    
        # q1 selects for everything in same chain and within 8A of r_pos
        q1 = f"(cname={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{r_pos[0]},{r_pos[1]},{r_pos[2]}}})"
        # q2 selects for everything within 12A of any other chain
        q2 = f"(12 <> [cname!={mol.QueryQuoteName(mdl_ch)}])"
        mdl_patch_one = self.modelmodel.CreateEmptyView()
        sel = repr_mdl.Select(" and ".join([q1, q2]))
        for r in sel.residues:
            mdl_r = self.modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
            mdl_patch_one.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
    
        # mdl_patch_two contains residues from all other chains. In detail:
        # the closest residue to r is identified in any other chain, and the
        # patch is filled with residues that are within 8A of that residue and
        # within 12A of chain from r
        sel = repr_mdl.Select(f"(cname!={mol.QueryQuoteName(mdl_ch)})")
        close_stuff = sel.FindWithin(r_pos, 8)
        min_pos = None
        min_dist = 42.0
        for close_at in close_stuff:
            dist = geom.Distance(r_pos, close_at.GetPos())
            if dist < min_dist:
                min_pos = close_at.GetPos()
                min_dist = dist
    
        # q1 selects for everything not in mdl_ch but within 8A of min_pos
        q1 = f"(cname!={mol.QueryQuoteName(mdl_ch)} and 8 <> {{{min_pos[0]},{min_pos[1]},{min_pos[2]}}})"
        # q2 selects for everything within 12A of mdl_ch
        q2 = f"(12 <> [cname={mol.QueryQuoteName(mdl_ch)}])"
        mdl_patch_two = self.modelmodel.CreateEmptyView()
        sel = repr_mdl.Select(" and ".join([q1, q2]))
        for r in sel.residues:
            mdl_r = self.modelmodel.FindResidue(r.GetChain().GetName(), r.GetNumber())
            mdl_patch_two.AddResidue(mdl_r, mol.ViewAddFlag.INCLUDE_ALL)
    
        # transfer mdl residues to trg
        flat_mapping = self.mappingmapping.GetFlatMapping(mdl_as_key=True)
        full_trg_coverage = True
        trg_patch_one = self.targettarget.CreateEmptyView()
        for r in mdl_patch_one.residues:
            trg_r = None
            mdl_cname = r.GetChain().GetName()
            if mdl_cname in flat_mapping:
                aln = self.mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
                for col in aln:
                    if col[0] != '-' and col[1] != '-':
                        if col.GetResidue(1).GetNumber() == r.GetNumber():
                            trg_r = col.GetResidue(0)
                            break
            if trg_r is not None:
                trg_patch_one.AddResidue(trg_r.handle,
                                         mol.ViewAddFlag.INCLUDE_ALL)
            else:
                full_trg_coverage = False
    
        trg_patch_two = self.targettarget.CreateEmptyView()
        for r in mdl_patch_two.residues:
            trg_r = None
            mdl_cname = r.GetChain().GetName()
            if mdl_cname in flat_mapping:
                aln = self.mappingmapping.alns[(flat_mapping[mdl_cname], mdl_cname)]
                for col in aln:
                    if col[0] != '-' and col[1] != '-':
                        if col.GetResidue(1).GetNumber() == r.GetNumber():
                            trg_r = col.GetResidue(0)
                            break
            if trg_r is not None:
                trg_patch_two.AddResidue(trg_r.handle,
                                         mol.ViewAddFlag.INCLUDE_ALL)
            else:
                full_trg_coverage = False
    
        return (full_trg_coverage, mdl_patch_one, mdl_patch_two,
                trg_patch_one, trg_patch_two)
 
    def _compute_patchqs_scores(self):
        LogScript("Computing patch QS-scores")
        self._patch_qs_patch_qs = dict()
        for cname, rnums in self.model_interface_residuesmodel_interface_residues.items():
            scores = list()
            for rnum in rnums:
                score = None
                r = self.modelmodel.FindResidue(cname, rnum)
                if r.IsValid() and r.GetChemType() == mol.ChemType.AMINOACIDS:
                    full_trg_coverage, mdl_patch_one, mdl_patch_two, \
                    trg_patch_one, trg_patch_two = \
                    self._get_interface_patches_get_interface_patches(cname, rnum)
                    if full_trg_coverage:
                        score = self._patchqs_patchqs(mdl_patch_one, mdl_patch_two,
                                              trg_patch_one, trg_patch_two)
                scores.append(score)
            self._patch_qs_patch_qs[cname] = scores
 
    def _compute_patchdockq_scores(self):
        LogScript("Computing patch DockQ scores")
        self._patch_dockq_patch_dockq = dict()
        for cname, rnums in self.model_interface_residuesmodel_interface_residues.items():
            scores = list()
            for rnum in rnums:
                score = None
                r = self.modelmodel.FindResidue(cname, rnum)
                if r.IsValid() and r.GetChemType() == mol.ChemType.AMINOACIDS:
                    full_trg_coverage, mdl_patch_one, mdl_patch_two, \
                    trg_patch_one, trg_patch_two = \
                    self._get_interface_patches_get_interface_patches(cname, rnum)
                    if full_trg_coverage:
                        score = self._patchdockq_patchdockq(mdl_patch_one, mdl_patch_two,
                                                 trg_patch_one, trg_patch_two)
                scores.append(score)
            self._patch_dockq_patch_dockq[cname] = scores
 
    def _patchqs(self, mdl_patch_one, mdl_patch_two, trg_patch_one, trg_patch_two):
        """ Score interface residue patches with QS-score
    
        In detail: Construct two entities with two chains each. First chain
        consists of residues from <x>_patch_one and second chain consists of
        <x>_patch_two. The returned score is the QS-score between the two
        entities
    
        :param mdl_patch_one: Interface patch representing scored residue
        :type mdl_patch_one: :class:`ost.mol.EntityView`
        :param mdl_patch_two: Interface patch representing scored residue
        :type mdl_patch_two: :class:`ost.mol.EntityView`
        :param trg_patch_one: Interface patch representing scored residue
        :type trg_patch_one: :class:`ost.mol.EntityView`
        :param trg_patch_two: Interface patch representing scored residue
        :type trg_patch_two: :class:`ost.mol.EntityView`
        :returns: PatchQS score
        """
        qs_ent_mdl = self._qs_ent_from_patches_qs_ent_from_patches(mdl_patch_one, mdl_patch_two)
        qs_ent_trg = self._qs_ent_from_patches_qs_ent_from_patches(trg_patch_one, trg_patch_two)
    
        alnA = seq.CreateAlignment()
        s = ''.join([r.one_letter_code for r in mdl_patch_one.residues])
        alnA.AddSequence(seq.CreateSequence("A", s))
        s = ''.join([r.one_letter_code for r in trg_patch_one.residues])
        alnA.AddSequence(seq.CreateSequence("A", s))
    
        alnB = seq.CreateAlignment()
        s = ''.join([r.one_letter_code for r in mdl_patch_two.residues])
        alnB.AddSequence(seq.CreateSequence("B", s))
        s = ''.join([r.one_letter_code for r in trg_patch_two.residues])
        alnB.AddSequence(seq.CreateSequence("B", s))
        alns = {("A", "A"): alnA, ("B", "B"): alnB}
    
        scorer = QSScorer(qs_ent_mdl, [["A"], ["B"]], qs_ent_trg, alns)
        score_result = scorer.Score([["A"], ["B"]])
    
        return score_result.QS_global
    
    def _patchdockq(self, mdl_patch_one, mdl_patch_two, trg_patch_one,
                    trg_patch_two):
        """ Score interface residue patches with DockQ
    
        In detail: Construct two entities with two chains each. First chain
        consists of residues from <x>_patch_one and second chain consists of
        <x>_patch_two. The returned score is the QS-score between the two
        entities
    
        :param mdl_patch_one: Interface patch representing scored residue
        :type mdl_patch_one: :class:`ost.mol.EntityView`
        :param mdl_patch_two: Interface patch representing scored residue
        :type mdl_patch_two: :class:`ost.mol.EntityView`
        :param trg_patch_one: Interface patch representing scored residue
        :type trg_patch_one: :class:`ost.mol.EntityView`
        :param trg_patch_two: Interface patch representing scored residue
        :type trg_patch_two: :class:`ost.mol.EntityView`
        :returns: DockQ score
        """
        m = self._qs_ent_from_patches_qs_ent_from_patches(mdl_patch_one, mdl_patch_two)
        t = self._qs_ent_from_patches_qs_ent_from_patches(trg_patch_one, trg_patch_two)
        dockq_result = dockq.DockQ(t, m, "A", "B", "A", "B")
        if dockq_result["nnat"] > 0:
            return dockq_result["DockQ"]
        return 0.0
 
    def _qs_ent_from_patches(self, patch_one, patch_two):
        """ Constructs Entity with two chains named "A" and "B""
    
        Blindly adds all residues from *patch_one* to chain A and residues from
        patch_two to chain B.
        """
        ent = mol.CreateEntity()
        ed = ent.EditXCS()
        added_ch = ed.InsertChain("A")
        for r in patch_one.residues:
            added_r = ed.AppendResidue(added_ch, r.GetName())
            added_r.SetChemClass(str(r.GetChemClass()))
            for a in r.atoms:
                ed.InsertAtom(added_r, a.handle)
        added_ch = ed.InsertChain("B")
        for r in patch_two.residues:
            added_r = ed.AppendResidue(added_ch, r.GetName())
            added_r.SetChemClass(str(r.GetChemClass()))
            for a in r.atoms:
                ed.InsertAtom(added_r, a.handle)
        return ent
 
    def _construct_custom_mapping(self, mapping):
        """ constructs a full blown MappingResult object from simple dict
 
        :param mapping: mapping with trg chains as key and mdl ch as values
        :type mapping: :class:`dict`
        """
        LogInfo("Setting custom chain mapping")
 
        chain_mapper = self.chain_mapperchain_mapper
        chem_mapping, chem_group_alns, mdl = \
        chain_mapper.GetChemMapping(self.modelmodel)
 
        # now that we have a chem mapping, lets do consistency checks
        # - check whether chain names are unique and available in structures
        # - check whether the mapped chains actually map to the same chem groups
        if len(mapping) != len(set(mapping.keys())):
            raise RuntimeError(f"Expect unique trg chain names in mapping. Got "
                               f"{mapping.keys()}")
        if len(mapping) != len(set(mapping.values())):
            raise RuntimeError(f"Expect unique mdl chain names in mapping. Got "
                               f"{mapping.values()}")
 
        trg_chains = set([ch.GetName() for ch in chain_mapper.target.chains])
        mdl_chains = set([ch.GetName() for ch in mdl.chains])
        for k,v in mapping.items():
            if k not in trg_chains:
                raise RuntimeError(f"Target chain \"{k}\" is not available "
                                   f"in target processed for chain mapping "
                                   f"({trg_chains})")
            if v not in mdl_chains:
                raise RuntimeError(f"Model chain \"{v}\" is not available "
                                   f"in model processed for chain mapping "
                                   f"({mdl_chains})")
 
        for trg_ch, mdl_ch in mapping.items():
            trg_group_idx = None
            mdl_group_idx = None
            for idx, group in enumerate(chain_mapper.chem_groups):
                if trg_ch in group:
                    trg_group_idx = idx
                    break
            for idx, group in enumerate(chem_mapping):
                if mdl_ch in group:
                    mdl_group_idx = idx
                    break
            if trg_group_idx is None or mdl_group_idx is None:
                raise RuntimeError("Could not establish a valid chem grouping "
                                   "of chain names provided in custom mapping.")
            
            if trg_group_idx != mdl_group_idx:
                raise RuntimeError(f"Chem group mismatch in custom mapping: "
                                   f"target chain \"{trg_ch}\" groups with the "
                                   f"following chemically equivalent target "
                                   f"chains: "
                                   f"{chain_mapper.chem_groups[trg_group_idx]} "
                                   f"but model chain \"{mdl_ch}\" maps to the "
                                   f"following target chains: "
                                   f"{chain_mapper.chem_groups[mdl_group_idx]}")
 
        pairs = set([(trg_ch, mdl_ch) for trg_ch, mdl_ch in mapping.items()])
        ref_mdl_alns =  \
        chain_mapping._GetRefMdlAlns(chain_mapper.chem_groups,
                                     chain_mapper.chem_group_alignments,
                                     chem_mapping,
                                     chem_group_alns,
                                     pairs = pairs)
 
        # translate mapping format
        final_mapping = list()
        for ref_chains in chain_mapper.chem_groups:
            mapped_mdl_chains = list()
            for ref_ch in ref_chains:
                if ref_ch in mapping:
                    mapped_mdl_chains.append(mapping[ref_ch])
                else:
                    mapped_mdl_chains.append(None)
            final_mapping.append(mapped_mdl_chains)
 
        alns = dict()
        for ref_group, mdl_group in zip(chain_mapper.chem_groups,
                                        final_mapping):
            for ref_ch, mdl_ch in zip(ref_group, mdl_group):
                if ref_ch is not None and mdl_ch is not None:
                    aln = ref_mdl_alns[(ref_ch, mdl_ch)]
                    trg_view = chain_mapper.target.Select(f"cname={mol.QueryQuoteName(ref_ch)}")
                    mdl_view = mdl.Select(f"cname={mol.QueryQuoteName(mdl_ch)}")
                    aln.AttachView(0, trg_view)
                    aln.AttachView(1, mdl_view)
                    alns[(ref_ch, mdl_ch)] = aln
 
        return chain_mapping.MappingResult(chain_mapper.target, mdl,
                                           chain_mapper.chem_groups,
                                           chem_mapping,
                                           final_mapping, alns)
 
    def _compute_tmscore(self):
        res = None
        if self.usalign_execusalign_exec is None:
            LogScript("Computing patch TM-score with USalign exectuable")
            if self.oumoum:
                flat_mapping = self.rigid_mappingrigid_mapping.GetFlatMapping()
                LogInfo("Overriding TM-score chain mapping")
                res = res = bindings.WrappedMMAlign(self.modelmodel, self.targettarget,
                                                    mapping=flat_mapping)
            else:
                res = bindings.WrappedMMAlign(self.modelmodel, self.targettarget)
        else:
            LogScript("Computing patch TM-score with built-in USalign")
            if self.oumoum:
                LogInfo("Overriding TM-score chain mapping")
                flat_mapping = self.rigid_mappingrigid_mapping.GetFlatMapping()
                res = tmtools.USAlign(self.modelmodel, self.targettarget,
                                      usalign = self.usalign_execusalign_exec,
                                      custom_chain_mapping = flat_mapping)
            else:
                res = tmtools.USAlign(self.modelmodel, self.targettarget,
                                      usalign = self.usalign_execusalign_exec)
 
        self._tm_score_tm_score = res.tm_score
        self._usalign_mapping_usalign_mapping = dict()
        for a,b in zip(res.ent1_mapped_chains, res.ent2_mapped_chains):
            self._usalign_mapping_usalign_mapping[b] = a
 
# specify public interface
__all__ = ('lDDTBSScorer', 'Scorer',)