libxc 0.1.4

//! Defines the layout of output buffers for libxc computations, and utilities
//! to compute them based on functional properties and requested derivative
//! levels.

use crate::prelude::*;

/* #region derivative flags */

/// Flags controlling which derivative levels to compute.
///
/// This struct is somehow cumbersome; but you can also use `usize` levels
/// (0..=4) which are converted to the corresponding flags via `From<usize>`.
#[derive(Debug, Clone, Copy)]
pub struct LibXCDerivativeFlags {
    pub do_exc: bool,
    pub do_vxc: bool,
    pub do_fxc: bool,
    pub do_kxc: bool,
    pub do_lxc: bool,
}

impl Default for LibXCDerivativeFlags {
    fn default() -> Self {
        Self { do_exc: true, do_vxc: true, do_fxc: false, do_kxc: false, do_lxc: false }
    }
}

impl From<usize> for LibXCDerivativeFlags {
    /// Converts a derivative level (0..=4) to the corresponding flags.
    ///
    /// - Level 0: EXC
    /// - Level 1: EXC + VXC
    /// - level 2: EXC + VXC + FXC
    /// - level 3: EXC + VXC + FXC + KXC
    /// - level 4: EXC + VXC + FXC + KXC + LXC
    fn from(level: usize) -> Self {
        match level {
            0 => Self { do_exc: true, do_vxc: false, do_fxc: false, do_kxc: false, do_lxc: false },
            1 => Self { do_exc: true, do_vxc: true, do_fxc: false, do_kxc: false, do_lxc: false },
            2 => Self { do_exc: true, do_vxc: true, do_fxc: true, do_kxc: false, do_lxc: false },
            3 => Self { do_exc: true, do_vxc: true, do_fxc: true, do_kxc: true, do_lxc: false },
            4 => Self { do_exc: true, do_vxc: true, do_fxc: true, do_kxc: true, do_lxc: true },
            _ => panic!("invalid derivative level {level}"),
        }
    }
}

/* #endregion */

/* #region Output label definitions (mirroring pylibxc output_labels) */

#[rustfmt::skip]
pub const LDA_OUTPUT_LABELS: [&str; 5] = [
    "zk",                                                               // 1
    "vrho",                                                             // 1
    "v2rho2",                                                           // 1
    "v3rho3",                                                           // 1
    "v4rho4",                                                           // 1
];
const LDA_EXC_END: usize = 1;
const LDA_VXC_END: usize = 2;
const LDA_FXC_END: usize = 3;
const LDA_KXC_END: usize = 4;
const LDA_LXC_END: usize = 5;

#[rustfmt::skip]
pub const GGA_OUTPUT_LABELS: [&str; 15] = [
    "zk",                                                               // 1
    "vrho", "vsigma",                                                   // 2
    "v2rho2", "v2rhosigma", "v2sigma2",                                 // 3
    "v3rho3", "v3rho2sigma", "v3rhosigma2", "v3sigma3",                 // 4
    "v4rho4", "v4rho3sigma", "v4rho2sigma2", "v4rhosigma3", "v4sigma4", // 5
];
const GGA_EXC_END: usize = 1;
const GGA_VXC_END: usize = 3;
const GGA_FXC_END: usize = 6;
const GGA_KXC_END: usize = 10;
const GGA_LXC_END: usize = 15;

#[rustfmt::skip]
pub const MGGA_OUTPUT_LABELS: [&str; 70] = [
    "zk",                                                               // 1
    "vrho", "vsigma", "vlapl", "vtau",                                  // 4
    "v2rho2", "v2rhosigma", "v2rholapl", "v2rhotau", "v2sigma2",        // 10
    "v2sigmalapl", "v2sigmatau", "v2lapl2", "v2lapltau", "v2tau2",
    "v3rho3", "v3rho2sigma", "v3rho2lapl", "v3rho2tau", "v3rhosigma2",  // 20
    "v3rhosigmalapl", "v3rhosigmatau", "v3rholapl2", "v3rholapltau",
    "v3rhotau2", "v3sigma3", "v3sigma2lapl", "v3sigma2tau",
    "v3sigmalapl2", "v3sigmalapltau", "v3sigmatau2", "v3lapl3",
    "v3lapl2tau", "v3lapltau2", "v3tau3",
    "v4rho4", "v4rho3sigma", "v4rho3lapl", "v4rho3tau", "v4rho2sigma2",  // 35
    "v4rho2sigmalapl", "v4rho2sigmatau", "v4rho2lapl2", "v4rho2lapltau",
    "v4rho2tau2", "v4rhosigma3", "v4rhosigma2lapl", "v4rhosigma2tau",
    "v4rhosigmalapl2", "v4rhosigmalapltau", "v4rhosigmatau2",
    "v4rholapl3", "v4rholapl2tau", "v4rholapltau2", "v4rhotau3",
    "v4sigma4", "v4sigma3lapl", "v4sigma3tau", "v4sigma2lapl2",
    "v4sigma2lapltau", "v4sigma2tau2", "v4sigmalapl3", "v4sigmalapl2tau",
    "v4sigmalapltau2", "v4sigmatau3", "v4lapl4", "v4lapl3tau",
    "v4lapl2tau2", "v4lapltau3", "v4tau4",
];
const MGGA_EXC_END: usize = 1;
const MGGA_VXC_END: usize = 5;
const MGGA_FXC_END: usize = 15;
const MGGA_KXC_END: usize = 35;
const MGGA_LXC_END: usize = 70;

