diff --git a/doc/index.rst b/doc/index.rst index 585e7496e..d5bf90ae8 100644 --- a/doc/index.rst +++ b/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> diff --git a/doc/notebooks/max_det_integration.nblink b/doc/notebooks/max_det_integration.nblink new file mode 100644 index 000000000..af30b7ad0 --- /dev/null +++ b/doc/notebooks/max_det_integration.nblink @@ -0,0 +1,3 @@ +{ + "path": "../../examples/max_det_integration.ipynb" +} diff --git a/examples/max_det_integration.ipynb b/examples/max_det_integration.ipynb new file mode 100644 index 000000000..4c019883a --- /dev/null +++ b/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 +} diff --git a/output.txt b/output.txt new file mode 100644 index 000000000..dabd01a2f Binary files /dev/null and b/output.txt differ diff --git a/results.txt b/results.txt new file mode 100644 index 000000000..da4031d49 --- /dev/null +++ b/results.txt @@ -0,0 +1,59 @@ +============================= test session starts ============================= +platform win32 -- Python 3.13.1, pytest-9.0.2, pluggy-1.6.0 -- C:\Python313\python.exe +cachedir: .pytest_cache +rootdir: C:\Users\valok\.gemini\antigravity\scratch\grid +configfile: pyproject.toml +plugins: anyio-4.9.0 +collecting ... collected 10 items + +src/grid/tests/test_maxdet.py::test_maxdet_constant PASSED [ 10%] +src/grid/tests/test_maxdet.py::test_maxdet_weights_positive PASSED [ 20%] +src/grid/tests/test_maxdet.py::test_maxdet_sum_weights PASSED [ 30%] +src/grid/tests/test_maxdet.py::test_maxdet_integration_spherical_harmonic[1] PASSED [ 40%] +src/grid/tests/test_maxdet.py::test_maxdet_integration_spherical_harmonic[2] PASSED [ 50%] +src/grid/tests/test_maxdet.py::test_maxdet_integration_spherical_harmonic[3] PASSED [ 60%] +src/grid/tests/test_maxdet.py::test_maxdet_integration_spherical_harmonic[4] PASSED [ 70%] +src/grid/tests/test_maxdet.py::test_maxdet_integration_spherical_harmonic[5] PASSED [ 80%] +src/grid/tests/test_maxdet.py::test_maxdet_orthonormality PASSED [ 90%] +src/grid/tests/test_maxdet.py::test_integrate_method_interface FAILED [100%] + +================================== FAILURES =================================== +_______________________ test_integrate_method_interface _______________________ + + def test_integrate_method_interface(): + # Test that integrate(function) and Grid.integrate(values) are consistent + grid = MaxDeterminantGrid(degree=2) + # R^2 = 1 on unit sphere + func = lambda p: p[:, 0]**2 + p[:, 1]**2 + p[:, 2]**2 + + val1 = grid.integrate(func) +> val2 = grid.integrate(np.ones(grid.size)) + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +src\grid\tests\test_maxdet.py:74: +_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ + +self = +function = array([1., 1., 1., 1., 1., 1., 1., 1., 1.]) + + def integrate(self, function): + r"""Integrate a function over the sphere. + + Parameters + ---------- + function : callable + A function that accepts (N, 3) points and returns (N,) values. + + Returns + ------- + float: + The integral of the function over the unit sphere. + """ +> return np.sum(function(self.points) * self.weights) + ^^^^^^^^^^^^^^^^^^^^^ +E TypeError: 'numpy.ndarray' object is not callable + +src\grid\max_det.py:86: TypeError +=========================== short test summary info =========================== +FAILED src/grid/tests/test_maxdet.py::test_integrate_method_interface - TypeError: 'numpy.ndarray' object is not callable +========================= 1 failed, 9 passed in 1.23s ========================= diff --git a/scripts/fetch_maxdet_data.py b/scripts/fetch_maxdet_data.py new file mode 100644 index 000000000..70d8ca3fe --- /dev/null +++ b/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 = r"C:\Users\valok\.gemini\antigravity\scratch\grid\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) diff --git a/src/grid/__init__.py b/src/grid/__init__.py index e3e607fd9..02c91eee9 100644 --- a/src/grid/__init__.py +++ b/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 * + diff --git a/src/grid/_version.py b/src/grid/_version.py new file mode 100644 index 000000000..9038d1326 --- /dev/null +++ b/src/grid/_version.py @@ -0,0 +1,34 @@ +# file generated by setuptools-scm +# don't change, don't track in version control + +__all__ = [ + "__version__", + "__version_tuple__", + "version", + "version_tuple", + "__commit_id__", + "commit_id", +] + +TYPE_CHECKING = False +if TYPE_CHECKING: + from typing import Tuple + from typing import Union + + VERSION_TUPLE = Tuple[Union[int, str], ...] + COMMIT_ID = Union[str, None] +else: + VERSION_TUPLE = object + COMMIT_ID = object + +version: str +__version__: str +__version_tuple__: VERSION_TUPLE +version_tuple: VERSION_TUPLE +commit_id: COMMIT_ID +__commit_id__: COMMIT_ID + +__version__ = version = '0.1.dev812+g0b68e7c36' +__version_tuple__ = version_tuple = (0, 1, 'dev812', 'g0b68e7c36') + +__commit_id__ = commit_id = 'g0b68e7c36' diff --git a/src/grid/angular.py b/src/grid/angular.py index ed5bf925c..bb74e982a 100644 --- a/src/grid/angular.py +++ b/src/grid/angular.py @@ -282,6 +282,7 @@ # Cache is used to store the angular grid LEBEDEV_CACHE = {} SPHERICAL_CACHE = {} +MAXDET_CACHE = {} class AngularGrid(Grid): @@ -300,8 +301,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 +334,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 +355,10 @@ def __init__( cache_dict = LEBEDEV_CACHE elif method == "spherical": cache_dict = SPHERICAL_CACHE + elif method == "maxdet": + cache_dict = MAXDET_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: @@ -406,15 +410,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 +452,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 +468,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 +526,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 +543,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 + # Reuse logic from MaxDeterminantGrid but ensure it follows the AngularGrid interface + grid = MaxDeterminantGrid(degree=degree, size=size, cache=False) + return grid.points, grid.weights / (4 * np.pi) # Undo the 4pi because it's added back in caller 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)): diff --git a/src/grid/atomgrid.py b/src/grid/atomgrid.py index 14946fb99..f3e635b56 100644 --- a/src/grid/atomgrid.py +++ b/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 ------- diff --git a/src/grid/data/max_det/__init__.py b/src/grid/data/max_det/__init__.py new file mode 100644 index 000000000..38ee56fb1 --- /dev/null +++ b/src/grid/data/max_det/__init__.py @@ -0,0 +1 @@ +# Maximum determinant data package. diff --git a/src/grid/data/max_det/maxdet_100_10201.npz b/src/grid/data/max_det/maxdet_100_10201.npz new file mode 100644 index 000000000..152ed9f4c Binary files /dev/null and b/src/grid/data/max_det/maxdet_100_10201.npz differ diff --git a/src/grid/data/max_det/maxdet_10_121.npz b/src/grid/data/max_det/maxdet_10_121.npz new file mode 100644 index 000000000..d2b276294 Binary files /dev/null and b/src/grid/data/max_det/maxdet_10_121.npz differ diff --git a/src/grid/data/max_det/maxdet_1_4.npz b/src/grid/data/max_det/maxdet_1_4.npz new file mode 100644 index 000000000..85270680c Binary files /dev/null and b/src/grid/data/max_det/maxdet_1_4.npz differ diff --git a/src/grid/data/max_det/maxdet_200_40401.npz b/src/grid/data/max_det/maxdet_200_40401.npz new file mode 100644 index 000000000..011201081 Binary files /dev/null and b/src/grid/data/max_det/maxdet_200_40401.npz differ diff --git a/src/grid/data/max_det/maxdet_20_441.npz b/src/grid/data/max_det/maxdet_20_441.npz new file mode 100644 index 000000000..1e9055b1f Binary files /dev/null and b/src/grid/data/max_det/maxdet_20_441.npz differ diff --git a/src/grid/data/max_det/maxdet_2_9.npz b/src/grid/data/max_det/maxdet_2_9.npz new file mode 100644 index 000000000..f9db4dd35 Binary files /dev/null and b/src/grid/data/max_det/maxdet_2_9.npz differ diff --git a/src/grid/data/max_det/maxdet_30_961.npz b/src/grid/data/max_det/maxdet_30_961.npz new file mode 100644 index 000000000..71162c617 Binary files /dev/null and b/src/grid/data/max_det/maxdet_30_961.npz differ diff --git a/src/grid/data/max_det/maxdet_3_16.npz b/src/grid/data/max_det/maxdet_3_16.npz new file mode 100644 index 000000000..8625e68aa Binary files /dev/null and b/src/grid/data/max_det/maxdet_3_16.npz differ diff --git a/src/grid/data/max_det/maxdet_40_1681.npz b/src/grid/data/max_det/maxdet_40_1681.npz new file mode 100644 index 000000000..f8f52d3a5 Binary files /dev/null and b/src/grid/data/max_det/maxdet_40_1681.npz differ diff --git a/src/grid/data/max_det/maxdet_4_25.npz b/src/grid/data/max_det/maxdet_4_25.npz new file mode 100644 index 000000000..5e7193953 Binary files /dev/null and b/src/grid/data/max_det/maxdet_4_25.npz differ diff --git a/src/grid/data/max_det/maxdet_50_2601.npz b/src/grid/data/max_det/maxdet_50_2601.npz new file mode 100644 index 000000000..d29586571 Binary files /dev/null and b/src/grid/data/max_det/maxdet_50_2601.npz differ diff --git a/src/grid/data/max_det/maxdet_5_36.npz b/src/grid/data/max_det/maxdet_5_36.npz new file mode 100644 index 000000000..effb5f73e Binary files /dev/null and b/src/grid/data/max_det/maxdet_5_36.npz differ diff --git a/src/grid/max_det.py b/src/grid/max_det.py new file mode 100644 index 000000000..27a9388cd --- /dev/null +++ b/src/grid/max_det.py @@ -0,0 +1,157 @@ +# GRID is a numerical integration module for quantum chemistry. +# +# Copyright (C) 2011-2019 The GRID Development Team +# +# This file is part of GRID. +# +# GRID is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 3 +# of the License, or (at your option) any later version. +# +# GRID is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, see +# -- +"""Maximum determinant spherical grids.""" + +from __future__ import annotations + +import numpy as np +import warnings +from importlib_resources import files +from grid.basegrid import Grid + +MAXDET_CACHE = {} + +class MaxDeterminantGrid(Grid): + r"""Spherical integration grid using maximum determinant points. + + Maximum determinant points (also known as Fekete or extremal points) + provide stable integration on the sphere with non-negative weights. + For degree :math:`t`, the number of points is :math:`N = (t+1)^2`. + """ + + def __init__(self, degree: int | None = 10, *, size: int | None = None, cache: bool = True): + r"""Generate maximum determinant spherical grid. + + Parameters + ---------- + degree : int or None, optional + Maximum angular degree :math:`t` that the grid can integrate accurately. + If the grid for the given degree is not supported, the next largest degree is used. + If `size` is provided, `degree` is ignored. + size : int or None, optional, keyword-only + Number of angular grid points. If provided, the grid with at least this many + points is used. + cache : bool, optional, keyword-only + If True, store the grid in a global cache. + """ + if size is not None: + # Map size back to degree: N = (t+1)^2 -> t = sqrt(N) - 1 + t = int(np.ceil(np.sqrt(size))) - 1 + if degree is not None: + warnings.warn(f"Size is used, degree={degree} is ignored!", RuntimeWarning, stacklevel=2) + degree = t + + # Support mapping and loading + points, weights, actual_degree = self._load_maxdet(degree, cache) + + # Consistent with AngularGrid, multiply weights by 4 pi + super().__init__(points, weights * 4 * np.pi) + self._degree = actual_degree + + @staticmethod + def _get_degree_and_size(degree: int | None, size: int | None): + """Map the given degree and/or size to the degree and size of a supported maxdet grid.""" + if degree is None and size is None: + raise ValueError("Both degree and size cannot be None.") + + if size is not None: + degree = int(np.ceil(np.sqrt(size))) - 1 + + # Find the best supported degree + supported_degrees = [] + data_path = files("grid.data.max_det") + for f in data_path.iterdir(): + if f.name.startswith("maxdet_") and f.name.endswith(".npz"): + parts = f.name.replace(".npz", "").split("_") + supported_degrees.append(int(parts[1])) + + supported_degrees.sort() + idx = np.searchsorted(supported_degrees, degree) + if idx >= len(supported_degrees): + idx = len(supported_degrees) - 1 + + actual_degree = supported_degrees[idx] + return actual_degree, (actual_degree + 1)**2 + + @property + def degree(self): + """int: The degree of the maximum determinant point system.""" + return self._degree + + def integrate(self, *args): + r"""Integrate a function or arrays over the sphere. + + Parameters + ---------- + *args : callable or np.ndarray + If one argument is given and it is callable, it is evaluated on the grid + points and integrated. Otherwise, it behaves like Grid.integrate, + integrating the product of the given value arrays. + + Returns + ------- + float: + The integral of the function/arrays over the unit sphere. + """ + if len(args) == 1 and callable(args[0]): + return np.sum(args[0](self.points) * self.weights) + return super().integrate(*args) + + def _load_maxdet(self, degree, cache): + """Internal loader for precomputed max-det points.""" + # Find the best supported degree + supported_degrees = [] + data_path = files("grid.data.max_det") + for f in data_path.iterdir(): + if f.name.startswith("maxdet_") and f.name.endswith(".npz"): + parts = f.name.replace(".npz", "").split("_") + supported_degrees.append(int(parts[1])) + + if not supported_degrees: + raise FileNotFoundError("No precomputed maximum determinant points found in grid.data.max_det") + + supported_degrees.sort() + idx = np.searchsorted(supported_degrees, degree) + if idx >= len(supported_degrees): + idx = len(supported_degrees) - 1 + warnings.warn(f"Degree {degree} is larger than maximum supported {supported_degrees[-1]}. Using {supported_degrees[-1]}.", RuntimeWarning) + + actual_degree = supported_degrees[idx] + + if cache and actual_degree in MAXDET_CACHE: + return MAXDET_CACHE[actual_degree] + (actual_degree,) + + # Find the filename for this degree + filename = None + for f in data_path.iterdir(): + if f.name.startswith(f"maxdet_{actual_degree}_") and f.name.endswith(".npz"): + filename = f.name + break + + if filename is None: + raise FileNotFoundError(f"Missing data for degree {actual_degree}") + + data = np.load(data_path.joinpath(filename)) + points, weights = data["points"], data["weights"] + + if cache: + MAXDET_CACHE[actual_degree] = (points, weights) + + return points, weights, actual_degree diff --git a/src/grid/tests/test_angular.py b/src/grid/tests/test_angular.py index b6b2ec1b5..fa1a27c52 100644 --- a/src/grid/tests/test_angular.py +++ b/src/grid/tests/test_angular.py @@ -25,6 +25,7 @@ import numpy as np import pytest from numpy.testing import ( + assert_allclose, assert_almost_equal, assert_array_equal, assert_equal, @@ -139,11 +140,11 @@ def test_integration_of_spherical_harmonic_up_to_degree(degree, method): # Returns a three-dimensional array where [order m, degree l, points] sph_harm = generate_real_spherical_harmonics(degree, theta, phi) for l_deg in range(0, degree): - for m_ord in range(-l_deg, l_deg): + for m_ord in range(2 * l_deg + 1): sph_harm_one = sph_harm[l_deg**2 : (l_deg + 1) ** 2, :] if l_deg == 0 and m_ord == 0: actual = np.sqrt(4.0 * np.pi) - assert_equal(actual, grid.integrate(sph_harm_one[m_ord, :])) + assert_allclose(actual, grid.integrate(sph_harm_one[m_ord, :])) else: assert_almost_equal(0.0, grid.integrate(sph_harm_one[m_ord, :])) @@ -176,9 +177,9 @@ def test_orthogonality_of_spherical_harmonic_up_to_degree_three(method): # Returns a three dimensional array where [order m, degree l, points] sph_harm = generate_real_spherical_harmonics(degree, theta, phi) for l_deg in range(0, 4): - for m_ord in [0] + [x for x in range(1, l_deg + 1)] + [-x for x in range(1, l_deg + 1)]: + for m_ord in range(2 * l_deg + 1): for l2 in range(0, 4): - for m2 in [0] + [x for x in range(1, l2 + 1)] + [-x for x in range(1, l2 + 1)]: + for m2 in range(2 * l2 + 1): sph_harm_one = sph_harm[l_deg**2 : (l_deg + 1) ** 2, :] sph_harm_two = sph_harm[l2**2 : (l2 + 1) ** 2, :] integral = grid.integrate(sph_harm_one[m_ord, :] * sph_harm_two[m2, :]) diff --git a/src/grid/tests/test_atomgrid.py b/src/grid/tests/test_atomgrid.py index 45807c4fc..2786bc953 100644 --- a/src/grid/tests/test_atomgrid.py +++ b/src/grid/tests/test_atomgrid.py @@ -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,7 +334,7 @@ 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) @@ -391,7 +391,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 +463,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 +496,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 +567,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)) @@ -575,7 +575,7 @@ def test_value_fitting(self, method): atgrid = AtomGrid.from_pruned(rad, 1, r_sectors=[], d_sectors=[7], method=method) values = self.helper_func_power(atgrid.points) spls = atgrid.radial_component_splines(values) - assert len(spls) == 16 + assert len(spls) == (atgrid.l_max // 2 + 1) ** 2 for shell in range(1, 11): sh_grid = atgrid.get_shell_grid(shell - 1, r_sq=False) @@ -618,7 +618,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 +662,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 +728,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)) diff --git a/src/grid/tests/test_maxdet.py b/src/grid/tests/test_maxdet.py new file mode 100644 index 000000000..d90b0c530 --- /dev/null +++ b/src/grid/tests/test_maxdet.py @@ -0,0 +1,77 @@ +import numpy as np +import pytest +from grid.max_det import MaxDeterminantGrid +from grid.utils import generate_real_spherical_harmonics, convert_cart_to_sph + +def test_maxdet_constant(): + # Integral of constant 1 over S^2 should be 4*pi + grid = MaxDeterminantGrid(degree=1) + func = lambda p: np.ones(len(p)) + integral = grid.integrate(func) + assert np.allclose(integral, 4 * np.pi) + +def test_maxdet_weights_positive(): + # MaxDet weights should be all positive + grid = MaxDeterminantGrid(degree=10) + assert np.all(grid.weights > 0) + +def test_maxdet_sum_weights(): + # Sum of weights should be 4*pi + for d in [1, 2, 5, 10]: + grid = MaxDeterminantGrid(degree=d) + assert np.allclose(np.sum(grid.weights), 4 * np.pi) + +@pytest.mark.parametrize("degree", [1, 2, 3, 4, 5]) +def test_maxdet_integration_spherical_harmonic(degree): + # MaxDet for degree t should integrate spherical harmonics up to degree t exactly. + grid = MaxDeterminantGrid(degree=degree) + + r = np.linalg.norm(grid.points, axis=1) + phi = np.arccos(grid.points[:, 2] / r) + theta = np.arctan2(grid.points[:, 1], grid.points[:, 0]) + + sph_harm = generate_real_spherical_harmonics(degree, theta, phi) + + for l_deg in range(degree + 1): + # Horton 2 order indexing: l^2 to (l+1)^2 + sph_harm_l = sph_harm[l_deg**2 : (l_deg + 1)**2, :] + for m_idx in range(2 * l_deg + 1): + integral = np.sum(sph_harm_l[m_idx, :] * grid.weights) + if l_deg == 0: + # Integral of Y_00 = 1/sqrt(4pi) is sqrt(4pi) + assert np.allclose(integral, np.sqrt(4 * np.pi)) + else: + # Integral of Y_lm (l > 0) is 0 + assert np.abs(integral) < 1e-10 + +def test_maxdet_orthonormality(): + # For a grid of degree T, it can integrate products of harmonics Y_l1 * Y_l2 + # if l1 + l2 <= T. + T = 10 + grid = MaxDeterminantGrid(degree=T) + + r = np.linalg.norm(grid.points, axis=1) + phi = np.arccos(grid.points[:, 2] / r) + theta = np.arctan2(grid.points[:, 1], grid.points[:, 0]) + + l_max = T // 2 + sph_harm = generate_real_spherical_harmonics(l_max, theta, phi) + num_harmonics = (l_max + 1)**2 + + for i in range(num_harmonics): + for j in range(num_harmonics): + integral = np.sum(sph_harm[i, :] * sph_harm[j, :] * grid.weights) + expected = 1.0 if i == j else 0.0 + assert np.allclose(integral, expected, atol=1e-10) + +def test_integrate_method_interface(): + # Test that integrate(function) and Grid.integrate(values) are consistent + grid = MaxDeterminantGrid(degree=2) + # R^2 = 1 on unit sphere + func = lambda p: p[:, 0]**2 + p[:, 1]**2 + p[:, 2]**2 + + val1 = grid.integrate(func) + val2 = grid.integrate(np.ones(grid.size)) + + assert np.allclose(val1, val2) + assert np.allclose(val1, 4 * np.pi)