AutoREACTER v0.3: Functional group detection refactor and reaction library expansion

AutoREACTER/__init__.py

@@ -3,7 +3,7 @@
 
 AutoREACTER is a tool for automated reaction-based molecular system generation.
 """
-__version__ = "0.2.3"
+__version__ = "0.3"
 
 __title__ = "AutoREACTER"
 __author__ = "Janitha Mahanthe, Jacob Gissinger"

AutoREACTER/arx_cli.py

@@ -14,6 +14,7 @@
 from contextlib import contextmanager
 import os
 from pathlib import Path
+import shutil
 import sys
 import threading
 from PIL import Image
@@ -30,6 +31,9 @@
 from AutoREACTER.reaction_preparation.ff_wrapper.REACTER_files_builder import REACTERFilesBuilder
 from AutoREACTER.sim_setup.simulation_setup import SimulationSetupManager
 
+class NoReactionGenerated(Exception):
+    """Custom exception raised when no reaction is generated in the pipeline.""" 
+    pass
 
 class ErrorHandler:
     """
@@ -98,7 +102,8 @@ def __init__(self, input: Path) -> None:
         self.img_dir = self.session.images_dir
         # with open(self.session.output_dir / "AutoREACTER.log", 'w') as f:
         #     f.write("--- Starting AutoREACTER Session ---\n")
-
+        # Save a copy of the input JSON to the output directory
+        self._save_input_json(abs_path)
         # Save an initial grid image of all monomers
         self._save_rdkit_img(
             InputParser().initial_molecules_image_grid(self.session),
@@ -314,6 +319,10 @@ def process(self):
     # ------------------------------------------------------------------
     # Internal helpers – lazy detection & image saving
     # ------------------------------------------------------------------
+    def _save_input_json(self, abs_path: Path): 
+            destination_file = "input.json"
+            destination_path = self.session.output_dir / destination_file
+            shutil.copy(abs_path, destination_path)
 
     def _ensure_fg_detected(self):
         """
@@ -384,7 +393,10 @@ def _save_rdkit_img(self, img, path: Path, is_non_reactant: bool = False):
         if img is None:
             if is_non_reactant:
                 return
-            raise ValueError("No image was generated. Cannot save molecule image.")
+            raise NoReactionGenerated(
+                "No reaction was generated. This is an error from AutoREACTER. "
+                "Please file an issue on https://github.com/NanoCIPHER-Lab/AutoREACTER/issues to improve the software."
+            )
 
         # Case 1: PIL image
         if hasattr(img, "save"):

AutoREACTER/detectors/functional_groups_detector.py

@@ -1,5 +1,6 @@
 from __future__ import annotations  
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, List
+
 """
 * Monomer Functionality Detection Module
 --------------------------------------
@@ -111,12 +112,14 @@
 from AutoREACTER.input_parser import MonomerEntry
 
 # Conditional import for FunctionalGroupsLibrary to support both installed and local usage.
-from .functional_groups_library import FunctionalGroupsLibrary
+from AutoREACTER.detectors.functional_groups_library.registry import FunctionalGroupsLibrary
 
 logger = logging.getLogger(__name__)  # Module-level logger for future diagnostics.
 if TYPE_CHECKING:
     from AutoREACTER.session import Session
-
+    from AutoREACTER.reaction_preparation.reaction_processor.reaction_progression import (
+        MonomerRoleforIndexBasedFGDetection,
+    )
 
 @dataclass(slots=True)
 class FunctionalGroupInfo:
@@ -127,19 +130,23 @@ class FunctionalGroupInfo:
         functionality_type (str): Type of functionality (e.g., 'vinyl', 'mono', 'di_identical', 'di_different').
         fg_name (str): Name of the functional group (e.g., 'acrylate').
         fg_smarts_1 (str): Primary SMARTS pattern for matching.
+        fg_1_indexes (Optional[Tuple[int, ...]]): Atom indices for matches of fg_smarts_1.
         fg_count_1 (int): Number of matches for fg_smarts_1.
         fg_smarts_2 (Optional[str]): Secondary SMARTS pattern (for 'di_different' types).
+        fg_2_indexes (Optional[Tuple[int, ...]]): Atom indices for matches of fg_smarts_2.
         fg_count_2 (Optional[int]): Number of matches for fg_smarts_2.
     """
     functionality_type: str
     fg_name: str
     fg_smarts_1: str
     fg_count_1: int
+    fg_1_indexes: Optional[Tuple[Tuple[int, ...], ...]] = None
     fg_smarts_2: Optional[str] = None
     fg_count_2: Optional[int] = None
+    fg_2_indexes: Optional[Tuple[Tuple[int, ...], ...]] = None
 
 
-@dataclass(slots=True, frozen=True)
+@dataclass(slots=True)
 class MonomerRole:
     """
     Immutable dataclass representing a monomer with its detected functional groups.
@@ -152,6 +159,10 @@ class MonomerRole:
     smiles: str
     name: str
     functionalities: Tuple[FunctionalGroupInfo, ...]  # Tuple of detected functionalities for the monomer
+    rdkit_mol: Optional[rdchem.Mol] = None  # Optional RDKit molecule object for the monomer
+    indexes_in_template: List[int] = None  # Optional list of atom indices in the template
+    is_monomer: bool = False  # Flag indicating if the monomer is eligible for polymerization
+    is_looped: bool = False  
 
 @dataclass(slots=True)
 class FunctionalGroupVisualization:
@@ -349,6 +360,7 @@ def functional_groups_detector(
                         smiles=smiles,
                         name=monomer.name,
                         functionalities=tuple(detected_functionalities),
+                        is_monomer=True
                     )
                 )
 
