openmm

openmm

OpenMM molecular mechanics backend.

Provides a full-featured MM engine using OpenMM for energy, minimization, Hessian, and frequency calculations. Supports both harmonic and MM3 functional forms with runtime parameter updates via :class:OpenMMHandle.

OpenMMHandle dataclass

OpenMMHandle(molecule: Q2MMMolecule, system: object, integrator: object, context: object, bond_force: object | None, angle_force: object | None, torsion_force: object | None, vdw_force: object | None, ub_force: object | None, cmap_force: object | None, bond_terms: list[_BondTerm], angle_terms: list[_AngleTerm], torsion_terms: list[_TorsionTerm], vdw_terms: list[_VdwTerm], ub_terms: list[_UBTerm] = list(), cmap_terms: list[_CmapTerm] = list(), exceptions_14: list[_Exception14] = list(), functional_form: FunctionalForm = MM3)

Reusable OpenMM system/context pair for fast parameter updates.

Attributes:

Name	Type	Description
molecule	Q2MMMolecule	Deep copy of the input molecule.
system	object	The openmm.System object.
integrator	object	The openmm.Integrator used by the context.
context	object	The openmm.Context for energy evaluation.
bond_force	object | None	The OpenMM bond force object, or None if no bonds.
angle_force	object | None	The OpenMM angle force object, or None if no angles.
torsion_force	object | None	The OpenMM torsion force object, or None if no torsions.
vdw_force	object | None	The OpenMM vdW force object, or None if no vdW terms.
ub_force	object | None	The OpenMM HarmonicBondForce for Urey-Bradley terms, or None.
cmap_force	object | None	The OpenMM CMAPTorsionForce, or None if no CMAP terms.
bond_terms	list[_BondTerm]	Mapping of molecule bonds to force indices.
angle_terms	list[_AngleTerm]	Mapping of molecule angles to force indices.
torsion_terms	list[_TorsionTerm]	Mapping of molecule torsions to force indices.
vdw_terms	list[_VdwTerm]	Mapping of atoms to vdW particle indices.
ub_terms	list[_UBTerm]	Mapping of Urey-Bradley 1-3 pairs to force indices.
cmap_terms	list[_CmapTerm]	Mapping of CMAP corrections to force indices.
exceptions_14	list[_Exception14]	1-4 nonbonded exceptions (harmonic form only).
functional_form	FunctionalForm	The functional form used when the handle was created.

OpenMMEngine

OpenMMEngine(platform_name: str | None = None, precision: str | None = None)

Bases: MMEngine

Molecular mechanics backend powered by OpenMM.

Supports both harmonic (AMBER-style) and MM3 functional forms. Provides reusable :class:OpenMMHandle objects for fast parameter updates during optimization loops.

Initialize the OpenMM engine.

Parameters:

Name	Type	Description	Default
platform_name	str | None	OpenMM platform to use (e.g. "CPU", "CUDA", "OpenCL"). When None, the fastest available platform is auto-detected via :func:detect_best_platform (CUDA > OpenCL > CPU > Reference). WSL2 is recommended for CUDA on modern GPUs (e.g. RTX 5090 Blackwell) when running on Windows hardware.	None
precision	str | None	Floating-point precision for GPU platforms ("single", "mixed", or "double"). Ignored for CPU/Reference platforms. Defaults to "mixed" when a GPU platform is selected.	None

Raises:

Type	Description
ImportError	If OpenMM is not installed.

Source code in q2mm/backends/mm/openmm.py

def __init__(
    self,
    platform_name: str | None = None,
    precision: str | None = None,
) -> None:
    """Initialize the OpenMM engine.

    Args:
        platform_name: OpenMM platform to use (e.g. ``"CPU"``,
            ``"CUDA"``, ``"OpenCL"``).  When ``None``, the fastest
            available platform is auto-detected via
            :func:`detect_best_platform` (CUDA > OpenCL > CPU >
            Reference).  WSL2 is recommended for CUDA on modern
            GPUs (e.g. RTX 5090 Blackwell) when running on Windows
            hardware.
        precision: Floating-point precision for GPU platforms
            (``"single"``, ``"mixed"``, or ``"double"``).  Ignored
            for CPU/Reference platforms.  Defaults to ``"mixed"``
            when a GPU platform is selected.

    Raises:
        ImportError: If OpenMM is not installed.

    """
    _ensure_openmm()
    if platform_name is None:
        platform_name = detect_best_platform()
    self._platform_name = platform_name

    _VALID_PRECISIONS = {"single", "mixed", "double"}
    if precision is not None:
        precision = precision.strip().lower()
        if precision not in _VALID_PRECISIONS:
            raise ValueError(
                f"Invalid precision {precision!r}. Allowed values: {', '.join(sorted(_VALID_PRECISIONS))}."
            )
    self._precision = precision
    logger.info("OpenMM platform: %s", self._platform_name)

name property

name: str

Human-readable engine name including the active platform.

Returns:

Name	Type	Description
str	str	e.g. "OpenMM (CUDA)".

supported_functional_forms

supported_functional_forms() -> frozenset[str]

Functional forms this engine can evaluate.

Returns:

Type	Description
frozenset[str]	frozenset[str]: {"harmonic", "mm3"}.

Source code in q2mm/backends/mm/openmm.py

def supported_functional_forms(self) -> frozenset[str]:
    """Functional forms this engine can evaluate.

    Returns:
        frozenset[str]: ``{"harmonic", "mm3"}``.

    """
    return frozenset({"harmonic", "mm3"})

is_available

is_available() -> bool

Check if OpenMM is installed.

Returns:

Name	Type	Description
bool	bool	True if the openmm package is importable.

Source code in q2mm/backends/mm/openmm.py

def is_available(self) -> bool:
    """Check if OpenMM is installed.

    Returns:
        bool: ``True`` if the ``openmm`` package is importable.

    """
    return _HAS_OPENMM

supports_runtime_params

supports_runtime_params() -> bool

Whether parameters can be updated without rebuilding the system.

Returns:

Name	Type	Description
bool	bool	Always True for OpenMM.

Source code in q2mm/backends/mm/openmm.py

