Reacter¶

The reacter module provides tools for chemical reactions and topology modification.

Base¶

Core Reacter implementation for chemical transformations.

This module defines the base Reacter class and ProductSet dataclass, providing the foundation for SMIRKS-style reaction semantics.

AnchorSelector `module-attribute` ¶

AnchorSelector = Callable[[Atomistic, Atom], Atom]

Select anchor atom given a port atom.

Port atoms are SMILES-marked connection sites (e.g. $, , <, >) stored on atoms via the "port" / "ports" attributes. An anchor* is the actual atom that participates in bond formation.

Parameters:

Name	Type	Description	Default
`assembly`		The Atomistic structure to select from	required
`port_atom`		The port atom entity (SMILES-marked atom)	required

Returns:

Name	Type	Description
`Atom`		The anchor atom where the new bond should be formed

BondFormer `module-attribute` ¶

BondFormer = Callable[[Atomistic, Atom, Atom], Bond | None]

Create or modify bonds between two anchor atoms in an assembly.

Parameters:

Name	Description	Default
`assembly`	The atomistic assembly to modify	required
`i`	First anchor atom	required
`j`	Second anchor atom	required

Side effects

Adds or updates bonds in assembly.links

LeavingSelector `module-attribute` ¶

LeavingSelector = Callable[[Atomistic, Atom], list[Atom]]

Select leaving group atoms given an anchor atom.

Parameters:

Name	Type	Description	Default
`assembly`		The Atomistic structure containing the atoms	required
`anchor_atom`		The anchor atom entity (actual reaction site)	required

Returns:

Type	Description
	List of atom entities to be removed

ProductInfo `dataclass` ¶

ProductInfo(product, anchor_L=None, anchor_R=None)

Information about the reaction product.

This captures the final product structure and the anchor atoms where the new bond was formed.

ReactantInfo `dataclass` ¶

ReactantInfo(merged_reactants, port_atom_L=None, port_atom_R=None)

Information about the reactants in a reaction.

Ports vs anchors: - Ports are SMILES-marked connection atoms (e.g. $, , <, >) - Anchors* are the actual atoms where bonds are formed

This dataclass tracks the merged reactants and the original port atoms on each side before the reaction is executed.

Reacter ¶

Reacter(name, anchor_selector_left, anchor_selector_right, leaving_selector_left, leaving_selector_right, bond_former)

Programmable chemical reaction executor.

A Reacter represents one specific chemical reaction type by composing: 1. Anchor selectors - map port atoms to anchor atoms 2. Leaving selectors - identify atoms to remove 3. Bond former - create new bonds between anchor atoms

The reaction is executed on copies of input monomers, ensuring original structures remain unchanged.

Port Selection Philosophy: Reacter does NOT handle port selection. The caller (e.g., MonomerLinker) must explicitly specify which ports to connect via port_L and port_R. Ports are marked directly on atoms using the "port" or "ports" attribute. This makes the reaction execution deterministic and explicit.

Attributes:

Name	Type	Description
`name`		Descriptive name for this reaction type
`anchor_selector_left`		Function to map left port atom to anchor atom
`anchor_selector_right`		Function to map right port atom to anchor atom
`leaving_selector_left`		Function to select left leaving group from left anchor
`leaving_selector_right`		Function to select right leaving group from right anchor
`bond_former`		Function to create bond between anchor atoms

Example

from molpy.reacter import Reacter, select_port_atom, select_one_hydrogen, form_single_bond from molpy import Atomistic

Mark ports on atoms¶

atom_a["port"] = "1" atom_b["port"] = "2"

cc_coupling = Reacter( ... name="C-C_coupling_with_H_loss", ... port_selector_left=select_port_atom, ... port_selector_right=select_port_atom, ... leaving_selector_left=select_one_hydrogen, ... leaving_selector_right=select_one_hydrogen, ... bond_former=form_single_bond, ... )

Port selection is explicit!¶

product = cc_coupling.run(structA, structB, port_L="1", port_R="2") print(product.removed_atoms) # [H1, H2]

Initialize a Reacter with reaction components.

Parameters:

Name	Type	Description	Default
`name`	`str`	Descriptive name for this reaction	required
`anchor_selector_left`	`AnchorSelector`	Selector that maps left port atom to anchor atom	required
`anchor_selector_right`	`AnchorSelector`	Selector that maps right port atom to anchor atom	required
`leaving_selector_left`	`LeavingSelector`	Selector for left leaving group (from left anchor)	required
`leaving_selector_right`	`LeavingSelector`	Selector for right leaving group (from right anchor)	required
`bond_former`	`BondFormer`	Function to create bond between anchor atoms	required

Source code in src/molpy/reacter/base.py

def __init__(
    self,
    name: str,
    anchor_selector_left: AnchorSelector,
    anchor_selector_right: AnchorSelector,
    leaving_selector_left: LeavingSelector,
    leaving_selector_right: LeavingSelector,
    bond_former: BondFormer,
):
    """
    Initialize a Reacter with reaction components.

    Args:
        name: Descriptive name for this reaction
        anchor_selector_left: Selector that maps left **port atom** to **anchor atom**
        anchor_selector_right: Selector that maps right **port atom** to **anchor atom**
        leaving_selector_left: Selector for left leaving group (from left anchor)
        leaving_selector_right: Selector for right leaving group (from right anchor)
        bond_former: Function to create bond between anchor atoms
    """
    self.name = name
    self.anchor_selector_left = anchor_selector_left
    self.anchor_selector_right = anchor_selector_right
    self.leaving_selector_left = leaving_selector_left
    self.leaving_selector_right = leaving_selector_right
    self.bond_former = bond_former

run ¶

run(left, right, port_atom_L, port_atom_R, compute_topology=True, record_intermediates=False, typifier=None)

Execute the reaction between two Atomistic structures.

IMPORTANT: port_atom_L and port_atom_R must be explicit Atom objects. Use find_port_atom() or find_port_atom_by_node() to get them first.

Workflow: 1. Transform port atoms to reaction sites via port selectors 2. Merge structures (or work on single copy for ring closure) 3. Select leaving groups from reaction sites 4. Create bond between reaction sites 5. Remove leaving groups 6. (Optional) Compute new angles/dihedrals 7. Return ReactionResult

Parameters:

Name	Type	Description	Default
`left`	`Atomistic`	Left reactant Atomistic structure	required
`right`	`Atomistic`	Right reactant Atomistic structure	required
`port_atom_L`	`Entity`	Port atom from left structure (the atom with port marker)	required
`port_atom_R`	`Entity`	Port atom from right structure (the atom with port marker)	required
`compute_topology`	`bool`	If True, compute new angles/dihedrals (default True)	`True`
`record_intermediates`	`bool`	If True, record intermediate states	`False`
`typifier`	`TypifierBase \| None`	Optional typifier for incremental retypification	`None`

Returns:

Type	Description
`ReactionResult`	ReactionResult containing product and metadata

Raises:

Type	Description
`ValueError`	If port atoms invalid

Source code in src/molpy/reacter/base.py

