build against NPY2

pyproject.toml

@@ -1,8 +1,9 @@
 [build-system]
 requires = [
     "Cython>=0.29.23",
-    # pin NumPy version used in the build
-    "oldest-supported-numpy",
+    # Building against NPY2 will support both NPY1 and NPY2
+    # https://numpy.org/devdocs/dev/depending_on_numpy.html#build-time-dependency
+    "numpy>=2.0.1",
     "setuptools>=65.0.0",
 ]
 build-backend = "setuptools.build_meta"
@@ -59,8 +60,8 @@ dependencies = [
     "matplotlib>=3.8",
     "monty>=2024.7.29",
     "networkx>=2.2",
-    "numpy>=1.25.0 ; platform_system != 'Windows'",
-    "numpy>=1.25.0,<2.0 ; platform_system == 'Windows'",
+    # NumPy documentation suggests pinning the current major version as the C API is used
+    # https://numpy.org/devdocs/dev/depending_on_numpy.html#runtime-dependency-version-ranges
     "palettable>=3.3.3",
     "pandas>=2",
     "plotly>=4.5.0",
@@ -73,6 +74,7 @@ dependencies = [
     "tabulate>=0.9",
     "tqdm>=4.60",
     "uncertainties>=3.1.4",
+    'numpy>=1.25.0,<3.0',
 ]
 version = "2024.7.18"
 

src/pymatgen/analysis/eos.py

@@ -15,6 +15,11 @@
 import numpy as np
 from scipy.optimize import leastsq, minimize
 
+try:
+    from numpy.exceptions import RankWarning  # NPY2
+except ImportError:
+    from numpy import RankWarning  # NPY1
+
 from pymatgen.core.units import FloatWithUnit
 from pymatgen.util.plotting import add_fig_kwargs, get_ax_fig, pretty_plot
 
@@ -420,7 +425,7 @@ def fit(self, min_ndata_factor=3, max_poly_order_factor=5, min_poly_order=2):
             min_poly_order (int): minimum order of the polynomial to be
                 considered for fitting.
         """
-        warnings.simplefilter("ignore", np.RankWarning)
+        warnings.simplefilter("ignore", RankWarning)
 
         def get_rms(x, y):
             return np.sqrt(np.sum((np.array(x) - np.array(y)) ** 2) / len(x))
@@ -490,7 +495,7 @@ def get_rms(x, y):
             norm += weight
             coeffs = np.array(val[0])
             # pad the coefficient array with zeros
-            coeffs = np.lib.pad(coeffs, (0, max(fit_poly_order - len(coeffs), 0)), "constant")
+            coeffs = np.pad(coeffs, (0, max(fit_poly_order - len(coeffs), 0)), "constant")
             weighted_avg_coeffs += weight * coeffs
 
         # normalization

src/pymatgen/optimization/linear_assignment.pyx

@@ -50,7 +50,8 @@ class LinearAssignment:
     """
 
     def __init__(self, costs, epsilon=1e-13):
-        self.orig_c = np.array(costs, dtype=np.float_, copy=False, order="C")
+        # https://numpy.org/devdocs/numpy_2_0_migration_guide.html#adapting-to-changes-in-the-copy-keyword
+        self.orig_c = np.asarray(costs, dtype=np.float64, order="C")
         self.nx, self.ny = self.orig_c.shape
         self.n = self.ny
 
@@ -63,7 +64,7 @@ class LinearAssignment:
         if self.nx == self.ny:
             self.c = self.orig_c
         else:
-            self.c = np.zeros((self.n, self.n), dtype=np.float_)
+            self.c = np.zeros((self.n, self.n), dtype=np.float64)
             self.c[:self.nx] = self.orig_c
 
         # initialize solution vectors
@@ -76,15 +77,15 @@ class LinearAssignment:
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-cdef np.float_t compute(int size, np.float_t[:, :] c, np.int_t[:] x, np.int_t[:] y, np.float_t eps) nogil:
+cdef np.float_t compute(int size, np.float_t[:, :] c, np.int64_t[:] x, np.int64_t[:] y, np.float_t eps) nogil:
 
     # augment
     cdef int i, j, k, i1, j1, f, f0, cnt, low, up, z, last, nrr
     cdef int n = size
     cdef bint b
-    cdef np.int_t * col = <np.int_t *> malloc(n * sizeof(np.int_t))
-    cdef np.int_t * fre = <np.int_t *> malloc(n * sizeof(np.int_t))
-    cdef np.int_t * pred = <np.int_t *> malloc(n * sizeof(np.int_t))
+    cdef np.int64_t * col = <np.int64_t *> malloc(n * sizeof(np.int64_t))
+    cdef np.int64_t * fre = <np.int64_t *> malloc(n * sizeof(np.int64_t))
+    cdef np.int64_t * pred = <np.int64_t *> malloc(n * sizeof(np.int64_t))
     cdef np.float_t * v = <np.float_t *> malloc(n * sizeof(np.float_t))
     cdef np.float_t * d = <np.float_t *> malloc(n * sizeof(np.float_t))
     cdef np.float_t h, m, u1, u2, cost

src/pymatgen/optimization/neighbors.pyx

@@ -45,7 +45,7 @@ def find_points_in_spheres(
         const double[:, ::1] all_coords,
         const double[:, ::1] center_coords,
         const double r,
-        const long[::1] pbc,
+        const np.int64_t[::1] pbc,
         const double[:, ::1] lattice,
         const double tol=1e-8,
         const double min_r=1.0):
@@ -57,7 +57,7 @@ def find_points_in_spheres(
             When periodic boundary is considered, this is all the points in the lattice.
         center_coords: (np.ndarray[double, dim=2]) all centering points
         r: (float) cutoff radius
-        pbc: (np.ndarray[long, dim=1]) whether to set periodic boundaries
+        pbc: (np.ndarray[np.int64_t, dim=1]) whether to set periodic boundaries
         lattice: (np.ndarray[double, dim=2]) 3x3 lattice matrix
         tol: (float) numerical tolerance
         min_r: (float) minimal cutoff to calculate the neighbor list
@@ -87,10 +87,10 @@ def find_points_in_spheres(
 
         int n_center = center_coords.shape[0]
         int n_total = all_coords.shape[0]
-        long nlattice = 1
+        np.int64_t nlattice = 1
 
-        long[3] max_bounds = [1, 1, 1]
-        long[3] min_bounds = [0, 0, 0]
+        np.int64_t[3] max_bounds = [1, 1, 1]
+        np.int64_t[3] min_bounds = [0, 0, 0]
         double [:, ::1] frac_coords =  <double[:n_center, :3]> safe_malloc(
             n_center * 3 * sizeof(double)
         )
@@ -114,23 +114,23 @@ def find_points_in_spheres(
         double *expanded_coords_p_temp = <double*> safe_malloc(
             n_atoms * 3 * sizeof(double)
         )
-        long *indices_p_temp = <long*> safe_malloc(n_atoms * sizeof(long))
+        np.int64_t *indices_p_temp = <np.int64_t*> safe_malloc(n_atoms * sizeof(np.int64_t))
         double coord_temp[3]
-        long ncube[3]
+        np.int64_t ncube[3]
 
-        long[:, ::1] center_indices3 = <long[:n_center, :3]> safe_malloc(
-            n_center*3*sizeof(long)
+        np.int64_t[:, ::1] center_indices3 = <np.int64_t[:n_center, :3]> safe_malloc(
+            n_center*3*sizeof(np.int64_t)
         )
-        long[::1] center_indices1 = <long[:n_center]> safe_malloc(n_center*sizeof(long))
+        np.int64_t[::1] center_indices1 = <np.int64_t[:n_center]> safe_malloc(n_center*sizeof(np.int64_t))
 
         int malloc_chunk = 10000  # size of memory chunks to re-allocate dynamically
         int failed_malloc = 0  # flag for failed reallocation within loops
-        long *index_1 = <long*> safe_malloc(malloc_chunk*sizeof(long))
-        long *index_2 = <long*> safe_malloc(malloc_chunk*sizeof(long))
+        np.int64_t *index_1 = <np.int64_t*> safe_malloc(malloc_chunk*sizeof(np.int64_t))
+        np.int64_t *index_2 = <np.int64_t*> safe_malloc(malloc_chunk*sizeof(np.int64_t))
         double *offset_final = <double*> safe_malloc(3*malloc_chunk*sizeof(double))
         double *distances = <double*> safe_malloc(malloc_chunk*sizeof(double))
-        long cube_index_temp
-        long link_index
+        np.int64_t cube_index_temp
+        np.int64_t link_index
         double d_temp2
         double r2 = r * r
 
@@ -201,8 +201,8 @@ def find_points_in_spheres(
                             expanded_coords_p_temp = <double*> realloc(
                                 expanded_coords_p_temp, n_atoms * 3 * sizeof(double)
                             )
-                            indices_p_temp = <long*> realloc(
-                                indices_p_temp, n_atoms * sizeof(long)
+                            indices_p_temp = <np.int64_t*> realloc(
+                                indices_p_temp, n_atoms * sizeof(np.int64_t)
                             )
                         if (
                                 offset_final is NULL or
@@ -256,38 +256,38 @@ def find_points_in_spheres(
         double *expanded_coords_p = <double*> safe_realloc(
             expanded_coords_p_temp, count * 3 * sizeof(double)
         )
-        long *indices_p = <long*> safe_realloc(
-            indices_p_temp, count * sizeof(long)
+        np.int64_t *indices_p = <np.int64_t*> safe_realloc(
+            indices_p_temp, count * sizeof(np.int64_t)
         )
 
         double[:, ::1] offsets = <double[:count, :3]> offsets_p
         double[:, ::1] expanded_coords = <double[:count, :3]> expanded_coords_p
-        long[::1] indices = <long[:count]> indices_p
+        np.int64_t[::1] indices = <np.int64_t[:count]> indices_p
 
         # Construct linked cell list
-        long[:, ::1] all_indices3 = <long[:n_atoms, :3]> safe_malloc(
-            n_atoms * 3 * sizeof(long)
+        np.int64_t[:, ::1] all_indices3 = <np.int64_t[:n_atoms, :3]> safe_malloc(
+            n_atoms * 3 * sizeof(np.int64_t)
         )
-        long[::1] all_indices1 = <long[:n_atoms]> safe_malloc(
-            n_atoms * sizeof(long)
+        np.int64_t[::1] all_indices1 = <np.int64_t[:n_atoms]> safe_malloc(
+            n_atoms * sizeof(np.int64_t)
         )
 
     for i in range(3):
-        ncube[i] = <long>(ceil((valid_max[i] - valid_min[i]) / ledge))
+        ncube[i] = <np.int64_t>(ceil((valid_max[i] - valid_min[i]) / ledge))
 
     compute_cube_index(expanded_coords, valid_min, ledge, all_indices3)
     three_to_one(all_indices3, ncube[1], ncube[2], all_indices1)
 
     cdef:
-        long nb_cubes = ncube[0] * ncube[1] * ncube[2]
-        long *head = <long*> safe_malloc(nb_cubes*sizeof(long))
-        long *atom_indices = <long*> safe_malloc(n_atoms*sizeof(long))
-        long[:, ::1] neighbor_map = <long[:nb_cubes, :27]> safe_malloc(
-            nb_cubes * 27 * sizeof(long)
+        np.int64_t nb_cubes = ncube[0] * ncube[1] * ncube[2]
+        np.int64_t *head = <np.int64_t*> safe_malloc(nb_cubes*sizeof(np.int64_t))
+        np.int64_t *atom_indices = <np.int64_t*> safe_malloc(n_atoms*sizeof(np.int64_t))
+        np.int64_t[:, ::1] neighbor_map = <np.int64_t[:nb_cubes, :27]> safe_malloc(
+            nb_cubes * 27 * sizeof(np.int64_t)
         )
 
-    memset(<void*>head, -1, nb_cubes*sizeof(long))
-    memset(<void*>atom_indices, -1, n_atoms*sizeof(long))
+    memset(<void*>head, -1, nb_cubes*sizeof(np.int64_t))
+    memset(<void*>atom_indices, -1, n_atoms*sizeof(np.int64_t))
 
     get_cube_neighbors(ncube, neighbor_map)
     for i in range(n_atoms):
@@ -321,8 +321,8 @@ def find_points_in_spheres(
                     # compared to using vectors in cpp
                     if count >= malloc_chunk:
                         malloc_chunk += malloc_chunk  # double the size
-                        index_1 = <long*> realloc(index_1, malloc_chunk * sizeof(long))
-                        index_2 = <long*> realloc(index_2, malloc_chunk*sizeof(long))
+                        index_1 = <np.int64_t*> realloc(index_1, malloc_chunk * sizeof(np.int64_t))
+                        index_2 = <np.int64_t*> realloc(index_2, malloc_chunk*sizeof(np.int64_t))
                         offset_final = <double*> realloc(
                             offset_final, 3*malloc_chunk*sizeof(double)
                         )
@@ -355,14 +355,14 @@ def find_points_in_spheres(
         py_distances = np.array([], dtype=float)
     else:
         # resize to the actual size
-        index_1 = <long*>safe_realloc(index_1, count * sizeof(long))
-        index_2 = <long*>safe_realloc(index_2, count*sizeof(long))
+        index_1 = <np.int64_t*>safe_realloc(index_1, count * sizeof(np.int64_t))
+        index_2 = <np.int64_t*>safe_realloc(index_2, count*sizeof(np.int64_t))
         offset_final = <double*>safe_realloc(offset_final, 3*count*sizeof(double))
         distances = <double*>safe_realloc(distances, count*sizeof(double))
 
         # convert to python objects
-        py_index_1 = np.array(<long[:count]>index_1)
-        py_index_2 = np.array(<long[:count]>index_2)
+        py_index_1 = np.array(<np.int64_t[:count]>index_1)
+        py_index_2 = np.array(<np.int64_t[:count]>index_2)
         py_offsets = np.array(<double[:count, :3]>offset_final)
         py_distances = np.array(<double[:count]>distances)
 
@@ -391,39 +391,39 @@ def find_points_in_spheres(
     return py_index_1, py_index_2, py_offsets, py_distances
 
 
-cdef void get_cube_neighbors(long[3] ncube, long[:, ::1] neighbor_map):
+cdef void get_cube_neighbors(np.int64_t[3] ncube, np.int64_t[:, ::1] neighbor_map):
     """
     Get {cube_index: cube_neighbor_indices} map
     """
     cdef:
         int i, j, k
         int count = 0
-        long ncubes = ncube[0] * ncube[1] * ncube[2]
-        long[::1] counts = <long[:ncubes]> safe_malloc(ncubes * sizeof(long))
-        long[:, ::1] cube_indices_3d = <long[:ncubes, :3]> safe_malloc(
-            ncubes*3*sizeof(long)
+        np.int64_t ncubes = ncube[0] * ncube[1] * ncube[2]
+        np.int64_t[::1] counts = <np.int64_t[:ncubes]> safe_malloc(ncubes * sizeof(np.int64_t))
+        np.int64_t[:, ::1] cube_indices_3d = <np.int64_t[:ncubes, :3]> safe_malloc(
+            ncubes*3*sizeof(np.int64_t)
         )
-        long[::1] cube_indices_1d = <long[:ncubes]> safe_malloc(ncubes*sizeof(long))
+        np.int64_t[::1] cube_indices_1d = <np.int64_t[:ncubes]> safe_malloc(ncubes*sizeof(np.int64_t))
 
         # creating the memviews of c-arrays once substantially improves speed
         # but for some reason it makes the runtime scaling with the number of
         # atoms worse
-        long[1][3] index3_arr
-        long[:, ::1] index3 = index3_arr
-        long[1] index1_arr
-        long[::1] index1 = index1_arr
+        np.int64_t[1][3] index3_arr
+        np.int64_t[:, ::1] index3 = index3_arr
+        np.int64_t[1] index1_arr
+        np.int64_t[::1] index1 = index1_arr
 
         int n = 1
-        long ntotal = (2 * n + 1) * (2 * n + 1) * (2 * n + 1)
-        long[:, ::1] ovectors
-        long *ovectors_p = <long *> safe_malloc(ntotal * 3 * sizeof(long))
+        np.int64_t ntotal = (2 * n + 1) * (2 * n + 1) * (2 * n + 1)
+        np.int64_t[:, ::1] ovectors
+        np.int64_t *ovectors_p = <np.int64_t *> safe_malloc(ntotal * 3 * sizeof(np.int64_t))
         int n_ovectors = compute_offset_vectors(ovectors_p, n)
 
     # now resize to the actual size
-    ovectors_p = <long *> safe_realloc(ovectors_p, n_ovectors * 3 * sizeof(long))
-    ovectors = <long[:n_ovectors, :3]>ovectors_p
+    ovectors_p = <np.int64_t *> safe_realloc(ovectors_p, n_ovectors * 3 * sizeof(np.int64_t))
+    ovectors = <np.int64_t[:n_ovectors, :3]>ovectors_p
 
-    memset(<void*>&neighbor_map[0, 0], -1, neighbor_map.shape[0] * 27 * sizeof(long))
+    memset(<void*>&neighbor_map[0, 0], -1, neighbor_map.shape[0] * 27 * sizeof(np.int64_t))
 
     for i in range(ncubes):
         counts[i] = 0
@@ -461,7 +461,7 @@ cdef void get_cube_neighbors(long[3] ncube, long[:, ::1] neighbor_map):
     free(ovectors_p)
 
 
-cdef int compute_offset_vectors(long* ovectors, long n) nogil:
+cdef int compute_offset_vectors(np.int64_t* ovectors, np.int64_t n) nogil:
     cdef:
         int i, j, k, ind
         int count = 0
@@ -492,9 +492,9 @@ cdef int compute_offset_vectors(long* ovectors, long n) nogil:
 cdef double distance2(
         const double[:, ::1] m1,
         const double[:, ::1] m2,
-        long index1,
-        long index2,
-        long size
+        np.int64_t index1,
+        np.int64_t index2,
+        np.int64_t size
     ) nogil:
     """Faster way to compute the distance squared by not using slice but providing indices
     in each matrix
@@ -511,9 +511,9 @@ cdef double distance2(
 cdef void get_bounds(
         const double[:, ::1] frac_coords,
         const double[3] maxr,
-        const long[3] pbc,
-        long[3] max_bounds,
-        long[3] min_bounds
+        const np.int64_t[3] pbc,
+        np.int64_t[3] max_bounds,
+        np.int64_t[3] min_bounds
     ) nogil:
     """
     Given the fractional coordinates and the number of repeation needed in each
@@ -532,8 +532,8 @@ cdef void get_bounds(
 
     for i in range(3):
         if pbc[i]:
-            min_bounds[i] = <long>(floor(min_fcoords[i] - maxr[i] - 1e-8))
-            max_bounds[i] = <long>(ceil(max_fcoords[i] + maxr[i] + 1e-8))
+            min_bounds[i] = <np.int64_t>(floor(min_fcoords[i] - maxr[i] - 1e-8))
+            max_bounds[i] = <np.int64_t>(ceil(max_fcoords[i] + maxr[i] + 1e-8))
 
 cdef void get_frac_coords(
         const double[:, ::1] lattice,
@@ -697,18 +697,18 @@ cdef void max_and_min(
 cdef void compute_cube_index(
         const double[:, ::1] coords,
         const double[3] global_min,
-        double radius, long[:, ::1] return_indices
+        double radius, np.int64_t[:, ::1] return_indices
     ) nogil:
     cdef int i, j
     for i in range(coords.shape[0]):
         for j in range(coords.shape[1]):
-            return_indices[i, j] = <long>(
+            return_indices[i, j] = <np.int64_t>(
                 floor((coords[i, j] - global_min[j] + 1e-8) / radius)
             )
 
 
 cdef void three_to_one(
-        const long[:, ::1] label3d, long ny, long nz, long[::1] label1d
+        const np.int64_t[:, ::1] label3d, np.int64_t ny, np.int64_t nz, np.int64_t[::1] label1d
     ) nogil:
     """
     3D vector representation to 1D
@@ -740,7 +740,7 @@ cdef bint distance_vertices(
 
 cdef void offset_cube(
         const double[8][3] center,
-        long n, long m, long l,
+        np.int64_t n, np.int64_t m, np.int64_t l,
         const double[8][3] (&offsetted)
     ) nogil:
     cdef int i, j, k