Integral tools

This module contains tools for generating integral dumps from AO integrals — both density-fitted and exact (non-DF) — and transforming them AO→MO.

ao_integrals(EC::ECInfo) -> Float64

Generate the exact (non-density-fitted) AO integrals of the current system and store them as files: the 2-e integrals <μν|ρσ> (physicists' notation, triangular (ρ,σ) packing) on the memory-mapped file "ao_int2", and the 1-e integrals "S_AA"/"h_AA" (see save_ao_1e_integrals!). This is the entry point behind the @ints macro.

The 2-e integrals are assembled batch-wise straight into the memory-mapped triangular int2[μ,ν,tri(ρ,σ)] = ⟨μν|ρσ⟩ (see eri_2e4idx_tri!); the full nao⁴ tensor is never materialized.

If EC.fd currently holds (MO/FCIDUMP) integrals, they are discarded with a warning — a fresh AO integral generation supersedes any previously loaded/generated MO dump.

Return the nuclear repulsion energy.

calc_system_df_integrals(EC::ECInfo)

Calculate 3-index integrals for the EC.system and store them in mmL file. The routine is intended to be used in a combination with FDump integrals.

contract_df_integrals!(EC::ECInfo)

Contract 3-index DF integrals from scratch file mmL into 4-index integrals and store them in EC.fd.int2.

The 3-index integrals $B_p^{qL}$ are stored in the scratch file mmL with shape (norb, norb, naux). The 4-index integrals are computed as: $v_{pr}^{qs} = \sum_L B_p^{qL} B_r^{sL}$ and stored in upper-triangular format for the last two indices (q ≤ s).

After contraction, EC.fd.df3idx is set to false.

delete_ao_integrals!(EC::ECInfo)

Delete all exact AO integral scratch files (the ± supermatrix store, "S_AA", "h_AA") if present. Called when the geometry or basis changes — then even the 2-e AO integrals (μν|ρσ) and the overlap are invalid. For a pure nuclear charge/dummy change use invalidate_ao_1e_integrals!, which keeps the (unchanged) 2-e integrals.

ensure_ao_integrals!(EC::ECInfo; method="@hf", alternative="@bohf")

Make sure the exact AO integral files for the current system exist: the 2-e integrals ("ao_int2") are generated with ao_integrals when the file is missing, and the cheap 1-e integrals ("S_AA"/"h_AA") are always refreshed. The files are invalidated on geometry/basis changes by @setupEC.

If EC.fd holds (MO/FCIDUMP) integrals, they are discarded with a warning: method (@hf/@uhf) runs on exact AO integrals, and a leftover MO dump would shadow the AO flow in subsequent correlated calculations. To run HF directly on FCIDUMP integrals, use alternative (@bohf/@bouhf) instead.

generate_3idx_integrals(EC::ECInfo, cMO::SpinMatrix, fitbasis="mpfit"; save3idx=true)

Generate $v_p^{qL}$ with $v_{pr}^{qs} = v_p^{qL} δ_{LL'} v_r^{sL'}$ and store in file mmL. If save3idx is false, no 3-index integrals are calculated, only save pseudo-square-root-inverse Cholesky decomposition.

generate_AO_DF_integrals(EC::ECInfo, fitbasis="mpfit"; save3idx=true)

Generate AO integrals using DF + Cholesky. If save3idx is true, save Cholesky-decomposed 3-index integrals, otherwise save pseudo-square-root-inverse Cholesky decomposition.

Return nuclear repulsion energy.

generate_DF_Fock(EC::ECInfo, cMO::SpinMatrix; check_diagonal=false)

Generate DF Fock matrix in MO basis. If check_diagonal is true, check the off-diagonal elements of the Fock matrix to be small. The Fock matrix is saved in files f_mm/f_MM and orbital energies in e_m/e_M.

Return reference energy.

generate_DF_integrals(EC::ECInfo, cMO::SpinMatrix; save3idx=true)

Generate $v_p^{qL}$ and $f_p^q$ with $v_{pr}^{qs} = v_p^{qL} δ_{LL'} v_r^{sL'}$. The $v_p^{qL}$ are generated using mpfit fitting basis, and the $f_p^q$ are generated using jkfit fitting basis. The integrals are stored in files mmL and f_mm.

Return reference energy (calculated using jkfit fitting basis).

generate_mo_dump(EC::ECInfo, cMO::AbstractMatrix; core::AbstractMatrix, persistent=false) -> FDump

Build an MO-basis FDump in EC.fd from the exact AO integral files ("ao_int2"/"h_AA", see ao_integrals) and the (restricted/closed-shell) MO coefficients cMO[μ,p]. This is the non-DF analogue of dfdump.

Exact (non-density-fitted) O(N⁵) four-index transformation, written straight onto a fresh memory-mapped scratch file ("mo_int2"): the full MO tensor is never materialized in memory. cMO may be rectangular (nao × nout with nout ≤ nao): only those nout orbitals are kept (e.g. deleted virtuals and frozen virtuals excluded).

By default the file is temporary (the dump is transient and re-derived on demand); persistent=true keeps it for the rest of the session, so that the dump survives the delete_temporary_files! at the end of a driver run (used by @moints).

NELEC is stored as the full (neutral) electron count (without frozen core) A frozen core is folded in the AO basis via core, so only the active orbitals are transformed.

generate_mo_dump(EC::ECInfo, cMO::SpinMatrix; core::SpinMatrix, persistent=false) -> FDump

Build an MO-basis FDump in EC.fd from the exact AO integral files and the MO coefficients cMO. For a restricted cMO this builds a closed-shell (RHF) dump (see the matrix method); for an unrestricted cMO it builds an unrestricted (UHF) dump with the spin blocks int2aa/int2bb (v_{pq}^{rs} in each spin, triangular) and int2ab (v_{pQ}^{rS}, full 4-index), and per-spin 1-e integrals — the exact-AO analogue of the rhf→uhf branch of transform_fcidump!. Both spins must have the same orbital count (a single NORB); each block may be rectangular (deleted / frozen-virtual orbitals dropped). A frozen core is folded in the AO basis via core, so only the active orbitals are transformed. persistent is as in the matrix method.

invalidate_ao_1e_integrals!(EC::ECInfo)

Invalidate only the 1-electron AO integral file "h_AA" (the core Hamiltonian T + V), which depends on the nuclear charges, so it is recomputed for the current system on demand. The exact 2-e integrals "ao_int2" and the overlap "S_AA" depend only on the basis functions (positions + basis), so they are kept — e.g. across a charge/@dummy change, where ghost atoms retain their basis functions and the ERIs are unchanged.

pm_transform(EC, pm, Tl, Tl2, Tr, Tr2, key; triangular=true) -> int2t
pm_transform(EC, pm, C, key; kw...) -> int2t                      # RHF shorthand (C,C,C,C)

Transform the ± AO store to the MO basis, int2t[p,q,rs] = Σ_μνρσ ⟨μν|ρσ⟩ Tl[μ,p] Tl2[ν,q] Tr[ρ,r] Tr2[σ,s], without ever reconstructing the joint nao⁴/2 int2 — the high-level verb that hides the pair-space-vs-N⁵ choice. Dispatches to the all-BLAS-3 pair-space pm_transform_int2 for a small target space and the batched N⁵ pm_transform_int2_n5 for a near-full one (pm_transform_worthwhile). triangular packs the ket pair tri(r,s) (needs size(Tr,2)==size(Tr2,2)); false writes the full dense [p,q,r,s] (the UHF αβ block).

save_ao_1e_integrals!(EC::ECInfo) -> BasisSet

Compute and store the 1-e integrals of the current system in the AO basis: the overlap $S_{μν}$ and the core Hamiltonian $h_{μν} = T_{μν} + V_{μν}$ under the standard "S_AA"/"h_AA" scratch keys. Cheap (2-index); (re)computed on demand by every consumer of the AO integral files. Return the AO basis.

transform_3idx!(EC::ECInfo, fname::String, U::AbstractMatrix)

Transform 3-index integrals in-place: $B_{pq}^{L} \leftarrow U^\dagger B U$. The integrals are memory-mapped from file fname.

transform_fcidump!(EC::ECInfo, fd::FDump, Tl::SpinMatrix, Tr::SpinMatrix)

Transform the integrals of fd in place to a new basis using Tl, Tr. If Tl/Tr are unrestricted, an RHF dump is turned into a UHF dump. The transformed 2-e integrals are written to memory-mapped scratch files (their size bounded by available_memory(EC)); the 1-e integrals are transformed in memory. Intended as a one-shot rotation — the scratch keys are the block names, which must not already back fd's current integrals (they never do for the in-memory / mo_*-backed dumps this is called on).

transform_fcidump!(fd::FDump, Tl::SpinMatrix, Tr::SpinMatrix)

Like transform_fcidump!(EC, fd, Tl, Tr) but keeps the transformed 2-e integrals in memory (no EC / scratch files needed) — convenient for ad-hoc FDump manipulation. For large dumps prefer the EC method, which memory-maps the result.

transform_int2(int2, Tl, Tl2, Tr, Tr2; membudget=available_memory()) -> int2t

In-memory triangular-output transform (allocates a zeros output and calls transform_int2!). For large integrals prefer the memory-mapped transform_int2(EC, …, key).

transform_int2(EC::ECInfo, int2, Tl, Tl2, Tr, Tr2, key) -> int2t

Triangular-output transform writing the result to a fresh memory-mapped scratch file named key (newmmap); the memory budget comes from available_memory(EC) (honoring @set mem budget/fraction). key must not name the scratch file currently backing int2.

transform_int2_Q(int2, Tl, Tl2, Tr, Tr2; membudget=available_memory()) -> int2t

In-memory full 4-index (dense) transform (allocates a zeros output). For large integrals prefer the memory-mapped transform_int2_Q(EC, …, key).

transform_int2_Q(EC::ECInfo, int2, Tl, Tl2, Tr, Tr2, key) -> int2t

Full 4-index (dense) transform writing the result to a fresh memory-mapped scratch file named key; budget from available_memory(EC). Accepts a triangular (Array{T,3}) or dense (Array{T,4}) input.

pm_integrals!(EC::ECInfo, bao; maxcols=0)

Fused generation of the ± supermatrix store straight from the ERI generator: the triangular ket-column blocks are assembled in a bounded RAM slab (shell-aligned σ-blocks, batches within a block generated in parallel), ±-folded with calc_tri_sym_antisym! and written as PM panels — the joint ao_int2 intermediate is never created (disk ≈ n⁴/4 throughout, no transient peak). maxcols bounds the block width (0 = the deterministic PMStore.pm_default_maxcols default). Shell quartets whose Cauchy–Schwarz bound falls below int.screen are skipped (see schwarz_bounds).

The generation computes exactly the stored lower block-triangle (rowcut, default on): bra pairs whose packed row falls below a block's floor are the conj-Hermitian mirror that PMStore.pm_write_block! never stores, so their quartets are skipped in the ERI kernel (bra shell pairs entirely below first(σblock)) and the ± fold is restricted to the kept rows — roughly half the quartets and fold flops, bit-identical store. rowcut=false restores the full computation (for testing).

pm_pairfold(pm, Ca, Cb, pairs) -> (Ds, Da)

±-fold of the coefficient outer products Ca[:,p]⊗Cb[:,q] over packed AO pairs, for the list of MO pairs pairs: Ds[tri(μν),j] = ½(Ca[μ,p]Cb[ν,q] + Ca[ν,p]Cb[μ,q]), Da = ½(… − …) (½ on the μ=ν diagonal, Da diagonal 0) — the pair-space transform matrix (the pm_matmul! "density" convention). Each is npp × length(pairs).

pm_transform_int2(EC, pm, Tl, Tl2, Tr, Tr2, key; triangular=true) -> int2t

On-the-fly AO→MO 4-index transform directly from the ± supermatrix store pm, entirely in BLAS-3 and without ever reconstructing the joint AO integrals — the analog of transform_int2 / transform_int2_Q. Writes int2t[p,q,rs] = Σ_μνρσ ⟨μν|ρσ⟩ Tl[μ,p] Tl2[ν,q] Tr[ρ,r] Tr2[σ,s] to a fresh mmap key via the pair-space identity ⟨pq|rs⟩ = 2·(Lsᵀ·Vs·Rs + Laᵀ·Va·Ra) (see the section comment). triangular=true packs the ket pair tri(r,s) (r≤s, requires size(Tr,2)==size(Tr2,2)); false writes the full dense [p,q,r,s] (the UHF αβ block). Cost O(nao⁴·nmo²) — intended for nmo ≪ nao (see pm_transform_worthwhile). The ket-output pairs are blocked to keep the working set within membudget (the bra-fold Ls,La is held once).

pm_transform_int2_n5(EC, pm, Tl, Tl2, Tr, Tr2, key; triangular=true) -> int2t

N⁵ AO→MO 4-index transform directly from the ± store, for a NEAR-FULL target space (where the pair-space pm_transform_int2 would be N⁶) — and, crucially, WITHOUT ever materializing the joint nao⁴/2 int2 on disk.

Batched index-by-index quarter-transforms, all BLAS-3. Per p-block: a chunked two-role sweep reconstructs the ± slabs into a bounded buffer and bra-transforms them into H2[p,q,ρ,σ] = Σ_μν ⟨μν|ρσ⟩ Tl[μ,p] Tl2[ν,q], then the ket pair is transformed Σ_ρσ H2 Tr[ρ,r] Tr2[σ,s]. Peak memory ≈ H2 (|pb|·nq·nao²) + the O(nao²·chunk) slab buffer; both blocked to membudget.

pm_transform_worthwhile(nout, nao, T, membudget) -> Bool

Whether the pair-space pm_transform_int2 is the right choice over the batched N⁵ pm_transform_int2_n5 (both transform directly from the ± store — neither materializes the joint int2). Two conditions:

• speed — 10·nout ≤ nao: the pair-space is O(nao⁴·nmo²), so it wins comfortably only while the target space is a fraction of the basis (an active space); as nout → nao it becomes N⁶ and the N⁵ slab transform is faster, so route near-full transforms there.
• memory — the bra-fold Ls,La (2·npp·nout²) must fit comfortably in membudget.

transform_int1(int1::AbstractArray, Tl::AbstractArray, Tr::AbstractArray) -> int1t

Transform 1-e integrals to a new basis: int1t[p,q] = int1[p',q'] Tl[p',p] Tr[q',q].

transform_int2!(int2t, int2::Array{T,3}, Tl, Tl2, Tr, Tr2; membudget) -> int2t

In-place triangular-output transform: write the transformed 2-e integrals into the preallocated int2t of size (nout, nout, nout*(nout+1)÷2) (in-memory or memory-mapped). Each packed output column is written exactly once, so int2t need not be zero-initialized. Streams the triangular integrals in blocks of the output partner index (via transform_int2_pqs_block!) and finishes with a matrix multiplication of the remaining index restricted to r ≤ s (so the triangular symmetry is exploited — only r ≤ s is formed). Requires the same-spin pattern Tl≡Tl2, Tr≡Tr2 (implicit in triangular storage).

transform_int2_Q!(int2t, int2::Array{T,3}, Tl, Tl2, Tr, Tr2; membudget) -> int2t

In-place full 4-index (dense) output transform from a triangular input. int2t must have size (size(Tl,2), size(Tl2,2), size(Tr,2), size(Tr2,2)) (may be rectangular / spin-mixed, e.g. the int2ab block). nin == size(int2,1) must match the rows of all four matrices.

transform_int2_Q!(int2t, int2::Array{T,4}, Tl, Tl2, Tr, Tr2; membudget) -> int2t

In-place full 4-index (dense) output transform from a dense 4-index input (int2t of size (norb, norb, norb, norb)).

transform_int2_blocksize(ns, nin, np, nq, elsize; membudget)

Choose the width of the output-4th-index block for the streaming 4-index integral transform so the per-block working set (the co-live intermediates Z(nin³) + W(np·nin²) + Z2(np·nq·nin), each of elsize bytes) stays within membudget bytes. Capped at ns (a single block ⇒ the input is read once) and floored at 1. Blocking only tiles a loop, so the numerical result is independent of the choice — this purely trades peak memory against the number of passes over the input: small blocks (many passes) on a memory-starved node, a single pass when membudget is ample (e.g. a fat node).

transform_int2_pqs_block!(buf, int2::Array{T,3}, Tl, Tl2, Tr2, tb) -> Z2

Shared kernel for the streaming 4-index transform from the triangular input int2[p',q',tri(r',s')] = <p'q'|r's'> (r'≤s'). For the output-4th-index block tb (via Tr2), narrow the packed partner index s' into the block (correctly desymmetrizing the joint packing <p'q'|r's'> = <q'p'|s'r'>) and transform the two bra indices p'→p, q'→q, leaving the third input index r' untransformed. Returns Z2[p,q,r',u] (allocated from the buffer arena buf); the caller applies Tr to the remaining r' index. All heavy steps are BLAS-3 gemms.