def supports_runtime_params(self) -> bool:
    """Whether parameters can be updated without rebuilding the system.

    Returns:
        bool: Always ``True`` for OpenMM.

    """
    return True

supports_analytical_gradients

supports_analytical_gradients() -> bool

Whether this engine provides analytical parameter gradients.

Both HARMONIC and MM3 functional forms use CustomBondForce, CustomAngleForce, and CustomTorsionForce with global parameters, so addEnergyParameterDerivative() provides exact dE/d(param) for bond, angle, and torsion parameters.

vdW parameters use per-particle values and are supplemented with central finite differences inside energy_and_param_grad, so the returned gradient is always complete.

Returns:

Name	Type	Description
bool	bool	Always True.

Source code in q2mm/backends/mm/openmm.py

def supports_analytical_gradients(self) -> bool:
    """Whether this engine provides analytical parameter gradients.

    Both HARMONIC and MM3 functional forms use ``CustomBondForce``,
    ``CustomAngleForce``, and ``CustomTorsionForce`` with global
    parameters, so ``addEnergyParameterDerivative()`` provides exact
    dE/d(param) for bond, angle, and torsion parameters.

    vdW parameters use per-particle values and are supplemented
    with central finite differences inside ``energy_and_param_grad``,
    so the returned gradient is always complete.

    Returns:
        bool: Always ``True``.

    """
    return True

create_context

create_context(structure: Q2MMMolecule | str | Path, forcefield: ForceField | None = None, *, precision: str | None = None) -> OpenMMHandle

Build an OpenMM system and context for a molecule.

Creates force objects (bond, angle, vdW) matching the force field's functional form and assigns per-term parameters from forcefield.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule | str | Path	A :class:~q2mm.models.molecule.Q2MMMolecule or path to an XYZ file.	required
forcefield	ForceField | None	Force field to apply. Auto-generated from the molecule if None.	None
precision	str | None	Override GPU precision ("single", "mixed", "double"). When None (default) uses the engine-level setting.	None

Returns:

Name	Type	Description
OpenMMHandle	OpenMMHandle	Reusable handle for energy evaluation and parameter updates.

Raises:

Type	Description
KeyError	If an atom's element has no defined mass.
ValueError	If no OpenMM terms could be created from the force field, or if a vdW parameter is missing for an atom.

Source code in q2mm/backends/mm/openmm.py