def run(
    self,
    left: Atomistic,
    right: Atomistic,
    port_atom_L: Entity,
    port_atom_R: Entity,
    compute_topology: bool = True,
    record_intermediates: bool = False,
    typifier: TypifierBase | None = None,
) -> ReactionResult:
    """
    Execute the reaction between two Atomistic structures.

    **IMPORTANT: port_atom_L and port_atom_R must be explicit Atom objects.**
    Use find_port_atom() or find_port_atom_by_node() to get them first.

    Workflow:
    1. Transform port atoms to reaction sites via port selectors
    2. Merge structures (or work on single copy for ring closure)
    3. Select leaving groups from reaction sites
    4. Create bond between reaction sites
    5. Remove leaving groups
    6. (Optional) Compute new angles/dihedrals
    7. Return ReactionResult

    Args:
        left: Left reactant Atomistic structure
        right: Right reactant Atomistic structure
        port_atom_L: Port atom from left structure (the atom with port marker)
        port_atom_R: Port atom from right structure (the atom with port marker)
        compute_topology: If True, compute new angles/dihedrals (default True)
        record_intermediates: If True, record intermediate states
        typifier: Optional typifier for incremental retypification

    Returns:
        ReactionResult containing product and metadata

    Raises:
        ValueError: If port atoms invalid
    """
    intermediates: list[dict] = []

    # Check if this is a ring closure (left and right are the same object)
    is_ring_closure = left is right

    if is_ring_closure:
        # For ring closure, work directly on a single copy of the structure
        merged = left.copy()
        for i, atom in enumerate(merged.atoms, start=1):
            atom["id"] = i

        # Build entity map: original -> copied
        original_atoms = list(left.atoms)
        copied_atoms = list(merged.atoms)
        entity_map: dict[Entity, Entity] = {}
        if len(original_atoms) == len(copied_atoms):
            for orig, copied in zip(original_atoms, copied_atoms):
                entity_map[orig] = copied

        left_entity_map = entity_map
        right_entity_map = entity_map

        # Map port atoms to their copies
        copied_port_L = entity_map.get(port_atom_L, port_atom_L)
        copied_port_R = entity_map.get(port_atom_R, port_atom_R)

        # Transform port atoms to anchors
        anchor_L = self.anchor_selector_left(merged, copied_port_L)
        anchor_R = self.anchor_selector_right(merged, copied_port_R)

        # Select leaving groups from anchors
        leaving_L = self.leaving_selector_left(merged, anchor_L)
        leaving_R = self.leaving_selector_right(merged, anchor_R)
    else:
        # Normal case: prepare reactants (creates copies)
        left_copy, right_copy, left_entity_map, right_entity_map = (
            self._prepare_reactants(left, right)
        )

        # Map port atoms to their copies
        copied_port_L = left_entity_map.get(port_atom_L, port_atom_L)
        copied_port_R = right_entity_map.get(port_atom_R, port_atom_R)

        # Transform port atoms to anchors
        anchor_L = self.anchor_selector_left(left_copy, copied_port_L)
        anchor_R = self.anchor_selector_right(right_copy, copied_port_R)

        # Select leaving groups from anchors
        leaving_L = self.leaving_selector_left(left_copy, anchor_L)
        leaving_R = self.leaving_selector_right(right_copy, anchor_R)

        # Merge structures
        merged = self._merge_structures(left_copy, right_copy)

    # Save merged reactants BEFORE reaction (for template generation)
    merged_reactants_before_reaction = merged.copy()
    merged_copy = merged.copy()

    if record_intermediates:
        if compute_topology:
            merged.get_topo(gen_angle=True, gen_dihe=True)
        intermediates.append(
            {
                "step": "reactants",
                "description": "Reactants before bond formation",
                "product": merged_copy,
            }
        )

    # Step 4: Execute reaction (bond formation and leaving group removal)
    new_bond, removed_atoms = self._execute_reaction(
        merged,
        anchor_L,
        anchor_R,
        leaving_L,
        leaving_R,
    )

    if record_intermediates:
        product_copy = merged_copy.copy()
        if compute_topology:
            product_copy.get_topo(gen_angle=True, gen_dihe=True)

        intermediates.append(
            {
                "step": "bond_formation",
                "description": "After forming new bond between port atoms",
                "product": product_copy,
                "new_bond": new_bond,
            }
        )

        intermediates.append(
            {
                "step": "remove_leaving",
                "description": f"After removing {len(removed_atoms)} leaving atoms",
                "product": product_copy,
                "removed_atoms": removed_atoms,
            }
        )

    # Step 5: Detect and update topology if requested
    new_angles: list[Angle] = []
    new_dihedrals: list[Dihedral] = []
    removed_angles: list[Angle] = []
    removed_dihedrals: list[Dihedral] = []

    if compute_topology and new_bond:
        new_angles, new_dihedrals, removed_angles, removed_dihedrals = (
            self._detect_and_update_topology(merged, [new_bond], removed_atoms)
        )

    # Step 6: Build entity maps
    final_entity_map = self._build_entity_maps(
        left_entity_map, right_entity_map, merged
    )

    # Step 7: Determine if retypification is needed
    requires_retype = bool(new_bond or removed_atoms)

    # Step 8: Build result structure
    # Use merged reactants saved BEFORE reaction (with all atoms including removed_atoms)
    reactant_info = ReactantInfo(
        merged_reactants=merged_reactants_before_reaction,
        port_atom_L=port_atom_L,
        port_atom_R=port_atom_R,
    )

    product_info = ProductInfo(product=merged, anchor_L=anchor_L, anchor_R=anchor_R)

    topology_changes = TopologyChanges(
        new_bonds=[new_bond] if new_bond else [],
        new_angles=new_angles,
        new_dihedrals=new_dihedrals,
        removed_angles=removed_angles,
        removed_dihedrals=removed_dihedrals,
        removed_atoms=removed_atoms,
        modified_atoms=({anchor_L, anchor_R} if anchor_L and anchor_R else set()),
    )

    metadata = ReactionMetadata(
        reaction_name=self.name,
        requires_retype=requires_retype,
        entity_maps=[final_entity_map],
        intermediates=intermediates,
    )

    result = ReactionResult(
        reactant_info=reactant_info,
        product_info=product_info,
        topology_changes=topology_changes,
        metadata=metadata,
    )

    # Step 9: Perform incremental typification if requested
    if typifier:
        self._incremental_typify(merged, result, typifier)

    return result

ReactionMetadata `dataclass` ¶

ReactionMetadata(reaction_name, requires_retype=False, entity_maps=list(), intermediates=list())

Metadata about the reaction.

ReactionResult `dataclass` ¶

ReactionResult(reactant_info, product_info, topology_changes, metadata)

Container for reaction products and metadata with organized structure.

This class organizes reaction information into logical groups:

reactant_info: Information about the reactants and their port atoms
product_info: Information about the product and the anchor atoms
topology_changes: All topology changes (bonds, angles, dihedrals)
metadata: Reaction metadata (name, retyping info, etc.)

TopologyChanges `dataclass` ¶

TopologyChanges(new_bonds=list(), new_angles=list(), new_dihedrals=list(), removed_angles=list(), removed_dihedrals=list(), removed_atoms=list(), modified_atoms=set())

Topology changes resulting from the reaction.

Connector¶

Connector interface for integrating Reacter with the linear polymer builder.

This module provides connectors that use Reacter for polymer assembly, allowing flexible reaction specification with default and specialized reactors.

MonomerLinker ¶

MonomerLinker(default_reaction, specialized_reactions=None)

Connector adapter that manages Reacter instances for polymer assembly.

This connector allows specifying a default reaction for most connections, with the ability to override specific monomer pair connections with specialized reacters.

Port Selection Philosophy: This connector does NOT handle port selection. The caller must explicitly provide port_L and port_R when calling connect(). Ports must be marked directly on atoms using the "port" or "ports" attribute. Port selection should be handled by the higher-level builder logic.

Example

from molpy.reacter import Reacter, MonomerLinker from molpy import Atomistic

Mark ports on atoms¶

atom_a["port"] = "1" atom_b["port"] = "2"

Default reaction for most connections¶

default_reacter = Reacter(...)

Special reaction for A-B connection¶

special_reacter = Reacter(...)

connector = MonomerLinker( ... default_reaction=default_reacter, ... specialized_reactions={('A', 'B'): special_reacter} ... )

Explicit port specification required¶