/// Maps an output label name to its dimension value from `xc_dimensions`.
pub fn get_dim(dim: &ffi::xc_dimensions, label: &str) -> i32 {
    match label {
        "zk" => dim.zk,
        "vrho" => dim.vrho,
        "vsigma" => dim.vsigma,
        "vlapl" => dim.vlapl,
        "vtau" => dim.vtau,
        "v2rho2" => dim.v2rho2,
        "v2rhosigma" => dim.v2rhosigma,
        "v2rholapl" => dim.v2rholapl,
        "v2rhotau" => dim.v2rhotau,
        "v2sigma2" => dim.v2sigma2,
        "v2sigmalapl" => dim.v2sigmalapl,
        "v2sigmatau" => dim.v2sigmatau,
        "v2lapl2" => dim.v2lapl2,
        "v2lapltau" => dim.v2lapltau,
        "v2tau2" => dim.v2tau2,
        "v3rho3" => dim.v3rho3,
        "v3rho2sigma" => dim.v3rho2sigma,
        "v3rho2lapl" => dim.v3rho2lapl,
        "v3rho2tau" => dim.v3rho2tau,
        "v3rhosigma2" => dim.v3rhosigma2,
        "v3rhosigmalapl" => dim.v3rhosigmalapl,
        "v3rhosigmatau" => dim.v3rhosigmatau,
        "v3rholapl2" => dim.v3rholapl2,
        "v3rholapltau" => dim.v3rholapltau,
        "v3rhotau2" => dim.v3rhotau2,
        "v3sigma3" => dim.v3sigma3,
        "v3sigma2lapl" => dim.v3sigma2lapl,
        "v3sigma2tau" => dim.v3sigma2tau,
        "v3sigmalapl2" => dim.v3sigmalapl2,
        "v3sigmalapltau" => dim.v3sigmalapltau,
        "v3sigmatau2" => dim.v3sigmatau2,
        "v3lapl3" => dim.v3lapl3,
        "v3lapl2tau" => dim.v3lapl2tau,
        "v3lapltau2" => dim.v3lapltau2,
        "v3tau3" => dim.v3tau3,
        "v4rho4" => dim.v4rho4,
        "v4rho3sigma" => dim.v4rho3sigma,
        "v4rho3lapl" => dim.v4rho3lapl,
        "v4rho3tau" => dim.v4rho3tau,
        "v4rho2sigma2" => dim.v4rho2sigma2,
        "v4rho2sigmalapl" => dim.v4rho2sigmalapl,
        "v4rho2sigmatau" => dim.v4rho2sigmatau,
        "v4rho2lapl2" => dim.v4rho2lapl2,
        "v4rho2lapltau" => dim.v4rho2lapltau,
        "v4rho2tau2" => dim.v4rho2tau2,
        "v4rhosigma3" => dim.v4rhosigma3,
        "v4rhosigma2lapl" => dim.v4rhosigma2lapl,
        "v4rhosigma2tau" => dim.v4rhosigma2tau,
        "v4rhosigmalapl2" => dim.v4rhosigmalapl2,
        "v4rhosigmalapltau" => dim.v4rhosigmalapltau,
        "v4rhosigmatau2" => dim.v4rhosigmatau2,
        "v4rholapl3" => dim.v4rholapl3,
        "v4rholapl2tau" => dim.v4rholapl2tau,
        "v4rholapltau2" => dim.v4rholapltau2,
        "v4rhotau3" => dim.v4rhotau3,
        "v4sigma4" => dim.v4sigma4,
        "v4sigma3lapl" => dim.v4sigma3lapl,
        "v4sigma3tau" => dim.v4sigma3tau,
        "v4sigma2lapl2" => dim.v4sigma2lapl2,
        "v4sigma2lapltau" => dim.v4sigma2lapltau,
        "v4sigma2tau2" => dim.v4sigma2tau2,
        "v4sigmalapl3" => dim.v4sigmalapl3,
        "v4sigmalapl2tau" => dim.v4sigmalapl2tau,
        "v4sigmalapltau2" => dim.v4sigmalapltau2,
        "v4sigmatau3" => dim.v4sigmatau3,
        "v4lapl4" => dim.v4lapl4,
        "v4lapl3tau" => dim.v4lapl3tau,
        "v4lapl2tau2" => dim.v4lapl2tau2,
        "v4lapltau3" => dim.v4lapltau3,
        "v4tau4" => dim.v4tau4,
        _ => 0,
    }
}

