± supermatrix and half-transformed AO stores

Persisted ± (plus/minus) supermatrix store for the exact AO integrals.

The two symmetric-pair-space supermatrices

Vs[tri(μν),tri(ρσ)] = ⟨μν|ρσ⟩ + ⟨νμ|ρσ⟩      (row/col diagonal μ=ν carries ×2)
Va[tri(μν),tri(ρσ)] = ⟨μν|ρσ⟩ − ⟨νμ|ρσ⟩      (row/col diagonals = 0)

(exactly the QMTensors.calc_tri_sym_antisym! convention) are each symmetric for real integrals and Hermitian for complex ones. Only the lower block-triangle is stored, as dense column panels — so kext consumer is a plain zero-copy GEMM ('N' + adjoint flag), never a packed-format special case. Total on disk ≈ n⁴/4.

Dimensions: μ,ν,ρ,σ are AO indices running 1:nao; tri(μν) packs an AO pair μ≤ν into 1:npp with npp = nao·(nao+1)/2. So each supermatrix is npp × npp in pair space, and the n⁴/4 size estimates use n = nao (npp² ≈ n⁴/4).

Files: ao_pm_s, ao_pm_a (1-D panel concatenations, ket-2-orbital σ-aligned blocks) and ao_pm_meta (the σ block breakpoints, so reader and writer never disagree on the layout).

THE CONJUGATION RULE: pair-internal swaps (μ↔ν, ρ↔σ) are the ± signs and NEVER conjugate; the block mirror (bra-pair↔ket-pair, the unstored triangle) conjugates with no sign. In code: mirror-role GEMMs use adjoint, mirror-role scatters use conj — both no-ops for real T (one code path).

Entry points: build with pm_from_joint! (or IntegralTools.pm_integrals! for the fused, no-joint build); open a built store with open_pm_store; contract with pm_matmul! or the FockFactory.ao_JK!(…, ::PMSupermatrices, …) methods. See PMSupermatrices for the handle type, its fields, and usage examples.

BothSlab{TF}

One ket-pair slab yielded by eachslab(pm; roles=:both): the reconstructed ⟨μν|ρσ⟩ (s.w.G) and its transpose (s.w.Gt) over the band (s.lo, s.hi], for ket pair (s.ρ, s.σ). Aliases the iterator scratch — valid only until the next iteration. Contract with pm_bra_half!.

HTStore{T}

Open handle to a half-transformed AO-integral store built by pm_half_trans: the mmapped map[m·nao, 2, npp] with map[(i,ν), 1, tri(ρσ)] = Σ_μ ⟨μν|ρσ⟩ C[μ,i] (ket order (ρ,σ)) and map[(i,ν), 2, tri(ρσ)] = Σ_μ ⟨μν|σρ⟩ C[μ,i] (ket order (σ,ρ)). Read dense σ-columns with ht_column!; the pair-packed layout never leaks to consumers.

PMElementReader(pm) -> reader

O(1) random access to the ± supermatrices: pm_vs_va(reader, B, K) returns (V_s[B,K], V_a[B,K]) for packed pair indices B (bra) and K (ket), resolving the block-triangle storage through the hermitian mirror when B lies above the ket panel's row floor. For consumers that genuinely need scattered elements (e.g. a pivoted Cholesky of the integral matrix) — sweeps should use eachslab.

PMSlab{TF}

One reconstructed ket-pair slab yielded by eachslab: Gslab(s)[μ,ν] = ⟨μν|ρσ⟩ for the ket pair (s.ρ, s.σ), with the native/mirror L-band lower bounds s.native_lo/s.mirror_lo (band = (lo, s.nao]). The slab aliases the iterator's single reusable scratch — it is valid ONLY until the next iteration; do not collect/retain it. Contract it with slab_bandmul! / slab_bandtmul! (native role) and slab_mirror! / slab_mirrort! (Hermitian mirror role).

PMSupermatrices{T}

Open handle to a memory-mapped ± supermatrix AO-integral store (see the module docstring for the on-disk format). In the usual flow you never touch this directly — ao_integrals(EC) builds the store and the consumers (FockFactory.ao_JK!, cc.pm_K2!, the dressing sweeps) dispatch on it. To use it by hand: open_pm_store it, contract, close_pm_store!.

Notation

