diff --git a/pymatgen/analysis/chemenv/connectivity/connected_components.py b/pymatgen/analysis/chemenv/connectivity/connected_components.py
index bc8a4ac198a..7b933c4c56b 100644
--- a/pymatgen/analysis/chemenv/connectivity/connected_components.py
+++ b/pymatgen/analysis/chemenv/connectivity/connected_components.py
@@ -458,7 +458,7 @@ def compute_periodicity_all_simple_paths_algorithm(self):
                     path.append(test_node)
                     # TODO: there are some paths that appears twice for cycles, and there are some paths that should
                     # probably not be considered
-                    this_path_deltas = [np.zeros(3, np.int_)]
+                    this_path_deltas = [np.zeros(3, int)]
                     for (node1, node2) in [(node1, path[inode1 + 1]) for inode1, node1 in enumerate(path[:-1])]:
                         this_path_deltas_new = []
                         for key, edge_data in self._connected_subgraph[node1][node2].items():
@@ -499,7 +499,7 @@ def compute_periodicity_cycle_basis(self):
         for cyc in cycles:
             mycyc = list(cyc)
             mycyc.append(cyc[0])
-            this_cycle_deltas = [np.zeros(3, np.int_)]
+            this_cycle_deltas = [np.zeros(3, int)]
             for (node1, node2) in [(node1, mycyc[inode1 + 1]) for inode1, node1 in enumerate(mycyc[:-1])]:
                 this_cycle_deltas_new = []
                 for key, edge_data in self._connected_subgraph[node1][node2].items():
@@ -778,7 +778,7 @@ def elastic_centered_graph(self, start_node=None):
                     logging.debug("          Edge outside the cell ... getting neighbor back inside")
                     if (0, 0, 0) in ddeltas:
                         ddeltas.remove((0, 0, 0))
-                    myddelta = np.array(ddeltas[0], np.int_)
+                    myddelta = np.array(ddeltas[0], int)
                     node_neighbor_edges = centered_connected_subgraph.edges(
                         nbunch=[node_neighbor], data=True, keys=True
                     )
@@ -797,11 +797,11 @@ def elastic_centered_graph(self, start_node=None):
                         ):
                             if edata["start"] == node_neighbor.isite and edata["end"] != node_neighbor.isite:
                                 centered_connected_subgraph[n1][n2][key]["delta"] = tuple(
-                                    np.array(edata["delta"], np.int_) + myddelta
+                                    np.array(edata["delta"], int) + myddelta
                                 )
                             elif edata["end"] == node_neighbor.isite:
                                 centered_connected_subgraph[n1][n2][key]["delta"] = tuple(
-                                    np.array(edata["delta"], np.int_) - myddelta
+                                    np.array(edata["delta"], int) - myddelta
                                 )
                             else:
                                 raise ValueError("DUHH")
