Algorithms¶

Variational algorithms built on the same IR, ansatz vocabulary and gradients as everything else — so an ansatz you registered, a gate you defined, or a backend you wrote works in all of them without any of them knowing about it.

Every one of these takes its ansatz, feature map or operator pool as an argument and must actually use it. tests/test_injection.py injects two of different sizes and asserts the parameter count follows, because a constructor that accepts ansatz= and silently ignores it looks identical from the outside.

`qmlkit.algorithms.vqe`¶

Ground-state energy by variational minimisation. Worked end to end in study 4, including the case where a too-shallow ansatz converges confidently to an energy 601 mHa wrong.

vqe ¶

The variational quantum eigensolver.

VQE is a loop, not an architecture: prepare a parameterised state, measure an energy, step downhill, repeat. Everything that makes one VQE differ from another — the ansatz, the optimiser, the gradient method, the shot budget — is therefore an argument, and the class itself is thin on purpose.

from qmlkit.algorithms import VQE, ising_hamiltonian

H = ising_hamiltonian(4, j=1.0, h=0.5)
result = VQE(H, n_qubits=4).run(seed=0)
print(result.energy, result.error_vs_exact)

VQEResult `dataclass` ¶

VQEResult(
    energy: float,
    theta: NDArray[Any],
    history: list[float] = list(),
    exact: float | None = None,
    n_evaluations: int = 0,
)

What a run produced, and how good it actually is.

error_vs_exact `property` ¶

error_vs_exact: float | None

Absolute error against dense diagonalisation, when that was computed.

VQE ¶

VQE(
    hamiltonian: Observable,
    ansatz: Ansatz | None = None,
    n_qubits: int | None = None,
    optimizer: str | Optimizer = "rotosolve",
    gradient: str = "auto",
    backend: BackendLike = None,
    shots: int | None = None,
)

Minimise <H> over a parameterised state.

Parameters:

Name	Type	Description	Default
`hamiltonian`	`Observable`	Any observable. :mod:`qmlkit.algorithms.hamiltonians` has constructors.	required
`ansatz`	`Ansatz \| None`	The trial state. Defaults to a hardware-efficient circuit wide enough for the Hamiltonian's support — replaceable with anything, including one you invented, because an `Ansatz` is the only contract.	`None`
`optimizer`	`str \| Optimizer`	A name from :data:`OPTIMIZERS` or any `fn(loss, theta0, **kw)` returning `(theta, history)`.	`'rotosolve'`
`gradient`	`str`	Passed through to :func:`qmlkit.grad`. Only consulted by gradient-based optimisers; Rotosolve and SPSA never ask for one.	`'auto'`

Source code in src/qmlkit/algorithms/vqe.py

def __init__(
    self,
    hamiltonian: Observable,
    ansatz: Ansatz | None = None,
    n_qubits: int | None = None,
    optimizer: str | Optimizer = "rotosolve",
    gradient: str = "auto",
    backend: BackendLike = None,
    shots: int | None = None,
) -> None:
    support = observable_support(hamiltonian)
    width = n_qubits or (max(support) + 1 if support else 1)
    if ansatz is not None and ansatz.n_qubits < width:
        raise ValueError(
            f"the ansatz has {ansatz.n_qubits} qubits but the Hamiltonian acts on {width}"
        )
    self.hamiltonian = hamiltonian
    self.ansatz = ansatz or hardware_efficient(width, n_layers=2)
    self.n_qubits = self.ansatz.n_qubits
    self.optimizer = optimizer
    self.gradient = gradient
    self.backend = backend
    self.shots = shots
    self._spec = self.ansatz.build()
    self.n_evaluations = 0

energy ¶

energy(theta: Sequence[float]) -> float

<H> at these parameters.

Source code in src/qmlkit/algorithms/vqe.py

def energy(self, theta: Sequence[float]) -> float:
    """``<H>`` at these parameters."""
    self.n_evaluations += 1
    value = expectation(
        self._spec,
        self.hamiltonian,
        theta=np.asarray(theta, dtype=float),
        shots=self.shots,
        backend=self.backend,
    )
    return float(value)