def create_context(
    self,
    structure: Q2MMMolecule | str | Path,
    forcefield: ForceField | None = None,
    *,
    precision: str | None = None,
) -> OpenMMHandle:
    """Build an OpenMM system and context for a molecule.

    Creates force objects (bond, angle, vdW) matching the force field's
    functional form and assigns per-term parameters from *forcefield*.

    Args:
        structure (Q2MMMolecule | str | Path): A
            :class:`~q2mm.models.molecule.Q2MMMolecule` or path to an
            XYZ file.
        forcefield: Force field to apply. Auto-generated from the
            molecule if ``None``.
        precision: Override GPU precision (``"single"``, ``"mixed"``,
            ``"double"``).  When ``None`` (default) uses the
            engine-level setting.

    Returns:
        OpenMMHandle: Reusable handle for energy evaluation and parameter
            updates.

    Raises:
        KeyError: If an atom's element has no defined mass.
        ValueError: If no OpenMM terms could be created from the force
            field, or if a vdW parameter is missing for an atom.

    """
    molecule = _as_molecule(structure)
    forcefield = _coerce_forcefield(forcefield, molecule)
    self._validate_forcefield(forcefield)

    # Default to MM3 for backward compatibility when functional_form is None
    ff_form = forcefield.functional_form or FunctionalForm.MM3
    use_harmonic = ff_form == FunctionalForm.HARMONIC

    system = mm.System()
    for symbol in molecule.symbols:
        if symbol not in MASSES:
            raise KeyError(f"No atomic mass is defined for element '{symbol}'.")
        system.addParticle(MASSES[symbol] * unit.dalton)

    # --- Create force objects based on functional form ---
    if use_harmonic:
        bond_force = mm.HarmonicBondForce()
    else:
        bond_force = mm.CustomBondForce(
            f"k*(10*(r-r0))^2*(1-c3*(10*(r-r0))+c4*(10*(r-r0))^2);c3={MM3_BOND_C3};c4={MM3_BOND_C4}"
        )
        bond_force.addPerBondParameter("k")
        bond_force.addPerBondParameter("r0")

    if use_harmonic:
        angle_force = mm.HarmonicAngleForce()
    else:
        angle_force = mm.CustomAngleForce(
            "k*(theta-theta0)^2*("
            "1+a3*((theta-theta0)*deg)"
            "+a4*((theta-theta0)*deg)^2"
            "+a5*((theta-theta0)*deg)^3"
            "+a6*((theta-theta0)*deg)^4"
            ");"
            f"a3={MM3_ANGLE_C3};"
            f"a4={MM3_ANGLE_C4};"
            f"a5={MM3_ANGLE_C5};"
            f"a6={MM3_ANGLE_C6};"
            f"deg={RAD_TO_DEG}"
        )
        angle_force.addPerAngleParameter("k")
        angle_force.addPerAngleParameter("theta0")

    # Torsion force: both harmonic (AMBER) and MM3 use PeriodicTorsionForce
    # E = k * (1 + cos(n*θ − phase))
    torsion_force = mm.PeriodicTorsionForce()

    if use_harmonic:
        vdw_force = mm.NonbondedForce()
        vdw_force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
    else:
        # MM3 Buckingham exp-6 with short-range repulsive wall.
        # Below r < 0.34·rv the attractive r^-6 dominates the exponential,
        # causing divergence to -∞.  Switch to a hard repulsive form
        # at that boundary using step() to prevent collapse.
        vdw_force = mm.CustomNonbondedForce(
            f"step(r - rc) * epsilon*(-{MM3_VDW_C}*(rv/r)^6 + {MM3_VDW_A}*exp(-{MM3_VDW_B}*r/rv))"
            f" + step(rc - r) * epsilon*{MM3_VDW_A}*exp(-{MM3_VDW_B}*rc/rv) * (rc/r)^12;"
            "rc=0.34*rv;"
            "rv=radius1+radius2;"
            "epsilon=sqrt(epsilon1*epsilon2)"
        )
        vdw_force.addPerParticleParameter("radius")
        vdw_force.addPerParticleParameter("epsilon")
        vdw_force.setNonbondedMethod(mm.CustomNonbondedForce.NoCutoff)

    # --- Assign bond parameters ---
    bond_terms: list[_BondTerm] = []
    for bond in molecule.bonds:
        param = _match_bond(forcefield, bond.elements, env_id=bond.env_id, ff_row=bond.ff_row)
        if param is None:
            continue
        if use_harmonic:
            force_index = bond_force.addBond(
                bond.atom_i,
                bond.atom_j,
                ang_to_nm(param.equilibrium),
                _bond_k_to_harmonic(param.force_constant),
            )
        else:
            force_index = bond_force.addBond(
                bond.atom_i,
                bond.atom_j,
                [_bond_k_to_openmm(param.force_constant), ang_to_nm(param.equilibrium)],
            )
        bond_terms.append(
            _BondTerm(
                force_index=force_index,
                atom_i=bond.atom_i,
                atom_j=bond.atom_j,
                elements=bond.elements,
                env_id=bond.env_id,
                ff_row=bond.ff_row,
            )
        )

    # --- Assign angle parameters ---
    angle_terms: list[_AngleTerm] = []
    for angle in molecule.angles:
        param = _match_angle(forcefield, angle.elements, env_id=angle.env_id, ff_row=angle.ff_row)
        if param is None:
            continue
        if use_harmonic:
            force_index = angle_force.addAngle(
                angle.atom_i,
                angle.atom_j,
                angle.atom_k,
                np.deg2rad(float(param.equilibrium)),
                _angle_k_to_harmonic(param.force_constant),
            )
        else:
            force_index = angle_force.addAngle(
                angle.atom_i,
                angle.atom_j,
                angle.atom_k,
                [_angle_k_to_openmm(param.force_constant), np.deg2rad(float(param.equilibrium))],
            )
        angle_terms.append(
            _AngleTerm(
                force_index=force_index,
                atom_i=angle.atom_i,
                atom_j=angle.atom_j,
                atom_k=angle.atom_k,
                elements=angle.elements,
                env_id=angle.env_id,
                ff_row=angle.ff_row,
            )
        )

    # --- Assign Urey-Bradley terms (1-3 distance for CHARMM angles) ---
    ub_force = None
    ub_terms: list[_UBTerm] = []
    for angle_term in angle_terms:
        param = _match_angle(forcefield, angle_term.elements, env_id=angle_term.env_id, ff_row=angle_term.ff_row)
        if param is None:
            continue
        # Require both UB parameters to be provided together. If neither is
        # set, there is simply no UB term for this angle; if only one is set,
        # this is a configuration error that should fail fast.
        if param.ub_force_constant is None and param.ub_equilibrium is None:
            continue
        if param.ub_force_constant is None or param.ub_equilibrium is None:
            raise ValueError(
                "Inconsistent Urey-Bradley parameters for angle "
                f"{angle_term.elements} (env_id={angle_term.env_id}, ff_row={angle_term.ff_row}): "
                "both 'ub_force_constant' and 'ub_equilibrium' must be set or both must be None."
            )
        if ub_force is None:
            ub_force = mm.HarmonicBondForce()
        # UB uses atom_i and atom_k (the two outer atoms of the angle)
        # Convert: k in kcal/(mol·Å²) → kJ/(mol·nm²) via _bond_k_to_harmonic
        # eq in Å → nm via ang_to_nm
        force_index = ub_force.addBond(
            angle_term.atom_i,
            angle_term.atom_k,
            ang_to_nm(param.ub_equilibrium),
            _bond_k_to_harmonic(param.ub_force_constant),
        )
        ub_terms.append(
            _UBTerm(
                force_index=force_index,
                atom_i=angle_term.atom_i,
                atom_k=angle_term.atom_k,
                elements=angle_term.elements,
                env_id=angle_term.env_id,
                ff_row=angle_term.ff_row,
            )
        )

    # --- Assign proper torsion parameters ---
    torsion_terms: list[_TorsionTerm] = []
    for torsion in molecule.torsions:
        params = _match_torsions(
            forcefield,
            torsion.element_quad,
            env_id=torsion.env_id,
            ff_row=torsion.ff_row,
            is_improper=False,
        )
        for param in params:
            force_index = torsion_force.addTorsion(
                torsion.atom_i,
                torsion.atom_j,
                torsion.atom_k,
                torsion.atom_l,
                param.periodicity,
                np.deg2rad(float(param.phase)),
                _torsion_k_to_openmm(param.force_constant),
            )
            torsion_terms.append(
                _TorsionTerm(
                    force_index=force_index,
                    atom_i=torsion.atom_i,
                    atom_j=torsion.atom_j,
                    atom_k=torsion.atom_k,
                    atom_l=torsion.atom_l,
                    elements=torsion.element_quad,
                    periodicity=param.periodicity,
                    env_id=torsion.env_id,
                    ff_row=param.ff_row,
                )
            )

    # --- Assign improper torsion parameters ---
    # Improper torsions are detected from trigonal centres in the
    # molecule topology and matched against FF params with
    # is_improper=True.
    for imp_torsion in molecule.improper_torsions:
        params = _match_torsions(
            forcefield,
            imp_torsion.element_quad,
            env_id=imp_torsion.env_id,
            ff_row=imp_torsion.ff_row,
            is_improper=True,
        )
        for param in params:
            force_index = torsion_force.addTorsion(
                imp_torsion.atom_i,
                imp_torsion.atom_j,
                imp_torsion.atom_k,
                imp_torsion.atom_l,
                param.periodicity,
                np.deg2rad(float(param.phase)),
                _torsion_k_to_openmm(param.force_constant),
            )
            torsion_terms.append(
                _TorsionTerm(
                    force_index=force_index,
                    atom_i=imp_torsion.atom_i,
                    atom_j=imp_torsion.atom_j,
                    atom_k=imp_torsion.atom_k,
                    atom_l=imp_torsion.atom_l,
                    elements=imp_torsion.element_quad,
                    periodicity=param.periodicity,
                    env_id=imp_torsion.env_id,
                    ff_row=param.ff_row,
                    is_improper=True,
                )
            )

    # --- Assign vdW parameters ---
    vdw_terms: list[_VdwTerm] = []
    if forcefield.vdws:
        for atom_index, (symbol, atom_type) in enumerate(zip(molecule.symbols, molecule.atom_types, strict=False)):
            param = _match_vdw(forcefield, atom_type=atom_type, element=symbol)
            if param is None:
                raise ValueError(f"Missing vdW parameter for atom {atom_index + 1} ({atom_type or symbol}).")
            if use_harmonic:
                vdw_force.addParticle(0.0, _vdw_sigma_nm(param.radius), _vdw_epsilon_to_openmm(param.epsilon))
            else:
                vdw_force.addParticle([_vdw_radius_to_openmm(param.radius), _vdw_epsilon_to_openmm(param.epsilon)])
            vdw_terms.append(
                _VdwTerm(
                    particle_index=atom_index,
                    atom_type=atom_type,
                    element=symbol,
                    ff_row=param.ff_row,
                )
            )
        if use_harmonic:
            # Exclude 1-2 and 1-3 pairs; scale 1-4 pairs by 1/2 (AMBER scnb=2.0)
            excluded_12: set[tuple[int, int]] = set()
            for bond in molecule.bonds:
                excluded_12.add((min(bond.atom_i, bond.atom_j), max(bond.atom_i, bond.atom_j)))

            excluded_13: set[tuple[int, int]] = set()
            for angle in molecule.angles:
                excluded_13.add((min(angle.atom_i, angle.atom_k), max(angle.atom_i, angle.atom_k)))
            excluded_13 -= excluded_12  # pure 1-3 only

            # Build adjacency for 1-4 detection
            neighbors: dict[int, set[int]] = {}
            for bond in molecule.bonds:
                neighbors.setdefault(bond.atom_i, set()).add(bond.atom_j)
                neighbors.setdefault(bond.atom_j, set()).add(bond.atom_i)

            pairs_14: set[tuple[int, int]] = set()
            for angle in molecule.angles:
                # 1-4 pairs: atoms bonded to angle endpoints but not in the angle
                for nb in neighbors.get(angle.atom_i, ()):
                    if nb != angle.atom_j and nb != angle.atom_k:
                        pairs_14.add((min(nb, angle.atom_k), max(nb, angle.atom_k)))
                for nb in neighbors.get(angle.atom_k, ()):
                    if nb != angle.atom_j and nb != angle.atom_i:
                        pairs_14.add((min(nb, angle.atom_i), max(nb, angle.atom_i)))
            pairs_14 -= excluded_12
            pairs_14 -= excluded_13

            # Fully exclude 1-2 and 1-3
            for p1, p2 in sorted(excluded_12 | excluded_13):
                vdw_force.addException(p1, p2, 0.0, 1.0, 0.0)

            # Scale 1-4 pairs (AMBER: scnb=2.0 → epsilon/2, scee=1.2 → no charges here)
            SCNB = 2.0
            exceptions_14: list[_Exception14] = []
            for p1, p2 in sorted(pairs_14):
                sig1, eps1 = vdw_force.getParticleParameters(p1)[1:]
                sig2, eps2 = vdw_force.getParticleParameters(p2)[1:]
                sig_14 = 0.5 * (sig1 + sig2)
                eps_14 = (eps1 * eps2) ** 0.5 / SCNB
                exc_idx = vdw_force.addException(p1, p2, 0.0, sig_14, eps_14)
                exceptions_14.append(_Exception14(exception_index=exc_idx, particle_i=p1, particle_j=p2))
        else:
            vdw_force.createExclusionsFromBonds([(bond.atom_i, bond.atom_j) for bond in molecule.bonds], 2)

    # --- Assign CMAP correction terms (CHARMM backbone corrections) ---
    cmap_force = None
    cmap_terms: list[_CmapTerm] = []
    if forcefield.has_cmap:
        cmap_force = mm.CMAPTorsionForce()
        # Build atom-type-to-index mapping for the molecule
        type_to_indices = _build_atom_type_index(molecule)

        for grid in forcefield.cmaps:
            # Add the 2D energy grid (convert kcal/mol → kJ/mol for OpenMM)
            energy_kj = [e * KCAL_TO_KJ for e in grid.energy]
            map_index = cmap_force.addMap(grid.resolution, energy_kj)

            # Find atom index quadruples matching the phi/psi types
            phi_matches = _find_dihedral_atoms(type_to_indices, grid.atom_types_phi, molecule)
            psi_matches = _find_dihedral_atoms(type_to_indices, grid.atom_types_psi, molecule)

            # Pair phi/psi dihedrals sharing 3 overlapping atoms.
            # Index psi by first 3 atoms for O(n_phi + n_psi) lookup.
            psi_index: dict[tuple[int, ...], list[tuple[int, int, int, int]]] = {}
            for psi_atoms in psi_matches:
                key = tuple(psi_atoms[:3])
                psi_index.setdefault(key, []).append(psi_atoms)

            for phi_atoms in phi_matches:
                key = tuple(phi_atoms[1:])
                for psi_atoms in psi_index.get(key, ()):
                    torsion_index = cmap_force.addTorsion(
                        map_index,
                        phi_atoms[0],
                        phi_atoms[1],
                        phi_atoms[2],
                        phi_atoms[3],
                        psi_atoms[0],
                        psi_atoms[1],
                        psi_atoms[2],
                        psi_atoms[3],
                    )
                    cmap_terms.append(
                        _CmapTerm(
                            torsion_index=torsion_index,
                            map_index=map_index,
                            phi_atoms=phi_atoms,
                            psi_atoms=psi_atoms,
                            phi_types=grid.atom_types_phi,
                            psi_types=grid.atom_types_psi,
                        )
                    )

        if cmap_terms:
            logger.info("Created %d CMAP correction term(s).", len(cmap_terms))
        else:
            logger.warning("CMAP grids present but no matching atom pairs found.")
            cmap_force = None

    if not bond_terms and not angle_terms and not torsion_terms and not vdw_terms and not cmap_terms:
        raise ValueError(
            "No OpenMM terms were created. Force field did not match any "
            "detected bonds, angles, torsions, vdW types, or CMAP pairs."
        )

    if bond_terms:
        system.addForce(bond_force)
    else:
        bond_force = None

    if angle_terms:
        system.addForce(angle_force)
    else:
        angle_force = None

    if torsion_terms:
        system.addForce(torsion_force)
    else:
        torsion_force = None

    if vdw_terms:
        system.addForce(vdw_force)
    else:
        vdw_force = None

    if ub_terms:
        system.addForce(ub_force)
    else:
        ub_force = None

    if cmap_terms:
        system.addForce(cmap_force)
    else:
        cmap_force = None

    integrator, context = self._create_context(system, precision=precision)
    context.setPositions(self._positions(molecule))

    return OpenMMHandle(
        molecule=copy.deepcopy(molecule),
        system=system,
        integrator=integrator,
        context=context,
        bond_force=bond_force,
        angle_force=angle_force,
        torsion_force=torsion_force,
        vdw_force=vdw_force,
        ub_force=ub_force,
        cmap_force=cmap_force,
        bond_terms=bond_terms,
        angle_terms=angle_terms,
        torsion_terms=torsion_terms,
        vdw_terms=vdw_terms,
        ub_terms=ub_terms,
        cmap_terms=cmap_terms,
        exceptions_14=exceptions_14 if use_harmonic and vdw_terms else [],
        functional_form=ff_form,
    )

