API reference

This page collects the public Python objects that are useful when following the calculation workflow from code. The narrative chapters explain the physical ideas; the reference below documents the classes and methods that expose those ideas programmatically.

Public package interface

The top-level pydft package re-exports the commonly used classes and helpers, including pydft.DFT, pydft.MolecularGrid, pydft.AtomicGrid, pydft.Molecule, and pydft.MoleculeBuilder.

Public package interface for PyDFT.

SCF driver

class pydft.DFT(mol: pyqint.Molecule, basis: str | list[pyqint.CGF] = 'sto3g', functional: str = 'svwn5', nshells: Mapping[str, int] | None = None, nangpts: Mapping[str, int] | None = None, lmax: Mapping[str, int] | None = None, fdpts: int = 7, normalize: bool = True)

Bases: object

Coordinate a complete Kohn-Sham density-functional theory calculation.

The DFT object is the main entry point for users. It owns the molecular geometry, the Gaussian basis functions, the numerical molecular grid, and the matrices used in the self-consistent field (SCF) cycle. The implementation intentionally keeps the SCF steps visible: one-electron matrix construction, density-matrix formation, grid-based Hartree and exchange-correlation terms, Fock matrix diagonalization, and convergence acceleration are all represented by small methods.

Notes

Calling scf() returns a dictionary with matrices, energies, orbitals, and timing data. This makes the object useful both for running calculations and for inspecting the intermediate objects that appear in a DFT textbook.

get_construction_times() → dict

Return construct times dictionary from MolecularGrid to get insights into the time to construct particular properties

Returns:: dictionary of construction times for the MolecularGrid objects
Return type:: dict

get_data() → dict

Return results of the SCF calculation.

This method is only valid after a successful SCF run. It returns a dictionary containing all relevant physical quantities, matrices, energies, and timing information. The returned data layout is shared between PyQInt and PyDFT calculations to ensure consistent post-processing.

Returns:: Dictionary containing SCF results and metadata.
Return type:: dict

Notes

The returned dictionary contains the following entries:

System information

mol : Molecular object defining geometry and atoms
nuclei : Nuclear positions and charges
cgfs : Contracted Gaussian basis functions

Energies

energy : Final total electronic energy
energies : SCF energy history
ekin : Electronic kinetic energy
enuc : Electron-nuclear attraction energy
enucrep : Nuclear-nuclear repulsion energy
ex : Exchange energy
ec : Correlation energy
exc : Exchange-correlation energy

Orbital quantities

orbc : Molecular orbital coefficient matrix
C : Compatibility alias for orbc
orbe : Molecular orbital eigenvalues

Matrices and operators

overlap : Overlap matrix
kinetic : Kinetic energy matrix
nuclear : Nuclear attraction matrix
hcore : Core Hamiltonian matrix (T + V)
density : Density matrix
P : Compatibility alias for density
fock : Fock matrix
hartree : Hartree (Coulomb) matrix
J : Compatibility alias for hartree
xc : Exchange-correlation matrix (DFT only, None for HF)
XC : Compatibility alias for xc

Timing information

time_statsDictionary with timing breakdowns
- construct : Grid and setup construction times
- scf : SCF iteration timings

get_density_at_points(spoints: numpy.ndarray) → numpy.ndarray

Get the electron density at provided points

Parameters:: spoints (ndarray) – Set of sampling points (\(X \times 3\) array)
Returns:: Electron density scalar field (\(N \times 1\) array)
Return type:: ndarray

get_gradient_at_points(spoints: numpy.ndarray) → numpy.ndarray

Get the electron density gradient at provided points

Parameters:: spoints (ndarray) – Set of sampling points (\(X \times 3\) array)
Returns:: Electron density gradient vector field (\(N \times 3\) array)
Return type:: ndarray

get_molgrid_copy() → MolecularGrid

Get a copy of the underlying molgrid object

Returns:: copy of the MolecularGrid object
Return type:: MolecularGrid

print_time_statistics(): Print a summary of time statistics

scf(tol: float = 1e-05, verbose: bool = False) → dict

Perform the self-consistent field procedure.

The SCF cycle repeatedly builds the density-dependent parts of the Kohn-Sham matrix from the current density matrix, diagonalizes the resulting Fock/Kohn-Sham matrix, and updates the density matrix until the total energy changes by less than tol. Early iterations use linear mixing; later iterations use DIIS extrapolation when possible.