run ¶

run(
    theta0: Sequence[float] | None = None,
    seed: int | None = None,
    compare_exact: bool | None = None,
    **optimizer_kwargs: Any,
) -> VQEResult

Optimise, and report against the exact answer when that is affordable.

Source code in src/qmlkit/algorithms/vqe.py

def run(
    self,
    theta0: Sequence[float] | None = None,
    seed: int | None = None,
    compare_exact: bool | None = None,
    **optimizer_kwargs: Any,
) -> VQEResult:
    """Optimise, and report against the exact answer when that is affordable."""
    start = (
        np.asarray(theta0, dtype=float)
        if theta0 is not None
        else self.ansatz.init("small", seed=seed)
    )
    fn = OPTIMIZERS[self.optimizer] if isinstance(self.optimizer, str) else self.optimizer
    # both gradient optimisers need the gradient injected; naming only one of them
    # here is what made optimizer="adam" raise a TypeError from every algorithm
    # that advertises it
    if fn in (_gradient_descent, _adam):
        optimizer_kwargs.setdefault("grad", self.gradient_of_energy)
    if fn is _spsa:
        optimizer_kwargs.setdefault("seed", seed)

    self.n_evaluations = 0
    theta, history = fn(self.energy, start, **optimizer_kwargs)

    # dense diagonalisation is exponential, so only offer it where it is cheap
    if compare_exact is None:
        compare_exact = self.n_qubits <= 12
    exact = exact_ground_energy(self.hamiltonian, self.n_qubits) if compare_exact else None

    return VQEResult(
        energy=float(history[-1]),
        theta=theta,
        history=list(history),
        exact=exact,
        n_evaluations=self.n_evaluations,
    )

`qmlkit.algorithms.adapt`¶

ADAPT-VQE: grow the ansatz one operator at a time, chosen by gradient magnitude.

The trap worth knowing. A molecular Hamiltonian conserves particle number, so any operator that does not has exactly zero gradient at Hartree–Fock — the generic pool grows an empty circuit and reports convergence. Use chemistry_operator_pool. This is physics, not a bug, and a test pins it.

adapt ¶

ADAPT-VQE — grow the ansatz instead of guessing it.

A fixed ansatz is a bet placed before you have seen the Hamiltonian. ADAPT-VQE (Grimsley et al. 2019) makes the circuit itself part of the optimisation: keep a pool of candidate generators, and at each iteration append the one whose gradient is largest, then re-optimise everything.

The gradient of appending :math:e^{-i\theta P/2} to the current state, evaluated at :math:\theta = 0, is

.. math:: \left.\frac{\partial E}{\partial\theta}\right|_0 = -i\langle\psi|[H, P]|\psi\rangle

so ranking the pool costs one commutator expectation per candidate — no re-training to find out which operator would have helped.

This is the algorithm that most depends on a circuit being data: growing an ansatz mid-optimisation is a list append here, not a rebuild.

AdaptVQE ¶

AdaptVQE(
    hamiltonian: Observable,
    n_qubits: int,
    pool: Sequence[PauliString] | None = None,
    optimizer: str | Optimizer = "gradient-descent",
    backend: BackendLike = None,
    reference: Sequence[int] | None = None,
)

Build the ansatz one operator at a time, largest gradient first.

Source code in src/qmlkit/algorithms/adapt.py

def __init__(
    self,
    hamiltonian: Observable,
    n_qubits: int,
    pool: Sequence[PauliString] | None = None,
    optimizer: str | Optimizer = "gradient-descent",
    backend: BackendLike = None,
    reference: Sequence[int] | None = None,
) -> None:
    self.hamiltonian = hamiltonian
    self.n_qubits = n_qubits
    self.pool = list(pool) if pool is not None else default_operator_pool(n_qubits)
    self.optimizer = optimizer
    self.backend = backend
    #: the starting state, as bits. Hartree-Fock plays this role in chemistry.
    self.reference = list(reference) if reference is not None else []

