systems

systems

Benchmark system configurations.

Each :class:SystemSpec declaratively describes one benchmark system — its molecule source, FF-assembly strategy, and metadata — and the :data:SYSTEMS registry maps a CLI key to that spec. The single public entry point is :func:load_system, which dispatches to the right molecule loader and the right FF strategy based on the spec.

Adding a new system = appending a :class:SystemSpec entry to :data:SYSTEMS; no new module-level function is required.

Usage::

from q2mm.systems import load_system, SYSTEMS

sys_data = load_system("rh-enamide", engine=engine)

The FF-assembly strategies live in :mod:q2mm.models.loaders and correspond to the published-FF workflows in Farrugia, Helquist, Norrby & Wiest 2025 (the QFUERZA paper — see AGENTS.md "Key Papers").

FFStrategy module-attribute

FFStrategy = Literal['qfuerza_fresh', 'published_opt', 'published_opt_composed']

Names of the FF-assembly strategies in :mod:q2mm.models.loaders.

StartingPoint module-attribute

StartingPoint = Literal['published', 'qfuerza']

Choice of starting-parameter values for the optimizer.

• "qfuerza" (default, canonical): take the FF skeleton from the .fld file (atom-type rows, OPT-substructure topology, frozen/ active partition, vdW, stretch-bend, Urey-Bradley) and overwrite the OPT bond and angle scalars with QFUERZA Hessian-derived estimates averaged across the training molecules (Farrugia 2025 protocol). This is the standard QFUERZA workflow as defined in the literature: the chemist provides the FF skeleton (no tool can automate the decisions of which atom types to use or which substructure rows need OPT parameters); QFUERZA fills in the Hessian-derivable scalars. Frozen MM3 backbone parameters are untouched; torsions are zeroed; OPT parameters that :func:qfuerza_into does not estimate (stretch-bends, vdW, Urey-Bradley) retain their literature values — the audit in :func:load_system records this explicitly.
• "published": use the literature OPT values from the .fld file(s) as-is, with no QFUERZA overwrite. This is the publication-baseline path used to reproduce historical convergence runs.

For qfuerza_fresh strategy (CH3F), "qfuerza" is a no-op because the FF is already QFUERZA-derived; the audit records this.

SystemData dataclass

SystemData(molecules: list[Q2MMMolecule], forcefield: ForceField, reference: ReferenceData, qm_freqs_per_mol: list[ndarray], metadata: dict[str, Any] = dict(), normal_modes: dict[str, ndarray] | None = None)

Loaded data for a benchmark system, ready for optimization.

Attributes:

Name	Type	Description
molecules	list[Q2MMMolecule]	One or more molecules with geometry (and optionally Hessians).
forcefield	ForceField	Template force field (from QFUERZA estimation or file).
reference	ReferenceData	Reference data for the objective function.
qm_freqs_per_mol	list[ndarray]	QM real frequencies per molecule (for reporting only). These frequencies are stored separately from reference.
metadata	dict[str, Any]	Extra info (level of theory, molecule name, etc.).
normal_modes	dict[str, ndarray] | None	Pre-computed normal modes for PES distortion analysis. None when not available.

ExternalDataRoots dataclass

ExternalDataRoots(rh_enamide: Path | None = None, supporting_info: Path | None = None, mm3_base: Path | None = None)

Explicit locations for scientific data that Q2MM does not distribute.

Attributes:

Name	Type	Description
rh_enamide	Path | None	Directory containing mm3.fld and the rh_enamide_training_set directory. Environment fallback: Q2MM_RH_ENAMIDE.
supporting_info	Path | None	Root of the extracted Wahlers/Rosales dissertation supporting information. Environment fallback: Q2MM_SUPPORTING_INFO.
mm3_base	Path | None	Path to the licensed MM3 base .fld file. Environment fallback: Q2MM_MM3_BASE.

from_environment classmethod

from_environment() -> ExternalDataRoots

Build roots from Q2MM's documented environment variables.

Source code in q2mm/systems.py

