M3b C++ Kernel — Handover for new chat

Date: 2026-05-22 | Status: Fully specified, ready for implementation

Design constraint: additive, keyword-activated

The new kernel must run alongside existing code — do NOT modify build_jk_gamma_molecular_limit or build_jk_gamma_molecular_limit_explicit. Add a new function with its own name, own pybind11 binding, own Python export. The existing J/K builders remain the default; the new kernel is called explicitly by the Python symmetry-orchestration layer.

What this is

M3b is the blocked milestone for two-electron Fock build symmetry reduction. The C++ substrate (explicit-cell kernels) exists, but the Python orbit-reconstruction layer fails because per-cell J(g)/K(g) blocks cannot be reconstructed via Wigner-D: J(g) involves an internal sum over g_λ that is lost when cells are restricted.

The solution is a new C++ kernel that works at the cell-PAIR level rather than the cell level.

What needs to be built

New C++ function

// In cpp/src/periodic_fock.cpp, declared in cpp/include/vibeqc/periodic_fock.hpp

struct PairJKContribution {
    int c_g;                     // index into the cell list
    int c_p;                     // index into the cell list
    Eigen::MatrixXd J_contrib;   // (nbf, nbf) J contribution from this pair
    Eigen::MatrixXd K_contrib;   // (nbf, nbf) K contribution from this pair
};

// Evaluate J and K contributions for each (c_g, c_p) pair in `pairs`.
// `cells` is the full cell list; `pairs[i] = {c_g, c_p}` indexes into it.
// `P_gamma` is the Γ-only density matrix (nbf, nbf).
// `omega` selects Coulomb kernel (0 = full 1/r, >0 = erfc-screened).
std::vector<PairJKContribution> build_jk_pair_contributions(
    const BasisSet& basis,
    const PeriodicSystem& system,
    const std::vector<LatticeCell>& cells,
    const std::vector<std::pair<int,int>>& pairs,
    const LatticeSumOptions& opts,
    const Eigen::MatrixXd& P_gamma,
    double omega = 0.0);

Algorithm (copy from build_jk_gamma_molecular_limit_explicit, lines 55-275 of periodic_fock.cpp)

For each pair (c_g, c_p):

  1. Get shells_g = shells_at[c_g], shells_p = shells_at[c_p]

  2. Pre-compute Schwarz screening as in the existing kernel

  3. For each shell quartet (s1, s2, s3, s4):

    • Compute ERI via libint

    • Compute J contribution: J(μ,ν) += P(λ,σ) * (μ0 νg | λp σp)

    • Compute K contribution: K(μ,ν) += P(λ,σ) * (μ0 λp | νg σp)

  4. Store J_contrib and K_contrib in the output vector

The existing build_jk_gamma_molecular_limit_explicit kernel does EXACTLY this, but it SUMS over all pairs into a single J and K. The new kernel just needs to NOT sum — store each pair’s contribution separately.

Key implementation note

The existing kernel already has the inner loop that computes J and K contributions for each pair. It uses thread-local accumulators. For the new kernel:

  • Keep the #pragma omp parallel for over pairs

  • Instead of accumulating into shared Jm_tls/Km_tls, write each pair’s result directly into a pre-allocated std::vector<PairJKContribution>

  • Each thread writes to a distinct slot (pair index), no locking needed

pybind11 binding

py::class_<PairJKContribution>(m, "PairJKContribution")
    .def_readonly("c_g", &PairJKContribution::c_g)
    .def_readonly("c_p", &PairJKContribution::c_p)
    .def_readonly("J_contrib", &PairJKContribution::J_contrib)
    .def_readonly("K_contrib", &PairJKContribution::K_contrib);

m.def("build_jk_pair_contributions", &build_jk_pair_contributions,
      py::arg("basis"), py::arg("system"), py::arg("cells"),
      py::arg("pairs"), py::arg("options"), py::arg("density"),
      py::arg("omega") = 0.0,
      py::call_guard<py::gil_scoped_release>(),
      "Per-pair J/K contributions for symmetry-reduced Fock build (M3b).");

Python orchestration (to be built after C++ kernel)

The new kernel is called explicitly by the Python symmetry layer — it does NOT replace build_jk_gamma_molecular_limit. In the GDF SCF driver, the call site would be guarded by a keyword like symmetry_reduce_fock=True.

# In symmetry_integrals_reduced.py or a new symmetry_fock_reduced.py module

def compute_jk_gamma_reduced(basis, system, opts, D_gamma, operations):
    """Build J(Γ) and K(Γ) with symmetry-reduced cell pairs."""
    from ._vibeqc_core import build_jk_pair_contributions, direct_lattice_cells
    from .symmetry_lattice_c import identify_atom_pair_orbits, _cell_tuple
    from .symmetry_integrals import symmorphic_operations
    from .symmetry_scf import build_ao_permutation_cache

    sym_ops = symmorphic_operations(operations)
    full_cells = list(direct_lattice_cells(system, opts.cutoff_bohr))
    orbits = identify_atom_pair_orbits(full_cells, sym_ops, system, require_closed=False)
    P_cache = build_ao_permutation_cache(system, basis, sym_ops)

    # Build representative pair list from orbit representatives
    # Each orbit rep is (a, b, h) — we need (c_g=h, c_p=any) pairs.
    # For Γ-only, the density is at g=0, so g_σ = g_λ, and cell pairs
    # are (g, g_λ) where both are the orbit rep cells.
    
    nbf = basis.nbasis
    J_gamma = np.zeros((nbf, nbf))
    K_gamma = np.zeros((nbf, nbf))

    # For each orbit of cell pairs, compute at rep, scatter via Wigner-D
    for orb in orbits.orbits:
        # ... compute representative pair contributions via C++ kernel
        # ... for each member, apply D_a @ J_rep @ D_b^T and accumulate
        pass

    return J_gamma, K_gamma

Validation

After the C++ kernel is built, verify:

  1. Calling with all pairs gives the same J and K as build_jk_gamma_molecular_limit

  2. Calling with a subset of pairs gives partial results

  3. The Python reconstruction produces J/K matching the full computation for a SINGLE-ATOM origin-fixed cell (Mg primitive, Pm-3m)

References

  • Existing kernel: cpp/src/periodic_fock.cpp lines 55-275 (build_jk_gamma_molecular_limit_explicit)

  • Header: cpp/include/vibeqc/periodic_fock.hpp

  • Bindings: cpp/src/bindings.cpp (search build_jk_gamma_molecular_limit_explicit)

  • Python exports: python/vibeqc/__init__.py

  • Status file: docs/symmetry_status.md