/* #endregion */

/* #region layout handling main part */

/// Mirrors pylibxc's `_check_arrays`: iterates output labels in
/// `labels[start..end]`, pushing those whose derivative level is
/// required and whose lapl/tau needs are satisfied.
#[allow(clippy::too_many_arguments)]
pub(crate) fn check_arrays(
    layout: &mut LibXCOutputLayout,
    labels: &[&'static str],
    start: usize,
    end: usize,
    dim: &ffi::xc_dimensions,
    required: bool,
    needs_lapl: bool,
    needs_tau: bool,
) {
    for &label in &labels[start..end] {
        let label_required = required
            && !(!needs_lapl && label.contains("lapl"))
            && !(!needs_tau && label.contains("tau"));
        if label_required {
            layout.push(label, get_dim(dim, label));
        }
    }
}

/// Describes the layout of a contiguous compute output buffer.
///
/// Each component (e.g. `"zk"`, `"vrho"`) occupies a contiguous range
/// `[offset .. offset + size)` within the buffer, stored in row-major order
/// `[n_comp, npoints]` where the last dimension is contiguous.
#[derive(Debug, Clone)]
pub struct LibXCOutputLayout {
    /// Total size of the buffer in f64 elements.
    pub total_size: usize,
    /// Number of grid points.
    pub npoints: usize,
    /// List of components in the layout (name, offset, n_comp).
    pub components: IndexMap<&'static str, (usize, usize)>,
}

impl LibXCOutputLayout {
    fn new(npoints: usize) -> Self {
        Self { total_size: 0, npoints, components: IndexMap::new() }
    }

    fn push(&mut self, name: &'static str, n_comp: i32) {
        if n_comp == 0 {
            return;
        }
        let size = self.npoints * (n_comp as usize);
        let offset = self.total_size;
        self.components.insert(name, (offset, n_comp as usize));
        self.total_size += size;
    }

    /// Returns the byte range `[offset .. offset + size)` for a named
    /// component.
    pub fn get(&self, name: &str) -> Option<Range<usize>> {
        self.components.get(name).map(|&(offset, n_comp)| offset..offset + n_comp * self.npoints)
    }

    /// Iterates over component names in order.
    pub fn component_names(&self) -> impl Iterator<Item = &str> {
        self.components.keys().copied()
    }

    /// Returns the size (in f64 elements) for a named component.
    pub fn component_size(&self, name: &str) -> Option<usize> {
        self.components.get(name).map(|&(_, n_comp)| n_comp * self.npoints)
    }

    /// Returns the number of components (not grid points) for a named
    /// component.
    pub fn component_dim(&self, name: &str) -> Option<usize> {
        self.components.get(name).map(|&(_, n_comp)| n_comp)
    }

    /// Iterates over components as (name, range) pairs.
    pub fn iter_to_range(&self) -> impl Iterator<Item = (&str, Range<usize>)> {
        self.components
            .iter()
            .map(|(&name, &(offset, n_comp))| (name, offset..offset + n_comp * self.npoints))
    }
}

impl core::fmt::Display for LibXCOutputLayout {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        writeln!(
            f,
            "LibXCOutputLayout (total_size={}, npoints={}):",
            self.total_size, self.npoints
        )?;
        for (name, (offset, n_comp)) in self.components.iter() {
            writeln!(
                f,
                "  {name}: [{offset}..{}], size={}, dim={}",
                offset + n_comp * self.npoints,
                n_comp * self.npoints,
                n_comp
            )?;
        }
        Ok(())
    }
}