@classmethod
def from_environment(cls) -> ExternalDataRoots:
    """Build roots from Q2MM's documented environment variables."""

    def optional_path(name: str) -> Path | None:
        value = os.environ.get(name)
        return Path(value).expanduser() if value else None

    return cls(
        rh_enamide=optional_path("Q2MM_RH_ENAMIDE"),
        supporting_info=optional_path("Q2MM_SUPPORTING_INFO"),
        mm3_base=optional_path("Q2MM_MM3_BASE"),
    )

SystemSpec dataclass

SystemSpec(key: str, name: str, molecule_loader: Callable[..., list[Q2MMMolecule]], ff_strategy: FFStrategy, uses_external_data_roots: bool = False, ff_paths: Mapping[str, Callable[..., Path]] = dict(), normal_modes_path: Callable[[Path | None], Path | None] | None = None, metadata: Mapping[str, Any] = dict(), metal: str | None = None, description: str = '', default_forms: tuple[str, ...] = ('mm3',))

Declarative spec for one benchmark system.

The :func:load_system dispatcher reads this spec to build a :class:SystemData for the system. Add new systems by appending to :data:SYSTEMS; do not write per-system loader functions.

Attributes:

Name	Type	Description
key	str	CLI key (e.g. "ch3f", "rh-enamide").
name	str	Human-readable system name.
molecule_loader	Callable[..., list[Q2MMMolecule]]	Zero-argument callable returning the training- set molecules (with QM Hessians for the published-FF systems).
uses_external_data_roots	bool	Whether :func:load_system passes the resolved :class:ExternalDataRoots to molecule_loader. Defaults to False so custom zero-argument loaders remain compatible.
ff_strategy	FFStrategy	Name of the FF-assembly strategy in :mod:q2mm.models.loaders. One of "qfuerza_fresh", "published_opt", "published_opt_composed".
ff_paths	Mapping[str, Callable[..., Path]]	Mapping of strategy-specific path keys to callables that return a :class:Path. When uses_external_data_roots is true, they receive :class:ExternalDataRoots; otherwise they remain zero-argument callbacks. Required keys per strategy: qfuerza_fresh: no keys. published_opt: "ff_path" → published .fld. published_opt_composed: "opt_path" → standalone Wahlers OPT-only .fld; "base_path" → MM3 base .fld.
normal_modes_path	Callable[[Path | None], Path | None] | None	Optional callable returning a path to a .npz file with pre-computed normal-mode eigendecomposition (used by PES distortion analysis). Accepts the same data_dir override the CLI may pass via :func:load_system's molecule_loader_kwargs so molecule, Hessian, and normal modes stay co-located. Return None (or a non-existent path) to signal "no modes available".
metadata	Mapping[str, Any]	Static metadata merged into the returned :class:SystemData.metadata (level of theory, publication, etc.).
metal	str | None	Optional element symbol for vdW injection during published_opt_composed (e.g. "PD").
description	str	One-line CLI description.
default_forms	tuple[str, ...]	Functional forms to benchmark by default.

load_rh_enamide_molecules

load_rh_enamide_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]

Load 9 rh-enamide structures with Jaguar Hessians.

This is the shared loader used by both the benchmark CLI and tests.

Returns:

Type	Description
list[Q2MMMolecule]	list[Q2MMMolecule]: 9 molecules with Hessian matrices.

Raises:

Type	Description
FileNotFoundError	If the rh-enamide dataset is not found.
ValueError	If the number of MacroModel structures doesn't match the number of Jaguar input files.

Source code in q2mm/systems.py

