Skip to content

Reacter

Chemical reaction modeling: anchor selection, leaving group removal, bond formation, template generation.

Quick reference

Symbol Summary Preferred for
Reacter Execute a reaction between two structures Polymer growth, bond formation
BondReactReacter Reacter + subgraph extraction for LAMMPS fix bond/react Reactive MD templates
find_port(struct, name) Find first atom with port marker Locating </>/$ ports
select_self Use port atom itself as reaction anchor Right-side anchor selection
select_neighbor(element) Select neighbor of given element Left-side anchor selection
select_hydrogens(n) Select n hydrogen leaving atoms Dehydration reactions
select_none Return empty list (no leaving group) Addition reactions
form_single_bond Create single bond between atoms Standard bond formation

Selector protocol

  • Anchor selectors: (struct: Atomistic, port_atom: Atom) -> Atom — map the port atom to the atom that actually forms the bond.
  • Leaving selectors: (struct: Atomistic, anchor: Atom) -> list[Atom] — return zero or more atoms to remove.

Canonical example

from molpy.core.atomistic import Atom, Atomistic, Bond
from molpy.reacter import (
    Reacter,
    find_port,
    form_single_bond,
    select_hydrogens,
    select_self,
)


def methane_fragment(port):
    struct = Atomistic()
    carbon = Atom(element="C")
    hydrogens = [Atom(element="H") for _ in range(3)]
    struct.add_entity(carbon, *hydrogens)
    for hydrogen in hydrogens:
        struct.add_link(Bond(carbon, hydrogen))
    carbon["port"] = port
    return struct


left = methane_fragment(">")
right = methane_fragment("<")

rxn = Reacter(
    name="cc_coupling",
    anchor_selector_left=select_self,
    anchor_selector_right=select_self,
    leaving_selector_left=select_hydrogens(1),
    leaving_selector_right=select_hydrogens(1),
    bond_former=form_single_bond,
)
result = rxn.run(
    left=left,
    right=right,
    port_atom_L=find_port(left, ">"),
    port_atom_R=find_port(right, "<"),
    compute_topology=True,
)
product = result.product
assert len(result.removed_atoms) == 2
assert len(list(product.atoms)) == 6

Reacter.run() never mutates the caller's inputs — the reaction runs on internal copies and the product is a fresh structure.


Full API

Base

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.

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 must explicitly specify which port atoms to connect via port_atom_L and port_atom_R (locate them with :func:molpy.reacter.find_port). 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.core.atomistic import Atom, Atomistic, Bond from molpy.reacter import ( ... Reacter, ... find_port, ... form_single_bond, ... select_hydrogens, ... select_self, ... ) def methane_fragment(port): ... struct = Atomistic() ... carbon = Atom(element="C") ... hydrogens = [Atom(element="H") for _ in range(3)] ... struct.add_entity(carbon, *hydrogens) ... for hydrogen in hydrogens: ... struct.add_link(Bond(carbon, hydrogen)) ... carbon["port"] = port ... return struct left = methane_fragment(">") right = methane_fragment("<") cc_coupling = Reacter( ... name="C-C_coupling_with_H_loss", ... anchor_selector_left=select_self, ... anchor_selector_right=select_self, ... leaving_selector_left=select_hydrogens(1), ... leaving_selector_right=select_hydrogens(1), ... bond_former=form_single_bond, ... ) result = cc_coupling.run( ... left, ... right, ... port_atom_L=find_port(left, ">"), ... port_atom_R=find_port(right, "<"), ... ) len(result.removed_atoms) 2 len(list(result.product.atoms)) # 8 atoms in, 2 H removed 6

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
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() 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/impropers 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/impropers (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

Copy semantics

Caller-owned left/right are copied once each and never mutated. With record_intermediates=False the merged assembly is not copied at all by the base Reacter (result.reactants is None); subclasses with _needs_reactants_snapshot = True (e.g. BondReactReacter) take exactly one pre-reaction snapshot so templates can be generated. record_intermediates=True adds copies for the recorded intermediate states.

ReactionResult dataclass

ReactionResult(product, anchor_L=None, anchor_R=None, reactants=None, port_atom_L=None, port_atom_R=None, new_bonds=list(), new_angles=list(), new_dihedrals=list(), new_impropers=list(), removed_angles=list(), removed_dihedrals=list(), removed_impropers=list(), removed_atoms=list(), modified_atoms=set(), reaction_name='', requires_retype=False, entity_maps=list(), intermediates=list())

Result of a reaction execution.

Flat structure containing the product, reactant snapshot, topology changes, and bookkeeping fields.

Selectors

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:

  1. Port Selectors (port_selector_left, port_selector_right):
  2. Signature: (struct: Atomistic, port_atom: Entity) -> Entity
  3. Transform a port atom to its actual reaction site
  4. Example: O port → C neighbor for dehydration

  5. Leaving Selectors (leaving_selector_left, leaving_selector_right):

  6. Signature: (struct: Atomistic, reaction_site: Entity) -> list[Entity]
  7. Identify atoms to remove from the reaction site
  8. 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

find_port(struct, port_name, *, node_id=None)

Find an atom carrying the given port marker.

This is a utility function, not a selector. Use it to locate port atoms ("<" / ">" / custom labels) before passing them to :meth:molpy.reacter.Reacter.run. Matches either the scalar port attribute or membership in a ports list.

Parameters:

Name Type Description Default
struct Atomistic

Atomistic structure containing ports.

required
port_name str

Name of port to find.

required
node_id int | None

If given, restrict the search to atoms whose monomer_node_id equals this value (useful when a structure holds several monomer units).

None

Returns:

Type Description
Atom

Atom with matching port.

Raises:

Type Description
ValueError

If no atom carries the port (for the node, when node_id is given).

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

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

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

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

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

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='*')