impl LibXCFunctional {
    /// Compute the output layout for LDA at a given number of grid points.
    #[rustfmt::skip]
    pub fn lda_output_layout(
        &self,
        npoints: usize,
        flags: LibXCDerivativeFlags,
    ) -> LibXCOutputLayout {
        let dim = self.dim();
        let mut layout = LibXCOutputLayout::new(npoints);
        check_arrays(&mut layout, &LDA_OUTPUT_LABELS, 0, LDA_EXC_END, dim, flags.do_exc, false, false);
        check_arrays(&mut layout, &LDA_OUTPUT_LABELS, LDA_EXC_END, LDA_VXC_END, dim, flags.do_vxc, false, false);
        check_arrays(&mut layout, &LDA_OUTPUT_LABELS, LDA_VXC_END, LDA_FXC_END, dim, flags.do_fxc, false, false);
        check_arrays(&mut layout, &LDA_OUTPUT_LABELS, LDA_FXC_END, LDA_KXC_END, dim, flags.do_kxc, false, false);
        check_arrays(&mut layout, &LDA_OUTPUT_LABELS, LDA_KXC_END, LDA_LXC_END, dim, flags.do_lxc, false, false);
        layout
    }

    /// Compute the output layout for GGA at a given number of grid points.
    #[rustfmt::skip]
    pub fn gga_output_layout(
        &self,
        npoints: usize,
        flags: LibXCDerivativeFlags,
    ) -> LibXCOutputLayout {
        let dim = self.dim();
        let mut layout = LibXCOutputLayout::new(npoints);
        check_arrays(&mut layout, &GGA_OUTPUT_LABELS, 0, GGA_EXC_END, dim, flags.do_exc, false, false);
        check_arrays(&mut layout, &GGA_OUTPUT_LABELS, GGA_EXC_END, GGA_VXC_END, dim, flags.do_vxc, false, false);
        check_arrays(&mut layout, &GGA_OUTPUT_LABELS, GGA_VXC_END, GGA_FXC_END, dim, flags.do_fxc, false, false);
        check_arrays(&mut layout, &GGA_OUTPUT_LABELS, GGA_FXC_END, GGA_KXC_END, dim, flags.do_kxc, false, false);
        check_arrays(&mut layout, &GGA_OUTPUT_LABELS, GGA_KXC_END, GGA_LXC_END, dim, flags.do_lxc, false, false);
        layout
    }

    /// Compute the output layout for MGGA at a given number of grid points.
    #[rustfmt::skip]
    pub fn mgga_output_layout(
        &self,
        npoints: usize,
        flags: LibXCDerivativeFlags,
    ) -> LibXCOutputLayout {
        let dim = self.dim();
        let needs_lapl = self.needs_laplacian();
        let needs_tau = self.needs_tau();
        let mut layout = LibXCOutputLayout::new(npoints);
        check_arrays(&mut layout, &MGGA_OUTPUT_LABELS, 0, MGGA_EXC_END, dim, flags.do_exc, needs_lapl, needs_tau);
        check_arrays(&mut layout, &MGGA_OUTPUT_LABELS, MGGA_EXC_END, MGGA_VXC_END, dim, flags.do_vxc, needs_lapl, needs_tau);
        check_arrays(&mut layout, &MGGA_OUTPUT_LABELS, MGGA_VXC_END, MGGA_FXC_END, dim, flags.do_fxc, needs_lapl, needs_tau);
        check_arrays(&mut layout, &MGGA_OUTPUT_LABELS, MGGA_FXC_END, MGGA_KXC_END, dim, flags.do_kxc, needs_lapl, needs_tau);
        check_arrays(&mut layout, &MGGA_OUTPUT_LABELS, MGGA_KXC_END, MGGA_LXC_END, dim, flags.do_lxc, needs_lapl, needs_tau);
        layout
    }

    /// Compute the output layout for this functional at a given number of grid
    /// points.
    pub fn output_layout(
        &self,
        npoints: usize,
        flags: impl Into<LibXCDerivativeFlags>,
    ) -> LibXCOutputLayout {
        use crate::prelude::libxc_enum_items::*;
        let flags = flags.into();
        match self.family() {
            LDA | HybLDA => self.lda_output_layout(npoints, flags),
            GGA | HybGGA => self.gga_output_layout(npoints, flags),
            MGGA | HybMGGA => self.mgga_output_layout(npoints, flags),
            OEP | LCA => unimplemented!("output layout for OEP/LCA is not recognized."),
        }
    }

