QM/MM — Region Partitioning, Embedding, Forces

Everything QM/MM in this library, in one place (issue #283): real region partitioning around a reactive bond, real coordinate slicing, a Diffuse2Seg-derived relevance propagation, and the real Hellmann-Feynman forces + Velocity-Verlet MD this repo already had (moved here from dashboard_core.qmmm, which still re-exports them for compatibility).

Region partitioning and propagation need the qmmm extra (pip install dense-evolution[qmmm], installs RDKit); forces/MD need no extra dependency beyond what dashboard_core already requires.

from rdkit import Chem
from dense_evolution.qmmm import partition_qm_mm_region, sliced_geometry

mol = Chem.AddHs(Chem.MolFromSmiles("OCCCCCC"))  # 1-hexanol
qm_atoms, boundary_pairs = partition_qm_mm_region(mol, {0, 1}, radius=1)
print(sorted(qm_atoms))       # [0, 1, 2, 7, 8, 9, 10, 11]
print(boundary_pairs)         # [(2, 3)] -- one bond crosses the boundary

partition_qm_mm_region walks outward from the reactive bond's atoms in heavy-atom hops. Hydrogens always follow their own heavy atom (never independently walked -- that spuriously cuts a terminal C-H bond and leaves an empty MM fragment). A boundary bond that would cut into an aromatic ring pulls the whole ring into the QM region first: cutting a lone aromatic atom out of its ring and capping it with hydrogen is not a valid molecule.

import numpy as np

atoms, geom_bohr = sliced_geometry(atomic_numbers, geom_bohr, qm_atoms, boundary_pairs)

sliced_geometry takes a coordinate SUBSET of one whole-molecule conformer -- never an independently re-embedded fragment, which gives unrelated 3D structures across fragments and was the real cause behind an MMFF94 mechanical correction that looked like it helped but turned out to be compensating for that geometry choice instead of truncation itself (see Dense-Evolution-Discovery's docs/qmmm_region_partitioning_mmff_correction.md for the full retraction).

from dense_evolution.qmmm import propagate_relevance

relevance = propagate_relevance(bond_order_matrix, seed_idx=[0, 1], n_nodes=n_atoms)
region = {i for i in range(n_atoms) if relevance[i] > threshold}

propagate_relevance is Algorithm 1 of Hümmer, Sicking, Hüger & Gottschalk 2026 ("Diffuse2Seg", arXiv:2609.06491) -- non-linear p-Laplacian graph propagation, solved by Gauss-Jacobi iteration, originally built to spread point prompts through a diffusion model's self-attention for image segmentation. Here "affinity" can be any node-graph weight, e.g. real Mayer/Wiberg bond order between atoms instead of self-attention between image patches. Measured, not assumed: on a real branched-aromatic molecule, this does NOT reproduce a naive "stronger bond propagates further" result at the paper's own calibrated hyperparameters (tuned for a dense multi-prompt image pipeline, not a single molecular seed) -- see Dense-Evolution-Discovery's docs/qmmm_utils.md for the full lambda sweep. Included here as a real, correctly-implemented primitive; not a proven QM/MM region-selection win yet.

from ase import Atoms
from dense_evolution.qmmm.ase_bridge import DenseEvolutionCalculator

h2 = Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.7414]])
h2.calc = DenseEvolutionCalculator(atomic_numbers=[1, 1], nuclear_charges=[1.0, 1.0],
                                    n_electrons=2, basis_name="sto-3g")
print(h2.get_potential_energy())  # -30.39 eV