pauli_rotation ¶

pauli_rotation(
    qc: QCircuit, term: PauliString, angle: Any
) -> None

Emit :math:e^{-i\theta P/2} for an arbitrary Pauli string P.

The standard construction: rotate each wire into the Z basis, run a CX ladder to collect the parity onto one wire, apply a single rz, then undo both. Built from ordinary registered gates, so it runs on every backend.

Source code in src/qmlkit/algorithms/adapt.py

def pauli_rotation(qc: QCircuit, term: PauliString, angle: Any) -> None:
    r"""Emit :math:`e^{-i\theta P/2}` for an arbitrary Pauli string ``P``.

    The standard construction: rotate each wire into the Z basis, run a CX ladder to
    collect the parity onto one wire, apply a single ``rz``, then undo both. Built
    from ordinary registered gates, so it runs on every backend.
    """
    wires = [q for q, p in term.paulis if p != "I"]
    if not wires:
        return
    for qubit, pauli in term.paulis:
        if pauli == "X":
            qc.h(qubit)
        elif pauli == "Y":
            qc.sdg(qubit)
            qc.h(qubit)
    for a, b in zip(wires, wires[1:], strict=False):
        qc.cx(a, b)
    qc.rz(wires[-1], angle)
    for a, b in reversed(list(zip(wires, wires[1:], strict=False))):
        qc.cx(a, b)
    for qubit, pauli in term.paulis:
        if pauli == "X":
            qc.h(qubit)
        elif pauli == "Y":
            qc.h(qubit)
            qc.s(qubit)

default_operator_pool ¶

default_operator_pool(n_qubits: int) -> list[PauliString]

Single-qubit Y and neighbouring YZ generators.

Deliberately all imaginary-valued generators: those are the ones that move a real starting state, which is what a real-amplitude ground state needs. The pool is an argument, so a chemistry-flavoured (UCCSD-style) pool drops straight in.

Source code in src/qmlkit/algorithms/adapt.py

def default_operator_pool(n_qubits: int) -> list[PauliString]:
    """Single-qubit ``Y`` and neighbouring ``YZ`` generators.

    Deliberately all imaginary-valued generators: those are the ones that move a real
    starting state, which is what a real-amplitude ground state needs. The pool is an
    argument, so a chemistry-flavoured (UCCSD-style) pool drops straight in.
    """
    # Note these do *not* conserve particle number, so they are useless on a
    # molecular Hamiltonian -- see chemistry_operator_pool.
    pool = [PauliString(((q, "Y"),), 1.0) for q in range(n_qubits)]
    pool += [PauliString(((q, "Y"), (q + 1, "Z")), 1.0) for q in range(n_qubits - 1)]
    pool += [PauliString(((q, "Y"), (q + 1, "X")), 1.0) for q in range(n_qubits - 1)]
    return pool

chemistry_operator_pool ¶

chemistry_operator_pool(n_qubits: int) -> list[PauliString]

A particle-number-conserving pool, for molecular Hamiltonians.

:func:default_operator_pool is generic and wrong for chemistry: a molecular Hamiltonian commutes with the number operator, so any generator that changes particle number has exactly zero gradient at the Hartree-Fock state. Measured on H\ :sub:2: every operator in the default pool scores 0.00e+00, ADAPT correctly concludes nothing helps, and returns an empty circuit.

This is the qubit-ADAPT pool of Tang et al. (2021) — the individual Pauli strings appearing in single and double excitations, which under Jordan-Wigner carry an odd number of Y factors. On H\ :sub:2 the winning operator is the double excitation Y0 X1 X2 X3, and one parameter is enough to reach the exact ground state.

Source code in src/qmlkit/algorithms/adapt.py

