Add maximum determinant angular grids

MANIFEST.in

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

src/grid/angular.py

@@ -274,14 +274,18 @@
     52328: 323,
     52978: 325,
 }
+MAX_DET_NPOINTS = {(degree + 1) ** 2: degree for degree in range(1, 200)}
+
 
 # Lebedev/Spherical dictionary of grid's degrees (keys) and numbers of points (values)
 LEBEDEV_DEGREES = dict([(v, k) for k, v in LEBEDEV_NPOINTS.items()])
 SPHERICAL_DEGREES = dict([(v, k) for k, v in SPHERICAL_NPOINTS.items()])
+MAX_DET_DEGREES = dict([(v, k) for k, v in MAX_DET_NPOINTS.items()])
 
 # Cache is used to store the angular grid
 LEBEDEV_CACHE = {}
 SPHERICAL_CACHE = {}
+MAX_DET_CACHE = {}
 
 
 class AngularGrid(Grid):
@@ -332,8 +336,9 @@ 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
+            grids).
 
         Returns
         -------
@@ -353,8 +358,12 @@ def __init__(
             cache_dict = LEBEDEV_CACHE
         elif method == "spherical":
             cache_dict = SPHERICAL_CACHE
+        elif method == "maxdet":
+            cache_dict = MAX_DET_CACHE
         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:
@@ -382,7 +391,11 @@ def __init__(
         self._degree = degree
         # Multiply weights by 4 pi, so that the spherical harmonics are orthonormal,
         #   etc. \int Y_l1 Y_l2 = \delta_{l1, l2}
-        super().__init__(points, weights * 4 * np.pi)
+        # (not necessary for maxdet, already normalized to 4 pi / N)
+        if method == "maxdet":
+            super().__init__(points, weights)
+        else:
+            super().__init__(points, weights * 4 * np.pi)
 
         if method == "lebedev" and np.any(weights < 0.0):
             # Lebedev degrees 13, 25, 27 have negative weights. Symmetric spherical t-design
@@ -448,8 +461,9 @@ 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
+            grids).
 
         Returns
         -------
@@ -464,8 +478,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":
+            dict_degrees, dict_npoints = MAX_DET_DEGREES, MAX_DET_NPOINTS
         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)):
@@ -535,8 +553,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":
+            dict_degrees, dict_npoints = MAX_DET_DEGREES, MAX_DET_NPOINTS
+            file_path = "grid.data.maxdet"
         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/data/maxdet/__init__.py

@@ -0,0 +1,23 @@
+r"""
+Maximum Determinant points on the sphere.
+
+These were obtained from `https://web.maths.unsw.edu.au/~rsw/Sphere/MaxDet/MaxDet1.html`.
+
+The weights are the ones in the data files. For this, the sum of the weights is 4*pi.
+
+References
+----------
+.. [1] I. H. Sloan and R. S. Womersley, "How good can polynomial interpolation on the sphere
+       be?", Advances in Computational Mathematics 14 (2001), 195-226.
+
+.. [2] I. H. Sloan and R. S. Womersley, "Extremal systems of points and numerical
+       integration on the sphere", Advances in Computational Mathematics 21 (2004), 107-125.
+
+.. [3] J. Marzo and J. Ortega-Cerda, "Equidistribution of Fekete Points on the Sphere",
+       Constructive Approximation 32 (2010), 513-521.
+       https://link.springer.com/article/10.1007/s00365-009-9051-5
+
+.. [4] S. V. Borodachov, D. P. Hardin, and E. B. Saff, Discrete Energy on Rectifiable Sets,
+       Springer Monographs in Mathematics, Springer, 2019.
+       https://www.springer.com/gp/book/9780387848075
+"""

src/grid/tests/test_angular.py

@@ -19,7 +19,6 @@
 # --
 """Test lebedev grid."""
 
-
 from unittest import TestCase
 
 import numpy as np
@@ -127,7 +126,7 @@ def test_errors_and_warnings(self):
 
 
 @pytest.mark.parametrize("degree", [5, 10, 100])
-@pytest.mark.parametrize("method", ["lebedev", "spherical"])
+@pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
 def test_integration_of_spherical_harmonic_up_to_degree(degree, method):
     r"""Test integration of spherical harmonic is accurate."""
     grid = AngularGrid(degree=degree, method=method)
@@ -148,7 +147,7 @@ def test_integration_of_spherical_harmonic_up_to_degree(degree, method):
                 assert_almost_equal(0.0, grid.integrate(sph_harm_one[m_ord, :]))
 
 
-@pytest.mark.parametrize("method", ["lebedev", "spherical"])
+@pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
 def test_integration_of_spherical_harmonic_not_accurate_beyond_degree(method):
     r"""Test integration of spherical harmonic of degree higher than grid is not accurate."""
     grid = AngularGrid(degree=3, method=method)
@@ -163,7 +162,7 @@ def test_integration_of_spherical_harmonic_not_accurate_beyond_degree(method):
     assert np.abs(grid.integrate(sph_harm[(6) ** 2, :])) > 1e-8
 
 
-@pytest.mark.parametrize("method", ["lebedev", "spherical"])
+@pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
 def test_orthogonality_of_spherical_harmonic_up_to_degree_three(method):
     r"""Test orthogonality of spherical harmonic up to degree 3 is accurate."""
     degree = 3

src/grid/tests/test_atomgrid.py

