Plan: Model Chemistry Step, Metadata Protocol, and MOPAC/xTB via MDI#

Date:

2026-06-20

Status:

Architecture proposed; metadata classmethod implemented (sketch) for MOPAC against its real existing metadata; xTB metadata sketched pending a new step package; lammps_step consumption logic sketched, not yet implemented.

Overview#

Replaces the ad hoc _forcefield / _OpenKIM_Potential / _pytorch_model global-variable scheme with a single new Model Chemistry Step, storing one dictionary (e.g. as the global variable _model_chemistry) that uniformly describes a classical forcefield, an MLFF, or a model-chemistry string (Program:Type@Method [/Basis[@Cutoff]]) – including, for model chemistries, enough information for a consuming driver step (lammps_step first, others later) to launch the right MDI engine.

This plan covers: the dict shape, the metadata protocol every program-step package implements to describe what it offers, a worked implementation of that protocol against mopac_step’s real existing metadata, a sketch of the same for xTB (which has no existing SEAMM step package at all), and how lammps_step would consume all of this to decide between a classical pair_style and an MDI connection.

The _model_chemistry dictionary#

# Classical forcefield (existing behavior, reshaped)
{
    "type": "forcefield",
    "name": "OPLS-AA",
}

# MLFF
{
    "type": "mlff",
    "name": "MACE-MP-0",
}

# Model chemistry -- MOPAC via MDI
{
    "type": "model_chemistry",
    "model": "MOPAC:SQM@PM6-ORG",   # the citable string
    "program": "MOPAC",
    "method_type": "SQM",
    "method": "PM6-ORG",
    "basis": None,
    "convergence": 1.0,             # mopactools relative-tolerance scale
    "max_iterations": None,
    "options": {},
}

# Model chemistry -- DFT with RI basis and explicit cutoff (VASP-style)
{
    "type": "model_chemistry",
    "model": "VASP:DFT@PBE/PAW@500eV",
    "program": "VASP",
    "method_type": "DFT",
    "method": "PBE",
    "basis": "PAW",
    "cutoff": "500eV",
    "convergence": None,
    "options": {},
}

Charge and multiplicity are deliberately not fields here – per existing SEAMM convention (and the SEAMM paper’s own statement that these are properties of the configuration, not of the calculation method), they are read from the configuration by whatever step launches the actual calculation, exactly as mopac_step/ gaussian_step already do.

The metadata protocol#

Every program-step package that wants to participate in the Model Chemistry Step (as a source of options in its dialog, and/or as an MDI-launchable engine for driver steps like lammps_step) implements one classmethod:

@classmethod
def get_model_chemistry_options(cls, periodic_only=False, mdi_only=False):
    """Return the model chemistries this step's program can provide.

    Parameters
    ----------
    periodic_only : bool
        Only return options valid for a periodic system (relevant
        for any driver running periodic MD, e.g. via MDI/LAMMPS).
    mdi_only : bool
        Only return options actually launchable as an MDI engine,
        as opposed to only available via this step's traditional
        batch/file-based path.

    Returns
    -------
    dict
        Keyed by the bare method name (e.g. "PM6-ORG"). Each value
        is a dict with at least:
            "model_chemistry" : str   -- full "Program:Type@Method" string
            "type" : str              -- "SQM", "DFT", "QC", "FF", "MLFF"
            "description" : str
            "periodic" : bool
            "mdi_capable" : bool
            "mdi_script" : str | None
            "mdi_method_arg" : str | None
    """

This is deliberately a thin classmethod, not a full instance method – the Model Chemistry Step’s dialog (and lammps_step, when deciding how to set up a calculation) need this information before any flowchart node for that program necessarily exists, so it has to be queryable from the class itself via whatever Stevedore already uses to discover the package, without requiring instantiation.

Worked implementation: MOPAC#

mopac_step already has essentially everything needed in its metadata["computational models"] structure – this just exposes it through the protocol above rather than inventing a new shape:

