Milestone: Add feature to lock-in phenyl rings

Author: bertrand-caron

Makefile

@@ -1,7 +1,7 @@
 PYTHON_BIN_DIR = /usr/local/python35/bin
 
 test_tautomers:
-	-@rm graphs/benzene* graphs/ethanal* graphs/propadiene*.pdf graphs/cyclooctyne*.pdf graphs/methylimidazole*.pdf graphs/porphyrin*.pdf
+	-@rm graphs/*tautomer*.pdf
 	python3 test_tautomers.py
 .PHONY: test_tautomers.py
 

helpers/molecule.py

@@ -12,7 +12,7 @@
 from fragment_capping.helpers.types_helpers import Fragment, ATB_Molid, Atom, FRAGMENT_CAPPING_DIR, Bond, ATOM_INDEX, MIN, MAX
 from fragment_capping.helpers.parameters import FULL_VALENCES, POSSIBLE_CHARGES,  Capping_Strategy, possible_bond_order_for_atom_pair, coordinates_n_angstroms_away_from, possible_charge_for_atom, ALL_ELEMENTS, electronegativity_spread, ELECTRONEGATIVITIES, VALENCE_ELECTRONS, MIN_ABSOLUTE_CHARGE, MAX_ABSOLUTE_CHARGE, MIN_BOND_ORDER, MAX_BOND_ORDER, MUST_BE_INT, MAX_NONBONDED_ELECTRONS, NO_CAP, ELECTRONS_PER_BOND, ALL_CAPPING_OPTIONS
 from fragment_capping.helpers.babel import energy_minimised_pdb
-from fragment_capping.helpers.rings import bonds_for_ring
+from fragment_capping.helpers.rings import bonds_for_ring, atom_indices_in_phenyl_rings_for
 from fragment_capping.helpers.graphs import unique_molecules
 from fragment_capping.helpers.tautomers import get_all_tautomers, get_all_tautomers_naive
 from fragment_capping.helpers.capping import get_best_capped_molecule_with_ILP, get_best_capped_molecule
@@ -67,6 +67,7 @@ def __init__(
         self.previously_uncapped = set()
         self.aromatic_bonds, self.aromatic_rings = (None, None)
         self.non_bonded_electrons = None
+        self.phenyl_atoms = set()
 
     def atom_desc(self, atom: Atom):
         if self.use_neighbour_valences:
@@ -924,6 +925,12 @@ def renumber_atoms(self) -> None:
                 for bond in self.aromatic_bonds
             }
 
+        if self.phenyl_atoms is not None:
+            self.phenyl_atoms = {
+                remap_atom(atom_index)
+                for atom_index in self.phenyl_atoms
+            }
+
         return None
 
     def remove_atoms_with_predicate(self, predicate: Callable[[Atom], bool]) -> None:
@@ -979,6 +986,9 @@ def remove_atoms_with_predicate(self, predicate: Callable[[Atom], bool]) -> None
         return self.renumber_atoms()
 
     def remove_all_hydrogens(self, mark_all_uncapped: bool = False) -> None:
+        # Store phenyl rings atoms before removing their hydrogens
+        self.phenyl_atoms = atom_indices_in_phenyl_rings_for(self)
+
         self.remove_atoms_with_predicate(
             lambda atom: atom.element == 'H',
         )

helpers/rings.py

@@ -11,7 +11,7 @@ def bonds_for_ring(ring: List[ATOM_INDEX]) -> List[Bond]:
         for atom_indices in zip(ring, ring[1:] + (ring[0],))
     ]
 
-def atoms_in_small_rings_for(molecule: 'Molecule') -> Set[Atom]:
+def atoms_in_small_rings_for(molecule: 'Molecule') -> List[Atom]:
     return [
         molecule.atoms[atom_index]
         for atom_index in reduce(
@@ -27,3 +27,23 @@ def bonds_in_small_rings_for(molecule: 'Molecule') -> Set[Bond]:
         [bonds_for_ring(ring) for ring in molecule.rings() if 0 < len(ring) < SMALL_RING],
         set(),
     )
+
+def atom_indices_in_phenyl_rings_for(molecule: 'Molecule') -> Set[ATOM_INDEX]:
+    return reduce(
+        lambda acc, e: acc | set(e),
+        [
+            ring
+            for ring in molecule.rings()
+            if (
+                len(ring) == 6
+                and
+                all(
+                    len([bond for bond in molecule.bonds if atom_index in bond]) == 3
+                    and
+                    molecule.atoms[atom_index].element == 'C'
+                    for atom_index in ring
+                )
+            )
+        ],
+        set(),
+    )

helpers/tautomers.py

@@ -237,9 +237,19 @@ def get_all_tautomers(
     max_tautomers: Optional[int] = 2000,
     disallow_triple_bond_in_small_rings: bool = True,
     disallow_allenes_in_small_rings: bool = True,
-    use_gurobi: bool = True,
+    lock_phenyl_rings: bool = True,
+    use_gurobi: bool = False,
     debug: Optional[TextIO] = None,
 ) -> 'Molecule':
+    '''
+    Args:
+        ``disallow_triple_bond_in_small_rings``: Disallow triple bonds in small rings (rings with size <= ``SMALL_RING``)
+        ``disallow_allenes_in_small_rings``: Disallow allenes (=C=) in small rings (rings with size <= ``SMALL_RING``)
+        ``lock_phenyl_rings``: Prevent phenyl rings (6 membered rings with all sp3 carbons) from being hydrogenised.
+
+    Returns:
+        A list of tautomers (Molecule).
+    '''
     if len([1 for atom in molecule.atoms.values() if atom.element == 'H']) > 0:
         molecule.remove_all_hydrogens(mark_all_uncapped=True)
 
@@ -254,31 +264,18 @@ def keep_capping_strategy_for_atom(capping_strategy: Capping_Strategy, atom: Ato
         else:
             return min_valence_for(atom) <= neighbour_counts[atom.index] + new_atom_for_capping_strategy(capping_strategy) <= max_valence_for(atom)
 
-    def possible_capping_strategies_for_atom(atom: Atom) -> List[Capping_Strategy]:
-        return [NO_CAP, H_CAP, H2_CAP, H3_CAP] + ([H4_CAP] if False else [])
-#        if debug is not None:
-#            write_to_debug(debug, '')
-#            write_to_debug(debug, atom)
-#            write_to_debug(debug, 'capping_strategy, new_atom_for_capping_strategy(), keep_capping_strategy_for_atom()')
-#            for capping_strategy in ALL_CAPPING_OPTIONS[molecule.atom_desc(atom)]:
-#                write_to_debug(
-#                    debug,
-#                    capping_strategy,
-#                    new_atom_for_capping_strategy(capping_strategy),
-#                    keep_capping_strategy_for_atom(capping_strategy, atom)
-#                )
-#        return [
-#            capping_strategy
-#            for capping_strategy in ALL_CAPPING_OPTIONS[molecule.atom_desc(atom)]
-#            if keep_capping_strategy_for_atom(capping_strategy, atom)
-#        ]
+    def possible_capping_strategies_for_atom(atom: Atom, is_phenyl_atom: bool = False) -> List[Capping_Strategy]:
+        if is_phenyl_atom and lock_phenyl_rings:
+            return [NO_CAP, H_CAP]
+        else:
+            return [NO_CAP, H_CAP, H2_CAP, H3_CAP] + ([H4_CAP] if False else [])
 
     atoms_need_capping = [atom for atom in molecule.sorted_atoms() if not atom.capped]
 
     write_to_debug(
         debug,
         [
-            possible_capping_strategies_for_atom(atom)
+            possible_capping_strategies_for_atom(atom, is_phenyl_atom=(atom.index in molecule.phenyl_atoms))
             for atom in atoms_need_capping
         ],
     )
@@ -291,7 +288,7 @@ def possible_capping_strategies_for_atom(atom: Atom) -> List[Capping_Strategy]:
     capping_atoms_for = {}
     new_bonds_sets = {}
     for uncapped_atom in atoms_need_capping:
-        possible_capping_strategies = possible_capping_strategies_for_atom(uncapped_atom)
+        possible_capping_strategies = possible_capping_strategies_for_atom(uncapped_atom, is_phenyl_atom=(uncapped_atom.index in molecule.phenyl_atoms))
         if len(possible_capping_strategies) == 0 or len(possible_capping_strategies) == 1 and possible_capping_strategies[0] == NO_CAP:
             pass
         else:

test_tautomers.py

@@ -5,6 +5,35 @@
 from fragment_capping.helpers.molecule import Molecule, molecule_from_pdb_str
 
 class Test_Capping(unittest.TestCase):
+    def test_warfarin(self, use_ILP: bool = True) -> None:
+        with open('pdbs/warfarin.pdb') as fh:
+            input_molecule = molecule_from_pdb_str(
+                fh.read(),
+                name='warfarin',
+            )
+        input_molecule.remove_all_hydrogens(mark_all_uncapped=True)
+        input_molecule.write_graph(
+            'input',
+            output_size=(1200, 1200),
+            graph_kwargs={'include_atom_index': False},
+        )
+        tautomers = input_molecule.get_all_tautomers(
+            net_charge=0,
+            total_number_hydrogens=16,
+            enforce_octet_rule=True,
+            allow_radicals=False,
+        )
+
+        assert len(tautomers) == 52, tautomers
+
+        for (n, molecule) in enumerate(tautomers, start=1):
+            molecule.name = input_molecule.name
+            molecule.write_graph(
+                '_tautomer_{0}'.format(n),
+                output_size=(1200, 1200),
+                graph_kwargs={'include_atom_index': False},
+            )
+
     def test_porphyrin(self, use_ILP: bool = True) -> None:
         with open('pdbs/tetraphenylporphyrin.pdb') as fh:
             input_molecule = molecule_from_pdb_str(
@@ -24,7 +53,7 @@ def test_porphyrin(self, use_ILP: bool = True) -> None:
             allow_radicals=False,
         )
 
-        assert len(tautomers) == 46, tautomers
+        assert len(tautomers) == 21, tautomers
 
         for (n, molecule) in enumerate(tautomers, start=1):
             molecule.name = input_molecule.name
@@ -178,4 +207,4 @@ def test_cyclooctyne(self):
 
 if __name__ == '__main__':
     print('Note: This test is very long (~15 minutes)')
-    unittest.main(warnings='ignore')
+    unittest.main(warnings='ignore', verbosity=2)