SCF Fock-build modes¶

vibe-qc offers two Fock-build modes for molecular SCF, the strategy that constructs the J (Coulomb) and K (exchange) matrices each iteration. Density fitting is a separate route: when density_fit=True, its RI-J or RIJCOSX kernels supersede the four-index path and the mode is ignored. Within that route, cosx=True selects COSX-K. The high-level run_job(cosx=True) convenience also enables density fitting automatically.

For a worked resource comparison, including the 170-basis-function H2O2/cc-pVQZ timing and memory measurements, see Direct SCF or in-core integrals?.

Three modes¶

Mode	Enum value	Memory	Best for
`AUTO`	`SCFMode.AUTO`	-	Default, picks based on system size
`CONVENTIONAL`	`SCFMode.CONVENTIONAL`	O(n_bf⁴)	Small systems whose dense tensor fits comfortably
`DIRECT`	`SCFMode.DIRECT`	O(n_shells² + n_bf²)	Large systems, bounded integral memory

All molecular SCF drivers (RHF, UHF, RKS, UKS) support the mode via identical scf_mode / scf_mode_auto_threshold fields on their options struct. Range-separated RKS and UKS are a deliberate exception to the requested mode: their attenuated-exchange kernel requires the direct builder, as described below.

AUTO, the default¶

from vibeqc import RHFOptions, SCFMode

opts = RHFOptions()
# opts.scf_mode defaults to SCFMode.AUTO
# opts.scf_mode_auto_threshold defaults to 140 BF

When scf_mode = AUTO (the default), vibe-qc counts the number of basis functions after the basis set is constructed and decides:

(n_{bf} \le 140) -> CONVENTIONAL, the in-core 4-index ERI tensor is built once and the per-iteration cost is a tensor contraction.
(n_{bf} > 140) -> DIRECT, the in-core tensor would be multi-GB; switch to on-the-fly Schwarz-screened libint quartet evaluation that stays O(n_shells² + n_bf²) in memory.

The threshold is a resource policy, not a physical boundary. At the default cutoff, a dense 140-function tensor is already 2.862 GiB before any other process memory. To cap that tensor below 1 GiB, use 107, the largest basis count satisfying that limit:

from vibeqc import RHFOptions

opts = RHFOptions()
opts.scf_mode_auto_threshold = 107  # switch to DIRECT above about 1 GiB

CONVENTIONAL, the in-core path¶

from vibeqc import RHFOptions, SCFMode

opts = RHFOptions()
opts.scf_mode = SCFMode.CONVENTIONAL

The full 4-index electron-repulsion integral tensor (μν|λσ) is computed once at the start and reused across SCF iterations. Fast for small systems, but the tensor grows as n_bf⁴:

n_bf	Full tensor size
24	2.53 MiB
95	621.42 MiB
107	1000.06 MiB
108	1.014 GiB
170	6.223 GiB
200	11.921 GiB
250	29.104 GiB
394	179.546 GiB

These are \(n_\mathrm{bf}^4\times8\) bytes for the dense tensor alone. Add Fock, density, DIIS, eigensolver, runtime, and requested-output workspaces when sizing a job. Lower the AUTO cutoff when the tensor exceeds the workflow’s budget, or set DIRECT explicitly.

DIRECT, the integral-driven path¶

from vibeqc import RHFOptions, SCFMode

opts = RHFOptions()
opts.scf_mode = SCFMode.DIRECT

Instead of storing the ERI tensor, the direct Fock build considers the unique shell quartets on each SCF iteration and evaluates only those that survive screening. An 8-fold-symmetric loop skips the vast majority of quartets using the strict Cauchy-Schwarz bound:

[ |\langle \mu\nu | \lambda\sigma \rangle| ;\le; Q(s_\mu, s_\nu) \cdot Q(s_\lambda, s_\sigma) ]

where (Q(s_a, s_b)) is precomputed from the diagonal shell-pair integrals at construction time and reused across iterations.