# Proposed addition to mopac_step (exact module/file location to be
# confirmed against the real package -- this is a clean-room sketch
# against the metadata structure retrieved from project knowledge,
# not a verified patch against the actual current file).

# mopactools' model_dict currently supports exactly these six --
# see mopac_mdi.py's --method choices, sourced from
# mopactools.api.MopacSystem.model_dict.
_MDI_CAPABLE_METHODS = {"PM7", "PM6-D3H4", "PM6-ORG", "PM6", "AM1", "RM1"}

@classmethod
def get_model_chemistry_options(cls, periodic_only=False, mdi_only=False):
    options = {}
    for theory_class, class_data in metadata["computational models"].items():
        for family, family_data in class_data["models"].items():
            for name, param in family_data["parameterizations"].items():
                is_periodic = param.get("periodic", False)
                if periodic_only and not is_periodic:
                    continue

                mdi_capable = name in _MDI_CAPABLE_METHODS
                if mdi_only and not mdi_capable:
                    continue

                # PM6-D3H4's periodic=False here is a confirmed real
                # limitation (MOPAC's docs explicitly say D3H4 does
                # not work under PBC) -- this filter correctly
                # excludes it from periodic_only results, not an
                # over-cautious guess. The same is NOT yet confirmed
                # one way or the other for PM6-DH2/DH2X/D3H4X --
                # treat those as unconfirmed for periodic+MDI use
                # until someone actually runs one and checks, per
                # the governing principle in the xTB section below.

                options[name] = {
                    "model_chemistry": f"MOPAC:SQM@{name}",
                    "type": "SQM",
                    "description": param.get("description", ""),
                    "periodic": is_periodic,
                    "elements": param.get("elements", ""),
                    "mdi_capable": mdi_capable,
                    "mdi_script": "mopac_mdi.py" if mdi_capable else None,
                    "mdi_method_arg": name if mdi_capable else None,
                }
    return options

Calling get_model_chemistry_options(periodic_only=True, mdi_only=True) against the real metadata retrieved this week would currently exclude PM6-D3H4 and everything outside the six mopactools-supported names, leaving PM7, PM6-ORG, PM6 as confirmed periodic+MDI options.

PM6-D3H4’s exclusion is now confirmed correct, not just cautious. Per Paul, MOPAC’s documentation explicitly states D3H4 does not work under periodic boundary conditions – this is a real, documented method limitation, not an oversight in mopac_step’s metadata. The filter above is doing exactly the right thing.

The same caution now extends to the other correction variants (PM6-DH2, PM6-DH2X, PM6-D3H4X, etc.) – whether they also have undocumented or differently-documented periodic limitations is not yet known and, per Paul, “cannot be solved now” without actually running them. Per the governing principle stated in the xTB section below, none of these should be treated as periodic+MDI-capable until someone actually runs one periodic via MDI and confirms it works – regardless of what any existing metadata flag says one way or the other. AM1 and RM1’s periodic flags were similarly not visible in the metadata excerpt available while drafting this plan and fall under the same default-to-unconfirmed treatment.

Worked implementation: xTB#

An xtb_step package now exists (added to project knowledge after this plan’s first draft) – this supersedes the earlier speculative sketch. Its real metadata["computational models"] wraps Grimme’s xtb binary via subprocess/file I/O (writes coord.xyz, parses xtbout.json) – a genuinely different code path from tblite_mdi.py, which calls tblite’s Python Calculator directly. Two important consequences fall out of this distinction rather than being assumed:

  • xtb_step’s own metadata flags every method "periodic": False, with an explicit comment that this is because the v1 subprocess wrapper “refuses periodic input.” This is correct for that specific path and does not describe a limitation of the methods themselves – periodic GFN2-xTB was independently validated this week via tblite’s library API directly. The metadata protocol therefore needs two separate periodicity flags, not one: periodic_native (this package’s subprocess path) and periodic_mdi (validated only for the specific methods actually run periodic via MDI).

  • The method sets only partially overlap. xtb_step exposes GFN2-xTB, GFN1-xTB, GFN0-xTB, GFN-FF (via the xtb binary). tblite.interface.Calculator (what tblite_mdi.py wraps) only covers GFN1-xTB, GFN2-xTB, IPEA1-xTB. The overlap is GFN1/GFN2 only – GFN0-xTB and GFN-FF are not reachable via MDI at all through this path; IPEA1-xTB is MDI-reachable but isn’t in xtb_step’s native list.

