removed downloaded plip

1 files changed, 0 insertions, 479 deletions

diff --git a/plip/structure/detection.py b/plip/structure/detection.py
deleted file mode 100644
index 71fec63..0000000
--- a/plip/structure/detection.py
+++ /dev/null
@@ -1,479 +0,0 @@
-import itertools
-from collections import defaultdict
-from collections import namedtuple
-
-import numpy as np
-from openbabel.openbabel import OBAtomAtomIter
-
-from plip.basic import config, logger
-from plip.basic.supplemental import vecangle, vector, euclidean3d, projection
-from plip.basic.supplemental import whichresnumber, whichrestype, whichchain
-
-logger = logger.get_logger()
-
-def filter_contacts(pairings):
-    """Filter interactions by two criteria:
-    1. No interactions between the same residue (important for intra mode).
-    2. No duplicate interactions (A with B and B with A, also important for intra mode)."""
-    if not config.INTRA:
-        return pairings
-    filtered1_pairings = [p for p in pairings if (p.resnr, p.reschain) != (p.resnr_l, p.reschain_l)]
-    already_considered = []
-    filtered2_pairings = []
-    for contact in filtered1_pairings:
-        try:
-            dist = 'D{}'.format(round(contact.distance, 2))
-        except AttributeError:
-            try:
-                dist = 'D{}'.format(round(contact.distance_ah, 2))
-            except AttributeError:
-                dist = 'D{}'.format(round(contact.distance_aw, 2))
-        res1, res2 = ''.join([str(contact.resnr), contact.reschain]), ''.join(
-            [str(contact.resnr_l), contact.reschain_l])
-        data = {res1, res2, dist}
-        if data not in already_considered:
-            filtered2_pairings.append(contact)
-            already_considered.append(data)
-    return filtered2_pairings
-
-
-##################################################
-# FUNCTIONS FOR DETECTION OF SPECIFIC INTERACTIONS
-##################################################
-
-def hydrophobic_interactions(atom_set_a, atom_set_b):
-    """Detection of hydrophobic pliprofiler between atom_set_a (binding site) and atom_set_b (ligand).
-    Definition: All pairs of qualified carbon atoms within a distance of HYDROPH_DIST_MAX
-    """
-    data = namedtuple('hydroph_interaction', 'bsatom bsatom_orig_idx ligatom ligatom_orig_idx '
-                                             'distance restype resnr reschain restype_l, resnr_l, reschain_l')
-    pairings = []
-    for a, b in itertools.product(atom_set_a, atom_set_b):
-        if a.orig_idx == b.orig_idx:
-            continue
-        e = euclidean3d(a.atom.coords, b.atom.coords)
-        if not config.MIN_DIST < e < config.HYDROPH_DIST_MAX:
-            continue
-        restype, resnr, reschain = whichrestype(a.atom), whichresnumber(a.atom), whichchain(a.atom)
-        restype_l, resnr_l, reschain_l = whichrestype(b.orig_atom), whichresnumber(b.orig_atom), whichchain(b.orig_atom)
-        contact = data(bsatom=a.atom, bsatom_orig_idx=a.orig_idx, ligatom=b.atom, ligatom_orig_idx=b.orig_idx,
-                       distance=e, restype=restype, resnr=resnr,
-                       reschain=reschain, restype_l=restype_l,
-                       resnr_l=resnr_l, reschain_l=reschain_l)
-        pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def hbonds(acceptors, donor_pairs, protisdon, typ):
-    """Detection of hydrogen bonds between sets of acceptors and donor pairs.
-    Definition: All pairs of hydrogen bond acceptor and donors with
-    donor hydrogens and acceptor showing a distance within HBOND DIST MIN and HBOND DIST MAX
-    and donor angles above HBOND_DON_ANGLE_MIN
-    """
-    data = namedtuple('hbond', 'a a_orig_idx d d_orig_idx h distance_ah distance_ad angle type protisdon resnr '
-                               'restype reschain resnr_l restype_l reschain_l sidechain atype dtype')
-    pairings = []
-    for acc, don in itertools.product(acceptors, donor_pairs):
-        if not typ == 'strong':
-            continue
-        # Regular (strong) hydrogen bonds
-        dist_ah = euclidean3d(acc.a.coords, don.h.coords)
-        dist_ad = euclidean3d(acc.a.coords, don.d.coords)
-        if not config.MIN_DIST < dist_ad < config.HBOND_DIST_MAX:
-            continue
-        vec1, vec2 = vector(don.h.coords, don.d.coords), vector(don.h.coords, acc.a.coords)
-        v = vecangle(vec1, vec2)
-        if not v > config.HBOND_DON_ANGLE_MIN:
-            continue
-        protatom = don.d.OBAtom if protisdon else acc.a.OBAtom
-        ligatom = don.d.OBAtom if not protisdon else acc.a.OBAtom
-        is_sidechain_hbond = protatom.GetResidue().GetAtomProperty(protatom, 8)  # Check if sidechain atom
-        resnr = whichresnumber(don.d) if protisdon else whichresnumber(acc.a)
-        resnr_l = whichresnumber(acc.a_orig_atom) if protisdon else whichresnumber(don.d_orig_atom)
-        restype = whichrestype(don.d) if protisdon else whichrestype(acc.a)
-        restype_l = whichrestype(acc.a_orig_atom) if protisdon else whichrestype(don.d_orig_atom)
-        reschain = whichchain(don.d) if protisdon else whichchain(acc.a)
-        rechain_l = whichchain(acc.a_orig_atom) if protisdon else whichchain(don.d_orig_atom)
-        # Next line prevents H-Bonds within amino acids in intermolecular interactions
-        if config.INTRA is not None and whichresnumber(don.d) == whichresnumber(acc.a):
-            continue
-        # Next line prevents backbone-backbone H-Bonds
-        if config.INTRA is not None and protatom.GetResidue().GetAtomProperty(protatom,
-                                                                              8) and ligatom.GetResidue().GetAtomProperty(
-                ligatom, 8):
-            continue
-        contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, d=don.d, d_orig_idx=don.d_orig_idx, h=don.h,
-                       distance_ah=dist_ah, distance_ad=dist_ad, angle=v, type=typ, protisdon=protisdon,
-                       resnr=resnr, restype=restype, reschain=reschain, resnr_l=resnr_l,
-                       restype_l=restype_l, reschain_l=rechain_l, sidechain=is_sidechain_hbond,
-                       atype=acc.a.type, dtype=don.d.type)
-        pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def pistacking(rings_bs, rings_lig):
-    """Return all pi-stackings between the given aromatic ring systems in receptor and ligand."""
-    data = namedtuple(
-        'pistack',
-        'proteinring ligandring distance angle offset type restype resnr reschain restype_l resnr_l reschain_l')
-    pairings = []
-    for r, l in itertools.product(rings_bs, rings_lig):
-        # DISTANCE AND RING ANGLE CALCULATION
-        d = euclidean3d(r.center, l.center)
-        b = vecangle(r.normal, l.normal)
-        a = min(b, 180 - b if not 180 - b < 0 else b)  # Smallest of two angles, depending on direction of normal
-
-        # RING CENTER OFFSET CALCULATION (project each ring center into the other ring)
-        proj1 = projection(l.normal, l.center, r.center)
-        proj2 = projection(r.normal, r.center, l.center)
-        offset = min(euclidean3d(proj1, l.center), euclidean3d(proj2, r.center))
-
-        # RECEPTOR DATA
-        resnr, restype, reschain = whichresnumber(r.atoms[0]), whichrestype(r.atoms[0]), whichchain(r.atoms[0])
-        resnr_l, restype_l, reschain_l = whichresnumber(l.orig_atoms[0]), whichrestype(
-            l.orig_atoms[0]), whichchain(l.orig_atoms[0])
-
-        # SELECTION BY DISTANCE, ANGLE AND OFFSET
-        passed = False
-        if not config.MIN_DIST < d < config.PISTACK_DIST_MAX:
-            continue
-        if 0 < a < config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
-            ptype = 'P'
-            passed = True
-        if 90 - config.PISTACK_ANG_DEV < a < 90 + config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX:
-            ptype = 'T'
-            passed = True
-        if passed:
-            contact = data(proteinring=r, ligandring=l, distance=d, angle=a, offset=offset,
-                           type=ptype, resnr=resnr, restype=restype, reschain=reschain,
-                           resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l)
-            pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def pication(rings, pos_charged, protcharged):
-    """Return all pi-Cation interaction between aromatic rings and positively charged groups.
-    For tertiary and quaternary amines, check also the angle between the ring and the nitrogen.
-    """
-    data = namedtuple(
-        'pication', 'ring charge distance offset type restype resnr reschain restype_l resnr_l reschain_l protcharged')
-    pairings = []
-    if len(rings) == 0 or len(pos_charged) == 0:
-        return pairings
-    for ring in rings:
-        c = ring.center
-        for p in pos_charged:
-            d = euclidean3d(c, p.center)
-            # Project the center of charge into the ring and measure distance to ring center
-            proj = projection(ring.normal, ring.center, p.center)
-            offset = euclidean3d(proj, ring.center)
-            if not config.MIN_DIST < d < config.PICATION_DIST_MAX or not offset < config.PISTACK_OFFSET_MAX:
-                continue
-            if type(p).__name__ == 'lcharge' and p.fgroup == 'tertamine':
-                # Special case here if the ligand has a tertiary amine, check an additional angle
-                # Otherwise, we might have have a pi-cation interaction 'through' the ligand
-                n_atoms = [a_neighbor for a_neighbor in OBAtomAtomIter(p.atoms[0].OBAtom)]
-                n_atoms_coords = [(a.x(), a.y(), a.z()) for a in n_atoms]
-                amine_normal = np.cross(vector(n_atoms_coords[0], n_atoms_coords[1]),
-                                        vector(n_atoms_coords[2], n_atoms_coords[0]))
-                b = vecangle(ring.normal, amine_normal)
-                # Smallest of two angles, depending on direction of normal
-                a = min(b, 180 - b if not 180 - b < 0 else b)
-                if not a > 30.0:
-                    resnr, restype = whichresnumber(ring.atoms[0]), whichrestype(ring.atoms[0])
-                    reschain = whichchain(ring.atoms[0])
-                    resnr_l, restype_l = whichresnumber(p.orig_atoms[0]), whichrestype(p.orig_atoms[0])
-                    reschain_l = whichchain(p.orig_atoms[0])
-                    contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular',
-                                   restype=restype, resnr=resnr, reschain=reschain,
-                                   restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l,
-                                   protcharged=protcharged)
-                    pairings.append(contact)
-                break
-            resnr = whichresnumber(p.atoms[0]) if protcharged else whichresnumber(ring.atoms[0])
-            resnr_l = whichresnumber(ring.orig_atoms[0]) if protcharged else whichresnumber(p.orig_atoms[0])
-            restype = whichrestype(p.atoms[0]) if protcharged else whichrestype(ring.atoms[0])
-            restype_l = whichrestype(ring.orig_atoms[0]) if protcharged else whichrestype(p.orig_atoms[0])
-            reschain = whichchain(p.atoms[0]) if protcharged else whichchain(ring.atoms[0])
-            reschain_l = whichchain(ring.orig_atoms[0]) if protcharged else whichchain(p.orig_atoms[0])
-            contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular', restype=restype,
-                           resnr=resnr, reschain=reschain, restype_l=restype_l, resnr_l=resnr_l,
-                           reschain_l=reschain_l, protcharged=protcharged)
-            pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def saltbridge(poscenter, negcenter, protispos):
-    """Detect all salt bridges (pliprofiler between centers of positive and negative charge)"""
-    data = namedtuple(
-        'saltbridge', 'positive negative distance protispos resnr restype reschain resnr_l restype_l reschain_l')
-    pairings = []
-    for pc, nc in itertools.product(poscenter, negcenter):
-        if not config.MIN_DIST < euclidean3d(pc.center, nc.center) < config.SALTBRIDGE_DIST_MAX:
-            continue
-        resnr = pc.resnr if protispos else nc.resnr
-        resnr_l = whichresnumber(nc.orig_atoms[0]) if protispos else whichresnumber(pc.orig_atoms[0])
-        restype = pc.restype if protispos else nc.restype
-        restype_l = whichrestype(nc.orig_atoms[0]) if protispos else whichrestype(pc.orig_atoms[0])
-        reschain = pc.reschain if protispos else nc.reschain
-        reschain_l = whichchain(nc.orig_atoms[0]) if protispos else whichchain(pc.orig_atoms[0])
-        contact = data(positive=pc, negative=nc, distance=euclidean3d(pc.center, nc.center), protispos=protispos,
-                       resnr=resnr, restype=restype, reschain=reschain, resnr_l=resnr_l, restype_l=restype_l,
-                       reschain_l=reschain_l)
-        pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def halogen(acceptor, donor):
-    """Detect all halogen bonds of the type Y-O...X-C"""
-    data = namedtuple('halogenbond', 'acc acc_orig_idx don don_orig_idx distance don_angle acc_angle restype '
-                                     'resnr reschain restype_l resnr_l reschain_l donortype acctype sidechain')
-    pairings = []
-    for acc, don in itertools.product(acceptor, donor):
-        dist = euclidean3d(acc.o.coords, don.x.coords)
-        if not config.MIN_DIST < dist < config.HALOGEN_DIST_MAX:
-            continue
-        vec1, vec2 = vector(acc.o.coords, acc.y.coords), vector(acc.o.coords, don.x.coords)
-        vec3, vec4 = vector(don.x.coords, acc.o.coords), vector(don.x.coords, don.c.coords)
-        acc_angle, don_angle = vecangle(vec1, vec2), vecangle(vec3, vec4)
-        is_sidechain_hal = acc.o.OBAtom.GetResidue().GetAtomProperty(acc.o.OBAtom, 8)  # Check if sidechain atom
-        if not config.HALOGEN_ACC_ANGLE - config.HALOGEN_ANGLE_DEV < acc_angle \
-               < config.HALOGEN_ACC_ANGLE + config.HALOGEN_ANGLE_DEV:
-            continue
-        if not config.HALOGEN_DON_ANGLE - config.HALOGEN_ANGLE_DEV < don_angle \
-               < config.HALOGEN_DON_ANGLE + config.HALOGEN_ANGLE_DEV:
-            continue
-        restype, reschain, resnr = whichrestype(acc.o), whichchain(acc.o), whichresnumber(acc.o)
-        restype_l, reschain_l, resnr_l = whichrestype(don.orig_x), whichchain(don.orig_x), whichresnumber(don.orig_x)
-        contact = data(acc=acc, acc_orig_idx=acc.o_orig_idx, don=don, don_orig_idx=don.x_orig_idx,
-                       distance=dist, don_angle=don_angle, acc_angle=acc_angle,
-                       restype=restype, resnr=resnr,
-                       reschain=reschain, restype_l=restype_l,
-                       reschain_l=reschain_l, resnr_l=resnr_l, donortype=don.x.OBAtom.GetType(), acctype=acc.o.type,
-                       sidechain=is_sidechain_hal)
-        pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water):
-    """Find water-bridged hydrogen bonds between ligand and protein. For now only considers bridged of first degree."""
-    data = namedtuple('waterbridge', 'a a_orig_idx atype d d_orig_idx dtype h water water_orig_idx distance_aw '
-                                     'distance_dw d_angle w_angle type resnr restype reschain resnr_l restype_l reschain_l protisdon')
-    pairings = []
-    # First find all acceptor-water pairs with distance within d
-    # and all donor-water pairs with distance within d and angle greater theta
-    lig_aw, prot_aw, lig_dw, prot_hw = [], [], [], []
-    for w in water:
-        for acc1 in lig_hba:
-            dist = euclidean3d(acc1.a.coords, w.oxy.coords)
-            if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST:
-                lig_aw.append((acc1, w, dist))
-        for acc2 in bs_hba:
-            dist = euclidean3d(acc2.a.coords, w.oxy.coords)
-            if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST:
-                prot_aw.append((acc2, w, dist))
-        for don1 in lig_hbd:
-            dist = euclidean3d(don1.d.coords, w.oxy.coords)
-            d_angle = vecangle(vector(don1.h.coords, don1.d.coords), vector(don1.h.coords, w.oxy.coords))
-            if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST \
-                    and d_angle > config.WATER_BRIDGE_THETA_MIN:
-                lig_dw.append((don1, w, dist, d_angle))
-        for don2 in bs_hbd:
-            dist = euclidean3d(don2.d.coords, w.oxy.coords)
-            d_angle = vecangle(vector(don2.h.coords, don2.d.coords), vector(don2.h.coords, w.oxy.coords))
-            if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST \
-                    and d_angle > config.WATER_BRIDGE_THETA_MIN:
-                prot_hw.append((don2, w, dist, d_angle))
-
-    for l, p in itertools.product(lig_aw, prot_hw):
-        acc, wl, distance_aw = l
-        don, wd, distance_dw, d_angle = p
-        if not wl.oxy == wd.oxy:
-            continue
-        # Same water molecule and angle within omega
-        w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords))
-        if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX:
-            continue
-        resnr, reschain, restype = whichresnumber(don.d), whichchain(don.d), whichrestype(don.d)
-        resnr_l, reschain_l, restype_l = whichresnumber(acc.a_orig_atom), whichchain(
-            acc.a_orig_atom), whichrestype(acc.a_orig_atom)
-        contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, atype=acc.a.type, d=don.d, d_orig_idx=don.d_orig_idx,
-                       dtype=don.d.type, h=don.h, water=wl.oxy, water_orig_idx=wl.oxy_orig_idx,
-                       distance_aw=distance_aw, distance_dw=distance_dw, d_angle=d_angle, w_angle=w_angle,
-                       type='first_deg', resnr=resnr, restype=restype,
-                       reschain=reschain, restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l, protisdon=True)
-        pairings.append(contact)
-    for p, l in itertools.product(prot_aw, lig_dw):
-        acc, wl, distance_aw = p
-        don, wd, distance_dw, d_angle = l
-        if not wl.oxy == wd.oxy:
-            continue
-        # Same water molecule and angle within omega
-        w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords))
-        if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX:
-            continue
-        resnr, reschain, restype = whichresnumber(acc.a), whichchain(acc.a), whichrestype(acc.a)
-        resnr_l, reschain_l, restype_l = whichresnumber(don.d_orig_atom), whichchain(
-            don.d_orig_atom), whichrestype(don.d_orig_atom)
-        contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, atype=acc.a.type, d=don.d, d_orig_idx=don.d_orig_idx,
-                       dtype=don.d.type, h=don.h, water=wl.oxy, water_orig_idx=wl.oxy_orig_idx,
-                       distance_aw=distance_aw, distance_dw=distance_dw,
-                       d_angle=d_angle, w_angle=w_angle, type='first_deg', resnr=resnr,
-                       restype=restype, reschain=reschain,
-                       restype_l=restype_l, reschain_l=reschain_l, resnr_l=resnr_l, protisdon=False)
-        pairings.append(contact)
-    return filter_contacts(pairings)
-
-
-def metal_complexation(metals, metal_binding_lig, metal_binding_bs):
-    """Find all metal complexes between metals and appropriate groups in both protein and ligand, as well as water"""
-    data = namedtuple('metal_complex', 'metal metal_orig_idx metal_type target target_orig_idx target_type '
-                                       'coordination_num distance resnr restype '
-                                       'reschain  restype_l reschain_l resnr_l location rms, geometry num_partners complexnum')
-    pairings_dict = {}
-    pairings = []
-    # #@todo Refactor
-    metal_to_id = {}
-    metal_to_orig_atom = {}
-    for metal, target in itertools.product(metals, metal_binding_lig + metal_binding_bs):
-        distance = euclidean3d(metal.m.coords, target.atom.coords)
-        if not distance < config.METAL_DIST_MAX:
-            continue
-        if metal.m not in pairings_dict:
-            pairings_dict[metal.m] = [(target, distance), ]
-            metal_to_id[metal.m] = metal.m_orig_idx
-            metal_to_orig_atom[metal.m] = metal.orig_m
-        else:
-            pairings_dict[metal.m].append((target, distance))
-    for cnum, metal in enumerate(pairings_dict):
-        rms = 0.0
-        excluded = []
-        # cnum +1 being the complex number
-        contact_pairs = pairings_dict[metal]
-        num_targets = len(contact_pairs)
-        vectors_dict = defaultdict(list)
-        for contact_pair in contact_pairs:
-            target, distance = contact_pair
-            vectors_dict[target.atom.idx].append(vector(metal.coords, target.atom.coords))
-
-        # Listing of coordination numbers and their geometries
-        configs = {2: ['linear', ],
-                   3: ['trigonal.planar', 'trigonal.pyramidal'],
-                   4: ['tetrahedral', 'square.planar'],
-                   5: ['trigonal.bipyramidal', 'square.pyramidal'],
-                   6: ['octahedral', ]}
-
-        # Angle signatures for each geometry (as seen from each target atom)
-        ideal_angles = {'linear': [[180.0]] * 2,
-                        'trigonal.planar': [[120.0, 120.0]] * 3,
-                        'trigonal.pyramidal': [[109.5, 109.5]] * 3,
-                        'tetrahedral': [[109.5, 109.5, 109.5, 109.5]] * 4,
-                        'square.planar': [[90.0, 90.0, 90.0, 90.0]] * 4,
-                        'trigonal.bipyramidal': [[120.0, 120.0, 90.0, 90.0]] * 3 + [[90.0, 90.0, 90.0, 180.0]] * 2,
-                        'square.pyramidal': [[90.0, 90.0, 90.0, 180.0]] * 4 + [[90.0, 90.0, 90.0, 90.0]],
-                        'octahedral': [[90.0, 90.0, 90.0, 90.0, 180.0]] * 6}
-        angles_dict = {}
-
-        for target in vectors_dict:
-            cur_vector = vectors_dict[target]
-            other_vectors = []
-            for t in vectors_dict:
-                if not t == target:
-                    [other_vectors.append(x) for x in vectors_dict[t]]
-            angles = [vecangle(pair[0], pair[1]) for pair in itertools.product(cur_vector, other_vectors)]
-            angles_dict[target] = angles
-
-        all_total = []  # Record fit information for each geometry tested
-        gdata = namedtuple('gdata', 'geometry rms coordination excluded diff_targets')  # Geometry Data
-        # Can't specify geometry with only one target
-        if num_targets == 1:
-            final_geom = 'NA'
-            final_coo = 1
-            excluded = []
-            rms = 0.0
-        else:
-            for coo in sorted(configs, reverse=True):  # Start with highest coordination number
-                geometries = configs[coo]
-                for geometry in geometries:
-                    signature = ideal_angles[geometry]  # Set of ideal angles for geometry, from each perspective
-                    geometry_total = 0
-                    geometry_scores = []  # All scores for one geometry (from all subsignatures)
-                    used_up_targets = []  # Use each target just once for a subsignature
-                    not_used = []
-                    coo_diff = num_targets - coo  # How many more observed targets are there?
-
-                    # Find best match for each subsignature
-                    for subsignature in signature:  # Ideal angles from one perspective
-                        best_target = None  # There's one best-matching target for each subsignature
-                        best_target_score = 999
-
-                        for k, target in enumerate(angles_dict):
-                            if target not in used_up_targets:
-                                observed_angles = angles_dict[target]  # Observed angles from perspective of one target
-                                single_target_scores = []
-                                used_up_observed_angles = []
-                                for i, ideal_angle in enumerate(subsignature):
-                                    # For each angle in the signature, find the best-matching observed angle
-                                    best_match = None
-                                    best_match_diff = 999
-                                    for j, observed_angle in enumerate(observed_angles):
-                                        if j not in used_up_observed_angles:
-                                            diff = abs(ideal_angle - observed_angle)
-                                            if diff < best_match_diff:
-                                                best_match_diff = diff
-                                                best_match = j
-                                    if best_match is not None:
-                                        used_up_observed_angles.append(best_match)
-                                        single_target_scores.append(best_match_diff)
-                                # Calculate RMS for target angles
-                                target_total = sum([x ** 2 for x in single_target_scores]) ** 0.5  # Tot. score targ/sig
-                                if target_total < best_target_score:
-                                    best_target_score = target_total
-                                    best_target = target
-
-                        used_up_targets.append(best_target)
-                        geometry_scores.append(best_target_score)
-                        # Total score is mean of RMS values
-                        geometry_total = np.mean(geometry_scores)
-                    # Record the targets not used for excluding them when deciding for a final geometry
-                    [not_used.append(target) for target in angles_dict if target not in used_up_targets]
-                    all_total.append(gdata(geometry=geometry, rms=geometry_total, coordination=coo,
-                                           excluded=not_used, diff_targets=coo_diff))
-
-        # Make a decision here. Starting with the geometry with lowest difference in ideal and observed partners ...
-        # Check if the difference between the RMS to the next best solution is not larger than 0.5
-        if not num_targets == 1:  # Can't decide for any geoemtry in that case
-            all_total = sorted(all_total, key=lambda x: abs(x.diff_targets))
-            for i, total in enumerate(all_total):
-                next_total = all_total[i + 1]
-                this_rms, next_rms = total.rms, next_total.rms
-                diff_to_next = next_rms - this_rms
-                if diff_to_next > 0.5:
-                    final_geom, final_coo, rms, excluded = total.geometry, total.coordination, total.rms, total.excluded
-                    break
-                elif next_total.rms < 3.5:
-                    final_geom, final_coo, = next_total.geometry, next_total.coordination
-                    rms, excluded = next_total.rms, next_total.excluded
-                    break
-                elif i == len(all_total) - 2:
-                    final_geom, final_coo, rms, excluded = "NA", "NA", float('nan'), []
-                    break
-
-        # Record all contact pairing, excluding those with targets superfluous for chosen geometry
-        only_water = set([x[0].location for x in contact_pairs]) == {'water'}
-        if not only_water:  # No complex if just with water as targets
-            logger.info(f'metal ion {metal.type} complexed with {final_geom} geometry (coo. number {final_coo}/ {num_targets} observed)')
-            for contact_pair in contact_pairs:
-                target, distance = contact_pair
-                if target.atom.idx not in excluded:
-                    metal_orig_atom = metal_to_orig_atom[metal]
-                    restype_l, reschain_l, resnr_l = whichrestype(metal_orig_atom), whichchain(
-                        metal_orig_atom), whichresnumber(metal_orig_atom)
-                    contact = data(metal=metal, metal_orig_idx=metal_to_id[metal], metal_type=metal.type,
-                                   target=target, target_orig_idx=target.atom_orig_idx, target_type=target.type,
-                                   coordination_num=final_coo, distance=distance, resnr=target.resnr,
-                                   restype=target.restype, reschain=target.reschain, location=target.location,
-                                   rms=rms, geometry=final_geom, num_partners=num_targets, complexnum=cnum + 1,
-                                   resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l)
-                    pairings.append(contact)
-    return filter_contacts(pairings)