def load_rh_enamide_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]:
    """Load 9 rh-enamide structures with Jaguar Hessians.

    This is the shared loader used by both the benchmark CLI and tests.

    Returns:
        list[Q2MMMolecule]: 9 molecules with Hessian matrices.

    Raises:
        FileNotFoundError: If the rh-enamide dataset is not found.
        ValueError: If the number of MacroModel structures doesn't match
            the number of Jaguar input files.

    """
    from q2mm.models.molecule import Q2MMMolecule
    from q2mm.io import JaguarIn, MacroModel

    rh_dir = _resolve_rh_enamide_dir(data_roots)
    training_set_dir = rh_dir / "rh_enamide_training_set"
    mmo_path = training_set_dir / "rh_enamide_training_set.mmo"
    jag_dir = training_set_dir / "jaguar_spe_freq_in_out"
    if not mmo_path.is_file():
        raise FileNotFoundError(f"Rh-enamide MacroModel training set not found: {mmo_path}")
    if not jag_dir.is_dir():
        raise FileNotFoundError(f"Rh-enamide Jaguar training set not found: {jag_dir}")

    mm = MacroModel(str(mmo_path))
    jag_files = sorted(jag_dir.glob("*.in"), key=_natural_sort_key)
    n_structures = len(mm.structures)
    n_jag = len(jag_files)
    if n_structures != n_jag:
        raise ValueError(
            f"Rh-enamide dataset inconsistent: {n_structures} MacroModel structures "
            f"but {n_jag} Jaguar .in files in {jag_dir}"
        )

    molecules = []
    for struct, jag_path in zip(mm.structures, jag_files):
        jag = JaguarIn(str(jag_path))
        hess = jag.get_hessian(len(struct.atoms))
        molecules.append(Q2MMMolecule.from_structure(struct, hessian=hess))
    return molecules

load_heck_relay_molecules

load_heck_relay_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]

Load the 23 Heck relay TS molecules from Gaussian logs.

Returns:

Type	Description
list[Q2MMMolecule]	List of Q2MMMolecule with Hessians and MM3 atom types assigned.

Raises:

Type	Description
FileNotFoundError	If the training set directory is absent.

Source code in q2mm/systems.py

def load_heck_relay_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]:
    """Load the 23 Heck relay TS molecules from Gaussian logs.

    Returns:
        List of Q2MMMolecule with Hessians and MM3 atom types assigned.

    Raises:
        FileNotFoundError: If the training set directory is absent.

    """
    si = _resolve_supporting_info_dir(data_roots)
    ts_dir = si / "rosales" / "Rosales_Anthony_Supporting_Information" / "Chapter3_Heck" / "TrainingSet"
    return _load_gaussian_molecules(ts_dir)

load_pd_allyl_molecules

load_pd_allyl_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]

Load the 21 Pd-allyl amination TS molecules (Wahlers Ch 3).

Source code in q2mm/systems.py

def load_pd_allyl_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]:
    """Load the 21 Pd-allyl amination TS molecules (Wahlers Ch 3)."""
    return _load_wahlers_molecules("Chapter 3", data_roots=data_roots)

load_pd_conjugate_molecules

load_pd_conjugate_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]

Load the 10 Pd 1,4-conjugate addition TS molecules (Wahlers Ch 5).

Source code in q2mm/systems.py

def load_pd_conjugate_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]:
    """Load the 10 Pd 1,4-conjugate addition TS molecules (Wahlers Ch 5)."""
    return _load_wahlers_molecules("Chapter 5", data_roots=data_roots)

load_rh_conjugate_molecules

load_rh_conjugate_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]

Load the 10 Rh 1,4-conjugate addition TS molecules (Wahlers Ch 6).

Source code in q2mm/systems.py

def load_rh_conjugate_molecules(*, data_roots: ExternalDataRoots | None = None) -> list[Q2MMMolecule]:
    """Load the 10 Rh 1,4-conjugate addition TS molecules (Wahlers Ch 6)."""
    return _load_wahlers_molecules("Chapter 6", data_roots=data_roots)

load_system

load_system(key: str, *, engine: Any | None = None, functional_form: str | None = None, molecule_loader_kwargs: dict[str, Any] | None = None, data_roots: ExternalDataRoots | None = None, starting_point: StartingPoint = 'qfuerza', qfuerza_replace_with: float = 1.0) -> SystemData

Build a :class:SystemData for one benchmark system.

The single loader entry point. Dispatches on the system's :class:SystemSpec to build the molecule list and the force field, then constructs the reference data and metadata.

Reference construction depends on the FF strategy:

• qfuerza_fresh (single-molecule, e.g. CH3F): a frequency-only reference is built from engine.frequencies(molecule, ff) compared against the QM frequencies derived from the Hessian. The engine must be provided in this case.
• published_opt / published_opt_composed (multi-molecule published-FF systems): an eigenmatrix-diagonal reference is built via :meth:ReferenceData.from_molecules. The engine is unused.

Parameters:

Name	Type	Description	Default
key	str	System key from :data:SYSTEMS.	required
engine	Any | None	MM engine instance; required for qfuerza_fresh systems (CH3F-style). Unused for other strategies.	None
functional_form	str | None	Override the FF's functional form ("harmonic" or "mm3").	None
molecule_loader_kwargs	dict[str, Any] | None	Optional kwargs forwarded to the system's molecule loader (e.g. data_dir for CH3F).	None
data_roots	ExternalDataRoots | None	Explicit locations for non-distributed scientific data. When omitted, the named Q2MM_* environment variables described by :class:ExternalDataRoots are used. Built-in CH3F/SN2 resources do not require external roots.	None
starting_point	StartingPoint	Where the starting OPT parameter values come from. See :data:StartingPoint. "qfuerza" is the canonical default: it overwrites OPT bond/angle values with multi-molecule QFUERZA estimates after FF assembly (per Farrugia 2025). "published" retains the literature OPT values verbatim — pass this to reproduce the historical publication-baseline runs. CH3F (qfuerza_fresh strategy) treats "qfuerza" as a no-op since the FF is already QFUERZA-derived.	'qfuerza'
qfuerza_replace_with	float	Replacement value (Hartree/Bohr²) for the most negative TS-Hessian eigenvalue during QFUERZA projection. Default 1.0 matches Limé & Norrby Method C and preserves historical numbers. Smaller values (e.g. 0.03, Method D's "natural" value) reduce the chance of negative angle force constants in the starting FF but produce a softer reaction-coordinate mode. Applied whenever QFUERZA invokes invert_ts_curvature — i.e. for the qfuerza_fresh strategy regardless of starting_point, and for the published-FF strategies only when starting_point="qfuerza" triggers a QFUERZA overwrite. Has no effect for published-FF strategies with starting_point="published" (no QFUERZA invocation occurs) or for any ground-state Hessian (no negative eigenvalue to replace).	1.0

Returns:

Type	Description
SystemData	A fully-populated :class:SystemData. The
SystemData	metadata["starting_point_audit"] block reports, for each
SystemData	parameter type, how many active scalars were QFUERZA-overwritten
SystemData	vs retained from the published FF vs frozen — see
SystemData	func:_audit_starting_point.

Raises:

Type	Description
KeyError	If key is not in :data:SYSTEMS.
TypeError	If engine is required but not provided.
ValueError	If starting_point is not one of "published" or "qfuerza".

Source code in q2mm/systems.py

def load_system(
    key: str,
    *,
    engine: Any | None = None,
    functional_form: str | None = None,
    molecule_loader_kwargs: dict[str, Any] | None = None,
    data_roots: ExternalDataRoots | None = None,
    starting_point: StartingPoint = "qfuerza",
    qfuerza_replace_with: float = 1.0,
) -> SystemData:
    """Build a :class:`SystemData` for one benchmark system.

    The single loader entry point.  Dispatches on the system's
    :class:`SystemSpec` to build the molecule list and the force
    field, then constructs the reference data and metadata.

    Reference construction depends on the FF strategy:

    - ``qfuerza_fresh`` (single-molecule, e.g. CH3F): a frequency-only
      reference is built from ``engine.frequencies(molecule, ff)``
      compared against the QM frequencies derived from the Hessian.
      The *engine* must be provided in this case.
    - ``published_opt`` / ``published_opt_composed`` (multi-molecule
      published-FF systems): an eigenmatrix-diagonal reference is built
      via :meth:`ReferenceData.from_molecules`.  The engine is unused.

    Args:
        key: System key from :data:`SYSTEMS`.
        engine: MM engine instance; required for ``qfuerza_fresh``
            systems (CH3F-style).  Unused for other strategies.
        functional_form: Override the FF's functional form
            (``"harmonic"`` or ``"mm3"``).
        molecule_loader_kwargs: Optional kwargs forwarded to the
            system's molecule loader (e.g. ``data_dir`` for CH3F).
        data_roots: Explicit locations for non-distributed scientific data.
            When omitted, the named ``Q2MM_*`` environment variables described
            by :class:`ExternalDataRoots` are used. Built-in CH3F/SN2 resources
            do not require external roots.
        starting_point: Where the starting OPT parameter values come
            from.  See :data:`StartingPoint`.  ``"qfuerza"`` is the
            canonical default: it overwrites OPT bond/angle values with
            multi-molecule QFUERZA estimates after FF assembly (per
            Farrugia 2025).  ``"published"`` retains the literature OPT
            values verbatim — pass this to reproduce the historical
            publication-baseline runs.  CH3F (``qfuerza_fresh``
            strategy) treats ``"qfuerza"`` as a no-op since the FF is
            already QFUERZA-derived.
        qfuerza_replace_with: Replacement value (Hartree/Bohr²) for the
            most negative TS-Hessian eigenvalue during QFUERZA
            projection.  Default ``1.0`` matches Limé & Norrby Method C
            and preserves historical numbers.  Smaller values (e.g.
            ``0.03``, Method D's "natural" value) reduce the chance of
            negative angle force constants in the starting FF but
            produce a softer reaction-coordinate mode.  Applied
            whenever QFUERZA invokes ``invert_ts_curvature`` — i.e.
            for the ``qfuerza_fresh`` strategy regardless of
            ``starting_point``, and for the published-FF strategies
            only when ``starting_point="qfuerza"`` triggers a QFUERZA
            overwrite.  Has no effect for published-FF strategies with
            ``starting_point="published"`` (no QFUERZA invocation
            occurs) or for any ground-state Hessian (no negative
            eigenvalue to replace).

    Returns:
        A fully-populated :class:`SystemData`.  The
        ``metadata["starting_point_audit"]`` block reports, for each
        parameter type, how many active scalars were QFUERZA-overwritten
        vs retained from the published FF vs frozen — see
        :func:`_audit_starting_point`.

    Raises:
        KeyError: If *key* is not in :data:`SYSTEMS`.
        TypeError: If *engine* is required but not provided.
        ValueError: If *starting_point* is not one of ``"published"`` or ``"qfuerza"``.

    """
    from q2mm.models import loaders
    from q2mm.models.forcefield import FunctionalForm
    from q2mm.optimizers.objective import ReferenceData

    if key not in SYSTEMS:
        raise KeyError(f"Unknown system {key!r}; available: {sorted(SYSTEMS)}")
    if starting_point not in ("published", "qfuerza"):
        raise ValueError(f"Unknown starting_point {starting_point!r}; must be one of: 'published', 'qfuerza'")
    spec = SYSTEMS[key]
    resolved_data_roots = _external_roots(data_roots)

    # 1. Molecules ---------------------------------------------------------
    loader_kwargs = dict(molecule_loader_kwargs or {})
    if spec.uses_external_data_roots:
        loader_kwargs["data_roots"] = resolved_data_roots
    molecules = spec.molecule_loader(**loader_kwargs)

    # 2. Force field via the named strategy --------------------------------
    def ff_path(name: str) -> Path:
        callback = spec.ff_paths[name]
        if spec.uses_external_data_roots:
            return callback(resolved_data_roots)
        return callback()

    if spec.ff_strategy == "qfuerza_fresh":
        if len(molecules) != 1:
            raise ValueError(
                f"qfuerza_fresh strategy requires exactly 1 molecule, got {len(molecules)} for system {key!r}"
            )
        ff = loaders.load_qfuerza_fresh(molecules[0], replace_with=qfuerza_replace_with)
    elif spec.ff_strategy == "published_opt":
        ff = loaders.load_published_opt(ff_path("ff_path"))
    elif spec.ff_strategy == "published_opt_composed":
        ff = loaders.load_published_opt_composed(
            ff_path("opt_path"),
            ff_path("base_path"),
            metal=spec.metal,
        )
    else:  # pragma: no cover — unreachable per FFStrategy literal
        raise ValueError(f"Unknown ff_strategy {spec.ff_strategy!r} for system {key!r}")

    # Functional form: explicit override > strategy default (MM3)
    if functional_form is not None:
        ff.functional_form = FunctionalForm(functional_form)

    # 2b. Optional QFUERZA overwrite of OPT parameter values --------------
    # For ``starting_point="qfuerza"`` we keep the published FF
    # topology + frozen/active partition but replace the active
    # OPT bond/angle values with multi-molecule QFUERZA estimates
    # (Farrugia 2025 §"Methods", per-molecule mean averaging).
    # ``qfuerza_into`` honours the frozen partition, so MM3 backbone
    # rows are never touched.  Active rows that QFUERZA cannot match
    # to any training molecule retain their published values — the
    # audit captures this so the leak is visible downstream.
    from q2mm.models.seminario import qfuerza_into

    before_vec: np.ndarray | None = None
    if starting_point == "qfuerza" and spec.ff_strategy != "qfuerza_fresh":
        before_vec = ff.get_param_vector().copy()
        qfuerza_into(
            ff,
            molecules,
            invert_ts_curvature=True,
            replace_with=qfuerza_replace_with,
        )
    starting_point_audit = _audit_starting_point(ff, before_vec=before_vec, starting_point=starting_point)

    # 3. Reference data ----------------------------------------------------
    if spec.ff_strategy == "qfuerza_fresh":
        # Frequency-only reference: requires the engine to compute MM
        # frequencies against the QM ones derived from the Hessian.
        if engine is None:
            raise TypeError(
                f"System {key!r} uses ff_strategy='qfuerza_fresh' which requires "
                "engine= to build the frequency-only reference."
            )
        molecule = molecules[0]
        qm_freqs_all = _qm_frequencies_from_hessian(molecule.hessian, molecule.symbols)
        mm_all = engine.frequencies(molecule, ff)
        reference, qm_real = _build_frequency_reference(qm_freqs_all, mm_all)
        qm_freqs_per_mol = [qm_real]
    else:
        reference = ReferenceData.from_molecules(molecules, eigenmatrix_diagonal_only=False)
        qm_freqs_per_mol = []
        for mol in molecules:
            qm_freqs = _qm_frequencies_from_hessian(mol.hessian, mol.symbols)
            qm_real = np.array(sorted(f for f in qm_freqs if f > 50.0))
            qm_freqs_per_mol.append(qm_real)

    # 4. Optional normal-modes side-load (declared per-spec).  Pass the
    #    CLI's data_dir override (if any) through so molecule, Hessian,
    #    and normal modes stay co-located.
    normal_modes: dict[str, np.ndarray] | None = None
    if spec.normal_modes_path is not None:
        data_dir_override = (molecule_loader_kwargs or {}).get("data_dir")
        modes_path = spec.normal_modes_path(data_dir_override)
        if modes_path is not None and modes_path.exists():
            from q2mm.diagnostics.pes_distortion import load_normal_modes

            normal_modes = load_normal_modes(modes_path)

    # 5. Wrap in SystemData ------------------------------------------------
    resolved_form = functional_form
    if resolved_form is None and ff.functional_form is not None:
        resolved_form = ff.functional_form.value
    metadata = {
        "molecule_name": spec.name,
        "n_molecules": len(molecules),
        "n_atoms_per_mol": [len(m.symbols) for m in molecules],
        "starting_point": starting_point,
        "starting_point_audit": starting_point_audit,
        **dict(spec.metadata),
        **({"functional_form": resolved_form} if resolved_form else {}),
    }
    return SystemData(
        molecules=molecules,
        forcefield=ff,
        reference=reference,
        qm_freqs_per_mol=qm_freqs_per_mol,
        metadata=metadata,
        normal_modes=normal_modes,
    )