Direct SCF options¶

All fields live on the per-method options struct (RHFOptions, UHFOptions, RKSOptions, UKSOptions):

Field	Default	Description
`schwarz_threshold`	1e-10	Per-quartet skip bound; lower = fewer skipped quartets
`schwarz_threshold_loose`	1e-7	Loose threshold for early SCF iterations (two-phase Schwarz screening)
`schwarz_threshold_tighten_at`	1e-3	Gradient-norm cutoff, switch from loose to tight below this
`incremental_fock`	`True`	Use a difference-density ΔP Fock build on the direct path
`incremental_fock_reset_freq`	8	Full-density rebuild every N iterations to limit recurrence error

Two-phase Schwarz tightening¶

By default, the direct Fock build starts with the looser schwarz_threshold_loose (1e-7) so the per-shell density envelope catches mid-SCF quartets aggressively. Once the SCF gradient norm drops below schwarz_threshold_tighten_at, the builder switches to the tight schwarz_threshold and discards the incremental cache. An energy difference across that threshold and cache change compares two different numerical Fock maps, so it is excluded from convergence testing.

Häser and Ahlrichs analyse iteration-dependent screening and show that the error in a screened difference-Fock update can inherit errors from earlier increments (Eqs. 30-35). That analysis motivates the coarse-to-tight refinement; it does not prescribe vibe-qc’s thresholds or transition policy. Set schwarz_threshold_loose <= schwarz_threshold to disable the coarsening and use a uniform tight threshold throughout.

Incremental ΔP Fock build¶

When incremental_fock = True, the builder caches the previous density D_prev and two-electron Fock matrix G_2e_prev. Each iteration computes ΔD = D − D_prev and returns:

[ G_{2e}[D] = G_{2e}[D_\text{prev}] + G_{2e}[\Delta D] ]

This recurrence is Eq. 20 of Almlöf, Fægri, and Korsell. Their discussion on p. 389 also recommends a final energy-oriented iteration after a non-energy monitor is satisfied. The exact certification policy below is vibe-qc’s implementation of that principle, not a prescription quoted from the paper.

The direct kernel applies the published density-weighted Schwarz test to the difference density. For an incremental update with (s = \max |\Delta D| < 1), vibe-qc first screens (\Delta D / s), builds that normalized contribution, and multiplies the result by (s). This is a vibe-qc relative-error policy: the omitted part of an update scales to zero with its amplitude, instead of each small update being allowed the same absolute omission and accumulating a persistent drift floor. It is derived from the linearity of (G_{2e}); it is not a scaling formula prescribed by Almlöf et al. or Häser and Ahlrichs.

Amplitude normalization deliberately prevents a uniform rescaling of (\Delta D) from making the retained quartet set artificially sparser. Incremental builds still benefit when the change density becomes spatially localized, while accuracy no longer depends on how one density change is split across SCF iterations.

A full rebuild happens every incremental_fock_reset_freq iterations (default 8) to limit error inherited from screened increments as well as ordinary arithmetic accumulation. The reset interval is a vibe-qc implementation policy, not a reset rule attributed to Almlöf et al. When the two-phase screen changes to its tight threshold, the driver also discards the incremental cache. Set incremental_fock=False for a controlled full-build comparison.

from vibeqc import RHFOptions

opts = RHFOptions()
opts.incremental_fock = True
opts.incremental_fock_reset_freq = 8

Final full-Fock certification¶

Loose-screened and incremental builds are acceleration phases. A molecular direct SCF calculation enters a certified full-Fock phase when either:

the orbital gradient first satisfies conv_tol_grad, making the row eligible for ordinary convergence; or
the restart monitor detects a high-gradient stall while a coarse Fock map is still active.

