OpenMM Engine¶

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, vdw_force: object | None, bond_terms: list[_BondTerm], angle_terms: list[_AngleTerm], vdw_terms: list[_VdwTerm], 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.
`vdw_force`	`object \| None`	The OpenMM vdW force object, or `None` if no vdW terms.
`bond_terms`	`list[_BondTerm]`	Mapping of molecule bonds to force indices.
`angle_terms`	`list[_AngleTerm]`	Mapping of molecule angles to force indices.
`vdw_terms`	`list[_VdwTerm]`	Mapping of atoms to vdW particle 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)

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"`). Auto-selected if `None`.	`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):
    """Initialize the OpenMM engine.

    Args:
        platform_name: OpenMM platform to use (e.g. ``"CPU"``,
            ``"CUDA"``). Auto-selected if ``None``.

    Raises:
        ImportError: If OpenMM is not installed.
    """
    _ensure_openmm()
    self._platform_name = platform_name

name `property` ¶

name: str

Human-readable engine name.

Returns:

Name	Type	Description
`str`	`str`	`"OpenMM"`.

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

create_context ¶

create_context(structure, forcefield: ForceField | 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`

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, forcefield: ForceField | 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``.

    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")

    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(
            "step(r - rc) * epsilon*(-2.25*(rv/r)^6 + 184000*exp(-12*r/rv))"
            " + step(rc - r) * epsilon*184000*exp(-12*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,
                float(param.equilibrium) * 0.1,
                _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), float(param.equilibrium) * 0.1],
            )
        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 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)

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

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

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

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

    integrator, context = self._create_context(system)
    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,
        vdw_force=vdw_force,
        bond_terms=bond_terms,
        angle_terms=angle_terms,
        vdw_terms=vdw_terms,
        exceptions_14=exceptions_14 if use_harmonic and vdw_terms else [],
        functional_form=ff_form,
    )

update_forcefield ¶

update_forcefield(handle: OpenMMHandle, forcefield: ForceField)

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):
    """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,
                    float(param.equilibrium) * 0.1,
                    _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), float(param.equilibrium) * 0.1],
                )
        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.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)

export_system_xml ¶

export_system_xml(path: str | Path, structure, 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,
    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 ¶

energy(structure, 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, 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, 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, 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, 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, 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 * KJMOLNM2_TO_HESSIAN_AU

frequencies ¶

frequencies(structure, 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, 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⁻¹.
    """
    handle = self._prepare_handle(structure, forcefield)
    hessian_au = self.hessian(handle)  # Hartree/Bohr²

    # Convert Hartree/Bohr² → J/m² (per molecule, not per mol)
    bohr_to_m = BOHR_TO_ANG * 1e-10
    hessian_si = hessian_au * HARTREE_TO_J / (bohr_to_m**2)

    masses = np.array([MASSES[symbol] * AMU_TO_KG for symbol in handle.molecule.symbols], dtype=float)
    mass_vector = np.repeat(masses, 3)
    mass_weighted = hessian_si / np.sqrt(np.outer(mass_vector, mass_vector))

    eigenvalues = np.linalg.eigvalsh(mass_weighted)
    frequencies = np.sign(eigenvalues) * np.sqrt(np.abs(eigenvalues))
    frequencies /= 2.0 * np.pi * SPEED_OF_LIGHT_MS * 100.0
    return frequencies.tolist()

OpenMM Engine¶

openmm ¶

OpenMMHandle dataclass ¶

OpenMMEngine ¶

name property ¶

supported_functional_forms ¶

is_available ¶

supports_runtime_params ¶

create_context ¶

update_forcefield ¶

export_system_xml ¶

load_system_xml staticmethod ¶

energy ¶

minimize ¶

hessian ¶

frequencies ¶

OpenMMHandle `dataclass` ¶

name `property` ¶

load_system_xml `staticmethod` ¶