seamm_mbe package#

Submodules#

seamm_mbe.algebra module#

The many-body increments and the MBE correction: energy, forces, virial.

For a fragment F (a set of molecules at definite images) computed in isolation at the high level and at a low level,

delta_F = E_high(F) - E_low(F)

and its increment is the F-body part, with every proper sub-fragment’s increment removed (Möbius inversion over the subsets):

dE_F = delta_F - sum over proper subsets S of F of dE_S

i.e. dE_i = delta_i, dE_ij = delta_ij - dE_i - dE_j, and so on. Forces go the same way, atom by atom through the slots. The correction is the sum of the selected fragments’ increments, each at its own low level.

Each increment’s virial, W_F = sum_a r_a (x) (f_a - <f>_F), uses its forces less their mean, which makes it independent of the origin (a true increment has zero net force; the residual is numerical but would otherwise be multiplied by an arbitrary origin), with r_a the fragment’s own coordinates.

Units: energies eV, forces eV/Å, virials eV (see seamm_mbe.units).

class seamm_mbe.algebra.Correction(energy: float, forces: ndarray, virial: ndarray, increments: dict = <factory>, per_body: dict = <factory>, per_level: dict = <factory>, max_increment_net_force: float = 0.0)[source]#

Bases: object

The MBE correction of a system.

Variables:

energy (float) – eV.
forces (numpy.ndarray) – (n_atoms, 3) eV/Å, on the system’s atoms.
virial (numpy.ndarray) – (3, 3) eV.
increments ({str: Increment}) – The selected fragments’ increments, each at its own level.
per_body ({int: {"energy": eV, "virial": (3, 3) eV, "forces": (n_atoms, 3)) – eV/Å, “count”: int}} The sums by order.
per_level ({str: {int: int}}) – How many increments of each order used each level.
max_increment_net_force (float) – The largest component of any increment’s net force (eV/Å).

energy: float#

forces: ndarray#

increments: dict#

max_increment_net_force: float = 0.0#

per_body: dict#

per_level: dict#

virial: ndarray#

class seamm_mbe.algebra.FragmentResult(energy: float, forces: ndarray)[source]#

Bases: object

One fragment calculation: energy (eV) and forces (eV/Å, the fragment’s atom order).

energy: float#

forces: ndarray#

class seamm_mbe.algebra.Increment(name: str, order: int, level: str, energy: float, forces: ndarray, net_force: ndarray, virial: ndarray, high_level: str = 'high', low_level: str = '')[source]#

Bases: object

One fragment’s many-body increment.

Variables:

name (str)
order (int)
level (str) – The low level it is referenced to.
energy (float) – eV.
forces (numpy.ndarray) – (n_atoms, 3) eV/Å in the fragment’s atom order.
net_force (numpy.ndarray) – (3,) eV/Å, the sum of the forces: ideally zero.
virial (numpy.ndarray) – (3, 3) eV, origin-independent (net force removed).

energy: float#

forces: ndarray#

high_level: str = 'high'#

level: str#

low_level: str = ''#

name: str#

net_force: ndarray#

order: int#

virial: ndarray#

exception seamm_mbe.algebra.MissingFragmentsError(missing)[source]#

Bases: ValueError

Results are missing or unusable, so the correction cannot be made.

Variables:: missing ([(str, str)]) – (fragment name, level) of each missing result.

seamm_mbe.algebra.increments(fragments, names, high, low)[source]#

The increments of the named fragments from high- and low-level results.

Parameters:

fragments (seamm_mbe.FragmentSet)
names ([str]) – Fragments to compute, closed under sub-fragments, ascending in order (as FragmentSet.calculations() gives them).
high, low ({str: FragmentResult or (energy, forces) or dict}) – Results by fragment name: energy in eV, forces in eV/Å.

Returns:

{str – Each fragment’s increment: energy (eV) and forces (eV/Å).

Return type:

(float, numpy.ndarray)}

seamm_mbe.algebra.mbe_correction(fragments, high, periodic=None, molecular=None, corrections=None, high_by_order=None, low_by_order=None, variants=None, variant_results=None)[source]#

The MBE correction from the fragment results.

Parameters:

fragments (seamm_mbe.FragmentSet) – With levels assigned (seamm_mbe.assign_levels()).
high ({str: result}) – High-level results by fragment name.
periodic, molecular ({str: result}) – Low-level results by fragment name, for the fragments FragmentSet.calculations() lists at each level. A result is a FragmentResult, an (energy, forces) pair or a dict with “energy” and “forces”: eV and eV/Å in the fragment’s atom order (its molecules in slot order, each molecule’s atoms ascending).
high_by_order ({int: {str: result}} or None) – High-level results for orders that have a high level of their own, e.g. {3: tz_results} with high at a larger basis for the monomers and pairs. Each order’s increments are built entirely from its own level’s results (the triple, its pairs and its monomers at that level), so the orders listed need their sub-fragments at that level too (FragmentSet.calculations() with high_levels). With equal results this gives exactly what high alone gives.
low_by_order ({int: {str: result}} or None) – Likewise molecular low-level results for orders with a low level of their own, e.g. {3: tz_results} with molecular at a larger basis for the pairs. They apply to that order’s fragments assigned the molecular low level, each increment built from the triple, its pairs and its monomers at that level (FragmentSet.calculations() with low_levels). An order’s increment is a difference within its own ladder, so its levels need only be consistent within it, never across orders. With equal results this gives exactly what molecular alone gives.
high_by_order, low_by_order (keyed “shell”) – Either map may also have the key "shell" (seamm_mbe.SHELL): the results for every selected fragment of two or more units holding an ion shell, which then takes that level instead of its order’s (e.g. the shell pairs at def2-TZVPPD with the solvent pairs at def2-QZVPPD). Each such increment is built entirely from that level’s results, so its sub-fragments, the shell’s monomer among them, are needed at that level too. The shell’s own 1-body increment stays at the order-1 level, which sets the absolute energy scale for every unit alike.
corrections ({str: result} or None) – Corrections, each a FragmentResult, an (energy, forces) pair or a dict with “energy” and “forces”, added to selected fragments’ increments in the sum only, never to the sub-fragment increments that higher fragments subtract (eV and eV/Å, as the results). This is how a pairwise counterpoise correction enters: with dE_ij^CP - dE_ij for each pair, the sum is E(1) + sum dE_ij^CP + sum dE_ijk, the triples still subtracting the uncorrected pairs, so a pair’s BSSE does not move into the 3-body terms.
variants ({str: (str, set of str)} or None) – Variants of levels for some selected fragments, by the level’s name, as for FragmentSet.calculations(): each listed fragment’s increment is built entirely from the variant’s results (the fragment and its sub-fragments at that level), as for a level of its own. The other fragments, the sub-fragments’ own increments among them, are unchanged.
variant_results ({str: {str: result}} or None) – The results of each variant level, by its name.

