diff --git a/gbasis/base.py b/gbasis/base.py
index 6a3ab65c..3e866a61 100644
--- a/gbasis/base.py
+++ b/gbasis/base.py
@@ -1,4 +1,5 @@
 """Base class for arrays that depend on one or more contracted Gaussians."""
+
 import abc
 
 from gbasis.contractions import GeneralizedContractionShell
diff --git a/gbasis/base_four_symm.py b/gbasis/base_four_symm.py
index 28f375a6..6791b876 100644
--- a/gbasis/base_four_symm.py
+++ b/gbasis/base_four_symm.py
@@ -1,4 +1,5 @@
 """Base class for arrays that depend on four contracted Gaussians."""
+
 import abc
 import itertools as it
 
diff --git a/gbasis/base_one.py b/gbasis/base_one.py
index 5da708c3..ce2202b1 100644
--- a/gbasis/base_one.py
+++ b/gbasis/base_one.py
@@ -1,4 +1,5 @@
 """Base class for arrays that depend on one contracted Gaussian."""
+
 import abc
 
 from gbasis.base import BaseGaussianRelatedArray
diff --git a/gbasis/base_two_asymm.py b/gbasis/base_two_asymm.py
index 5d5c2e69..4ca3024e 100644
--- a/gbasis/base_two_asymm.py
+++ b/gbasis/base_two_asymm.py
@@ -1,4 +1,5 @@
 """Base class for arrays that depend on two contracted Gaussians."""
+
 import abc
 
 from gbasis.base import BaseGaussianRelatedArray
diff --git a/gbasis/evals/_deriv.py b/gbasis/evals/_deriv.py
index c357a82b..930fb901 100644
--- a/gbasis/evals/_deriv.py
+++ b/gbasis/evals/_deriv.py
@@ -1,4 +1,5 @@
 """Derivative of a Gaussian Contraction."""
+
 import numpy as np
 from scipy.special import comb, eval_hermite, perm
 
