!> @brief Cross-iteration cache of the collocation matrix Phi (AO values and !> derivatives on the quadrature grid) used by the DFT XC build. !> !> @details The collocation block Phi[mu,i] = phi_mu(r_i) (and its grid !> derivatives), together with the grid weights and the per-slice !> significant-AO pruning metadata, depends ONLY on geometry + basis + grid !> + integration thresholds -- NOT on the electron density. It is therefore !> identical on every SCF iteration, yet the native grid loop recomputes it !> (compAOs / pruneAOs) on every Fock build. This module stores the !> post-pruning Phi block per grid slice once per geometry and replays it on !> subsequent iterations, mirroring the incremental reuse OpenQP already does !> for the J/K (HF) part of the Fock matrix. !> !> The cache is OPT-IN (env var OQP_XC_PHI_CACHE) and is only populated by the !> repeated SCF energy/Fock build (mod_dft_gridint_energy::dmatd_blk); one-shot !> consumers (gradients, response) never set the opt-in flag. Validity is keyed !> by a geometry hash + derivative order + DFT threshold, so any change !> transparently rebuilds the cache. The replayed Phi is bit-for-bit identical !> to the recomputed Phi, so the converged energy/gradient are unchanged. !> !> @author Claude (Anthropic), 2026 module mod_dft_gridint_phi_cache use precision, only: fp, i8b implicit none private !> @brief Cached collocation data for a single grid slice type :: phi_slice_t logical :: stored = .false. !< this slice was populated during the build pass logical :: skip = .false. !< slice contributes nothing (empty or fully pruned) integer :: numPts = 0 !< number of significant grid points in the slice integer :: numAOs_p = 0 !< effective leading AO dimension (pruned or full) logical :: skip_p = .true. !< .true. if AO pruning was skipped (dense slice) integer, allocatable :: indices_p(:) !< significant AO indices (numAOs_p) real(fp), allocatable :: ao(:) !< compacted Phi block (numAOs_p*numPts*numAOVecs) real(fp), allocatable :: wts(:) !< grid weights (numPts) end type !> @brief Geometry-keyed cache of the collocation matrix across SCF iterations type, public :: phi_cache_t logical :: active = .false. !< this run participates in caching logical :: replay = .false. !< .true. => load cached Phi; .false. => build it logical :: ready = .false. !< cache is fully populated and valid integer :: nSlices = 0 integer :: numAOVecs = 0 integer :: numAOs = 0 ! validity signature integer(i8b) :: sig_geom = 0_i8b integer :: sig_vecs = 0 integer :: sig_natoms = 0 integer :: sig_nmolpts = 0 !< grid identity: total quadrature points real(fp) :: sig_thr = -1.0_fp integer(i8b) :: nbytes = 0_i8b type(phi_slice_t), allocatable :: slices(:) contains procedure :: begin_run procedure :: finish_run procedure :: store procedure :: get_meta procedure :: get_bulk procedure :: free end type !> Module-global singleton: persists across run_xc calls / SCF iterations type(phi_cache_t), public, save :: g_phi_cache public :: phi_cache_geom_hash public :: phi_cache_env_enabled contains !> @brief Deterministic 64-bit hash of the molecular geometry (FNV-1a over the !> raw IEEE bits of every coordinate). Any coordinate change invalidates. pure function phi_cache_geom_hash(xyz) result(h) real(fp), intent(in) :: xyz(:,:) integer(i8b) :: h integer(i8b) :: bits integer :: i, j h = -3750763034362895579_i8b ! FNV-1a 64-bit offset basis (1469598103934665603 as signed) do j = 1, size(xyz, 2) do i = 1, size(xyz, 1) bits = transfer(real(xyz(i,j), fp), 1_i8b) h = ieor(h, bits) h = h * 1099511628211_i8b ! FNV-1a 64-bit prime (modular wraparound is fine) end do end do end function !> @brief Is the Phi cache enabled via the environment? (cached after first read) function phi_cache_env_enabled() result(en) logical :: en logical, save :: known = .false. logical, save :: value = .false. character(len=32) :: s integer :: st, ln if (.not. known) then call get_environment_variable('OQP_XC_PHI_CACHE', s, length=ln, status=st) if (st == 0 .and. ln > 0) then s = adjustl(s) value = (trim(s) == '1' .or. trim(s) == 'true' .or. trim(s) == 'TRUE' & .or. trim(s) == 'on' .or. trim(s) == 'ON' .or. trim(s) == 'yes' & .or. trim(s) == 'YES' .or. trim(s) == 'T' .or. trim(s) == 't') end if known = .true. end if en = value end function !> @brief Decide build/replay mode for the current run and (re)allocate as needed. !> @param[in] enable master on/off (env AND opt-in caller) !> @param[in] nSlices number of grid slices !> @param[in] nMolPts total quadrature points (grid-identity discriminator) !> @param[in] numAOs number of basis functions !> @param[in] numAOVecs AO vectors per point (1=val, 4=+grad, 10=+2nd der) !> @param[in] natoms number of atoms !> @param[in] geom_hash geometry hash (see phi_cache_geom_hash) !> @param[in] thr DFT integration threshold (affects nonzero-point set) !> !> @note nMolPts is an explicit grid-identity key so a mid-run grid swap (e.g. !> the coarse->fine XC ramp, OQP_XC_C2F) can never be mistaken for a cache hit: !> with a user-fixed dft_threshold the thr/AO/atom/geom signatures are identical !> across the coarse and fine grids, leaving nSlices as the only other guard; !> the two grids' point counts differ even if their slice counts collide. subroutine begin_run(self, enable, nSlices, nMolPts, numAOs, numAOVecs, natoms, geom_hash, thr) class(phi_cache_t), intent(inout) :: self logical, intent(in) :: enable integer, intent(in) :: nSlices, nMolPts, numAOs, numAOVecs, natoms integer(i8b), intent(in) :: geom_hash real(fp), intent(in) :: thr logical :: match self%active = enable self%replay = .false. if (.not. enable) return match = self%ready & .and. self%nSlices == nSlices & .and. self%sig_nmolpts == nMolPts & .and. self%numAOs == numAOs & .and. self%numAOVecs == numAOVecs & .and. self%sig_natoms == natoms & .and. self%sig_geom == geom_hash & .and. self%sig_thr == thr if (match) then self%replay = .true. return end if ! (Re)build: drop any stale cache and prepare fresh per-slice storage. call self%free() self%nSlices = nSlices self%sig_nmolpts = nMolPts self%numAOs = numAOs self%numAOVecs = numAOVecs self%sig_natoms = natoms self%sig_geom = geom_hash self%sig_thr = thr allocate(self%slices(nSlices)) self%ready = .false. self%replay = .false. end subroutine !> @brief Finalize a build pass: mark the cache ready and tally its footprint. subroutine finish_run(self) class(phi_cache_t), intent(inout) :: self integer :: i if (.not. self%active) return if (self%replay) return self%nbytes = 0_i8b do i = 1, self%nSlices if (allocated(self%slices(i)%ao)) self%nbytes = self%nbytes + 8_i8b*size(self%slices(i)%ao, kind=i8b) if (allocated(self%slices(i)%wts)) self%nbytes = self%nbytes + 8_i8b*size(self%slices(i)%wts, kind=i8b) end do self%ready = .true. end subroutine !> @brief Store one slice's post-pruning collocation block (build pass). !> @param[in] iSlice slice index !> @param[in] skip whether the slice contributes nothing !> @param[in] numPts significant points in the slice !> @param[in] numAOs_p effective leading AO dimension !> @param[in] skip_p whether pruning was skipped !> @param[in] indices_p significant AO indices (length >= numAOs_p) !> @param[in] ao flat AO memory (xce%aoMem_); first numAOs_p*numPts*numAOVecs used !> @param[in] wts grid weights (length >= numPts) subroutine store(self, iSlice, skip, numPts, numAOs_p, skip_p, indices_p, ao, wts) class(phi_cache_t), intent(inout) :: self integer, intent(in) :: iSlice, numPts, numAOs_p logical, intent(in) :: skip, skip_p integer, intent(in) :: indices_p(:) real(fp), intent(in) :: ao(:) real(fp), intent(in) :: wts(:) integer :: n associate (s => self%slices(iSlice)) s%stored = .true. s%skip = skip if (skip) return s%numPts = numPts s%numAOs_p = numAOs_p s%skip_p = skip_p n = numAOs_p*numPts*self%numAOVecs ! Significant-AO indices are only consumed by the pruned (sparse) scatter ! path; the dense path adds the full block, so skip storing them there. if (.not. skip_p) s%indices_p = indices_p(1:numAOs_p) s%ao = ao(1:n) s%wts = wts(1:numPts) end associate end subroutine !> @brief Read a slice's scalar metadata (replay pass, phase 1). subroutine get_meta(self, iSlice, skip, numPts, numAOs_p, skip_p) class(phi_cache_t), intent(in) :: self integer, intent(in) :: iSlice logical, intent(out) :: skip, skip_p integer, intent(out) :: numPts, numAOs_p associate (s => self%slices(iSlice)) if (.not. s%stored) then ! defensive: treat unpopulated slice as empty skip = .true.; numPts = 0; numAOs_p = 0; skip_p = .true. return end if skip = s%skip numPts = s%numPts numAOs_p = s%numAOs_p skip_p = s%skip_p end associate end subroutine !> @brief Copy a slice's bulk data back into the engine arrays (replay, phase 2). !> @param[out] indices_p receives significant AO indices (first numAOs_p) !> @param[out] ao receives the Phi block into xce%aoMem_ (first n) !> @param[out] wts receives grid weights into xce%xyzw(:,4) (first numPts) subroutine get_bulk(self, iSlice, indices_p, ao, wts) class(phi_cache_t), intent(in) :: self integer, intent(in) :: iSlice integer, intent(out) :: indices_p(:) real(fp), intent(out) :: ao(:) real(fp), intent(out) :: wts(:) integer :: n, np associate (s => self%slices(iSlice)) np = s%numAOs_p n = np*s%numPts*self%numAOVecs if (.not. s%skip_p) indices_p(1:np) = s%indices_p(1:np) ao(1:n) = s%ao(1:n) wts(1:s%numPts) = s%wts(1:s%numPts) end associate end subroutine !> @brief Release all cached storage and reset to the empty/invalid state. subroutine free(self) class(phi_cache_t), intent(inout) :: self if (allocated(self%slices)) deallocate(self%slices) self%ready = .false. self%replay = .false. self%nSlices = 0 self%nbytes = 0_i8b end subroutine end module mod_dft_gridint_phi_cache