Native Hartree-Fock (dense_evolution.native_hf)

A from-scratch, JAX-vectorized ab-initio Hartree-Fock engine — overlap, kinetic, nuclear-attraction, and electron-repulsion integrals over s/p Cartesian Gaussian shells via the Obara-Saika recursion (Obara & Saika, J. Chem. Phys. 84, 3963, 1986), each shell-pair/quartet batched with jax.lax.scan/jax.vmap and compiled with jax.jit instead of looping in the Python interpreter.

It exists because PennyLane's own differentiable Hartree-Fock solver (qml.qchem, method="dhf") builds these same integrals through a Python-level loop wrapped in its autograd-tracing numpy layer — correct, but profiled directly at 482 of 483 total seconds for Si2/STO-3G, almost entirely per-scalar-op tracer overhead rather than real FLOPs. This module only replaces the integral/SCF stage; the converged result still goes to PennyLane's own fermionic_observable + jordan_wigner for the qubit mapping, since that stage is already fast (under 2 seconds) and well-tested. Basis-set parameters come from the basis_set_exchange package, so any element it has STO-3G data for is reachable, not just PennyLane's bundled H–Ne table. An element needing d-orbitals or higher (e.g. Fe) fails with a clear NotImplementedError naming the real limitation, not a silent wrong energy for an incomplete basis.

Design and algorithm structure were informed by studying lowdanie/hartree-fock-solver ("slaterform", Apache-2.0) as a reference for structuring the Obara-Saika recursion with jax.lax.scan, and by PennyLane's own white paper (Delgado et al., "Differentiable quantum computational chemistry with PennyLane", arXiv:2111.09967) — no source code from either project is copied here. Verified element-wise against an independent JAX Hartree-Fock implementation (slaterform) to machine precision on individual integrals and to 10 significant figures on Si2/STO-3G's full SCF energy.

dashboard_core.hamiltonians calls this engine automatically — bridge.build_qubit_hamiltonian — whenever a requested molecule uses an element outside PennyLane's own STO-3G table; existing molecules (H2/HeH+/H3+/LiH/H2O) are unaffected and keep using PennyLane's dhf pipeline directly. See Dashboard Core — Hamiltonians for the dispatch logic and the Si2 catalog entry this engine backs.

bridge

Bridge from our native Hartree-Fock result to a PennyLane qubit Hamiltonian.