Return type:

Correction

Raises:

MissingFragmentsError – If any needed result is missing: an incomplete configuration never gives a correction.

seamm_mbe.algebra.missing_results(fragments, high, periodic=None, molecular=None, high_by_order=None, low_by_order=None, variants=None, variant_results=None)[source]#: The (name, level) of every result the correction needs but lacks. Results that are None count as missing (a failed calculation). high_by_order, low_by_order, variants and variant_results as for mbe_correction(); their levels are named “high:<order>” and “molecular:<order>” (“high:shell” and “molecular:shell” for the key “shell”), and the variants by their own names.

seamm_mbe.assemble module#

Assembling the labels of a configuration: the cell’s own terms plus the MBE correction, the energy reference, stress and pressures.

E = sum of cell terms + dE_MBE (+ per-molecule-type energy offsets) F = sum of cell terms’ forces + F_MBE W = sum of cell terms’ virials + W_MBE

The cell terms are whatever the step computed on the whole cell: the periodic low level, and its dispersion add-on (dftd4 r2SCAN-D4 for VASP) as a separate term. This library never runs a calculation.

Units: eV, eV/Å, virial eV, stress eV/Å³ (sigma = -P), pressures atm.

class seamm_mbe.assemble.CellTerm(label: str, energy: float, forces: ndarray, virial: ndarray | None = None)[source]#

Bases: object

A whole-cell contribution.

Variables:

label (str) – e.g. “VASP r2SCAN”, “D4”.
energy (float) – eV.
forces (numpy.ndarray) – (n_atoms, 3) eV/Å, the system’s atom order.
virial (numpy.ndarray or None) – (3, 3) eV (W = V * P; see seamm_mbe.units for converting a code’s stress); None for a cluster.

energy: float#

forces: ndarray#

label: str#

virial: ndarray | None = None#

class seamm_mbe.assemble.Labels(energy: float, reference_energy: float, forces: ndarray, virial: ndarray | None, stress: ndarray | None, pressure: float | None, molecular_pressure: float | None, breakdown: dict = <factory>, per_body: dict = <factory>, per_level: dict = <factory>, max_increment_net_force: float = 0.0, net_force: float = 0.0)[source]#

Bases: object

The labels of one configuration.

Variables:

energy (float) – Total energy, eV, on the high level’s absolute scale.
reference_energy (float) – energy plus the per-type offsets (eV): e.g. the formation-energy scale of the training sets. Equal to energy without offsets.
forces (numpy.ndarray) – (n_atoms, 3) eV/Å.
virial (numpy.ndarray or None) – (3, 3) eV.
stress (numpy.ndarray or None) – (3, 3) eV/Å³, sigma = -W/V (ASE/xnn).
pressure (float or None) – Atomic configurational pressure tr(W)/3V, atm.
molecular_pressure (float or None) – tr(W - W_intra)/3V, atm (W_intra about each molecule’s centre of mass).
breakdown ({str: {"energy": eV, "pressure": atm, "molecular pressure": atm}}) – Per cell term, and “MBE” for the correction.
per_body ({int: {"energy": eV, "pressure": atm, "molecular pressure": atm,) – “count”: int}} The correction by order. Each order’s molecular pressure is its virial less its own intramolecular part, so the orders add up to the MBE correction’s, and the monomers’ is zero (their increments have no net force on any molecule).
per_level ({str: {int: int}}) – Increments per low level and order.
max_increment_net_force (float) – eV/Å, the largest component of any increment’s net force.
net_force (float) – eV/Å, the norm of the total net force.

breakdown: dict#

energy: float#

forces: ndarray#

max_increment_net_force: float = 0.0#

molecular_pressure: float | None#

net_force: float = 0.0#

per_body: dict#

per_level: dict#

pressure: float | None#

reference_energy: float#

stress: ndarray | None#

virial: ndarray | None#

seamm_mbe.assemble.assemble(system, correction, cell_terms, offsets=None)[source]#

The labels of a configuration.

Parameters:

system (seamm_mbe.System)
correction (seamm_mbe.Correction) – From seamm_mbe.mbe_correction().
cell_terms ([CellTerm]) – The whole-cell calculations (eV, eV/Å, virial eV).
offsets ({str: float} or None) – Energy offset per molecule, by type name (eV).

Return type:

Labels

seamm_mbe.assemble.energy_offset(system, offsets)[source]#: The total offset (eV) for a system from per-type offsets {type: eV per molecule}; every type present needs one.

seamm_mbe.assemble.intramolecular_virial(system, forces)[source]#

W_intra = sum over molecules, sum over their atoms, (r_a - R_com) (x) f_a, with each molecule whole (eV, from forces in eV/Å). The molecular virial is W - W_intra; monomer increments drop out of it exactly.

Applied to a part of the correction (a body order), W_intra uses the raw forces while each increment’s virial uses forces centred on its fragment (net force removed). The two differ by sum over molecules of (sum_a (r_a - R_com)) (x) <f>, nonzero only because R_com is mass-weighted; it scales with the increments’ net forces (max_increment_net_force), so it is numerical noise. The fragment frame (each molecule at the image it interacts at) is the right one for the increments’ virial.

seamm_mbe.catalog module#

Molecule types: identification by formula and bond-graph topology.

