Getting Started#

seamm_mbe is the bookkeeping behind SEAMM’s many-body correction (MBE). It estimates high-level energies, forces and stress for a periodic cell or a large cluster as a cheap calculation of the whole system plus many-body increments of [high − low], each computed on a small isolated fragment:

E(cell) ≈ E_low(cell) + Σ_F dE_F

The increments come from monomers, selected pairs and selected triples (designed for general order). The library does no calculations itself. The mbe_step plug-in runs them through the Model Chemistry batch contract and passes the results here.

Installing#

pip install seamm_mbe

It needs numpy, and molsystem for System.from_configuration().

Units#

Everything is in eV, Å, eV/Å (forces, not gradients) and eV for virials (W = Σ r ⊗ f, so P = W/V). Stress σ is in eV/Å³ with σ = −P, the ASE/xnn convention. seamm_mbe.units converts from SEAMM’s analyze_task results (kJ/mol, kJ/mol/Å gradients, GPa) and back. Every external stress needs an explicit convention: "pressure" (VASP, MDI) or "stress" (ASE, xnn).

Basic usage#

import seamm_mbe

# The system: from arrays (bonds define the molecules) ...
system = seamm_mbe.System(symbols, xyz, cell, bonds=bonds)
# ... or from a molsystem configuration with bonds
system = seamm_mbe.System.from_configuration(configuration)

# The fragments: the prototype's water rules are the default
rules = seamm_mbe.SelectionRules(
    max_order=3,
    criterion="designated",            # water O-O
    cutoffs={2: 4.5, 3: 3.5},          # or tables by type pair
    rules={3: "connected"},
)
fragments = seamm_mbe.enumerate_fragments(system, rules)

# The mixed low level: monomers and pairs < 3.5 Å periodic, the rest molecular
seamm_mbe.assign_levels(fragments, {1: True, 2: 3.5})
todo = fragments.calculations()   # {"high": [...], "periodic": [...], "molecular": [...]}

# ... run them (mbe_step's job), giving {name: (energy_eV, forces_eV_per_A)} ...

correction = seamm_mbe.mbe_correction(fragments, high, periodic, molecular)
labels = seamm_mbe.assemble(
    system,
    correction,
    [
        seamm_mbe.CellTerm("VASP r2SCAN", e, f, seamm_mbe.units.vasp_stress_to_virial(kB, V)),
        seamm_mbe.CellTerm("D4", e4, f4, seamm_mbe.units.dftd4_virial_to_virial(dEde)),
    ],
    offsets={"water": 2074.69325},     # eV per molecule, to the training sets' scale
)
labels.reference_energy, labels.forces, labels.stress, labels.molecular_pressure

labels.per_body breaks the correction down by body order, and labels.breakdown by cell term and “MBE”. Each entry gives the energy, its share of the atomic pressure and its share of the molecular pressure. The molecular share is the term’s virial less its own intramolecular part, so the shares add up to labels.molecular_pressure. The monomers’ molecular share is zero, since a monomer increment has no net force on its molecule. This is the number to compare with the virial errors of a cluster validation.

Molecules and types#

The molecules are the connected parts of the bond graph, made whole across the cell boundary. Each is typed by its formula and a Weisfeiler–Lehman hash of its bond graph. The built-in catalog has water, EC, FEC, DMC, EMC, Li⁺, BF₄⁻, PF₆⁻ and the monatomic ions Na⁺, K⁺, F⁻, Cl⁻, Br⁻ and I⁻, each with its charge and designated atom (water’s O, the carbonyl C, the ion’s central atom). Other molecules are named by their formula.

Charges come from the structure’s per-atom formal charges when it has any (e.g. SDF M CHG). They are checked against the catalog, and a mismatch names the molecule. Otherwise the catalog’s charges apply. An electron-count parity check catches a charge put on the wrong molecule. Explicit charges by type override both.

Fragment names#

Names are the prototype’s (gen_stage2_frame.py), so stored results match:

m00: a monomer.
d03_41: a pair. When the second molecule is not at its minimum image, the image follows, e.g. d03_41_x0p0.
t00_17_30 and q..: selected triples and 4-bodies.

A selected triple’s or 4-body’s sub-fragments that are not selected themselves (e.g. the “outer pair” of a connected triple) are computed too, with in_sum False. Higher-order ones carry their images in their names.

In a cell narrower than the cutoffs need for uniqueness, the same molecules can form several different fragments through different images: in a 15.2 Å cell of 32 ethylene carbonates, molecules 0, 1 and 25 form three different connected triples. SelectionRules.check() then warns, and each of these fragments is enumerated and computed on its own. The one with every molecule at its minimum image keeps the plain name (t00_01_25); the others carry their images (t00_01_25_xm00m00), as do their sub-pairs. A fragment never holds a molecule twice: when molecule a is bonded to two images of b, the triple (b, a, b’) is skipped, while (a, b, k) and (a, b’, k) are both enumerated. Sub-fragments that equal another fragment up to whole cells share its key and are computed once.

The neighbour search looks as many image layers out as the cutoffs need, so a cell narrower than the cutoff plus the molecules’ size (as with large units such as ion shells) is handled. One case is refused, since no fragment list would be right: a molecule within the pair cutoff of its own image, a real interaction that no fragment can hold, which skipping would silently drop.

Corrections to selected increments#

mbe_correction(..., corrections={name: (dE, dF)}) adds a correction to a selected fragment’s increment in the sum only, never to the sub-fragment increments that higher fragments subtract. Pairwise counterpoise works this way: with dE_ij^CP − dE_ij for each pair, the result is E(1) + Σ dE_ij^CP + Σ dE_ijk. The triples still subtract the uncorrected pairs, so a pair’s basis-set superposition error does not move into the 3-body terms.

