Select the Embedded Region with the EmbeddedAtomsSelector
When you build a projective-embedding ground-state problem or an even-handed reaction-path problem, you must decide which atoms belong to the embedded region — the chemically important part treated with a high-level wavefunction method — and which atoms form the environment, treated at the mean-field (DFT) level. See Projection-Based Wavefunction-in-DFT Embedding for the theory behind this partition.
Goal
Automatically select the atom indices to include in the embedded region of a Projection-Based Wavefunction-in-DFT (WF-in-DFT) calculation, instead of specifying them by hand.
Specifying the Embedded Region Manually
Before reaching for automation, remember that you can always specify the embedded region explicitly using the
embedded_atoms=[...] keyword (indices of atoms to embed).
For a single molecular configuration:
import qrunch as qc
molecular_configuration = qc.build_molecular_configuration(
[
("O", 0.0, 0.0, 0.11779),
("H", 0.0, 0.75545, -0.47116),
("H", 0.0, -0.75545, -0.47116),
("O", water_separation, 0.0, 0.11779),
("H", water_separation, 0.75545, -0.47116),
("H", water_separation, -0.75545, -0.47116),
],
basis_set="sto3g",
embedded_atoms=[0, 1, 2], # indices of atoms in the embedded region
)
or for a reaction configuration:
import qrunch as qc
reaction_configuration = qc.build_reaction_configuration(
reaction=[
[("H", 0.0, 0.0, 0.0), ("H", 0.72, 0.0, 0.0)],
[("H", 0.0, 0.0, 0.0), ("H", 0.74, 0.0, 0.0)],
[("H", 0.0, 0.0, 0.0), ("H", 0.76, 0.0, 0.0)],
],
basis_set="sto3g",
embedded_atoms=[0],
)
See Create a MolecularConfiguration and Create a ReactionConfiguration for more on manual setup.
Automating the Choice: the EmbeddedAtomsSelector
Choosing embedded_atoms by hand is fine for small systems, but it becomes tedious and error-prone for large
molecules or long reaction paths. The EmbeddedAtomsSelector is Qrunch’s effort to automate this choice: given a
reaction (a sequence of molecular geometries, or images), it selects the atom indices whose local chemistry changes
most along the path, and returns them as the embedded region.
The EmbeddedAtomsSelector combines two kinds of components:
Selectors identify the core active atoms — those that genuinely participate in the reaction.
Extenders grow that core selection with additional context atoms (neighbors, whole molecules, aromatic rings, …), so that the embedded region is chemically complete.
Building an EmbeddedAtomsSelector
The EmbeddedAtomsSelector is created with the qc.embedded_atoms_selector_creator() builder, which follows
Qrunch’s fluent builder pattern: each .with_… method configures one
selector, extender, or option and returns the builder, so calls can be chained and terminated with .create().
import qrunch as qc
embedded_atoms_selector = (
qc.embedded_atoms_selector_creator()
.with_coordination_change_selector()
.with_eeqbc_change_selector()
.with_electronic_structure_change_selector()
.with_bond_length_change_selector()
.with_hydrogen_bond_change_selector()
.with_options(options=qc.options.EmbeddedAtomsSelectorOptions(max_number_of_embedded_atoms=8))
.with_distance_extender(distance_threshold=2.0)
.with_nearest_neighbor_extender()
.with_full_molecule_extender()
.with_aromatic_ring_extender()
.with_embedded_atoms_connection_extender()
.create()
)
Each line is explained below.
Selectors: Identifying the Core Active Atoms
Selectors examine the reaction images and flag the atoms whose chemistry changes. You can add any combination of them; their selections are combined (unioned).
.with_coordination_change_selector()
Include a coordination change selector. This selector identifies atoms whose coordination environment changes across reaction images. It builds a connectivity graph for each image and detects atoms whose sorted neighbor atomic numbers change between any pair of images.
The coordination change selector highlighting the atoms whose bonding changes during an SN2 reaction.
.with_eeqbc_change_selector()
Include an EEQBC atomic charge change selector. This selector identifies atoms whose partial charges change significantly across reaction images, using the bond capacity electronegativity equilibration charge model (EEQBC) [FMullerHG25]. It is useful for detecting charge-transfer reactions where the bond connectivity may not change. Compared to ordinary QEq/EEQ, EEQBC introduces a distance-dependent bond capacitance that strongly suppresses artificial long-range charge transfer between distant fragments, so it is less likely to select chemically irrelevant, distant spectator atoms.
.with_electronic_structure_change_selector()
Include an electronic structure change selector. This selector identifies atoms whose local electronic structure (measured via IAO density matrices) changes significantly across reaction images. It supports both restricted and unrestricted calculations.
Note
Unlike the other selectors, with_electronic_structure_change_selector requires access to molecular orbital
information (see Using the EmbeddedAtomsSelector on Its Own below).
The EEQBC and electronic structure change selectors responding to charge redistribution in nitroaniline.
.with_bond_length_change_selector()
Include a bond-length change selector. This selector identifies atoms involved in bonds whose length changes most across reaction images. It is a purely geometric, rigid-motion-invariant selector that also detects symmetric homonuclear bond stretches (e.g. H2 dissociation) that electronic descriptors cannot see.
The bond-length change selector tracking the bonds that stretch and contract in glyoxal.
.with_hydrogen_bond_change_selector()
Include a hydrogen-bond change selector. This selector identifies atoms in hydrogen bonds (donor–hydrogen···acceptor, e.g. C–H···O) whose geometry — the donor–hydrogen distance, the hydrogen–acceptor distance and the donor–hydrogen–acceptor angle — changes most across reaction images. It is a purely geometric, rigid-motion-invariant selector that detects reactions driven by hydrogen-bond formation, breaking or rearrangement.
The hydrogen-bond change selector following the hydrogen bonds in a formic acid dimer.
Options: Controlling the Selection
.with_options()
Set additional options on the selector via an
EmbeddedAtomsSelectorOptions instance. The most important
option is max_number_of_embedded_atoms, which limits the size of the embedded region: if the selectors and extenders
together propose more atoms than this cap, only the highest-scoring atoms are kept (primary atoms always outrank
extender-added context atoms).
from qrunch.options import EmbeddedAtomsSelectorOptions
.with_options(options=EmbeddedAtomsSelectorOptions(max_number_of_embedded_atoms=8))
Extenders: Growing the Selection
Extenders start from the primary atoms selected by the selectors and add further atoms in various ways, so that the embedded region contains enough chemical context to be treated accurately.
.with_nearest_neighbor_extender()
Include a nearest neighbor extender. This extender includes atoms directly bonded to any active atom in any image.
.with_distance_extender(distance_threshold)
Include a distance-based extender. This extender includes atoms within distance_threshold (in Angstrom) from any
active atom.
.with_full_molecule_extender()
Include a full molecule extender. This extender includes all atoms in connected fragments that contain active atoms. It is useful for ensuring entire molecules are included when any part of them is active — for example a reaction solvated in water, where you want to keep each involved solvent molecule whole.
.with_aromatic_ring_extender()
Include an aromatic ring extender. This extender includes the complete aromatic ring (e.g. a benzene ring) whenever any of its atoms is active, so that the whole conjugated ring is embedded together. Only the ring atoms are added — hydrogens and side chains attached to the ring are not included. Aromatic rings are detected geometrically as planar rings of aromatic-eligible elements, which handles both benzene and heteroaromatic rings (e.g. pyridine, furan).
.with_embedded_atoms_connection_extender()
Include a connection extender that chemically repairs the active fragment. This extender adds the minimum number of backbone atoms needed to link the disconnected pieces of the active selection into a single connected fragment (or as few fragments as possible — atoms in genuinely different molecules cannot be connected). Afterward it applies an aromatic ring extender so that connecting paths do not cut through aromatic (or conjugated planar) ring systems.
Using the EmbeddedAtomsSelector in a Problem Builder
The EmbeddedAtomsSelector can be attached to the even-handed reaction-path problem builder, which will use it to determine the embedded region consistently along the reaction path (see Construct a Reaction-Path Problem):
import qrunch as qc
even_handed_reaction_problem_builder = (
qc.problem_builder_creator()
.reaction_path()
.even_handed()
.with_embedded_atoms_selector(embedded_atoms_selector)
.create()
)
Using the EmbeddedAtomsSelector on Its Own
The EmbeddedAtomsSelector can also be used directly to obtain the embedded atom indices for a reaction:
embedded_atoms = embedded_atoms_selector.select_embedded_atoms(
reaction=reaction_configuration.reaction
)
Warning
The with_electronic_structure_change_selector selector will perform a small Hartree-Fock calculation
on each molecular configuration, because it requires access to molecular orbital information.
This may be expensive and the purely geometric and charge-based selectors
(coordination, EEQBC, bond-length and hydrogen-bond change), may be preferable for large systems.
Visualizing the Selection
You can write the selected atoms to a multi-frame XYZ file for inspection using
write_embedded_atoms_reaction_xyz(). The
output file contains one XYZ frame per reaction image, which can be viewed as an animation in molecular visualization
tools:
from pathlib import Path
EmbeddedAtomsSelector.write_embedded_atoms_reaction_xyz(
reaction=reaction_configuration.reaction,
embedded_atoms_indices=embedded_atoms,
file_path=Path("embedded_atoms.xyz"),
)