Add MaxDeterminantGrid with full AngularGrid and AtomGrid integration

.gitignore

@@ -126,3 +126,6 @@ nosetests.xml
 coverage.xml
 *,cover
 .pytest_cache/
+output.txt
+results.txt
+_version.py

MANIFEST.in

@@ -5,3 +5,5 @@ include grid.data.lebedev/*.npz
 include grid.data.spherical_design/*.npz
 include grid.data.prune_grid/*.npz
 include grid.data.proatoms/*.npz
+include grid.data.max_det/*.npz
+

doc/index.rst

@@ -91,6 +91,7 @@ Modules
     ./notebooks/quickstart
     One-Dimensional Grids <./notebooks/one_dimensional_grids>
     Angular Grids <./notebooks/angular_grid>
+    Max-Det Integration <./notebooks/max_det_integration>
     Atom Grid Construction <./notebooks/atom_grid_construction>
     Atom Grid Application <./notebooks/atom_grid>
     Molecular Grid Construction <./notebooks/molecular_grid_construction>

doc/notebooks/max_det_integration.nblink

@@ -0,0 +1,3 @@
+{
+    "path": "../../examples/max_det_integration.ipynb"
+}

examples/max_det_integration.ipynb

@@ -0,0 +1,95 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Spherical Integration with Maximum Determinant Points\n",
+    "\n",
+    "This example demonstrates how to use the `MaxDeterminantGrid` class for numerical integration on the unit sphere.\n",
+    "\n",
+    "## 1. Introduction\n",
+    "Maximum determinant points (Fekete points) provide a stable set of points for integration on the sphere with positive weights. \n",
+    "For a given degree $t$, there are $N = (t+1)^2$ points."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from grid.max_det import MaxDeterminantGrid\n",
+    "from grid.utils import generate_real_spherical_harmonics\n",
+    "\n",
+    "# Create a Max-Det grid of degree 5\n",
+    "grid = MaxDeterminantGrid(degree=5)\n",
+    "print(f\"Grid Size: {grid.size}\")\n",
+    "print(f\"Degrees Supported: {grid.degree}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 2. Integration of a simple function\n",
+    "Let's integrate $f(x, y, z) = x^2 + y^2 + z^2$ over the unit sphere. \n",
+    "Since $x^2 + y^2 + z^2 = 1$ on the sphere, the integral should be $4\\pi$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def func(p):\n",
+    "    return p[:, 0]**2 + p[:, 1]**2 + p[:, 2]**2\n",
+    "\n",
+    "integral = grid.integrate(func)\n",
+    "print(f\"Calculated Integral: {integral}\")\n",
+    "print(f\"Expected Integral:   {4 * np.pi}\")\n",
+    "print(f\"Error:               {np.abs(integral - 4 * np.pi)}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## 3. Integration of Spherical Harmonics\n",
+    "Spherical harmonics of degree $l > 0$ should integrate to 0."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r = np.linalg.norm(grid.points, axis=1)\n",
+    "phi = np.arccos(grid.points[:, 2] / r)\n",
+    "theta = np.arctan2(grid.points[:, 1], grid.points[:, 0])\n",
+    "\n",
+    "l = 2; m = 0\n",
+    "sph_harm = generate_real_spherical_harmonics(l, theta, phi)\n",
+    "Y_lm = sph_harm[l**2, :] # Horton 2 index for (l=2, m=0)\n",
+    "\n",
+    "integral_Y = grid.integrate(Y_lm)\n",
+    "print(f\"Integral of Y_{l}{m}: {integral_Y}\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

scripts/fetch_maxdet_data.py

@@ -0,0 +1,43 @@
+import urllib.request
+import os
+import numpy as np
+
+def download_maxdet(degree, npts):
+    filename = f"md{degree:03d}.{npts:05d}"
+    url = f"https://web.maths.unsw.edu.au/~rsw/Sphere/S2Pts/MD/{filename}"
+    print(f"Downloading {url}...")
+    try:
+        with urllib.request.urlopen(url) as f:
+            data = f.read().decode().splitlines()
+        points_weights = []
+        for line in data:
+            if line.strip():
+                points_weights.append([float(x) for x in line.split()])
+        pw = np.array(points_weights)
+        points = pw[:, :3]
+        # Sum of weights in the file is already 4*pi
+        # Let's normalize to 1 to stay consistent with other npz files
+        weights = pw[:, 3] / (4 * np.pi)
+        return points, weights
+    except Exception as e:
+        print(f"Error downloading {url}: {e}")
+        return None, None
+
+def save_as_npz(degree, npts, points, weights, target_dir):
+    filename = f"maxdet_{degree}_{npts}.npz"
+    filepath = os.path.join(target_dir, filename)
+    print(f"Saving to {filepath}...")
+    np.savez(filepath, degree=degree, size=npts, points=points, weights=weights)
+
+if __name__ == "__main__":
+    target_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "src", "grid", "data", "max_det"))
+    if not os.path.exists(target_dir):
+        os.makedirs(target_dir)
+
+    # Degrees and corresponding number of points
+    degrees = [1, 2, 3, 4, 5, 10, 20, 30] 
+    for d in degrees:
+        npts = (d + 1) ** 2
+        p, w = download_maxdet(d, npts)
+        if p is not None:
+            save_as_npz(d, npts, p, w, target_dir)

src/grid/__init__.py

@@ -34,4 +34,6 @@
 from grid.onedgrid import *
 from grid.periodicgrid import *
 from grid.rtransform import *
+from grid.max_det import *
 from grid.ngrid import *
+

src/grid/angular.py

@@ -300,8 +300,9 @@ class AngularGrid(Grid):
     :math:`4\pi` normalization factor present in the original quadrature scheme
     is included in the integration weights.
 
-    Two types of angular grids are supported: Lebedev-Laikov grid and symmetric
-    spherical t-design. Specifically, for spherical t-design, the weights are constant
+    Two types of angular grids are supported: Lebedev-Laikov grid, symmetric
+    spherical t-design, and maximum determinant points. 
+    Specifically, for spherical t-design, the weights are constant
     value of :math:`4 \pi / N`\, where :math:`N` is the number of points in the grid.
     The weights are chosen so that the spherical harmonics are normalized.
 
@@ -332,8 +333,8 @@ def __init__(
             If True, then store the points and weights of the AngularGrid in cache
             to avoid duplicate grids that have the same `degree`\.
         method: str, optional, keyword-only
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Returns
         -------
@@ -353,8 +354,10 @@ def __init__(
             cache_dict = LEBEDEV_CACHE
         elif method == "spherical":
             cache_dict = SPHERICAL_CACHE
+        elif method == "maxdet":
+            cache_dict = None
         else:
-            raise ValueError(f"Method {method} is not supported, choose 'lebedev' or 'spherical'")
+            raise ValueError(f"Method {method} is not supported, choose 'lebedev', 'spherical' or 'maxdet'")
 
         # allow only one of degree or size to be given
         if size is not None:
@@ -369,16 +372,16 @@ def __init__(
         # map degree and size to the supported (i.e., pre-computed) degree and size
         degree, size = self._get_degree_and_size(degree=degree, size=size, method=method)
         # load pre-computed angular points & weights and make angular grid
-        if degree not in cache_dict:
+        if cache_dict is not None and degree in cache_dict:
+            points, weights = cache_dict[degree]
+        else:
             points, weights = self._load_precomputed_angular_grid(degree, size, method)
             # store the points and weights in cache_dict which updates the global cache
             # dictionary for the given method (i.e., LEBEDEV_CACHE or SPHERICAL_CACHE)
             # in this case, if another instance of AngularGrid is created with the same degree,
             # the points and weights are not recomputed.
-            if cache:
+            if cache and cache_dict is not None:
                 cache_dict[degree] = points, weights
-        else:
-            points, weights = cache_dict[degree]
         self._degree = degree
         # Multiply weights by 4 pi, so that the spherical harmonics are orthonormal,
         #   etc. \int Y_l1 Y_l2 = \delta_{l1, l2}
@@ -406,15 +409,15 @@ def method(self):
     @staticmethod
     def convert_angular_sizes_to_degrees(sizes: np.ndarray, method: str):
         """
-        Convert given Lebedev/Spherical design grid sizes to degrees.
+        Convert given Lebedev/Spherical/Max-Det design grid sizes to degrees.
 
         Parameters
         ----------
         sizes : ndarray[int]
             Sequence of angular grid sizes (e.g., number of points for each atomic shell).
         method: str
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Returns
         -------
@@ -448,8 +451,8 @@ def _get_degree_and_size(degree: int | None, size: int | None, method: str):
             degree is used for constructing the grid. Use None, if `degree` is given. If both
             `degree` and `size` are given `degree` is used for constructing the grid.
         method: str
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Returns
         -------
@@ -464,8 +467,12 @@ def _get_degree_and_size(degree: int | None, size: int | None, method: str):
             dict_degrees, dict_npoints = LEBEDEV_DEGREES, LEBEDEV_NPOINTS
         elif method == "spherical":
             dict_degrees, dict_npoints = SPHERICAL_DEGREES, SPHERICAL_NPOINTS
+        elif method == "maxdet":
+            # For maxdet, N = (t+1)^2
+            from grid.max_det import MaxDeterminantGrid
+            return MaxDeterminantGrid._get_degree_and_size(degree, size)
         else:
-            raise ValueError(f"Method {method} is not supported, choose 'lebedev' or 'spherical'")
+            raise ValueError(f"Method {method} is not supported, choose 'lebedev', 'spherical' or 'maxdet'")
 
         # check whether degree and size are valid
         if not (degree is None or (isinstance(degree, (int, np.integer)) and degree >= 0)):
@@ -518,8 +525,8 @@ def _load_precomputed_angular_grid(degree: int, size: int, method: str):
         size : int
             Number of angular grid points.
         method: str
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Returns
         -------
@@ -535,8 +542,13 @@ def _load_precomputed_angular_grid(degree: int, size: int, method: str):
         elif method == "spherical":
             dict_degrees, dict_npoints = SPHERICAL_DEGREES, SPHERICAL_NPOINTS
             file_path = "grid.data.spherical_design"
+        elif method == "maxdet":
+            from grid.max_det import MaxDeterminantGrid
+            # MaxDeterminantGrid stores raw (unnormalized) weights
+            grid = MaxDeterminantGrid(degree=degree, size=size, cache=False)
+            return grid.points, grid.weights
         else:
-            raise ValueError(f"Method {method} is not supported, choose 'lebedev' or 'spherical'")
+            raise ValueError(f"Method {method} is not supported, choose 'lebedev', 'spherical' or 'maxdet'")
 
         # check whether degree and size are valid
         if not (degree is None or (isinstance(degree, (int, np.integer)) and degree >= 0)):

src/grid/atomgrid.py

@@ -89,8 +89,8 @@ def __init__(
             spherical grids at each radial grid point. If the integer is zero, then no rotate
             is used.
         method: str, optional, keyword-only
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         """
         # check stage, if center is None, set to (0., 0., 0.)
@@ -172,8 +172,8 @@ def from_preset(
             Integer used as a seed for generating random rotation matrices to rotate the angular
             spherical grids at each radial grid point. If 0, then no rotate is made.
         method: str, optional
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Notes
         -----
@@ -279,8 +279,8 @@ def from_pruned(
             spherical grids at each radial grid point. If the integer is zero, then no rotate
             is used.
         method: str, optional
-            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov)
-            and "spherical" (for symmetric spherical t-design).
+            Method for constructing the angular grid. Options are "lebedev" (for Lebedev-Laikov),
+            "spherical" (for symmetric spherical t-design), and "maxdet" (for maximum determinant).
 
         Returns
         -------

src/grid/data/max_det/__init__.py

@@ -0,0 +1 @@
+# Maximum determinant data package.