Dashboard Core — QM/MM

Moved to dense_evolution.qmmm — this page documents the backward-compatible re-export only; dashboard_core.qmmm still works unchanged, but the real implementation and any new QM/MM work (region partitioning, embedding) now live in the library itself, alongside real Hellmann-Feynman nuclear forces and a real Velocity-Verlet MD step for this project's molecule catalog — no fabricated geometry, no placeholder forces. Backs Composer's QM/MM force and MD-trajectory panels.

from dashboard_core.qmmm import compute_hellmann_feynman_forces, run_md_trajectory
from dashboard_core.hamiltonians import MOLECULE_CATALOG

# MOLECULE_CATALOG keys are descriptive strings, not bare element symbols:
h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]

forces = compute_hellmann_feynman_forces(h2)
print(forces["energy_hartree"])  # -1.1373 (Hartree-Fock ground state)
print(forces["force_norm"])      # ~0.0154 Hartree/Angstrom -- small residual at equilibrium

trajectory = run_md_trajectory(h2, n_steps=3, dt_fs=0.5)
print(trajectory["force_norm"])  # decreasing as the bond relaxes toward equilibrium

qmmm

Backward-compatibility re-export. The real implementation moved to dense_evolution.qmmm.forces (Dense-Evolution issue #283 -- qmmm now lives in the library, alongside real QM/MM region partitioning/embedding, instead of being buried inside the dashboard tool). Nothing in dashboard code needed to change: from dashboard_core.qmmm import ... still works.

compute_hellmann_feynman_forces

compute_hellmann_feynman_forces(
    name: str,
    statevector=None,
    mapping: str = "jordan_wigner",
    geometry=None,
    fd_step_angstrom: float = 0.001,
)

Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of MOLECULE_CATALOG[name]: F = -d/dR, with H(R) this project's own real Hamiltonian (build_molecular_hamiltonian) and psi held fixed. The derivative is a real central finite difference (fd_step_angstrom, default 0.001 A -- verified converged against 0.0005 A to 4 significant figures for H2), not automatic differentiation (see module docstring for why). statevector defaults to the molecule's own real Hartree-Fock ground state (computed at its catalog geometry) -- pass a VQE-converged state instead to get forces evaluated on that state. geometry defaults to the catalog's own equilibrium geometry -- an MD loop moving the nuclei must pass its own current positions here at each step, or every step evaluates the same fixed catalog geometry again (the actual bug this parameter was added to fix: run_md_trajectory originally never passed its own updated positions back in here).

Cost (prog.txt, dashboard_core audit point 4a): the central-difference derivative evaluates energy_at 6n_atoms+1 times (H2's 2 atoms -> 13 Hamiltonian builds per call, each at a genuinely different geometry). This looks like it should be cacheable -- build_molecular_hamiltonian already caches by exact geometry -- but it isn't in practice: every one of the 13 geometries differs by fd_step_angstrom, so every call is a cache miss. Measured directly on H2 (dhf/PennyLane path) before deciding not to add a "cache the Pauli-term basis" layer here: HF + fermion-to-qubit mapping took 0.130s, Pauli-term extraction 0.0005s, dense matrix assembly 0.0031s -- the Hartree-Fock solve itself is 96%+ of the cost, not the bookkeeping after it, so caching the Pauli-term structure would save a few percent at best, not the 6n_atoms multiplier prog.txt's framing suggests. A real geometry change requires a real HF re-solve regardless of how its output gets packaged afterward -- see the module docstring's own account of why analytic differentiation (which WOULD avoid re-solving HF this many times) was tried and dropped for a real cross-platform PennyLane/ autograd bug, not reattempted here.

Parameters:

Name	Type	Description	Default
name	str	A key from MOLECULE_CATALOG (e.g. one of list(MOLECULE_CATALOG) -- these are descriptive strings like "H2 (Idrogeno) - R = 0.7414 A [equilibrio reale]", not bare element symbols like "H2").	required
statevector	see description above.		None
mapping	see description above.		None
geometry	see description above.		None
fd_step_angstrom	see description above.		None

Returns:

Type	Description
dict	name, symbols, energy_hartree, positions_angstrom, forces_hartree_per_angstrom, force_norm.

Examples:

>>> from dense_evolution.qmmm import compute_hellmann_feynman_forces
>>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
>>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
>>> result = compute_hellmann_feynman_forces(h2)
>>> round(result['energy_hartree'], 4)
-1.1373
>>> 0.01 < result['force_norm'] < 0.02  # small residual at equilibrium, not exactly zero
True

Source code in dense_evolution/qmmm/forces.py

def compute_hellmann_feynman_forces(name: str, statevector=None, mapping: str = "jordan_wigner",
                                     geometry=None, fd_step_angstrom: float = 0.001):
    """Real Hellmann-Feynman forces (Hartree/Angstrom) on every nucleus of
    MOLECULE_CATALOG[name]: F = -d<psi|H(R)|psi>/dR, with H(R) this
    project's own real Hamiltonian (build_molecular_hamiltonian) and psi
    held fixed. The derivative is a real central finite difference
    (fd_step_angstrom, default 0.001 A -- verified converged against
    0.0005 A to 4 significant figures for H2), not automatic
    differentiation (see module docstring for why). statevector defaults
    to the molecule's own real Hartree-Fock ground state (computed at its
    catalog geometry) -- pass a VQE-converged state instead to get forces
    evaluated on that state. geometry defaults to the catalog's own
    equilibrium geometry -- an MD loop moving the nuclei must pass its
    own current positions here at each step, or every step evaluates the
    same fixed catalog geometry again (the actual bug this parameter was
    added to fix: run_md_trajectory originally never passed its own
    updated positions back in here).

    Cost (prog.txt, dashboard_core audit point 4a): the central-difference
    derivative evaluates energy_at 6*n_atoms+1 times (H2's 2 atoms -> 13
    Hamiltonian builds per call, each at a genuinely different geometry).
    This looks like it should be cacheable -- build_molecular_hamiltonian
    already caches by exact geometry -- but it isn't in practice: every
    one of the 13 geometries differs by fd_step_angstrom, so every call
    is a cache miss. Measured directly on H2 (dhf/PennyLane path) before
    deciding not to add a "cache the Pauli-term basis" layer here: HF +
    fermion-to-qubit mapping took 0.130s, Pauli-term extraction 0.0005s,
    dense matrix assembly 0.0031s -- the Hartree-Fock solve itself is
    96%+ of the cost, not the bookkeeping after it, so caching the
    Pauli-term structure would save a few percent at best, not the
    6*n_atoms multiplier prog.txt's framing suggests. A real geometry
    change requires a real HF re-solve regardless of how its output gets
    packaged afterward -- see the module docstring's own account of why
    analytic differentiation (which WOULD avoid re-solving HF this many
    times) was tried and dropped for a real cross-platform PennyLane/
    autograd bug, not reattempted here.

    Parameters
    ----------
    name : str
        A key from `MOLECULE_CATALOG` (e.g. one of `list(MOLECULE_CATALOG)`
        -- these are descriptive strings like `"H2 (Idrogeno) - R = 0.7414
        A [equilibrio reale]"`, not bare element symbols like `"H2"`).
    statevector, mapping, geometry, fd_step_angstrom : see description above.

    Returns
    -------
    dict
        `name`, `symbols`, `energy_hartree`, `positions_angstrom`,
        `forces_hartree_per_angstrom`, `force_norm`.

    Examples
    --------
    >>> from dense_evolution.qmmm import compute_hellmann_feynman_forces
    >>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
    >>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
    >>> result = compute_hellmann_feynman_forces(h2)
    >>> round(result['energy_hartree'], 4)
    -1.1373
    >>> 0.01 < result['force_norm'] < 0.02  # small residual at equilibrium, not exactly zero
    True
    """
    MOLECULE_CATALOG, build_molecular_hamiltonian = _import_dashboard_hamiltonians()
    if name not in MOLECULE_CATALOG:
        raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
    spec = MOLECULE_CATALOG[name]
    symbols = spec["symbols"]
    if geometry is None:
        geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
    charge = spec["charge"]
    # BUG FIX (prog.txt, dashboard_core audit point 3a): these two were
    # never read from spec at all, so any catalog entry needing active-
    # space reduction (Si2 is the reason it's in MOLECULE_CATALOG) built
    # the FULL Hamiltonian instead of the reduced one here -- for Si2
    # specifically, 36 qubits instead of the intended 8, which
    # SafeMemoryGuard correctly refuses to allocate. Both
    # _reference_ground_state and energy_at below need these to build
    # the same, correctly-reduced Hamiltonian this molecule's other
    # dashboard panels (VQE, energy scan) already use.
    active_electrons = spec.get("active_electrons")
    active_orbitals = spec.get("active_orbitals")

    unknown = [s for s in symbols if s not in ATOMIC_MASSES_AMU]
    if unknown:
        raise ValueError(f"no real atomic mass on file for {unknown} -- add to "
                          f"ATOMIC_MASSES_AMU before using this molecule here")

    if statevector is None:
        statevector, _gs_energy, _n_qubits = _reference_ground_state(
            symbols, geometry, charge, mapping, active_electrons, active_orbitals)
    sv = np.asarray(statevector, dtype=np.complex128)
    geometry = np.asarray(geometry, dtype=np.float64)

    def energy_at(geom):
        h_matrix, _n_qubits = build_molecular_hamiltonian(
            symbols, geom, charge, mapping, active_electrons, active_orbitals)
        return float(np.real(np.vdot(sv, h_matrix @ sv)))

    energy = energy_at(geometry)
    forces = np.zeros_like(geometry)
    h = fd_step_angstrom
    for i in range(geometry.shape[0]):
        for j in range(3):
            geom_plus = geometry.copy()
            geom_plus[i, j] += h
            geom_minus = geometry.copy()
            geom_minus[i, j] -= h
            forces[i, j] = -(energy_at(geom_plus) - energy_at(geom_minus)) / (2 * h)

    return {
        "name": name,
        "symbols": symbols,
        "energy_hartree": energy,
        "positions_angstrom": geometry.tolist(),
        "forces_hartree_per_angstrom": forces.tolist(),
        "force_norm": float(np.linalg.norm(forces)),
    }

md_step

md_step(
    positions_angstrom,
    velocities_angstrom_per_fs,
    forces_hartree_per_angstrom,
    symbols,
    dt_fs: float = 0.5,
)

One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)dt/2, r(t+dt) = r(t) + v(t+dt/2)dt) using the real Hellmann-Feynman forces above and each atom's real atomic mass -- ordinary classical Newtonian mechanics (F=ma), nothing invented. Positions in Angstrom, velocities in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds.