A high level per order#

The high level can differ by order, e.g. revDSD at def2-QZVPPD for the monomers and pairs and at def2-TZVPPD for the triples. Each order’s increments are then built entirely at that order’s level: a triple’s increment uses the triple, its pairs and its monomers all at the triples’ level, so those pairs and monomers are computed at both levels:

todo = fragments.calculations(high_levels={3: "high:3"})
# todo["high"]: the monomers and pairs (and their sub-fragments)
# todo["high:3"]: the triples with all their sub-pairs and monomers
correction = seamm_mbe.mbe_correction(
    fragments, high, periodic, molecular, high_by_order={3: high_tz}
)

Each increment records the high level it used (Increment.high_level). With the same results for both levels the correction is exactly the single-level one.

The molecular low level can differ by order in the same way, e.g. r2SCAN-D4 at def2-QZVPPD for the pairs and at def2-TZVPPD for the triples. It applies to the fragments of that order assigned the molecular low level:

todo = fragments.calculations(
    high_levels={3: "high:3"}, low_levels={3: "molecular:3"}
)
# todo["molecular"]: the molecular monomers and pairs
# todo["molecular:3"]: the molecular triples with all their sub-pairs and monomers
correction = seamm_mbe.mbe_correction(
    fragments, high, periodic, molecular,
    high_by_order={3: high_tz}, low_by_order={3: low_tz},
)

Each increment records its low level too (Increment.low_level).

Different levels for different orders are not a shortcut: they follow from the decomposition. Each order’s increment is a difference within its own ladder (a triple minus its pairs plus its monomers, all at the same levels), so an order’s levels need only be consistent within that ladder, never across orders. A single low level at a smaller basis would instead put the basis error of every pair into the sum.

Ion shells#

Around Li⁺, the expansion over single molecules converges slowly. Treating the ion and its first shell as one unit fixes it (campaign 2026-10-09). seamm_mbe.ion_shells() returns a seamm_mbe.UnitSystem: a view of a seamm_mbe.System whose molecules are the units, each shell (whole around its ion) and every other molecule. Fragments are then enumerated from it as usual:

system = seamm_mbe.System.from_configuration(configuration)
units = seamm_mbe.ion_shells(system)  # Li-O and Li-F 2.6 Å, nearest ion, <= 5
fragments = seamm_mbe.enumerate_fragments(units, rules)

How the units are built:

Membership: seamm_mbe.IonShellRules sets the cutoffs per ion element and partner element. A molecule within reach of two ions joins the nearer one. max_members keeps a crowded shell’s nearest molecules.
Types: a shell’s type is named by its contents, e.g. Li+[DMC3 PF6-], and its charge is the sum of its members’.
What stays per molecule: the real molecules are on units.parent, for what is per molecule (energy offsets and the molecular virial).
Truncation: SelectionRules.shell_max_order = 2 selects a shell’s pairs but none of its triples.
Criterion: shells need a contact criterion (contact or heavy contact). A shell’s reference point is its ion, so a point criterion would miss partners near its members, and enumerate_fragments refuses it.
Levels of their own: the fragments holding a shell can have their own high and molecular low levels, with the key "shell" (seamm_mbe.SHELL) in the same maps as the per-order levels above. For example, the shell fragments can use revDSD and r2SCAN-D4 at def2-TZVPPD while the solvent pairs stay at def2-QZVPPD:
```
todo = fragments.calculations(
    high_levels={3: "high:3", "shell": "high:shell"},
    low_levels={3: "molecular:3", "shell": "molecular:shell"},
)
# todo["high:shell"]: every selected fragment holding a shell, with its
# sub-fragments (the shell, and the partner monomers, at this level too)
correction = seamm_mbe.mbe_correction(
    fragments, high, periodic, molecular,
    high_by_order={3: high_tz, "shell": high_shell},
    low_by_order={3: low_tz, "shell": low_shell},
)
```
The shell key applies to the fragments of two or more units and takes precedence over the order, so a shell triple (without truncation) uses the shell level. A shell’s own monomer keeps the 1-body level: the 1-body terms set the labels’ absolute energy scale, and its energy offsets, so they must be at one level for every unit. (Its monomer is also computed at the shell level, inside the shell pairs’ increments.) As for the orders, each such increment is a difference within its own ladder, and with the same results for every level the correction is exactly the single-level one. seamm_mbe.level_key() gives the entry of a level map a fragment uses.

Variants of a level for single fragments#

A fragment whose calculation cannot be done at its level as the others are – e.g. a 104-atom Li⁺ shell pair for which ORCA’s AutoAux auxiliary basis is linearly dependent – can have its increment built at a variant of that level instead, such as the same method with explicit auxiliary basis sets. variants maps a level’s name to the variant’s name and the fragments that use it:

variants = {"high:shell": ("high:shell@aux", {"d01_09"})}
todo = fragments.calculations(high_levels=..., variants=variants)
# todo["high:shell@aux"]: d01_09 and its sub-fragments; todo["high:shell"]
# no longer has d01_09 itself, only what the other increments need
correction = seamm_mbe.mbe_correction(
    fragments, high, periodic, molecular,
    high_by_order={"shell": high_shell},
    variants=variants,
    variant_results={"high:shell@aux": high_shell_aux},
)

The increment of each listed fragment is built entirely from the variant’s results (the fragment and its sub-fragments), as for a level of its own, so one auxiliary basis is used throughout it. Every other increment, the sub-fragments’ own among them, is unchanged. Low levels take variants the same way.