A molecule’s signature is its Hill formula plus a Weisfeiler–Lehman hash of its element-labelled bond graph, so isomers with the same formula are told apart (in practice; WL refinement does not separate every pair of non-isomorphic graphs, which is irrelevant for the small molecules here). The same refinement labels each atom by its symmetry class within the molecule, which is how a type’s designated atom (the atom whose position defines the molecule’s position for the distance criterion: water’s O, a carbonate’s carbonyl C, an ion’s central atom) is found in every molecule of that type.

The built-in CATALOG knows the molecules of the water/electrolyte campaign: water, ethylene carbonate (EC), fluoroethylene carbonate (FEC), dimethyl carbonate (DMC), ethyl methyl carbonate (EMC), Li⁺, BF₄⁻ and PF₆⁻, plus the common monatomic ions Na⁺, K⁺, F⁻, Cl⁻, Br⁻ and I⁻. Anything else is typed automatically, named by its formula.

seamm_mbe.catalog.CATALOG = {'BF4:cc066cca103b98c9': MoleculeType(name='BF4-', formula='BF4', signature='BF4:cc066cca103b98c9', charge=-1, multiplicity=1, designated='fbebf833ae16648d'), 'Br:1d6664f594d8f8ab': MoleculeType(name='Br-', formula='Br', signature='Br:1d6664f594d8f8ab', charge=-1, multiplicity=1, designated='681ebfaf1c1b0b85'), 'C3H3FO3:bc2862c3458977e4': MoleculeType(name='FEC', formula='C3H3FO3', signature='C3H3FO3:bc2862c3458977e4', charge=0, multiplicity=1, designated='a74d3322c641e935'), 'C3H4O3:404341a3a83e87f3': MoleculeType(name='EC', formula='C3H4O3', signature='C3H4O3:404341a3a83e87f3', charge=0, multiplicity=1, designated='4252ca90f028ff0b'), 'C3H6O3:7cb88292f5e76c55': MoleculeType(name='DMC', formula='C3H6O3', signature='C3H6O3:7cb88292f5e76c55', charge=0, multiplicity=1, designated='4252ca90f028ff0b'), 'C4H8O3:70775b073c6b4e26': MoleculeType(name='EMC', formula='C4H8O3', signature='C4H8O3:70775b073c6b4e26', charge=0, multiplicity=1, designated='1b8054939067fdfd'), 'Cl:2c1871835f6d5cf7': MoleculeType(name='Cl-', formula='Cl', signature='Cl:2c1871835f6d5cf7', charge=-1, multiplicity=1, designated='335432670dfa4ef4'), 'F6P:43cec40088a0d55f': MoleculeType(name='PF6-', formula='F6P', signature='F6P:43cec40088a0d55f', charge=-1, multiplicity=1, designated='134da7ef5a3ff162'), 'F:51f2c960f6108a4c': MoleculeType(name='F-', formula='F', signature='F:51f2c960f6108a4c', charge=-1, multiplicity=1, designated='2a45c6ebf2a286eb'), 'H2O:e62d91c53f6ca05c': MoleculeType(name='water', formula='H2O', signature='H2O:e62d91c53f6ca05c', charge=0, multiplicity=1, designated='9a250af442bfb3b3'), 'I:9fad73a86f97286b': MoleculeType(name='I-', formula='I', signature='I:9fad73a86f97286b', charge=-1, multiplicity=1, designated='3187cf6c5fa3c155'), 'K:c263daef2ca58b22': MoleculeType(name='K+', formula='K', signature='K:c263daef2ca58b22', charge=1, multiplicity=1, designated='b48d3ad274f282a2'), 'Li:17fb4b60b90c6635': MoleculeType(name='Li+', formula='Li', signature='Li:17fb4b60b90c6635', charge=1, multiplicity=1, designated='ee6e9886b2ca901e'), 'Na:88f2174d333c6278': MoleculeType(name='Na+', formula='Na', signature='Na:88f2174d333c6278', charge=1, multiplicity=1, designated='772aac09e5178690')}#: The built-in molecule types, keyed by signature

class seamm_mbe.catalog.MoleculeType(name: str, formula: str, signature: str, charge: int | None = 0, multiplicity: int = 1, designated: str | None = None)[source]#

Bases: object

A kind of molecule.

Variables:

name (str) – The type’s name, e.g. “water”, “EC”, “Li+”, or the formula for an automatically typed molecule.
formula (str) – Hill formula.
signature (str) – Formula + WL graph hash (see signature()).
charge (int or None) – The molecule’s total charge; None if the type does not fix it.
multiplicity (int) – Spin multiplicity, 2S + 1.
designated (str or None) – The WL label of the designated atom, or None to use the heavy atom nearest the centre of mass (set when the type is first seen).

charge: int | None = 0#

designated: str | None = None#

formula: str#

multiplicity: int = 1#

name: str#

signature: str#

seamm_mbe.catalog.define_type(name, symbols, bonds, charge=0, multiplicity=1, designated=None)[source]#

A MoleculeType from a template graph.

Parameters:

name (str) – The type’s name.
symbols ([str]) – Element symbols of the template.
bonds ([(int, int)]) – The template’s bonds (0-based).
charge, multiplicity (int) – Defaults for the type.
designated (int or None) – Index of the designated atom in the template.

seamm_mbe.catalog.signature(symbols, bonds)[source]#: The topology signature of a molecule: formula + WL graph hash.

seamm_mbe.catalog.wl_labels(symbols, bonds)[source]#

The Weisfeiler–Lehman atom labels of an element-labelled graph.

Parameters:

symbols ([str]) – Element symbols of the molecule’s atoms.
bonds ([(int, int)]) – Bonds as pairs of 0-based indices into symbols.

Returns:

A label per atom; atoms with equal labels are (WL-)equivalent.

Return type:

[str]

seamm_mbe.elements module#

Element symbols, atomic numbers and standard atomic weights.

The weights are the IUPAC abridged standard atomic weights (conventional values for elements with an interval), used only for molecular centres of mass. They are kept here, for H to Ba, so that common systems need nothing beyond numpy; heavier elements fall back to molsystem’s table.