diff --git a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
index 73283a2540a..474da317263 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
@@ -1700,7 +1700,7 @@ def coordination_geometry_symmetry_measures_separation_plane_optim(
                         continue
                     if sep not in nb_set.separations:
                         nb_set.separations[sep] = {}
-                    mysep = [np.array(ss, dtype=np.int8) for ss in separation]
+                    mysep = [np.array(ss, dtype=int) for ss in separation]
                     nb_set.separations[sep][separation] = (plane, mysep)
                     if sep == separation_plane_algo.separation:
                         new_seps.append(mysep)
diff --git a/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py b/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
index 966afc460df..10f5cc9ef64 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
@@ -1648,7 +1648,7 @@ def from_structure_environments(cls, strategy, structure_environments, valences=
                             raise ValueError(
                                 "Weird, differences between one site in a periodic image cell is not integer ..."
                             )
-                        nb_image_cell = np.array(rounddiff, np.int_)
+                        nb_image_cell = np.array(rounddiff, int)
                         nb_allnbs_sites_index = len(_all_nbs_sites)
                         _all_nbs_sites.append(
                             {
@@ -2129,13 +2129,13 @@ def from_dict(cls, d):
         for nb_site in d["all_nbs_sites"]:
             site = dec.process_decoded(nb_site["site"])
             if "image_cell" in nb_site:
-                image_cell = np.array(nb_site["image_cell"], np.int_)
+                image_cell = np.array(nb_site["image_cell"], int)
             else:
                 diff = site.frac_coords - structure[nb_site["index"]].frac_coords
                 rounddiff = np.round(diff)
                 if not np.allclose(diff, rounddiff):
                     raise ValueError("Weird, differences between one site in a periodic image cell is not integer ...")
-                image_cell = np.array(rounddiff, np.int_)
+                image_cell = np.array(rounddiff, int)
             all_nbs_sites.append({"site": site, "index": nb_site["index"], "image_cell": image_cell})
         neighbors_sets = [
             [
diff --git a/pymatgen/analysis/chemenv/utils/graph_utils.py b/pymatgen/analysis/chemenv/utils/graph_utils.py
index 65d7fc1abcd..f58f5d84b54 100644
--- a/pymatgen/analysis/chemenv/utils/graph_utils.py
+++ b/pymatgen/analysis/chemenv/utils/graph_utils.py
@@ -21,9 +21,9 @@ def get_delta(node1, node2, edge_data):
     :return:
     """
     if node1.isite == edge_data["start"] and node2.isite == edge_data["end"]:
-        return np.array(edge_data["delta"], dtype=np.int_)
+        return np.array(edge_data["delta"], dtype=int)
     if node2.isite == edge_data["start"] and node1.isite == edge_data["end"]:
-        return -np.array(edge_data["delta"], dtype=np.int_)
+        return -np.array(edge_data["delta"], dtype=int)
     raise ValueError("Trying to find a delta between two nodes with an edge that seems not to link these nodes.")
 
 
@@ -490,7 +490,7 @@ def get_all_elementary_cycles(graph):
         all_edges_dict[(n2, n1)] = nedges
         index2edge.append((n1, n2))
         nedges += 1
-    cycles_matrix = np.zeros(shape=(len(cycle_basis), nedges), dtype=np.bool)
+    cycles_matrix = np.zeros(shape=(len(cycle_basis), nedges), dtype=bool)
     for icycle, cycle in enumerate(cycle_basis):
         for in1, n1 in enumerate(cycle):
             n2 = cycle[(in1 + 1) % len(cycle)]
@@ -502,7 +502,7 @@ def get_all_elementary_cycles(graph):
 
     for ncycles in range(1, len(cycle_basis) + 1):
         for cycles_combination in itertools.combinations(cycles_matrix, ncycles):
-            edges_counts = np.array(np.mod(np.sum(cycles_combination, axis=0), 2), dtype=np.bool)
+            edges_counts = np.array(np.mod(np.sum(cycles_combination, axis=0), 2), dtype=bool)
             myedges = [edge for iedge, edge in enumerate(index2edge) if edges_counts[iedge]]
             # print(myedges)
             try:
diff --git a/pymatgen/analysis/defects/dilute_solution_model.py b/pymatgen/analysis/defects/dilute_solution_model.py
index 935fd2922db..a88540e76f4 100644
--- a/pymatgen/analysis/defects/dilute_solution_model.py
+++ b/pymatgen/analysis/defects/dilute_solution_model.py
@@ -139,7 +139,7 @@ def dilute_solution_model(structure, e0, vac_defs, antisite_defs, T, trial_chem_
                 dC[i,k,k] = 0 due to no effect on ith type atom
                 dC[i,j,k] = 0 if i!=j!=k
     """
-    dC = np.zeros((n, n, n), dtype=np.int_)
+    dC = np.zeros((n, n, n), dtype=int)
     for i in range(n):
         for j in range(n):
             for k in range(n):
@@ -852,7 +852,7 @@ def solute_site_preference_finder(
                 dC[i,k,k] = 0 due to no effect on ith type atom
                 dC[i,j,k] = 0 if i!=j!=k
     """
-    dC = np.zeros((n + 1, n + 1, n), dtype=np.int_)
+    dC = np.zeros((n + 1, n + 1, n), dtype=int)
     for i in range(n):
         for j in range(n):
             for k in range(n):
diff --git a/pymatgen/analysis/gb/grain.py b/pymatgen/analysis/gb/grain.py
index 660d02ca864..489578fd13c 100644
--- a/pymatgen/analysis/gb/grain.py
+++ b/pymatgen/analysis/gb/grain.py
@@ -1256,7 +1256,7 @@ def get_trans_mat(
         r_matrix = np.dot(np.dot(np.linalg.inv(trans_cry.T), r_matrix), trans_cry.T)
         # set one vector of the basis to the rotation axis direction, and
         # obtain the corresponding transform matrix
-        eye = np.eye(3, dtype=np.int_)
+        eye = np.eye(3, dtype=int)
         for h in range(3):
             if abs(r_axis[h]) != 0:
                 eye[h] = np.array(r_axis)
@@ -2107,7 +2107,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals
         # quickly generate a supercell, normal is not work in this way
         if quick_gen:
             scale_factor = []
-            eye = np.eye(3, dtype=np.int_)
+            eye = np.eye(3, dtype=int)
             for i, j in enumerate(miller):
                 if j == 0:
                     scale_factor.append(eye[i])
diff --git a/pymatgen/analysis/graphs.py b/pymatgen/analysis/graphs.py
index dae81c8daab..cd1aa484122 100644
--- a/pymatgen/analysis/graphs.py
+++ b/pymatgen/analysis/graphs.py
@@ -1135,9 +1135,9 @@ def __mul__(self, scaling_matrix):
         # easier to extend to a general 3x3 scaling matrix.
 
         # code adapted from Structure.__mul__
-        scale_matrix = np.array(scaling_matrix, np.int16)
+        scale_matrix = np.array(scaling_matrix, int)
         if scale_matrix.shape != (3, 3):
-            scale_matrix = np.array(scale_matrix * np.eye(3), np.int16)
+            scale_matrix = np.array(scale_matrix * np.eye(3), int)
         else:
             # TODO: test __mul__ with full 3x3 scaling matrices
             raise NotImplementedError("Not tested with 3x3 scaling matrices yet.")
diff --git a/pymatgen/analysis/local_env.py b/pymatgen/analysis/local_env.py
index 8b38cfba95c..2554b035713 100644
--- a/pymatgen/analysis/local_env.py
+++ b/pymatgen/analysis/local_env.py
@@ -778,7 +778,7 @@ def get_all_voronoi_polyhedra(self, structure):
             indices.extend([(x[2],) + x[3] for x in neighs])
 
         # Get the non-duplicates (using the site indices for numerical stability)
-        indices = np.array(indices, dtype=np.int_)
+        indices = np.array(indices, dtype=int)
         indices, uniq_inds = np.unique(indices, return_index=True, axis=0)
         sites = [sites[i] for i in uniq_inds]
 
diff --git a/pymatgen/analysis/structure_matcher.py b/pymatgen/analysis/structure_matcher.py
index a3db9a44952..0d3aca1713d 100644
--- a/pymatgen/analysis/structure_matcher.py
+++ b/pymatgen/analysis/structure_matcher.py
@@ -553,7 +553,7 @@ def _get_mask(self, struct1, struct2, fu, s1_supercell):
         Returns:
         mask, struct1 translation indices, struct2 translation index
         """
-        mask = np.zeros((len(struct2), len(struct1), fu), dtype=np.bool)
+        mask = np.zeros((len(struct2), len(struct1), fu), dtype=bool)
 
         inner = []
         for sp2, i in itertools.groupby(enumerate(struct2.species_and_occu), key=lambda x: x[1]):
@@ -580,7 +580,7 @@ def _get_mask(self, struct1, struct2, fu, s1_supercell):
         if s1_supercell:
             # remove the symmetrically equivalent s1 indices
             inds = inds[::fu]
-        return np.array(mask, dtype=np.int_), inds, i
+        return np.array(mask, dtype=int), inds, i
 
     def fit(self, struct1, struct2, symmetric=False):
         """
diff --git a/pymatgen/core/lattice.py b/pymatgen/core/lattice.py
index a349d35c25b..273ca9d01d0 100644
--- a/pymatgen/core/lattice.py
+++ b/pymatgen/core/lattice.py
@@ -186,7 +186,7 @@ class in `pymatgen.analysis.ferroelectricity.polarization`.
         Returns:
             Lattice coordinates.
         """
-        return self.lengths * self.get_fractional_coords(cart_coords)
+        return self.lengths * self.get_fractional_coords(cart_coords)  # type: ignore
 
     def d_hkl(self, miller_index: ArrayLike) -> float:
         """
@@ -1028,7 +1028,7 @@ def find_all_mappings(
         # this can't be broadcast because they're different lengths
         inds = [np.logical_and(dist / l < 1 + ltol, dist / l > 1 / (1 + ltol)) for l in lengths]  # type: ignore
         c_a, c_b, c_c = (cart[i] for i in inds)
-        f_a, f_b, f_c = (frac[i] for i in inds)
+        f_a, f_b, f_c = (frac[i] for i in inds)  # type: ignore
         l_a, l_b, l_c = (np.sum(c ** 2, axis=-1) ** 0.5 for c in (c_a, c_b, c_c))
 
         def get_angles(v1, v2, l1, l2):
@@ -1045,7 +1045,7 @@ def get_angles(v1, v2, l1, l2):
         for i, all_j in enumerate(gammab):
             inds = np.logical_and(all_j[:, None], np.logical_and(alphab, betab[i][None, :]))
             for j, k in np.argwhere(inds):
-                scale_m = np.array((f_a[i], f_b[j], f_c[k]), dtype=np.int_)  # type: ignore
+                scale_m = np.array((f_a[i], f_b[j], f_c[k]), dtype=int)  # type: ignore
                 if abs(np.linalg.det(scale_m)) < 1e-8:
                     continue
 
@@ -1123,7 +1123,7 @@ def _calculate_lll(self, delta: float = 0.75) -> Tuple[np.ndarray, np.ndarray]:
         """
         # Transpose the lattice matrix first so that basis vectors are columns.
         # Makes life easier.
-        # pylint: disable=E1136,E1137
+        # pylint: disable=E1136,E1137,E1126
         a = self._matrix.copy().T
 
         b = np.zeros((3, 3))  # Vectors after the Gram-Schmidt process
@@ -1179,7 +1179,8 @@ def _calculate_lll(self, delta: float = 0.75) -> Tuple[np.ndarray, np.ndarray]:
                 else:
                     # We have to do p/q, so do lstsq(q.T, p.T).T instead.
                     p = dot(a[:, k:3].T, b[:, (k - 2) : k])
-                    q = np.diag(m[(k - 2) : k])
+                    q = np.diag(m[(k - 2) : k])  # type: ignore
+                    # pylint: disable=E1101
                     result = np.linalg.lstsq(q.T, p.T, rcond=None)[0].T  # type: ignore
                     u[k:3, (k - 2) : k] = result
 
@@ -1436,7 +1437,7 @@ def norm(self, coords: ArrayLike, frac_coords: bool = True) -> float:
         Returns:
             one-dimensional `numpy` array.
         """
-        return np.sqrt(self.dot(coords, coords, frac_coords=frac_coords))
+        return np.sqrt(self.dot(coords, coords, frac_coords=frac_coords))  # type: ignore
 
     def get_points_in_sphere(
         self,
@@ -1630,12 +1631,12 @@ def get_points_in_sphere_old(
         # Generate all possible images that could be within `r` of `center`
         mins = np.floor(pcoords - nmax)
         maxes = np.ceil(pcoords + nmax)
-        arange = np.arange(start=mins[0], stop=maxes[0], dtype=np.int_)
-        brange = np.arange(start=mins[1], stop=maxes[1], dtype=np.int_)
-        crange = np.arange(start=mins[2], stop=maxes[2], dtype=np.int_)
-        arange = arange[:, None] * np.array([1, 0, 0], dtype=np.int_)[None, :]
-        brange = brange[:, None] * np.array([0, 1, 0], dtype=np.int_)[None, :]
-        crange = crange[:, None] * np.array([0, 0, 1], dtype=np.int_)[None, :]
+        arange = np.arange(start=mins[0], stop=maxes[0], dtype=int)
+        brange = np.arange(start=mins[1], stop=maxes[1], dtype=int)
+        crange = np.arange(start=mins[2], stop=maxes[2], dtype=int)
+        arange = arange[:, None] * np.array([1, 0, 0], dtype=int)[None, :]
+        brange = brange[:, None] * np.array([0, 1, 0], dtype=int)[None, :]
+        crange = crange[:, None] * np.array([0, 0, 1], dtype=int)[None, :]
         images = arange[:, None, None] + brange[None, :, None] + crange[None, None, :]
 
         # Generate the coordinates of all atoms within these images
@@ -1644,7 +1645,7 @@ def get_points_in_sphere_old(
         # Determine distance from `center`
         cart_coords = self.get_cartesian_coords(fcoords)
         cart_images = self.get_cartesian_coords(images)
-        coords = cart_coords[:, None, None, None, :] + cart_images[None, :, :, :, :]
+        coords = cart_coords[:, None, None, None, :] + cart_images[None, :, :, :, :]  # pylint: disable=E1126
         coords -= center[None, None, None, None, :]
         coords **= 2
         d_2 = np.sum(coords, axis=4)
@@ -1747,13 +1748,13 @@ def get_distance_and_image(
         """
         if jimage is None:
             v, d2 = pbc_shortest_vectors(self, frac_coords1, frac_coords2, return_d2=True)
-            fc = self.get_fractional_coords(v[0][0]) + frac_coords1 - frac_coords2
-            fc = np.array(np.round(fc), dtype=np.int_)
+            fc = self.get_fractional_coords(v[0][0]) + frac_coords1 - frac_coords2  # type: ignore
+            fc = np.array(np.round(fc), dtype=int)
             return np.sqrt(d2[0, 0]), fc
 
         jimage = np.array(jimage)
-        mapped_vec = self.get_cartesian_coords(jimage + frac_coords2 - frac_coords1)
-        return np.linalg.norm(mapped_vec), jimage
+        mapped_vec = self.get_cartesian_coords(jimage + frac_coords2 - frac_coords1)  # type: ignore
+        return np.linalg.norm(mapped_vec), jimage  # type: ignore
 
     def get_miller_index_from_coords(
         self,
@@ -1846,14 +1847,14 @@ def get_integer_index(miller_index: Sequence[float], round_dp: int = 4, verbose:
     # deal with the case we have nice fractions
     md = [Fraction(n).limit_denominator(12).denominator for n in mi]
     mi *= reduce(lambda x, y: x * y, md)
-    int_miller_index = np.int_(np.round(mi, 1))
+    int_miller_index = np.round(mi, 1).astype(int)
     mi /= np.abs(reduce(math.gcd, int_miller_index))
 
     # round to a reasonable precision
     mi = np.array([round(h, round_dp) for h in mi])
 
     # need to recalculate this after rounding as values may have changed
-    int_miller_index = np.int_(np.round(mi, 1))
+    int_miller_index = np.round(mi, 1).astype(int)
     if np.any(np.abs(mi - int_miller_index) > 1e-6) and verbose:
         warnings.warn("Non-integer encountered in Miller index")
     else:
@@ -1902,7 +1903,7 @@ def get_points_in_spheres(
     """
     if isinstance(pbc, bool):
         pbc = [pbc] * 3
-    pbc = np.array(pbc, dtype=np.bool_)  # type: ignore
+    pbc = np.array(pbc, dtype=bool)  # type: ignore
     if return_fcoords and lattice is None:
         raise ValueError("Lattice needs to be supplied to compute fractional coordinates")
     center_coords_min = np.min(center_coords, axis=0)
@@ -1959,7 +1960,7 @@ def get_points_in_spheres(
     else:
         valid_coords = all_coords  # type: ignore
         valid_images = [[0, 0, 0]] * len(valid_coords)
-        valid_indices = np.arange(len(valid_coords))
+        valid_indices = np.arange(len(valid_coords))  # type: ignore
 
     # Divide the valid 3D space into cubes and compute the cube ids
     all_cube_index = _compute_cube_index(valid_coords, global_min, r)  # type: ignore
diff --git a/pymatgen/core/structure.py b/pymatgen/core/structure.py
index fcafb6d45d4..8972327cd31 100644
--- a/pymatgen/core/structure.py
+++ b/pymatgen/core/structure.py
@@ -1096,9 +1096,9 @@ def __mul__(self, scaling_matrix: Union[int, Sequence[int], Sequence[Sequence[in
             you prefer a subclass to return its own type, you need to override
             this method in the subclass.
         """
-        scale_matrix = np.array(scaling_matrix, np.int16)
+        scale_matrix = np.array(scaling_matrix, int)
         if scale_matrix.shape != (3, 3):
-            scale_matrix = np.array(scale_matrix * np.eye(3), np.int16)
+            scale_matrix = np.array(scale_matrix * np.eye(3), int)
         new_lattice = Lattice(np.dot(scale_matrix, self._lattice.matrix))
 
         f_lat = lattice_points_in_supercell(scale_matrix)
diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py
index ec9e38ea671..bb2940f2811 100644
--- a/pymatgen/core/surface.py
+++ b/pymatgen/core/surface.py
@@ -844,7 +844,7 @@ def __init__(
 
         slab_scale_factor = []
         non_orth_ind = []
-        eye = np.eye(3, dtype=np.int_)
+        eye = np.eye(3, dtype=int)
         for i, j in enumerate(miller_index):
             if j == 0:
                 # Lattice vector is perpendicular to surface normal, i.e.,
diff --git a/pymatgen/core/tests/test_units.py b/pymatgen/core/tests/test_units.py
index 67a213f6fb2..435a40c485d 100644
--- a/pymatgen/core/tests/test_units.py
+++ b/pymatgen/core/tests/test_units.py
@@ -184,7 +184,7 @@ def test_time(self):
         Check whether EnergyArray and FloatWithUnit have same behavior.
         """
         # here there's a minor difference because we have a ndarray with
-        # dtype=np.int_.
+        # dtype=int.
         a = TimeArray(20, "h")
         self.assertAlmostEqual(a.to("s"), 3600 * 20)
         # Test left and right multiplication.
diff --git a/pymatgen/io/abinit/abiobjects.py b/pymatgen/io/abinit/abiobjects.py
index c05707316d8..d8b90dc6200 100644
--- a/pymatgen/io/abinit/abiobjects.py
+++ b/pymatgen/io/abinit/abiobjects.py
@@ -159,7 +159,7 @@ def structure_from_abivars(cls=None, *args, **kwargs):
         raise ValueError(f"len(typat) != len(coords):\ntypat: {typat}\ncoords: {coords}")
 
     # Note conversion to int and Fortran --> C indexing
-    typat = np.array(typat, dtype=np.int_)
+    typat = np.array(typat, dtype=int)
     species = [znucl_type[typ - 1] for typ in typat]
 
     return cls(
@@ -244,7 +244,7 @@ def structure_to_abivars(structure, enforce_znucl=None, enforce_typat=None, **kw
 
         # [ntypat] list
         znucl_type = [specie.number for specie in types_of_specie]
-        typat = np.zeros(natom, np.int_)
+        typat = np.zeros(natom, int)
         for atm_idx, site in enumerate(structure):
             typat[atm_idx] = types_of_specie.index(site.specie) + 1
     else:
diff --git a/pymatgen/io/vasp/outputs.py b/pymatgen/io/vasp/outputs.py
index d0934e7bae6..b9d3357eb22 100644
--- a/pymatgen/io/vasp/outputs.py
+++ b/pymatgen/io/vasp/outputs.py
@@ -4812,7 +4812,7 @@ def __init__(self, filename="WAVECAR", verbose=False, precision="normal", vasp_t
         self._C = 0.262465831
         with open(self.filename, "rb") as f:
             # read the header information
-            recl, spin, rtag = np.fromfile(f, dtype=np.float64, count=3).astype(np.int_)
+            recl, spin, rtag = np.fromfile(f, dtype=np.float64, count=3).astype(int)
             if verbose:
                 print(f"recl={recl}, spin={spin}, rtag={rtag}")
             recl8 = int(recl / 8)
@@ -4831,7 +4831,7 @@ def __init__(self, filename="WAVECAR", verbose=False, precision="normal", vasp_t
             np.fromfile(f, dtype=np.float64, count=(recl8 - 3))
 
             # extract kpoint, bands, energy, and lattice information
-            self.nk, self.nb, self.encut = np.fromfile(f, dtype=np.float64, count=3).astype(np.int_)
+            self.nk, self.nb, self.encut = np.fromfile(f, dtype=np.float64, count=3).astype(int)
             self.a = np.fromfile(f, dtype=np.float64, count=9).reshape((3, 3))
             self.efermi = np.fromfile(f, dtype=np.float64, count=1)[0]
             if verbose:
@@ -5007,7 +5007,7 @@ def _generate_nbmax(self) -> None:
         nbmaxC[2] /= np.abs(np.sin(phi23))
         nbmaxC += 1
 
-        self._nbmax = np.max([nbmaxA, nbmaxB, nbmaxC], axis=0).astype(np.int_)
+        self._nbmax = np.max([nbmaxA, nbmaxB, nbmaxC], axis=0).astype(int)
 
     def _generate_G_points(self, kpoint: np.ndarray, gamma: bool = False) -> Tuple[List, List, List]:
         """
@@ -5131,7 +5131,7 @@ def fft_mesh(self, kpoint: int, band: int, spin: int = 0, spinor: int = 0, shift
 
         mesh = np.zeros(tuple(self.ng), dtype=np.complex_)
         for gp, coeff in zip(self.Gpoints[kpoint], tcoeffs):
-            t = tuple(gp.astype(np.int_) + (self.ng / 2).astype(np.int_))
+            t = tuple(gp.astype(int) + (self.ng / 2).astype(int))
             mesh[t] = coeff
 
         if shift:
@@ -5523,9 +5523,9 @@ def readData(dtype):
                 np.array of given dtype."""
                 data = b""
                 while True:
-                    prefix = np.fromfile(fp, dtype=np.int32, count=1)[0]
+                    prefix = np.fromfile(fp, dtype=int, count=1)[0]
                     data += fp.read(abs(prefix))
-                    suffix = np.fromfile(fp, dtype=np.int32, count=1)[0]
+                    suffix = np.fromfile(fp, dtype=int, count=1)[0]
                     if abs(prefix) - abs(suffix):
                         raise RuntimeError(
                             "Read wrong amount of bytes.\n"
@@ -5535,7 +5535,7 @@ def readData(dtype):
                         break
                 return np.frombuffer(data, dtype=dtype)
 
-            nbands, nelect, nk, ispin = readData(np.int32)
+            nbands, nelect, nk, ispin = readData(int)
             _ = readData(np.float_)  # nodes_in_dielectric_function
             _ = readData(np.float_)  # wplasmon
             if gamma_only:
diff --git a/pymatgen/io/vasp/tests/test_outputs.py b/pymatgen/io/vasp/tests/test_outputs.py
index 6114fc54b01..e95b4ffb5f1 100644
--- a/pymatgen/io/vasp/tests/test_outputs.py
+++ b/pymatgen/io/vasp/tests/test_outputs.py
@@ -1865,7 +1865,7 @@ def test_fft_mesh_basic(self):
         mesh = self.w.fft_mesh(0, 5)
         ind = np.argmax(np.abs(mesh))
         self.assertEqual(np.unravel_index(ind, mesh.shape), (14, 1, 1))
-        self.assertEqual(mesh[tuple((self.w.ng / 2).astype(np.int_))], 0j)
+        self.assertEqual(mesh[tuple((self.w.ng / 2).astype(int))], 0j)
         mesh = self.w.fft_mesh(0, 5, shift=False)
         ind = np.argmax(np.abs(mesh))
         self.assertEqual(np.unravel_index(ind, mesh.shape), (6, 8, 8))
diff --git a/pymatgen/io/wannier90.py b/pymatgen/io/wannier90.py
index 97b07e1a809..134d6548fc5 100644
--- a/pymatgen/io/wannier90.py
+++ b/pymatgen/io/wannier90.py
@@ -149,7 +149,7 @@ def write_file(self, filename: str) -> None:
                 YY is 1 or 2 for the spin index or NC if noncollinear
         """
         with FortranFile(filename, "w") as f:
-            f.write_record(np.array([*self.ng, self.ik, self.nbnd], dtype=np.int32))
+            f.write_record(np.array([*self.ng, self.ik, self.nbnd], dtype=int))
             for ib in range(self.nbnd):
                 if self.is_noncollinear:
                     f.write_record(self.data[ib, 0].flatten("F"))
diff --git a/pymatgen/optimization/linear_assignment_numpy.py b/pymatgen/optimization/linear_assignment_numpy.py
index 9b02f35d139..56a3b759793 100644
--- a/pymatgen/optimization/linear_assignment_numpy.py
+++ b/pymatgen/optimization/linear_assignment_numpy.py
@@ -71,7 +71,7 @@ def __init__(self, costs, epsilon=1e-6):
             self.c[: self.nx] = self.orig_c
 
         # initialize solution vectors
-        self._x = np.zeros(self.n, dtype=np.int_) - 1
+        self._x = np.zeros(self.n, dtype=int) - 1
         self._y = self._x.copy()
 
         # if column reduction doesn't find a solution, augment with shortest
@@ -195,7 +195,7 @@ def _build_tree(self):
 
         # compute distances
         self._d = self.c[istar] - self._v
-        _pred = np.zeros(self.n, dtype=np.int_) + istar
+        _pred = np.zeros(self.n, dtype=int) + istar
 
         # initialize sets
         # READY: set of nodes visited and in the path (whose price gets
@@ -203,9 +203,9 @@ def _build_tree(self):
         # SCAN: set of nodes at the bottom of the tree, which we need to
         # look at
         # T0DO: unvisited nodes
-        _ready = np.zeros(self.n, dtype=np.bool)
-        _scan = np.zeros(self.n, dtype=np.bool)
-        _todo = np.zeros(self.n, dtype=np.bool) + True
+        _ready = np.zeros(self.n, dtype=bool)
+        _scan = np.zeros(self.n, dtype=bool)
+        _todo = np.zeros(self.n, dtype=bool) + True
 
         while True:
             # populate scan with minimum reduced distances
diff --git a/pymatgen/optimization/tests/test_linear_assignment.py b/pymatgen/optimization/tests/test_linear_assignment.py
index c674281efed..4a73c792427 100644
--- a/pymatgen/optimization/tests/test_linear_assignment.py
+++ b/pymatgen/optimization/tests/test_linear_assignment.py
@@ -259,7 +259,7 @@ def test_small_range(self):
         self.assertAlmostEqual(LinearAssignment(x).min_cost, 48)
 
     def test_boolean_inputs(self):
-        w = np.ones((135, 135), dtype=np.bool)
+        w = np.ones((135, 135), dtype=bool)
         np.fill_diagonal(w, False)
         la = LinearAssignment(w)
         # if the input doesn't get converted to a float, the masking
diff --git a/pymatgen/util/coord.py b/pymatgen/util/coord.py
index 5f34581cd12..de35e5f4c31 100644
--- a/pymatgen/util/coord.py
+++ b/pymatgen/util/coord.py
@@ -258,9 +258,9 @@ def is_coord_subset_pbc(subset, superset, atol=1e-8, mask=None):
     c1 = np.array(subset, dtype=np.float64)
     c2 = np.array(superset, dtype=np.float64)
     if mask is not None:
-        m = np.array(mask, dtype=np.int_)
+        m = np.array(mask, dtype=int)
     else:
-        m = np.zeros((len(subset), len(superset)), dtype=np.int_)
+        m = np.zeros((len(subset), len(superset)), dtype=int)
     atol = np.zeros(3, dtype=np.float64) + atol
     return cuc.is_coord_subset_pbc(c1, c2, atol, m)