cpp_hf — Hartree–Fock on k‑grids (C++ reimplementation of jax_hf)

cpp_hf is a C++ reimplementation of jax_hf: two Hartree–Fock solvers on 2D k‑meshes that share the common kernel, config, solver, continuation, and result surface, with the JAX dependency replaced by a pybind11 extension built on FFTW + Eigen.

• Direct minimisation (default) — preconditioned Riemannian CG on Stiefel × capped simplex with Cayley retraction, or trust-region Newton. Trial free energies are evaluated explicitly, so backtracking and Newton Hessian-vector products may require multiple Fock builds per outer iteration.
• Self-consistent field — safeguarded full-density-residual DIIS and Anderson acceleration, plus the legacy linear, commutator-DIIS and optimal-damping methods. rdiis is the recommended fixed-point solver for finite-temperature workloads.

Both solvers run in double precision throughout (complex128 / float64).

Install

Requires FFTW (double precision: fftw3) and a C++17 compiler.

# macOS
brew install fftw eigen
pip install -e .

# Linux
sudo apt install libfftw3-dev libeigen3-dev
pip install -e .

Minimal example

import numpy as np
import cpp_hf

# Build a HartreeFockKernel: precomputes the FFT of the interaction kernel,
# the Hartree matrix, etc., ready for the solver.
kernel = cpp_hf.HartreeFockKernel(
    weights=weights,          # (nk1, nk2) k-point weights
    hamiltonian=hamiltonian,  # (nk1, nk2, nb, nb) single-particle Hamiltonian
    coulomb_q=coulomb_q,      # (nk1, nk2, 1, 1) scalar or (nk1, nk2, nb, nb) layer-resolved
    T=0.1,
    include_hartree=False,    # set True for Hartree; also pass reference_density + hartree_matrix
    include_exchange=True,
)

# Safeguarded acceleration with explicit stationarity tolerances.
result = cpp_hf.solve_scf(
    kernel, P0=np.zeros_like(hamiltonian), n_electrons=N,
    config=cpp_hf.SCFConfig(acceleration="rdiis", diis_size=8,
                           density_tol=1e-6, comm_tol=1e-6, max_iter=500),
)
print(result.energy, result.converged, result.iterations)
print(result.history["fock_builds"])
# result.density_matrix, result.fock_matrix, result.chemical_potential

# Or use the hybrid with Newton recovery and the direct-solver result surface.
result_dm = cpp_hf.solve(
    kernel, P0=np.zeros_like(hamiltonian), n_electrons=N,
    config=cpp_hf.SolverConfig(optimizer="auto", tol_grad=1e-6, max_iter=500),
)
# result_dm.density, result_dm.fock, result_dm.Q, result_dm.p, result_dm.mu

For a local minimum with a stability check and bounded saddle recovery:

minimum = cpp_hf.solve_minimum(kernel, np.zeros_like(hamiltonian), N)
print(minimum.converged, minimum.solution.energy, minimum.escapes)
print(minimum.fock_builds)  # includes stability and recovery cost

A converged residual alone does not establish that a solution is a minimum. solve_minimum starts with safeguarded RDIIS, checks negative curvature, and tries finite feasible descent perturbations when needed. Escape trials first use a fixed-particle Fermi-density map along the interaction response of the unstable mode; this preserves the tangent near saturated occupations. Projected density perturbations remain a fallback. Both require actual free-energy descent. For phase comparisons, compare free energies across seeds. The standalone diagnostic is cpp_hf.stability_analysis(kernel, result.density_matrix, N). These are local checks within the supplied basis; they do not certify a global minimum. Existing calls to solve without a config still select the CG optimizer, but its stopping test and objective changed in 1.2.0 (see CHANGELOG.md).

Architecture

Layer	Where it lives	What it does
C++ core	cpp/include/cpp_hf/*.hpp + cpp/cpp_hf_native.cpp	FFT‑based exchange (selfenergy_fft), batched Hermitian eigendecomposition, contact‑term Fock construction, SCF main loop, direct‑minimisation main loop (preconditioned Riemannian CG with spectral Cayley line search), k‑grid resampling — all in double precision.
Native extension	cpp_hf._native (compiled .so inside the package)	pybind11 wrapper exposing the C++ entry points; the GIL is released for the duration of every solver call.
Python surface	src/cpp_hf/*.py	Public dataclasses (SolverConfig, SCFConfig, SolveResult, SCFResult, ContinuationResult), the HartreeFockKernel constructor (input validation + kernel precomputation), the symmetry projector framework (passed to the C++ solver as a Python callback), and the continuation driver (composes two C++ solver calls + resample_kgrid).

The Python surface is intentionally thin: it validates inputs, packs them into the dict shape the C++ binding expects, then hands off. The expensive Fock and iteration kernels run in C++; stability and escape orchestration use Python/SciPy. Importing the package without the compiled extension raises a clear RuntimeError from native_required().

Public API

Name	Purpose
HartreeFockKernel	Problem definition + precomputed arrays
solve (alias solve_direct_minimization), SolverConfig, SolveResult	Primary solver
solve_scf, SCFConfig, SCFResult	Safeguarded full-residual and legacy SCF methods
solve_minimum, MinimumResult	RDIIS with local stability checks and bounded saddle recovery
stability_analysis, StabilityResult	Matrix-free local stability test of an unconstrained finite-T solution
build_fock, hf_energy, free_energy, occupation_entropy	HF objective building blocks
solve_continuation, ContinuationResult, resample_kgrid	Coarse → fine multigrid driver + k‑grid resampler
cpp_hf.symmetry.make_project_fn	Symmetry projector builder (unitary / spatial / time‑reversal)
cpp_hf.linalg.eigh	Block‑diagonal Hermitian eigh with optional structure check

The common API closely follows jax_hf, and ordinary built-in-kernel scripts can generally switch by changing the import. Backend-specific features remain explicit: cpp_hf does not accept JAX tracing/differentiation or the fock_build_fn= callback, and typed config objects belong to the backend that created them. Use a plain config dictionary in higher-level code that needs to route settings across both backends.

Coarse → fine continuation

from cpp_hf import HartreeFockKernel, SolverConfig, SCFConfig, solve_continuation

coarse = HartreeFockKernel(weights_c, h_c, Vq_c, T=0.1)
fine   = HartreeFockKernel(weights_f, h_f, Vq_f, T=0.1)

result = solve_continuation(
    coarse, fine, P0_coarse=np.zeros_like(h_c),
    n_electrons_coarse=N, n_electrons_fine=N,
    coarse_config=SCFConfig(max_iter=50, mixing=0.5),
    fine_config=SolverConfig(max_iter=200, tol_E=1e-8),
)
# result.coarse, result.fine, result.P0_fine

Tests

pip install -e .
python -m pytest tests/

The suite is a port of jax_hf/tests/; it covers the chemical‑potential solver, block‑diagonal eigh and self‑energy, the symmetry projectors, direct minimisation (basic convergence, multi‑k, contact terms, edge cases, Cayley spectral identities), SCF, continuation, and package‑import contracts. The bilayer regression test (test_bilayer_regression.py, ported from jax_hf) is not included by default since it requires contimod.

License

GPLv2+ — see LICENSE.

Acknowledgments

• FFTW — GPLv2+ (www.fftw.org)
• Eigen — MPL2
• pybind11 — BSD‑style