    pub(crate) fn validate_flags(
        &self,
        flags: impl Into<LibXCDerivativeFlags>,
    ) -> Result<(), LibXCError> {
        let flags = flags.into();
        for (flag, has_cap, name) in [
            (flags.do_exc, self.has_exc(), "EXC"),
            (flags.do_vxc, self.has_vxc(), "VXC"),
            (flags.do_fxc, self.has_fxc(), "FXC"),
            (flags.do_kxc, self.has_kxc(), "KXC"),
            (flags.do_lxc, self.has_lxc(), "LXC"),
        ] {
            if flag && !has_cap {
                return Err(LibXCError::ComputeError(format!(
                    "functional '{}' does not have {name} capabilities",
                    self.identifier(),
                )));
            }
        }
        Ok(())
    }
}

/* #endregion */

/* #region Shared unsafe FFI call functions (used by both CPU and CUDA paths) */

/// Helper to get a pointer from the validated map, defaulting to null.
pub(crate) fn ptr_of(ptrs: &HashMap<&'static str, *mut f64>, key: &str) -> *mut f64 {
    ptrs.get(key).copied().unwrap_or(std::ptr::null_mut())
}

/// Invoke `xc_lda` FFI call with a contiguous output buffer and layout.
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// `rho_ptr` points to at least `npoints * nspin` readable `f64`s, and
/// `output_base` points to at least `layout.total_size` writable `f64`s.
pub(crate) unsafe fn xc_lda_call(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    output_base: *mut f64,
    layout: &LibXCOutputLayout,
) {
    let ptr_for = |name: &str| -> *mut f64 {
        match layout.get(name) {
            Some(range) => output_base.add(range.start),
            None => std::ptr::null_mut::<f64>(),
        }
    };
    ffi::xc_lda(
        func_ptr,
        npoints,
        rho_ptr,
        ptr_for("zk"),
        ptr_for("vrho"),
        ptr_for("v2rho2"),
        ptr_for("v3rho3"),
        ptr_for("v4rho4"),
    );
}

/// Invoke `xc_lda` FFI call with named per-component output pointers.
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// `rho_ptr` points to at least `npoints * nspin` readable `f64`s, and each
/// non-null pointer in `ptrs` points to a buffer of the correct size.
pub(crate) unsafe fn xc_lda_call_with_output(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    ptrs: &HashMap<&'static str, *mut f64>,
) {
    ffi::xc_lda(
        func_ptr,
        npoints,
        rho_ptr,
        ptr_of(ptrs, "zk"),
        ptr_of(ptrs, "vrho"),
        ptr_of(ptrs, "v2rho2"),
        ptr_of(ptrs, "v3rho3"),
        ptr_of(ptrs, "v4rho4"),
    );
}

/// Invoke `xc_gga` FFI call with a contiguous output buffer and layout.
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// input pointers are valid for their expected sizes, and `output_base` points
/// to at least `layout.total_size` writable `f64`s.
pub(crate) unsafe fn xc_gga_call(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    sigma_ptr: *const f64,
    output_base: *mut f64,
    layout: &LibXCOutputLayout,
) {
    let ptr_for = |name: &str| -> *mut f64 {
        match layout.get(name) {
            Some(range) => output_base.add(range.start),
            None => std::ptr::null_mut::<f64>(),
        }
    };
    ffi::xc_gga(
        func_ptr,
        npoints,
        rho_ptr,
        sigma_ptr as *mut f64,
        ptr_for("zk"),
        ptr_for("vrho"),
        ptr_for("vsigma"),
        ptr_for("v2rho2"),
        ptr_for("v2rhosigma"),
        ptr_for("v2sigma2"),
        ptr_for("v3rho3"),
        ptr_for("v3rho2sigma"),
        ptr_for("v3rhosigma2"),
        ptr_for("v3sigma3"),
        ptr_for("v4rho4"),
        ptr_for("v4rho3sigma"),
        ptr_for("v4rho2sigma2"),
        ptr_for("v4rhosigma3"),
        ptr_for("v4sigma4"),
    );
}

/// Invoke `xc_gga` FFI call with named per-component output pointers.
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// input pointers are valid for their expected sizes, and each non-null
/// pointer in `ptrs` points to a buffer of the correct size.
pub(crate) unsafe fn xc_gga_call_with_output(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    sigma_ptr: *const f64,
    ptrs: &HashMap<&'static str, *mut f64>,
) {
    ffi::xc_gga(
        func_ptr,
        npoints,
        rho_ptr,
        sigma_ptr as *mut f64,
        ptr_of(ptrs, "zk"),
        ptr_of(ptrs, "vrho"),
        ptr_of(ptrs, "vsigma"),
        ptr_of(ptrs, "v2rho2"),
        ptr_of(ptrs, "v2rhosigma"),
        ptr_of(ptrs, "v2sigma2"),
        ptr_of(ptrs, "v3rho3"),
        ptr_of(ptrs, "v3rho2sigma"),
        ptr_of(ptrs, "v3rhosigma2"),
        ptr_of(ptrs, "v3sigma3"),
        ptr_of(ptrs, "v4rho4"),
        ptr_of(ptrs, "v4rho3sigma"),
        ptr_of(ptrs, "v4rho2sigma2"),
        ptr_of(ptrs, "v4rhosigma3"),
        ptr_of(ptrs, "v4sigma4"),
    );
}