@@ -129,9 +130,7 @@ def _eval_deriv_contractions(coords, orders, center, angmom_comps, alphas, prim_
         # to evaluate multiple orders at the same time. Creating/finding a better function for
         # evaluating the hermite polynomial at different orders (in sequence) may be nice in the
         # future.
-        hermite = np.sum(
-            coeffs * eval_hermite(indices_herm, alphas**0.5 * nonzero_coords), axis=0
-        )
+        hermite = np.sum(coeffs * eval_hermite(indices_herm, alphas**0.5 * nonzero_coords), axis=0)
         hermite = np.prod(hermite, axis=1)
 
         # NOTE: `hermite` now has axis 0 for primitives, 1 for angular momentum vector, and axis 2
diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py
index 4d850975..e1c80b10 100644
--- a/gbasis/evals/electrostatic_potential.py
+++ b/gbasis/evals/electrostatic_potential.py
@@ -1,4 +1,5 @@
 """Module for computing electrostatic potential integrals."""
+
 from gbasis.integrals.point_charge import point_charge_integral
 import numpy as np
 
@@ -121,9 +122,11 @@ def electrostatic_potential(
     elif isinstance(coord_type, (list, tuple)):
         if (
             sum(
-                cont.num_sph * cont.num_seg_cont
-                if j == "spherical"
-                else cont.num_cart * cont.num_seg_cont
+                (
+                    cont.num_sph * cont.num_seg_cont
+                    if j == "spherical"
+                    else cont.num_cart * cont.num_seg_cont
+                )
                 for cont, j in zip(basis, coord_type)
             )
             != one_density_matrix.shape[0]
diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py
index 6a33ac51..24e6537e 100644
--- a/gbasis/evals/stress_tensor.py
+++ b/gbasis/evals/stress_tensor.py
@@ -1,4 +1,5 @@
 """Module for computing properties related to the stress tensor."""
+
 from gbasis.evals.density import (
     evaluate_density_laplacian,
     evaluate_deriv_density,
diff --git a/gbasis/integrals/__init__.py b/gbasis/integrals/__init__.py
index 62abae86..f583fb37 100644
--- a/gbasis/integrals/__init__.py
+++ b/gbasis/integrals/__init__.py
@@ -1 +1,4 @@
 """Collection of modules that compute different integrals of the contractions."""
+
+from .density import compute_intracule
+from .density import compute_extracule
diff --git a/gbasis/integrals/_diff_operator_int.py b/gbasis/integrals/_diff_operator_int.py
index d590d8bf..83f9a14f 100644
--- a/gbasis/integrals/_diff_operator_int.py
+++ b/gbasis/integrals/_diff_operator_int.py
@@ -1,4 +1,5 @@
 """Integrals over differential operator involving contracted Cartesian Gaussians."""
+
 from gbasis.integrals._moment_int import (
     _cleanup_intermediate_integrals,
     _compute_multipole_moment_integrals_intermediate,
diff --git a/gbasis/integrals/_moment_int.py b/gbasis/integrals/_moment_int.py
index 0fae05a6..a5b31a95 100644
--- a/gbasis/integrals/_moment_int.py
+++ b/gbasis/integrals/_moment_int.py
@@ -1,4 +1,5 @@
 """Multipole moment integrals involving Contracted Cartesian Gaussians."""
+
 import numpy as np
 
 
diff --git a/gbasis/integrals/_one_elec_int.py b/gbasis/integrals/_one_elec_int.py
index 688d7d9d..7a3506a4 100644
--- a/gbasis/integrals/_one_elec_int.py
+++ b/gbasis/integrals/_one_elec_int.py
@@ -1,4 +1,5 @@
 """One-electron integrals involving Contracted Cartesian Gaussians."""
+
 import numpy as np
 from gbasis.utils import factorial2
 
diff --git a/gbasis/integrals/_two_elec_int.py b/gbasis/integrals/_two_elec_int.py
index 6ff2a651..49d4ab14 100644
--- a/gbasis/integrals/_two_elec_int.py
+++ b/gbasis/integrals/_two_elec_int.py
@@ -1,4 +1,5 @@
 """Two-electron integrals involving Contracted Cartesian Gaussians."""
+
 import numpy as np
 from gbasis.utils import factorial2
 
diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py
new file mode 100644
index 00000000..b6a5514c
--- /dev/null
+++ b/gbasis/integrals/density.py
@@ -0,0 +1,137 @@
+"""
+Density API built on top of the arbitrary-order Gaussian overlap engine.
+
+This module provides high-level routines for constructing matrices of
+Gaussian overlap integrals between shells. These matrices serve as
+building blocks in density-related calculations such as intracule and
+extracule analysis.
+
+The overlap integrals are evaluated using the arbitrary-order Gaussian
+overlap engine implemented in ``arbitrary_order_overlap``.
+"""
+
+import numpy as np
+
+from .overlap_n import arbitrary_order_overlap
+
+
+def compute_intracule(shells):
+    """
+    Compute pairwise intracule overlap matrix between Gaussian shells.
+
+    This function constructs a matrix whose elements correspond to
+    two-center overlap integrals between Gaussian basis functions.
+    Each matrix element represents the overlap integral
+
+    .. math::
+
+        S_{ij} = \\int \\phi_i(\\mathbf{r}) \\phi_j(\\mathbf{r}) \\, d\\mathbf{r}
+
+    where :math:`\\phi_i` and :math:`\\phi_j` are Cartesian Gaussian basis
+    functions centered on different atoms or centers.
+
+    The overlap integrals are evaluated using the
+    :func:`arbitrary_order_overlap` engine, which computes Gaussian
+    overlap tensors of arbitrary order.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+        List of Gaussian shells defining the basis functions.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array of shape ``(n, n)`` containing pairwise Gaussian
+        overlap integrals between shells.
+
+    Notes
+    -----
+    The intracule matrix is constructed by evaluating pairwise
+    two-center overlap integrals between Gaussian shells. Each
+    integral is extracted from the sparse tensor returned by the
+    arbitrary-order overlap engine.
+
+    References
+    ----------
+    Helgaker, T., Jørgensen, P., Olsen, J.
+    *Molecular Electronic-Structure Theory*, Wiley (2000).
+
+    Szabo, A., Ostlund, N. S.
+    *Modern Quantum Chemistry*, Dover (1996).
+    """
+    n = len(shells)
+
+    result = np.zeros((n, n))
+
+    for i in range(n):
+        for j in range(n):
+
+            tensor = arbitrary_order_overlap([shells[i], shells[j]])
+
+            # Extract scalar overlap value from sparse tensor returned by
+            # the arbitrary-order overlap engine.
+            value = tensor.data[0] if tensor.nnz > 0 else 0.0
+
+            result[i, j] = value
+
+    return result
+
+
+def compute_extracule(shells):
+    """
+    Compute pairwise extracule overlap matrix between Gaussian shells.
+
+    This function constructs a matrix of pairwise overlap integrals
+    between Gaussian shells using the arbitrary-order overlap engine.
+
+    Each matrix element corresponds to the two-center Gaussian
+    overlap integral
+
+    .. math::
+
+        S_{ij} = \\int \\phi_i(\\mathbf{r}) \\phi_j(\\mathbf{r}) \\, d\\mathbf{r}
+
+    where :math:`\\phi_i` and :math:`\\phi_j` are Cartesian Gaussian
+    basis functions.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+        List of Gaussian shells defining the basis functions.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array of shape ``(n, n)`` containing pairwise Gaussian
+        overlap integrals between shells.
+
+    Notes
+    -----
+    This function provides an API wrapper around the
+    :func:`arbitrary_order_overlap` routine for density-related
+    calculations involving Gaussian basis functions.
+    This function currently constructs the same overlap matrix as
+    :func:`compute_intracule`, but is provided as a separate API
+    entry point for extracule-related density calculations.
+
+    References
+    ----------
+    Helgaker, T., Jørgensen, P., Olsen, J.
+    *Molecular Electronic-Structure Theory*, Wiley (2000).
+    """
+    n = len(shells)
+
+    result = np.zeros((n, n))
+
+    for i in range(n):
+        for j in range(n):
+            # evaluate two-center Gaussian overlap tensor between shells i and j
+            tensor = arbitrary_order_overlap([shells[i], shells[j]])
+
+            # extract scalar overlap value from sparse tensor
+            value = tensor.data[0] if tensor.nnz > 0 else 0.0
+
+            result[i, j] = value
+
+    return result
diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py
index a6df60ff..1aa423ce 100644
--- a/gbasis/integrals/libcint.py
+++ b/gbasis/integrals/libcint.py
@@ -21,7 +21,6 @@
 
 from gbasis.utils import factorial2
 
-
 __all__ = [
     "LIBCINT",
     "CBasis",
@@ -193,6 +192,7 @@ def from_param(cls, obj):
 
 class PairData(Structure):
     r"""``libcint`` ``PairData`` class."""
+
     _fields_ = [
         ("rij", c_double * 3),
         ("eij", c_double),
@@ -202,6 +202,7 @@ class PairData(Structure):
 
 class CINTOpt(Structure):
     r"""``libcint`` ``CINTOpt`` class."""
+
     _fields_ = [
         ("index_xyz_array", POINTER(POINTER(c_int))),
         ("non0ctr", POINTER(POINTER(c_int))),
diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py
index 00e2241e..ec7c4d6b 100644
--- a/gbasis/integrals/nuclear_electron_attraction.py
+++ b/gbasis/integrals/nuclear_electron_attraction.py
@@ -1,4 +1,5 @@
 """Module for computing the nuclear electron attraction."""
+
 from gbasis.integrals.point_charge import point_charge_integral
 import numpy as np
 
diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py
new file mode 100644
index 00000000..62322631
--- /dev/null
+++ b/gbasis/integrals/overlap_n.py
@@ -0,0 +1,520 @@
+"""
+Arbitrary-order N-center Gaussian overlap integrals.
+
+This module implements algorithms for computing overlap integrals
+between Cartesian Gaussian basis functions located on different
+centers. The implementation supports overlap integrals involving
+an arbitrary number of Gaussian shells.
+
+The evaluation follows the Gaussian product theorem and Hermite
+Gaussian expansion techniques commonly used in quantum chemistry
+integral engines.
+
+References
+----------
+Helgaker, T., Jørgensen, P., Olsen, J.
+Molecular Electronic-Structure Theory, Wiley (2000).
+
+Szabo, A., Ostlund, N. S.
+Modern Quantum Chemistry, Dover (1996).
+
+Obara, S., Saika, A. (1986)
+Efficient recursive computation of molecular integrals over
+Cartesian Gaussian functions.
+Journal of Chemical Physics.
+"""
+
+import numpy as np
+from scipy.sparse import coo_matrix
+from itertools import product
+
+
+class PrimitiveNEngine:
+
+    # Gaussian collapse
+    @staticmethod
+    def collapse_gaussians(alphas, centers):
+        """
+        Collapse multiple primitive Gaussian functions into a single
+        effective Gaussian using the Gaussian product theorem.
+
+        Given primitive Gaussian functions
+
+            exp(-α_i |r - A_i|²)
+
+        their product can be written as
+
+            exp(-α_tot |r - P|²) * exp(term2 - term1)
+
+        where
+
+            α_tot = Σ α_i
+            P = (Σ α_i A_i) / α_tot
+
+        This transformation simplifies the evaluation of Gaussian
+        integrals by reducing multiple Gaussians into a single
+        Gaussian centered at P.
+
+        Parameters
+        ----------
+        alphas : array_like
+            Exponents of the primitive Gaussian functions.
+
+        centers : array_like
+            Cartesian coordinates of the Gaussian centers.
+
+        Returns
+        -------
+        alpha_tot : float
+            Total Gaussian exponent.
+
+        P : ndarray
+            Product center of the collapsed Gaussian.
+
+        prefactor : float
+            Exponential prefactor arising from Gaussian collapse.
+            Raises:
+            ------
+            ValueError:
+            If the total Gaussian exponent is non-positive.
+        """
+        alphas = np.asarray(alphas, dtype=np.float64)
+        centers = np.asarray(centers, dtype=np.float64)
+
+        alpha_tot = np.sum(alphas)
+
+        if alpha_tot <= 0.0:
+            raise ValueError("Total Gaussian exponent must be positive.")
+
+        P = np.sum(alphas[:, None] * centers, axis=0) / alpha_tot
+
+        term1 = np.sum(alphas * np.sum(centers**2, axis=1))
+        term2 = alpha_tot * np.dot(P, P)
+
+        exponent = term2 - term1
+        prefactor = np.exp(exponent)
+
+        return alpha_tot, P, prefactor
+
+    # Pure binomial Hermite shift
+
+    @staticmethod
+    def hermite_coefficients(l, PA):
+        """
+        Compute Hermite translation coefficients.
+
+        This function expands a polynomial centered at A in terms of
+        powers centered at the Gaussian product center P:
+
+        .. math::
+
+            (x - A)^l = \\sum_t E_t (x - P)^t
+
+        These coefficients are used in the Hermite Gaussian formalism
+        for evaluating molecular integrals.
+
+        Parameters
+        ----------
+        l : int
+            Angular momentum order.
+
+        PA : float
+            Distance between the product center P and the original
+            Gaussian center A.
+
+        Returns
+        -------
+        numpy.ndarray
+            Array containing the Hermite expansion coefficients.
+        """
+        E = np.zeros(l + 1, dtype=np.float64)
+        E[0] = 1.0
+
+        for i in range(l):
+            E_new = np.zeros(l + 1, dtype=np.float64)
+
+            for t in range(i + 1):
+                E_new[t] += PA * E[t]
+                E_new[t + 1] += E[t]
+
+            E = E_new
+
+        return E
+
+    # Gaussian moments
+
+    @staticmethod
+    def gaussian_moments(alpha, max_order):
+        """
+        Compute Gaussian moment integrals.
+
+        Evaluates integrals of the form
+
+        .. math::
+
+            \\int_{-\\infty}^{\\infty} (x-P)^k e^{-\\alpha (x-P)^2} dx
+
+        Only even moments are non-zero. Higher moments are computed
+        recursively from the zeroth moment.
+
+        Parameters
+        ----------
+        alpha : float
+            Gaussian exponent.
+
+        max_order : int
+            Maximum moment order to compute.
+
+        Returns
+        -------
+        numpy.ndarray
+            Array containing Gaussian moments up to ``max_order``.
+        """
+        moments = np.zeros(max_order + 1, dtype=np.float64)
+
+        # zeroth moment
+        moments[0] = np.sqrt(np.pi / alpha)
+
+        # only even moments survive
+        for k in range(0, max_order - 1, 2):
+            moments[k + 2] = (k + 1) / (2.0 * alpha) * moments[k]
+
+        return moments
+
+    # Full primitive N-center overlap
+    @staticmethod
+    def primitive_overlap(alphas, centers, angmoms):
+        """
+        Compute primitive N-center Gaussian overlap integral.
+
+        This function evaluates the overlap integral between N primitive
+        Gaussian basis functions using the Hermite Gaussian formalism.
+
+        The primitive overlap integral is
+
+        .. math::
+
+            S = \\int \\prod_i \\phi_i(\\mathbf{r}) \\, d\\mathbf{r}
+
+        where each primitive Gaussian basis function has the form
+
+        .. math::
+
+            \\phi(\\mathbf{r}) =
+            (x-A_x)^l (y-A_y)^m (z-A_z)^n
+            e^{-\\alpha |\\mathbf{r}-A|^2}
+
+        Parameters
+        ----------
+        alphas : list
+            Gaussian exponents.
+
+        centers : list
+            Cartesian coordinates of Gaussian centers.
+
+        angmoms : list
+            Angular momentum tuples ``(l, m, n)``.
+
+        Returns
+        -------
+        float
+            Value of the primitive N-center overlap integral.
+
+        Notes
+        -----
+        The overlap integral factorizes into independent Cartesian
+        components (x, y, z). Each component is evaluated using
+        Gaussian moment integrals obtained from the Hermite expansion.
+        """
+        alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians(alphas, centers)
+
+        result = prefactor
+
+        # factorize into x, y, z
+        # Evaluate overlap as product of x, y, z integrals
+
+        for axis in range(3):
+
+            # build total polynomial via convolution
+            E_total = np.array([1.0], dtype=np.float64)
+
+            for i in range(len(alphas)):
+                l = angmoms[i][axis]
+                PA = P[axis] - centers[i][axis]
+
+                E = PrimitiveNEngine.hermite_coefficients(l, PA)
+
+                E_total = np.convolve(E_total, E)
+
+            moments = PrimitiveNEngine.gaussian_moments(alpha_tot, len(E_total) - 1)
+
+            axis_integral = np.dot(E_total, moments[: len(E_total)])
+
+            result *= axis_integral
+
+        return result
+
+
+# Screening Function
+def is_n_shell_overlap_screened(shells, tol=1e-12):
+    """
+    Determine whether an N-center overlap integral can be skipped
+    using a conservative exponential bound.
+
+    This screening estimates an upper bound for the magnitude of
+    the contracted overlap integral using Gaussian decay between
+    shell centers. If the bound is below the specified tolerance,
+    the integral is considered negligible.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+        Shells involved in the overlap integral.
+
+    tol : float
+        Screening tolerance.
+
+    Returns
+    -------
+    bool
+        True if the integral is negligible and can be skipped.
+    """
+
+    alpha_mins = [np.min(shell.exps) for shell in shells]
+    centers = [shell.coord for shell in shells]
+
+    alpha_tot = sum(alpha_mins)
+
+    if alpha_tot <= 0.0:
+        return True
+
+    # Exponential decay from Gaussian collapse
+    decay_sum = 0.0
+    N = len(shells)
+
+    for i in range(N):
+        for j in range(i + 1, N):
+            Rij = centers[i] - centers[j]
+            Rij2 = np.dot(Rij, Rij)
+            decay_sum += alpha_mins[i] * alpha_mins[j] * Rij2
+
+    D = decay_sum / alpha_tot
+
+    # Contraction-level magnitude bound
+    coeff_bound = 1.0
+    norm_bound = 1.0
+
+    for shell in shells:
+        coeff_bound *= np.max(np.abs(shell.coeffs))
+        norm_bound *= np.max(np.abs(shell.norm_prim_cart))
+
+    volume_bound = (np.pi / alpha_tot) ** 1.5
+
+    bound = coeff_bound * norm_bound * volume_bound * np.exp(-D)
+
+    return bound < tol
+
+
+def contracted_n_overlap(shells):
+    """
+    Compute contracted N-center overlap for a list of
+    GeneralizedContractionShell objects.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+
+    Returns
+    -------
+    np.ndarray
+        N-dimensional array over segmented contractions
+        and Cartesian angular components.
+
+        Shape:
+        (M1, L1, M2, L2, ..., MN, LN)
+    """
+
+    N = len(shells)
+
+    # Build shape for final tensor
+    shape = []
+    for shell in shells:
+        shape.append(shell.num_seg_cont)
+        shape.append(shell.num_cart)
+
+    result = np.zeros(shape, dtype=np.float64)
+
+    # Primitive exponent index ranges
+    prim_ranges = [range(len(shell.exps)) for shell in shells]
+
+    # Segmented contraction indices
+    seg_ranges = [range(shell.num_seg_cont) for shell in shells]
+
+    # Cartesian angular component indices
+    cart_ranges = [range(shell.num_cart) for shell in shells]
+
+    for seg_indices in product(*seg_ranges):
+        for cart_indices in product(*cart_ranges):
+
+            total_value = 0.0
+            for prim_indices in product(*prim_ranges):
+
+                alphas = []
+                centers = []
+                angmoms = []
+                coeff_prod = 1.0
+                norm_prod = 1.0
+
+                for i, shell in enumerate(shells):
+
+                    p = prim_indices[i]
+                    m = seg_indices[i]
+                    c = cart_indices[i]
+
+                    alpha = shell.exps[p]
+                    coeff = shell.coeffs[
+                        p, m
+                    ]  # contraction coefficient for primitive p in segment m
+                    norm = shell.norm_prim_cart[
+                        c, p
+                    ]  # normalization factor for primitive p and Cartesian component c
+                    angmom = tuple(
+                        shell.angmom_components_cart[c]
+                    )  # Cartesian angular momentum (lx, ly, lz)
+
+                    alphas.append(alpha)
+                    centers.append(shell.coord)
+                    angmoms.append(angmom)
+
+                    coeff_prod *= coeff
+                    norm_prod *= norm
+
+                prim_val = PrimitiveNEngine.primitive_overlap(alphas, centers, angmoms)
+
+                total_value += coeff_prod * norm_prod * prim_val
+
+            index = []
+            for i in range(N):
+                index.append(seg_indices[i])
+                index.append(cart_indices[i])
+
+            result[tuple(index)] = total_value
+
+    return result
+
+
+def build_n_overlap_tensor(shells, tol=1e-12):
+    """
+    Construct sparse N-center overlap tensor.
+
+    This function builds the full N-center overlap tensor over
+    Gaussian basis functions. Each tensor element corresponds
+    to the overlap integral between a combination of contracted
+    Gaussian shells.
+
+    To reduce computational cost, shell-level screening is
+    applied before evaluating the contracted integrals.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+        List of Gaussian shells participating in the overlap
+        integral.
+
+    tol : float
+        Screening tolerance used to discard negligible
+        contributions.
+
+    Returns
+    -------
+    scipy.sparse.coo_matrix
+        Sparse tensor containing the N-center overlap integrals.
+        The tensor is flattened so that each N-index integral
+        maps to a single row index.
+    """
+
+    # Total AO dimension
+    shell_sizes = [shell.num_seg_cont * shell.num_cart for shell in shells]
+
+    total_ao = sum(shell_sizes)
+    N = len(shells)
+
+    data = []
+    rows = []
+
+    # compute AO offsets per shell
+    offsets = []
+    acc = 0
+    for size in shell_sizes:
+        offsets.append(acc)
+        acc += size
+
+    for shell_indices in product(range(len(shells)), repeat=N):
+
+        shell_tuple = [shells[i] for i in shell_indices]
+
+        # Screening
+        if is_n_shell_overlap_screened(shell_tuple, tol=tol):
+            continue
+
+        block = contracted_n_overlap(shell_tuple)
+
+        block_flat = block.reshape(-1)
+
+        local_sizes = [shells[i].num_seg_cont * shells[i].num_cart for i in shell_indices]
+
+        local_offsets = [offsets[i] for i in shell_indices]
+
+        for local_idx, value in enumerate(block_flat):
+
+            if abs(value) < tol:
+                continue
+
+            # convert local multi-index to global index
+            multi = []
+            tmp = local_idx
+
+            for size in reversed(local_sizes):
+                multi.append(tmp % size)
+                tmp //= size
+
+            multi = list(reversed(multi))
+
+            global_index = 0
+            for k in range(N):
+                global_index = global_index * total_ao + local_offsets[k] + multi[k]
+
+            rows.append(global_index)
+            data.append(value)
+
+    shape = (total_ao**N, 1)
+
+    return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape)
+
+
+# Public API function (Week 3 deliverable)
+def arbitrary_order_overlap(shells, tol=1e-12):
+    """
+    Compute arbitrary-order Gaussian overlap tensor.
+
+    This is the main public API for evaluating N-center overlap
+    integrals between Gaussian basis functions.
+
+    Internally this function constructs the overlap tensor using
+    ``build_n_overlap_tensor``.
+
+    Parameters
+    ----------
+    shells : list[GeneralizedContractionShell]
+        List of Gaussian shells defining the basis functions.
+
+    tol : float
+        Screening tolerance.
+
+    Returns
+    -------
+    scipy.sparse.coo_matrix
+        Sparse representation of the N-center overlap tensor.
+    """
+
+    return build_n_overlap_tensor(shells, tol)
diff --git a/gbasis/screening.py b/gbasis/screening.py
index 162aec54..6ae61544 100644
--- a/gbasis/screening.py
+++ b/gbasis/screening.py
@@ -158,6 +158,7 @@ def compute_primitive_cutoff_radius(c, alpha, angm, deriv_order, tol_screen):
     lambert_input_value = -2 * alpha * (eff_tol_screen / (c * n)) ** (2 / eff_angm) / eff_angm
     return np.sqrt(-(eff_angm / (2 * alpha)) * lambertw(lambert_input_value, k=-1).real)
 
+
 # TODO: Fix this, it fails for some reason, it is needed for screening of 1rdms
 def compute_contraction_upper_bond(contractions, deriv_order):
     r"""Compute an upper bound for a contraction or its derivatives for any point.
diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py
index f915e26c..a4445fc8 100644
--- a/tests/test_angular_momentum.py
+++ b/tests/test_angular_momentum.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.angular_momentum."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import (
     _compute_differential_operator_integrals_intermediate,
diff --git a/tests/test_base.py b/tests/test_base.py
index cb4d5497..47a1658e 100644
--- a/tests/test_base.py
+++ b/tests/test_base.py
@@ -1,4 +1,5 @@
 """Test gbasis.base.BaseGuassianRelatedArray."""
+
 from gbasis.base import BaseGaussianRelatedArray
 from gbasis.contractions import GeneralizedContractionShell
 import numpy as np
diff --git a/tests/test_base_four_symm.py b/tests/test_base_four_symm.py
index 374c6155..4aa0b371 100644
--- a/tests/test_base_four_symm.py
+++ b/tests/test_base_four_symm.py
@@ -1,4 +1,5 @@
 """Test gbasis.base_four_symmetric."""
+
 from gbasis.base_four_symm import BaseFourIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.spherical import generate_transformation
diff --git a/tests/test_base_one.py b/tests/test_base_one.py
index dd2975f1..b1cae29a 100644
--- a/tests/test_base_one.py
+++ b/tests/test_base_one.py
@@ -1,4 +1,5 @@
 """Test gbasis.base_one."""
+
 from gbasis.base_one import BaseOneIndex
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.spherical import generate_transformation
diff --git a/tests/test_base_two_asymm.py b/tests/test_base_two_asymm.py
index 4b609c00..7bdb3cac 100644
--- a/tests/test_base_two_asymm.py
+++ b/tests/test_base_two_asymm.py
@@ -1,4 +1,5 @@
 """Test gbasis.base_two_asymm."""
+
 from gbasis.base_two_asymm import BaseTwoIndexAsymmetric
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.spherical import generate_transformation
diff --git a/tests/test_base_two_symm.py b/tests/test_base_two_symm.py
index ac6be7a3..9a30bf91 100644
--- a/tests/test_base_two_symm.py
+++ b/tests/test_base_two_symm.py
@@ -1,4 +1,5 @@
 """Test gbasis.base_two_symm."""
+
 from gbasis.base_two_asymm import BaseTwoIndexAsymmetric
 from gbasis.base_two_symm import BaseTwoIndexSymmetric
 from gbasis.contractions import GeneralizedContractionShell
diff --git a/tests/test_contractions.py b/tests/test_contractions.py
index 0d446b06..93851862 100644
--- a/tests/test_contractions.py
+++ b/tests/test_contractions.py
@@ -1,4 +1,5 @@
 """Test gbasis.contractions."""
+
 from gbasis.contractions import GeneralizedContractionShell
 import numpy as np
 import pytest
diff --git a/tests/test_deriv.py b/tests/test_deriv.py
index 0d725fd7..cc30647c 100644
--- a/tests/test_deriv.py
+++ b/tests/test_deriv.py
@@ -1,4 +1,5 @@
 """Test gbasis.evals._deriv."""
+
 import itertools as it
 
 from gbasis.evals._deriv import _eval_deriv_contractions
diff --git a/tests/test_diff_operator_int.py b/tests/test_diff_operator_int.py
index c5a40c56..cb08f70b 100644
--- a/tests/test_diff_operator_int.py
+++ b/tests/test_diff_operator_int.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals._diff_operator_int."""
+
 import itertools as it
 
 from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals
diff --git a/tests/test_electrostatic_potential.py b/tests/test_electrostatic_potential.py
index 520349c8..089b76d5 100644
--- a/tests/test_electrostatic_potential.py
+++ b/tests/test_electrostatic_potential.py
@@ -1,4 +1,5 @@
 """Tests for gbasis.evals.electrostatic_potential."""
+
 from gbasis.evals.electrostatic_potential import electrostatic_potential
 from gbasis.parsers import make_contractions, parse_nwchem
 import numpy as np
diff --git a/tests/test_eval_deriv.py b/tests/test_eval_deriv.py
index 4f318170..22c10f35 100644
--- a/tests/test_eval_deriv.py
+++ b/tests/test_eval_deriv.py
@@ -1,4 +1,5 @@
 """Test gbasis.evals.evaluate_deriv."""
+
 import itertools as it
 
 from gbasis.contractions import GeneralizedContractionShell
diff --git a/tests/test_extracule.py b/tests/test_extracule.py
new file mode 100644
index 00000000..07f7d32a
--- /dev/null
+++ b/tests/test_extracule.py
@@ -0,0 +1,32 @@
+import numpy as np
+
+from gbasis.integrals.density import compute_extracule
+
+
+class DummyTensor:
+    def __init__(self):
+        self.data = np.array([1.0])
+        self.nnz = 1
+
+
+class DummyShell:
+    pass
+
+
+import gbasis.integrals.density as density
+
+
+def dummy_overlap(shells):
+    return DummyTensor()
+
+
+density.arbitrary_order_overlap = dummy_overlap
+
+
+def test_extracule_shape():
+
+    shells = [DummyShell(), DummyShell()]
+
+    result = compute_extracule(shells)
+
+    assert result.shape == (2, 2)
diff --git a/tests/test_extracule_real.py b/tests/test_extracule_real.py
new file mode 100644
index 00000000..75505338
--- /dev/null
+++ b/tests/test_extracule_real.py
@@ -0,0 +1,52 @@
+import numpy as np
+
+from gbasis.parsers import parse_gbs
+from gbasis.contractions import GeneralizedContractionShell
+from gbasis.integrals.density import compute_extracule
+
+
+def convert_to_shells(parsed_shells):
+    """
+    Convert parsed GBS tuples into GeneralizedContractionShell objects.
+    """
+
+    origin = np.array([0.0, 0.0, 0.0], dtype=float)
+
+    shells = []
+
+    for angmom, exps, coeffs in parsed_shells:
+
+        shell = GeneralizedContractionShell(
+            angmom,  # angular momentum
+            origin,  # center coordinate
+            np.array(coeffs, dtype=float),  # contraction coefficients
+            np.array(exps, dtype=float),  # exponents
+            "cartesian",  # coordinate type
+        )
+
+        shells.append(shell)
+
+    return shells
+
+
+def test_extracule_real_basis():
+    """
+    Validate compute_extracule using real GBasis basis data.
+    """
+
+    basis_dict = parse_gbs("tests/data_631g.gbs")
+
+    parsed_shells = basis_dict["H"]
+
+    shells = convert_to_shells(parsed_shells)
+
+    result = compute_extracule(shells)
+
+    # shape validation
+    assert result.shape == (len(shells), len(shells))
+
+    # symmetry validation
+    assert np.allclose(result, result.T)
+
+    # diagonal must be non-negative
+    assert np.all(result.diagonal() >= 0)
diff --git a/tests/test_integrals_density.py b/tests/test_integrals_density.py
new file mode 100644
index 00000000..91a18505
--- /dev/null
+++ b/tests/test_integrals_density.py
@@ -0,0 +1,37 @@
+import numpy as np
+
+from gbasis.integrals.density import (
+    compute_intracule,
+    compute_extracule,
+)
+
+
+class DummyShell:
+
+    def __init__(self):
+
+        self.exps = np.array([1.0])
+        self.coeffs = np.array([[1.0]])
+        self.norm_prim_cart = np.array([[1.0]])
+        self.angmom_components_cart = [(0, 0, 0)]
+        self.coord = np.array([0.0, 0.0, 0.0])
+        self.num_seg_cont = 1
+        self.num_cart = 1
+
+
+def test_compute_intracule():
+
+    shells = [DummyShell(), DummyShell()]
+
+    result = compute_intracule(shells)
+
+    assert result is not None
+
+
+def test_compute_extracule():
+
+    shells = [DummyShell(), DummyShell()]
+
+    result = compute_extracule(shells)
+
+    assert result is not None
diff --git a/tests/test_intracule.py b/tests/test_intracule.py
new file mode 100644
index 00000000..3fd89698
--- /dev/null
+++ b/tests/test_intracule.py
@@ -0,0 +1,51 @@
+import numpy as np
+
+from gbasis.integrals.density import compute_intracule
+
+
+class DummyTensor:
+    def __init__(self):
+        self.data = np.array([1.0])
+        self.nnz = 1
+
+
+class DummyShell:
+    pass
+
+
+# Monkeypatch arbitrary_order_overlap
+import gbasis.integrals.density as density
+
+
+def dummy_overlap(shells):
+    return DummyTensor()
+
+
+density.arbitrary_order_overlap = dummy_overlap
+
+
+def test_intracule_shape():
+
+    shells = [DummyShell(), DummyShell(), DummyShell()]
+
+    result = compute_intracule(shells)
+
+    assert result.shape == (3, 3)
+
+
+def test_intracule_empty():
+
+    shells = []
+
+    result = compute_intracule(shells)
+
+    assert result.shape == (0, 0)
+
+
+def test_intracule_symmetry():
+
+    shells = [DummyShell(), DummyShell()]
+
+    result = compute_intracule(shells)
+
+    assert np.allclose(result, result.T)
diff --git a/tests/test_intracule_real.py b/tests/test_intracule_real.py
new file mode 100644
index 00000000..fb7a89cc
--- /dev/null
+++ b/tests/test_intracule_real.py
@@ -0,0 +1,42 @@
+import numpy as np
+from gbasis.parsers import parse_gbs
+from gbasis.contractions import GeneralizedContractionShell
+from gbasis.integrals.density import compute_intracule
+
+
+def convert_to_shells(parsed_shells):
+
+    origin = np.array([0.0, 0.0, 0.0], dtype=float)
+
+    shells = []
+
+    for angmom, exps, coeffs in parsed_shells:
+
+        shell = GeneralizedContractionShell(
+            angmom,  # 1
+            origin,  # 2
+            np.array(coeffs, dtype=float),  # 3
+            np.array(exps, dtype=float),  # 4
+            "cartesian",  # 5
+        )
+
+        shells.append(shell)
+
+    return shells
+
+
+def test_intracule_real_basis():
+
+    basis_dict = parse_gbs("tests/data_631g.gbs")
+
+    parsed_shells = basis_dict["H"]
+
+    shells = convert_to_shells(parsed_shells)
+
+    result = compute_intracule(shells)
+
+    assert result.shape == (len(shells), len(shells))
+
+    assert np.allclose(result, result.T)
+
+    assert np.all(result.diagonal() >= 0)
diff --git a/tests/test_kinetic_energy.py b/tests/test_kinetic_energy.py
index 949b975e..b58c18eb 100644
--- a/tests/test_kinetic_energy.py
+++ b/tests/test_kinetic_energy.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.kinetic_energy."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals
 from gbasis.integrals.kinetic_energy import kinetic_energy_integral, KineticEnergyIntegral
@@ -250,10 +251,11 @@ def test_kinetic_energy_integral_horton_anorcc_bec():
     ke_int_values = kinetic_energy_integral(basis, screen_basis=False)
     assert np.allclose(ke_int_values, horton_kinetic_energy_integral)
 
+
 @pytest.mark.parametrize("precision", [1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8])
 def test_kinetic_screening_accuracy(precision):
     """Test kinetic energy screening."""
-    
+
     basis_dict = parse_gbs(find_datafile("data_631g.gbs"))
     atsymbols = ["H", "C", "Kr"]
     atcoords = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
diff --git a/tests/test_libcint.py b/tests/test_libcint.py
index 1d1deddb..d279a897 100644
--- a/tests/test_libcint.py
+++ b/tests/test_libcint.py
@@ -29,7 +29,6 @@
 
 from utils import find_datafile
 
-
 TEST_BASIS_SETS = [
     pytest.param("data_sto6g.nwchem", id="STO-6G"),
     pytest.param("data_631g.nwchem", id="6-31G"),
@@ -66,6 +65,7 @@
     pytest.param("moment", id="Moment"),
 ]
 
+
 @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows")
 @pytest.mark.skipif(
     len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0,
@@ -175,11 +175,10 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral):
     npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol)
 
 
-
 TEST_SYSTEMS_IODATA = [
     pytest.param("h2o_hf_ccpv5z_cart.fchk", ["O", "H", "H"], "Cartesian", id="h2o_cart"),
     pytest.param("h2o_hf_ccpv5z_sph.fchk", ["O", "H", "H"], "Spherical", id="h2o_sph"),
-    ]
+]
 
 TEST_COORD_TRANSFORM = [
     pytest.param(False, id="no-transform"),
@@ -192,10 +191,14 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral):
     pytest.param("nuclear_attraction", id="NuclearAttraction"),
     pytest.param("momentum", id="Momentum"),
     pytest.param("angular_momentum", id="AngularMomentum"),
-    pytest.param("electron_repulsion", marks=pytest.mark.skip(reason='TOO SLOW'), id="ElectronRepulsion"),
+    pytest.param(
+        "electron_repulsion", marks=pytest.mark.skip(reason="TOO SLOW"), id="ElectronRepulsion"
+    ),
     pytest.param("point_charge", id="PointCharge"),
     pytest.param("moment", id="Moment"),
 ]
+
+
 @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows")
 @pytest.mark.skipif(
     len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0,
@@ -211,8 +214,8 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform):
 
     atol, rtol = 1e-6, 1e-6
 
-    mol=load_one(find_datafile(fname))
-    py_basis=from_iodata(mol)
+    mol = load_one(find_datafile(fname))
+    py_basis = from_iodata(mol)
 
     lc_basis = CBasis(py_basis, elements, mol.atcoords, coord_type=coord_type)
 
@@ -230,7 +233,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform):
 
     elif integral == "kinetic_energy":
         if transform:
-            py_int = kinetic_energy_integral(py_basis, transform=mol.mo.coeffs.T, screen_basis=False)
+            py_int = kinetic_energy_integral(
+                py_basis, transform=mol.mo.coeffs.T, screen_basis=False
+            )
             npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn))
             lc_int = lc_basis.kinetic_energy_integral(transform=mol.mo.coeffs.T)
             npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn))