• nao — number of AOs (basis size); AO indices μ,ν,ρ,σ run 1:nao.
• npp — number of packed AO pairs, nao·(nao+1)÷2. This is the pair-space dimension: each supermatrix Vs,Va is conceptually npp × npp, its rows/cols indexed by packed pairs.
• tri(μν) — the packed index of AO pair (μ,ν), μ≤ν (QMTensors.uppertriangular_index); runs 1:npp. pair_luts(nao) gives the inverse (packed index → μ, ν).
• panel / block J — the store keeps only the lower block-triangle of Vs/Va, cut into column blocks. Panel J is a dense matrix: the columns c_J (a contiguous ket-pair range) for every bra-pair row ≥ first(c_J) (the rows above are the unstored Hermitian mirror).

Fields

• nao, npp — as above.
• σblocks[J] — the ket-2 orbital (σ) range that defines block J.
• pairblocks[J] — c_J, the contiguous packed pair range panel J owns as its columns.
• offsets[J] — 1-based start of panel J inside each flat memory map.
• sio,smap / aio,amap — open IOStream + mmapped flat vector for Vs / Va; the panels are reshaped views into these (via spanel / apanel).

Panel J has shape (npp − first(c_J) + 1) × |c_J| (column-major): its top |c_J| rows are the diagonal tile V[c_J, c_J] (diagtile), the rest the sub-diagonal V[last(c_J)+1:npp, c_J] (subpanel).

Example — kext contraction K2[tri(pq)] = Σ_rs ⟨pq|rs⟩ D2[rs]

Its pair-space form is K2 = Vs·Ds + Va·Da, where Ds,Da are the ± fold of the density over ket pairs (both npp × m, m right-hand sides — e.g. nocc² for the CC 4-external term):

pm = open_pm_store(EC)                        # T inferred from the stored files
Ks = zeros(eltype(pm), pm.npp, m)             # eltype(pm) === T
Ka = zeros(eltype(pm), pm.npp, m)
pm_matmul!(Ks, pm, :s, Ds)                    # Ks = Vs · Ds   (zero-copy panel GEMMs)
pm_matmul!(Ka, pm, :a, Da)                    # Ka = Va · Da
K2 = Ks .+ Ka                                 # unpack tri(pq) → K2[p,q] with pair_luts(pm.nao)
close_pm_store!(EC, pm)

Example — low-level panel access (zero-copy views into the mmap)

for J in 1:pm_nblocks(pm)
  S = spanel(pm, J)                           # Vs panel: (npp-first(c_J)+1) × |c_J|
  Sdiag = diagtile(S, pm, J)                  #   top |c_J| rows = Vs[c_J, c_J]
  Ssub  = subpanel(S, pm, J)                  #   remaining rows = Vs[below, c_J]
end

PMWriter{T}

Write handle for building a ± store block by block: the panel layout (from the σ-block breakpoints) plus the two open output mmaps. Fill with pm_write_block! (any block order), finalize with pm_close_writer!.

SlabWork{TF}(pm::PMSupermatrices) -> SlabWork{TF}

Reusable scratch for a pm_slab_sweep!: the reconstructed nao×nao slab G, two length-nao buffers for the mirror-role GEMVs, and the packed-pair lookup tables. TF is the working element type — take TF = promote_type(eltype(pm), eltype(data)) so a real store still contracts against complex data (a real store gives TF = eltype(pm) = T).

add_mirror_row!(out, i, w, D, j, lo, transposed)

Mirror-role GEMV: out[i,:] += Σ conj(G)·D[j,·] over the L-band (lo, nao] — the current slab read as the bra ⟨ρσ|··⟩. Runs as conj(band · conj(D[j,:])) so the BLAS happens between the conjs (both no-ops for real TF); transposed selects band_tmul! over band_mul!. Uses the w.tv/w.tw scratch of the SlabWork.

Antisymmetric-panel J: dense V_a[first(c_J):npp, c_J] view (no copy).

band_htrans!(H, L, G, lo, hi)

Half-transform over the L-band: H[i,y] = Σ_x G[x,y] L[x,i] restricted to max(x,y) ∈ (lo,hi] (two rectangle GEMMs, no zero-padding flops; columns > hi of H are zeroed). Plain contraction — no conjugation of L (the coefficient convention of the dressing half-transforms). The reconstructed slab G over the panel's L-band (via the Fock builders' reconstruct_slab!) is the input.

band_mul!(y, G, lo, hi, x)   /   band_tmul!(y, G, lo, hi, x)