update_forcefield

update_forcefield(handle: OpenMMHandle, forcefield: ForceField) -> None

Update per-term parameters in an existing OpenMM Context.

Modifies bond, angle, and vdW parameters in-place, then pushes changes to the OpenMM context. Much faster than rebuilding the system from scratch.

Parameters:

Name	Type	Description	Default
handle	OpenMMHandle	An existing :class:OpenMMHandle to update.	required
forcefield	ForceField	New force field parameters to apply.	required

Raises:

Type	Description
ValueError	If the force field's functional form does not match the handle's form, or if a required parameter is missing.

Source code in q2mm/backends/mm/openmm.py

def update_forcefield(self, handle: OpenMMHandle, forcefield: ForceField) -> None:
    """Update per-term parameters in an existing OpenMM Context.

    Modifies bond, angle, and vdW parameters in-place, then pushes
    changes to the OpenMM context.  Much faster than rebuilding the
    system from scratch.

    Args:
        handle: An existing :class:`OpenMMHandle` to update.
        forcefield: New force field parameters to apply.

    Raises:
        ValueError: If the force field's functional form does not match
            the handle's form, or if a required parameter is missing.

    """
    incoming_form = forcefield.functional_form
    if incoming_form is not None and incoming_form != handle.functional_form:
        raise ValueError(
            f"Force field functional form {incoming_form!r} does not match "
            f"the handle's form {handle.functional_form!r}. "
            f"Create a new context instead of reusing this handle."
        )
    use_harmonic = handle.functional_form == FunctionalForm.HARMONIC

    if handle.bond_force is not None:
        for term in handle.bond_terms:
            param = _match_bond(forcefield, term.elements, env_id=term.env_id, ff_row=term.ff_row)
            if param is None:
                raise ValueError(f"Updated force field is missing bond parameter for {term.elements}.")
            if use_harmonic:
                handle.bond_force.setBondParameters(
                    term.force_index,
                    term.atom_i,
                    term.atom_j,
                    ang_to_nm(param.equilibrium),
                    _bond_k_to_harmonic(param.force_constant),
                )
            else:
                handle.bond_force.setBondParameters(
                    term.force_index,
                    term.atom_i,
                    term.atom_j,
                    [_bond_k_to_openmm(param.force_constant), ang_to_nm(param.equilibrium)],
                )
        handle.bond_force.updateParametersInContext(handle.context)

    if handle.angle_force is not None:
        for term in handle.angle_terms:
            param = _match_angle(forcefield, term.elements, env_id=term.env_id, ff_row=term.ff_row)
            if param is None:
                raise ValueError(f"Updated force field is missing angle parameter for {term.elements}.")
            if use_harmonic:
                handle.angle_force.setAngleParameters(
                    term.force_index,
                    term.atom_i,
                    term.atom_j,
                    term.atom_k,
                    np.deg2rad(float(param.equilibrium)),
                    _angle_k_to_harmonic(param.force_constant),
                )
            else:
                handle.angle_force.setAngleParameters(
                    term.force_index,
                    term.atom_i,
                    term.atom_j,
                    term.atom_k,
                    [_angle_k_to_openmm(param.force_constant), np.deg2rad(float(param.equilibrium))],
                )
        handle.angle_force.updateParametersInContext(handle.context)

    if handle.torsion_force is not None:
        for term in handle.torsion_terms:
            params = _match_torsions(
                forcefield, term.elements, env_id=term.env_id, ff_row=term.ff_row, is_improper=term.is_improper
            )
            matched = [p for p in params if p.periodicity == term.periodicity]
            if not matched:
                raise ValueError(
                    f"Updated force field is missing torsion parameter for "
                    f"{term.elements} periodicity={term.periodicity}."
                )
            param = matched[0]
            handle.torsion_force.setTorsionParameters(
                term.force_index,
                term.atom_i,
                term.atom_j,
                term.atom_k,
                term.atom_l,
                param.periodicity,
                np.deg2rad(float(param.phase)),
                _torsion_k_to_openmm(param.force_constant),
            )
        handle.torsion_force.updateParametersInContext(handle.context)

    if handle.vdw_force is not None:
        for term in handle.vdw_terms:
            param = _match_vdw(forcefield, atom_type=term.atom_type, element=term.element, ff_row=term.ff_row)
            if param is None:
                raise ValueError(
                    f"Updated force field is missing vdW parameter for {term.atom_type or term.element}."
                )
            if use_harmonic:
                handle.vdw_force.setParticleParameters(
                    term.particle_index,
                    0.0,
                    _vdw_sigma_nm(param.radius),
                    _vdw_epsilon_to_openmm(param.epsilon),
                )
            else:
                handle.vdw_force.setParticleParameters(
                    term.particle_index,
                    [_vdw_radius_to_openmm(param.radius), _vdw_epsilon_to_openmm(param.epsilon)],
                )

        # Recompute 1-4 exception params from updated particle params
        if use_harmonic and handle.exceptions_14:
            SCNB = 2.0
            for exc in handle.exceptions_14:
                _, sig1, eps1 = handle.vdw_force.getParticleParameters(exc.particle_i)
                _, sig2, eps2 = handle.vdw_force.getParticleParameters(exc.particle_j)
                sig_14 = 0.5 * (sig1 + sig2)
                eps_14 = (eps1 * eps2) ** 0.5 / SCNB
                handle.vdw_force.setExceptionParameters(
                    exc.exception_index, exc.particle_i, exc.particle_j, 0.0, sig_14, eps_14
                )

        handle.vdw_force.updateParametersInContext(handle.context)

    if handle.ub_force is not None:
        for term in handle.ub_terms:
            param = _match_angle(forcefield, term.elements, env_id=term.env_id, ff_row=term.ff_row)
            if param is None:
                raise ValueError(f"Updated force field is missing UB parameter for angle {term.elements}.")
            if param.ub_force_constant is None and param.ub_equilibrium is None:
                raise ValueError(f"Updated force field is missing UB parameter for angle {term.elements}.")
            if param.ub_force_constant is None or param.ub_equilibrium is None:
                raise ValueError(
                    "Inconsistent Urey-Bradley parameters for angle "
                    f"{term.elements} (env_id={term.env_id}, ff_row={term.ff_row}): "
                    "both 'ub_force_constant' and 'ub_equilibrium' must be set or both must be None."
                )
            handle.ub_force.setBondParameters(
                term.force_index,
                term.atom_i,
                term.atom_k,
                ang_to_nm(param.ub_equilibrium),
                _bond_k_to_harmonic(param.ub_force_constant),
            )
        handle.ub_force.updateParametersInContext(handle.context)