Keep the conventional values (O 15.999, C 12.011, Li 6.94): do not replace them with molsystem’s interval midpoints (O 15.9995, C 12.0105, Li 6.9675). The molecular pressure depends on the centres of mass, and the prototype that the regression test reproduces (tests/data) used the conventional values; the midpoints move P_mol of the pilot frame.

seamm_mbe.elements.IONIC_ELEMENTS = frozenset({'Ba', 'Be', 'Ca', 'Cs', 'K', 'Li', 'Mg', 'Na', 'Rb', 'Sr'})#

never bonded (as molsystem)

Type:: Elements that are ions in molecular systems

seamm_mbe.elements.MONATOMIC = frozenset({'Ar', 'Ba', 'Be', 'Br', 'Ca', 'Cl', 'Cs', 'F', 'He', 'I', 'K', 'Kr', 'Li', 'Mg', 'Na', 'Ne', 'Rb', 'Rn', 'Sr', 'Xe'})#

those ions, the halide ions and the noble gases

Type:: Elements that can be a molecule on their own

seamm_mbe.elements.atomic_number(symbol)[source]#: The atomic number of an element symbol.

seamm_mbe.elements.hill_formula(symbols)[source]#: The Hill-order formula, e.g. ‘H2O’, ‘C3H4O3’, ‘BF4’, ‘Li’.

seamm_mbe.elements.mass(symbol)[source]#

The standard atomic weight (g/mol) of an element symbol.

The table here (H–Ba) uses IUPAC’s conventional values, as the prototype did (O 15.999); heavier elements come from molsystem’s table.

seamm_mbe.fragments module#

Fragment enumeration with periodic images and canonical keys.

A fragment is a set of whole molecules, each at a definite periodic image. Its key is canonical: the sorted molecule indices plus the integer images of molecules 2..n relative to molecule 1, so the same physical fragment reached two ways (a pair selected for itself and as part of a triple) has one key and is computed once.

Its name is the prototype’s (gen_stage2_frame.py), so results stored under the prototype’s names match:

monomers m00;
pairs d00_17, with _x and one character per cell axis (m, 0, p for -1, 0, +1) appended when the image of the second molecule is not its minimum image relative to the first, e.g. d03_41_x0p0 (images beyond +-1, never needed so far, are written out: _x2_0_-1);
selected fragments of order 3 and 4 t00_17_30, q00_17_30_41 – without images, which is unique because SelectionRules.check() guarantees the same molecules cannot form two selected fragments;
auxiliary fragments of order >= 3 (sub-fragments of a selected fragment that are not selected themselves) carry the images of molecules 2..n relative to molecule 1 in the same _x form in a periodic system (in a cluster there are no images, and no suffix).

Every selected fragment’s proper sub-fragments are computed too, because its increment needs them (the prototype’s “outer pair of a connected triple”): they are in the set with in_sum = False unless selected themselves.

Each fragment has a frame: the images its molecules are placed at for the calculation. Monomers are at home; pairs and auxiliary fragments have their first molecule at home; selected fragments of order >= 3 have their hub (the first molecule, in index order, bonded to the most others) at home, as the prototype did. Frames differ only by lattice translations, which change nothing physical.

class seamm_mbe.fragments.Fragment(name: str, key: tuple, molecules: tuple, images: tuple, in_sum: bool, coordinates: ndarray, atoms: ndarray, slot_atoms: list, distances: dict, charge: int, multiplicities: tuple, subfragments: list = <factory>, level: str | None = None)[source]#

Bases: object

A fragment of whole molecules.

Variables:

name (str) – The prototype-compatible name (see the module docstring), used as the task key.
key (tuple) – The canonical key (see canonical_key()).
molecules (tuple of int) – The molecules, ascending: the fragment’s slot order.
images (tuple of (int, int, int)) – The image of each molecule in the fragment’s frame.
in_sum (bool) – Whether its increment enters the energy (selected), or it is only a sub-fragment of selected ones.
coordinates (numpy.ndarray) – (n_atoms, 3) Cartesian coordinates in Å in the frame; the atoms are the molecules’ atoms, molecule by molecule in slot order.
atoms (numpy.ndarray of int) – The system’s atom index of each fragment atom.
slot_atoms ([numpy.ndarray]) – Per slot, the positions (rows) of that molecule’s atoms in the fragment – with atoms, what a counterpoise ghost set needs.
distances ({(int, int): float}) – The criterion distance (Å) between each pair of slots, in the frame.
charge (int) – Total charge.
multiplicities (tuple of int) – Each molecule’s multiplicity.
subfragments ([(str, tuple of int)]) – Every proper sub-fragment (all non-empty proper subsets of the molecules) by name, with the slots of this fragment it occupies.
level (str or None) – For a selected fragment: which low level its increment uses, “periodic” or “molecular” (see seamm_mbe.assign_levels()).

atoms: ndarray#

charge: int#

coordinates: ndarray#

distances: dict#

images: tuple#

in_sum: bool#

key: tuple#

level: str | None = None#

property max_distance#: The largest criterion distance (Å) between two members (0 for a monomer).

molecules: tuple#

multiplicities: tuple#

property multiplicity#

a molecule’s own, or 1 when every molecule is closed-shell. Open-shell molecules in a fragment of several are not supported (their coupling is ambiguous).

Type:: The fragment’s multiplicity

property n_atoms#

name: str#

property order#: The number of molecules.

slot_atoms: list#

subfragments: list#

symbols(system)[source]#: The element symbols of the fragment’s atoms.

class seamm_mbe.fragments.FragmentSet(system, rules, fragments)[source]#

Bases: object

The fragments of a system, in a stable order: by order, then by key (the prototype’s order for monomers, pairs and selected triples).

Use enumerate_fragments() to make one.

by_order(order, in_sum=None)[source]#: The fragments of an order; only selected (True) or only auxiliary (False) ones if in_sum is given.

calculations(molecular_everywhere=False, high_levels=None, low_levels=None, variants=None)[source]#

What must be computed, by level.

Parameters:

molecular_everywhere (bool) – Put every fragment in “molecular”, which costs more but allows a consistency check of the two low levels.
high_levels ({int or str: str} or None) – A high level other than “high” for some orders, e.g. {3: "high:3"} for triples at their own level (a smaller basis than the pairs’). Each order’s increments are built entirely at its own high level, so that level is computed on those fragments and all their sub-fragments: the monomers and sub-pairs of the triples are computed at both levels. The key "shell" (SHELL) gives the level of every selected fragment of two or more units holding an ion shell, e.g. {"shell": "high:shell"} for the shell pairs at a smaller basis; the shell’s own monomer keeps the 1-body level (see level_key()).
low_levels ({int or str: str} or None) – Likewise a molecular low level other than “molecular” for some orders, or for the fragments holding an ion shell ("shell"), e.g. {3: "molecular:3"} for triples at a smaller basis than the pairs’ low level. It applies to those fragments assigned the molecular low level, and is computed on them and all their sub-fragments.
variants ({str: (str, set of str)} or None) – For some selected fragments, a variant of one of their levels, by the level’s name: {"high:shell": ("high:shell@aux", {"d01_09"})} builds the increment of d01_09 at the variant “high:shell@aux” instead of “high:shell” (e.g. the same method with explicit auxiliary basis sets where AutoAux failed). The variant level is computed on those fragments and all their sub-fragments; the base level no longer needs the fragments themselves, only what its other increments need.

Returns:

{str – Fragment names for “periodic” and “molecular” (the selected fragments assigned to that low level, with all their sub-fragments: each increment is built from sub-fragments at its own level), for “high” (what the increments of the orders without a level of their own need) and for each level in high_levels and low_levels.

Return type:

[str]}

closure(fragments)[source]#: The names of the given fragments and all their sub-fragments, in set order.

counts()[source]#: {order: {“selected”: n, “auxiliary”: n}}.

property names#

selected()[source]#: The fragments whose increments enter the sum.

property shells#: The units that are ion shells (empty without ion shells).

seamm_mbe.fragments.MAX_IMAGE_LAYERS = 3#: The most layers of periodic images the neighbour search looks through ((2k + 1)^3 cells); a cell needing more is refused

seamm_mbe.fragments.SHELL = 'shell'#: The key of a level map (high_levels, low_levels, high_by_order, low_by_order) for the fragments holding an ion shell.

seamm_mbe.fragments.canonical_key(placement)[source]#

The canonical key of a fragment from its placement.

Parameters:: placement ({int: (int, int, int)}) – The image of each molecule, in any common frame.
Returns:: (molecule indices ascending, images of molecules 2..n relative to molecule 1).
Return type:: tuple

seamm_mbe.fragments.check_variants(fragments, variants)[source]#: The level variants (see FragmentSet.calculations()), checked: each names selected fragments and a variant level of a different name.

seamm_mbe.fragments.enumerate_fragments(system, rules=None)[source]#

Enumerate the fragments of a system.

Parameters:

system (seamm_mbe.System) – The molecules, with or without a cell.
rules (seamm_mbe.SelectionRules) – The selection; default the prototype’s (pairs < 4.5 Å, connected triples < 3.5 Å, designated-atom distances).

Return type:

FragmentSet

Raises:

SelectionError – If the selection is not well defined for the cell (see SelectionRules.check()).

seamm_mbe.fragments.level_key(fragment, levels, default, shells=frozenset({}))[source]#

Which entry of a level map a selected fragment’s increment uses.

Parameters:

fragment (Fragment)
levels ({int or str: object} or None) – By order, and for the fragments of two or more units holding an ion shell by SHELL, which takes precedence over the order. A shell’s own monomer keeps its order’s level: the 1-body terms set the labels’ absolute energy scale (and its energy offsets), so they stay at one level for every unit.
default (object) – For a fragment the map does not cover.
shells (frozenset of int) – The units that are ion shells (FragmentSet.shells).

Return type:

object

seamm_mbe.fragments.low_level_of(fragment, low_levels=None, shells=frozenset({}))[source]#: The low level a selected fragment’s increment uses: its assigned level, or for a molecular one its own molecular level, if it has one (low_levels, by order, e.g. {3: "molecular:3"}, or by SHELL for the fragments holding one of shells).

seamm_mbe.fragments.variant_of(fragment, level, variants=None)[source]#: The level a selected fragment uses: level, or its variant if variants (see FragmentSet.calculations()) lists the fragment for that level.

seamm_mbe.levels module#

Which low level each selected fragment’s increment uses.

The cell’s low level is usually a periodic plane-wave code, and the fragments’ low level need not be the same code. The prototype’s mixed scheme references monomers and compact pairs (O–O < 3.5 Å) to the periodic code (VASP, in registered boxes), cancelling the cell’s own low-level error best, and everything else to a molecular code (ORCA r2SCAN-D4), because extended fragments pick up image interactions in affordable periodic boxes.

Each increment is built entirely at its own level: a triple referenced to the molecular level subtracts the molecular-level increments of its pairs, even when those pairs are themselves referenced to the periodic level in the sum. So the periodic level is computed on the periodic fragments and their sub-fragments, and the molecular level on the molecular fragments and theirs (FragmentSet.calculations()).

seamm_mbe.levels.assign_levels(fragments, periodic=None)[source]#

Assign each selected fragment a low level.

Parameters:

fragments (seamm_mbe.FragmentSet) – The fragments; their level attributes are set.
periodic ({int: bool or float} or None) – Per order, whether selected fragments use the periodic low level: True for all of that order, a distance r (Å) for those whose members are all closer than r (by the selection’s criterion), False or absent for none. Only real booleans mean all or none: 1 is a distance of 1 Å. None (the default) puts everything at the molecular level. The prototype’s mixed scheme is {1: True, 2: 3.5}.

Returns:

{str – The number of selected fragments per level and order.

Return type:

{int: int}}

seamm_mbe.selection module#

Which fragments are selected: the distance criterion, the cutoffs and the rules for each order.

seamm_mbe.selection.CRITERIA = ('designated', 'com', 'cog', 'contact', 'heavy contact')#: How the distance between two molecules is measured