def chemistry_operator_pool(n_qubits: int) -> list[PauliString]:
    r"""A particle-number-conserving pool, for molecular Hamiltonians.

    :func:`default_operator_pool` is generic and *wrong for chemistry*: a molecular
    Hamiltonian commutes with the number operator, so any generator that changes
    particle number has exactly zero gradient at the Hartree-Fock state. Measured on
    H\ :sub:`2`: every operator in the default pool scores ``0.00e+00``, ADAPT
    correctly concludes nothing helps, and returns an empty circuit.

    This is the qubit-ADAPT pool of Tang et al. (2021) — the individual Pauli strings
    appearing in single and double excitations, which under Jordan-Wigner carry an
    **odd number of Y factors**. On H\ :sub:`2` the winning operator is the double
    excitation ``Y0 X1 X2 X3``, and one parameter is enough to reach the exact ground
    state.
    """
    pool: list[PauliString] = []
    for p_, q_ in itertools.combinations(range(n_qubits), 2):  # singles
        pool.append(PauliString(((p_, "Y"), (q_, "X")), 1.0))
        pool.append(PauliString(((p_, "X"), (q_, "Y")), 1.0))
    for quad in itertools.combinations(range(n_qubits), 4):  # doubles
        for y_at in range(4):
            paulis = tuple((wire, "Y" if index == y_at else "X") for index, wire in enumerate(quad))
            pool.append(PauliString(paulis, 1.0))
    return pool

`qmlkit.algorithms.qaoa`¶

Quantum approximate optimisation.

Rotosolve is not valid here. QAOA's cost angle drives one rz per edge, and those do not compose into a single sinusoid — measured: five frequencies. Rotosolve's three-point fit then converges instantly to the wrong point and reports it as a result. Check with qmlkit.optim.supports_rotosolve before trusting it.

qaoa ¶

QAOA as a solver, not just an ansatz.

The ansatz has been in the zoo since Phase 3. What was missing is the part that makes it an algorithm: turn a combinatorial problem into a cost Hamiltonian, optimise the angles, then sample the state and read off a bitstring you can act on.

from qmlkit.algorithms import QAOA

edges = [(0, 1), (1, 2), (2, 3), (3, 0)]
result = QAOA(edges, p=2).run(seed=0)
print(result.bitstring, result.cut_value)

The cost Hamiltonian is an argument, so anything expressible as a Pauli sum — MaxCut, Max-2-SAT, a weighted graph, a portfolio constraint — is the same call.

QAOAResult `dataclass` ¶

QAOAResult(
    energy: float,
    theta: NDArray[Any],
    bitstring: str,
    probability: float,
    history: list[float] = list(),
    exact: float | None = None,
    top: list[tuple[str, float]] = list(),
)

The angles, and — more usefully — the bitstring they point at.

cut_value `property` ¶

cut_value: float

For a MaxCut cost, the number of edges the returned bitstring cuts.

approximation_ratio `property` ¶

approximation_ratio: float | None

Energy reached over the best possible, when the exact answer is known.

QAOA ¶

QAOA(
    problem: Observable | Sequence[tuple[int, int]],
    p: int = 1,
    n_qubits: int | None = None,
    mixer: str = "x",
    ansatz: Ansatz | None = None,
    optimizer: str | Optimizer = "gradient-descent",
    backend: BackendLike = None,
    shots: int | None = None,
)

Optimise QAOA angles, then sample a solution out of the state.

Parameters:

Name	Type	Description	Default
`problem`	`Observable \| Sequence[tuple[int, int]]`	Either an edge list (treated as MaxCut) or any cost `Observable`.	required
`p`	`int`	Rounds. Two angles per round regardless of problem size — which is the whole appeal, and also why more rounds is the only way to improve.	`1`
`mixer`	`str`	The structure. `mixer` is passed to :func:`qaoa_ansatz`; pass `ansatz` directly to use a warm-started or otherwise non-standard construction.	`'x'`
`ansatz`	`str`	The structure. `mixer` is passed to :func:`qaoa_ansatz`; pass `ansatz` directly to use a warm-started or otherwise non-standard construction.	`'x'`

Source code in src/qmlkit/algorithms/qaoa.py

