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

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. SelectionRules.check() refuses a cell too small for the cutoffs, since the same molecules could otherwise form two different fragments.

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.