@@ -31,7 +31,7 @@
 from scipy.spatial.transform import Rotation as R
 from scipy.integrate import trapezoid
 
-from grid.angular import LEBEDEV_DEGREES, AngularGrid
+from grid.angular import LEBEDEV_DEGREES, SPHERICAL_DEGREES, MAX_DET_DEGREES, AngularGrid
 from grid.atomgrid import AtomGrid, _get_rgrid_size
 from grid.basegrid import Grid, OneDGrid
 from grid.onedgrid import GaussLaguerre, GaussLegendre, UniformInteger
@@ -280,7 +280,7 @@ def test_atomic_rotate(self):
             rot_mt = R.random(random_state=atgrid2.rotate + i).as_matrix()
             assert_allclose(rot_shell, non_rot_shell @ rot_mt)
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_get_shell_grid(self, method):
         """Test angular grid get from get_shell_grid function."""
         rad_pts = np.array([0.1, 0.5, 1])
@@ -334,22 +334,28 @@ def test_convert_points_to_sph(self):
         assert_allclose(np.array([r, theta, phi]).reshape(-1, 3), calc_sph)
 
     @pytest.mark.filterwarnings("ignore:Lebedev weights are negative which can*")
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_spherical_complete(self, method):
         """Test atomic grid consistence for spherical integral."""
-        num_pts = len(LEBEDEV_DEGREES)
-        pts = UniformInteger(num_pts)
+        degree_maps = {
+            "lebedev": LEBEDEV_DEGREES,
+            "spherical": SPHERICAL_DEGREES,
+            "maxdet": MAX_DET_DEGREES,
+        }
+        degree_map = degree_maps[method]
+        num_shells = len(degree_map)
+        pts = UniformInteger(num_shells)
         for _ in range(10):
             start = np.random.rand() * 1e-5
             end = np.random.rand() * 10 + 10
             tf = PowerRTransform(start, end)
             rad_grid = tf.transform_1d_grid(pts)
             atgrid = AtomGrid(
                 rad_grid,
-                degrees=list(LEBEDEV_DEGREES.keys()),
+                degrees=list(degree_map.keys()),
                 method=method,
             )
-            values = np.random.rand(len(LEBEDEV_DEGREES))
+            values = np.random.rand(len(degree_map))
             pt_val = np.zeros(atgrid.size)
             for index, value in enumerate(values):
                 pt_val[atgrid._indices[index] : atgrid._indices[index + 1]] = value
@@ -391,7 +397,7 @@ def helper_func_power_deriv(self, points):
         dzf = 8 * points[:, 2] * points[:, 2] / r
         return dxf + dyf + dzf
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_integrating_angular_components_spherical(self, method):
         """Test integrating angular components of a spherical harmonics of maximum degree 3."""
         odg = OneDGrid(np.array([0.0, 1e-16, 1e-8, 1e-4, 1e-2]), np.ones(5), (0, np.inf))
@@ -463,7 +469,7 @@ def test_integrating_angular_components_with_offcentered_gaussian(self, numb_rad
         assert np.all(integrals[2, :] < 1e-5)  # px, py-orbital should be zero
         assert np.all(integrals[3, :] < 1e-5)
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_fitting_spherical_harmonics(self, method):
         r"""Test fitting the radial components of spherical harmonics is just a one, rest zeros."""
         max_degree = 10  # Maximum degree
@@ -496,7 +502,7 @@ def test_fitting_spherical_harmonics(self, method):
                         assert_almost_equal(radial_comp(radial_pts), 0.0, decimal=8)
                 i += 1
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_fitting_product_of_spherical_harmonic(self, method):
         r"""Test fitting the radial components of r**2 times spherical harmonic."""
         max_degree = 7  # Maximum degree
@@ -567,7 +573,7 @@ def func(sph_points):
                 assert_almost_equal(fit[counter](radial_pts, 2), 0.0)
                 counter += 1
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_value_fitting(self, method):
         """Test spline projection the same as spherical harmonics."""
         odg = OneDGrid(np.arange(10) + 1, np.ones(10), (0, np.inf))
@@ -618,7 +624,7 @@ def test_interpolation_of_gaussian(self, centers):
             self.helper_func_gauss(input_points, centers), interfunc(input_points), atol=1e-3
         )
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_cubicspline_and_interp_mol(self, method):
         """Test cubicspline interpolation values."""
         odg = OneDGrid(np.arange(10) + 1, np.ones(10), (0, np.inf))
@@ -662,7 +668,7 @@ def test_cubicspline_and_interp(self):
                 interp_func = atgrid.interpolate(values)
                 assert_allclose(interp_func(xyz), ref_value)
 
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_spherical_average_of_gaussian(self, method):
         r"""Test spherical average of a Gaussian (radial) function is itself and its integral."""
 
@@ -728,7 +734,7 @@ def func(sph_points):
         assert_allclose(spherical_avg2, 0.0, atol=1e-4)
 
     @pytest.mark.filterwarnings("ignore:Lebedev weights are negative which can*")
-    @pytest.mark.parametrize("method", ["lebedev", "spherical"])
+    @pytest.mark.parametrize("method", ["lebedev", "spherical", "maxdet"])
     def test_interpolate_and_its_derivatives_on_polynomial(self, method):
         """Test interpolation of derivative of polynomial function."""
         odg = OneDGrid(np.linspace(0.01, 10, num=50), np.ones(50), (0, np.inf))