Merge pull request #276 from rwxayheee/reduce_prepare_overhead

Remove from class Ring the unused attributes to reduce overhead in molecule preparation

meeko/macrocycle.py

@@ -83,9 +83,7 @@ def collect_rings(self, setup):
             setup.rings.keys()
         ):  # ring_id are the atom indices in each ring
             size = len(ring_id)
-            if setup.rings[ring_id].is_aromatic:
-                rigid_rings.append(ring_id)
-            elif size < self._min_ring_size:
+            if size < self._min_ring_size:
                 rigid_rings.append(ring_id)
                 # do not add rings > _max_ring_size to rigid_rings
                 # because bonds in rigid rings will not be breakable

meeko/molsetup.py

@@ -48,9 +48,6 @@
 
 DEFAULT_BOND_ROTATABLE = False
 
-DEFAULT_RING_CORNER_FLIP = False
-DEFAULT_RING_GRAPH = []
-DEFAULT_RING_IS_AROMATIC = False
 DEFAULT_RING_CLOSURE_BONDS_REMOVED = []
 DEFAULT_RING_CLOSURE_PSEUDOS_BY_ATOM = defaultdict
 # endregion
@@ -348,9 +345,6 @@ def from_json(obj: dict):
 @dataclass
 class Ring:
     ring_id: tuple
-    corner_flip: bool = DEFAULT_RING_CORNER_FLIP
-    graph: dict = field(default_factory=list)
-    is_aromatic: bool = DEFAULT_RING_IS_AROMATIC
 
     @staticmethod
     def from_json(obj: dict):
@@ -374,16 +368,13 @@ def from_json(obj: dict):
             return obj
 
         # Check that all the keys we expect are in the object dictionary as a safety measure
-        expected_json_keys = {"ring_id", "corner_flip", "graph", "is_aromatic"}
+        expected_json_keys = {"ring_id"}
         if set(obj.keys()) != expected_json_keys:
             return obj
 
         # Constructs a Ring object from the provided keys.
         ring_id = string_to_tuple(obj["ring_id"], int)
-        corner_flip = obj["corner_flip"]
-        graph = obj["graph"]
-        is_aromatic = obj["is_aromatic"]
-        output_ring = Ring(ring_id, corner_flip, graph, is_aromatic)
+        output_ring = Ring(ring_id)
         return output_ring
 
 
@@ -1524,8 +1515,6 @@ class RDKitMoleculeSetup(MoleculeSetup, MoleculeSetupExternalToolkit):
         a mapping from tuples of atom indices to dihedral labels
     atom_to_ring_id: dict()
         mapping of atom index to ring id of each atom belonging to the ring
-    ring_corners: dict()
-        unclear what this is a mapping of, but is used to store corner flexibility for the rings
     rmsd_symmetry_indices: tuple
         Tuples of the indices of the molecule's atoms that match a substructure query. needs info.
 
@@ -1543,7 +1532,6 @@ def __init__(self, name: str = None, is_sidechain: bool = False):
         self.dihedral_partaking_atoms: dict = {}
         self.dihedral_labels: dict = {}
         self.atom_to_ring_id = {}
-        self.ring_corners = {}
         self.rmsd_symmetry_indices = ()
 
     def copy(self):
@@ -1622,7 +1610,7 @@ def from_mol(
         molsetup.init_atom(compute_gasteiger_charges, read_charges_from_prop, coords)
         molsetup.init_bond()
         molsetup.perceive_rings(keep_chorded_rings, keep_equivalent_rings)
-        molsetup.rmsd_symmetry_indices = cls.get_symmetries_for_rmsd(mol)
+        # molsetup.rmsd_symmetry_indices = cls.get_symmetries_for_rmsd(mol)
 
         # to store sets of coordinates, e.g. docked poses, as dictionaries indexed by
         # the atom index, because not all atoms need to have new coordinates specified
@@ -1974,13 +1962,6 @@ def perceive_rings(self, keep_chorded_rings: bool, keep_equivalent_rings: bool):
         rings = hjk_ring_detection.scan(keep_chorded_rings, keep_equivalent_rings)
         for ring_atom_indices in rings:
             ring_to_add = Ring(ring_atom_indices)
-            if self._is_ring_aromatic(ring_atom_indices):
-                ring_to_add.is_aromatic = True
-            for atom_idx in ring_atom_indices:
-                # TODO: add it to some sort of atom to ring id tracking thing -> add to atom data structure
-                ring_to_add.graph = self._recursive_graph_walk(
-                    atom_idx, collected=[], exclude=list(ring_atom_indices)
-                )
             self.rings[ring_atom_indices] = ring_to_add
         return
 
@@ -2152,7 +2133,6 @@ def from_json(obj):
             "dihedral_partaking_atoms",
             "dihedral_labels",
             "atom_to_ring_id",
-            "ring_corners",
             "rmsd_symmetry_indices",
         }
         for key in expected_molsetup_keys:
@@ -2179,7 +2159,6 @@ def from_json(obj):
             int(k): [string_to_tuple(t) for t in v]
             for k, v in obj["atom_to_ring_id"].items()
         }
-        rdkit_molsetup.ring_corners = obj["ring_corners"]
         rdkit_molsetup.rmsd_symmetry_indices = [
             string_to_tuple(v) for v in obj["rmsd_symmetry_indices"]
         ]
@@ -2280,9 +2259,6 @@ def default(self, obj):
         if isinstance(obj, Ring):
             return {
                 "ring_id": tuple_to_string(obj.ring_id),
-                "corner_flip": obj.corner_flip,
-                "graph": obj.graph,
-                "is_aromatic": obj.is_aromatic,
             }
         return json.JSONEncoder.default(self, obj)
 
@@ -2387,7 +2363,6 @@ def default(self, obj):
             output_dict["dihedral_partaking_atoms"] = {tuple_to_string(k): v for k,v in obj.dihedral_partaking_atoms.items()}
             output_dict["dihedral_labels"] = {tuple_to_string(k): v for k,v in obj.dihedral_labels.items()}
             output_dict["atom_to_ring_id"] = obj.atom_to_ring_id
-            output_dict["ring_corners"] = obj.ring_corners
             output_dict["rmsd_symmetry_indices"] = obj.rmsd_symmetry_indices
         # If nothing is in the dict, then none of the possible object types for this encoder matched and we should
         # return the default JSON encoder.

test/json_serialization_test.py

@@ -525,13 +525,6 @@ def check_ring_equality(decoded_obj: Ring, starting_obj: Ring):
     """
     assert isinstance(decoded_obj.ring_id, tuple)
     assert decoded_obj.ring_id == starting_obj.ring_id
-    assert isinstance(decoded_obj.corner_flip, bool)
-    assert decoded_obj.corner_flip == starting_obj.corner_flip
-    assert len(decoded_obj.graph) == len(starting_obj.graph)
-    for idx, val in enumerate(starting_obj.graph):
-        assert decoded_obj.graph[idx] == val
-    assert isinstance(decoded_obj.is_aromatic, bool)
-    assert decoded_obj.is_aromatic == starting_obj.is_aromatic
     return
 
 

test/macrocycle_test.py

@@ -62,6 +62,30 @@ def test_external_ring_closure():
             mindist = min(mindist, abs(x - xcoord))
         assert(mindist < 1e-3) # loose matching
 
+def test_aromatic_rings_unbroken():
+    write_pdbqt = PDBQTWriterLegacy.write_string
+    mol = Chem.MolFromMolFile(filenames["macrocycle3"], removeHs=False)
+    mk_prep = MoleculePreparation(min_ring_size=5)
+    setups = mk_prep(mol)
+    assert len(setups) == 1
+    setup = setups[0]
+    nr_rot = sum([b.rotatable for _, b in setup.bond_info.items()])
+    assert nr_rot == 9
+    pdbqt_string, is_ok, error_msg = write_pdbqt(setup)
+    assert pdbqt_string.count("ENDBRANCH") == 8
+    # now we will check that a 6-member ring that is not aromatic does break
+    # when we set the min ring size to six or less
+    mol2 = Chem.MolFromMolFile(filenames["macrocycle5"], removeHs=False)
+    setups_broken = mk_prep(mol2)
+    mk_prep = MoleculePreparation(min_ring_size=7) 
+    setups_unbroken = mk_prep(mol2)
+    assert len(setups_broken) == 1
+    assert len(setups_unbroken) == 1
+    pdbqt_broken, is_ok, err = write_pdbqt(setups_broken[0])
+    pdbqt_unbroken, is_ok, err = write_pdbqt(setups_unbroken[0])
+    assert pdbqt_broken.count("ENDBRANCH") > pdbqt_unbroken.count("ENDBRANCH")
+
+
 def run(molname):
     filename = filenames[molname]
     mol = Chem.MolFromMolFile(filename, removeHs=False)
@@ -76,3 +100,4 @@ def run(molname):
 def test_all():
     for molname in num_cycle_breaks:
         run(molname)
+