variable bond length based on element types

src/DockQ/DockQ.py

@@ -18,13 +18,15 @@
 try:
     from .operations import residue_distances, get_fnat_stats
     from .parsers import PDBParser, MMCIFParser
+    from .constants import *
 except ImportError:
     warnings.warn(
         """WARNING: It looks like cython is not working,
          falling back on native python. This will make DockQ slower"""
     )
     from operations_nocy import residue_distances, get_fnat_stats
     from parsers import PDBParser, MMCIFParser
+    from constants import *
 
 
 def parse_args():
@@ -197,8 +199,11 @@ def calc_sym_corrected_lrmsd(
         get_atoms_per_residue((aligned_ref_receptor, aligned_sample_receptor), what="receptor")
     )
 
-    sample_ligand_atoms_ids = [atom.id for atom in sample_ligand.get_atoms() if atom.id != "H"]
-    ref_ligand_atoms_ids = [atom.id for atom in ref_ligand.get_atoms() if atom.id != "H"]
+    sample_ligand_atoms_ids = [atom.id for atom in sample_ligand.get_atoms()]
+    sample_ligand_atoms_ele = [atom.element for atom in sample_ligand.get_atoms()]
+
+    ref_ligand_atoms_ids = [atom.id for atom in ref_ligand.get_atoms()]
+    ref_ligand_atoms_ele = [atom.element for atom in ref_ligand.get_atoms()]
 
     sample_ligand_atoms = np.array([atom.coord for atom in sample_ligand.get_atoms() if atom.id in ref_ligand_atoms_ids])
     ref_ligand_atoms = np.array([atom.coord for atom in ref_ligand.get_atoms() if atom.id in sample_ligand_atoms_ids])
@@ -211,13 +216,13 @@ def calc_sym_corrected_lrmsd(
 
     sample_rotated_ligand_atoms = np.dot(sample_ligand_atoms, rot) + tran
 
-    sample_graph = create_graph(sample_ligand_atoms)
-    ref_graph = create_graph(ref_ligand_atoms)
+    sample_graph = create_graph(sample_ligand_atoms, sample_ligand_atoms_ele)
+    ref_graph = create_graph(ref_ligand_atoms, ref_ligand_atoms_ele)
 
     min_Lrms = float("inf")
     best_mapping = None
+    print(sample_graph, ref_graph)
     for isomorphism in nx.vf2pp_all_isomorphisms(sample_graph, ref_graph):
-
         model_i = list(isomorphism.keys())
         native_i = list(isomorphism.values())
 
@@ -593,14 +598,20 @@ def run_on_chains(
     return info
 
 
-def create_graph(atom_list, threshold=2.0):
+def create_graph(atom_list, atom_ids):
     import networkx as nx
     G = nx.Graph()
+
     for i, atom_i in enumerate(atom_list):
+        cr_i = covalent_radius[atom_ids[i]]
         for j, atom_j in enumerate(atom_list):
+            cr_j = covalent_radius[atom_ids[j]]
             distance = np.linalg.norm(atom_i - atom_j)
+            threshold = (cr_i + cr_j + bond_tolerance) if i != j else 0.0
+            #print(atom_ids[i], atom_ids[j], threshold)
             if distance < threshold:  # Adjust threshold as needed
                 G.add_edge(i, j)
+
     return G
 
 

src/DockQ/constants.py

@@ -0,0 +1,124 @@
+from typing import Dict
+
+bond_tolerance: float = 0.4
+
+covalent_radius: Dict = {
+    "H": 0.31,
+    "HE": 0.28,
+    "LI": 1.28,
+    "BE": 0.96,
+    "B": 0.84,
+    "C": 0.76,
+    "N": 0.71,
+    "O": 0.66,
+    "F": 0.57,
+    "NE": 0.58,
+    "NA": 1.66,
+    "MG": 1.41,
+    "AL": 1.21,
+    "SI": 1.11,
+    "P": 1.07,
+    "S": 1.05,
+    "CL": 1.02,
+    "AR": 1.06,
+    "K": 2.03,
+    "CA": 1.76,
+    "SC": 1.7,
+    "TI": 1.6,
+    "V": 1.53,
+    "CR": 1.39,
+    "MN": 1.5,
+    "FE": 1.42,
+    "CO": 1.38,
+    "NI": 1.24,
+    "CU": 1.32,
+    "ZN": 1.22,
+    "GA": 1.22,
+    "GE": 1.2,
+    "AS": 1.19,
+    "SE": 1.2,
+    "BR": 1.2,
+    "KR": 1.16,
+    "RB": 2.2,
+    "SR": 1.95,
+    "Y": 1.9,
+    "ZR": 1.75,
+    "NB": 1.64,
+    "MO": 1.54,
+    "TC": 1.47,
+    "RU": 1.46,
+    "RH": 1.42,
+    "PD": 1.39,
+    "AG": 1.45,
+    "CD": 1.44,
+    "IN": 1.42,
+    "SN": 1.39,
+    "SB": 1.39,
+    "TE": 1.38,
+    "I": 1.39,
+    "XE": 1.4,
+    "CS": 2.44,
+    "BA": 2.15,
+    "LA": 2.07,
+    "CE": 2.04,
+    "PR": 2.03,
+    "ND": 2.01,
+    "PM": 1.99,
+    "SM": 1.98,
+    "EU": 1.98,
+    "GD": 1.96,
+    "TB": 1.94,
+    "DY": 1.92,
+    "HO": 1.92,
+    "ER": 1.89,
+    "TM": 1.9,
+    "YB": 1.87,
+    "LU": 1.87,
+    "HF": 1.75,
+    "TA": 1.7,
+    "W": 1.62,
+    "RE": 1.51,
+    "OS": 1.44,
+    "IR": 1.41,
+    "PT": 1.36,
+    "AU": 1.36,
+    "HG": 1.32,
+    "TL": 1.45,
+    "PB": 1.46,
+    "BI": 1.48,
+    "PO": 1.4,
+    "AT": 1.5,
+    "RN": 1.5,
+    "FR": 2.6,
+    "RA": 2.21,
+    "AC": 2.15,
+    "TH": 2.06,
+    "PA": 2.0,
+    "U": 1.96,
+    "NP": 1.9,
+    "PU": 1.87,
+    "AM": 1.8,
+    "CM": 1.69,
+    "BK": 2.0,
+    "CF": 2.0,
+    "ES": 2.0,
+    "FM": 2.0,
+    "MD": 2.0,
+    "NO": 2.0,
+    "LR": 2.0,
+    "RF": 2.0,
+    "DB": 2.0,
+    "SG": 2.0,
+    "BH": 2.0,
+    "HS": 2.0,
+    "MT": 2.0,
+    "DS": 2.0,
+    "RG": 2.0,
+    "CN": 2.0,
+    "UUT": 2.0,
+    "UUQ": 2.0,
+    "UUP": 2.0,
+    "UUH": 2.0,
+    "UUS": 2.0,
+    "UUO": 2.0,
+}