@@ -242,8 +247,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform):
 
     elif integral == "nuclear_attraction":
         if transform:
-            py_int = nuclear_electron_attraction_integral(py_basis, mol.atcoords,
-                                                          mol.atnums, transform=mol.mo.coeffs.T)
+            py_int = nuclear_electron_attraction_integral(
+                py_basis, mol.atcoords, mol.atnums, transform=mol.mo.coeffs.T
+            )
             npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn))
             lc_int = lc_basis.nuclear_attraction_integral(transform=mol.mo.coeffs.T)
             npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn))
@@ -298,11 +304,13 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform):
         charges = np.asarray([1.0, 0.666, -3.1415926])
         if transform:
             for i in range(1, len(charges) + 1):
-                py_int = point_charge_integral(py_basis, charge_coords[:i],
-                                               charges[:i], transform=mol.mo.coeffs.T)
+                py_int = point_charge_integral(
+                    py_basis, charge_coords[:i], charges[:i], transform=mol.mo.coeffs.T
+                )
                 npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i))
-                lc_int = lc_basis.point_charge_integral(charge_coords[:i],
-                                                        charges[:i], transform=mol.mo.coeffs.T)
+                lc_int = lc_basis.point_charge_integral(
+                    charge_coords[:i], charges[:i], transform=mol.mo.coeffs.T
+                )
                 npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i))
 
         else:
