WoC: Final combined implementation (Weeks 1–5)

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

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
 

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

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

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

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]

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,

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

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,

gbasis/integrals/_moment_int.py

@@ -1,4 +1,5 @@
 """Multipole moment integrals involving Contracted Cartesian Gaussians."""
+
 import numpy as np
 
 

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
 

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
 

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

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))),

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