product = connector.connect( ... left=struct_a, right=struct_b, ... left_type='A', right_type='B', ... port_L='1', port_R='2' # REQUIRED ... )

Initialize connector with default and specialized reacters.

Parameters:

Name	Type	Description	Default
`default_reaction`	`Reacter`	Default Reacter used for most connections	required
`specialized_reactions`	`dict[tuple[str, str], Reacter] \| None`	Dict mapping (left_type, right_type) -> specialized Reacter	`None`

Source code in src/molpy/reacter/connector.py

def __init__(
    self,
    default_reaction: Reacter,
    specialized_reactions: dict[tuple[str, str], Reacter] | None = None,
):
    """
    Initialize connector with default and specialized reacters.

    Args:
        default_reaction: Default Reacter used for most connections
        specialized_reactions: Dict mapping (left_type, right_type) -> specialized Reacter
    """
    self.default_reaction = default_reaction
    self.specialized_reactions = specialized_reactions or {}
    self._history: list[ReactionResult] = []

clear_history ¶

clear_history()

Clear reaction history.

Source code in src/molpy/reacter/connector.py

def clear_history(self):
    """Clear reaction history."""
    self._history.clear()

connect ¶

connect(left, right, port_L, port_R, left_type=None, right_type=None)

Connect two Atomistic structures using appropriate reacter.

IMPORTANT: port_L and port_R must be explicitly specified. Ports must be marked on atoms using the "port" or "ports" attribute. No automatic port selection or fallback is performed.

This connector is responsible for locating port atoms by name. The underlying :class:Reacter only operates on anchors and receives explicit port atoms which it converts to anchors via its anchor_selector_left / anchor_selector_right callables.

Parameters:

Name	Type	Description	Default
`left`	`Atomistic`	Left Atomistic structure	required
`right`	`Atomistic`	Right Atomistic structure	required
`port_L`	`str`	Port name on left structure (REQUIRED)	required
`port_R`	`str`	Port name on right structure (REQUIRED)	required
`left_type`	`str \| None`	Type label for left structure (e.g., 'A', 'B')	`None`
`right_type`	`str \| None`	Type label for right structure	`None`

Returns:

Type	Description
`Atomistic`	Connected Atomistic assembly

Raises:

Type	Description
`ValueError`	If ports not found on structures

Source code in src/molpy/reacter/connector.py

def connect(
    self,
    left: Atomistic,
    right: Atomistic,
    port_L: str,
    port_R: str,
    left_type: str | None = None,
    right_type: str | None = None,
) -> Atomistic:
    """
    Connect two Atomistic structures using appropriate reacter.

    **IMPORTANT: port_L and port_R must be explicitly specified.**
    Ports must be marked on atoms using the "port" or "ports" attribute.
    No automatic port selection or fallback is performed.

    This connector is responsible for locating **port atoms** by name.
    The underlying :class:`Reacter` only operates on **anchors** and
    receives explicit port atoms which it converts to anchors via its
    ``anchor_selector_left`` / ``anchor_selector_right`` callables.

    Args:
        left: Left Atomistic structure
        right: Right Atomistic structure
        port_L: Port name on left structure (REQUIRED)
        port_R: Port name on right structure (REQUIRED)
        left_type: Type label for left structure (e.g., 'A', 'B')
        right_type: Type label for right structure

    Returns:
        Connected Atomistic assembly

    Raises:
        ValueError: If ports not found on structures
    """
    # Validate ports exist by checking if any atom has the port marker
    left_has_port = any(
        atom.get("port") == port_L
        or (isinstance(atom.get("ports"), list) and port_L in atom.get("ports", []))
        for atom in left.atoms
    )
    if not left_has_port:
        raise ValueError(
            f"Port '{port_L}' not found on left structure. "
            f"Atoms must have 'port' or 'ports' attribute set to '{port_L}'"
        )

    right_has_port = any(
        atom.get("port") == port_R
        or (isinstance(atom.get("ports"), list) and port_R in atom.get("ports", []))
        for atom in right.atoms
    )
    if not right_has_port:
        raise ValueError(
            f"Port '{port_R}' not found on right structure. "
            f"Atoms must have 'port' or 'ports' attribute set to '{port_R}'"
        )

    # Select reacter based on monomer types
    reacter = self._select_reacter(left_type, right_type)

    # Locate concrete port atoms – Reacter itself only works with anchors
    from molpy.reacter.selectors import find_port_atom

    port_atom_L = find_port_atom(left, port_L)
    port_atom_R = find_port_atom(right, port_R)

    # Execute reaction on explicit port atoms
    result: ReactionResult = reacter.run(
        left,
        right,
        port_atom_L=port_atom_L,
        port_atom_R=port_atom_R,
    )

    # Store in history for retypification
    self._history.append(result)

    return result.product_info.product

get_all_modified_atoms ¶

get_all_modified_atoms()

Get all atoms that have been modified across all reactions.

Returns:

Type	Description
`set[Entity]`	Set of all atoms that need retypification

Source code in src/molpy/reacter/connector.py

def get_all_modified_atoms(self) -> set[Entity]:
    """
    Get all atoms that have been modified across all reactions.

    Returns:
        Set of all atoms that need retypification
    """
    modified: set[Entity] = set()
    for result in self._history:
        modified.update(result.topology_changes.modified_atoms)
    return modified

get_history ¶

get_history()

Get history of all reactions performed.

Useful for batch retypification after polymer assembly.

Returns:

Type	Description
`list[ReactionResult]`	List of ReactionResult for each connection made

Source code in src/molpy/reacter/connector.py

def get_history(self) -> list[ReactionResult]:
    """
    Get history of all reactions performed.

    Useful for batch retypification after polymer assembly.

    Returns:
        List of ReactionResult for each connection made
    """
    return self._history.copy()

needs_retypification ¶

needs_retypification()

Check if any reactions require retypification.

Returns:

Type	Description
`bool`	True if retypification needed

Source code in src/molpy/reacter/connector.py

def needs_retypification(self) -> bool:
    """
    Check if any reactions require retypification.

    Returns:
        True if retypification needed
    """
    return any(r.metadata.requires_retype for r in self._history)

Selectors¶

Selector functions for identifying reaction sites and leaving groups.

Selectors are composable functions that identify specific atoms in a structure for use in reactions.

Selector Types:

Port Selectors (port_selector_left, port_selector_right):
Signature: (struct: Atomistic, port_atom: Entity) -> Entity
Transform a port atom to its actual reaction site
Example: O port → C neighbor for dehydration
Leaving Selectors (leaving_selector_left, leaving_selector_right):
Signature: (struct: Atomistic, reaction_site: Entity) -> list[Entity]
Identify atoms to remove from the reaction site
Example: C → [O, H] for hydroxyl removal

Naming: - struct: The Atomistic structure containing atoms - port_atom: The specific atom with port marker - reaction_site: The atom where the bond will form

find_port_atom ¶

find_port_atom(struct, port_name)

Find an atom with the specified port marker.

This is a utility function, not a selector. Use it to find port atoms before passing them to selectors.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure containing ports	required
`port_name`	`str`	Name of port to find	required

Returns:

Type	Description
`Atom`	Atom with matching port

Raises:

Type	Description
`ValueError`	If port not found

Source code in src/molpy/reacter/selectors.py

def find_port_atom(struct: Atomistic, port_name: str) -> Atom:
    """
    Find an atom with the specified port marker.

    This is a utility function, not a selector. Use it to find
    port atoms before passing them to selectors.

    Args:
        struct: Atomistic structure containing ports
        port_name: Name of port to find

    Returns:
        Atom with matching port

    Raises:
        ValueError: If port not found
    """
    for atom in struct.atoms:
        if atom.get("port") == port_name:
            return atom  # type: ignore[return-value]
        ports = atom.get("ports")
        if isinstance(ports, list) and port_name in ports:
            return atom

    raise ValueError(f"Port '{port_name}' not found")

find_port_atom_by_node ¶

find_port_atom_by_node(struct, port_name, node_id)

Find port atom for a specific node ID.

Useful when structure has multiple nodes and you need to find the port for a specific one.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`port_name`	`str`	Name of port	required
`node_id`	`int`	The monomer_node_id to match	required

Returns:

Type	Description
`Atom`	Atom with matching port and node_id

Raises:

Type	Description
`ValueError`	If port not found for the specified node

Source code in src/molpy/reacter/selectors.py

def find_port_atom_by_node(struct: Atomistic, port_name: str, node_id: int) -> Atom:
    """
    Find port atom for a specific node ID.

    Useful when structure has multiple nodes and you need
    to find the port for a specific one.

    Args:
        struct: Atomistic structure
        port_name: Name of port
        node_id: The monomer_node_id to match

    Returns:
        Atom with matching port and node_id

    Raises:
        ValueError: If port not found for the specified node
    """
    for atom in struct.atoms:
        if atom.get("monomer_node_id") != node_id:
            continue
        if atom.get("port") == port_name:
            return atom  # type: ignore[return-value]
        ports = atom.get("ports")
        if isinstance(ports, list) and port_name in ports:
            return atom

    raise ValueError(f"Port '{port_name}' not found for node {node_id}")

select_all_hydrogens ¶

select_all_hydrogens(struct, reaction_site)

Select all hydrogen atoms bonded to the reaction site.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`reaction_site`	`Atom`	Atom to find H neighbors of	required

Returns:

Type	Description
`list[Atom]`	List of all H atoms bonded to reaction site

Source code in src/molpy/reacter/selectors.py

def select_all_hydrogens(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select all hydrogen atoms bonded to the reaction site.

    Args:
        struct: Atomistic structure
        reaction_site: Atom to find H neighbors of

    Returns:
        List of all H atoms bonded to reaction site
    """
    return [
        a
        for a in find_neighbors(struct, reaction_site, element="H")
        if isinstance(a, Atom)
    ]

select_c_neighbor ¶

select_c_neighbor(struct, port_atom)

Select C neighbor of the port atom as reaction site.

Useful when port is on O but reaction site is adjacent C.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`port_atom`	`Atom`	Port atom (typically O)	required

Returns:

Type	Description
`Atom`	First C neighbor of port_atom

Source code in src/molpy/reacter/selectors.py

def select_c_neighbor(struct: Atomistic, port_atom: Atom) -> Atom:
    """
    Select C neighbor of the port atom as reaction site.

    Useful when port is on O but reaction site is adjacent C.

    Args:
        struct: Atomistic structure
        port_atom: Port atom (typically O)

    Returns:
        First C neighbor of port_atom
    """
    c_neighbors = [
        a for a in find_neighbors(struct, port_atom, element="C") if isinstance(a, Atom)
    ]
    return c_neighbors[0]

select_dehydration_left ¶

select_dehydration_left(struct, port_atom)

Left-side selector for dehydration reactions.

Returns C as reaction site, handling ports on either O or C: - If port on O: returns C neighbor - If port on C: returns C itself

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`port_atom`	`Atom`	Port atom (O or C)	required

Returns:

Type	Description
`Atom`	C atom as reaction site

Source code in src/molpy/reacter/selectors.py

def select_dehydration_left(struct: Atomistic, port_atom: Atom) -> Atom:
    """
    Left-side selector for dehydration reactions.

    Returns C as reaction site, handling ports on either O or C:
    - If port on O: returns C neighbor
    - If port on C: returns C itself

    Args:
        struct: Atomistic structure
        port_atom: Port atom (O or C)

    Returns:
        C atom as reaction site
    """
    symbol = port_atom.get("symbol")

    if symbol == "O":
        c_neighbors = [
            a
            for a in find_neighbors(struct, port_atom, element="C")
            if isinstance(a, Atom)
        ]
        return c_neighbors[0]
    else:  # C
        return port_atom

select_dehydration_right ¶

select_dehydration_right(struct, port_atom)

Right-side selector for dehydration reactions.

Returns O as reaction site, handling ports on either O or C: - If port on O: returns O itself - If port on C: returns O neighbor

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`port_atom`	`Atom`	Port atom (O or C)	required

Returns:

Type	Description
`Atom`	O atom as reaction site

Source code in src/molpy/reacter/selectors.py

def select_dehydration_right(struct: Atomistic, port_atom: Atom) -> Atom:
    """
    Right-side selector for dehydration reactions.

    Returns O as reaction site, handling ports on either O or C:
    - If port on O: returns O itself
    - If port on C: returns O neighbor

    Args:
        struct: Atomistic structure
        port_atom: Port atom (O or C)

    Returns:
        O atom as reaction site
    """
    symbol = port_atom.get("symbol")

    if symbol == "O":
        return port_atom
    else:  # C
        o_neighbors = [
            a
            for a in find_neighbors(struct, port_atom, element="O")
            if isinstance(a, Atom)
        ]
        return o_neighbors[0]

select_dummy_atoms ¶

select_dummy_atoms(struct, reaction_site)

Select dummy atoms (*) bonded to the reaction site.

Useful for BigSMILES-style reactions where * marks connection points.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`reaction_site`	`Atom`	Atom to find dummy neighbors of	required

Returns:

Type	Description
`list[Atom]`	List of dummy atoms (symbol='*')

Source code in src/molpy/reacter/selectors.py

def select_dummy_atoms(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select dummy atoms (*) bonded to the reaction site.

    Useful for BigSMILES-style reactions where * marks connection points.

    Args:
        struct: Atomistic structure
        reaction_site: Atom to find dummy neighbors of

    Returns:
        List of dummy atoms (symbol='*')
    """
    neighbors = [
        a for a in find_neighbors(struct, reaction_site) if isinstance(a, Atom)
    ]
    return [n for n in neighbors if n.get("symbol") == "*"]

select_hydroxyl_group ¶

select_hydroxyl_group(struct, reaction_site)

Select hydroxyl group (-OH) bonded to the reaction site.

Finds single-bonded O neighbor, then finds H bonded to that O. Used for esterification and dehydration reactions.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`reaction_site`	`Atom`	Atom (typically C) with -OH attached	required

Returns:

Type	Description
`list[Atom]`	[O, H] atoms from hydroxyl group

Source code in src/molpy/reacter/selectors.py

def select_hydroxyl_group(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select hydroxyl group (-OH) bonded to the reaction site.

    Finds single-bonded O neighbor, then finds H bonded to that O.
    Used for esterification and dehydration reactions.

    Args:
        struct: Atomistic structure
        reaction_site: Atom (typically C) with -OH attached

    Returns:
        [O, H] atoms from hydroxyl group
    """

    o_neighbors = [
        a
        for a in find_neighbors(struct, reaction_site, element="O")
        if isinstance(a, Atom)
    ]

    if not o_neighbors:
        raise ValueError("No oxygen neighbors found for reaction site")

    # Find single-bonded oxygen (hydroxyl, not carbonyl)
    hydroxyl_o = None
    for o in o_neighbors:
        bond = get_bond_between(struct, reaction_site, o)
        if bond and bond.get("order") == 1:
            hydroxyl_o = o
            break

    if hydroxyl_o is None:
        raise ValueError("No single-bonded oxygen (hydroxyl) found for reaction site")

    h_neighbors = [
        a
        for a in find_neighbors(struct, hydroxyl_o, element="H")
        if isinstance(a, Atom)
    ]

    if not h_neighbors:
        raise ValueError("No hydrogen found bonded to hydroxyl oxygen")

    return [hydroxyl_o, h_neighbors[0]]

select_hydroxyl_h_only ¶

select_hydroxyl_h_only(struct, reaction_site)

Select only the H from hydroxyl group bonded to reaction site.

The O remains (becomes part of the new bond).

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`reaction_site`	`Atom`	Atom with -OH attached (typically O itself)	required

Returns:

Type	Description
`list[Atom]`	[H] if found, otherwise []

Source code in src/molpy/reacter/selectors.py

def select_hydroxyl_h_only(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select only the H from hydroxyl group bonded to reaction site.

    The O remains (becomes part of the new bond).

    Args:
        struct: Atomistic structure
        reaction_site: Atom with -OH attached (typically O itself)

    Returns:
        [H] if found, otherwise []
    """
    # If reaction_site is O, look for H directly bonded
    if reaction_site.get("symbol") == "O":
        h_neighbors = [
            a
            for a in find_neighbors(struct, reaction_site, element="H")
            if isinstance(a, Atom)
        ]
        if h_neighbors:
            return [h_neighbors[0]]
        return []

    # If reaction_site is C, look for O neighbor then H
    from .utils import get_bond_between

    o_neighbors = [
        a
        for a in find_neighbors(struct, reaction_site, element="O")
        if isinstance(a, Atom)
    ]
    for o in o_neighbors:
        bond = get_bond_between(struct, reaction_site, o)
        if not bond or bond.get("order") != 1:
            continue
        h_neighbors = [
            a for a in find_neighbors(struct, o, element="H") if isinstance(a, Atom)
        ]
        if h_neighbors:
            return [h_neighbors[0]]
    return []

select_none ¶

select_none(struct, reaction_site)

Select no leaving group - returns empty list.

Useful for addition reactions where nothing is eliminated.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure (ignored)	required
`reaction_site`	`Atom`	Reaction site atom (ignored)	required

Returns:

Type	Description
`list[Atom]`	Empty list

Source code in src/molpy/reacter/selectors.py

def select_none(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select no leaving group - returns empty list.

    Useful for addition reactions where nothing is eliminated.

    Args:
        struct: Atomistic structure (ignored)
        reaction_site: Reaction site atom (ignored)

    Returns:
        Empty list
    """
    return []

select_o_neighbor ¶

select_o_neighbor(struct, port_atom)

Select O neighbor of the port atom as reaction site.

Useful when port is on C but reaction site is adjacent O.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`port_atom`	`Atom`	Port atom (typically C)	required

Returns:

Type	Description
`Atom`	First O neighbor of port_atom

Source code in src/molpy/reacter/selectors.py

def select_o_neighbor(struct: Atomistic, port_atom: Atom) -> Atom:
    """
    Select O neighbor of the port atom as reaction site.

    Useful when port is on C but reaction site is adjacent O.

    Args:
        struct: Atomistic structure
        port_atom: Port atom (typically C)

    Returns:
        First O neighbor of port_atom
    """
    o_neighbors = [
        a for a in find_neighbors(struct, port_atom, element="O") if isinstance(a, Atom)
    ]
    return o_neighbors[0]

select_one_hydrogen ¶

select_one_hydrogen(struct, reaction_site)

Select one hydrogen atom bonded to the reaction site.

Useful for condensation reactions where H is eliminated.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure	required
`reaction_site`	`Atom`	Atom to find H neighbor of	required

Returns:

Type	Description
`list[Atom]`	List with one H atom, or empty list if none

Source code in src/molpy/reacter/selectors.py

def select_one_hydrogen(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """
    Select one hydrogen atom bonded to the reaction site.

    Useful for condensation reactions where H is eliminated.

    Args:
        struct: Atomistic structure
        reaction_site: Atom to find H neighbor of

    Returns:
        List with one H atom, or empty list if none
    """
    h_neighbors = [
        a
        for a in find_neighbors(struct, reaction_site, element="H")
        if isinstance(a, Atom)
    ]
    if h_neighbors:
        return [h_neighbors[0]]
    return []

select_port ¶

select_port(struct, port_atom)

Port selector - returns the port atom as the reaction site.

Use this when the port atom itself is the reaction site.

Parameters:

Name	Type	Description	Default
`struct`	`Atomistic`	Atomistic structure containing the atoms	required
`port_atom`	`Atom`	The atom with port marker	required

Returns:

Type	Description
`Atom`	The same port_atom

Source code in src/molpy/reacter/selectors.py

def select_port(struct: Atomistic, port_atom: Atom) -> Atom:
    """
    Port selector - returns the port atom as the reaction site.

    Use this when the port atom itself is the reaction site.

    Args:
        struct: Atomistic structure containing the atoms
        port_atom: The atom with port marker

    Returns:
        The same port_atom
    """
    return port_atom

Template¶

TemplateReacter: Specialized Reacter for LAMMPS fix bond/react template generation.

This module provides a wrapper around the base Reacter that handles: 1. react_id assignment for atom tracking across reactions 2. Pre/post template extraction with correct atom/bond references 3. Writing LAMMPS molecule and map files

TemplateReacter ¶

TemplateReacter(name, anchor_selector_left, anchor_selector_right, leaving_selector_left, leaving_selector_right, bond_former, radius=4)

Reacter wrapper for LAMMPS fix bond/react template generation.

Composes a Reacter internally and adds: 1. react_id assignment before reaction 2. Subgraph extraction for templates 3. Correct atom/bond reference handling

The function signature matches Reacter, but run_with_template() returns both ReactionResult and TemplateResult.

Example

from molpy.reacter import select_port, select_one_hydrogen, form_single_bond template_reacter = TemplateReacter( ... name="C-C_coupling",

... anchor_selector_left=select_port, ... anchor_selector_right=select_port, ... leaving_selector_left=select_one_hydrogen, ... leaving_selector_right=select_one_hydrogen, ... bond_former=form_single_bond, ... radius=4 ... ) >>> result, template = template_reacter.run_with_template( ... left, right, port_atom_L, port_atom_R ... ) >>> template.pre # Pre-reaction subgraph >>> template.post # Post-reaction subgraph

Initialize TemplateReacter.

Parameters:

Name	Type	Description	Default
`name`	`str`	Descriptive name for this reaction	required
`anchor_selector_left`	`AnchorSelector`	Selector that maps left port atom to anchor atom	required
`anchor_selector_right`	`AnchorSelector`	Selector that maps right port atom to anchor atom	required
`leaving_selector_left`	`LeavingSelector`	Selector for left leaving group (from left anchor)	required
`leaving_selector_right`	`LeavingSelector`	Selector for right leaving group (from right anchor)	required
`bond_former`	`BondFormer`	Function to create bond between anchor atoms	required
`radius`	`int`	Topological radius for subgraph extraction	`4`

Source code in src/molpy/reacter/template.py

def __init__(
    self,
    name: str,
    anchor_selector_left: AnchorSelector,
    anchor_selector_right: AnchorSelector,
    leaving_selector_left: LeavingSelector,
    leaving_selector_right: LeavingSelector,
    bond_former: BondFormer,
    radius: int = 4,
):
    """
    Initialize TemplateReacter.

    Args:
        name: Descriptive name for this reaction
        anchor_selector_left: Selector that maps left port atom to anchor atom
        anchor_selector_right: Selector that maps right port atom to anchor atom
        leaving_selector_left: Selector for left leaving group (from left anchor)
        leaving_selector_right: Selector for right leaving group (from right anchor)
        bond_former: Function to create bond between anchor atoms
        radius: Topological radius for subgraph extraction
    """
    # Create internal Reacter instance
    self.reacter = Reacter(
        name=name,
        anchor_selector_left=anchor_selector_left,
        anchor_selector_right=anchor_selector_right,
        leaving_selector_left=leaving_selector_left,
        leaving_selector_right=leaving_selector_right,
        bond_former=bond_former,
    )
    self.radius = radius
    self._react_id_counter = 0

run_with_template ¶

run_with_template(left, right, port_atom_L, port_atom_R, compute_topology=True, record_intermediates=False, typifier=None)

Run reaction and generate templates.

Function signature matches Reacter.run() but returns both ReactionResult and TemplateResult.

Parameters:

Name	Type	Description	Default
`left`	`Atomistic`	Left reactant structure	required
`right`	`Atomistic`	Right reactant structure	required
`port_atom_L`	`Entity`	Port atom in left structure	required
`port_atom_R`	`Entity`	Port atom in right structure	required
`compute_topology`	`bool`	If True, compute new angles/dihedrals (default True)	`True`
`record_intermediates`	`bool`	If True, record intermediate states	`False`
`typifier`	`TypifierBase \| None`	Optional typifier for incremental retypification	`None`

Returns:

Type	Description
`tuple[ReactionResult, TemplateResult]`	Tuple of (ReactionResult, TemplateResult)

Source code in src/molpy/reacter/template.py

def run_with_template(
    self,
    left: Atomistic,
    right: Atomistic,
    port_atom_L: Entity,
    port_atom_R: Entity,
    compute_topology: bool = True,
    record_intermediates: bool = False,
    typifier: TypifierBase | None = None,
) -> tuple[ReactionResult, TemplateResult]:
    """
    Run reaction and generate templates.

    Function signature matches Reacter.run() but returns both
    ReactionResult and TemplateResult.

    Args:
        left: Left reactant structure
        right: Right reactant structure
        port_atom_L: Port atom in left structure
        port_atom_R: Port atom in right structure
        compute_topology: If True, compute new angles/dihedrals (default True)
        record_intermediates: If True, record intermediate states
        typifier: Optional typifier for incremental retypification

    Returns:
        Tuple of (ReactionResult, TemplateResult)
    """
    # Step 1: Assign react_id to ALL atoms BEFORE reaction
    self._assign_react_ids(left)
    self._assign_react_ids(right)

    # Verify port atoms have react_id
    if "react_id" not in port_atom_L.data:
        raise ValueError("port_atom_L missing react_id after assignment!")
    if "react_id" not in port_atom_R.data:
        raise ValueError("port_atom_R missing react_id after assignment!")

    # Step 2: Run the underlying reaction
    result = self.reacter.run(
        left,
        right,
        port_atom_L,
        port_atom_R,
        compute_topology=compute_topology,
        record_intermediates=record_intermediates,
        typifier=typifier,
    )

    # Step 3: Generate templates from result
    template = self._generate_template(result, port_atom_L, port_atom_R, typifier)

    return result, template

TemplateResult `dataclass` ¶

TemplateResult(pre, post, init_atoms, edge_atoms, removed_atoms, pre_react_id_to_atom, post_react_id_to_atom)

Result of template generation.

write_template_files ¶

write_template_files(base_path, template, typifier=None)

Write LAMMPS fix bond/react template files.

Parameters:

Name	Type	Description	Default
`base_path`	`Path`	Base path for files (e.g., "rxn1" -> rxn1_pre.mol, rxn1_post.mol, rxn1.map)	required
`template`	`TemplateResult`	TemplateResult from TemplateReacter	required
`typifier`		Optional typifier to ensure types are set	`None`

Source code in src/molpy/reacter/template.py

def write_template_files(
    base_path: Path, template: TemplateResult, typifier=None
) -> None:
    """
    Write LAMMPS fix bond/react template files.

    Args:
        base_path: Base path for files (e.g., "rxn1" -> rxn1_pre.mol, rxn1_post.mol, rxn1.map)
        template: TemplateResult from TemplateReacter
        typifier: Optional typifier to ensure types are set
    """
    pre_path = Path(f"{base_path}_pre.mol")
    post_path = Path(f"{base_path}_post.mol")
    map_path = Path(f"{base_path}.map")

    pre = template.pre
    post = template.post

    # Note: Topology and typification should already be done in TemplateReacter._generate_template
    # This is just a safety check for backward compatibility
    if typifier:
        _ensure_typified(pre, typifier)
        _ensure_typified(post, typifier)

    # Build react_id to index mappings and set atom IDs
    pre_react_id_to_idx = {}
    for i, atom in enumerate(pre.atoms, start=1):
        if "react_id" not in atom.data:
            raise ValueError(f"Atom {i} in pre missing react_id!")
        pre_react_id_to_idx[atom["react_id"]] = i
        atom["id"] = i

    post_react_id_to_idx = {}
    for i, atom in enumerate(post.atoms, start=1):
        if "react_id" not in atom.data:
            raise ValueError(f"Atom {i} in post missing react_id!")
        post_react_id_to_idx[atom["react_id"]] = i
        atom["id"] = i

    # Verify pre and post have same atoms
    pre_react_ids = set(pre_react_id_to_idx.keys())
    post_react_ids = set(post_react_id_to_idx.keys())

    if pre_react_ids != post_react_ids:
        raise ValueError(
            f"Pre and post have different atoms!\n"
            f"  Missing in post: {pre_react_ids - post_react_ids}\n"
            f"  Missing in pre: {post_react_ids - pre_react_ids}"
        )

    # Build equivalences
    equiv = [
        (pre_react_id_to_idx[rid], post_react_id_to_idx[rid]) for rid in pre_react_ids
    ]

    # Get special atom indices
    initiator_ids = []
    for a in template.init_atoms:
        rid = a.get("react_id")
        if rid and rid in pre_react_id_to_idx:
            initiator_ids.append(pre_react_id_to_idx[rid])

    # Edge atoms: atoms in subgraph that have neighbors outside subgraph
    # Exclude initiator atoms from edge atoms (they are reaction sites, not boundaries)
    edge_ids = []
    initiator_rids = {a.get("react_id") for a in template.init_atoms}
    for a in template.edge_atoms:
        rid = a.get("react_id")
        # Don't mark initiator atoms as edge atoms
        if rid and rid in pre_react_id_to_idx and rid not in initiator_rids:
            edge_ids.append(pre_react_id_to_idx[rid])

    deleted_ids = []
    for a in template.removed_atoms:
        rid = a.get("react_id")
        # Mark all removed atoms for deletion, even if they are also initiators.
        # This is necessary for reactions like dehydration where:
        # - Initiator atoms are port atoms (two O atoms)
        # - Anchor atoms are the actual reaction sites (left: C neighbor, right: O port)
        # - Left port O is both an initiator AND part of the leaving group (O+H)
        # In such cases, the port O must be deleted even though it's an initiator
        if rid and rid in pre_react_id_to_idx:
            deleted_ids.append(pre_react_id_to_idx[rid])

    # Write map file
    with map_path.open("w", encoding="utf-8") as f:
        f.write("# auto-generated map file for fix bond/react\n\n")
        f.write(f"{len(equiv)} equivalences\n")
        f.write(f"{len(edge_ids)} edgeIDs\n")
        f.write(f"{len(deleted_ids)} deleteIDs\n\n")
        f.write("InitiatorIDs\n\n")
        for idx in initiator_ids:
            f.write(f"{idx}\n")
        f.write("\nEdgeIDs\n\n")
        for idx in edge_ids:
            f.write(f"{idx}\n")
        f.write("\nDeleteIDs\n\n")
        for idx in deleted_ids:
            f.write(f"{idx}\n")
        f.write("\nEquivalences\n\n")
        for pre_id, post_id in sorted(equiv):
            f.write(f"{pre_id}   {post_id}\n")

    # Write mol files with UNIFIED type mappings
    # CRITICAL: Pre and post must use the same type-to-ID mapping so LAMMPS
    # sees consistent types. Otherwise edge atoms may appear to have type changes.
    pre_frame = pre.to_frame()
    post_frame = post.to_frame()

    if "charge" in pre_frame["atoms"]:
        pre_frame["atoms"]["q"] = pre_frame["atoms"]["charge"]
    if "charge" in post_frame["atoms"]:
        post_frame["atoms"]["q"] = post_frame["atoms"]["charge"]

    # Create unified type mappings from BOTH pre and post
    _unify_type_mappings(pre_frame, post_frame)

    # Write with pre-converted types (won't trigger _convert_types_to_ids again)
    mp.io.write_lammps_molecule(pre_path, pre_frame)
    mp.io.write_lammps_molecule(post_path, post_frame)

    print(f"✅ Written: {pre_path.name}, {post_path.name}, {map_path.name}")

Topology Detector¶

Topology detection and update for chemical reactions.

This module provides intelligent detection and updating of angles and dihedrals after bond formation in reactions. It identifies affected old topology items and replaces them with new ones based on the updated bond structure.

TopologyDetector ¶

Detects and updates angles and dihedrals after reaction bond formation.

This class implements intelligent topology detection that: 1. Identifies atoms affected by new bonds 2. Removes old topology items (angles/dihedrals) involving affected atoms 3. Generates new topology items based on current bond structure

detect_and_update_topology `classmethod` ¶

detect_and_update_topology(assembly, new_bonds, removed_atoms)

Detect and update topology structure after reaction.

This method: 1. Removes topology items involving removed atoms 2. Generates new topology items based on new bonds 3. Adds new topology items (deduplicated with existing)

Note: We do NOT remove topology involving "affected atoms" (bond endpoints and neighbors) because those angles/dihedrals don't need to change - they still exist with the same atoms. We only need to add NEW topology created by the new bond.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	The Atomistic structure (will be modified)	required
`new_bonds`	`list[Bond]`	List of newly formed bonds	required
`removed_atoms`	`Collection[Entity]`	List of atoms that were removed from the structure	required

Returns:

Type	Description
`tuple[list[Angle], list[Dihedral], list[Angle], list[Dihedral]]`	Tuple of (new_angles, new_dihedrals, removed_angles, removed_dihedrals)

Source code in src/molpy/reacter/topology_detector.py

@classmethod
def detect_and_update_topology(
    cls,
    assembly: Atomistic,
    new_bonds: list[Bond],
    removed_atoms: Collection[Entity],
) -> tuple[list[Angle], list[Dihedral], list[Angle], list[Dihedral]]:
    """
    Detect and update topology structure after reaction.

    This method:
    1. Removes topology items involving removed atoms
    2. Generates new topology items based on new bonds
    3. Adds new topology items (deduplicated with existing)

    Note: We do NOT remove topology involving "affected atoms" (bond endpoints
    and neighbors) because those angles/dihedrals don't need to change -
    they still exist with the same atoms. We only need to add NEW topology
    created by the new bond.

    Args:
        assembly: The Atomistic structure (will be modified)
        new_bonds: List of newly formed bonds
        removed_atoms: List of atoms that were removed from the structure

    Returns:
        Tuple of (new_angles, new_dihedrals, removed_angles, removed_dihedrals)
    """
    # Step 1: Remove topology items involving removed atoms ONLY
    removed_angles, removed_dihedrals = cls._remove_topology_with_removed_atoms(
        assembly, removed_atoms
    )

    # Step 2: Generate new topology items based on new bonds
    # Get existing topology items for deduplication
    existing_angles = list(assembly.links.bucket(Angle))
    existing_dihedrals = list(assembly.links.bucket(Dihedral))

    new_angles_candidates = []
    new_dihedrals_candidates = []

    # Generate angles and dihedrals around each new bond
    for bond in new_bonds:
        # Angles around the bond
        angles_around_bond = cls._generate_angles_around_bond(assembly, bond)
        new_angles_candidates.extend(angles_around_bond)

        # Dihedrals through the bond
        dihedrals_through_bond = cls._generate_dihedrals_through_bond(
            assembly, bond
        )
        new_dihedrals_candidates.extend(dihedrals_through_bond)

        # Dihedrals continuing from the bond
        dihedrals_continuing = cls._generate_dihedrals_continuing_from_bond(
            assembly, bond
        )
        new_dihedrals_candidates.extend(dihedrals_continuing)

    # Step 3: Deduplicate and add new topology items
    unique_new_angles = cls._deduplicate_angles(
        new_angles_candidates, existing_angles
    )
    unique_new_dihedrals = cls._deduplicate_dihedrals(
        new_dihedrals_candidates, existing_dihedrals
    )

    # Add new items to assembly
    if unique_new_angles:
        assembly.add_link(*unique_new_angles)
    if unique_new_dihedrals:
        assembly.add_link(*unique_new_dihedrals)

    return (
        unique_new_angles,
        unique_new_dihedrals,
        removed_angles,
        removed_dihedrals,
    )

Transformers¶

Bond-making transformer functions for reactions.

Transformers create or modify bonds between anchor atoms in the product assembly.

break_bond ¶

break_bond(assembly, i, j)

Remove existing bond between two atoms.

Opposite of bond makers - used for bond-breaking reactions.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct containing the bond	required
`i`	`Entity`	First atom	required
`j`	`Entity`	Second atom	required

Side effects

Removes bond from assembly.links

Example

break_bond(assembly, carbon1, oxygen1)

Breaks C-O bond¶

Source code in src/molpy/reacter/transformers.py

def break_bond(assembly: Atomistic, i: Entity, j: Entity) -> None:
    """
    Remove existing bond between two atoms.

    Opposite of bond makers - used for bond-breaking reactions.

    Args:
        assembly: Struct containing the bond
        i: First atom
        j: Second atom

    Side effects:
        Removes bond from assembly.links

    Example:
        >>> break_bond(assembly, carbon1, oxygen1)
        >>> # Breaks C-O bond
    """
    existing = get_bond_between(assembly, i, j)
    if existing is not None:
        assembly.remove_link(existing)

create_bond_former ¶

create_bond_former(order)

Factory function to create bond former with specific order.

Parameters:

Name	Type	Description	Default
`order`	`int`	Bond order (1, 2, 3, or 1.5 for aromatic)	required

Returns:

Type	Description
`BondFormer`	BondFormer function that creates bonds with specified order

Example

double_bond_former = create_bond_former(2) reacter = Reacter( ... bond_former=double_bond_former, ... ... ... )

Source code in src/molpy/reacter/transformers.py

def create_bond_former(order: int) -> BondFormer:
    """
    Factory function to create bond former with specific order.

    Args:
        order: Bond order (1, 2, 3, or 1.5 for aromatic)

    Returns:
        BondFormer function that creates bonds with specified order

    Example:
        >>> double_bond_former = create_bond_former(2)
        >>> reacter = Reacter(
        ...     bond_former=double_bond_former,
        ...     ...
        ... )
    """

    def bond_former(assembly: Atomistic, i: Entity, j: Entity) -> Bond | None:
        existing = get_bond_between(assembly, i, j)

        # Determine kind symbol
        kind_map = {1: "-", 2: "=", 3: "#", 1.5: ":"}
        kind = kind_map.get(order, "-")

        if existing is not None:
            existing["order"] = order
            existing["kind"] = kind
            if order == 1.5:
                existing["aromatic"] = True
            return existing
        else:
            attrs = {"order": order, "kind": kind}
            if order == 1.5:
                attrs["aromatic"] = True
            bond = Bond(cast(Atom, i), cast(Atom, j), **attrs)
            assembly.add_link(bond, include_endpoints=False)
            return bond

    return bond_former

form_aromatic_bond ¶

form_aromatic_bond(assembly, i, j)

Create an aromatic bond between two atoms.

If a bond already exists, updates it to aromatic (order=1.5 by convention).

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct to add bond to	required
`i`	`Entity`	First atom	required
`j`	`Entity`	Second atom	required

Returns:

Type	Description
`Bond \| None`	The created or updated Bond object

Side effects

Adds Bond(i, j, order=1.5, kind=':') to assembly.links

Source code in src/molpy/reacter/transformers.py

def form_aromatic_bond(assembly: Atomistic, i: Entity, j: Entity) -> Bond | None:
    """
    Create an aromatic bond between two atoms.

    If a bond already exists, updates it to aromatic (order=1.5 by convention).

    Args:
        assembly: Struct to add bond to
        i: First atom
        j: Second atom

    Returns:
        The created or updated Bond object

    Side effects:
        Adds Bond(i, j, order=1.5, kind=':') to assembly.links
    """
    existing = get_bond_between(assembly, i, j)

    if existing is not None:
        existing["order"] = 1.5
        existing["kind"] = ":"
        existing["aromatic"] = True
        return existing
    else:
        bond = Bond(cast(Atom, i), cast(Atom, j), order=1.5, kind=":", aromatic=True)
        assembly.add_link(bond, include_endpoints=False)
        return bond

form_double_bond ¶

form_double_bond(assembly, i, j)

Create a double bond between two atoms.

If a bond already exists, updates it to double bond (order=2).

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct to add bond to	required
`i`	`Entity`	First atom	required
`j`	`Entity`	Second atom	required

Returns:

Type	Description
`Bond \| None`	The created or updated Bond object

Side effects

Adds Bond(i, j, order=2) to assembly.links

Source code in src/molpy/reacter/transformers.py

def form_double_bond(assembly: Atomistic, i: Entity, j: Entity) -> Bond | None:
    """
    Create a double bond between two atoms.

    If a bond already exists, updates it to double bond (order=2).

    Args:
        assembly: Struct to add bond to
        i: First atom
        j: Second atom

    Returns:
        The created or updated Bond object

    Side effects:
        Adds Bond(i, j, order=2) to assembly.links
    """
    existing = get_bond_between(assembly, i, j)

    if existing is not None:
        existing["order"] = 2
        existing["kind"] = "="
        return existing
    else:
        bond = Bond(cast(Atom, i), cast(Atom, j), order=2, kind="=")
        assembly.add_link(bond, include_endpoints=False)
        return bond

form_single_bond ¶

form_single_bond(assembly, i, j)

Create a single bond between two atoms.

If a bond already exists, updates it to single bond (order=1).

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct to add bond to	required
`i`	`Entity`	First atom	required
`j`	`Entity`	Second atom	required

Returns:

Type	Description
`Bond \| None`	The created or updated Bond object

Side effects

Adds Bond(i, j, order=1) to assembly.links

Example

bond = form_single_bond(merged, carbon1, carbon2)

Source code in src/molpy/reacter/transformers.py

def form_single_bond(assembly: Atomistic, i: Entity, j: Entity) -> Bond | None:
    """
    Create a single bond between two atoms.

    If a bond already exists, updates it to single bond (order=1).

    Args:
        assembly: Struct to add bond to
        i: First atom
        j: Second atom

    Returns:
        The created or updated Bond object

    Side effects:
        Adds Bond(i, j, order=1) to assembly.links

    Example:
        >>> bond = form_single_bond(merged, carbon1, carbon2)
    """
    # Avoid creating self-bonds (loop edges) which break topology algorithms
    if i is j:
        return None

    # Check if bond exists
    existing = get_bond_between(assembly, i, j)

    if existing is not None:
        # Update existing bond
        existing["order"] = 1
        existing["kind"] = "-"
        return existing
    else:
        # Create new bond
        bond = Bond(cast(Atom, i), cast(Atom, j), order=1, kind="-")
        assembly.add_link(bond, include_endpoints=False)
        return bond

form_triple_bond ¶

form_triple_bond(assembly, i, j)

Create a triple bond between two atoms.

If a bond already exists, updates it to triple bond (order=3).

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct to add bond to	required
`i`	`Entity`	First atom	required
`j`	`Entity`	Second atom	required

Returns:

Type	Description
`Bond \| None`	The created or updated Bond object

Side effects

Adds Bond(i, j, order=3) to assembly.links

Source code in src/molpy/reacter/transformers.py

def form_triple_bond(assembly: Atomistic, i: Entity, j: Entity) -> Bond | None:
    """
    Create a triple bond between two atoms.

    If a bond already exists, updates it to triple bond (order=3).

    Args:
        assembly: Struct to add bond to
        i: First atom
        j: Second atom

    Returns:
        The created or updated Bond object

    Side effects:
        Adds Bond(i, j, order=3) to assembly.links
    """
    existing = get_bond_between(assembly, i, j)

    if existing is not None:
        existing["order"] = 3
        existing["kind"] = "#"
        return existing
    else:
        bond = Bond(cast(Atom, i), cast(Atom, j), order=3, kind="#")
        assembly.add_link(bond, include_endpoints=False)
        return bond

skip_bond_formation ¶

skip_bond_formation(assembly, i, j)

Do not create any bond.

Useful for reactions that only remove atoms without forming new bonds.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Struct (ignored)	required
`i`	`Entity`	First atom (ignored)	required
`j`	`Entity`	Second atom (ignored)	required

Side effects

None

Example

reacter = Reacter( ... bond_former=skip_bond_formation, # Just remove leaving groups ... ... ... )

Source code in src/molpy/reacter/transformers.py

def skip_bond_formation(assembly: Atomistic, i: Entity, j: Entity) -> None:
    """
    Do not create any bond.

    Useful for reactions that only remove atoms without forming new bonds.

    Args:
        assembly: Struct (ignored)
        i: First atom (ignored)
        j: Second atom (ignored)

    Side effects:
        None

    Example:
        >>> reacter = Reacter(
        ...     bond_former=skip_bond_formation,  # Just remove leaving groups
        ...     ...
        ... )
    """
    return None

Utils¶

Utility functions for assembly manipulation in reactions.

This module provides helper functions for finding neighbors and other common operations needed for reactions.

count_bonds ¶

count_bonds(assembly, atom)

Count the number of bonds connected to an atom.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Atomistic assembly containing the atom	required
`atom`	`Entity`	Atom entity to count bonds for	required

Returns:

Type	Description
`int`	Number of bonds

Example

valence = count_bonds(asm, carbon_atom) print(f"Carbon has {valence} bonds")

Source code in src/molpy/reacter/utils.py

def count_bonds(assembly: Atomistic, atom: Entity) -> int:
    """
    Count the number of bonds connected to an atom.

    Args:
        assembly: Atomistic assembly containing the atom
        atom: Atom entity to count bonds for

    Returns:
        Number of bonds

    Example:
        >>> valence = count_bonds(asm, carbon_atom)
        >>> print(f"Carbon has {valence} bonds")
    """
    count = 0
    for bond in assembly.links.bucket(Bond):
        # Use identity check (is) not equality check (==)
        if any(ep is atom for ep in bond.endpoints):
            count += 1
    return count

create_atom_mapping ¶

create_atom_mapping(pre_atoms, post_atoms)

Create atom mapping between pre-reaction and post-reaction states.

Creates a mapping of atom IDs from pre-reaction template to post-reaction template. This is used to generate map files for LAMMPS fix bond/react.

Parameters:

Name	Type	Description	Default
`pre_atoms`	`list`	List of atoms in pre-reaction state	required
`post_atoms`	`list`	List of atoms in post-reaction state	required

Returns:

Type	Description
`dict[int, int]`	Dictionary mapping pre-reaction atom indices to post-reaction indices
`dict[int, int]`	(1-indexed for LAMMPS)

Note

Atoms that are deleted during the reaction will not appear in the mapping. Atoms are matched by their 'id' attribute if present, otherwise by position.

Example

from molpy.reacter.utils import create_atom_mapping mapping = create_atom_mapping(pre_atoms, post_atoms)

Write to map file¶

with open("reaction.map", 'w') as f: ... for pre_id, post_id in sorted(mapping.items()): ... f.write(f"{pre_id} {post_id}\n")

Source code in src/molpy/reacter/utils.py

def create_atom_mapping(pre_atoms: list, post_atoms: list) -> dict[int, int]:
    """Create atom mapping between pre-reaction and post-reaction states.

    Creates a mapping of atom IDs from pre-reaction template to post-reaction
    template. This is used to generate map files for LAMMPS fix bond/react.

    Args:
        pre_atoms: List of atoms in pre-reaction state
        post_atoms: List of atoms in post-reaction state

    Returns:
        Dictionary mapping pre-reaction atom indices to post-reaction indices
        (1-indexed for LAMMPS)

    Note:
        Atoms that are deleted during the reaction will not appear in the mapping.
        Atoms are matched by their 'id' attribute if present, otherwise by position.

    Example:
        >>> from molpy.reacter.utils import create_atom_mapping
        >>> mapping = create_atom_mapping(pre_atoms, post_atoms)
        >>> # Write to map file
        >>> with open("reaction.map", 'w') as f:
        ...     for pre_id, post_id in sorted(mapping.items()):
        ...         f.write(f"{pre_id} {post_id}\\n")
    """
    mapping = {}

    # Try matching by atom 'id' first
    pre_id_map = {}
    for i, atom in enumerate(pre_atoms):
        atom_id = atom.get("id")
        if atom_id is not None:
            pre_id_map[atom_id] = i + 1

    post_id_map = {}
    for i, atom in enumerate(post_atoms):
        atom_id = atom.get("id")
        if atom_id is not None:
            post_id_map[atom_id] = i + 1

    # Match atoms by 'id' attribute
    for atom_id, pre_idx in pre_id_map.items():
        if atom_id in post_id_map:
            mapping[pre_idx] = post_id_map[atom_id]

    # If no matches by ID, try matching by position (for atoms without IDs)
    if not mapping:
        import numpy as np

        for i, pre_atom in enumerate(pre_atoms):
            pre_pos = np.array(
                [
                    pre_atom["x"],
                    pre_atom["y"],
                    pre_atom["z"],
                ]
            )

            for j, post_atom in enumerate(post_atoms):
                post_pos = np.array(
                    [
                        post_atom["x"],
                        post_atom["y"],
                        post_atom["z"],
                    ]
                )

                # Check if positions are close (within 0.01 Angstrom)
                dist_sq = np.sum((pre_pos - post_pos) ** 2)
                if dist_sq < 1e-4:
                    mapping[i + 1] = j + 1
                    break

    return mapping

find_neighbors ¶

find_neighbors(assembly, atom, *, element=None)

Find neighboring atoms of a given atom.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Atomistic assembly containing the atom	required
`atom`	`Entity`	Atom entity to find neighbors of	required
`element`	`str \| None`	Optional element symbol to filter by (e.g., 'H', 'C')	`None`

Returns:

Type	Description
`list[Entity]`	List of neighboring atom entities

Example

h_neighbors = find_neighbors(asm, carbon_atom, element='H') all_neighbors = find_neighbors(asm, carbon_atom)

Source code in src/molpy/reacter/utils.py

def find_neighbors(
    assembly: Atomistic,
    atom: Entity,
    *,
    element: str | None = None,
) -> list[Entity]:
    """
    Find neighboring atoms of a given atom.

    Args:
        assembly: Atomistic assembly containing the atom
        atom: Atom entity to find neighbors of
        element: Optional element symbol to filter by (e.g., 'H', 'C')

    Returns:
        List of neighboring atom entities

    Example:
        >>> h_neighbors = find_neighbors(asm, carbon_atom, element='H')
        >>> all_neighbors = find_neighbors(asm, carbon_atom)
    """
    neighbors: list[Entity] = []

    # Look through all bonds
    for bond in assembly.links.bucket(Bond):
        # Use identity check (is) not equality check (==)
        if any(ep is atom for ep in bond.endpoints):
            # Found a bond involving this atom
            for endpoint in bond.endpoints:
                if endpoint is not atom:
                    # Filter by element if specified
                    if element is None or endpoint.get("symbol") == element:
                        neighbors.append(endpoint)

    return neighbors

get_bond_between ¶

get_bond_between(assembly, i, j)

Find existing bond between two atoms.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Atomistic assembly containing the atoms	required
`i`	`Entity`	First atom entity	required
`j`	`Entity`	Second atom entity	required

Returns:

Type	Description
`Bond \| None`	Bond entity if found, None otherwise

Example

bond = get_bond_between(asm, itom, jtom) if bond: ... print(f"Bond order: {bond.get('order', 1)}")

Source code in src/molpy/reacter/utils.py

def get_bond_between(
    assembly: Atomistic,
    i: Entity,
    j: Entity,
) -> Bond | None:
    """
    Find existing bond between two atoms.

    Args:
        assembly: Atomistic assembly containing the atoms
        i: First atom entity
        j: Second atom entity

    Returns:
        Bond entity if found, None otherwise

    Example:
        >>> bond = get_bond_between(asm, itom, jtom)
        >>> if bond:
        ...     print(f"Bond order: {bond.get('order', 1)}")
    """
    for bond in assembly.links.bucket(Bond):
        endpoints = bond.endpoints
        # Use identity check (is) not equality check (==)
        if any(ep is i for ep in endpoints) and any(ep is j for ep in endpoints):
            return bond
    return None

remove_dummy_atoms ¶

remove_dummy_atoms(assembly)

Remove all dummy atoms (element '' or symbol '') from assembly.

Parameters:

Name	Type	Description	Default
`assembly`	`Atomistic`	Atomistic assembly to clean	required

Returns:

Type	Description
`list[Entity]`	List of removed dummy atoms

Example

removed = remove_dummy_atoms(merged_asm) print(f"Removed {len(removed)} dummy atoms")

Source code in src/molpy/reacter/utils.py

def remove_dummy_atoms(assembly: Atomistic) -> list[Entity]:
    """
    Remove all dummy atoms (element '*' or symbol '*') from assembly.

    Args:
        assembly: Atomistic assembly to clean

    Returns:
        List of removed dummy atoms

    Example:
        >>> removed = remove_dummy_atoms(merged_asm)
        >>> print(f"Removed {len(removed)} dummy atoms")
    """
    dummy_atoms: list[Entity] = []

    for atom in assembly.entities.bucket(Atom):
        symbol = atom.get("symbol", "")
        element = atom.get("element", "")
        if symbol == "*" or element == "*":
            dummy_atoms.append(atom)

    if dummy_atoms:
        assembly.remove_entity(*dummy_atoms, drop_incident_links=True)

    return dummy_atoms

Reacter¶

Base¶

AnchorSelector module-attribute ¶

BondFormer module-attribute ¶

LeavingSelector module-attribute ¶

ProductInfo dataclass ¶

ReactantInfo dataclass ¶

Reacter ¶

Mark ports on atoms¶

Port selection is explicit!¶

run ¶

ReactionMetadata dataclass ¶

ReactionResult dataclass ¶

TopologyChanges dataclass ¶

Connector¶

MonomerLinker ¶

Mark ports on atoms¶

Default reaction for most connections¶

Special reaction for A-B connection¶

Explicit port specification required¶

clear_history ¶

connect ¶

get_all_modified_atoms ¶

get_history ¶

needs_retypification ¶

Selectors¶

find_port_atom ¶

find_port_atom_by_node ¶

select_all_hydrogens ¶

select_c_neighbor ¶

select_dehydration_left ¶

select_dehydration_right ¶

select_dummy_atoms ¶

select_hydroxyl_group ¶

select_hydroxyl_h_only ¶

select_none ¶

select_o_neighbor ¶

select_one_hydrogen ¶

select_port ¶

Template¶

TemplateReacter ¶

run_with_template ¶

TemplateResult dataclass ¶

write_template_files ¶

Topology Detector¶

TopologyDetector ¶

detect_and_update_topology classmethod ¶

Transformers¶

break_bond ¶

Breaks C-O bond¶

create_bond_former ¶

form_aromatic_bond ¶

form_double_bond ¶

form_single_bond ¶

form_triple_bond ¶

skip_bond_formation ¶

Utils¶

count_bonds ¶

create_atom_mapping ¶

Write to map file¶

find_neighbors ¶

get_bond_between ¶

remove_dummy_atoms ¶

AnchorSelector `module-attribute` ¶

BondFormer `module-attribute` ¶

LeavingSelector `module-attribute` ¶

ProductInfo `dataclass` ¶

ReactantInfo `dataclass` ¶

ReactionMetadata `dataclass` ¶

ReactionResult `dataclass` ¶

TopologyChanges `dataclass` ¶

TemplateResult `dataclass` ¶

detect_and_update_topology `classmethod` ¶