qrunch.chemistry.embedded_atoms.embedded_atoms_selector

Active site selection based on coordination environment changes and partial charge changes during a reaction.

Classes

AromaticRingExtender

Extends the embedded atoms to include the complete aromatic ring whenever any ring atom is active.

BondLengthChangeSelector

Selects atoms involved in the bonds whose length changes most across reaction images.

CoordinationChangeSelector

Selects atoms whose coordination environment changes across reaction images.

DistanceExtender

Extends the embedded atoms to include atoms within a distance threshold from any active atom.

ElectronicStructureChangeSelector

Selects atoms whose local electronic structure changes significantly across reaction images.

EmbeddedAtomsConnectionExtender

Chemically "repairs" the active selection by connecting fragments into complete pieces.

EmbeddedAtomsSelector

Composes primary embedded atoms selectors and extenders to identify active atoms in a reaction.

EmbeddedAtomsSelectorOptions

Options controlling EmbeddedAtomsSelector.

FullMoleculeExtender

Extends the embedded atoms to include all atoms in connected fragments containing active atoms.

HydrogenBondChangeSelector

Selects atoms in hydrogen bonds whose geometry changes most across reaction images.

NearestNeighborExtender

Extends the embedded atoms to include atoms directly bonded to any active atom in any image.
class AromaticRingExtender

Bases: EmbeddedAtomsExtender

Extends the embedded atoms to include the complete aromatic ring whenever any ring atom is active.

If the current selection includes part of an aromatic ring (e.g. a benzene ring), this extender adds the remaining atoms of that ring, so that the whole conjugated ring system is embedded together. Only the ring atoms themselves are added - hydrogens and side chains attached to the ring are not included.

Aromatic rings are detected geometrically: a ring is treated as aromatic when it is planar and composed of aromatic-eligible elements (see _detect_aromatic_rings()). This handles benzene as well as heteroaromatic rings (e.g. pyridine, furan) without needing to perceive bond orders. A ring is included if any of its atoms is active in any reaction image.

__init__(ring_sizes: tuple[int, ...] = (5, 6), planarity_tolerance: float = 0.1) → None
Parameters:

  • ring_sizes (tuple[int, ...])

  • planarity_tolerance (float)

Return type:

None

extend(active_indices: set[int], reaction: Reaction) → set[int]

Return the atom indices completing any aromatic ring that overlaps the active set.

Parameters:

  • active_indices (set[int]) – The current set of active atom indices.

  • reaction (Reaction) – The reaction being analyzed.

Return type:

set[int]

planarity_tolerance: float = 0.1

Maximum out-of-plane deviation (in Angstrom) for a ring to be considered planar/aromatic.

ring_sizes: tuple[int, ...] = (5, 6)

5 and 6).

Type:

Ring sizes (number of atoms) to treat as candidate aromatic rings (default

class BondLengthChangeSelector

Bases: PrimaryEmbeddedAtomsSelector

Selects atoms involved in the bonds whose length changes most across reaction images.

For each bond (taken as the union of the connectivity graphs of all images) the change in length is the spread max - min of its length across the images. Each atom is scored by the largest change among its incident bonds, and atoms are separated into active and inactive groups either by an unsupervised clustering algorithm (Jenks natural breaks) or a user-provided threshold.

This is a purely geometric selector: it is invariant to rigid-body motion and, unlike the per-atom electronic descriptor used by ElectronicStructureChangeSelector, it detects symmetric homonuclear bond stretches (e.g. H2 dissociation). The two selectors are intended to be combined for a robust selection.

__init__(threshold: float | None = None) → None
Parameters:

threshold (float | None)

Return type:

None

scores(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → ndarray[tuple[Any, ...], dtype[float64]]

Score atoms by their largest incident bond-length change.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[float64]]

select(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → set[int]

Select atoms involved in the bonds whose length changes most across reaction images.

Parameters:
Return type:

set[int]

threshold: float | None = None

Optional explicit threshold (in Angstrom) for active atom selection.

When None (default), the threshold is determined automatically using Jenks natural breaks clustering. When set, atoms whose largest incident bond-length change is at or above this value are considered active.

class CoordinationChangeSelector

Bases: PrimaryEmbeddedAtomsSelector

Selects atoms whose coordination environment changes across reaction images.

Builds a connectivity graph for each image using RDKit’s hybridization-aware bond determination and detects atoms whose coordination environment changes between any pair of images. The coordination environment is the sorted tuple of atomic numbers of an atom’s neighbors, so it distinguishes both changes in coordination number and changes in the types of bonded atoms (e.g., a hydrogen switching from being bonded to another hydrogen to being bonded to a carbon).

__init__() → None
Return type:

None

scores(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → ndarray[tuple[Any, ...], dtype[float64]]

Score atoms by whether their coordination environment changed.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[float64]]

select(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → set[int]

Select atoms whose coordination environment changed across the reaction images.

Parameters:
Return type:

set[int]

class DistanceExtender

Bases: EmbeddedAtomsExtender

Extends the embedded atoms to include atoms within a distance threshold from any active atom.

__init__(distance_threshold: float) → None
Parameters:

distance_threshold (float)

Return type:

None

distance_threshold: float

Maximum distance (in Angstrom) for inclusion.

extend(active_indices: set[int], reaction: Reaction) → set[int]

Return atom indices within the distance threshold from active atoms.

Parameters:

  • active_indices (set[int]) – The current set of active atom indices.

  • reaction (Reaction) – The reaction being analyzed.

Return type:

set[int]

class ElectronicStructureChangeSelector

Bases: PrimaryEmbeddedAtomsSelector

Selects atoms whose local electronic structure changes significantly across reaction images.

For each image the total electronic density is symmetrically (Lowdin) orthogonalized, \(D' = S^{1/2} D S^{1/2}\), and partitioned into atom-local diagonal blocks in the full AO basis. Each atom is described by the sorted eigenvalues of its block, i.e. its local Lowdin orbital occupations. Atoms are ranked by the maximum change of this descriptor across all pairs of images.

The descriptor is invariant under rigid-body translation and rotation of any fragment (an atom’s block transforms as an orthogonal similarity \(B \\to O B O^T\), leaving its eigenvalues unchanged), so spectator fragments that merely move between images score near zero. Genuine electronic rearrangements - charge transfer, changes in hybridization, polarization, or bond order that redistribute local occupation - do change the eigenvalues and are detected. This directly targets local electronic-structure change without relying on connectivity, which is handled separately by CoordinationChangeSelector.

Either an unsupervised clustering algorithm (Jenks natural breaks) or a user-provided threshold separates atoms into active and inactive groups.

Supports both restricted and unrestricted prepared data. For unrestricted calculations, the total (alpha + beta) density matrix is used.

Notes that this selector cannot detect homonuclear bond stretches (e.g. H2 dissociation) because the local occupation spectra of the two atoms remain identical.

Requires prepared_data and molecular_configurations to be provided to select(), or they will be automatically generated using a minimal basis (STO-3G) unrestricted Hartree-Fock calculation followed by Pipek-Mezey localization.

__init__(threshold: float | None = None) → None
Parameters:

threshold (float | None)

Return type:

None

scores(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → ndarray[tuple[Any, ...], dtype[float64]]

Score atoms by their local electronic-structure change.

Parameters:

  • reaction (Reaction) – The reaction to analyze.

  • prepared_data (Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None) – Sequence of prepared data (one per reaction image) from the mean-field and localization stages. Both restricted and unrestricted data are supported.

  • molecular_configurations (Sequence[MolecularConfiguration] | None) – Sequence of molecular configurations (one per reaction image), providing PySCF mol and overlap matrix for the Lowdin orthogonalization.

Return type:

ndarray[tuple[Any, …], dtype[float64]]

select(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → set[int]

Select atoms whose local electronic structure changes significantly across reaction images.

Parameters:

  • reaction (Reaction) – The reaction to analyze.

  • prepared_data (Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None) – Sequence of prepared data (one per reaction image) from the mean-field and localization stages. Both restricted and unrestricted data are supported.

  • molecular_configurations (Sequence[MolecularConfiguration] | None) – Sequence of molecular configurations (one per reaction image), providing PySCF mol and overlap matrix for the Lowdin orthogonalization.

Return type:

set[int]

threshold: float | None = None

Optional explicit threshold for active atom selection.

When None (default), the threshold is determined automatically using Jenks natural breaks clustering. When set, atoms with a change score at or above this value are considered active.

class EmbeddedAtomsConnectionExtender

Bases: EmbeddedAtomsExtender

Chemically “repairs” the active selection by connecting fragments into complete pieces.

Electronic-structure changes are often localized on a few atoms, but wavefunction methods generally perform better on chemically complete fragments. This extender therefore adds a minimum number of backbone atoms (typically carbons) needed to link the disconnected pieces of the active selection into a single connected fragment - or into as few fragments as possible.

Atoms are connected component by component: within each connected molecule of the system, the active atoms are linked using a minimum-weight Steiner tree of the connectivity graph, which adds the fewest intermediate atoms. Active atoms sitting in genuinely different molecules (for example a solute and a hydrogen-bonded water) cannot be connected and are left in separate fragments - the extender only links what is physically bonded.

After connecting the fragments, an AromaticRingExtender is applied so that whenever the connecting path enters an aromatic (or conjugated planar) ring, the whole ring is included rather than cutting through it.

__init__(aromatic_ring_extender: AromaticRingExtender = AromaticRingExtender(ring_sizes=(5, 6), planarity_tolerance=0.1)) → None
Parameters:

aromatic_ring_extender (AromaticRingExtender)

Return type:

None

aromatic_ring_extender: AromaticRingExtender = AromaticRingExtender(ring_sizes=(5, 6), planarity_tolerance=0.1)

Applied after connecting fragments to complete any aromatic ring that is partially included.

extend(active_indices: set[int], reaction: Reaction) → set[int]

Return the backbone atoms needed to connect the active fragments (plus completed rings).

Parameters:

  • active_indices (set[int]) – The current set of active atom indices.

  • reaction (Reaction) – The reaction being analyzed.

Return type:

set[int]

class EmbeddedAtomsSelector

Bases: object

Composes primary embedded atoms selectors and extenders to identify active atoms in a reaction.

Primary selectors identify the core active atoms (e.g., atoms with changed coordination or significant charge redistribution). Extenders then expand the selection (e.g., by including nearest neighbors, atoms within a distance threshold, or entire connected fragments).

Extenders are not chained: each extender is applied to the primary selection only, never to the atoms proposed by a preceding extender. The final selection is the union of the primary selection with every extender’s proposals. As a result the ordering of extenders does not affect the outcome, and an extender cannot cascade off atoms added by another extender (for example, a nearest-neighbor extender followed by a second nearest-neighbor extender does not grow the region by two bonds). To grow the embedded region past the immediate neighborhood of the primary atoms, use an extender that reaches the desired distance directly (e.g. DistanceExtender) rather than composing several short-range extenders.

__init__(primary_selectors: list[PrimaryEmbeddedAtomsSelector], extenders: list[EmbeddedAtomsExtender] | None = None, options: EmbeddedAtomsSelectorOptions | None = None) → None

Initialize the EmbeddedAtomsSelector.

Parameters:
Return type:

None

select_embedded_atoms(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → list[int]

Identify embedded atoms atoms by running all selectors and extenders.

The active atoms are formed by unioning each selector’s select() indices, then applying the extenders. Only if the result exceeds max_number_of_embedded_atoms are the selectors’ scores() consulted to rank the atoms and keep the highest-scoring ones.

Parameters:

  • reaction (Reaction) – The reaction to analyze.

  • prepared_data (Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None) – Optional sequence of prepared data (one per reaction image) from the mean-field and localization stages. Required by selectors that depend on electronic structure information.

  • molecular_configurations (Sequence[MolecularConfiguration] | None) – Optional sequence of molecular configurations (one per reaction image). Required by selectors that need basis set information.

Return type:

list[int]

static write_embedded_atoms_reaction_xyz(reaction: Reaction, embedded_atoms_indices: list[int], file_path: Path, *, comment: str = '') → None

Write a reaction XYZ file containing only the specified embedded_atoms atoms.

The output file contains multiple XYZ frames (one per reaction image), which can be viewed as an animation in molecular visualization tools.

Parameters:

  • reaction (Reaction) – The reaction to write.

  • embedded_atoms_indices (list[int]) – The atom indices to include in the output.

  • file_path (Path) – Path to the output XYZ file.

  • comment (str) – Optional comment to include in each frame’s comment line.

Return type:

None

class EmbeddedAtomsSelectorOptions

Bases: DataclassPublicAPI

Options controlling EmbeddedAtomsSelector.

All fields are immutable (frozen=True) so an instance can be safely reused.

Parameters:

max_number_of_embedded_atoms – Optional hard cap on the number of embedded atoms. When set and the selection (after extenders) exceeds it, only the highest-scoring atoms are kept. Primary atoms always outrank extender-added context atoms. Must be a positive integer when provided. (default=None)

__init__(*, max_number_of_embedded_atoms: int | None = None) → None
Parameters:

max_number_of_embedded_atoms (int | None)

Return type:

None

max_number_of_embedded_atoms: int | None = None
class FullMoleculeExtender

Bases: EmbeddedAtomsExtender

Extends the embedded atoms to include all atoms in connected fragments containing active atoms.

This is useful for ensuring that entire molecules are included when any part of them is active, which can be important for preserving chemical context in reactions involving many small molecules, like a reaction solvated in water where only a few water molecules are directly involved but the entire solvent environment is relevant.

__init__() → None
Return type:

None

extend(active_indices: set[int], reaction: Reaction) → set[int]

Return atom indices belonging to the same molecule as active atoms.

This assumes the reaction consist of several disconnected molecules (e.g., reactants, products, solvent molecules) and that we want to include entire molecules if any part of them is active.

Parameters:

  • active_indices (set[int]) – The current set of active atom indices.

  • reaction (Reaction) – The reaction being analyzed.

Return type:

set[int]

class HydrogenBondChangeSelector

Bases: PrimaryEmbeddedAtomsSelector

Selects atoms in hydrogen bonds whose geometry changes most across reaction images.

A hydrogen bond D-H...A is characterized by three atoms: a hydrogen atom, its covalently bonded heavy donor D and a nearby electronegative acceptor A (e.g. a C-H…O bond). A pair is treated as a hydrogen bond when the hydrogen-acceptor distance is at most hydrogen_acceptor_max_distance and the donor-hydrogen-acceptor angle is at least minimum_angle.

Each hydrogen bond is described by three geometric features - the donor-hydrogen distance, the hydrogen-acceptor distance and the donor-hydrogen-acceptor angle. The change of each feature (spread max - min across images) is made dimensionless via distance_change_scale and angle_change_scale and combined by their Euclidean norm. Each of the three participating atoms is scored by the largest change among its incident hydrogen bonds, and atoms are separated into active and inactive groups either by an unsupervised clustering algorithm (Jenks natural breaks) or a user-provided threshold.

This is a purely geometric selector detecting reactions driven by hydrogen-bond formation, breaking or rearrangement without a change in covalent connectivity.

__init__(acceptor_atomic_numbers: tuple[int, ...] = (7, 8, 9, 15, 16, 17), hydrogen_acceptor_max_distance: float = 2.5, minimum_angle: float = 120.0, distance_change_scale: float = 0.3, angle_change_scale: float = 30.0, threshold: float | None = None) → None
Parameters:

  • acceptor_atomic_numbers (tuple[int, ...])

  • hydrogen_acceptor_max_distance (float)

  • minimum_angle (float)

  • distance_change_scale (float)

  • angle_change_scale (float)

  • threshold (float | None)

Return type:

None

acceptor_atomic_numbers: tuple[int, ...] = (7, 8, 9, 15, 16, 17)

N, O, F, P, S, Cl).

Type:

Atomic numbers eligible to act as hydrogen-bond acceptors (default

angle_change_scale: float = 30.0

Angle change (in degrees) that maps to a unit contribution in the combined score.

distance_change_scale: float = 0.3

Distance change (in Angstrom) that maps to a unit contribution in the combined score.

hydrogen_acceptor_max_distance: float = 2.5

Maximum hydrogen-acceptor distance (in Angstrom) for a pair to be a hydrogen bond.

minimum_angle: float = 120.0

Minimum donor-hydrogen-acceptor angle (in degrees) for a pair to be a hydrogen bond.

scores(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → ndarray[tuple[Any, ...], dtype[float64]]

Score atoms by their largest incident hydrogen-bond change.

Parameters:
Return type:

ndarray[tuple[Any, …], dtype[float64]]

select(reaction: Reaction, prepared_data: Sequence[PreparedRestrictedData] | Sequence[PreparedUnrestrictedData] | None = None, molecular_configurations: Sequence[MolecularConfiguration] | None = None) → set[int]

Select atoms in the hydrogen bonds whose geometry changes most across reaction images.

Parameters:
Return type:

set[int]

threshold: float | None = None

Optional explicit threshold for active atom selection.

When None (default), the threshold is determined automatically using Jenks natural breaks clustering. When set, atoms whose largest incident hydrogen-bond change is at or above this value are considered active.

class NearestNeighborExtender

Bases: EmbeddedAtomsExtender

Extends the embedded atoms to include atoms directly bonded to any active atom in any image.

__init__() → None
Return type:

None

extend(active_indices: set[int], reaction: Reaction) → set[int]

Return atom indices that are nearest neighbors of active atoms.

Parameters:

  • active_indices (set[int]) – The current set of active atom indices.

  • reaction (Reaction) – The reaction being analyzed.

Return type:

set[int]