seamm_mbe.selection.RULES = ('none', 'connected', 'hub', 'compact')#: The rules for fragments of order 3 and higher

exception seamm_mbe.selection.SelectionError[source]#

Bases: ValueError

The selection is not well defined for this system.

class seamm_mbe.selection.SelectionRules(max_order: int = 3, criterion: str = 'designated', cutoffs: dict = <factory>, rules: dict = <factory>, shell_max_order: int | None = None)[source]#

Bases: object

The rules selecting the fragments.

Variables:

max_order (int) – The highest order of fragment (1 = monomers only, 2 = pairs, …).
criterion (str) – The intermolecular distance (Å): “designated” (between the types’ designated atoms, e.g. water O–O), “com” or “cog” (centres of mass or geometry), “contact” (closest atom–atom contact) or “heavy contact” (closest contact between non-hydrogen atoms).
cutoffs ({int: float or {(str, str): float}}) – Per order, the cutoff (Å): one value, or a table by molecule-type pair. A table’s keys are (type, type) tuples in either order; “*” matches any type, and the most specific entry wins. For pairs the cutoff selects the pairs; for order n >= 3 it defines the bonds of the connectivity graph that the order’s rule tests.
rules ({int: str}) – Per order n >= 3: “none” (no fragments of that order), “connected” (the n molecules form a connected graph), “hub” (one molecule is bonded to all the others) or “compact” (all of them bonded to each other). For n = 3, “connected” and “hub” are the same.
shell_max_order (int or None) – The highest order of a selected fragment holding an ion shell (see seamm_mbe.shells); None for the same as every fragment. With 2, a shell’s pairs are selected but none of its triples.

check(system)[source]#

Warn when the same molecules may form several different fragments.

The neighbour search looks as many image layers out as the cutoffs need, so a cell narrower than a cutoff (plus, for the contact criteria, a molecule’s diameter) is handled; only a molecule within the pair cutoff of its own image is refused, by the enumeration, since no fragment can hold that interaction. Narrower than the bound below, the same molecules can form two different fragments through different images. Those are genuinely different fragments of the lattice and the enumeration keeps both, naming all but the minimum-image one with its images, so this only warns.

The bound. Let a selected fragment of order n have, under its rule, radius at most R (some member A is within R of every member) and diameter at most D (no two members further apart than D), measured with the criterion’s distance. Suppose two different selected fragments F1 and F2 contain the same molecules. Put both with molecule A, F1’s centre, at the same position. Some member k then sits at x in F1 and at x + T in F2, T a non-zero lattice vector, so:

length(T) <= d1(A, k) + d2(A, k) <= R + D.

Every non-zero lattice vector is at least as long as the smallest perpendicular width of the cell, L_min. So L_min > R + D makes the collision impossible. With c the order’s largest cutoff:

pairs (one bond): R = c, D = c -> 2c
“compact” (all bonded), any n: R = c, D = c -> 2c
“hub” (one bonded to all): R = c, D = 2c -> 3c
“connected”, n molecules: a spanning tree has a centre within floor(n/2) bonds of every member and a diameter of at most n - 1 bonds: R = floor(n/2) c, D = (n-1) c -> 3c for n = 3, 5c for n = 4.

For the contact criteria the bound is on the reference points (centres of geometry), so each bond length becomes c + 2 r_max, r_max being the largest distance of an atom from its molecule’s centre of geometry. Clusters (no cell) are always fine.

criterion: str = 'designated'#

cutoff(order, type_a, type_b)[source]#: The cutoff (Å) of an order for a pair of molecule types.

cutoffs: dict#

max_order: int = 3#

reach(order)[source]#: The bound on the radius R and diameter D of a selected fragment of an order, in units of its cutoff (see check()).

rule(order)[source]#: The rule for an order: “pair” for 2, else the configured rule.

rules: dict#

shell_max_order: int | None = None#

seamm_mbe.shells module#

Ion shells: an ion and its first-shell molecules as one unit of the expansion.

Around a strongly polarizing ion (Li⁺) the per-molecule expansion converges slowly: in the C.3 cluster (Li⁺/BF₄⁻/6 H₂O/4 EC, 2026-10-09) the Li⁺ 3- and 4-body [revDSD − r2SCAN] increments summed to +19 to +23 and −5 to −18 kJ/mol. With Li⁺ and its first shell as one unit the error at three units fell to +1–2 kJ/mol and the forces to 3.6 meV/Å RMS (seamm_mbe campaign 2026-10-09).

IonShellRules says which ions get a shell and how far it reaches. ion_shells() returns a UnitSystem, a view of a System whose molecules are the units: each shell (the ion and its members, placed whole around the ion) and every other molecule by itself. enumerate_fragments() works on it unchanged, so the fragments are built from units. The real molecules stay on UnitSystem.parent, for what is per molecule (energy offsets, the molecular virial).

Membership. A molecule is a candidate for an ion’s shell when one of its atoms of a listed element lies within that element’s cutoff of the ion (e.g. Li–O 2.6 Å, Li–F 2.6 Å). Ions never join a shell (only molecules with listed atoms can), and anions get none of their own, so a contact ion pair puts the anion in the cation’s shell. A molecule near two ions joins the nearer one (shared = "nearest"), so units never overlap; with max_members a crowded shell keeps its nearest members and the rest stay units of their own. The units are rebuilt from each frame’s geometry.

class seamm_mbe.shells.IonShellRules(cutoffs: dict = <factory>, shared: str = 'nearest', max_members: int | None = 5)[source]#

Bases: object

Which ions get a shell, and its reach.

Variables:

cutoffs ({str: {str: float}}) – Per ion element, the shell cutoff (Å) to each partner element: a molecule joins the shell when one of its atoms of that element is within the cutoff of the ion. Default Li⁺ with O and F at 2.6 Å, just inside the first minima of g(r) in carbonate electrolytes.
shared (str) – What happens to a molecule within reach of two ions: “nearest” joins the nearer one.
max_members (int or None) – The most molecules in one shell (the nearest are kept); None for no limit.

cutoffs: dict#

max_members: int | None = 5#

shared: str = 'nearest'#