@@ -329,7 +337,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform):
             ]
         )
         if transform:
-            py_int = moment_integral(py_basis, origin, orders, transform=mol.mo.coeffs.T, screen_basis=False)
+            py_int = moment_integral(
+                py_basis, origin, orders, transform=mol.mo.coeffs.T, screen_basis=False
+            )
             npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders)))
             lc_int = lc_basis.moment_integral(orders, origin=origin, transform=mol.mo.coeffs.T)
             npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders)))
diff --git a/tests/test_moment.py b/tests/test_moment.py
index 6d920fab..b68b6e21 100644
--- a/tests/test_moment.py
+++ b/tests/test_moment.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.moment."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._moment_int import _compute_multipole_moment_integrals
 from gbasis.integrals.moment import Moment, moment_integral
diff --git a/tests/test_moment_int.py b/tests/test_moment_int.py
index 13b1348b..01ac4531 100644
--- a/tests/test_moment_int.py
+++ b/tests/test_moment_int.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals._moment_int."""
+
 import itertools as it
 
 from gbasis.integrals._moment_int import _compute_multipole_moment_integrals
diff --git a/tests/test_momentum.py b/tests/test_momentum.py
index 718e6580..979588ad 100644
--- a/tests/test_momentum.py
+++ b/tests/test_momentum.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.momentum."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals
 from gbasis.integrals.momentum import momentum_integral, MomentumIntegral
diff --git a/tests/test_n.py b/tests/test_n.py
new file mode 100644
index 00000000..b9af0754
--- /dev/null
+++ b/tests/test_n.py
@@ -0,0 +1,50 @@
+import numpy as np
+import pytest
+
+from gbasis.parsers import parse_gbs, make_contractions
+from gbasis.integrals.overlap_n import contracted_n_overlap
+
+# =====================================================
+# Helper: Build hydrogen basis
+# =====================================================
+
+
+def build_h2_basis():
+
+    atoms = ["H", "H"]
+
+    coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+
+    basis_dict = parse_gbs("tests/data_631g.gbs")
+
+    basis = make_contractions(basis_dict, atoms, coords, coord_types="c")
+
+    return basis
+
+
+# =====================================================
+# Test arbitrary orders N=1 to N=6
+# =====================================================
+
+
+@pytest.mark.parametrize("N", [1, 2, 3, 4, 5, 6])
+def test_arbitrary_n_orders(N):
+
+    basis = build_h2_basis()
+
+    shells = []
+
+    for i in range(N):
+        shells.append(basis[i % len(basis)])
+
+    result = contracted_n_overlap(shells)
+
+    print(f"N={N}, max abs value:", np.max(np.abs(result)))
+
+    assert result is not None
+
+    assert not np.isnan(result).any()
+
+    assert not np.isinf(result).any()
+
+    assert np.all(np.isfinite(result))
diff --git a/tests/test_n_engine_stress.py b/tests/test_n_engine_stress.py
new file mode 100644
index 00000000..ecce90cd
--- /dev/null
+++ b/tests/test_n_engine_stress.py
@@ -0,0 +1,111 @@
+import numpy as np
+import pytest
+
+from gbasis.parsers import parse_gbs, make_contractions
+from gbasis.integrals.overlap import overlap_integral
+from gbasis.integrals.overlap_n import contracted_n_overlap
+
+
+# =====================================================
+# Helper: Build small hydrogen basis
+# =====================================================
+def build_h2_basis():
+    atoms = ["H", "H"]
+    coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+
+    basis_dict = parse_gbs("tests/data_631g.gbs")
+
+    basis = make_contractions(basis_dict, atoms, coords, coord_types="c")
+
+    return basis
+
+
+# =====================================================
+# 1️⃣ N=2 FULL AO MATRIX VALIDATION
+# =====================================================
+
+
+def test_n2_matches_gbasis():
+
+    basis = build_h2_basis()
+
+    ref = overlap_integral(basis, screen_basis=False)
+
+    n_shells = len(basis)
+    total_ao = ref.shape[0]
+    mine = np.zeros((total_ao, total_ao))
+
+    ao_i = 0
+    for shell_i in basis:
+        size_i = shell_i.num_seg_cont * shell_i.num_cart
+
+        ao_j = 0
+        for shell_j in basis:
+            size_j = shell_j.num_seg_cont * shell_j.num_cart
+
+            block = contracted_n_overlap([shell_i, shell_j])
+            block = block.reshape(size_i, size_j)
+
+            mine[ao_i : ao_i + size_i, ao_j : ao_j + size_j] = block
+
+            ao_j += size_j
+
+        ao_i += size_i
+
+    diff = np.max(np.abs(ref - mine))
+
+    print("N=2 max difference:", diff)
+
+    assert np.allclose(ref, mine, atol=1e-10)
+
+
+# =====================================================
+# 2️⃣ N=3 SYMMETRY TEST
+# =====================================================
+
+
+def test_n3_symmetry():
+
+    basis = build_h2_basis()
+
+    shell1 = basis[0]
+    shell2 = basis[1]
+    shell3 = basis[0]
+
+    S123 = contracted_n_overlap([shell1, shell2, shell3])
+    S321 = contracted_n_overlap([shell3, shell2, shell1])
+
+    diff = np.max(np.abs(S123 - S321))
+
+    print("N=3 symmetry difference:", diff)
+
+    assert np.allclose(S123, S321, atol=1e-10)
+
+
+# =====================================================
+# 3️⃣ HIGH ANGULAR MOMENTUM + DIFFUSE TEST
+# =====================================================
+
+
+def test_high_angmom_diffuse():
+
+    basis = build_h2_basis()
+
+    # choose a higher angular momentum shell if available
+    # try to find p-shell
+    p_shells = [shell for shell in basis if shell.angmom == 1]
+
+    if not p_shells:
+        pytest.skip("No p-shell available in basis.")
+
+    shell = p_shells[0]
+
+    # artificially make exponents diffuse
+    shell.exps = shell.exps * 0.01
+
+    S = contracted_n_overlap([shell, shell])
+
+    print("Diffuse shell max abs value:", np.max(np.abs(S)))
+
+    assert not np.isnan(S).any()
+    assert not np.isinf(S).any()
diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py
index 0e67871b..59fb1eca 100644
--- a/tests/test_nuclear_electron_attraction.py
+++ b/tests/test_nuclear_electron_attraction.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.nuclear_electron_attraction."""
+
 from gbasis.integrals.nuclear_electron_attraction import nuclear_electron_attraction_integral
 from gbasis.integrals.point_charge import point_charge_integral
 from gbasis.parsers import make_contractions, parse_nwchem