Governing principle, stated explicitly: a (program, method) combination is only marked periodic-and-MDI-capable if it was actually run periodic via MDI and observed to work – not because some metadata flag says “periodic: True” and not by structural inference from “the underlying library probably supports it.” This week that means GFN2-xTB only – GFN1-xTB and IPEA1-xTB use the identical tblite.interface.Calculator path and are plausibly fine, but were never actually tested periodic, so they default to unconfirmed rather than assumed-working.

# Proposed addition to xtb_step/metadata.py, alongside the existing
# metadata["computational models"] dict (real, from project knowledge).

# tblite.interface.Calculator's scope -- see tblite_mdi.py --method choices.
_MDI_CAPABLE_METHODS = {"GFN1-xTB", "GFN2-xTB", "IPEA1-xTB"}

# Only methods actually run periodic via MDI this week. Everything
# else defaults to unconfirmed regardless of what any other flag says.
_MDI_PERIODIC_VALIDATED = {"GFN2-xTB"}

@classmethod
def get_model_chemistry_options(cls, periodic_only=False, mdi_only=False):
    options = {}
    for theory_class, class_data in metadata["computational models"].items():
        for family, family_data in class_data["models"].items():
            for name, param in family_data["parameterizations"].items():
                # GFN-FF is a force field, not SQM -- Paul's own example
                # of why the method-type tag matters.
                method_type = "FF" if name == "GFN-FF" else "SQM"

                mdi_capable = name in _MDI_CAPABLE_METHODS
                periodic_mdi = name in _MDI_PERIODIC_VALIDATED

                if periodic_only and not periodic_mdi:
                    continue
                if mdi_only and not mdi_capable:
                    continue

                options[name] = {
                    "model_chemistry": f"xTB:{method_type}@{name}",
                    "type": method_type,
                    "description": "",  # pull from param/docs as available
                    "periodic_native": param.get("periodic", False),
                    "periodic_mdi": periodic_mdi,
                    "elements": param.get("elements", ""),
                    "mdi_capable": mdi_capable,
                    "mdi_script": "tblite_mdi.py" if mdi_capable else None,
                    "mdi_method_arg": name if mdi_capable else None,
                }
    return options

Program tag: settled as xTB: , not tblite:. Per Paul: xTB and DFTB+ both use tblite under the hood for GFN-xTB methods, but users recognize the program names, not the underlying library – matching how the rest of this nomenclature already prioritizes user-recognizable program names. (This does mean the same implementation-precision question raised for B3LYP/Gaussian-vs-Psi4 technically applies here too – xtb_step’s binary path and tblite_mdi.py’s library path could in principle disagree slightly even though both report “GFN2-xTB” – but the decision is made deliberately, not by oversight.)

Related, not solved here: dftbplus_step already exists and, per Paul, also uses tblite for GFN-xTB internally alongside its own native DFTB parameterizations. It will eventually want the same two-path (native vs. MDI) treatment in its own metadata. Not pursued in this plan.

Open question, raised but not decided: should tblite_mdi.py (and mopac_mdi.py) actually move into their respective step packages (e.g. xtb_step/mdi_engine.py, mopac_step/mdi_engine.py) now that get_model_chemistry_options() is proposed as a classmethod on the step class? That would put the metadata and the thing it describes in the same package, rather than a standalone script elsewhere referencing the package only by string ("mopac_mdi.py") with no actual import-time connection between them.

