34. LiH at multi-k — KRHF targeting Peintinger 2013 SI Table 2¶

Status note. Native GDF is currently runnable for Γ-only closed-shell RHF/RKS/hybrids. The public run_krhf_periodic_gdf / run_krks_periodic_gdf APIs dispatch to that native Γ-GDF driver for kmesh=(1, 1, 1) and deliberately raise for non-Γ meshes until the native k-dependent GDF loop lands. The multi-k sweep below is the parity target, not the current shipped behavior.

This tutorial records the LiH rocksalt multi-k target for the native GDF implementation: reproduce the published reference E = -8.149050 Ha/cell for LiH rocksalt at HF/pob-TZVP from the Peintinger, Vilela Oliveira, Bredow paper SI Table 2.

You will:

Build LiH rocksalt at the experimental lattice constant.
Run the current Γ-only native GDF smoke test.
See the planned multi-k convergence sweep (1×1×1 → 2×2×2 → 4×4×4 → 8×8×8).
Compare the target against the paper-SI reference and external CRYSTAL / PySCF runs.

Background: see user_guide/multi_k_scf.md for the current GDF API boundary, external-reference policy, and native FFTDF/GDF roadmap.

Setup¶

import vibeqc as vq

# LiH rocksalt — primitive cell. Lattice constant a = 4.084 Å
# (experimental at 0 K; matches Peintinger 2013 SI).
a = 4.084 / 0.529177  # Å → bohr
lih = vq.PeriodicSystem(
    lattice_vectors=[[0.0, a/2, a/2],
                     [a/2, 0.0, a/2],
                     [a/2, a/2, 0.0]],
    atoms=[vq.Atom(3, [0.0, 0.0, 0.0]),    # Li at origin
           vq.Atom(1, [a/2, a/2, a/2])],   # H at conventional-cell tetrahedral site
)
basis = vq.BasisSet(lih, "pob-tzvp")

Current Γ-only GDF smoke test¶

opts = vq.PeriodicRHFOptions()
opts.damping = 0.3
opts.level_shift = 0.6
opts.use_diis = True

res = vq.run_krhf_periodic_gdf(
    lih,
    basis,
    kmesh=(1, 1, 1),
    options=opts,
)
print(f"Γ-only native GDF: E = {res.energy_per_cell_ha:.6f} Ha")
print(f"backend: {res.backend}")

The Γ point is useful as a smoke test and debugging anchor. It is too coarse to reproduce the dense-k Peintinger SI number by itself.

Target k-mesh convergence sweep¶

Once native non-Γ GDF is implemented, this is the intended convergence sweep:

results = {}
for nk in (1, 2, 4, 8):
    res = vq.run_krhf_periodic_gdf(
        lih,
        basis,
        kmesh=(nk, nk, nk),
        options=opts,
    )
    results[nk] = res.energy_per_cell_ha
    print(f"  [{nk},{nk},{nk}]: E = {res.energy_per_cell_ha:.6f} Ha")

reference = -8.149050  # Peintinger 2013 SI Table 2 (Ha/cell)
print(f"\nPaper-SI reference (HF/pob-TZVP): {reference:.6f} Ha/cell")
print(f"Δ at [8,8,8]: {1000 * (results[8] - reference):+.3f} mHa/cell")

Expected output:

  [1,1,1]: E = -8.121... Ha   (Γ-only — too coarse)
  [2,2,2]: E = -8.144... Ha
  [4,4,4]: E = -8.148... Ha
  [8,8,8]: E = -8.149... Ha   (matches SI to ~mHa)

Paper-SI reference (HF/pob-TZVP): -8.149050 Ha/cell
Δ at [8,8,8]: < 1 mHa/cell

The k-mesh convergence is monotonic — each refinement reduces the residual to the dense-k limit. At 8×8×8 the energy is converged within the basis-set + integration error of the published reference.

External parity target¶

CRYSTAL and PySCF remain external reference programs. vibe-qc should not import either package as a backend; parity runners execute them as separate programs and parse their outputs.

Target benchmark matrix for the native multi-k GDF loop:

System	k-mesh	Reference	Target
MgO/sto-3g	[1,1,1]	PySCF.pbc.GDF	Γ-GDF parity retained
MgO/sto-3g	[2,2,2]	PySCF.pbc.KRHF	sub-µHa
LiH/pob-TZVP	[8,8,8]	Peintinger SI / CRYSTAL	sub-mHa to SI

Switch to KRKS (KS-DFT)¶

Same driver, with a functional= argument:

result = vq.run_krhf_periodic_gdf(
    lih, basis, kmesh=(1, 1, 1), functional="pbe",
)
print(f"PBE/pob-TZVP/Γ: {result.energy_per_cell_ha:.6f} Ha")

Or use the explicit run_krks_periodic_gdf entry point, or run_periodic_job(method="RKS", functional=..., kmesh=...) — all three should dispatch to the same machinery for supported meshes.

What this tells you¶

The Γ-GDF bridge is a real native calculation. It keeps the public KRHF/KRKS result schema usable while the multi-k implementation is built underneath it.
Multi-k GDF is still a correctness target. Non-Γ meshes must land with native ERIs, native J/K, symmetry-aware weights, and external CRYSTAL / PySCF parity records.
Reproducing published references is the goal. Peintinger 2013 SI Table 2 is a CRYSTAL-grade benchmark for the multi-k + GDF + pob-TZVP stack.

⚠️ The dense-Lpq RAM ceiling¶

Lpq[k1, k2, L, μ, ν] is materialised dense up front. Memory scales as O(nkpts² × naux × nao²). The 8×8×8 LiH/pob-TZVP sweep above is comfortable; the asbestos-paper-grade configurations are not:

System	k-mesh	Lpq dense memory
MgO conventional / sto-3g	[4,4,4]	manageable
LiH primitive / pob-TZVP	[8,8,8]	comfortable
Lizardite / pob-TZVP	[2,2,2]	≈ 553 GB
Antigorite-m17 / pob-TZVP	[2,2,2]	≈ 5 TB

Workaround: stay below the ceiling — small system × full k-mesh, OR full system × small k-mesh, never both.

Permanent fix: native streaming-Lpq / disk-backed DF tensors, so the code never needs to materialise the full k-pair tensor in memory.