Examples:

>>> import numpy as np
>>> from dense_evolution.qmmm import md_step
>>> positions = np.array([[0, 0, 0.0], [0, 0, 0.7414]])   # H2 at equilibrium
>>> velocities = np.zeros_like(positions)
>>> forces = np.array([[0, 0, 0.0109], [0, 0, -0.0109]])  # restoring force, pulling atoms together
>>> new_pos, new_vel, accel = md_step(positions, velocities, forces, ['H', 'H'], dt_fs=0.5)
>>> bool(new_pos[1, 2] < 0.7414)  # the second atom moved toward the first
True

Source code in dense_evolution/qmmm/forces.py

def md_step(positions_angstrom, velocities_angstrom_per_fs, forces_hartree_per_angstrom,
            symbols, dt_fs: float = 0.5):
    """One real Velocity-Verlet half-step (v(t+dt/2) = v(t) + a(t)*dt/2,
    r(t+dt) = r(t) + v(t+dt/2)*dt) using the real Hellmann-Feynman forces
    above and each atom's real atomic mass -- ordinary classical Newtonian
    mechanics (F=ma), nothing invented. Positions in Angstrom, velocities
    in Angstrom/fs, forces in Hartree/Angstrom, dt in femtoseconds.

    Examples
    --------
    >>> import numpy as np
    >>> from dense_evolution.qmmm import md_step
    >>> positions = np.array([[0, 0, 0.0], [0, 0, 0.7414]])   # H2 at equilibrium
    >>> velocities = np.zeros_like(positions)
    >>> forces = np.array([[0, 0, 0.0109], [0, 0, -0.0109]])  # restoring force, pulling atoms together
    >>> new_pos, new_vel, accel = md_step(positions, velocities, forces, ['H', 'H'], dt_fs=0.5)
    >>> bool(new_pos[1, 2] < 0.7414)  # the second atom moved toward the first
    True
    """
    positions = np.asarray(positions_angstrom, dtype=np.float64)
    velocities = np.asarray(velocities_angstrom_per_fs, dtype=np.float64)
    forces = np.asarray(forces_hartree_per_angstrom, dtype=np.float64)
    masses = np.array([ATOMIC_MASSES_AMU[s] for s in symbols], dtype=np.float64)

    accel = ACCEL_CONVERSION * forces / masses[:, None]
    velocities_half = velocities + 0.5 * accel * dt_fs
    positions_new = positions + velocities_half * dt_fs
    return positions_new, velocities_half, accel

