add: make (parts) of sympleint testable

add: basic Boys-function implementation for testing
add: int2c-driver for testing
add: tests against pyscf/libcint
add: cart2sph module from pysisyphus to remove dependecy on it
add: improved generation of spherical integrals by using more precise
transformation factors generated using mpmath
add: sympleints.defs.overlap.gen_prefactor_shell for GDMA

ressources/cgto_normalization.ipynb

@@ -661,7 +661,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,

sympleints/Functions.py

@@ -26,7 +26,7 @@ class Functions:
     header: str = ""
     doc_func: Optional[Callable] = None
     comment: str = ""
-    boys: bool = False
+    boys_func: Optional[str] = None
     ncomponents: int = 0
     with_ref_center: bool = True
     full_name: Optional["str"] = None
@@ -141,7 +141,8 @@ def ls(self):
 
     @property
     def nbfs(self):
-        return len(self.coeffs)
+        # Return number of total angular momenta
+        return len(self.ls_exprs[0][0])
 
     @property
     def ls_name_map(self):
@@ -170,8 +171,3 @@ def ndim_act(self):
     @property
     def cartesian(self):
         return not self.spherical
-
-    # def numba_func_type(self):
-    # pass
-
-    # def numba_driver_type(self):

sympleints/PythonRenderer.py

@@ -81,6 +81,7 @@ def render_equi_function(
         name,
         equi_name,
         equi_inds,
+        reshape,
         from_axes,
         to_axes,
     ):
@@ -96,6 +97,7 @@ def render_equi_function(
             args=args,
             from_axes=from_axes,
             to_axes=to_axes,
+            reshape=reshape,
         )
         return rendered
 
@@ -110,7 +112,7 @@ def render_module(self, functions, rendered_funcs, **tpl_kwargs):
         _tpl_kwargs = {
             "header": functions.header,
             "comment": functions.comment,
-            "boys": functions.boys,
+            "boys": functions.boys_func,
             "funcs": rendered_funcs,
             "func_dict": func_dict,
         }

sympleints/Renderer.py

@@ -64,6 +64,7 @@ def render_equi_function(
         name,
         act_name,
         equi_inds,
+        L_tots,
         from_axes,
         to_axes,
     ):
@@ -133,16 +134,21 @@ def render_functions(self, functions: Functions):
                     from_axes = tuple(from_axes)
 
                 more_components = functions.ncomponents > 1
-                print(
-                    f"{more_components=}, {functions.ncomponents=}, {self._drop_dim=}"
-                )
                 add_axis = more_components or (not self._drop_dim and nbfs == 2)
-                print(f"{add_axis=}")  # , {inds_missing=}")
                 if add_axis:
                     from_axes = tuple([0] + [fa + 1 for fa in from_axes])
                     to_axes = tuple([0] + [ta + 1 for ta in to_axes])