The driver then forces the configured tight Schwarz threshold, disables incremental updates, resets the Fock builder, and clears extrapolation and mixing state. The first tight, nonincremental full-density row has an energy difference relative to the preceding coarse-map row. That cross-map value is reported but cannot certify convergence. The following row provides an energy difference between two comparable tight, nonincremental full-density builds. Only then can the calculation converge, and both the original conv_tol_energy and conv_tol_grad must pass unchanged.

Consequently, certification requires two comparable same-map rows. It does not relax either tolerance or infer an empirical energy-noise floor. A standard diagnostic identifies a gradient-triggered or stall-triggered move into this phase.

When density_fit or cosx is active¶

At the options-struct level, density_fit=True selects the DF J/K infrastructure and causes it to ignore scf_mode. Inside that route, cosx=True replaces fitted exchange with COSX-K. A low-level options object with cosx=True but density_fit=False does not activate COSX; set both and provide aux_basis. The high-level run_job(cosx=True) argument performs those two steps and auto-selects an auxiliary basis when possible.

Range-separated hybrids force DIRECT¶

RKS and UKS range-separated hybrids need an erf-attenuated exchange build. Only the direct builder currently implements that kernel, so these jobs force DIRECT when density_fit=False, even if AUTO would choose conventional or the options request CONVENTIONAL. Combining a range-separated hybrid with density_fit=True currently raises an explicit unsupported-route error because the required attenuated three-center integrals are not implemented.

Performance¶

On compute-small (i9-10885H, 4 threads), RHF direct / def2-SVP, same-box ORCA comparison:

This sweep uses def2-SVP throughout, so its H2O row has 24 basis functions. The H2O/6-31G* examples use 18 basis functions instead.

System	n_bf	vibe-qc (s)	ORCA (s)	Ratio
H₂O	24	2.0	15.0	0.13
glycine	95	5.0	13.0	0.39
naphthalene	180	31.1	45.1	0.69
n-decane	250	226.3	67.2	3.37
n-hexadecane	394	639.2	185.3	3.45

This historical sweep placed the speed crossover against ORCA near 200 BF on that host. The current AUTO policy switches at 140 BF because the CPU-time crossover is not a safe memory cutoff: the dense tensor is already 2.862 GiB at 140 BF. Energy agreement against ORCA was ≤ 2.4 µHa across all converged cells.

That wall-time crossover does not define a safe memory crossover. A matched one-core H2O2/cc-pVQZ case at 170 basis functions used 45.550 s wall time, 45.032 s active CPU time, and 6611.05 MiB peak RSS with vibe-qc’s conventional path. ORCA 6.1.1 used 52.669 s wall time, 51.232 s active CPU time, and 126.25 MiB peak RSS in the same allocation. The energy difference was \(1.53\times10^{-10}\) Ha. The worked tutorial derives the \(170^4\) allocation, separates the different-host direct validation, and defines a safe crossover sweep.

Component timing¶

Set component_timing = True on RHFOptions or RKSOptions to record where a molecular restricted SCF spends its wall time. The option is off by default. Enabling it does not change energies, densities or iteration counts.

opts = vq.RHFOptions()
opts.density_fit = True
opts.aux_basis = "def2-universal-jkfit"
opts.component_timing = True
result = vq.run_rhf(mol, basis, opts)
ct = result.component_timing
for phase in ct.phases:
    print(phase.path, phase.inclusive_s, phase.exclusive_s, phase.calls)
print(ct.wall_s, ct.unclassified_s)

Phase paths join nested names with /. The inclusive time of a phase contains its nested phases, while its exclusive time does not. The exclusive times of all phases plus unclassified_s equal wall_s, which runs from driver entry to return.