exception seamm_mbe.shells.ShellError[source]#

Bases: StructureError

The ion-shell rules cannot be applied to this system.

class seamm_mbe.shells.UnitSystem(parent, members, placements, shell_types)[source]#

Bases: object

A System seen as units: shells and other molecules.

It has the attributes and methods of a System that fragment enumeration uses, with molecules and types describing the units. A unit’s atoms are system atom indices, ascending, and its coordinates are whole around the ion.

Variables:

parent (seamm_mbe.System) – The real system.
members ([tuple of int]) – Per unit, the parent’s molecules it holds (the ion first for a shell).
shells (frozenset of int) – The units that are shells.

formula()#: The system’s Hill formula.

property masses#: Atomic masses (g/mol), looked up when first needed.

minimum_image(vector)#: The integer image n making vector + n @ cell shortest, and that shortest vector. Searches the 27 images around the fractional rounding, which is exact for any vector shorter than half the cell’s smallest width (bonds, and everything within a valid cutoff) and for longer vectors in reasonably shaped cells; a strongly skewed (non-reduced) cell can miss the shortest image of a long vector. For a cluster, n = (0, 0, 0).

property periodic#: Whether the system is a periodic cell (False for a cluster).

reference_point(molecule, criterion)#: The point (Å, home coordinates) that locates a molecule for a distance criterion: “designated”, “com” or “cog”. The contact criteria use the centre of geometry as their reference point.

shell_sizes()[source]#: {unit index: number of molecules besides the ion} for the shells.

shift(image)#: The Cartesian translation (Å) of an integer image (n_a, n_b, n_c).

type_counts()#: {type name: number of molecules}.

property volume#: The cell volume in Å³ (None for a cluster).

property widths#

the spacing of the planes spanned by each pair of lattice vectors. The smallest is L_min, the bound on any cutoff. None for a cluster.

Type:: The three perpendicular widths of the cell (Å)

seamm_mbe.shells.ion_shells(system, rules=None)[source]#

The units of a system under ion-shell rules.

Parameters:

system (seamm_mbe.System)
rules (IonShellRules) – Default Li⁺ with O and F at 2.6 Å, nearest ion, at most 5 members.

Returns:

Units numbered by their lowest molecule index, so a system without any shell keeps its molecule numbering.

Return type:

UnitSystem

seamm_mbe.system module#

The system: atoms, cell, whole molecules and their types.

A System is built from arrays (symbols, Cartesian coordinates in Å, the cell, the bonds) or from a molsystem configuration (System.from_configuration()). It finds the molecules as the connected components of the bond graph, makes each molecule whole, types the molecules by formula and topology, and settles each type’s charge and multiplicity.

Molecules are numbered in the order of their first atom, and each molecule’s atoms are in ascending atom order. A molecule’s home coordinates keep its first atom where it is and place the others by following the bonds with the minimum image, so a molecule that is already whole keeps exactly the coordinates it was given.

class seamm_mbe.system.Molecule(index: int, atoms: ndarray, type: str, coordinates: ndarray, designated: int)[source]#

Bases: object

One whole molecule.

Variables:

index (int) – 0-based molecule number.
atoms (numpy.ndarray of int) – The atoms (0-based indices in the system), ascending.
type (str) – The name of its MoleculeType.
coordinates (numpy.ndarray) – (n, 3) home coordinates in Å, whole.
designated (int) – Index into atoms of the designated atom.

atoms: ndarray#

coordinates: ndarray#

designated: int#

index: int#

type: str#

exception seamm_mbe.system.StructureError[source]#

Bases: ValueError

The structure cannot be used as given (bad charges, missing bonds, …).

class seamm_mbe.system.System(symbols, coordinates, cell=None, *, bonds=(), formal_charges=None, catalog=None, charges=None, multiplicities=None, masses=None)[source]#

Bases: object

Atoms, an optional cell, and the molecules.

Parameters:

symbols ([str]) – Element symbols.
coordinates (array-like) – (n, 3) Cartesian coordinates in Å.
cell (array-like or None) – (3, 3) lattice vectors in Å, one per row; None for a cluster.
bonds ([(int, int)]) – Bonds as 0-based atom index pairs. They define the molecules; atoms with no bonds are one-atom molecules (ions).
formal_charges ([int] or None) – Per-atom formal charges, e.g. from an SDF file’s M CHG lines. When given and not all zero they define each molecule’s charge, and a catalog type whose charge disagrees is an error. When absent or all zero (bare xyz input) the catalog’s charges are used.
catalog ({str: MoleculeType} or None) – Known types keyed by signature; default seamm_mbe.CATALOG.
charges ({str: int} or None) – Charges by type name, overriding both the formal charges and the catalog.
multiplicities ({str: int} or None) – Multiplicities by type name, overriding the catalog (default 1).
masses ([float] or None) – Atomic masses (g/mol); default the standard atomic weights.

formula()[source]#: The system’s Hill formula.

classmethod from_configuration(configuration, **kwargs)[source]#

A System from a molsystem configuration.

The configuration must have its bonds (call its perceive_bonds() first if it has none): they define the molecules. Per-atom formal charges, if the configuration carries them, set the molecules’ charges (see System). Keyword arguments are passed to System.

property masses#: Atomic masses (g/mol), looked up when first needed.

minimum_image(vector)[source]#: The integer image n making vector + n @ cell shortest, and that shortest vector. Searches the 27 images around the fractional rounding, which is exact for any vector shorter than half the cell’s smallest width (bonds, and everything within a valid cutoff) and for longer vectors in reasonably shaped cells; a strongly skewed (non-reduced) cell can miss the shortest image of a long vector. For a cluster, n = (0, 0, 0).

property periodic#: Whether the system is a periodic cell (False for a cluster).

reference_point(molecule, criterion)[source]#: The point (Å, home coordinates) that locates a molecule for a distance criterion: “designated”, “com” or “cog”. The contact criteria use the centre of geometry as their reference point.

shift(image)[source]#: The Cartesian translation (Å) of an integer image (n_a, n_b, n_c).

type_counts()[source]#: {type name: number of molecules}.