@@ -401,6 +413,140 @@ def _functional_groups_detector_for_visualization(
                     )
                 )
         return monomer_roles_visualization
+
+    def _detect_functional_groups_by_index(
+        self,
+        mol: Chem.Mol,
+        smarts: str,
+        atom_indices: list[int],
+    ) -> bool:
+        """Return True when any SMARTS match overlaps the supplied atom indices."""
+        target_indices = set(atom_indices)
+
+        patt = Chem.MolFromSmarts(smarts)
+        if patt is None:
+            logger.warning("Invalid SMARTS pattern: %s", smarts)
+            return False
+
+        matches = mol.GetSubstructMatches(patt, uniquify=True)
+        return any(target_indices.intersection(match) for match in matches)
+
+    def index_based_functional_groups_detector(
+        self,
+        monomer_roles_in: list[MonomerRoleforIndexBasedFGDetection],
+    ) -> list[MonomerRole] | bool:
+        """
+        Detect functional groups across a list of monomers and categorize them into roles,
+        restricted to a given set of atom indices per monomer.
+
+        Iterates over predefined monomer_types, matches each against the monomer's
+        rdkit_mol, and keeps only matches that overlap with the monomer's
+        `indexes_in_template`. Prints matches for debugging/user feedback.
+
+        Args:
+            monomer_roles_in (list[MonomerRoleforIndexBasedFGDetection]): List of monomer
+                roles to process, each carrying the atom indices of interest.
+
+        Returns:
+            list[MonomerRole] | bool: List of MonomerRole objects with index-filtered
+            functionalities, or False if none detected.
+
+        Notes:
+            - Index-based rule: at least ONE match overlapping the given indices is
+            enough to qualify, regardless of functionality_type. This intentionally
+            breaks the whole-molecule 'di_identical' (>=2 matches) rule, since here
+            we only care whether the given index sits inside a valid functional group,
+            not how many total sites exist on the monomer.
+        """
+
+        monomer_roles_out = []
+
+        for monomer in monomer_roles_in:
+            if monomer.is_looped:
+                continue  # Skip already processed monomers
+
+            mol = monomer.rdkit_mol
+
+            target_indices = set(monomer.indexes_in_template or [])
+            detected_functionalities = []
+            all_matches = []
+
+            # Check against each predefined functional group type.
+            for functional_group in self.monomer_types.values():
+                ftype = functional_group["functionality_type"]
+                smarts_1 = functional_group["smarts_1"]
+                smarts_2 = functional_group.get("smarts_2")
+
+                patt1 = Chem.MolFromSmarts(smarts_1)
+                if patt1 is None:
+                    logger.warning(f"Invalid primary SMARTS: {smarts_1}")
+                    continue
+
+                matches1 = mol.GetSubstructMatches(patt1, uniquify=True)
+                # Index-based filter: keep only matches touching at least one target index.
+                matches1_hit = [m for m in matches1 if target_indices.intersection(m)]
+                count_1 = len(matches1_hit)
+
+                count_2 = None
+                matches2_hit = []
+
+                if smarts_2:
+                    patt2 = Chem.MolFromSmarts(smarts_2)
+                    if patt2 is None:
+                        logger.warning(f"Invalid secondary SMARTS: {smarts_2}")
+                        continue
+
+                    matches2 = mol.GetSubstructMatches(patt2, uniquify=True)
+                    matches2_hit = [m for m in matches2 if target_indices.intersection(m)]
+                    count_2 = len(matches2_hit)
+
+                    # di_different: still need one overlapping hit on EACH pattern.
+                    functionality_count = 2 if (count_1 >= 1 and count_2 >= 1) else 0
+                else:
+                    # vinyl / mono / di_identical: ONE overlapping match is enough.
+                    # (Breaks the normal di_identical >=2 rule on purpose for index-based detection.)
+                    functionality_count = 1 if count_1 >= 1 else 0
+
+                if functionality_count > 0:
+                    functional_matches = tuple(matches1_hit) + tuple(matches2_hit)
+                    all_matches.extend(functional_matches)
+
+                    # Log detected functionality for debugging/user feedback.
+                    # print(f"{monomer.smiles} has functionality: {functional_group['group_name']}") 
+
+                    detected_functionalities.append(
+                        FunctionalGroupInfo(
+                            functionality_type=ftype,
+                            fg_name=functional_group["group_name"],
+                            fg_smarts_1=smarts_1,
+                            fg_count_1=count_1,
+                            fg_1_indexes=tuple(matches1_hit) if matches1_hit else None,
+                            fg_smarts_2=smarts_2,
+                            fg_count_2=count_2,
+                            fg_2_indexes=tuple(matches2_hit) if matches2_hit else None,
+                        )
+                    )
+
+            # Add to roles if any functionalities detected.
+            if detected_functionalities:
+                monomer_roles_out.append(
+                    MonomerRole(
+                        smiles=monomer.smiles,
+                        name=monomer.name,
+                        rdkit_mol=monomer.rdkit_mol,
+                        functionalities=tuple(detected_functionalities),
+                        is_monomer=False,  # This is a product, not an input monomer
+                        is_looped=False,   # Yet to be processed in the loop
+                        indexes_in_template=monomer.indexes_in_template,
+                    )
+                )
+
+        # Store results for potential downstream use.
+        if not monomer_roles_out:
+            return False  # No functional groups detected; handle as needed
+        # first break condition: if no monomer roles are detected, return False to indicate no further processing is needed.
+
+        return monomer_roles_out  # Return list of MonomerRole; visualization not considered here.
 
     def functional_group_highlighted_molecules_image_grid(self, session: Session) -> Image:
         """Convert monomer roles with detected functionalities into visualizations.