diff --git a/tests/test_one_elec_int.py b/tests/test_one_elec_int.py
index ea6d9bfa..451defc8 100644
--- a/tests/test_one_elec_int.py
+++ b/tests/test_one_elec_int.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals._one_elec_int."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals._one_elec_int import _compute_one_elec_integrals
 import numpy as np
diff --git a/tests/test_overlap_asymm.py b/tests/test_overlap_asymm.py
index 1a1167ad..6072925c 100644
--- a/tests/test_overlap_asymm.py
+++ b/tests/test_overlap_asymm.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.overlap_asymm."""
+
 from gbasis.integrals.overlap import overlap_integral
 from gbasis.integrals.overlap_asymm import overlap_integral_asymmetric
 from gbasis.parsers import make_contractions, parse_nwchem, parse_gbs
@@ -119,6 +120,10 @@ def test_overlap_asymmetric_screening_accuracy(precision):
 
     #  the screening tolerance needs to be 1e-4 times the desired precision
     tol_screen = precision * 1e-4
-    overlap_asymmetric = overlap_integral_asymmetric(contraction, contraction, tol_screen=tol_screen)
-    overlap_asymmetric_no_screen = overlap_integral_asymmetric(contraction, contraction, screen_basis=False)
+    overlap_asymmetric = overlap_integral_asymmetric(
+        contraction, contraction, tol_screen=tol_screen
+    )
+    overlap_asymmetric_no_screen = overlap_integral_asymmetric(
+        contraction, contraction, screen_basis=False
+    )
     assert np.allclose(overlap_asymmetric, overlap_asymmetric_no_screen, atol=precision)
diff --git a/tests/test_overlap_n.py b/tests/test_overlap_n.py
new file mode 100644
index 00000000..670a9d21
--- /dev/null
+++ b/tests/test_overlap_n.py
@@ -0,0 +1,42 @@
+import numpy as np
+import pytest
+
+from gbasis.integrals.overlap import Overlap
+from gbasis.integrals.overlap_n import contracted_n_overlap, PrimitiveNEngine
+from gbasis.contractions import GeneralizedContractionShell
+
+
+def build_test_shell(center):
+    ...
+    return GeneralizedContractionShell(...)
+
+
+def test_two_center_matches_existing_overlap():
+    ...
+    assert ...
+
+
+def test_two_center_ss_analytic():
+    ...
+    assert ...
+
+
+def test_primitive_px_py_pz():
+    A = np.array([0.0, 0.0, 0.0])
+    B = np.array([0.3, 0.1, -0.2])
+
+    alpha = 0.5
+
+    print(
+        "px-px:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(1, 0, 0), (1, 0, 0)])
+    )
+
+    print(
+        "py-py:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 1, 0), (0, 1, 0)])
+    )
+
+    print(
+        "pz-pz:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 0, 1), (0, 0, 1)])
+    )
+
+    assert True  # temporary so pytest runs it
diff --git a/tests/test_parsers.py b/tests/test_parsers.py
index 5a42726f..e97ccb84 100644
--- a/tests/test_parsers.py
+++ b/tests/test_parsers.py
@@ -1,4 +1,5 @@
 """Test gbasis.parsers."""
+
 from gbasis.parsers import make_contractions, parse_gbs, parse_nwchem
 import numpy as np
 import pytest
diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py
index 0b4fd6a7..5be39409 100644
--- a/tests/test_point_charge.py
+++ b/tests/test_point_charge.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals.point_charge."""
+
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.integrals.point_charge import point_charge_integral, PointChargeIntegral
 from gbasis.parsers import make_contractions, parse_nwchem