Parameters:

tol (float, optional) – Electronic energy convergence criterion in Hartree. The default is 1e-5.
verbose (bool, optional) – If True, log one line per SCF iteration showing the iteration number, total energy, energy change, and elapsed time. Configure the logging module to display these messages.

Returns:

Result dictionary as returned by get_data(). The final total energy is available as result['energy'].

Return type:

dict

Numerical grids

class pydft.MolecularGrid(atoms: list[Tuple[str, numpy.typing.NDArray.numpy.float64]], cgfs: list[pyqint.CGF], nshells: Mapping[str, int] | int | None = None, nangpts: Mapping[str, int] | int | None = None, lmax: Mapping[str, int] | int | None = None, fdpts: int = 7, functional: str = 'svwn5', interpolation_method='CubicStencil')

Bases: object

Combine atom-centered grids into a molecular integration grid.

MolecularGrid is the numerical workhorse behind the SCF driver. It constructs Becke fuzzy-cell weights, caches basis-function amplitudes and gradients on all grid points, evaluates electron densities from density matrices, builds the Hartree potential from atom-centered spherical-harmonic expansions, and assembles exchange-correlation matrix contributions.

The object is intentionally more transparent than a production quantum chemistry grid engine. Many getters expose intermediate arrays so examples and documentation can visualize the quadrature grid, Becke weights, electron density, molecular orbitals, and local potentials.

build_density(P: numpy.ndarray, normalize: bool = True)

Build the electron density from the density matrix

Parameters:

P (np.ndarray) – density matrix
normalize (bool, optional) – whether to normalize the density matrix based on the total number of electrons, by default True

calculate_correlation()

Build the correlation matrix and correlation energy for the molecule.

The same local-density and gradient-corrected structure used for exchange is applied to the correlation functional. The returned matrix is added to the exchange matrix by the SCF driver to form the full exchange-correlation contribution.

Returns:: Correlation matrix and correlation energy
Return type:: np.ndarray,float

calculate_coulomb_potential_at_points(pts: numpy.ndarray) → numpy.ndarray

Calculate the Hartree (Coulomb) potential at arbitrary points.