export_system_xml

export_system_xml(path: str | Path, structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None) -> Path

Serialize the OpenMM System to XML.

Produces a topology-specific XML file containing the force objects (HarmonicBondForce/CustomBondForce, etc. depending on the functional form) with all per-term parameters. The file can be loaded back with openmm.XmlSerializer.deserialize().

Parameters:

Name	Type	Description	Default
path	str | Path	Output file path.	required
structure	Q2MMMolecule | str | Path | OpenMMHandle	A :class:~q2mm.models.molecule.Q2MMMolecule, path to an XYZ file, or an existing :class:OpenMMHandle.	required
forcefield	ForceField | None	Force field to apply. When structure is not an :class:OpenMMHandle, this is used to build the OpenMM system. When structure is an existing :class:OpenMMHandle, providing a non-None forcefield updates the per-term parameters of that handle; if forcefield is None, the handle's current parameters are used unchanged.	None

Returns:

Name	Type	Description
Path	Path	The resolved output path.

Source code in q2mm/backends/mm/openmm.py

def export_system_xml(
    self,
    path: str | Path,
    structure: Q2MMMolecule | str | Path | OpenMMHandle,
    forcefield: ForceField | None = None,
) -> Path:
    """Serialize the OpenMM System to XML.

    Produces a topology-specific XML file containing the force objects
    (``HarmonicBondForce``/``CustomBondForce``, etc. depending on the
    functional form) with all per-term parameters.  The file can be
    loaded back with ``openmm.XmlSerializer.deserialize()``.

    Args:
        path: Output file path.
        structure (Q2MMMolecule | str | Path | OpenMMHandle): A
            :class:`~q2mm.models.molecule.Q2MMMolecule`, path to an XYZ
            file, or an existing :class:`OpenMMHandle`.
        forcefield: Force field to apply.  When *structure* is not an
            :class:`OpenMMHandle`, this is used to build the OpenMM
            system.  When *structure* is an existing
            :class:`OpenMMHandle`, providing a non-None *forcefield*
            updates the per-term parameters of that handle; if
            *forcefield* is ``None``, the handle's current parameters
            are used unchanged.

    Returns:
        Path: The resolved output path.

    """
    handle = self._prepare_handle(structure, forcefield)
    xml_string = mm.XmlSerializer.serialize(handle.system)
    output = Path(path)
    output.write_text(xml_string, encoding="utf-8")
    return output