/// Extra output pointers for mGGA labels absent from the layout or pointer
/// map.
///
/// libxc 6.2.x unconditionally requires every tau-family output buffer of a
/// requested derivative order to be non-NULL, even when the functional does
/// not use tau (the `XC_FLAGS_NEEDS_TAU` flag was only introduced in v7.0),
/// and it also reads the tau input array for every meta-GGA functional.
/// Passing null pointers there makes libxc abort the whole process
/// (`xc_mgga_sanity_check` calls `exit(1)`).
///
/// The CPU wrapper therefore allocates zeroed scratch buffers for those
/// labels and passes their pointers here. Device-side (CUDA) callers can
/// pass an empty map: libxc 6.2.x has no GPU support at all, and v7.0+ only
/// requires tau-family buffers when the functional actually needs them.
pub(crate) type MggaExtraPtrs = HashMap<&'static str, *mut f64>;

/// Build zeroed scratch buffers for the tau-family output labels of every
/// requested derivative order.
///
/// `present` reports whether an output label will be passed to libxc (either
/// inside the contiguous output buffer or as a dedicated pointer). A
/// tau-family label is only scratched when the primary output of its
/// derivative order is requested; libxc only validates tau-family buffers
/// for the orders that are actually computed.
fn mgga_tau_scratch_impl(
    present: &impl Fn(&str) -> bool,
    dim: &ffi::xc_dimensions,
    npoints: usize,
) -> (Vec<Vec<f64>>, MggaExtraPtrs) {
    const MGGA_ORDER_PRIMARY: [(&str, usize, usize); 5] = [
        ("zk", 0, MGGA_EXC_END),
        ("vrho", MGGA_EXC_END, MGGA_VXC_END),
        ("v2rho2", MGGA_VXC_END, MGGA_FXC_END),
        ("v3rho3", MGGA_FXC_END, MGGA_KXC_END),
        ("v4rho4", MGGA_KXC_END, MGGA_LXC_END),
    ];
    let mut bufs = Vec::new();
    let mut ptrs = HashMap::new();
    for (primary, start, end) in MGGA_ORDER_PRIMARY {
        if !present(primary) {
            continue;
        }
        for &label in &MGGA_OUTPUT_LABELS[start..end] {
            if label.contains("tau") && !present(label) {
                let d = get_dim(dim, label);
                if d > 0 {
                    bufs.push(vec![0.0f64; (d as usize) * npoints]);
                    ptrs.insert(label, bufs.last().unwrap().as_ptr() as *mut f64);
                }
            }
        }
    }
    (bufs, ptrs)
}

/// Build zeroed scratch buffers for tau-family outputs missing from a
/// contiguous output buffer layout.
pub(crate) fn mgga_tau_scratch(
    layout: &LibXCOutputLayout,
    dim: &ffi::xc_dimensions,
    npoints: usize,
) -> (Vec<Vec<f64>>, MggaExtraPtrs) {
    mgga_tau_scratch_impl(&|label| layout.get(label).is_some(), dim, npoints)
}

/// Build zeroed scratch buffers for tau-family outputs missing from a
/// named-pointer map.
///
/// The map may contain null pointers for labels the user did not provide
/// (see `validate_output_ptrs`); only non-null entries count as present.
pub(crate) fn mgga_tau_scratch_from_ptrs(
    ptrs: &HashMap<&'static str, *mut f64>,
    dim: &ffi::xc_dimensions,
    npoints: usize,
) -> (Vec<Vec<f64>>, MggaExtraPtrs) {
    mgga_tau_scratch_impl(&|label| ptrs.get(label).is_some_and(|p| !p.is_null()), dim, npoints)
}

