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, "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)[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#

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)[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).
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.

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)[source]#: The (name, level) of every result the correction needs but lacks. Results that are None count as missing (a failed calculation).

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}}) – Per cell term, and “MBE” for the correction.
per_body ({int: {"energy": eV, "pressure": atm, "count": int}}) – The correction by order.
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.

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)[source]#

What must be computed, by level.

Returns:: {str – Fragment names for “high” (every fragment any increment needs), “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). molecular_everywhere puts every fragment in “molecular”, which costs more but allows a consistency check of the two 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.

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.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.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>)[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.

check(system)[source]#

Raise SelectionError unless the selection is well defined.

A fragment is named by its molecules (and, for a pair, the image), so two different physical fragments made of the same molecules would collide and one would be silently lost. This never warns: it refuses.

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.

The same bound for pairs (2c) also means a pair has at most one image within the cutoff. 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#

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/.