load_system_xml staticmethod

load_system_xml(path: str | Path) -> object

Deserialize an OpenMM System from XML.

Parameters:

Name	Type	Description	Default
path	str | Path	Path to the XML file.	required

Returns:

Name	Type	Description
object	object	An openmm.System object.

Source code in q2mm/backends/mm/openmm.py

@staticmethod
def load_system_xml(path: str | Path) -> object:
    """Deserialize an OpenMM System from XML.

    Args:
        path: Path to the XML file.

    Returns:
        object: An ``openmm.System`` object.

    """
    _ensure_openmm()
    xml_string = Path(path).read_text(encoding="utf-8")
    return mm.XmlSerializer.deserialize(xml_string)

energy_and_param_grad

energy_and_param_grad(structure: Q2MMMolecule, forcefield: ForceField) -> tuple[float, ndarray]

Compute energy and analytical gradient w.r.t. FF parameters.

Uses OpenMM's addEnergyParameterDerivative() on CustomForce objects to get exact dE/d(param) for bond, angle, and torsion parameters. vdW parameters use per-particle values that cannot be differentiated via global parameters, so their gradients are computed via central finite differences automatically.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule	Molecular structure.	required
forcefield	ForceField	Force field parameters.	required

Returns:

Type	Description
tuple[float, ndarray]	tuple[float, np.ndarray]: (energy, grad) where energy is in kcal/mol and grad has the same length as forcefield.get_param_vector().

Source code in q2mm/backends/mm/openmm.py