diff --git a/tests/test_spherical.py b/tests/test_spherical.py
index d6d23ad4..ad9d9ec7 100644
--- a/tests/test_spherical.py
+++ b/tests/test_spherical.py
@@ -1,4 +1,5 @@
 """Test gbasis.spherical."""
+
 import itertools as it
 
 from gbasis.spherical import (
diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py
index a9ef3f56..3d9c784d 100644
--- a/tests/test_stress_tensor.py
+++ b/tests/test_stress_tensor.py
@@ -1,4 +1,5 @@
 """Test gbasis.evals.stress_tensor."""
+
 from gbasis.evals.density import (
     evaluate_density_laplacian,
     evaluate_deriv_density,
diff --git a/tests/test_two_elec_int.py b/tests/test_two_elec_int.py
index 8d98617a..fb65538e 100644
--- a/tests/test_two_elec_int.py
+++ b/tests/test_two_elec_int.py
@@ -1,4 +1,5 @@
 """Test gbasis.integrals._two_elec_int."""
+
 from gbasis.integrals._two_elec_int import (
     _compute_two_elec_integrals,
     _compute_two_elec_integrals_angmom_zero,
diff --git a/tests/utils.py b/tests/utils.py
index 7ae48318..da6657d7 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -1,4 +1,5 @@
 """Utility functions for running tests."""
+
 import itertools as it
 import os
 
diff --git a/tools/auto_intor_modify.py b/tools/auto_intor_modify.py
index 881aefce..69f28fb3 100644
--- a/tools/auto_intor_modify.py
+++ b/tools/auto_intor_modify.py
@@ -1,6 +1,5 @@
 import sys
 
-
 # Original string
 old_intor1 = '(gen-cint "intor1.c"'
 
diff --git a/tools/benchmark_density.py b/tools/benchmark_density.py
new file mode 100644
index 00000000..87eafd90
--- /dev/null
+++ b/tools/benchmark_density.py
@@ -0,0 +1,71 @@
+"""
+Benchmark script for density-related integrals in gbasis.
+
+Measures execution time of:
+- compute_intracule
+- compute_extracule
+"""
+
+import time
+import numpy as np
+
+from gbasis.integrals.density import compute_intracule, compute_extracule
+
+
+class DummyShell:
+    """
+    Minimal shell object for benchmarking.
+    """
+
+    def __init__(self):
+        self.exps = np.array([1.0])
+        self.coeffs = np.array([[1.0]])
+        self.norm_prim_cart = np.array([[1.0]])
+
+        self.angmom_components_cart = [(0, 0, 0)]
+        self.angmom_components = [(0, 0, 0)]
+
+        self.coord = np.array([0.0, 0.0, 0.0])
+
+        self.num_seg_cont = 1
+        self.num_cart = 1
+
+
+def benchmark_intracule(shells):
+    start = time.perf_counter()
+    compute_intracule(shells)
+    elapsed = time.perf_counter() - start
+    print(f"intracule: {elapsed:.6f} seconds")
+
+
+def benchmark_extracule(shells):
+    start = time.perf_counter()
+    compute_extracule(shells)
+    elapsed = time.perf_counter() - start
+    print(f"extracule: {elapsed:.6f} seconds")
+
+
+def main(n_shells=20):
+    shells = [DummyShell() for _ in range(n_shells)]
+
+    print(f"Benchmarking with {n_shells} shells\n")
+
+    benchmark_intracule(shells)
+    benchmark_extracule(shells)
+
+
+def main(n_shells=20):
+    shells = [DummyShell() for _ in range(n_shells)]
+
+    print(f"Benchmarking with {n_shells} shells\n")
+
+    start = time.perf_counter()
+    compute_intracule(shells)
+    print(f"intracule: {time.perf_counter() - start:.6f} seconds")
+
+    start = time.perf_counter()
+    compute_extracule(shells)
+    print(f"extracule: {time.perf_counter() - start:.6f} seconds")
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file
diff --git a/tools/benchmark_intracule.py b/tools/benchmark_intracule.py
new file mode 100644
index 00000000..779942cb
--- /dev/null
+++ b/tools/benchmark_intracule.py
@@ -0,0 +1,31 @@
+import time
+import numpy as np
+from gbasis.integrals.density import compute_intracule
+
+
+class DummyShell:
+
+    def __init__(self):
+        self.exps = np.array([1.0])
+        self.coeffs = np.array([[1.0]])
+        self.norm_prim_cart = np.array([[1.0]])
+        self.angmom_components_cart = [(0, 0, 0)]
+        self.angmom_components = [(0, 0, 0)]
+        self.coord = np.array([0.0, 0.0, 0.0])
+        self.num_seg_cont = 1
+        self.num_cart = 1
+
+
+def benchmark_intracule(n_shells=20):
+
+    shells = [DummyShell() for _ in range(n_shells)]
+
+    start = time.perf_counter()
+    compute_intracule(shells)
+    elapsed = time.perf_counter() - start
+
+    print(f"compute_intracule benchmark: {elapsed:.6f} seconds for {n_shells} shells")
+
+
+if __name__ == "__main__":
+    benchmark_intracule()
\ No newline at end of file