select_hydrogens

select_hydrogens(n=None)

Factory that returns a leaving selector picking n hydrogen neighbors.

Parameters:

Name Type Description Default
n int | None

Number of hydrogens to select. None means all hydrogens bonded to the site.

None

Returns:

Type Description
Callable[[Atomistic, Atom], list[Atom]]

A leaving selector (struct, site) -> list[Atom]

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

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 []

select_neighbor

select_neighbor(element)

Factory that returns an anchor selector picking a neighbor of a given element.

Parameters:

Name Type Description Default
element str

Element symbol to look for (e.g. "C", "O").

required

Returns:

Type Description
Callable[[Atomistic, Atom], Atom]

An anchor selector (struct, port_atom) -> Atom

Raises:

Type Description
ValueError

If no neighbor with the requested element is found.

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

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

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

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

select_self

select_self(struct, port_atom)

Identity anchor selector -- returns the port atom as the reaction site.

This is the simplest anchor selector: the port atom itself is where the new bond will form.

Parameters:

Name Type Description Default
struct Atomistic

Atomistic structure (unused, kept for interface compatibility)

required
port_atom Atom

The port atom entity

required

Returns:

Type Description
Atom

The same port atom

Topology Detector

topology_detector

Topology detection and update for chemical reactions.

This module provides intelligent detection and updating of angles, dihedrals, and impropers 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 bonded topology after reaction bond formation.

Implements the topology bookkeeping required by REACTER-style template generation (Gissinger et al., Polymer 128 (2017) 211-217, DOI: 10.1016/j.polymer.2017.06.038; Gissinger et al., Macromolecules 53 (2020) 9953-9961, DOI: 10.1021/acs.macromol.0c02012; LAMMPS fix bond/react: https://docs.lammps.org/fix_bond_react.html):

  1. Identifies atoms affected by new bonds (endpoints + first shell)
  2. Removes topology items (angles/dihedrals/impropers) involving removed atoms
  3. Generates new angles/dihedrals around new bonds and improper candidates at affected sp2-like centers (exactly 3 bonded neighbors), deduplicated against existing topology

Impropers matter physically: OPLS-AA and GAFF use them to keep sp2 centers planar, so a post-reaction template that lost impropers would let planar groups pyramidalize.

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 (angles/dihedrals/impropers) involving removed atoms 2. Generates new angles/dihedrals around new bonds and improper candidates at affected atoms with exactly 3 bonded neighbors 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
list[Angle]

Tuple of (new_angles, new_dihedrals, new_impropers,

list[Dihedral]

removed_angles, removed_dihedrals, removed_impropers)

Utils

utils

Utility functions and type aliases for reactions.

AdjacencyMap module-attribute

AdjacencyMap = dict[Entity, list[tuple[Entity, Bond]]]

Atom → list of (neighbor, connecting bond) pairs; see build_adjacency.

AnchorSelector module-attribute

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

Select anchor atom given a port atom.

BondFormer module-attribute

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

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

LeavingSelector module-attribute

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

Select leaving group atoms given an anchor atom.

break_bond

break_bond(assembly, i, j)

Remove existing bond between two atoms.

build_adjacency

build_adjacency(assembly)

Build an atom → neighbors adjacency map in one pass over bonds.

Complexity: O(bonds) to build; afterwards neighbor queries through :func:find_neighbors / :func:get_bond_between / :func:count_bonds with adjacency= are O(degree) instead of O(bonds) per call.

Parameters:

Name Type Description Default
assembly Atomistic

Atomistic assembly to index.

required

Returns:

Type Description
AdjacencyMap

Mapping from each bonded atom to its (neighbor, bond) pairs.

AdjacencyMap

Atoms without bonds are absent (queries fall back to empty).

count_bonds

count_bonds(assembly, atom, *, adjacency=None)

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
adjacency AdjacencyMap | None

Optional prebuilt adjacency map for O(degree) lookup

None

Returns:

Type Description
int

Number of bonds

Example

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

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")

create_bond_former

create_bond_former(order)

Factory function to create bond former with specific order.

find_neighbors

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

Find neighboring atoms of a given atom.

Complexity: O(bonds) full scan by default; O(degree) when a prebuilt adjacency map (see :func:build_adjacency) is given.

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
adjacency AdjacencyMap | None

Optional prebuilt adjacency map for O(degree) lookup

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)

form_aromatic_bond

form_aromatic_bond(assembly, i, j)

Create an aromatic bond between two atoms.

form_double_bond

form_double_bond(assembly, i, j)

Create a double bond between two atoms.

form_single_bond

form_single_bond(assembly, i, j)

Create a single bond between two atoms.

form_triple_bond

form_triple_bond(assembly, i, j)

Create a triple bond between two atoms.

get_bond_between

get_bond_between(assembly, i, j, *, adjacency=None)

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
adjacency AdjacencyMap | None

Optional prebuilt adjacency map for O(degree) lookup

None

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)}")

remove_dummy_atoms

remove_dummy_atoms(assembly)

Remove all dummy atoms (element '*') 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")

skip_bond_formation

skip_bond_formation(assembly, i, j)

Do not create any bond. For reactions that only remove atoms.