DenseEvolutionCalculator (issue #288, needs the ase extra: pip install dense-evolution[ase]) is an ASE Calculator backed by native_hf's own differentiable energy (build_energy_fn) -- real Obara-Saika integrals and SCF, not a stub, for interop with ASE's optimizers/MD drivers and other engines' Atoms representations.

h2.calc = DenseEvolutionCalculator(atomic_numbers=[1, 1], nuclear_charges=[1.0, 1.0],
                                    n_electrons=2, basis_name="6-31g*")
print(h2.get_potential_energy())  # -30.66 eV

basis_name is a plain string -- "sto-3g", "6-31g", "6-31g*", anything basis_set_exchange has data for -- passed straight through to native_hf, which already supported arbitrary bases before this bridge existed (nothing here adds new basis-set capability; this is purely the ASE-interop layer). Swapping the basis on the SAME geometry is the whole point: 6-31G gives a lower (better, more variational freedom) energy than STO-3G for real hydrogen, exactly as physics requires. 6-31G and 6-31G* give the IDENTICAL energy for H2 specifically -- not a bug: * adds polarization d-functions to heavy atoms only, and hydrogen has none to add here.

Only energy is implemented (implemented_properties = ["energy"]) -- native_hf's own forces come from a separate, already-real implementation, compute_hellmann_feynman_forces above, with its own finite-difference derivative and its own calling convention (a molecule name from MOLECULE_CATALOG, not a bare ASE Atoms object). This bridge does not wrap that here, so ASE's gradient-based optimizers (BFGS, FIRE, ...) cannot be driven by it yet -- only single-point energies at any geometry/basis you construct directly.

qmmm

compute_hellmann_feynman_forces

compute_hellmann_feynman_forces(
    name: str,
    statevector=None,
    mapping: str = "jordan_wigner",
    geometry=None,
    fd_step_angstrom: float = 0.001,
)

Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of MOLECULE_CATALOG[name]: F = -d/dR, with H(R) this project's own real Hamiltonian (build_molecular_hamiltonian) and psi held fixed. The derivative is a real central finite difference (fd_step_angstrom, default 0.001 A -- verified converged against 0.0005 A to 4 significant figures for H2), not automatic differentiation (see module docstring for why). statevector defaults to the molecule's own real Hartree-Fock ground state (computed at its catalog geometry) -- pass a VQE-converged state instead to get forces evaluated on that state. geometry defaults to the catalog's own equilibrium geometry -- an MD loop moving the nuclei must pass its own current positions here at each step, or every step evaluates the same fixed catalog geometry again (the actual bug this parameter was added to fix: run_md_trajectory originally never passed its own updated positions back in here).

Cost (prog.txt, dashboard_core audit point 4a): the central-difference derivative evaluates energy_at 6n_atoms+1 times (H2's 2 atoms -> 13 Hamiltonian builds per call, each at a genuinely different geometry). This looks like it should be cacheable -- build_molecular_hamiltonian already caches by exact geometry -- but it isn't in practice: every one of the 13 geometries differs by fd_step_angstrom, so every call is a cache miss. Measured directly on H2 (dhf/PennyLane path) before deciding not to add a "cache the Pauli-term basis" layer here: HF + fermion-to-qubit mapping took 0.130s, Pauli-term extraction 0.0005s, dense matrix assembly 0.0031s -- the Hartree-Fock solve itself is 96%+ of the cost, not the bookkeeping after it, so caching the Pauli-term structure would save a few percent at best, not the 6n_atoms multiplier prog.txt's framing suggests. A real geometry change requires a real HF re-solve regardless of how its output gets packaged afterward -- see the module docstring's own account of why analytic differentiation (which WOULD avoid re-solving HF this many times) was tried and dropped for a real cross-platform PennyLane/ autograd bug, not reattempted here.

Parameters:

Name	Type	Description	Default
name	str	A key from MOLECULE_CATALOG (e.g. one of list(MOLECULE_CATALOG) -- these are descriptive strings like "H2 (Idrogeno) - R = 0.7414 A [equilibrio reale]", not bare element symbols like "H2").	required
statevector	see description above.		None
mapping	see description above.		None
geometry	see description above.		None
fd_step_angstrom	see description above.		None

Returns:

Type	Description
dict	name, symbols, energy_hartree, positions_angstrom, forces_hartree_per_angstrom, force_norm.

Examples:

>>> from dense_evolution.qmmm import compute_hellmann_feynman_forces
>>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
>>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
>>> result = compute_hellmann_feynman_forces(h2)
>>> round(result['energy_hartree'], 4)
-1.1373
>>> 0.01 < result['force_norm'] < 0.02  # small residual at equilibrium, not exactly zero
True

Source code in dense_evolution/qmmm/forces.py

def compute_hellmann_feynman_forces(name: str, statevector=None, mapping: str = "jordan_wigner",
                                     geometry=None, fd_step_angstrom: float = 0.001):
    """Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of
    MOLECULE_CATALOG[name]: F = -d<psi|H(R)|psi>/dR, with H(R) this
    project's own real Hamiltonian (build_molecular_hamiltonian) and psi
    held fixed. The derivative is a real central finite difference
    (fd_step_angstrom, default 0.001 A -- verified converged against
    0.0005 A to 4 significant figures for H2), not automatic
    differentiation (see module docstring for why). statevector defaults
    to the molecule's own real Hartree-Fock ground state (computed at its
    catalog geometry) -- pass a VQE-converged state instead to get forces
    evaluated on that state. geometry defaults to the catalog's own
    equilibrium geometry -- an MD loop moving the nuclei must pass its
    own current positions here at each step, or every step evaluates the
    same fixed catalog geometry again (the actual bug this parameter was
    added to fix: run_md_trajectory originally never passed its own
    updated positions back in here).

    Cost (prog.txt, dashboard_core audit point 4a): the central-difference
    derivative evaluates energy_at 6*n_atoms+1 times (H2's 2 atoms -> 13
    Hamiltonian builds per call, each at a genuinely different geometry).
    This looks like it should be cacheable -- build_molecular_hamiltonian
    already caches by exact geometry -- but it isn't in practice: every
    one of the 13 geometries differs by fd_step_angstrom, so every call
    is a cache miss. Measured directly on H2 (dhf/PennyLane path) before
    deciding not to add a "cache the Pauli-term basis" layer here: HF +
    fermion-to-qubit mapping took 0.130s, Pauli-term extraction 0.0005s,
    dense matrix assembly 0.0031s -- the Hartree-Fock solve itself is
    96%+ of the cost, not the bookkeeping after it, so caching the
    Pauli-term structure would save a few percent at best, not the
    6*n_atoms multiplier prog.txt's framing suggests. A real geometry
    change requires a real HF re-solve regardless of how its output gets
    packaged afterward -- see the module docstring's own account of why
    analytic differentiation (which WOULD avoid re-solving HF this many
    times) was tried and dropped for a real cross-platform PennyLane/
    autograd bug, not reattempted here.

    Parameters
    ----------
    name : str
        A key from `MOLECULE_CATALOG` (e.g. one of `list(MOLECULE_CATALOG)`
        -- these are descriptive strings like `"H2 (Idrogeno) - R = 0.7414
        A [equilibrio reale]"`, not bare element symbols like `"H2"`).
    statevector, mapping, geometry, fd_step_angstrom : see description above.

    Returns
    -------
    dict
        `name`, `symbols`, `energy_hartree`, `positions_angstrom`,
        `forces_hartree_per_angstrom`, `force_norm`.

    Examples
    --------
    >>> from dense_evolution.qmmm import compute_hellmann_feynman_forces
    >>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
    >>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
    >>> result = compute_hellmann_feynman_forces(h2)
    >>> round(result['energy_hartree'], 4)
    -1.1373
    >>> 0.01 < result['force_norm'] < 0.02  # small residual at equilibrium, not exactly zero
    True
    """
    MOLECULE_CATALOG, build_molecular_hamiltonian = _import_dashboard_hamiltonians()
    if name not in MOLECULE_CATALOG:
        raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
    spec = MOLECULE_CATALOG[name]
    symbols = spec["symbols"]
    if geometry is None:
        geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
    charge = spec["charge"]
    # BUG FIX (prog.txt, dashboard_core audit point 3a): these two were
    # never read from spec at all, so any catalog entry needing active-
    # space reduction (Si2 is the reason it's in MOLECULE_CATALOG) built
    # the FULL Hamiltonian instead of the reduced one here -- for Si2
    # specifically, 36 qubits instead of the intended 8, which
    # SafeMemoryGuard correctly refuses to allocate. Both
    # _reference_ground_state and energy_at below need these to build
    # the same, correctly-reduced Hamiltonian this molecule's other
    # dashboard panels (VQE, energy scan) already use.
    active_electrons = spec.get("active_electrons")
    active_orbitals = spec.get("active_orbitals")

    unknown = [s for s in symbols if s not in ATOMIC_MASSES_AMU]
    if unknown:
        raise ValueError(f"no real atomic mass on file for {unknown} -- add to "
                          f"ATOMIC_MASSES_AMU before using this molecule here")

    if statevector is None:
        statevector, _gs_energy, _n_qubits = _reference_ground_state(
            symbols, geometry, charge, mapping, active_electrons, active_orbitals)
    sv = np.asarray(statevector, dtype=np.complex128)
    geometry = np.asarray(geometry, dtype=np.float64)

    def energy_at(geom):
        h_matrix, _n_qubits = build_molecular_hamiltonian(
            symbols, geom, charge, mapping, active_electrons, active_orbitals)
        return float(np.real(np.vdot(sv, h_matrix @ sv)))

    energy = energy_at(geometry)
    forces = np.zeros_like(geometry)
    h = fd_step_angstrom
    for i in range(geometry.shape[0]):
        for j in range(3):
            geom_plus = geometry.copy()
            geom_plus[i, j] += h
            geom_minus = geometry.copy()
            geom_minus[i, j] -= h
            forces[i, j] = -(energy_at(geom_plus) - energy_at(geom_minus)) / (2 * h)

    return {
        "name": name,
        "symbols": symbols,
        "energy_hartree": energy,
        "positions_angstrom": geometry.tolist(),
        "forces_hartree_per_angstrom": forces.tolist(),
        "force_norm": float(np.linalg.norm(forces)),
    }

md_step

md_step(
    positions_angstrom,
    velocities_angstrom_per_fs,
    forces_hartree_per_angstrom,
    symbols,
    dt_fs: float = 0.5,
)

One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)dt/2, r(t+dt) = r(t) + v(t+dt/2)dt) using the real Hellmann-Feynman forces above and each atom's real atomic mass -- ordinary classical Newtonian mechanics (F=ma), nothing invented. Positions in Angstrom, velocities in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds.