def __init__(
    self,
    problem: Observable | Sequence[tuple[int, int]],
    p: int = 1,
    n_qubits: int | None = None,
    mixer: str = "x",
    ansatz: Ansatz | None = None,
    optimizer: str | Optimizer = "gradient-descent",
    backend: BackendLike = None,
    shots: int | None = None,
) -> None:
    if isinstance(problem, (list, tuple)) and not problem:
        raise ValueError(
            "QAOA needs a problem: an edge list for MaxCut, or a cost observable. "
            "An empty edge list defines nothing to optimise."
        )
    if isinstance(problem, (list, tuple)) and isinstance(problem[0], tuple):
        edges = [(int(a), int(b)) for a, b in problem]
        width = n_qubits or max(max(e) for e in edges) + 1
        self.cost: Observable = max_cut_hamiltonian(edges)
        self.edges: list[tuple[int, int]] | None = edges
    else:
        self.cost = problem  # type: ignore[assignment]
        support = observable_support(self.cost)
        width = n_qubits or (max(support) + 1 if support else 1)
        self.edges = None

    self.n_qubits = width
    self.p = p
    self.ansatz = ansatz or qaoa_ansatz(width, edges=self.edges, p=p, mixer=mixer)
    self.optimizer = optimizer
    self.backend = backend
    self.shots = shots
    self._spec = self.ansatz.build()
    self.n_evaluations = 0

distribution ¶

distribution(theta: ArrayLike) -> NDArray[Any]

Outcome probabilities of the optimised state.

Source code in src/qmlkit/algorithms/qaoa.py

def distribution(self, theta: ArrayLike) -> npt.NDArray[Any]:
    """Outcome probabilities of the optimised state."""
    return probabilities(self._spec.bind(np.asarray(theta, dtype=float)), backend=self.backend)

cut_size ¶

cut_size(bitstring: str) -> int

Edges cut by an assignment — the classical objective, computed classically.

Source code in src/qmlkit/algorithms/qaoa.py

def cut_size(self, bitstring: str) -> int:
    """Edges cut by an assignment — the classical objective, computed classically."""
    if self.edges is None:
        raise ValueError("cut_size only applies when the problem was given as edges")
    return sum(1 for a, b in self.edges if bitstring[a] != bitstring[b])

`qmlkit.algorithms.molecule` and `qmlkit.algorithms.chemistry`¶

Molecular Hamiltonians. Two routes, deliberately: from_integrals is the general one and takes PySCF or OpenFermion output for any molecule, while the built-in SCF handles s-orbital elements only. qmlkit is not a quantum chemistry package and does not try to become one.

molecule ¶

Molecular Hamiltonians for any molecule, by two routes.

Route one — bring your own integrals. This is the general one, and the one to reach for past a couple of light atoms. Anything that can produce one- and two-electron integrals in a molecular-orbital basis — PySCF, OpenFermion, Psi4 — hands them to :func:from_integrals and gets a qubit Hamiltonian back::

from pyscf import gto, scf, ao2mo
import numpy as np

mol = gto.M(atom="Li 0 0 0; H 0 0 1.6", basis="sto-3g")
mf = scf.RHF(mol).run()
c = mf.mo_coeff
h1 = c.T @ mf.get_hcore() @ c
h2 = ao2mo.restore(1, ao2mo.kernel(mol, c), c.shape[1])

hamiltonian, info = from_integrals(h1, h2, n_electrons=mol.nelectron,
                                   nuclear_repulsion=mol.energy_nuc())

That decoupling is deliberate. A quantum ML library should not also be a quantum chemistry package, and pretending otherwise would mean shipping a worse version of software that already exists.

Route two — the built-in SCF. For molecules built only from s-orbital atoms (hydrogen and helium in STO-3G) the integrals are computed here, with a real restricted Hartree–Fock loop rather than a symmetry shortcut. That covers the systems VQE is usually benchmarked on — H\ :sub:2, H\ :sub:3\ :sup:+, H\ :sub:4 chains and rings, HeH\ :sup:+ — at arbitrary geometry::