+
+                reshape = shape if functions.ncomponents > 1 else shape[1:]
+
                 func_equi = self.render_equi_function(
-                    functions, name, name_equi, equi_inds, from_axes, to_axes
+                    functions,
+                    name,
+                    name_equi,
+                    equi_inds,
+                    reshape,
+                    from_axes,
+                    to_axes,
                 )
                 rendered_funcs.append(
                     RenderedFunction(name=name_equi, Ls=L_tots_equi, text=func_equi)

sympleints/boys.py

@@ -0,0 +1,147 @@
+from typing import Callable
+
+import numpy as np
+import scipy as sp
+from scipy.special import factorial, factorial2
+
+
+def boys_quad(n, x, err=1e-13):
+    """Boys function from quadrature."""
+
+    def inner(t):
+        return t ** (2 * n) * np.exp(-x * t**2)
+
+    y, err = sp.integrate.quad(inner, 0.0, 1.0, epsabs=err, epsrel=err)
+    return y
+
+
+def get_table(nmax, order, dx, xmax):
+    nsteps = (xmax / dx) + 1
+    # assert that mantissa of nsteps is 0.0
+    assert (
+        float(int(nsteps)) == nsteps
+    ), f"{xmax=} must be an integer multiple of {dx=:.8e} but it isn't!"
+    nsteps = int(nsteps)
+    xs = np.linspace(0.0, xmax, num=nsteps)
+    table = np.empty_like(xs)
+    for i, x in enumerate(xs):
+        table[i] = boys_quad(nmax + order, x)
+    return xs, table
+
+
+def boys_x0(n):
+    """Boys-function when x equals 0.0."""
+    return 1 / (2 * n + 1)
+
+
+def boys_xlarge(n, x):
+    """Boys-function for large x values (x >= 30 is sensible)."""
+    _2n = 2 * n
+    f2 = factorial2(_2n - 1) if n > 0 else 1.0
+    return f2 / 2 ** (n + 1) * np.sqrt(np.pi / x ** (_2n + 1))
+
+
+def get_boys_func(
+    nmax: int, order: int = 6, dx: float = 0.1, xlarge: float = 30.0
+) -> Callable[[int, float], float]:
+    """Wrapper that returns Boys function object.
+
+               1
+               /
+              |
+              |          2
+    F_n(x) =  |   2*n  -t x
+              |  t    e      dt
+              |
+             /
+             0
+
+     Parameters
+     ----------
+     nmax
+         Positive integer Maximum n value for which the Boys function can be evaluated.
+     order
+         Order of the Tayor-expansion that is used to evaluate the Boys-function, defaults
+         to 6.
+     dx
+        Positive float; Distance between points on the grid.
+     xlarge
+        Cutoff value. For 0 < x < xlarge a Taylor expansion will be used to evaluate
+        the boys function. For x >= xlarge the Boys-function is evaluated using the
+        incomplete gamma function. For all n and x = 0.0 the integral is calculated
+        analytically.
+
+    Returns
+    -------
+    boys
+        Callable F(n, x) that accepts two arguments: an integer 0 <= n <= nmax and a
+        float 0.0 <= x. It returns the values of the Boys-function F_n at the given x.
+    """
+    # Pretabulated Boys-values on a grid for n = nmax + order
+    xs, table = get_table(nmax, order, dx, xlarge)
+    nxs = len(xs)
+
+    # Build up full table for all n = {0, 1, ... nmax + order} via Downward recursion
+    table_full = np.empty((nmax + order + 1, nxs))
+    table_full[-1] = table
+    exp_minus_xs = np.exp(-xs)
+    _2xs = 2 * xs
+    # Loop over all n-values
+    for n in range(nmax + order - 1, -1, -1):
+        table_full[n] = (_2xs * table_full[n + 1] + exp_minus_xs) / (2 * n + 1)
+    # 1 / k!
+    _1_factorials = 1 / factorial(np.arange(order, dtype=int))
+
+    # Prefactors for Boys-values for large x > xlarge
+    xlarge_prefact = np.empty(nmax + 1)
+    for n in range(nmax + 1):
+        xlarge_prefact[n] = (
+            (factorial2(2 * n - 1) if n > 0 else 1.0) / 2 ** (n + 1) * np.sqrt(np.pi)
+        )
+
+    def boys_xlarge_precalc(n, x):
+        return xlarge_prefact[n] * x ** (-n - 0.5)
+
+    def boys_taylor(n, x):
+        # Rounding to the nearest integer always yields a grid point before (or on) x
+        # so we interpolate and don't extrapolate
+        factor = int(1 / dx)
+        ind0 = int(x * factor)
+        # Closest grid point (on/before) x
+        xg = ind0 * dx
+
+        # Step length for Taylor expansion; depends on closest grid point
+        dx_taylor = x - xg
+        result = 0.0
+        # The loop is somehow faster than the vectorized expression
+        for k in range(order):
+            result += table_full[n + k, ind0] * (-dx_taylor) ** k * _1_factorials[k]
+        return result
+
+    def boys(n: int, xs: float) -> float:
+        def func(n, x):
+            if x == 0.0:
+                # Take values directly from table; should correspond to 1 / (2n + 1)
+                return table_full[n, 0]
+            elif x < xlarge:
+                return boys_taylor(n, x)
+            else:
+                return boys_xlarge_precalc(n, x)
+
+        if isinstance(xs, np.ndarray):
+            boys_list = list()
+
+            with np.nditer(xs) as it:
+                for x in it:
+                    boys_list.append(func(n, x))
+            boys_table = np.reshape(boys_list, xs.shape)
+            return boys_table
+        else:
+            return func(n, xs)
+
+    return boys
+
+
+# Calculation should take ~ 20 ms
+# python -m timeit -n 25 "from sympleints.boys import boys"
+boys = get_boys_func(nmax=8)