Examples:

>>> import numpy as np
>>> from dense_evolution.qmmm import md_step
>>> positions = np.array([[0, 0, 0.0], [0, 0, 0.7414]])   # H2 at equilibrium
>>> velocities = np.zeros_like(positions)
>>> forces = np.array([[0, 0, 0.0109], [0, 0, -0.0109]])  # restoring force, pulling atoms together
>>> new_pos, new_vel, accel = md_step(positions, velocities, forces, ['H', 'H'], dt_fs=0.5)
>>> bool(new_pos[1, 2] < 0.7414)  # the second atom moved toward the first
True

Source code in dense_evolution/qmmm/forces.py

def md_step(positions_angstrom, velocities_angstrom_per_fs, forces_hartree_per_angstrom,
            symbols, dt_fs: float = 0.5):
    """One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)*dt/2,
    r(t+dt) = r(t) + v(t+dt/2)*dt) using the real Hellmann-Feynman forces
    above and each atom's real atomic mass -- ordinary classical Newtonian
    mechanics (F=ma), nothing invented. Positions in Angstrom, velocities
    in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds.

    Examples
    --------
    >>> import numpy as np
    >>> from dense_evolution.qmmm import md_step
    >>> positions = np.array([[0, 0, 0.0], [0, 0, 0.7414]])   # H2 at equilibrium
    >>> velocities = np.zeros_like(positions)
    >>> forces = np.array([[0, 0, 0.0109], [0, 0, -0.0109]])  # restoring force, pulling atoms together
    >>> new_pos, new_vel, accel = md_step(positions, velocities, forces, ['H', 'H'], dt_fs=0.5)
    >>> bool(new_pos[1, 2] < 0.7414)  # the second atom moved toward the first
    True
    """
    positions = np.asarray(positions_angstrom, dtype=np.float64)
    velocities = np.asarray(velocities_angstrom_per_fs, dtype=np.float64)
    forces = np.asarray(forces_hartree_per_angstrom, dtype=np.float64)
    masses = np.array([ATOMIC_MASSES_AMU[s] for s in symbols], dtype=np.float64)

    accel = ACCEL_CONVERSION * forces / masses[:, None]
    velocities_half = velocities + 0.5 * accel * dt_fs
    positions_new = positions + velocities_half * dt_fs
    return positions_new, velocities_half, accel