Parameters:: pts (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
Returns:: Coulomb potential at grid points \((N \times 1)\) with \(N\) the number of grid points.
Return type:: np.array

calculate_coulombic_matrix(calculate_timestats: bool = False) → numpy.array

Build the Coulomb matrix \(\mathbf{J}\) from the Hartree potential.

The electron density is represented on each atomic grid as spherical-harmonic coefficients. The corresponding Hartree-potential coefficients are interpolated from every atomic center to every point of the molecular Becke grid. The final matrix element is the quadrature of the Hartree potential multiplied by a pair of basis functions.

Returns:: Coulomb matrix \(\mathbf{J}\)
Return type:: np.array

calculate_dfa_coulomb() → float

Get density functional approximation for the coulombic repulsion between electrons

Returns:: Total coulombic repulsion (Hartree)
Return type:: float

calculate_dfa_coulomb_no_interpolation() → float

Calculate the Coulombic self-repulsion in the limit of the sum of atomic centers, without any interpolation of the contribution of the other atoms

Returns:: Approximate coulombic repulsion (Hartree)
Return type:: float

calculate_dfa_exchange()

Get density functional approximation for the exchange energy

Returns:: Total exchange energy (Hartree)
Return type:: float

calculate_dfa_kinetic() → float

Get density functional approximation for the kinetic energy

Returns:: Total kinetic energy (Hartree)
Return type:: float

calculate_dfa_nuclear_attraction_full() → numpy.array

Get density functional approximation for the nuclear attraction energy

Returns:: Total nuclear attraction energy (Hartree)
Return type:: float

calculate_dfa_nuclear_attraction_local() → float

Get density functional approximation for the nuclear attraction energy using only the local potential for each atomic grid and ignoring the influence of atomic attraction due to other atoms

Returns:: Approximate nuclear attraction energy (Hartree)
Return type:: float

Notes

This function will of course not give an accurate result and is only meant for educational purposes

calculate_exchange()

Build the exchange matrix and exchange energy for the molecule.

For LDA functionals, the matrix contribution is local in the density. For GGA functionals such as PBE, the derivative with respect to \(\sigma = |\nabla\rho|^2\) adds a second term involving basis function gradients.

Returns:: Exchange matrix \(\mathbf{X}\) and exchange energy
Return type:: np.ndarray,float

calculate_weights_at_points(points: numpy.ndarray, k: int = 3) → numpy.ndarray

Custom function that (re-)calculates weights at given points

Parameters:

points (np.ndarray) – \(N \times 3\) array of grid points
k (int, optional) – smoothing value, by default 3

Returns:

\(N_{a} \times N\) array with \(N_{a}\) the number of atoms and \(N\) the number of grid points

Return type:

np.ndarray

count_electrons() → float

Count number of electrons in the system

Returns:: Total number of electrons
Return type:: float

count_electrons_from_rho_lm() → float

Count number of electrons from the spherical harmonic expansion; this function mainly acts to verify that the spherical harmonic expansion works correctly

Returns:: Total number of electrons
Return type:: float

Notes

The Becke weights are already used when generating the density coefficients and hence for the back-transformation, we should not multiply with the Becke weights

get_amplitude_at_points(spoints: numpy.ndarray, c: numpy.ndarray) → numpy.ndarray

Calculate the wave function amplitude at the points as given by spoints using solution vector c

Parameters:

spoints (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
c (np.ndarray) – Coefficient vector \(\vec{c}\) of a single molecular orbital

Returns:

Molecular orbital amplitude specified at grid points \((N \times 1)\) with \(N\) the number of grid points.

Return type:

np.array

get_atomic_grid(i: int) → AtomicGrid

Get an atomic grid of atom \(i\)

Parameters:: i (int) – atom index
Returns:: pydft.AtomicGrid of atom \(i\)
Return type:: AtomicGrid

get_atomic_relative_grid_coordinates() → list[numpy.array]

Get the molecular grid coordinates relative to their corresponding atomic centers

Returns:: \((N \times 3)\) array with \(N\) number of grid points
Return type:: np.array

get_becke_weights() → numpy.ndarray

Get the Becke coefficients

Returns:: \((N \times G \times 3)\) array with \(N\) number of atoms and \(G\) number of grid points
Return type:: np.array

get_cached_amplitudes() → list[numpy.array]: Get the pre-cached basis function amplitudes for all grid points per atom

get_correlation_potential_at_points(pts: numpy.ndarray, P: numpy.ndarray) → numpy.ndarray

Get the correlation potential at points pts

Parameters:

pts (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
P (np.ndarray) – Density matrix \((K \times K)\) with \(K\) the number of basis functions

Returns:

Correlation potential \(\nu_{c}\) at grid points \((N \times 1)\) with \(N\) the number of grid points

Return type:

np.ndarray

get_densities() → list[numpy.array]

Get the densities at the grid points

Returns:: list over atoms with for each atom a \((G \times 1)\) array of electron densities where \(G\) is the number of grid points per atom
Return type:: list[np.array]

get_density_at_points(spoints: numpy.ndarray, P: numpy.ndarray) → numpy.array

Calculate the electron density at the points as given by spoints using density matrix P

Parameters:

spoints (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
P (np.ndarray) – Density matrix \((K \times K)\) with \(K\) the number of basis functions

Returns:

Electron density specified at grid points \((N \times 1)\) with \(N\) the number of grid points.

Return type:

np.array

get_exchange_potential_at_points(pts: numpy.ndarray, P: numpy.ndarray) → numpy.ndarray

Get the exchange potential at points pts

Parameters:

pts (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
P (np.ndarray) – Density matrix \((K \times K)\) with \(K\) the number of basis functions

Returns:

Exchange potential \(\nu_{x}\) at grid points \((N \times 1)\) with \(N\) the number of grid points

Return type:

np.ndarray

get_gradient_at_points(spoints: numpy.ndarray, P: numpy.ndarray)

Calculate the gradient of the electron density at the points as given by spoints using density matrix \(\mathbf{P}\)

Parameters:

spoints (np.ndarray) – Array of grid points \((N \times 3)\) with \(N\) the number of grid points
P (np.ndarray) – Density matrix \((K \times K)\) with \(K\) the number of basis functions

Returns:

Electron density gradients at grid points \((N \times 3)\) with \(N\) the number of grid points

Return type:

np.array

get_gradients() → list[numpy.array]

Get the gradient magnitude of the electron density field

Returns:: list over atoms with for each atom a \((G \times 1)\) array of electron gradient magnitudes where \(G\) is the number of grid points per atom
Return type:: list[np.array]

get_grid_coordinates() → list[numpy.array]

Get the grid coordinates

Returns:: list over atoms with for each atom a \((G \times 1)\) array of electron densities where \(G\) is the number of grid points per atom
Return type:: list[np.array]

get_hartree_potential() → numpy.array

Get the Hartree potential at each grid cell

Returns:: \((G \times 1)\) array with \(G\) number of molecular grid points
Return type:: np.array

get_molecular_grid_coordinates() → list[numpy.array]

Get the molecular grid coordinates

Returns:: \((N \times 3)\) array with \(N\) number of grid points
Return type:: np.array

get_rgridpoints() → list[numpy.array]

Get the radii relative to their corresponding atomic centers

Returns:: \((N \times 3)\) array with \(N\) number of grid points
Return type:: np.array

get_rho_lm_atoms() → np.ndarray | list[np.ndarray]

Get the electron density projected onto spherical harmonics per atom

Returns:: Per-atom arrays of \(\rho_{n,lm}\) coefficients. A single stacked array is returned when all atomic grids have the same shape; otherwise a list is returned.
Return type:: np.ndarray or list[np.ndarray]

get_spherical_harmonic_expansion_of_amplitude(c: numpy.ndarray, radial_factor=False) → list[numpy.ndarray]

Get the spherical harmonic expansion representation of a wavefunction amplitude using solution vector c

Parameters:

c (np.ndarray) – Wave function expansion coefficient
radial_factor (bool, optional) – whether to multiply result by \(r^{2}\), by default False

Returns:

List of spherical harmonic expansions per atom stored as \((N_{r} \times N_{lm})\) arrays where \(N_{r}\) is the number of radial points and \(N_{lm}\) the set of spherical harmonics used in the expansion.

Return type:

list[np.ndarray]

initialize()

Initialize the molecular grid

This is a relatively expensive procedure and thus the user can delay the initialization of the molecular grid until.

class pydft.AtomicGrid(at, nshells: int = 32, nangpts: int = 110, lmax: int = 8, fdpts: int = 7)

Bases: object

Represent the numerical integration grid attached to one atom.

An atomic grid combines a radial Gauss-Chebychev grid with a Lebedev angular grid. During a molecular calculation, each AtomicGrid stores the portion of the electron density assigned to its fuzzy Becke cell, projects that density onto real spherical harmonics, and solves the radial Poisson equations used to reconstruct the Hartree potential.

The class deliberately exposes several intermediate quantities, such as quadrature weights, spherical-harmonic coefficients, and Hartree-potential coefficients, because these are useful when studying how the molecular grid calculation is assembled from atom-centered pieces.

build_cubic_stencil_evaluation(r_eval)

Precompute cubic interpolation stencils for requested radii.

The atomic radial grid is stored from large to small radius. For interpolation, the grid is reversed into ascending order, and every molecular-grid radius is associated with a local four-point Lagrange stencil. The cached indices and weights are reused when the Hartree potential coefficients change during the SCF cycle.

Parameters:: r_eval (array_like) – Radii, measured relative to this atom, at which U_lm should be evaluated.

build_hartree_potential(): Construct the Hartree potential at the grid points including cubic interpolation to obtain the Hartree potential expansion coefficients at other positions

calculate_coulomb_energy()

Calculate the electron-electron repulsion for this atomic grid

This function is only for educational purposes, for any evaluation of the coulomb energy, extensive interpolation over all the atomic grids is necessary

calculate_coulomb_energy_interpolation()

Calculate the electron-electron repulsion for this atomic grid using cubic spline interpolation

This function is only for educational purposes, for any evaluation of the coulomb energy, extensive interpolation over all the atomic grids is necessary

calculate_interpolated_ulm(rr)

Calculate interpolated Hartree potential expansion coefficients for all points r in the array rr using CubicSpline interpolation

Note that CubicSpline expect that the r values are ascending; in the default grid, they are descending, which is why both rr and ulm are reversed.

calculate_interpolated_ulm_at_points(rpoints)

Interpolate Hartree-potential coefficients at arbitrary radii.

This method is a compatibility wrapper around calculate_interpolated_ulm(), which builds cubic splines directly for the requested radii. For repeated evaluations on the fixed molecular grid, calculate_interpolated_ulm_stencil() is faster.

calculate_interpolated_ulm_stencil()

Interpolate Hartree-potential coefficients at cached molecular radii.

build_cubic_stencil_evaluation() first precomputes, for each requested radius, the four neighboring radial grid indices and their cubic interpolation weights. This method applies those cached stencils to every spherical-harmonic channel of U_lm. It is equivalent in purpose to calculate_interpolated_ulm(), but avoids rebuilding a SciPy spline object for every SCF iteration.

Returns:: Interpolated Hartree-potential expansion coefficients with shape (N_lm, N_points).
Return type:: np.ndarray

calculate_nuclear_attraction(): Calculate the nuclear attraction given the electron density

count_electrons(): Sum over the electron density to obtain the number of electrons for this atomic cell

get_becke_weights(): Return the fuzzy cell weights (Becke weights)

get_bragg_slater_radius() → float

Get the Bragg-Slater radius

Returns:: Bragg-Slater radius
Return type:: float

get_charge(): Get the atomic charge

get_density(): Get the number of electrons for this atomic grid

get_dfa_exchange(): Get density functional approximation of the exchange energy for this atomic cell

get_dfa_kinetic(): Get density functional approximation of the kinetic energy for this atomic cell

get_dfa_nuclear_local()

Get the local density approximation to nuclear attraction for this cell.

Only the nucleus at the center of this atomic grid is included. The molecular-grid method adds the attraction from all nuclei when the full electron-nuclear energy is needed.

get_full_grid(): Get all the grid points as a list of positions (Nx3) matrix

get_gradient(): Get the gradient of the electron density for this atomic grid

get_gradient_squared(): Get the modulus squared of the gradient for this atomic grid

get_gridpoints(): Return list of grid points

get_lmax(): Get lmax value

get_local_hartree_potential()

Return the Hartree potential generated by this atomic grid alone.

This educational helper ignores interactions with the other atomic grids. It is valid only after __htpot has been generated by calculate_coulomb_energy() or calculate_coulomb_energy_interpolation().

get_nr_pts(): Get total number of sampling points for this atom

get_radial_grid(): Return radial grid

get_rho_lm(): Get electron density coefficients

get_ulm(): Get Hartree potential coefficients

get_weights(): Get the weights on the atomic grid for the quadrature

get_weights_angular_points(): Get the weights for the angular points only

get_ylm(): Get values for the spherical harmonics

get_ylm_atom(): Grab the spherical harmonic values at the angular points

perform_spherical_harmonic_expansion(f)

Calculate the expansion coefficients using a basis set composed of spherical harmonics given a atomic orbital

The expansion is performed over the function f which needs to be a matrix of Nk x Nj elements where Nk is the number of radial points and Nj is the number of angular points.

set_density(density): Set the density at the grid points

set_gradient(gradient): Set the electron-density gradient at the grid points.

set_molecular_weights(mweights): Set the molecular weights (Becke weights) at the grid points

Spherical harmonics

Real spherical harmonics used for density and Hartree-potential expansions.

class pydft.spherical_harmonics.SphericalHarmonicsCache

Cache for real spherical harmonics Y_lm evaluated on Lebedev angular grids.

Cached objects:: (l, nangpts) -> ndarray of shape (2*l + 1, nangpts)

Assumptions: - For a given nangpts, the Lebedev angular grid is deterministic. - Angular points are generated internally and never part of the cache key.

classmethod clear_cache(): Clear all cached spherical harmonics and angular grids.

classmethod freeze(): Prevent new cache entries from being created.

classmethod get_l(l, nangpts)

Return Y_lm for fixed l and all m in [-l, l].

Shape: (2*l + 1, nangpts)

classmethod get_ylm(lmax, nangpts)

Return all real spherical harmonics up to lmax on a Lebedev grid.

The rows are ordered by increasing l and, within each l, by m=-l, ..., l. The resulting array has shape ((lmax + 1)**2, nangpts).

classmethod precache_bulk(requests, max_workers=None)

Precompute spherical harmonics for multiple (lmax, nangpts) requests.

Parameters:

requests (iterable of (lmax, nangpts))
max_workers (int or None) – Passed to ThreadPoolExecutor

pydft.spherical_harmonics.job_compute_l(args): Helper function for the ThreadPoolExecutor

pydft.spherical_harmonics.real_sph_harm_l_legendre(l, theta, phi)

Real spherical harmonics matching your sph_harm_y-based definition, but computed via associated Legendre polynomials. This function seems to perform best.

Parameters:

l (int) – Angular momentum quantum number
theta (array_like) – Azimuthal angles
phi (array_like) – Polar angles

pydft.spherical_harmonics.real_sph_harm_l_scipy(l, theta, phi)

Evaluate real Spherical Harmonics for all m-values corresponding to a single l-value. This function is faster than looping over spherical_harmonic, yet remains slower than real_sph_harm_l_legendre, which is the preferred function for this task.

Parameters:

l (int) – Angular momentum quantum number
theta (array_like) – Azimuthal angles
phi (array_like) – Polar angles

pydft.spherical_harmonics.spherical_harmonic(l, m, theta, phi)

Calculate value of spherical harmonic function

l: angular quantum number m: magnetic quantum number theta: azimuthal angle in radians phi: polar angle in radians

pydft.spherical_harmonics.spherical_harmonic_cart(l, m, p)

Calculate the value of the spherical harmonic depending on the position p in Cartesian coordinates. This function assumes that the position p lies on the unit sphere

l: angular quantum number m: magnetic quantum number p: position three-vector on the unit sphere

Exchange-correlation functionals

Exchange-correlation energy densities and potentials for PyDFT.

class pydft.xcfunctionals.Functionals(functional='svwn5')

Class holding all exchange-correlation functionals

Note that all the functions ending in _deriv are the derivative of the energy per particle with respect to rho.

To calculate the potential, nu, one has to calculate deltaf/deltarho = df/drho * rho + f

calc_c(rho, grad=None, deriv_sigma=False)

Evaluate correlation energy density and potential terms.

Parameters and return values follow calc_x(), but for the correlation part of the selected exchange-correlation functional.

calc_x(rho, grad=None, deriv_sigma=False)

Evaluate exchange energy density and potential terms.

Parameters:

rho (array_like) – Electron density values on the numerical grid.
grad (array_like, optional) – GGA input \(\sigma = |\nabla\rho|^2\). Required for PBE and ignored for LDA.
deriv_sigma (bool, optional) – If True, also return the derivative with respect to sigma.

Returns:

(epsilon_x, v_x) for LDA-style use, or (epsilon_x, v_x, dF_x/dsigma) when deriv_sigma is true. Here epsilon_x is the exchange energy per particle and v_x = rho * d epsilon_x / d rho + epsilon_x is the local potential contribution.

Return type:

tuple

is_gga(): Returns whether the loaded XC functional is of GGA type

Internal numerical helpers

These helpers are included because they support educational inspection of the numerical algorithms. They are not usually needed for ordinary calculations.

Lebedev angular quadrature support for atom-centered integration grids.

class pydft.angulargrid.AngularGrid

Load and cache Lebedev quadrature rules by number of angular points.

PyDFT uses Lebedev grids to integrate functions over the angular part of an atom-centered spherical grid. This small wrapper translates the educational input used throughout PyDFT, nangpts, into the Lebedev order expected by pylebedev and stores the resulting angles, Cartesian unit-sphere points, and quadrature weights.

get_coefficients(numpoints: int): Get the coefficients for Lebedev quadrature on the unit sphere by number of points of that order

get_dataset_sizes(): Get the number of points on the unit sphere for each of the Lebedev orders

Numba-accelerated cubic stencil interpolation kernels.

pydft.stencil.stencil_interp(vals, x, weights, out)

Evaluate four-point cubic interpolation stencils for many channels.

Parameters:

vals (ndarray) – Source values with shape (radial_points, channels).
x (ndarray) – Base stencil indices for each evaluation point. Each index j uses source rows j-1 through j+2.
weights (ndarray) – Cubic interpolation weights with shape (points, 4).
out (ndarray) – Output array with shape (channels, points).