/// Invoke `xc_mgga` FFI call with a contiguous output buffer and layout.
///
/// `extra` supplies output pointers for labels not covered by `layout`
/// (scratch buffers for libxc 6.2.x; see [`MggaExtraPtrs`]).
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// input pointers are valid for their expected sizes, and `output_base` points
/// to at least `layout.total_size` writable `f64`s. Pointers in `extra` must
/// be valid for `get_dim(dim, label) * npoints` writable `f64`s and outlive
/// the call.
#[allow(clippy::too_many_arguments)]
pub(crate) unsafe fn xc_mgga_call(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    sigma_ptr: *const f64,
    lapl_ptr: *const f64,
    tau_ptr: *const f64,
    output_base: *mut f64,
    layout: &LibXCOutputLayout,
    extra: &MggaExtraPtrs,
) {
    let ptr_for = |name: &str| -> *mut f64 {
        if let Some(range) = layout.get(name) {
            return output_base.add(range.start);
        }
        match extra.get(name) {
            Some(&ptr) => ptr,
            None => std::ptr::null_mut::<f64>(),
        }
    };

    ffi::xc_mgga(
        func_ptr,
        npoints,
        rho_ptr,
        sigma_ptr as *mut f64,
        lapl_ptr as *mut f64,
        tau_ptr as *mut f64,
        ptr_for("zk"),
        ptr_for("vrho"),
        ptr_for("vsigma"),
        ptr_for("vlapl"),
        ptr_for("vtau"),
        ptr_for("v2rho2"),
        ptr_for("v2rhosigma"),
        ptr_for("v2rholapl"),
        ptr_for("v2rhotau"),
        ptr_for("v2sigma2"),
        ptr_for("v2sigmalapl"),
        ptr_for("v2sigmatau"),
        ptr_for("v2lapl2"),
        ptr_for("v2lapltau"),
        ptr_for("v2tau2"),
        ptr_for("v3rho3"),
        ptr_for("v3rho2sigma"),
        ptr_for("v3rho2lapl"),
        ptr_for("v3rho2tau"),
        ptr_for("v3rhosigma2"),
        ptr_for("v3rhosigmalapl"),
        ptr_for("v3rhosigmatau"),
        ptr_for("v3rholapl2"),
        ptr_for("v3rholapltau"),
        ptr_for("v3rhotau2"),
        ptr_for("v3sigma3"),
        ptr_for("v3sigma2lapl"),
        ptr_for("v3sigma2tau"),
        ptr_for("v3sigmalapl2"),
        ptr_for("v3sigmalapltau"),
        ptr_for("v3sigmatau2"),
        ptr_for("v3lapl3"),
        ptr_for("v3lapl2tau"),
        ptr_for("v3lapltau2"),
        ptr_for("v3tau3"),
        ptr_for("v4rho4"),
        ptr_for("v4rho3sigma"),
        ptr_for("v4rho3lapl"),
        ptr_for("v4rho3tau"),
        ptr_for("v4rho2sigma2"),
        ptr_for("v4rho2sigmalapl"),
        ptr_for("v4rho2sigmatau"),
        ptr_for("v4rho2lapl2"),
        ptr_for("v4rho2lapltau"),
        ptr_for("v4rho2tau2"),
        ptr_for("v4rhosigma3"),
        ptr_for("v4rhosigma2lapl"),
        ptr_for("v4rhosigma2tau"),
        ptr_for("v4rhosigmalapl2"),
        ptr_for("v4rhosigmalapltau"),
        ptr_for("v4rhosigmatau2"),
        ptr_for("v4rholapl3"),
        ptr_for("v4rholapl2tau"),
        ptr_for("v4rholapltau2"),
        ptr_for("v4rhotau3"),
        ptr_for("v4sigma4"),
        ptr_for("v4sigma3lapl"),
        ptr_for("v4sigma3tau"),
        ptr_for("v4sigma2lapl2"),
        ptr_for("v4sigma2lapltau"),
        ptr_for("v4sigma2tau2"),
        ptr_for("v4sigmalapl3"),
        ptr_for("v4sigmalapl2tau"),
        ptr_for("v4sigmalapltau2"),
        ptr_for("v4sigmatau3"),
        ptr_for("v4lapl4"),
        ptr_for("v4lapl3tau"),
        ptr_for("v4lapl2tau2"),
        ptr_for("v4lapltau3"),
        ptr_for("v4tau4"),
    );
}