Phase	Covers
`setup.one_electron`, `setup.grids`, `setup.jk_construction`, `setup.guess`, `setup.scf_init`	Integrals, XC grid, J/K builder, DF tensor and COSX caches, initial guess and orthogonalization
`fock.jk_fused`	One RHF `build_g_rhf` call. Direct and conventional builders fuse J and K, so no separate J/K times exist for them
`df.j`, `df.k`, `cosx.k`	Separate DF and COSX kernels, nested under the driver phase that called them
`fock.j`, `fock.k`, `xc`	RKS Coulomb, exchange and XC builds. `xc` includes AO collocation, density/XC evaluation and Vxc assembly, which run batched together. PBE has no exchange phase
`accelerator`, `diagonalization`, `second_order_step`, `density`	DIIS-family extrapolation, Fock diagonalization or a second-order orbital step, and density construction
`final_rebuild`, `stability`	The post-convergence Fock rebuild and the stability analysis

ct.iterations repeats the phase table per SCF iteration. An iteration runs from one loop start to the next, so it includes that iteration’s diagonalization. For staged COSX grids, segment numbers the grid stages. ct.memory holds tracked DF container sizes at the end of construction: df.live_tracked_bytes is the retained tensors and df.construction_tracked_peak_bytes adds the temporary ones alive at that point. Neither includes BLAS or solver scratch, and neither is process RSS.

Timers run only on the thread that called the driver. The G0 CPU profiling harness in examples/regression/gpu_profile.py combines these records with runtime provenance and process RSS; see the G0 plan.

Bounded CPU profiling campaigns¶

The G0 harness can pin stability_check=False on every timed route with submit --no-stability-check. Add --measure-stability for one separate fresh-process driver run per case at that system’s largest thread count. That run requests stability explicitly and reports the native stability_wall_s separately from its total driver time. Unsupported stability responses remain refusals in batch_status.json; they are not zero-cost measurements or successful stability verdicts. These options do not change the production driver’s default stability policy.

Use --system-settings with an external JSON file mapping each system to its scheduling entries. Entries inherit --warmup, --repeats and --per-run-budget-s; threads is required. For example, the large-system campaign can use:

{
  "n-hexadecane": [
    {"threads": 4, "warmup": 1, "repeats": 3, "per_run_budget_s": 3600},
    {"threads": 1, "warmup": 0, "repeats": 1, "per_run_budget_s": 7200}
  ]
}

--compare-timing-off --timing-off-cases CASE_ID,... limits paired runs to selected cases, such as one naphthalene method per resolved route. Without --timing-off-cases, the comparison applies to every selected case. AUTO deduplication occurs inside the measured runtime; the worker writes batch_resolved.json with the actual cases, aliases and run schedule.

Submit with --python pointing to an immutable runtime launcher, and pair --runtime-root with --runtime-artifact to pin its deployed vq bundle. --mem-mb reserves the job’s memory through vq. Installed wheels without Git metadata are identified by the deployed bundle’s source revision, manifest digest and verification receipt. The harness checks the installed package inventory, native module and interpreter against the manifest both before and after each calculation, outside the measured SCF interval. Raw provenance stays in the external evidence root.

After fetching the job, summarize with --batch-json pointing to its batch_resolved.json, plus --expected-source-sha and --expected-native-sha256. Missing planned runs remain incomplete; separate stability runs cannot enter the SCF repeat or instrumentation comparison. Keep all campaign inputs, raw output and operational configuration outside the checkout. Publish only the redacted, measured summary after checking its completeness and identity results.

Fixed-density CPU noise¶

examples/regression/gpu_noise.py measures CPU arithmetic noise before a replacement kernel is compared. Stage it with gpu_profile.py and a direct case from the G0 manifest in an external vq workspace. Use the same pinned immutable runtime for fixture export and replay. The worker commands are:

"$RUNTIME_PYTHON" gpu_noise.py export case.json --out-dir evidence/fixture
"$RUNTIME_PYTHON" gpu_noise.py replay evidence/fixture/fixture.json \
  --route df --threads 1 --repeats 3 --out-dir evidence/replay
"$RUNTIME_PYTHON" gpu_noise.py replay evidence/fixture/fixture.json \
  --route df --threads 2 --repeats 3 --out-dir evidence/replay