def energy_and_param_grad(self, structure: Q2MMMolecule, forcefield: ForceField) -> tuple[float, np.ndarray]:
    """Compute energy and analytical gradient w.r.t. FF parameters.

    Uses OpenMM's ``addEnergyParameterDerivative()`` on ``CustomForce``
    objects to get exact dE/d(param) for bond, angle, and torsion
    parameters.  vdW parameters use per-particle values that cannot
    be differentiated via global parameters, so their gradients are
    computed via central finite differences automatically.

    Args:
        structure (Q2MMMolecule): Molecular structure.
        forcefield (ForceField): Force field parameters.

    Returns:
        tuple[float, np.ndarray]: ``(energy, grad)`` where ``energy``
            is in kcal/mol and ``grad`` has the same length as
            ``forcefield.get_param_vector()``.

    """
    molecule = _as_molecule(structure)
    diff = self._create_diff_handle(molecule, forcefield)

    state = diff.context.getState(getEnergy=True, getParameterDerivatives=True)
    energy = float(state.getPotentialEnergy().value_in_unit(unit.kilocalories_per_mole))
    derivs = state.getEnergyParameterDerivatives()

    param_vector = forcefield.get_param_vector()
    grad = np.zeros(len(param_vector))

    for name, pv_idx, unit_factor in zip(
        diff.param_names, diff.param_vector_indices, diff.grad_unit_factors, strict=True
    ):
        deriv_openmm = derivs[name]  # dE_kJ/dp_openmm
        grad[pv_idx] = kj_to_kcal(deriv_openmm * unit_factor)

    # vdW parameters use per-particle values without global-parameter
    # derivatives.  Supplement with central finite differences.
    # Reuse a single OpenMMHandle to avoid rebuilding the OpenMM
    # context for each perturbation.  Use double precision on GPU
    # so the finite differences are not lost to float32 rounding.
    if forcefield.vdws:
        vdw_start = 2 * len(forcefield.bonds) + 2 * len(forcefield.angles) + len(forcefield.torsions)
        vdw_end = vdw_start + 2 * len(forcefield.vdws)
        step = 1e-4
        handle = self.create_context(molecule, forcefield, precision="double")
        for i in range(vdw_start, vdw_end):
            pv_plus = param_vector.copy()
            pv_plus[i] += step
            pv_minus = param_vector.copy()
            pv_minus[i] -= step
            e_plus = self.energy(handle, forcefield.with_params(pv_plus))
            e_minus = self.energy(handle, forcefield.with_params(pv_minus))
            grad[i] = (e_plus - e_minus) / (2.0 * step)

    return energy, grad

energy

energy(structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None) -> float

Calculate MM energy in kcal/mol.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule | str | Path | OpenMMHandle	Molecule, XYZ path, or :class:OpenMMHandle.	required
forcefield	ForceField | None	Force field to apply. Auto-generated if None.	None

Returns:

Name	Type	Description
float	float	Potential energy in kcal/mol.

Source code in q2mm/backends/mm/openmm.py

def energy(
    self, structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None
) -> float:
    """Calculate MM energy in kcal/mol.

    Args:
        structure (Q2MMMolecule | str | Path | OpenMMHandle): Molecule,
            XYZ path, or :class:`OpenMMHandle`.
        forcefield: Force field to apply. Auto-generated if ``None``.

    Returns:
        float: Potential energy in kcal/mol.

    """
    handle = self._prepare_handle(structure, forcefield)
    state = handle.context.getState(getEnergy=True)
    return float(state.getPotentialEnergy().value_in_unit(unit.kilocalories_per_mole))

minimize

minimize(structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None, tolerance: float = 1.0, max_iterations: int = 200) -> tuple

Energy-minimize structure using L-BFGS.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule | str | Path | OpenMMHandle	Molecule, XYZ path, or :class:OpenMMHandle.	required
forcefield	ForceField | None	Force field to apply. Auto-generated if None.	None
tolerance	float	Energy convergence tolerance in kJ/mol.	1.0
max_iterations	int	Maximum minimization steps.	200

Returns:

Type	Description
tuple	tuple[float, list[str], np.ndarray]: (energy, atoms, coords) where energy is in kcal/mol and coords are in Å.

Source code in q2mm/backends/mm/openmm.py

def minimize(
    self,
    structure: Q2MMMolecule | str | Path | OpenMMHandle,
    forcefield: ForceField | None = None,
    tolerance: float = 1.0,
    max_iterations: int = 200,
) -> tuple:
    """Energy-minimize structure using L-BFGS.

    Args:
        structure (Q2MMMolecule | str | Path | OpenMMHandle): Molecule,
            XYZ path, or :class:`OpenMMHandle`.
        forcefield: Force field to apply. Auto-generated if ``None``.
        tolerance: Energy convergence tolerance in kJ/mol.
        max_iterations: Maximum minimization steps.

    Returns:
        tuple[float, list[str], np.ndarray]: ``(energy, atoms, coords)``
            where energy is in kcal/mol and coords are in Å.

    """
    handle = self._prepare_handle(structure, forcefield)
    mm.LocalEnergyMinimizer.minimize(handle.context, tolerance, max_iterations)
    state = handle.context.getState(getEnergy=True, getPositions=True)
    coords = np.array(state.getPositions(asNumpy=True).value_in_unit(unit.angstrom))
    handle.molecule.geometry = coords
    energy = float(state.getPotentialEnergy().value_in_unit(unit.kilocalories_per_mole))
    return energy, list(handle.molecule.symbols), coords

hessian

hessian(structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None, step: float = 0.0001) -> ndarray

Finite-difference Hessian in canonical units (Hartree/Bohr²).

Internally computed in kJ/mol/nm² (OpenMM native) and converted to Hartree/Bohr² before returning, matching the canonical unit contract defined in :class:~q2mm.backends.base.MMEngine.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule | str | Path | OpenMMHandle	Molecule, XYZ path, or :class:OpenMMHandle.	required
forcefield	ForceField | None	Force field to apply. Auto-generated if None.	None
step	float	Finite-difference displacement in nm.	0.0001

Returns:

Type	Description
ndarray	np.ndarray: Shape (3N, 3N) Hessian in Hartree/Bohr².