y += L·x (resp. transpose(L)·x) where L is the L-band max(row,col) ∈ (lo,hi] of the reconstructed slab G: two rectangle GEMVs (bottom band + right band) covering the band exactly — no zero-padding flops, and the unstored corner (row,col ≤ lo) is never read. transpose (not adjoint): swapping the ket pair is exchange, conjugation-free.

Transposed band product y[1:hi] += G[lo+1:hi, 1:hi]' * x[lo+1:hi] (the mirror rows' contribution).

Close a ± store handle from open_pm_store (a no-op for the cached handle, which release_pm_store! owns).

Delete a half-transformed store (data + metadata).

Diagonal tile V[c_J, c_J] of a panel P (its top |c_J| rows).

eachslab(pm; TF=eltype(pm)) -> iterator of PMSlab

The iterator form of pm_slab_sweep!: for slab in eachslab(pm) … reconstructs each stored ket-pair slab ⟨μν|ρσ⟩ into shared scratch and yields a PMSlab, in the exact same visit order (native columns of a σ-block in turn). Lets a custom AO-store contraction be written as an ordinary loop with the band helpers (slab_bandmul! etc.) instead of a pm_slab_sweep! callback. The yielded slab aliases one reusable buffer — single pass, do not retain.

roles=:both instead yields ALL ket pairs (native columns + the Hermitian-mirror sub-panel rows) as BothSlabs carrying G and Gt, for pm_bra_half! — the dressing sweep.

Example — AO Coulomb build J[μ,ρ] = Σ_νσ ⟨μν|ρσ⟩ D[ν,σ] (= FockFactory.add_coulomb!(J, s, D))

J = zeros(eltype(pm), pm.nao, pm.nao)
for s in eachslab(pm)
  slab_bandmul!(view(J,:,s.ρ), s, view(D,:,s.σ))             # native: ⟨··|ρσ⟩ · D[:,σ]
  s.ρ < s.σ && slab_bandtmul!(view(J,:,s.σ), s, view(D,:,s.ρ))
  slab_mirror!(J, s.ρ, s, D, s.σ)                           # mirror: ⟨ρσ|··⟩ · D[σ,:]
  s.ρ < s.σ && slab_mirrort!(J, s.σ, s, D, s.ρ)
end

ht_build_dress!(EC, pm, C; key="ht_oAAA", ookey="ht_ooAA")

Build the persistent half-transformed store for the AO dressing (called ONCE per orbital set in ao_cc_setup!): key = Σ_μ⟨μν|ρσ⟩C[μ,i] in both ket orders (pm_half_trans) plus the T1-independent ookey = v_ooAA[i,j,ρ,σ] = Σ_μν⟨μν|ρσ⟩C[μ,i]C[ν,j] (ht_build_oo!). ao_dressed_ints then reads these each iteration instead of re-streaming the ± store — one ± sweep per CC run, not per iteration. C is the (T1-independent) active occupied bra coefficients. The key/ookey arguments let the unrestricted build reuse this per spin.

ht_build_dress_unrestricted!(EC, pm, La_o, Lb_o)

Unrestricted analogue of ht_build_dress!: build the per-spin half-transformed stores plus the three T1-independent doubly-occ blocks the open-shell dressing reads each iteration. The bra-occ transforms La_o/Lb_o are T1-independent, so this runs ONCE per orbital set in ao_cc_setup!:

• "ht_oAAA_a" / "ht_oAAA_b" : the α / β half-transformed stores (pm_half_trans);
• "ht_ooAA_a" = v_ooAA(αα) , "ht_ooAA_b" = v_ooAA(ββ) : the same-spin doubly-occ blocks;
• "ht_oOAA" = v_oOAA[i,J,ρ,σ] = Σ_μν⟨μν|ρσ⟩La_o[μ,i]Lb_o[ν,J] : the opposite-spin cross block (built off the α store's A-role columns × Lb_o, so no separate ± sweep).

ht_build_oo!(EC, ht, C2, ookey)

Build the doubly-bra-occ intermediate v[i,j,ρ,σ] = Σ_μν⟨μν|ρσ⟩ C1[μ,i] C2[ν,j] (T1-independent) from an open half-transformed store ht (whose bra-1 transform used C1) and a second occupied bra C2, saved (mmapped) under ookey with shape (ht.m, size(C2,2), nao, nao). Reads each A-role σ-column (ht_column_A! — the B slab is not needed, so this touches each stored element once) and contracts the free ν slot with C2. For the closed shell C2 == C1 (→ v_ooAA); for the opposite-spin cross block C1 = La_o, C2 = Lb_o (→ v_oOAA).

ht_column!(Aσ, Bσ, ht, σ) -> (Aσ, Bσ)

Fill the dense σ-column slabs (both [m·nao, nao]) for all ρ:

• Aσ[(i,ν), ρ] = Σ_μ ⟨μν|ρσ⟩ C[μ,i] (A-role, μ→occ) — feeds v_oAoA/v_ooAA;
• Bσ[(i,μ), ρ] = Σ_ν ⟨μν|ρσ⟩ C[ν,i] (B-role, ν→occ, via exchange symmetry) — feeds v_AooA.

The two roles read the SAME pairs with swapped ket orders (Aσ: o1 for the contiguous σ-column ρ≤σ, o2 for the row picks ρ>σ; Bσ: the opposite). The pair-packing is invisible to the caller.

ht_column_A!(Aσ, ht, σ) -> Aσ

A-role half of ht_column!: the slab Aσ[(i,ν), ρ] = Σ_μ⟨μν|ρσ⟩C[μ,i] alone.

Each stored element is the A entry of exactly ONE column, so a full σ-sweep with this reader touches the store exactly once — half the bytes of the two-slab ht_column!. Consumers that need the other ket order do NOT need the B slab: by the particle-exchange symmetry ⟨μν|ρσ⟩ = ⟨νμ|σρ⟩ the B-role block is the A-role block with the two ket coefficients swapped and the output permuted (2,1,4,3) (see ht_mo_block). Only the per-column dressing sweeps, which contract both roles into the SAME output, genuinely need both slabs.

ht_column_B!(Bσ, ht, σ) -> Bσ

B-role half of ht_column!: the slab Bσ[(i,μ), ρ] = Σ_ν⟨μν|ρσ⟩C[ν,i] alone. Like ht_column_A! it touches each stored element exactly once, so a full σ-sweep reads the store once rather than twice.

A block builder does NOT need this: the B ordering is reachable from the A slab by particle exchange (see ht_mo_block), and a builder is free to choose its output permutation. It is for consumers whose output layout is fixed by someone else and for which the B order is the one that writes contiguously — save_ao_AAAo! streams ⟨μν|ρi⟩ into [μ,ν,ρ,i] and gets a contiguous [:,ν,:,:] slice per column from the B slab, where the A slab would leave it writing single strided elements.

True if the half-transformed store key (data + key*"_meta") exists on file.

ht_jk_columns!(t, ht, DJ, DK) -> t

Contract every σ-column's B- and A-role slabs with the AO densities DJ / DK over both of their free AO slots, for all σ in one go (t is [m, nao] and is overwritten):

t[i,σ] = Σ_xy ( B_σ[i,x,y]·DJ[x,y] − A_σ[i,x,y]·DK[x,y] )
       = Σ_xy ( ⟨x i|y σ⟩·DJ[x,y]  − ⟨i x|y σ⟩·DK[x,y] )

with the slabs exactly as ht_column! defines them (σ = ket-2 in both). That is the [bra-occ, AO-ket] block of a generalized J(DJ) − K(DK) Fock — the contraction the Λ residual's dD1 generalized-Fock term needs (CoupledCluster.dD1_ht_vo), which wants only that block and so never has to build an nao×nao Fock.

One sequential pass, each stored element read EXACTLY once — half of what a full ht_column! sweep reads (that gathers every pair twice: once as a column entry, once as a row pick) and with no dense per-column slab materialised at all. The four contributions of a stored pair p = tri(a,b), a ≤ b, are read straight off ht_column!'s two cases:

map[:,1,p] = A_{σ=b}[:,y=a] = B_{σ=a}[:,y=b]     (ket order (a,b))
map[:,2,p] = A_{σ=a}[:,y=b] = B_{σ=b}[:,y=a]     (ket order (b,a))

so each map[:,o,p] — a contiguous m×nao block, hence a BLAS gemv — is used twice, against one column of DJ/DK. On the pair diagonal a == b the two ket orders coincide with ht_column!'s single ρ ≤ σ case, which is why only the σ=b pair of contributions is applied there.

Threaded over contiguous b-ranges holding equal numbers of pairs (so every task still streams the map sequentially), each with a private accumulator: race-free and deterministic. DJ/DK need not be symmetric, and no conjugation is applied anywhere, so this is element-type generic.

ht_occ_early(EC, Ro; key="ht_oAAA", ookey="ht_ooAA") -> (v_ooAA, v_AooA, v_oAoA)

The occ-early intermediates read from a prebuilt half-transformed store (ht_build_dress!): v_ooAA is loaded (T1-independent); the T1-dependent ket contractions v_oAoA[i,ν,j,σ] = Σ_ρ Aσ[(i,ν),ρ]·Ro[ρ,j] and v_AooA[μ,i,j,σ] = Σ_ρ Bσ[(i,μ),ρ]·Ro[ρ,j] are one BLAS GEMM per σ off the A-/B-role slabs from ht_column!. Bit-identical to pm_occ_early but without re-streaming the ± store. key/ookey select the store (per-spin in the unrestricted path).

ht_occ_early_unrestricted(EC, Ra_o, Rb_o) -> (ssa_in, ssb_in, os_in)

All the open-shell occ-early intermediates the unrestricted dressing needs, read from the prebuilt per-spin half-transformed stores (ht_build_dress_unrestricted!) in ONE pass per store — the fused replacement for two ht_occ_early calls plus a pm_os_sweep. Reading each store once (its A-/B-role slabs both come from one ht_column! call) and reusing the shared products removes the redundant re-reads/re-GEMMs of the separate passes (6 GEMMs + 2 reads per σ-column vs 8 + 4). The doubly-occ blocks (v_ooAA(αα)/v_ooAA(ββ)/v_oOAA) are T1-independent and loaded.

Per α column c (Aa/Ba = A-/B-role), particle symmetry ⟨μν|ρσ⟩=⟨νμ|σρ⟩ routing every kept-AO to the fixed second-ket column: Aa·Ra_o→v_oAoA(αα) (shared by ssa AND os); Ba·Ra_o→v_AooA(αα); Ba·Rb_o→v_oAAO. Per β column: Ab·Rb_o→v_oAoA(ββ) (ssb) AND v_AOAO (os); Bb·Rb_o→ v_AooA(ββ); Bb·Ra_o→v_AOoA. Returns the tuples for ao_ss_finish(α), ao_ss_finish(β), ao_os_finish. Ra_o/Rb_o are the T1-dependent occupied ket coefficients.

open_ht_store(EC, key) -> HTStore    — memory-map an existing half-transformed store (read-only).

open_pm_store(EC) -> PMSupermatrices

Memory-map an existing ± store (read-only) and rebuild the layout from ao_pm_meta. The handle is cached and reused across calls for the same scratch directory — see the note above and release_pm_store!; close_pm_store! therefore does NOT unmap a cached handle.

Row-pair decode: lutμ[tri(μ,ν)] = μ, lutν[tri(μ,ν)] = ν for μ ≤ ν.

pm_bra_half!(hν, hμ, s::BothSlab, C) -> (hν, hμ)

Half-transform the bra pair of slab s with coefficients C ([nao, m]), filling hν[i,ν] = Σ_μ ⟨μν|ρσ⟩ C[μ,i] (μ→i, ν kept AO) and hμ[i,μ] = Σ_ν ⟨μν|ρσ⟩ C[ν,i] (ν→i, μ kept AO) — the two band GEMMs (band_htrans! on G and Gt) that the fused dressing outputs SHARE (compute once, reuse for every output). hν,hμ are caller-owned [m, nao] buffers (allocate once outside the loop). Cheap second-stage contractions (hν/hμ with the remaining coeffs) then form the outputs.

Flush both files and store the layout metadata (ao_pm_meta).

pm_from_joint!(EC; maxcols=pm_default_maxcols(nao))

One-time build of the ± store from the joint triangular ao_int2. Streams the packed ket columns block by block (each c_J is a contiguous packed range ⇒ sequential mmap read), folds each nao×nao slab into its ±-symmetrized bra pairs with calc_tri_sym_antisym!, and writes the panels via a PMWriter. Sequential read, sequential write, one pass. (For generation without a joint intermediate see IntegralTools.pm_integrals!.)

pm_half_trans(EC, pm, C, key) -> key

Build the half-transformed store key (data + key*"_meta") from the ± store pm and bra coefficients C[nao, m]: Σ_μ ⟨μν|ρσ⟩ C[μ,i] for every ket pair, in both ket orders. Written by an out-of-core gather — one output pair-block J at a time: the native contribution from panel J's own columns (sequential mmap), plus the Hermitian-mirror contribution of block-J pairs as sub-panel rows of every earlier panel, read via pread+WILLNEED (contiguous per-column runs). Each ± element is read ≈ twice; only band_htrans!-transformed blocks are held in RAM, each output block written once. Reuses the exact eachslab(:both) reconstruction, reordered by output block (so every reconstruct_*_both!+band_htrans! is self-contained).

pm_matmul!(out, pm, which, X) -> out

Left-multiply by a stored ± supermatrix: out .= V · X, V = Vs (which===:s) or Va (which===:a), reconstructed on the fly from its lower block-triangle. X, out are npp × m. Per panel: the diagonal tile once, the sub-diagonal in its own role ('N') and — by hermiticity — its mirror role via adjoint (BLAS 'C'/'T', no copy; a no-op transpose for real T). Every stored element fuels exactly its own and (sub-diagonal) its mirror product, so out = V·X for the full Hermitian V. Zero-copy panel GEMMs.

pm_matvec!(out, pm, which, X) -> out

Pair-space matvec out = Vs·X (which=:s) or out = Va·X (which=:a) over packed AO pairs — the verb name for the primitive pm_matmul!.

pm_occ_early(pm::PMSupermatrices, Lo, Ro) -> (v_ooAA, v_AooA, v_oAoA)

Occ-early half-transform on the persisted ± supermatrix store — three intermediates at half the integral streaming (each stored element read once, ≈ n⁴/4; flop parity). One eachslab(pm; roles=:both) pass does the shared bra half-transform pair (pm_bra_half! → hν[i,ν] = Σ_μ ⟨μν|ρσ⟩ Lo[μ,i] (μ→occ), hμ[i,μ] = Σ_ν ⟨μν|ρσ⟩ Lo[ν,i] (ν→occ)) and, fused in-cache, the v_ooAA GEMM; it also ACCUMULATES the half-transforms into A[(i,ν),ρ,σ] / B[(μ,i),ρ,σ] (the bra index fused LEADING so each σ-slice is a BLAS-contiguous matrix). The ket-contracted v_oAoA/v_AooA are then ONE BLAS GEMM per σ off A/B — batching the ket ρ-sum instead of a per-slab rank-1 update (memory-bound at BLAS-2). This is ≈1.1–1.6× the old per-slab outer-product sweep and grows with nao. A/B cost 2·nocc·nao³ in RAM (fine up to a few hundred AOs; an mmap fallback for larger systems is a follow-up).

pm_os_sweep(pm::PMSupermatrices, La_o, Ra_o, Lb_o, Rb_o) -> (v_oOAA, v_AOoA, v_oAoA, v_oAAO, v_AOAO)

The opposite-spin occ-early sweep on the ± store (the five intermediates of ao_os_blocks at half the streaming), as a eachslab(pm; roles=:both) sweep like pm_occ_early but with TWO coefficient sets: pm_bra_half! gives the shared half-transforms of La_o (hνa/hμa) and Lb_o (hνb/hμb).

pm_slab_sweep!(f, pm, w)

Visit each stored ket-pair column of the ± store pm exactly once: reconstruct its dense slab ⟨μν|ρσ⟩ into w.G (a SlabWork) and call f(ρ, σ, native_lo, mirror_lo). (ρ,σ) is the ket pair; native_lo/mirror_lo are the L-band lower bounds for the native/mirror roles (see the section comment above). f contracts w.G with band_mul!/band_tmul! (native role) and/or add_mirror_row! (mirror role).

Example — an AO Coulomb build J[μ,ρ] = Σ_νσ ⟨μν|ρσ⟩ D[ν,σ]

w = SlabWork{eltype(pm)}(pm)                 # real store: TF = eltype(pm) = T
J = zeros(eltype(pm), pm.nao, pm.nao)
pm_slab_sweep!(pm, w) do ρ, σ, native_lo, mirror_lo
  band_mul!(view(J,:,ρ), w.G, native_lo, pm.nao, view(D,:,σ))   # native: ⟨··|ρσ⟩ · D[:,σ]
  ρ < σ && band_tmul!(view(J,:,σ), w.G, native_lo, pm.nao, view(D,:,ρ))
  add_mirror_row!(J, ρ, w, D, σ, mirror_lo, false)              # mirror: ⟨ρσ|··⟩ · D[σ,:]
  ρ < σ && add_mirror_row!(J, σ, w, D, ρ, mirror_lo, true)
end

(FockFactory.ao_JK! wraps exactly this — Coulomb + exchange — into add_coulomb!/add_exchange!.)

pm_to_joint!(EC)

Inverse of pm_from_joint!: reconstruct the joint triangular "ao_int2" file (int2[μ,ν,tri(ρσ)] = ⟨μν|ρσ⟩) from the ± supermatrix store. One panel-major pass; each stored element writes its two native slab positions (± inversion) and — sub-panel rows — its two Hermitian-mirror positions (conj; the mirrors of diagonal-tile elements are their own stored partners). Used to serve joint-format consumers (e.g. the AO→MO transform) when only the ± store is on disk.

pm_write_block!(w::PMWriter, J, S, A)

Write the ±-folded block J: S/A hold the full-height npp × |c_J| symmetric/ antisymmetric combinations for the ket columns c_J; only the lower rows ≥ first(c_J) are stored (the upper rows are the Hermitian mirror owned by earlier panels).

Open a ± store for writing with the given σ-block breakpoints (any valid blocking).

reconstruct_slab!(w, Ps, Pa, jc, r0)

Fill w.G[μ,ν] = ⟨μν|ρσ⟩ for stored column jc (ket pair (ρ,σ)) of the ± panels Ps = spanel(pm,J), Pa = apanel(pm,J) (r0 = first(pairblocks[J])), via the inversion ⟨μν|ρσ⟩=(Vs+Va)/2, ⟨νμ|ρσ⟩=(Vs−Va)/2 (the stored ×2 row diagonals make μ=ν come out right with no special-casing). Only the panel's L-band bra pairs are written; the untouched corner is never read by the band GEMVs. Normally called for you by pm_slab_sweep!.

release_pm_store!()

Close and forget the cached ± store handle (see open_pm_store). Must be called before the store files are rewritten or deleted; a no-op when nothing is cached.

slab_bandmul!(y, s::PMSlab, x) — native role: `y += band(⟨··|ρσ⟩)·x` (ket order ρσ).

slab_bandtmul!(y, s::PMSlab, x) — native role, transposed: `y += band(⟨··|ρσ⟩)ᵀ·x` (ket order σρ).

slab_mirror!(out, i, s::PMSlab, D, j) — Hermitian mirror role: `out[i,:] += conj(band)·conj(D[j,:])`.

slab_mirrort!(out, i, s::PMSlab, D, j) — Hermitian mirror role, transposed band.

Symmetric-panel J: dense V_s[first(c_J):npp, c_J] view (no copy).

Sub-diagonal block V[last(c_J)+1:npp, c_J] of a panel P (its remaining rows).

GEMM-friendly σ-block column width (see pm_default_maxcols).

Gslab(s::PMSlab) -> Matrix — the reconstructed dense slab `⟨μν|ρσ⟩` (aliases scratch).

pm_breakpoints(nao; maxcols) -> Vector{Int}

σ-aligned block breakpoints: greedily grow a block until its column count Σσ would exceed maxcols (each σ contributes σ columns). maxcols ≥ nao guarantees every single σ fits, so no block is empty.

pm_default_maxcols(nao) -> Int

Default σ-block column ceiling for the ± store layout: clamp(PM_BLOCK_COLS, nao, npp).

Deterministic — a function of nao only, NOT of available memory. The layout is persisted in ao_pm_meta and rebuilt from it by every reader (open_pm_store), so a memory-dependent width would make the on-disk block structure non-reproducible (it would differ between the writing machine/run and any other) and size the build buffer to volatile free memory. A fixed width is reproducible and keeps the build buffer (2·npp·maxcols·sizeof(T)) predictable. The clamp bounds: ≥ nao is required (a single σ contributes up to nao columns), and ≤ npp means a tiny store (npp < PM_BLOCK_COLS, i.e. nao ≤ 31) is a single block (the full square) — the block-triangle saving needs ≥2 blocks anyway. Between the two, PM_BLOCK_COLS batches enough columns per panel for efficient BLAS-3.

pm_layout(nao, breakpoints) -> (σblocks, pairblocks, offsets, totlen)

Reconstruct the full panel layout from the σ-block breakpoints (the last σ of each block; breakpoints[end] == nao). Pure/closed-form: writer and reader call it identically.