from qmlkit.algorithms import Molecule, molecular_hamiltonian

h4 = Molecule([("H", (0, 0, 0)), ("H", (0, 0, 0.9)),
               ("H", (0, 0, 1.8)), ("H", (0, 0, 2.7))])
hamiltonian, info = molecular_hamiltonian(h4)

Anything with p orbitals needs route one. That boundary is stated rather than papered over.

Molecule `dataclass` ¶

Molecule(
    atoms: list[tuple[str, tuple[float, float, float]]],
    charge: int = 0,
)

Atoms and where they are. Positions in angstrom.

n_orbitals `property` ¶

n_orbitals: int

One 1s function per atom, in this basis.

MolecularInfo `dataclass` ¶

MolecularInfo(
    n_qubits: int,
    n_electrons: int,
    n_orbitals: int,
    n_terms: int,
    nuclear_repulsion: float,
    hartree_fock_energy: float | None = None,
    hartree_fock_occupation: list[int] = list(),
    active_space: tuple[int, ...] | None = None,
)

Everything that went into the Hamiltonian, so the result can be audited.

hydrogen_chain ¶

hydrogen_chain(n: int, spacing: float = 0.74) -> Molecule

n hydrogens in a line — the standard scaling benchmark for VQE.

Source code in src/qmlkit/algorithms/molecule.py

def hydrogen_chain(n: int, spacing: float = 0.74) -> Molecule:
    """``n`` hydrogens in a line — the standard scaling benchmark for VQE."""
    return Molecule([("H", (0.0, 0.0, i * spacing)) for i in range(n)])

hydrogen_ring ¶

hydrogen_ring(n: int, radius: float = 1.0) -> Molecule

n hydrogens on a circle; frustrated, and harder than the chain.

Source code in src/qmlkit/algorithms/molecule.py

def hydrogen_ring(n: int, radius: float = 1.0) -> Molecule:
    """``n`` hydrogens on a circle; frustrated, and harder than the chain."""
    return Molecule(
        [
            ("H", (radius * np.cos(2 * np.pi * i / n), radius * np.sin(2 * np.pi * i / n), 0.0))
            for i in range(n)
        ]
    )

from_integrals ¶

from_integrals(
    one_body: NDArray[Any],
    two_body: NDArray[Any],
    n_electrons: int,
    nuclear_repulsion: float = 0.0,
    active_space: tuple[int, ...] | None = None,
    tol: float = 1e-10,
) -> tuple[PauliSum, MolecularInfo]

A qubit Hamiltonian from molecular-orbital integrals.

Parameters:

Name	Type	Description	Default
`one_body`	`NDArray[Any]`	`h[p, q]`, the one-electron integrals in the MO basis.	required
`two_body`	`NDArray[Any]`	`g[p, q, r, s]` in chemist notation `(pq\|rs)`, which is what PySCF's `ao2mo.restore(1, ...)` returns.	required
`active_space`	`tuple[int, ...] \| None`	Spatial orbitals to keep. Everything else is dropped, which is the usual way to fit a molecule onto a machine you actually have — `2 * len(active_space)` qubits instead of `2 * n_orbitals`.	`None`

Notes

This is the general entry point: it never asks where the integrals came from.

Source code in src/qmlkit/algorithms/molecule.py

