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)