run_md_trajectory

run_md_trajectory(
    name: str,
    n_steps: int,
    dt_fs: float = 0.5,
    mapping: str = "jordan_wigner",
    recompute_electronic_state: bool = False,
    fd_step_angstrom: float = 0.001,
)

Real, minimal ab-initio-forces MD trajectory: at each step, real Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the real nuclear positions/velocities via real Velocity-Verlet (md_step). Starts from rest (zero initial velocities) at the catalog's real equilibrium geometry.

fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at every step -- previously not exposed here at all, silently using that function's own default (0.001 A) with no way for a caller to ask for a different finite-difference step (e.g. a molecule with an unusually steep energy landscape, where the default step size isn't the one already verified converged for H2).

recompute_electronic_state=False (default) holds the electronic state fixed at the initial Hartree-Fock reference through the whole trajectory -- forces stay exact only close to the starting geometry (a real, explicitly-stated approximation, not a fabricated one). True ab-initio MD (re-solving Hartree-Fock at every step's new geometry) is available by setting this True, at real, substantial extra cost per step.

Examples:

>>> from dense_evolution.qmmm import run_md_trajectory
>>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
>>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
>>> traj = run_md_trajectory(h2, n_steps=3, dt_fs=0.5)
>>> traj['step']
[0, 1, 2]
>>> len(traj['force_norm'])
3

Source code in dense_evolution/qmmm/forces.py

def run_md_trajectory(name: str, n_steps: int, dt_fs: float = 0.5, mapping: str = "jordan_wigner",
                       recompute_electronic_state: bool = False, fd_step_angstrom: float = 0.001):
    """Real, minimal ab-initio-forces MD trajectory: at each step, real
    Hellmann-Feynman forces (compute_hellmann_feynman_forces) move the
    real nuclear positions/velocities via real Velocity-Verlet (md_step).
    Starts from rest (zero initial velocities) at the catalog's real
    equilibrium geometry.

    fd_step_angstrom: forwarded to compute_hellmann_feynman_forces at
    every step -- previously not exposed here at all, silently using
    that function's own default (0.001 A) with no way for a caller to
    ask for a different finite-difference step (e.g. a molecule with an
    unusually steep energy landscape, where the default step size isn't
    the one already verified converged for H2).

    recompute_electronic_state=False (default) holds the electronic state
    fixed at the initial Hartree-Fock reference through the whole
    trajectory -- forces stay exact only close to the starting geometry
    (a real, explicitly-stated approximation, not a fabricated one).
    True ab-initio MD (re-solving Hartree-Fock at every step's new
    geometry) is available by setting this True, at real, substantial
    extra cost per step.

    Examples
    --------
    >>> from dense_evolution.qmmm import run_md_trajectory
    >>> from dashboard_core.hamiltonians import MOLECULE_CATALOG
    >>> h2 = [k for k in MOLECULE_CATALOG if k.startswith("H2 ")][0]
    >>> traj = run_md_trajectory(h2, n_steps=3, dt_fs=0.5)
    >>> traj['step']
    [0, 1, 2]
    >>> len(traj['force_norm'])
    3
    """
    MOLECULE_CATALOG, _build_molecular_hamiltonian = _import_dashboard_hamiltonians()
    if name not in MOLECULE_CATALOG:
        raise ValueError(f"unknown molecule {name!r}; available: {sorted(MOLECULE_CATALOG)}")
    spec = MOLECULE_CATALOG[name]
    symbols = spec["symbols"]
    geometry = spec["geometry"]() if callable(spec["geometry"]) else spec["geometry"]
    charge = spec["charge"]
    active_electrons = spec.get("active_electrons")
    active_orbitals = spec.get("active_orbitals")

    statevector, _gs_energy, _n_qubits = _reference_ground_state(
        symbols, geometry, charge, mapping, active_electrons, active_orbitals)
    positions = np.asarray(geometry, dtype=np.float64)
    velocities = np.zeros_like(positions)

    trajectory = {"step": [], "time_fs": [], "positions_angstrom": [], "energy_hartree": [], "force_norm": []}
    for step in range(n_steps):
        result = compute_hellmann_feynman_forces(name, statevector, mapping=mapping, geometry=positions,
                                                  fd_step_angstrom=fd_step_angstrom)
        forces = np.asarray(result["forces_hartree_per_angstrom"])
        trajectory["step"].append(step)
        trajectory["time_fs"].append(step * dt_fs)
        trajectory["positions_angstrom"].append(positions.tolist())
        trajectory["energy_hartree"].append(result["energy_hartree"])
        trajectory["force_norm"].append(result["force_norm"])

        positions, velocities, _accel = md_step(positions, velocities, forces, symbols, dt_fs=dt_fs)
        _assert_no_nuclear_collision(positions, step, dt_fs)

        if recompute_electronic_state:
            statevector, _gs_energy, _n_qubits = _reference_ground_state(
                symbols, positions, charge, mapping, active_electrons, active_orbitals)

    return trajectory

---

See also: dashboard_core.hamiltonians for MOLECULE_CATALOG and build_molecular_hamiltonian, which this module's forces are derived from.