How lammps_step consumes this#

Sketch of the decision logic, not yet implemented:

1. Read the global _model_chemistry dict (or a local override, if
   lammps_step's own GUI is later given a "use global" switch).

2. If type == "forcefield": existing behavior, unchanged.

3. If type == "model_chemistry":
     a. Parse "program" from the dict (already split out, no need
        to re-parse the string).
     b. Look up that program's step-package metadata classmethod
        via whatever Stevedore-based registry already locates
        plug-ins, calling get_model_chemistry_options(
            periodic_only=(configuration.periodicity == 3),
            mdi_only=True,
        ) to confirm the requested method is actually valid for
        this configuration and reachable via MDI.
     c. Pull charge/multiplicity from the configuration (existing
        convention) and pass them as CLI args when launching the
        engine subprocess.
     d. Build the "elements" list from the configuration's atom
        types, in LAMMPS type order.
     e. Emit "comm_modify cutoff 2.0" and "fix mdi/qm [virial yes]
        elements ..." instead of pair_style/pair_coeff.
     f. Launch "{mdi_script} -mdi '...' --method {mdi_method_arg}
        --charge {q} --multiplicity {mult}" as the engine process
        (TCP, given this week's validated approach), then launch
        LAMMPS as the driver.

This split deliberately keeps step (3.b) – “resolve a model chemistry to an engine launch spec” – as something any future driver step could call into identically, while steps (3.d)-(3.f) are LAMMPS-specific input-generation logic that stays inside lammps_step. Whether (3.b) should be a shared utility function (in seamm_util or a new small package) or just duplicated logic in each driver step that needs it is an open implementation decision, not resolved here – with only one driver (LAMMPS) existing today, it’s reasonable to write it inside lammps_step first and extract it if/when a second driver actually needs it, rather than build the abstraction speculatively.

Phased plan#

  • Phase 1 (this document): protocol defined; MOPAC metadata classmethod sketched against real existing metadata.

  • Phase 2: confirm the AM1/RM1 periodic flags and the PM6-D3H4 periodic/MDI conflict directly against mopac_step source; land the classmethod for real.

  • Phase 3: new minimal xtb_step package (cookiecutter), hosting the metadata classmethod above plus whatever minimal standalone capability is wanted.

  • Phase 4: the Model Chemistry Step package itself – dialog, _model_chemistry variable, dynamic field exposure driven by the metadata protocol.

  • Phase 5: lammps_step consumption logic (decision tree above) – the actual MDI integration this whole effort is for.

  • Phase 6 (later, optional): migrate gaussian_step/ mopac_step/etc. to have a “use global model chemistry” GUI switch, per the earlier design discussion.

Open questions carried into this plan#

  • PM6-D3H4 periodic-via-mopactools: resolved – confirmed real limitation, correctly excluded.

  • PM6-DH2/DH2X/D3H4X periodic reliability, and AM1/RM1’s periodic flags: still genuinely unconfirmed – per the governing principle, default to excluded from periodic+MDI use until actually tested.

  • GFN1-xTB and IPEA1-xTB periodic reliability via MDI: structurally plausible (same tblite.interface.Calculator as the validated GFN2-xTB path) but not actually tested – same default-to-excluded treatment.

  • xTB: as the program tag: settled – matches user-recognizable program names over the underlying tblite library, consistent with how DFTB+ will eventually need the same treatment.

  • Whether tblite_mdi.py/mopac_mdi.py should move into their respective step packages now that the metadata classmethod is proposed as living on the step class – raised, not decided.

  • Whether (3.b) in the lammps_step consumption section becomes a shared utility now or only once a second driver step exists – deferred deliberately.

  • Whether the Model Chemistry Step replaces or coexists indefinitely with the Forcefield Step’s existing variables during the transition (Paul’s stated intent is “leave Forcefield Step for compatibility for the time being” – the exact deprecation timeline is not fixed here).

References#