Source code in q2mm/backends/mm/openmm.py

def hessian(
    self,
    structure: Q2MMMolecule | str | Path | OpenMMHandle,
    forcefield: ForceField | None = None,
    step: float = 1.0e-4,
) -> np.ndarray:
    """Finite-difference Hessian in canonical units (Hartree/Bohr²).

    Internally computed in kJ/mol/nm² (OpenMM native) and converted
    to Hartree/Bohr² before returning, matching the canonical unit
    contract defined in :class:`~q2mm.backends.base.MMEngine`.

    Args:
        structure (Q2MMMolecule | str | Path | OpenMMHandle): Molecule,
            XYZ path, or :class:`OpenMMHandle`.
        forcefield: Force field to apply. Auto-generated if ``None``.
        step: Finite-difference displacement in nm.

    Returns:
        np.ndarray: Shape ``(3N, 3N)`` Hessian in Hartree/Bohr².

    """
    handle = self._prepare_handle(structure, forcefield)
    positions = np.array(
        handle.context.getState(getPositions=True).getPositions(asNumpy=True).value_in_unit(unit.nanometer)
    )
    n_atoms = positions.shape[0]
    hessian = np.zeros((3 * n_atoms, 3 * n_atoms))

    for atom_index in range(n_atoms):
        for coord_index in range(3):
            column = 3 * atom_index + coord_index

            displaced_plus = positions.copy()
            displaced_minus = positions.copy()
            displaced_plus[atom_index, coord_index] += step
            displaced_minus[atom_index, coord_index] -= step

            handle.context.setPositions(displaced_plus * unit.nanometer)
            forces_plus = np.array(
                handle.context.getState(getForces=True)
                .getForces(asNumpy=True)
                .value_in_unit(unit.kilojoule_per_mole / unit.nanometer)
            )

            handle.context.setPositions(displaced_minus * unit.nanometer)
            forces_minus = np.array(
                handle.context.getState(getForces=True)
                .getForces(asNumpy=True)
                .value_in_unit(unit.kilojoule_per_mole / unit.nanometer)
            )

            hessian[:, column] = -((forces_plus - forces_minus) / (2.0 * step)).reshape(-1)

    handle.context.setPositions(positions * unit.nanometer)
    hessian_symmetric = 0.5 * (hessian + hessian.T)

    # Convert from OpenMM native kJ/mol/nm² to canonical Hartree/Bohr²
    return hessian_symmetric * hessian_kjmolnm2_to_au(1.0)

frequencies

frequencies(structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None) -> list[float]

Approximate harmonic frequencies in cm⁻¹ from the numerical Hessian.

Parameters:

Name	Type	Description	Default
structure	Q2MMMolecule | str | Path | OpenMMHandle	Molecule, XYZ path, or :class:OpenMMHandle.	required
forcefield	ForceField | None	Force field to apply. Auto-generated if None.	None

Returns:

Type	Description
list[float]	list[float]: Vibrational frequencies in cm⁻¹.

Source code in q2mm/backends/mm/openmm.py

def frequencies(
    self, structure: Q2MMMolecule | str | Path | OpenMMHandle, forcefield: ForceField | None = None
) -> list[float]:
    """Approximate harmonic frequencies in cm⁻¹ from the numerical Hessian.

    Args:
        structure (Q2MMMolecule | str | Path | OpenMMHandle): Molecule,
            XYZ path, or :class:`OpenMMHandle`.
        forcefield: Force field to apply. Auto-generated if ``None``.

    Returns:
        list[float]: Vibrational frequencies in cm⁻¹.

    """
    from q2mm.models.hessian import hessian_to_frequencies

    handle = self._prepare_handle(structure, forcefield)
    hessian_au = self.hessian(handle)  # Hartree/Bohr²
    return hessian_to_frequencies(hessian_au, list(handle.molecule.symbols))

detect_best_platform

detect_best_platform() -> str

Return the name of the fastest available OpenMM platform.

Platform preference order: CUDA > OpenCL > CPU > Reference.

Logs a warning when CUDA is unavailable and the function falls back to OpenCL on a system with an NVIDIA GPU — OpenCL on modern NVIDIA GPUs gives very poor utilisation (~14%). The warning is suppressed on non-NVIDIA systems where OpenCL may be the intended backend.

Returns:

Name	Type	Description
str	str	Name of the best available platform.

Raises:

Type	Description
ImportError	If OpenMM is not installed.

Source code in q2mm/backends/mm/openmm.py

def detect_best_platform() -> str:
    """Return the name of the fastest available OpenMM platform.

    Platform preference order: CUDA > OpenCL > CPU > Reference.

    Logs a warning when CUDA is unavailable and the function falls back
    to OpenCL on a system with an NVIDIA GPU — OpenCL on modern NVIDIA
    GPUs gives very poor utilisation (~14%).  The warning is suppressed
    on non-NVIDIA systems where OpenCL may be the intended backend.

    Returns:
        str: Name of the best available platform.

    Raises:
        ImportError: If OpenMM is not installed.

    """
    _ensure_openmm()
    available = {mm.Platform.getPlatform(i).getName() for i in range(mm.Platform.getNumPlatforms())}
    for name in _PLATFORM_PRIORITY:
        if name in available:
            if name == "OpenCL" and "CUDA" not in available:
                # Only warn on NVIDIA GPUs where CUDA should be available
                _nvidia_present = False
                try:
                    import subprocess

                    result = subprocess.run(
                        ["nvidia-smi", "--query-gpu=name", "--format=csv,noheader"],
                        capture_output=True,
                        text=True,
                        timeout=5,
                    )
                    _nvidia_present = result.returncode == 0 and bool(result.stdout.strip())
                except Exception:
                    pass
                if _nvidia_present:
                    logger.warning(
                        "CUDA platform not available, falling back to OpenCL. "
                        "GPU utilization will be poor (~14%%). "
                        "Consider installing OpenMM-CUDA-12 or using WSL2."
                    )
            return name
    # Fallback — shouldn't happen since OpenMM always has Reference
    return mm.Platform.getPlatform(0).getName()  # pragma: no cover