[PR1] feat: add Boys function for 2 electron and validation tests

gbasis/integrals/boys_functions.py

@@ -0,0 +1,135 @@
+"""Boys functions for two-electron interaction potentials.
+
+This module implements the Boys function, which is the starting point for
+computing two-electron integrals using the Obara-Saika recursion relations.
+
+The Boys function is defined as:
+
+    F_m(T) = integral from 0 to 1 of t^(2m) * exp(-T*t^2) dt
+
+References:
+- Helgaker, T.; Jorgensen, P.; Olsen, J. "Molecular Electronic-Structure Theory"
+  (eq. 9.8.39 for hyp1f1 representation)
+- Ahlrichs, R. "A simple algebraic derivation of the Obara-Saika scheme for
+  general two-electron interaction potentials."
+  Phys. Chem. Chem. Phys., 2006, 8, 3072-3077.
+"""
+
+import numpy as np
+from scipy.special import hyp1f1
+
+
+def boys_function_standard(orders, weighted_dist):
+    r"""Compute standard Boys function for Coulomb potential (1/r12).
+
+    The Coulombic Boys function is defined as:
+
+    .. math::
+
+        F_m(T) = \int_0^1 t^{2m} e^{-T t^2} dt
+
+    This can be expressed in terms of the Kummer confluent hypergeometric
+    function (hyp1f1) as shown in Helgaker (eq. 9.8.39).
+
+    Parameters
+    ----------
+    orders : np.ndarray
+        Differentiation order of the Boys function (m values).
+        Shape can be (M,) or (M, 1, 1, ...) for broadcasting.
+    weighted_dist : np.ndarray
+        Weighted interatomic distances (T values).
+        T = rho * |P - Q|^2 where rho = zeta*eta/(zeta+eta).
+
+    Returns
+    -------
+    boys_eval : np.ndarray
+        Boys function values F_m(T) with shape determined by broadcasting
+        of orders and weighted_dist.
+
+    Notes
+    -----
+    For the Coulomb potential g(r) = 1/r:
+    - G_0(rho, T) = (2*pi/rho) * F_0(T)
+
+    This is the standard case used for electron repulsion integrals.
+
+    There's some documented instability for hyp1f1, mainly for large values.
+    For typical quantum chemistry calculations (m < 20), the values are stable.
+    """
+    return hyp1f1(orders + 0.5, orders + 1.5, -weighted_dist) / (2 * orders + 1)
+
+
+def boys_function_all_orders(m_max, weighted_dist):
+    r"""Compute Boys function for all orders from 0 to m_max.
+
+    Returns F_m(T) for m = 0, 1, ..., m_max simultaneously using
+    scipy's hyp1f1, which is numerically stable for all practical T
+    and m values encountered in quantum chemistry.
+
+    Parameters
+    ----------
+    m_max : int
+        Maximum order of the Boys function needed.
+    weighted_dist : np.ndarray
+        Weighted interatomic distances (T values).
+
+    Returns
+    -------
+    boys_all : np.ndarray
+        Boys function values for all orders from 0 to m_max.
+        Shape: (m_max + 1, *weighted_dist.shape)
+    """
+    # Vectorize hyp1f1 across all orders with broadcasting to avoid
+    # per-order Python loops while keeping numerical behaviour identical
+    # to the standard expression. hyp1f1 is stable for the practical
+    # range m < 20 and T encountered here.
+
+    T = np.asarray(weighted_dist)
+    orders = np.arange(m_max + 1, dtype=np.result_type(T, np.float64))
+
+    # Reshape orders to broadcast over T's dimensions: (m_max+1, 1, 1, ...)
+    orders_shape = (m_max + 1,) + (1,) * T.ndim
+    orders_b = orders.reshape(orders_shape)
+
+    boys_vals = hyp1f1(orders_b + 0.5, orders_b + 1.5, -T) / (2 * orders_b + 1)
+    return boys_vals
+
+
+def get_boys_function(potential="coulomb", omega=None):
+    """Get the appropriate Boys function for a given potential type.
+
+    Parameters
+    ----------
+    potential : str
+        Type of two-electron potential. Options:
+        - "coulomb" or "standard" or "1/r": Standard 1/r12 Coulomb potential
+    omega : float, optional
+        Range-separation or damping parameter. Reserved for future use
+        with screened and damped potentials.
+
+    Returns
+    -------
+    boys_func : callable
+        Boys function with signature boys_func(orders, weighted_dist, rho=None).
+
+    Raises
+    ------
+    ValueError
+        If potential type is not recognized.
+
+    Examples
+    --------
+    >>> boys = get_boys_function("coulomb")
+    >>> F0 = boys(np.array([0]), np.array([1.0]))
+    """
+    potential = potential.lower()
+
+    if potential in ["coulomb", "standard", "1/r"]:
+
+        def boys_func(orders, weighted_dist, rho=None):
+            return boys_function_standard(orders, weighted_dist)
+
+        return boys_func
+
+    else:
+        raise ValueError(f"Unknown potential type: {potential}. " f"Valid options: coulomb")