def from_integrals(
    one_body: npt.NDArray[Any],
    two_body: npt.NDArray[Any],
    n_electrons: int,
    nuclear_repulsion: float = 0.0,
    active_space: tuple[int, ...] | None = None,
    tol: float = 1e-10,
) -> tuple[PauliSum, MolecularInfo]:
    r"""A qubit Hamiltonian from molecular-orbital integrals.

    Parameters
    ----------
    one_body
        ``h[p, q]``, the one-electron integrals in the MO basis.
    two_body
        ``g[p, q, r, s]`` in **chemist notation** ``(pq|rs)``, which is what PySCF's
        ``ao2mo.restore(1, ...)`` returns.
    active_space
        Spatial orbitals to keep. Everything else is dropped, which is the usual way
        to fit a molecule onto a machine you actually have — ``2 * len(active_space)``
        qubits instead of ``2 * n_orbitals``.

    Notes
    -----
    This is the general entry point: it never asks where the integrals came from.
    """
    h = np.asarray(one_body, dtype=float)
    g = np.asarray(two_body, dtype=float)
    if h.ndim != 2 or h.shape[0] != h.shape[1]:
        raise ValueError(f"one_body must be square, got shape {h.shape}")
    if g.shape != (h.shape[0],) * 4:
        raise ValueError(f"two_body must have shape {(h.shape[0],) * 4}, got {g.shape}")

    if active_space is not None:
        keep = list(active_space)
        h = h[np.ix_(keep, keep)]
        g = g[np.ix_(keep, keep, keep, keep)]

    n_spatial = h.shape[0]
    n_spin = 2 * n_spatial
    if n_spin > 12:
        raise ValueError(
            f"{n_spin} qubits means a {2**n_spin}-dimensional matrix. Use active_space "
            "to pick the orbitals that matter."
        )

    a = [_annihilator(p, n_spin) for p in range(n_spin)]
    adag = [x.conj().T for x in a]
    spin = [p % 2 for p in range(n_spin)]
    spatial = [p // 2 for p in range(n_spin)]

    matrix = np.zeros((2**n_spin, 2**n_spin), dtype=complex)
    for p, q in itertools.product(range(n_spin), repeat=2):
        if spin[p] == spin[q]:
            matrix += h[spatial[p], spatial[q]] * (adag[p] @ a[q])
    for p, q, r, s in itertools.product(range(n_spin), repeat=4):
        if spin[p] == spin[q] and spin[r] == spin[s]:
            matrix += (
                0.5
                * g[spatial[p], spatial[q], spatial[r], spatial[s]]
                * (adag[p] @ adag[r] @ a[s] @ a[q])
            )
    matrix += float(nuclear_repulsion) * np.eye(2**n_spin)

    terms = _decompose(matrix, n_spin, tol)
    occupation = [1 if i < n_electrons else 0 for i in range(n_spin)]
    info = MolecularInfo(
        n_qubits=n_spin,
        n_electrons=n_electrons,
        n_orbitals=n_spatial,
        n_terms=len(terms),
        nuclear_repulsion=float(nuclear_repulsion),
        hartree_fock_occupation=occupation,
        active_space=active_space,
    )
    return PauliSum(tuple(terms)), info

molecular_hamiltonian ¶

molecular_hamiltonian(
    molecule: Molecule,
    active_space: tuple[int, ...] | None = None,
    tol: float = 1e-10,
) -> tuple[PauliSum, MolecularInfo]

Compute the integrals here, then hand them to :func:from_integrals.

Restricted to s-orbital elements — see the module docstring for why, and for the route to take when that is not enough.

Source code in src/qmlkit/algorithms/molecule.py

def molecular_hamiltonian(
    molecule: Molecule,
    active_space: tuple[int, ...] | None = None,
    tol: float = 1e-10,
) -> tuple[PauliSum, MolecularInfo]:
    """Compute the integrals here, then hand them to :func:`from_integrals`.

    Restricted to s-orbital elements — see the module docstring for why, and for the
    route to take when that is not enough.
    """
    overlap, core, eri, repulsion = _ao_integrals(molecule)
    coefficients, electronic = _rhf(overlap, core, eri, molecule.n_electrons)
    h = coefficients.T @ core @ coefficients
    g = np.einsum(
        "pi,qj,rk,sl,pqrs->ijkl",
        coefficients,
        coefficients,
        coefficients,
        coefficients,
        eri,
        optimize=True,
    )
    hamiltonian, info = from_integrals(
        h, g, molecule.n_electrons, float(repulsion), active_space=active_space, tol=tol
    )
    info.hartree_fock_energy = float(electronic + repulsion)
    return hamiltonian, info

chemistry ¶

Molecular Hamiltonians, computed here rather than quoted.

VQE's canonical demonstration is the ground-state energy of H\ :sub:2, and most tutorials get the Hamiltonian by importing coefficients from a chemistry package or copying a table out of a paper. This module computes it: STO-3G integrals over Gaussian primitives, symmetry-adapted molecular orbitals, second quantisation, and a Jordan–Wigner map to four qubits.

That matters for a library whose whole argument is that you should be able to see the cost of what you run. It is also checkable — the curve below reproduces the published FCI/STO-3G result to five decimals:

>>> from qmlkit.algorithms.chemistry import h2_hamiltonian
>>> from qmlkit.algorithms import exact_ground_energy
>>> h, info = h2_hamiltonian(0.735)
>>> round(exact_ground_energy(h, 4), 5)
-1.13731

Minimal basis only, and two centres only. Anything larger wants PySCF or OpenFermion, and the point here is transparency rather than coverage.

h2_hamiltonian ¶

h2_hamiltonian(
    bond_length: float = 0.735, tol: float = 1e-10
) -> tuple[PauliSum, dict[str, Any]]

The H2 qubit Hamiltonian at a given bond length in angstrom.

Returns the observable and a dictionary of what went into it. The Pauli coefficients come from projecting the dense matrix, c_P = Tr(P H) / 2^n, which needs no symbolic algebra and is trivially checkable in the other direction with :func:~qmlkit.algorithms.hamiltonian_matrix.

Source code in src/qmlkit/algorithms/chemistry.py

def h2_hamiltonian(
    bond_length: float = 0.735, tol: float = 1e-10
) -> tuple[PauliSum, dict[str, Any]]:
    """The H2 qubit Hamiltonian at a given bond length in angstrom.

    Returns the observable and a dictionary of what went into it. The Pauli
    coefficients come from projecting the dense matrix, ``c_P = Tr(P H) / 2^n``,
    which needs no symbolic algebra and is trivially checkable in the other
    direction with :func:`~qmlkit.algorithms.hamiltonian_matrix`.
    """
    matrix, repulsion = _matrix(bond_length)
    terms: list[PauliString] = []
    for letters in itertools.product("IXYZ", repeat=4):
        operator = _kron(*[_PAULI[c] for c in letters])
        coefficient = float(np.real(np.trace(operator @ matrix)) / 16)
        if abs(coefficient) > tol:
            paulis = tuple((q, c) for q, c in enumerate(letters) if c != "I")
            terms.append(PauliString(paulis, coefficient))
    info = {
        "bond_length": bond_length,
        "n_qubits": 4,
        "n_terms": len(terms),
        "nuclear_repulsion": repulsion,
        "hartree_fock_occupation": [1, 1, 0, 0],  # both electrons in sigma_g
    }
    return PauliSum(tuple(terms)), info

h2_curve ¶

h2_curve(
    bond_lengths: NDArray[Any] | list[float],
) -> list[tuple[float, PauliSum]]

(bond_length, hamiltonian) pairs — the dissociation curve as input data.

Source code in src/qmlkit/algorithms/chemistry.py

def h2_curve(bond_lengths: npt.NDArray[Any] | list[float]) -> list[tuple[float, PauliSum]]:
    """``(bond_length, hamiltonian)`` pairs — the dissociation curve as input data."""
    return [(float(r), h2_hamiltonian(float(r))[0]) for r in bond_lengths]

`qmlkit.algorithms.hamiltonians`¶

Standard model Hamiltonians — Ising, Heisenberg, and the rest — as PauliSums.

hamiltonians ¶

Hamiltonians to hand to VQE, and an exact answer to check it against.

A Hamiltonian here is just a :class:~qmlkit.core.observables.PauliSum — the same type an expectation value takes — so nothing new has to learn about it. These are constructors, not a new class hierarchy.

:func:exact_ground_energy diagonalises the dense matrix. That is exponential and useless past ~14 qubits, which is exactly the point: it is the oracle a variational result gets checked against on small systems, not a method to compete with.