The expensive part (Hartree-Fock: integrals + SCF, everything in this package) is entirely ours. Second quantization and the Jordan-Wigner mapping are cheap (profiled at under 2 seconds even for Si2 -- see dense_evolution/native_hf/init.py's module docstring) and PennyLane already does them well via public functions, so we call those directly instead of reimplementing them.

build_qubit_hamiltonian

build_qubit_hamiltonian(
    atomic_numbers: list[int],
    geometry_angstrom: ndarray,
    n_electrons: int,
    active_electrons: int = None,
    active_orbitals: int = None,
    basis_name: str = "sto-3g",
    cutoff: float = 1e-12,
) -> tuple[qml.Hamiltonian, int, HFResult]

Runs native Hartree-Fock, then hands the result to PennyLane for second quantization + Jordan-Wigner mapping.

Returns:

Type	Description
Hamiltonian	tuple[qml.Hamiltonian, int, HFResult]: the mapped qubit
int	Hamiltonian, the qubit count, and the native HFResult -- the
HFResult	latter useful for e.g. reporting the SCF energy alongside the
tuple[Hamiltonian, int, HFResult]	post-mapping ground-state energy.

Source code in dense_evolution/native_hf/bridge.py

def build_qubit_hamiltonian(
    atomic_numbers: list[int],
    geometry_angstrom: np.ndarray,
    n_electrons: int,
    active_electrons: int = None,
    active_orbitals: int = None,
    basis_name: str = "sto-3g",
    cutoff: float = 1e-12,
) -> tuple["qml.Hamiltonian", int, HFResult]:
    """Runs native Hartree-Fock, then hands the result to PennyLane for
    second quantization + Jordan-Wigner mapping.

    Returns:
        tuple[qml.Hamiltonian, int, HFResult]: the mapped qubit
        Hamiltonian, the qubit count, and the native HFResult -- the
        latter useful for e.g. reporting the SCF energy alongside the
        post-mapping ground-state energy.
    """
    geometry_bohr = np.asarray(geometry_angstrom) * _BOHR_PER_ANGSTROM
    shells = build_molecule_shells(atomic_numbers, geometry_bohr, basis_name)

    S = build_overlap_matrix(shells)
    H_core = build_core_hamiltonian(shells, [float(z) for z in atomic_numbers], geometry_bohr)
    repulsion = build_repulsion_tensor(shells)

    hf_result = run_scf(S, H_core, repulsion, n_electrons, [float(z) for z in atomic_numbers], geometry_bohr)

    one, two = _ao_to_mo(H_core, repulsion, hf_result.orbital_coefficients)
    n_orbitals = one.shape[0]

    if active_electrons is None and active_orbitals is None:
        core_idx, active_idx = [], list(range(n_orbitals))
    else:
        core_idx, active_idx = qml.qchem.active_space(
            n_electrons, n_orbitals, active_electrons=active_electrons, active_orbitals=active_orbitals
        )

    core_constant, one_active, two_active = _apply_active_space(
        hf_result.nuclear_repulsion_energy, one, two, core_idx, active_idx
    )

    fermi_op = _pl_observable.fermionic_observable(
        np.array([core_constant]), one_active, two_active, cutoff
    )
    qubit_hamiltonian = qml.jordan_wigner(fermi_op)
    n_qubits = 2 * len(active_idx)

    return qubit_hamiltonian, n_qubits, hf_result

scf

Restricted Hartree-Fock self-consistent field loop (Roothaan-Hall).

Standard textbook algorithm (e.g. Szabo & Ostlund, "Modern Quantum Chemistry", ch. 3): orthogonalize the AO basis via S^(-1/2), build the Fock matrix F = H_core + 2J - K from the current density, diagonalize in the orthogonal basis, form a new density, repeat to convergence. This part is genuinely simple compared to the integral evaluation and doesn't need vectorizing -- a closed-shell molecule's SCF loop is a few dozen matrix multiplies on an N x N matrix where N is a few tens at most for STO-3G, nowhere near where PennyLane's implementation loses its time (which is entirely in building H_core/repulsion tensor, done once in assembly.py, not in this loop).

BUG FOUND (Si2, minimal 4-electron/4-orbital active space, R=2.184 A): plain (undamped) density substitution never converged for this system -- 100/100 iterations, still oscillating -- because two pairs of orbitals near the active-space boundary are numerically degenerate (HOMO-1/HOMO and LUMO/LUMO+1 each split by <1e-9 Ha), so each iteration flips which member of a near-tied pair gets occupied, and the density never settles. Confirmed this is a real oscillation, not just slow convergence: three separate machines/runs of the undamped loop each hit the iteration cap at a DIFFERENT total energy (-571.63, -570.69, -571.02 Ha), all physically meaningless artifacts of whatever step the loop happened to be on. Fixed with standard linear density damping (P_next = alphaP_new + (1-alpha)P_old) -- textbook remedy for exactly this oscillation failure mode (Szabo & Ostlund ch. 3.4.9). Verified the damped loop converges to the SAME energy (-570.874032094871 Ha, agreeing to 10 significant figures) across alpha in {0.1, 0.2, 0.3, 0.5, 0.7} -- the fix's correctness doesn't hinge on the specific alpha chosen, only that damping is applied at all.

basis

Loading contracted basis-set shells for a molecule.

Basis-set parameters (exponents, contraction coefficients) are fetched from the Basis Set Exchange (the basis_set_exchange PyPI package, data-only, BSD-licensed -- https://www.basissetexchange.org), so this module works for any element the basis has data for, not just the handful PennyLane's own bundled STO-3G table covers (which stops at Ne and is why Silicon needed a hand-written patch earlier in this project).

The coefficients BSE reports are for unnormalized primitives, so each primitive's contraction coefficient must be rescaled by its own L2 norm before it can be summed against another shell's primitives. We get that norm for free by reusing our own overlap_3d on the primitive against itself: N = 1/sqrt().

ContractedShell dataclass

ContractedShell(
    atom_index: int,
    center: Array,
    degree: int,
    exponents: Array,
    coefficients: Array,
)

One angular-momentum shell (s, p, ...) of a contracted GTO.

build_molecule_shells

build_molecule_shells(
    atomic_numbers: list[int],
    geometry_bohr: ndarray,
    basis_name: str,
) -> list[ContractedShell]

geometry_bohr: shape (n_atoms, 3), atomic units.

Source code in dense_evolution/native_hf/basis.py

def build_molecule_shells(atomic_numbers: list[int], geometry_bohr: np.ndarray, basis_name: str) -> list[ContractedShell]:
    """geometry_bohr: shape (n_atoms, 3), atomic units."""
    shells = []
    for i, (z, r) in enumerate(zip(atomic_numbers, geometry_bohr)):
        shells.extend(load_element_shells(basis_name, z, jnp.asarray(r), i))
    return shells

---

See also: Dashboard Core — Hamiltonians for the production entry point (MOLECULE_CATALOG's Si2 entry), and dashboard_core.vqe for the ansatz circuits optimized against Hamiltonians this engine can build.