tests/test_boys_functions.py

@@ -0,0 +1,208 @@
+"""Test gbasis.integrals.boys_functions module."""
+
+import numpy as np
+import pytest
+from scipy.integrate import quad
+from scipy.special import hyp1f1
+
+from gbasis.integrals.boys_functions import (
+    boys_function_all_orders,
+    boys_function_standard,
+    get_boys_function,
+)
+
+
+class TestBoysStandard:
+    """Tests for the standard Boys function (Coulomb potential)."""
+
+    def test_boys_m0_zero(self):
+        """Test F_0(0) = 1."""
+        result = boys_function_standard(np.array([0]), np.array([0.0]))
+        assert np.allclose(result, 1.0)
+
+    def test_boys_m0_limit(self):
+        """Test F_0(T) -> sqrt(pi)/(2*sqrt(T)) for large T."""
+        T = np.array([100.0])
+        result = boys_function_standard(np.array([0]), T)
+        expected = np.sqrt(np.pi) / (2 * np.sqrt(T))
+        assert np.allclose(result, expected, rtol=0.01)
+
+    def test_boys_m1_zero(self):
+        """Test F_1(0) = 1/3."""
+        result = boys_function_standard(np.array([1]), np.array([0.0]))
+        assert np.allclose(result, 1.0 / 3.0)
+
+    def test_boys_general_zero(self):
+        """Test F_m(0) = 1/(2m+1)."""
+        for m in range(10):
+            result = boys_function_standard(np.array([m]), np.array([0.0]))
+            expected = 1.0 / (2 * m + 1)
+            assert np.allclose(result, expected), f"Failed for m={m}"
+
+    def test_boys_consistency_with_hyp1f1(self):
+        """Test that our implementation matches scipy's hyp1f1."""
+        orders = np.arange(5)
+        T_values = np.array([0.1, 0.5, 1.0, 2.0, 5.0])
+
+        for m in orders:
+            for T in T_values:
+                result = boys_function_standard(np.array([m]), np.array([T]))
+                expected = hyp1f1(m + 0.5, m + 1.5, -T) / (2 * m + 1)
+                assert np.allclose(result, expected), f"Failed for m={m}, T={T}"
+
+    def test_boys_broadcasting(self):
+        """Test that Boys function handles broadcasting correctly (shape + values)."""
+        orders = np.array([0, 1, 2])[:, None]
+        T = np.array([0.5, 1.0, 2.0])[None, :]
+        result = boys_function_standard(orders, T)
+        assert result.shape == (3, 3)
+
+        # Verify each element matches individual scalar calls
+        for i, m in enumerate([0, 1, 2]):
+            for j, t in enumerate([0.5, 1.0, 2.0]):
+                expected = boys_function_standard(np.array([m]), np.array([t]))[0]
+                assert np.isclose(result[i, j], expected), \
+                    f"Broadcasting value mismatch at m={m}, T={t}"
+
+    def test_boys_recurrence(self):
+        """Test the recurrence relation: (2m+1)*F_m(T) = 2T*F_{m+1}(T) + exp(-T)."""
+        T = np.array([1.0, 2.0, 3.0])
+        for m in range(5):
+            fm = boys_function_standard(np.array([m]), T)
+            fm1 = boys_function_standard(np.array([m + 1]), T)
+            lhs = (2 * m + 1) * fm
+            rhs = 2 * T * fm1 + np.exp(-T)
+            assert np.allclose(lhs, rhs), f"Recurrence failed for m={m}"
+
+
+class TestBoysAllOrders:
+    """Tests for the all-orders Boys function."""
+
+    def test_matches_standard_small_T(self):
+        """Test that all-orders function matches standard for small T."""
+        m_max = 5
+        T = np.array([0.1, 0.5, 1.0, 2.0])
+        result = boys_function_all_orders(m_max, T)
+
+        for m in range(m_max + 1):
+            expected = boys_function_standard(np.array([m]), T)
+            assert np.allclose(result[m], expected), f"Failed for m={m}"
+
+    def test_matches_standard_large_T(self):
+        """Test that all-orders function matches standard for large T (>= 25)."""
+        m_max = 5
+        T = np.array([30.0, 50.0, 100.0])
+        result = boys_function_all_orders(m_max, T)
+
+        for m in range(m_max + 1):
+            expected = boys_function_standard(np.array([m]), T)
+            np.testing.assert_allclose(result[m], expected.flatten(), rtol=1e-8,
+                err_msg=f"All-orders doesn't match standard for m={m}, large T")
+
+    def test_output_shape(self):
+        """Test that all-orders function returns correct shape (m_max+1, *T.shape)."""
+        m_max = 4
+        T = np.array([0.5, 1.0, 2.0, 5.0, 30.0])
+        result = boys_function_all_orders(m_max, T)
+        assert result.shape == (m_max + 1, T.shape[0]), \
+            f"Expected shape {(m_max + 1, T.shape[0])}, got {result.shape}"
+
+    def test_small_T_series_branch(self):
+        """Exercise the small-T series path and compare with reference."""
+        m_max = 6
+        T = np.array([0.0, 1e-12, 1e-9])
+        result = boys_function_all_orders(m_max, T)
+
+        for m in range(m_max + 1):
+            expected = boys_function_standard(np.array([m]), T)
+            np.testing.assert_allclose(
+                result[m], expected,
+                rtol=1e-12, atol=1e-14,
+                err_msg=f"Small-T mismatch at m={m}")
+
+    def test_multidim_broadcasting(self):
+        """Check broadcasting for multidimensional T arrays."""
+        m_max = 3
+        T = np.array([[0.2, 1.0, 5.0], [0.4, 2.0, 10.0]])[:, None, :]
+        result = boys_function_all_orders(m_max, T)
+
+        assert result.shape == (m_max + 1, *T.shape)
+
+        # Verify values against scalar reference
+        for idx, t_val in np.ndenumerate(T):
+            for m in range(m_max + 1):
+                expected = boys_function_standard(np.array([m]), np.array([t_val]))[0]
+                assert np.isclose(result[m][idx], expected), f"Mismatch at m={m}, idx={idx}"
+
+    def test_recursion_identity(self):
+        """Verify recurrence holds on all-orders output."""
+        m_max = 8
+        T = np.array([0.3, 1.5, 4.0])
+        all_vals = boys_function_all_orders(m_max, T)
+
+        for m in range(m_max):
+            lhs = (2 * m + 1) * all_vals[m]
+            rhs = 2 * T * all_vals[m + 1] + np.exp(-T)
+            np.testing.assert_allclose(lhs, rhs, rtol=1e-11, atol=1e-14,
+                err_msg=f"Recurrence failed for m={m}")
+
+    def test_random_values_against_reference(self):
+        """Random spot-check across orders and T values."""
+        rng = np.random.default_rng(0)
+        T = 10 ** rng.uniform(-10, 4, size=20)
+        m_max = 10
+        result = boys_function_all_orders(m_max, T)
+
+        for m in range(m_max + 1):
+            expected = boys_function_standard(np.array([m]), T)
+            np.testing.assert_allclose(result[m], expected, rtol=1e-11, atol=1e-14,
+                err_msg=f"Random check mismatch at m={m}")
+
+
+class TestGetBoysFunction:
+    """Tests for the get_boys_function factory."""
+
+    def test_coulomb_potential(self):
+        """Test that Coulomb potential returns standard Boys function."""
+        boys = get_boys_function("coulomb")
+        result = boys(np.array([0]), np.array([1.0]))
+        expected = boys_function_standard(np.array([0]), np.array([1.0]))
+        assert np.allclose(result, expected)
+
+    def test_aliases(self):
+        """Test that aliases work correctly."""
+        for alias in ["coulomb", "standard", "1/r"]:
+            boys = get_boys_function(alias)
+            result = boys(np.array([0]), np.array([1.0]))
+            assert result is not None
+
+    def test_unknown_potential(self):
+        """Test that unknown potential raises error."""
+        with pytest.raises(ValueError):
+            get_boys_function("unknown_potential")
+
+
+class TestBoysNumericalIntegration:
+    """Compare Boys function with numerical integration."""
+
+    def test_numerical_integration_m0(self):
+        """Test F_0(T) against numerical integration."""
+
+        def boys_integrand(t, T, m):
+            return t ** (2 * m) * np.exp(-T * t ** 2)
+
+        for T in [0.5, 1.0, 2.0, 5.0]:
+            result_analytic = boys_function_standard(np.array([0]), np.array([T]))[0]
+            result_numeric, _ = quad(boys_integrand, 0, 1, args=(T, 0))
+            assert np.allclose(result_analytic, result_numeric, rtol=1e-6)
+
+    def test_numerical_integration_m2(self):
+        """Test F_2(T) against numerical integration."""
+
+        def boys_integrand(t, T, m):
+            return t ** (2 * m) * np.exp(-T * t ** 2)
+
+        for T in [0.5, 1.0, 2.0]:
+            result_analytic = boys_function_standard(np.array([2]), np.array([T]))[0]
+            result_numeric, _ = quad(boys_integrand, 0, 1, args=(T, 2))
+            assert np.allclose(result_analytic, result_numeric, rtol=1e-6)