"$RUNTIME_PYTHON" gpu_noise.py summarize evidence --threads 1,2 \
  --expected-source-sha "$SOURCE_SHA" --expected-native-sha256 "$NATIVE_SHA256" \
  --json-out noise-summary.json

Set OMP_NUM_THREADS for each process and set OPENBLAS_NUM_THREADS, MKL_NUM_THREADS, VECLIB_MAXIMUM_THREADS and BLIS_NUM_THREADS to one before launching it. Reserve N CPUs and sufficient memory through vq; each replay invocation is a fresh process. Apply an explicit queue and subprocess wall budget appropriate to the case, as for the timing campaign. The accepted routes are direct, df, rijcosx and resolved auto.

Fixture export runs the direct SCF with the G0 recipe and stability disabled. It requires convergence and the pinned residual target, then captures the returned density of a second run capped at half the converged iteration count. It also exports SAD and two small opposite symmetric perturbations of the converged density. The perturbations preserve Tr(D S) and have Frobenius norm 1e-7 * max(1, ||D||_F); they are not projected onto positive semidefinite matrices. The native restart normalizer prepares the frozen inputs once; its changes are recorded. The density path is retained throughout replay. Archive and matrix hashes keep these exact inputs fixed across routes and thread counts. No GPU dependencies are needed.

Direct replay disables incremental caching and fixes the tight Schwarz threshold. COSX uses the FITTED final grid of the pinned level; its K is the iterated exchange operator, excluding the post-SCF one-center correction. The latter still needs a separate noise qualification. PBE has no K sample. Vxc uses the native evaluate_rks_xc_potential interface when available. It returns the production restricted XC matrix, nonlinear energy and observed batch worker count at the supplied spin-summed density, without normalization or a restart check. Finite symmetric indefinite perturbations are accepted; the production pointwise density cutoffs still apply. VV10, when present, uses the supplied grid rather than a separate coarse grid. The G0 cases use PBE or PBE0 and do not include VV10.

Older runtimes fall back to the unpolarised RKS zero-iteration path with zero H, J and K, so the returned Fock matrix is Vxc without subtracting large terms. The harness checks that the supplied density is unchanged. The SCF restart interface refuses indefinite perturbations; these remain explicit incomplete rows until replayed with a runtime containing the direct interface. The recorded evaluator identity prevents combining the two paths as one noise measurement.

The summary reports worst pairwise matrix max-absolute errors in Hartree, relative Frobenius errors, symmetry defects, and energy differences within each thread count and between OMP=1 and N. J and K use Tr(D J)/2 and -a_x Tr(D K)/4, with RHF a_x=1, PBE0 a_x=1/4 and PBE a_x=0. For Vxc, Tr(D delta Vxc) is a separate diagnostic; the energy comparison uses the nonlinear native E_xc. A zero reference norm yields relative error zero only if both matrices vanish, otherwise JSON null.

MEASURED means the recorded input has at least three valid repeats at both requested thread counts. It does not certify coverage of other systems, methods, routes, spin conventions or hardware, or set a scientific tolerance. Missing thread samples remain INCOMPLETE; mismatched input, runtime, settings or sleep evidence are refused. Publish only redacted summaries. The G0 operation contracts must combine these measured floors with explicit scientific budgets before admitting a replacement implementation.

References¶

Almlöf, Fægri, Korsell, J. Comput. Chem. 3, 385 (1982), doi:10.1002/jcc.540030314. Difference-Fock recurrence (Eq. 20) and the final energy-oriented iteration recommendation (p. 389).
Häser & Ahlrichs, J. Comput. Chem. 10, 104 (1989), doi:10.1002/jcc.540100111. Cauchy-Schwarz screening and inherited screened-increment errors (Eqs. 30-35).

These references are cited automatically in the .bibtex and .references output files when the requested mode is DIRECT or when AUTO resolves to DIRECT.