property volume#: The cell volume in Å³ (None for a cluster).

property widths#

the spacing of the planes spanned by each pair of lattice vectors. The smallest is L_min, the bound on any cutoff. None for a cluster.

Type:: The three perpendicular widths of the cell (Å)

seamm_mbe.units module#

Units and sign conventions.

The library works throughout in

energy: eV
length: Å
forces: eV/Å (forces, not gradients: F = -dE/dr)
virial: eV, the tensor W = sum_a r_a (x) f_a, so that for a configurational (no kinetic) pressure tensor P = W / V
pressure tensor P: eV/Å³, positive = the system pushes outward
stress sigma: eV/Å³ in the ASE/xnn convention sigma = (1/V) dE/d(strain) = -P

SEAMM’s Model Chemistry batch contract (analyze_task) and the step’s store_results use kJ/mol, kJ/mol/Å gradients and GPa; the functions here convert both ways. The conversion layer is where unit bugs live, so every function states its units and the tests check them against pint.

Codes disagree on the sign of the “stress” they print, so nothing here guesses: every stress or pressure from outside comes with an explicit convention, either "pressure" (positive = pushes outward: VASP’s in kB line, MDI’s <STRESS) or "stress" (sigma = -P: ASE, xnn).

seamm_mbe.units.EV_PER_A3_TO_ATM = 1581225.3974833456#: 1 eV/Å³ in atm (1 atm = 101325 Pa)

seamm_mbe.units.EV_PER_A3_TO_GPA = 160.2176634#

1.602176634e-19 J / 1e-30 m³ = 1.602176634e11 Pa)

Type:: 1 eV/Å³ in GPa (exact

seamm_mbe.units.EV_TO_KJ_PER_MOL = 96.48533212331002#

e * N_A / 1000, exact)

Type:: 1 eV in kJ/mol (CODATA 2018

seamm_mbe.units.KBAR_TO_EV_PER_A3 = 0.0006241509074460764#: 1 kbar (VASP’s “kB”) in eV/Å³

seamm_mbe.units.as_tensor(value)[source]#: A (3, 3) array from a (3, 3) array, 9 values (row-major) or Voigt 6 values in the order xx yy zz yz xz xy. Units are unchanged.

seamm_mbe.units.dftd4_virial_to_virial(dE_dstrain)[source]#: The virial (eV) from the dftd4 library’s “virial”, which is dE/d(strain) (convert it from hartree to eV first): W = -dE/d(strain).

seamm_mbe.units.ev_to_kj_per_mol(energy)[source]#: Energy in eV -> kJ/mol.

seamm_mbe.units.forces_to_gradients(forces)[source]#: Forces in eV/Å -> gradients in kJ/mol/Å (note the sign: g = -F).

seamm_mbe.units.from_analyze_task(result, *, volume=None, stress_convention=None)[source]#

Convert a Model Chemistry analyze_task result to the library’s units.

Parameters:

result (dict) – {“energy”: kJ/mol, “gradients”: (n, 3) kJ/mol/Å, and optionally “stress”: GPa as the program gives it}.
volume (float) – The cell volume in Å³; needed only with a stress.
stress_convention (str) – “pressure” or “stress” (see virial_from_tensor()); needed only with a stress, because the contract does not fix the sign.

Returns:

{“energy”: eV, “forces”: (n, 3) eV/Å, and “virial”: (3, 3) eV if the result had a stress}.

Return type:

dict

seamm_mbe.units.gradients_to_forces(gradients)[source]#: Gradients in kJ/mol/Å -> forces in eV/Å (note the sign: F = -dE/dr).

seamm_mbe.units.kj_per_mol_to_ev(energy)[source]#: Energy in kJ/mol -> eV.

seamm_mbe.units.pressure(virial, volume, units='atm')[source]#: The scalar pressure, tr(P)/3, from the virial (eV) and volume (Å³), in units (“atm”, “GPa”, “kbar” or “eV/Å^3”).

seamm_mbe.units.tensor_from_virial(virial, volume, *, convention, units='GPa')[source]#: A stress or pressure tensor from the virial W (eV): the inverse of virial_from_tensor(), with the same arguments. Returns (3, 3) in units.

seamm_mbe.units.to_seamm(energy=None, forces=None, virial=None, volume=None)[source]#

The library’s numbers in SEAMM’s store_results conventions.

energy (eV) -> kJ/mol; forces (eV/Å) -> gradients (kJ/mol/Å); virial (eV) with the volume (Å³) -> stress in GPa, sigma = -P (the ASE/xnn stress, the training-label convention). Returns a dict with the keys given.

seamm_mbe.units.vasp_stress_to_virial(stress_kB, volume)[source]#: The virial (eV) from VASP’s in kB line: XX YY ZZ XY YZ ZX in kbar, positive = pushing outward (a pressure). Note VASP’s order differs from Voigt’s.

seamm_mbe.units.virial_from_tensor(value, volume, *, convention, units='GPa')[source]#

The virial W (eV) from a stress or pressure tensor.

Parameters:

value (array-like) – (3, 3), 9 or Voigt 6 values (xx yy zz yz xz xy).
volume (float) – The cell volume in Å³.
convention (str) – “pressure” (positive = pushes outward; VASP, MDI) or “stress” (sigma = -P; ASE, xnn). Required: there is no default.
units (str) – “GPa”, “kbar” (VASP’s kB), “eV/Å^3” or “atm”.

Returns:

(3, 3) virial in eV, W = V * P.

Return type:

numpy.ndarray

Module contents#

Many-body expansion (MBE) corrections for periodic cells and clusters.

The energy, forces and stress of a system at a high level of theory are estimated as a cheap calculation on the whole system plus many-body increments of [high - low] computed on small isolated fragments (monomers, selected pairs and triples, …). This library holds the bookkeeping, with no calculations and no SEAMM GUI: molecule typing, fragment enumeration with periodic images and canonical keys, the mixed low-level assignment, the increment algebra for energy, forces and virial, and the assembled labels. The mbe_step plug-in runs the calculations. See docs/developer_guide/campaigns/2026-10-03/.