/// Invoke `xc_mgga` FFI call with named per-component output pointers.
///
/// `extra` supplies output pointers for labels not covered by `ptrs`
/// (scratch buffers for libxc 6.2.x; see [`MggaExtraPtrs`]).
///
/// # Safety
///
/// The caller must ensure that `func_ptr` is a valid `xc_func_type` pointer,
/// input pointers are valid for their expected sizes, and each non-null
/// pointer in `ptrs` points to a buffer of the correct size. Pointers in
/// `extra` must be valid for `get_dim(dim, label) * npoints` writable `f64`s
/// and outlive the call.
#[allow(clippy::too_many_arguments)]
pub(crate) unsafe fn xc_mgga_call_with_output(
    func_ptr: *mut ffi::xc_func_type,
    npoints: usize,
    rho_ptr: *const f64,
    sigma_ptr: *const f64,
    lapl_ptr: *const f64,
    tau_ptr: *const f64,
    ptrs: &HashMap<&'static str, *mut f64>,
    extra: &MggaExtraPtrs,
) {
    let ptr_of_either = |name: &str| -> *mut f64 {
        match ptrs.get(name) {
            Some(&ptr) if !ptr.is_null() => ptr,
            _ => extra.get(name).copied().unwrap_or(std::ptr::null_mut()),
        }
    };
    ffi::xc_mgga(
        func_ptr,
        npoints,
        rho_ptr,
        sigma_ptr as *mut f64,
        lapl_ptr as *mut f64,
        tau_ptr as *mut f64,
        ptr_of_either("zk"),
        ptr_of_either("vrho"),
        ptr_of_either("vsigma"),
        ptr_of_either("vlapl"),
        ptr_of_either("vtau"),
        ptr_of_either("v2rho2"),
        ptr_of_either("v2rhosigma"),
        ptr_of_either("v2rholapl"),
        ptr_of_either("v2rhotau"),
        ptr_of_either("v2sigma2"),
        ptr_of_either("v2sigmalapl"),
        ptr_of_either("v2sigmatau"),
        ptr_of_either("v2lapl2"),
        ptr_of_either("v2lapltau"),
        ptr_of_either("v2tau2"),
        ptr_of_either("v3rho3"),
        ptr_of_either("v3rho2sigma"),
        ptr_of_either("v3rho2lapl"),
        ptr_of_either("v3rho2tau"),
        ptr_of_either("v3rhosigma2"),
        ptr_of_either("v3rhosigmalapl"),
        ptr_of_either("v3rhosigmatau"),
        ptr_of_either("v3rholapl2"),
        ptr_of_either("v3rholapltau"),
        ptr_of_either("v3rhotau2"),
        ptr_of_either("v3sigma3"),
        ptr_of_either("v3sigma2lapl"),
        ptr_of_either("v3sigma2tau"),
        ptr_of_either("v3sigmalapl2"),
        ptr_of_either("v3sigmalapltau"),
        ptr_of_either("v3sigmatau2"),
        ptr_of_either("v3lapl3"),
        ptr_of_either("v3lapl2tau"),
        ptr_of_either("v3lapltau2"),
        ptr_of_either("v3tau3"),
        ptr_of_either("v4rho4"),
        ptr_of_either("v4rho3sigma"),
        ptr_of_either("v4rho3lapl"),
        ptr_of_either("v4rho3tau"),
        ptr_of_either("v4rho2sigma2"),
        ptr_of_either("v4rho2sigmalapl"),
        ptr_of_either("v4rho2sigmatau"),
        ptr_of_either("v4rho2lapl2"),
        ptr_of_either("v4rho2lapltau"),
        ptr_of_either("v4rho2tau2"),
        ptr_of_either("v4rhosigma3"),
        ptr_of_either("v4rhosigma2lapl"),
        ptr_of_either("v4rhosigma2tau"),
        ptr_of_either("v4rhosigmalapl2"),
        ptr_of_either("v4rhosigmalapltau"),
        ptr_of_either("v4rhosigmatau2"),
        ptr_of_either("v4rholapl3"),
        ptr_of_either("v4rholapl2tau"),
        ptr_of_either("v4rholapltau2"),
        ptr_of_either("v4rhotau3"),
        ptr_of_either("v4sigma4"),
        ptr_of_either("v4sigma3lapl"),
        ptr_of_either("v4sigma3tau"),
        ptr_of_either("v4sigma2lapl2"),
        ptr_of_either("v4sigma2lapltau"),
        ptr_of_either("v4sigma2tau2"),
        ptr_of_either("v4sigmalapl3"),
        ptr_of_either("v4sigmalapl2tau"),
        ptr_of_either("v4sigmalapltau2"),
        ptr_of_either("v4sigmatau3"),
        ptr_of_either("v4lapl4"),
        ptr_of_either("v4lapl3tau"),
        ptr_of_either("v4lapl2tau2"),
        ptr_of_either("v4lapltau3"),
        ptr_of_either("v4tau4"),
    );
}

/* #endregion */