run_md_trajectory

run_md_trajectory(
    name: str,
    n_steps: int,
    dt_fs: float = 0.5,
    mapping: str = "jordan_wigner",
    recompute_electronic_state: bool = False,
    fd_step_angstrom: float = 0.001,
)

Real, minimal ab-initio-forces MD trajectory: at each step, real Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the real nuclear positions/velocities via real Velocity-Verlet (md_step). Starts from rest (zero initial velocities) at the catalog's real equilibrium geometry.

fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at every step -- previously not exposed here at all, silently using that function's own default (0.001 A) with no way for a caller to ask for a different finite-difference step (e.g. a molecule with an unusually steep energy landscape, where the default step size isn't the one already verified converged for H2).

recompute_electronic_state=False (default) holds the electronic state fixed at the initial Hartree-Fock reference through the whole trajectory -- forces stay exact only close to the starting geometry (a real, explicitly-stated approximation, not a fabricated one). True ab-initio MD (re-solving Hartree-Fock at every step's new geometry) is available by setting this True, at real, substantial extra cost per step.

Examples:

>>> from dense_evolution.qmmm import run_md_trajectory
>>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
>>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
>>> traj = run_md_trajectory(h2, n_steps=3, dt_fs=0.5)
>>> traj['step']
[0, 1, 2]
>>> len(traj['force_norm'])
3

Source code in dense_evolution/qmmm/forces.py

def run_md_trajectory(name: str, n_steps: int, dt_fs: float = 0.5, mapping: str = "jordan_wigner",
                       recompute_electronic_state: bool = False, fd_step_angstrom: float = 0.001):
    """Real, minimal ab-initio-forces MD trajectory: at each step, real
    Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the
    real nuclear positions/velocities via real Velocity-Verlet (md_step).
    Starts from rest (zero initial velocities) at the catalog's real
    equilibrium geometry.

    fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at
    every step -- previously not exposed here at all, silently using
    that function's own default (0.001 A) with no way for a caller to
    ask for a different finite-difference step (e.g. a molecule with an
    unusually steep energy landscape, where the default step size isn't
    the one already verified converged for H2).

    recompute_electronic_state=False (default) holds the electronic state
    fixed at the initial Hartree-Fock reference through the whole
    trajectory -- forces stay exact only close to the starting geometry
    (a real, explicitly-stated approximation, not a fabricated one).
    True ab-initio MD (re-solving Hartree-Fock at every step's new
    geometry) is available by setting this True, at real, substantial
    extra cost per step.

    Examples
    --------
    >>> from dense_evolution.qmmm import run_md_trajectory
    >>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
    >>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
    >>> traj = run_md_trajectory(h2, n_steps=3, dt_fs=0.5)
    >>> traj['step']
    [0, 1, 2]
    >>> len(traj['force_norm'])
    3
    """
    MOLECULE_CATALOG, _build_molecular_hamiltonian = _import_dashboard_hamiltonians()
    if name not in MOLECULE_CATALOG:
        raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
    spec = MOLECULE_CATALOG[name]
    symbols = spec["symbols"]
    geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
    charge = spec["charge"]
    active_electrons = spec.get("active_electrons")
    active_orbitals = spec.get("active_orbitals")

    statevector, _gs_energy, _n_qubits = _reference_ground_state(
        symbols, geometry, charge, mapping, active_electrons, active_orbitals)
    positions = np.asarray(geometry, dtype=np.float64)
    velocities = np.zeros_like(positions)

    trajectory = {"step": [], "time_fs": [], "positions_angstrom": [], "energy_hartree": [], "force_norm": []}
    for step in range(n_steps):
        result = compute_hellmann_feynman_forces(name, statevector, mapping=mapping, geometry=positions,
                                                  fd_step_angstrom=fd_step_angstrom)
        forces = np.asarray(result["forces_hartree_per_angstrom"])
        trajectory["step"].append(step)
        trajectory["time_fs"].append(step * dt_fs)
        trajectory["positions_angstrom"].append(positions.tolist())
        trajectory["energy_hartree"].append(result["energy_hartree"])
        trajectory["force_norm"].append(result["force_norm"])

        positions, velocities, _accel = md_step(positions, velocities, forces, symbols, dt_fs=dt_fs)
        _assert_no_nuclear_collision(positions, step, dt_fs)

        if recompute_electronic_state:
            statevector, _gs_energy, _n_qubits = _reference_ground_state(
                symbols, positions, charge, mapping, active_electrons, active_orbitals)

    return trajectory

partition_qm_mm_region

partition_qm_mm_region(mol, seed_heavy_atoms, radius)

BFS outward from seed_heavy_atoms (a set/list of atom indices) across radius heavy-atom hops. Hydrogens always follow their own heavy atom afterward (never independently expanded -- see module docstring). Any boundary bond that would cut into an aromatic ring pulls that whole ring into the QM region first, iterating until stable, before hydrogens are added and the boundary is finalized.

Returns (qm_atoms: set[int], boundary_pairs: list[(kept_idx, cut_idx)]), one pair per bond crossing the QM/MM boundary, kept_idx on the QM side.

Examples:

>>> from rdkit import Chem
>>> from dense_evolution.qmmm import partition_qm_mm_region
>>> mol = Chem.AddHs(Chem.MolFromSmiles("OCCCCCC"))
>>> qm_atoms, boundary_pairs = partition_qm_mm_region(mol, {0, 1}, radius=1)
>>> len(boundary_pairs)
1

Source code in dense_evolution/qmmm/region.py

def partition_qm_mm_region(mol, seed_heavy_atoms, radius):
    """BFS outward from `seed_heavy_atoms` (a set/list of atom indices)
    across `radius` heavy-atom hops. Hydrogens always follow their own
    heavy atom afterward (never independently expanded -- see module
    docstring). Any boundary bond that would cut into an aromatic ring
    pulls that whole ring into the QM region first, iterating until
    stable, before hydrogens are added and the boundary is finalized.

    Returns (qm_atoms: set[int], boundary_pairs: list[(kept_idx, cut_idx)]),
    one pair per bond crossing the QM/MM boundary, `kept_idx` on the QM
    side.

    Examples
    --------
    >>> from rdkit import Chem
    >>> from dense_evolution.qmmm import partition_qm_mm_region
    >>> mol = Chem.AddHs(Chem.MolFromSmiles("OCCCCCC"))
    >>> qm_atoms, boundary_pairs = partition_qm_mm_region(mol, {0, 1}, radius=1)
    >>> len(boundary_pairs)
    1
    """
    qm_heavy = set(seed_heavy_atoms)
    frontier = set(seed_heavy_atoms)
    for _ in range(radius):
        new_frontier = set()
        for idx in frontier:
            for nbr in mol.GetAtomWithIdx(idx).GetNeighbors():
                if nbr.GetAtomicNum() > 1 and nbr.GetIdx() not in qm_heavy:
                    new_frontier.add(nbr.GetIdx())
        qm_heavy |= new_frontier
        frontier = new_frontier

    ring_info = mol.GetRingInfo()
    changed = True
    while changed:
        changed = False
        for bond in mol.GetBonds():
            if not bond.GetIsAromatic():
                continue
            a, b = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
            if (a in qm_heavy) != (b in qm_heavy):
                for ring in ring_info.AtomRings():
                    if a in ring or b in ring:
                        newly = set(ring) - qm_heavy
                        if newly:
                            qm_heavy |= newly
                            changed = True

    qm_atoms = set(qm_heavy)
    for idx in qm_heavy:
        for nbr in mol.GetAtomWithIdx(idx).GetNeighbors():
            if nbr.GetAtomicNum() == 1:
                qm_atoms.add(nbr.GetIdx())

    raw_boundary = [(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in mol.GetBonds()
                    if (b.GetBeginAtomIdx() in qm_atoms) != (b.GetEndAtomIdx() in qm_atoms)]
    boundary_pairs = [(i, j) if i in qm_atoms else (j, i) for i, j in raw_boundary]
    return qm_atoms, boundary_pairs

sliced_geometry

sliced_geometry(
    atomic_numbers, geom_bohr, keep_idx, boundary_pairs
)

A coordinate SUBSET of one whole-molecule conformer (never an independently re-embedded fragment -- that would give unrelated 3D structures across fragments, and was the real cause of the MMFF94 correction's apparent benefit turning out to be a geometry artifact, see module docstring). Boundary bonds get a capping H placed along the kept->cut bond direction at a standard C-H bond length.

Examples:

>>> import numpy as np
>>> from dense_evolution.qmmm import sliced_geometry
>>> numbers = [8, 6, 6]
>>> geom = np.array([[0.0, 0.0, 0.0], [1.4, 0.0, 0.0], [2.8, 0.0, 0.0]])
>>> new_numbers, new_geom = sliced_geometry(numbers, geom, {0, 1}, [(1, 2)])
>>> new_numbers
[8, 6, 1]

Source code in dense_evolution/qmmm/region.py

def sliced_geometry(atomic_numbers, geom_bohr, keep_idx, boundary_pairs):
    """A coordinate SUBSET of one whole-molecule conformer (never an
    independently re-embedded fragment -- that would give unrelated 3D
    structures across fragments, and was the real cause of the MMFF94
    correction's apparent benefit turning out to be a geometry artifact,
    see module docstring). Boundary bonds get a capping H placed along
    the kept->cut bond direction at a standard C-H bond length.

    Examples
    --------
    >>> import numpy as np
    >>> from dense_evolution.qmmm import sliced_geometry
    >>> numbers = [8, 6, 6]
    >>> geom = np.array([[0.0, 0.0, 0.0], [1.4, 0.0, 0.0], [2.8, 0.0, 0.0]])
    >>> new_numbers, new_geom = sliced_geometry(numbers, geom, {0, 1}, [(1, 2)])
    >>> new_numbers
    [8, 6, 1]
    """
    keep_idx = sorted(keep_idx)
    new_numbers = [atomic_numbers[i] for i in keep_idx]
    new_geom = [geom_bohr[i] for i in keep_idx]
    for kept, cut in boundary_pairs:
        vec = geom_bohr[cut] - geom_bohr[kept]
        vec = vec / np.linalg.norm(vec)
        new_geom.append(geom_bohr[kept] + vec * CH_BOND_BOHR)
        new_numbers.append(1)
    return new_numbers, np.array(new_geom)

propagate_relevance

propagate_relevance(
    affinity,
    seed_idx,
    n_nodes,
    p=1.6,
    lam=1e-05,
    tau_prop=0.0001,
    max_iter=500,
)

Algorithm 1 of Hummer, Sicking, Huger & Gottschalk 2026 (arXiv:2609.06491, "Diffuse2Seg", Sec. 3.4/A.1), read directly from the paper and implemented verbatim.

Non-linear p-Laplacian graph-regularized smoothing (Elmoataz et al. 2008), solved by Gauss-Jacobi iteration: propagates a one-hot seed vector over any node-affinity graph affinity (self-attention in the original paper; Mayer/Wiberg bond order in Dense-Evolution-Discovery's QM/MM experiments) into a soft relevance map that stays smooth within high-affinity regions and is throttled across low-affinity (edge) ones. p, lam, tau_prop are the paper's own final values (Sec. 4.2).

Measured, not assumed (Dense-Evolution-Discovery's qmmm_diffuse2seg_propagation_lambda_sweep.py, on the real Mayer bond-order graph of a branched-aromatic molecule, OCC(c1ccccc1)CCC): at the paper's own lam=1e-5 -- tuned for a dense grid of prompts later merged together, not a single isolated seed -- a stronger real bond (aromatic ring, bond order 1.412) actually propagates LESS relevance than a weaker one (alkyl chain, bond order 0.991), the opposite of the naive expectation. The ratio crosses 1.0 only around lam~0.5-1, and only reaches a large expected-direction differentiation (ratio 1.49) at lam=10, two orders of magnitude above the paper's own calibrated value. This is reported as the real, measured lam-dependence of this algorithm in a single-seed molecular setting, distinct from the many-prompt image setting it was designed and calibrated for -- lam is not silently retuned to whatever value looks best.

Degenerate-case handling: at any node i where g_i = sqrt(sum_j A_ij(f_j-f_i)^2) is exactly 0 (every affinity-neighbor already equals f_i), the formula's g_i^(p-2) term diverges for p<2. An earlier version clamped g_i away from 0 with an epsilon -- the same category of shortcut already rejected elsewhere in this project for degenerate eigenvalues (see dense_evolution.physics.spectral, Kato's divided- difference formula) in favor of the real mathematical limit. Worked out directly here: as g_i -> 0, every A_ij-connected f_j equals f_i by definition of g_i=0, so the g_i^(p-2)-weighted terms in both the numerator and denominator of the update come to dominate and cancel to exactly f_i -- i.e. a node already consistent with its whole affinity-neighborhood is unchanged by an edge-preserving smoothing step, exactly as expected. Implemented as an explicit special case below, not an epsilon.

Examples:

>>> import numpy as np
>>> from dense_evolution.qmmm import propagate_relevance
>>> affinity = np.array([[0, 1.0, 0], [1.0, 0, 1.0], [0, 1.0, 0]])
>>> rel = propagate_relevance(affinity, [0], 3, lam=1.0)
>>> bool(rel[0] > rel[1] > rel[2])
True

Source code in dense_evolution/qmmm/propagation.py

def propagate_relevance(affinity, seed_idx, n_nodes, p=1.6, lam=1e-5, tau_prop=1e-4, max_iter=500):
    """Algorithm 1 of Hummer, Sicking, Huger & Gottschalk 2026
    (arXiv:2609.06491, "Diffuse2Seg", Sec. 3.4/A.1), read directly from
    the paper and implemented verbatim.

    Non-linear p-Laplacian graph-regularized smoothing (Elmoataz et al.
    2008), solved by Gauss-Jacobi iteration: propagates a one-hot seed
    vector over any node-affinity graph `affinity` (self-attention in the
    original paper; Mayer/Wiberg bond order in Dense-Evolution-Discovery's
    QM/MM experiments) into a soft relevance map that stays smooth within
    high-affinity regions and is throttled across low-affinity (edge)
    ones. `p`, `lam`, `tau_prop` are the paper's own final values
    (Sec. 4.2).

    Measured, not assumed (Dense-Evolution-Discovery's
    qmmm_diffuse2seg_propagation_lambda_sweep.py, on the real Mayer
    bond-order graph of a branched-aromatic molecule,
    OCC(c1ccccc1)CCC): at the paper's own lam=1e-5 -- tuned for a dense
    grid of prompts later merged together, not a single isolated seed --
    a stronger real bond (aromatic ring, bond order 1.412) actually
    propagates LESS relevance than a weaker one (alkyl chain, bond order
    0.991), the opposite of the naive expectation. The ratio crosses 1.0
    only around lam~0.5-1, and only reaches a large expected-direction
    differentiation (ratio 1.49) at lam=10, two orders of magnitude above
    the paper's own calibrated value. This is reported as the real,
    measured lam-dependence of this algorithm in a single-seed molecular
    setting, distinct from the many-prompt image setting it was designed
    and calibrated for -- `lam` is not silently retuned to whatever value
    looks best.

    Degenerate-case handling: at any node i where g_i = sqrt(sum_j
    A_ij(f_j-f_i)^2) is exactly 0 (every affinity-neighbor already equals
    f_i), the formula's g_i^(p-2) term diverges for p<2. An earlier
    version clamped g_i away from 0 with an epsilon -- the same category
    of shortcut already rejected elsewhere in this project for degenerate
    eigenvalues (see dense_evolution.physics.spectral, Kato's divided-
    difference formula) in favor of the real mathematical limit. Worked
    out directly here: as g_i -> 0, every A_ij-connected f_j equals f_i
    by definition of g_i=0, so the g_i^(p-2)-weighted terms in both the
    numerator and denominator of the update come to dominate and cancel
    to exactly f_i -- i.e. a node already consistent with its whole
    affinity-neighborhood is unchanged by an edge-preserving smoothing
    step, exactly as expected. Implemented as an explicit special case
    below, not an epsilon.

    Examples
    --------
    >>> import numpy as np
    >>> from dense_evolution.qmmm import propagate_relevance
    >>> affinity = np.array([[0, 1.0, 0], [1.0, 0, 1.0], [0, 1.0, 0]])
    >>> rel = propagate_relevance(affinity, [0], 3, lam=1.0)
    >>> bool(rel[0] > rel[1] > rel[2])
    True
    """
    f0 = np.zeros(n_nodes)
    f0[seed_idx] = 1.0
    f = f0.copy()
    for _ in range(max_iter):
        diff = f[None, :] - f[:, None]
        g = np.sqrt(np.sum(affinity * diff ** 2, axis=1))
        degenerate = g == 0.0
        gp = np.zeros_like(g)
        gp[~degenerate] = g[~degenerate] ** (p - 2)
        gamma = affinity * (gp[:, None] + gp[None, :])
        numerator = lam * f0 + (gamma * f[None, :]).sum(axis=1)
        denominator = lam + gamma.sum(axis=1)
        f_new = np.where(degenerate, f, numerator / denominator)
        if np.sum((f_new - f) ** 2) <= tau_prop:
            f = f_new
            break
        f = f_new
    return f

---

See also: dense_evolution.native_hf for the Hartree-Fock engine these region functions feed into; dashboard_core.hamiltonians for MOLECULE_CATALOG and build_molecular_hamiltonian, which compute_hellmann_feynman_forces/run_md_trajectory are built from.