diff --git a/docs_rst/conf.py b/docs_rst/conf.py
index 9fd2d71..dce520f 100644
--- a/docs_rst/conf.py
+++ b/docs_rst/conf.py
@@ -181,7 +181,7 @@
 
 # -- Options for LaTeX output ------------------------------------------------
 
-latex_elements = {
+latex_elements: dict = {
     # The paper size ('letterpaper' or 'a4paper').
     #'papersize': 'letterpaper',
     # The font size ('10pt', '11pt' or '12pt').
diff --git a/pymatgen/analysis/diffusion/aimd/clustering.py b/pymatgen/analysis/diffusion/aimd/clustering.py
index 80935a4..1b8e497 100644
--- a/pymatgen/analysis/diffusion/aimd/clustering.py
+++ b/pymatgen/analysis/diffusion/aimd/clustering.py
@@ -24,14 +24,14 @@
 class Kmeans:
     """Simple kmeans clustering."""
 
-    def __init__(self, max_iterations: int = 1000):
+    def __init__(self, max_iterations: int = 1000) -> None:
         """
         Args:
             max_iterations (int): Maximum number of iterations to run KMeans algo.
         """
         self.max_iterations = max_iterations
 
-    def cluster(self, points, k, initial_centroids=None):
+    def cluster(self, points: np.ndarray, k: int, initial_centroids: np.ndarray = None) -> tuple:
         """
         Args:
             points (ndarray): Data points as a mxn ndarray, where m is the
@@ -67,7 +67,7 @@ def cluster(self, points, k, initial_centroids=None):
         return centroids, labels, ss
 
     @staticmethod
-    def get_labels(points, centroids):
+    def get_labels(points: np.ndarray, centroids: np.ndarray) -> tuple:
         """
         For each element in the dataset, chose the closest centroid.
         Make that centroid the element's label.
@@ -81,7 +81,7 @@ def get_labels(points, centroids):
         return np.where(dists == min_dists[:, None])[1], np.sum(min_dists**2)
 
     @staticmethod
-    def get_centroids(points, labels, k, centroids):
+    def get_centroids(points: np.ndarray, labels: np.ndarray, k: int, centroids: np.ndarray) -> np.ndarray:
         """
         Each centroid is the geometric mean of the points that
         have that centroid's label. Important: If a centroid is empty (no
@@ -94,16 +94,16 @@ def get_centroids(points, labels, k, centroids):
             centroids: List of centroids
         """
         labels = np.array(labels)
-        centroids = []
+        _centroids = []
         for i in range(k):
             ind = np.where(labels == i)[0]
             if len(ind) > 0:
-                centroids.append(np.average(points[ind, :], axis=0))
+                _centroids.append(np.average(points[ind, :], axis=0))
             else:
-                centroids.append(get_random_centroid(points))
-        return np.array(centroids)
+                _centroids.append(get_random_centroid(points))
+        return np.array(_centroids)
 
-    def should_stop(self, old_centroids, centroids, iterations):
+    def should_stop(self, old_centroids: np.ndarray | None, centroids: np.ndarray, iterations: int) -> bool:
         """
         Check for stopping conditions.
 
@@ -127,7 +127,7 @@ class KmeansPBC(Kmeans):
     fractional coordinates.
     """
 
-    def __init__(self, lattice, max_iterations=1000):
+    def __init__(self, lattice: np.ndarray, max_iterations: int = 1000) -> None:
         """
         Args:
             lattice: Lattice
@@ -136,7 +136,7 @@ def __init__(self, lattice, max_iterations=1000):
         self.lattice = lattice
         self.max_iterations = max_iterations
 
-    def get_labels(self, points, centroids):
+    def get_labels(self, points, centroids):  # noqa: ANN001,ANN201
         """
         For each element in the dataset, chose the closest centroid.
         Make that centroid the element's label.
@@ -149,7 +149,7 @@ def get_labels(self, points, centroids):
         min_dists = np.min(dists, axis=1)
         return np.where(dists == min_dists[:, None])[1], np.sum(min_dists**2)
 
-    def get_centroids(self, points, labels, k, centroids):
+    def get_centroids(self, points, labels, k, centroids):  # noqa: ANN001,ANN201
         """
         Each centroid is the geometric mean of the points that
         have that centroid's label. Important: If a centroid is empty (no
@@ -179,7 +179,7 @@ def get_centroids(self, points, labels, k, centroids):
             new_centroids.append(c)
         return np.array(new_centroids)
 
-    def should_stop(self, old_centroids, centroids, iterations):
+    def should_stop(self, old_centroids: np.ndarray | None, centroids: np.ndarray, iterations: int) -> bool:
         """
         Check for stopping conditions.
 
@@ -196,7 +196,7 @@ def should_stop(self, old_centroids, centroids, iterations):
         return all(np.allclose(pbc_diff(c1, c2), [0, 0, 0]) for c1, c2 in zip(old_centroids, centroids))
 
 
-def get_random_centroid(points):
+def get_random_centroid(points: np.ndarray) -> np.ndarray:
     """
     Generate a random centroid based on points.
 
@@ -209,7 +209,7 @@ def get_random_centroid(points):
     return np.array([random.uniform(mind[i], maxd[i]) for i in range(n)])
 
 
-def get_random_centroids(points, k):
+def get_random_centroids(points: np.ndarray, k: int) -> np.ndarray:
     """
     Generate k random centroids based on points.
 
diff --git a/pymatgen/analysis/diffusion/aimd/pathway.py b/pymatgen/analysis/diffusion/aimd/pathway.py
index 9ddc21a..380ccae 100644
--- a/pymatgen/analysis/diffusion/aimd/pathway.py
+++ b/pymatgen/analysis/diffusion/aimd/pathway.py
@@ -4,11 +4,19 @@
 
 import itertools
 from collections import Counter
+from typing import TYPE_CHECKING
 
 import numpy as np
 from scipy.cluster.hierarchy import fcluster, linkage
 from scipy.spatial.distance import squareform
 
+if TYPE_CHECKING:
+    from collections.abc import Sequence
+
+    from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
+    from pymatgen.core.structure import Structure
+    from pymatgen.util.typing import PathLike, SpeciesLike
+
 
 class ProbabilityDensityAnalysis:
     r"""
@@ -25,7 +33,13 @@ class ProbabilityDensityAnalysis:
     Conductor". Chem. Mater. (2015), 27, pp 8318-8325.
     """
 
-    def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")):
+    def __init__(
+        self,
+        structure: Structure,
+        trajectories: np.ndarray,
+        interval: float = 0.5,
+        species: Sequence[SpeciesLike] = ("Li", "Na"),
+    ) -> None:
         """
         Initialization.
 
@@ -66,7 +80,7 @@ def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")):
         grid = agrid[:, None, None] + bgrid[None, :, None] + cgrid[None, None, :]
 
         # Calculate time-averaged probability density function distribution Pr
-        count = Counter()
+        count: Counter = Counter()
         Pr = np.zeros(ngrid, dtype=np.double)
 
         for it in range(nsteps):
@@ -113,10 +127,12 @@ def __init__(self, structure, trajectories, interval=0.5, species=("Li", "Na")):
         self.lens = lens
         self.Pr = Pr
         self.species = species
-        self.stable_sites = None
+        self.stable_sites: np.ndarray | None = None
 
     @classmethod
-    def from_diffusion_analyzer(cls, diffusion_analyzer, interval=0.5, species=("Li", "Na")):
+    def from_diffusion_analyzer(
+        cls, diffusion_analyzer: DiffusionAnalyzer, interval: float = 0.5, species: Sequence[SpeciesLike] = ("Li", "Na")
+    ) -> ProbabilityDensityAnalysis:
         """
         Create a ProbabilityDensityAnalysis from a diffusion_analyzer object.
 
@@ -128,16 +144,16 @@ def from_diffusion_analyzer(cls, diffusion_analyzer, interval=0.5, species=("Li"
             species(list of str): list of species that are of interest
         """
         structure = diffusion_analyzer.structure
-        trajectories = []
+        _trajectories = []
 
         for _i, s in enumerate(diffusion_analyzer.get_drift_corrected_structures()):
-            trajectories.append(s.frac_coords)
+            _trajectories.append(s.frac_coords)
 
-        trajectories = np.array(trajectories)
+        trajectories = np.array(_trajectories)
 
         return ProbabilityDensityAnalysis(structure, trajectories, interval=interval, species=species)
 
-    def generate_stable_sites(self, p_ratio=0.25, d_cutoff=1.0):
+    def generate_stable_sites(self, p_ratio: float = 0.25, d_cutoff: float = 1.0) -> None:
         """
         Obtain a set of low-energy sites from probability density function with
         given probability threshold 'p_ratio'. The set of grid points with
@@ -157,14 +173,14 @@ def generate_stable_sites(self, p_ratio=0.25, d_cutoff=1.0):
             as a Nx3 numpy array.
         """
         # Set of grid points with high probability density.
-        grid_fcoords = []
+        _grid_fcoords = []
         indices = np.where(self.Pr > self.Pr.max() * p_ratio)
         lattice = self.structure.lattice
 
         for x, y, z in zip(indices[0], indices[1], indices[2]):
-            grid_fcoords.append([x / self.lens[0], y / self.lens[1], z / self.lens[2]])
+            _grid_fcoords.append([x / self.lens[0], y / self.lens[1], z / self.lens[2]])
 
-        grid_fcoords = np.array(grid_fcoords)
+        grid_fcoords = np.array(_grid_fcoords)
         dist_matrix = np.array(lattice.get_all_distances(grid_fcoords, grid_fcoords))
         np.fill_diagonal(dist_matrix, 0)
 
@@ -191,34 +207,35 @@ def generate_stable_sites(self, p_ratio=0.25, d_cutoff=1.0):
         stable_sites = []
 
         for i in set(cluster_indices):
-            indices = np.where(cluster_indices == i)[0]
+            _indices = np.where(cluster_indices == i)[0]
 
-            if len(indices) == 1:
-                stable_sites.append(grid_fcoords[indices[0]])
+            if len(_indices) == 1:
+                stable_sites.append(grid_fcoords[_indices[0]])
                 continue
 
             # Consider periodic boundary condition
-            members = grid_fcoords[indices] - grid_fcoords[indices[0]]
+            members = grid_fcoords[_indices] - grid_fcoords[_indices[0]]
             members = np.where(members > 0.5, members - 1.0, members)
             members = np.where(members < -0.5, members + 1.0, members)
-            members += grid_fcoords[indices[0]]
+            members += grid_fcoords[_indices[0]]
 
             stable_sites.append(np.mean(members, axis=0))
 
         self.stable_sites = np.array(stable_sites)
 
-    def get_full_structure(self):
+    def get_full_structure(self) -> Structure:
         """
         Generate the structure with the low-energy sites included. In the end, a
         pymatgen Structure object will be returned.
         """
         full_structure = self.structure.copy()
+        assert self.stable_sites is not None, "Please run generate_stable_sites() first!"
         for fcoord in self.stable_sites:
             full_structure.append("X", fcoord)
 
         return full_structure
 
-    def to_chgcar(self, filename="CHGCAR.vasp"):
+    def to_chgcar(self, filename: PathLike = "CHGCAR.vasp") -> None:
         """
         Save the probability density distribution in the format of CHGCAR,
         which can be visualized by VESTA.
@@ -282,7 +299,13 @@ class SiteOccupancyAnalyzer:
 
     """
 
-    def __init__(self, structure, coords_ref, trajectories, species=("Li", "Na")):
+    def __init__(
+        self,
+        structure: Structure,
+        coords_ref: np.ndarray | Sequence[Sequence[float]],
+        trajectories: np.ndarray | Sequence[Sequence[Sequence[float]]],
+        species: Sequence[SpeciesLike] = ("Li", "Na"),
+    ) -> None:
         """
         Args:
             structure (pmg_structure): Initial structure.
@@ -296,7 +319,7 @@ def __init__(self, structure, coords_ref, trajectories, species=("Li", "Na")):
         lattice = structure.lattice
         coords_ref = np.array(coords_ref)
         trajectories = np.array(trajectories)
-        count = Counter()
+        count: Counter = Counter()
 
         indices = [i for i, site in enumerate(structure) if site.specie.symbol in species]
 
@@ -318,12 +341,17 @@ def __init__(self, structure, coords_ref, trajectories, species=("Li", "Na")):
         self.nsteps = len(trajectories)
         self.site_occ = site_occ
 
-    def get_average_site_occupancy(self, indices):
+    def get_average_site_occupancy(self, indices: list) -> float:
         """Get the average site occupancy over a subset of reference sites."""
         return np.sum(self.site_occ[indices]) / len(indices)
 
     @classmethod
-    def from_diffusion_analyzer(cls, coords_ref, diffusion_analyzer, species=("Li", "Na")):
+    def from_diffusion_analyzer(
+        cls,
+        coords_ref: np.ndarray | Sequence[Sequence[float]],
+        diffusion_analyzer: DiffusionAnalyzer,
+        species: Sequence[SpeciesLike] = ("Li", "Na"),
+    ) -> SiteOccupancyAnalyzer:
         """
         Create a SiteOccupancyAnalyzer object using a diffusion_analyzer object.
 
diff --git a/pymatgen/analysis/diffusion/aimd/rdf.py b/pymatgen/analysis/diffusion/aimd/rdf.py
index da4d264..42b6902 100644
--- a/pymatgen/analysis/diffusion/aimd/rdf.py
+++ b/pymatgen/analysis/diffusion/aimd/rdf.py
@@ -8,6 +8,7 @@
 from collections import Counter
 from math import ceil
 from multiprocessing import cpu_count
+from typing import TYPE_CHECKING
 
 import numpy as np
 from joblib import Parallel, delayed
@@ -18,6 +19,11 @@
 from pymatgen.core import Structure
 from pymatgen.util.plotting import pretty_plot
 
+if TYPE_CHECKING:
+    from collections.abc import Sequence
+
+    from pymatgen.util.typing import SpeciesLike
+
 
 class RadialDistributionFunction:
     """
@@ -34,7 +40,7 @@ def __init__(
         rmax: float = 10.0,
         cell_range: int = 1,
         sigma: float = 0.1,
-    ):
+    ) -> None:
         """
         Args:
             structures ([Structure]): list of structure
@@ -145,14 +151,14 @@ def __init__(
     @classmethod
     def from_species(
         cls,
-        structures: list,
+        structures: list[Structure],
         ngrid: int = 101,
         rmax: float = 10.0,
         cell_range: int = 1,
         sigma: float = 0.1,
-        species: tuple | list = ("Li", "Na"),
-        reference_species: tuple | list | None = None,
-    ):
+        species: Sequence[SpeciesLike] = ("Li", "Na"),
+        reference_species: Sequence[SpeciesLike] | None = None,
+    ) -> RadialDistributionFunction:
         """
         Initialize using species.
 
@@ -191,7 +197,7 @@ def from_species(
         )
 
     @property
-    def coordination_number(self):
+    def coordination_number(self) -> np.ndarray:
         """
         returns running coordination number.
 
@@ -208,7 +214,7 @@ def get_rdf_plot(
         xlim: tuple = (0.0, 8.0),
         ylim: tuple = (-0.005, 3.0),
         loc_peak: bool = False,
-    ):
+    ) -> np.ndarray:
         """
         Plot the average RDF function.
 
@@ -248,7 +254,7 @@ def get_rdf_plot(
 
         return ax
 
-    def export_rdf(self, filename: str):
+    def export_rdf(self, filename: str) -> None:
         """
         Output RDF data to a csv file.
 
@@ -280,8 +286,8 @@ def __init__(
         rmax: float = 10.0,
         ngrid: float = 101,
         sigma: float = 0.0,
-        n_jobs=None,
-    ):
+        n_jobs: int | None = None,
+    ) -> None:
         """
         This method calculates rdf on `np.linspace(rmin, rmax, ngrid)` points.
 
@@ -350,7 +356,7 @@ def __init__(
         self.n_structures = len(self.structures)
         self.sigma = ceil(sigma / self.dr)
 
-    def _dist_to_counts(self, d):
+    def _dist_to_counts(self, d: np.ndarray) -> np.ndarray:
         """
         Convert a distance array for counts in the bin.
 
@@ -360,7 +366,11 @@ def _dist_to_counts(self, d):
         Returns:
             1D array of counts in the bins centered on self.r
         """
-        counts = np.zeros((self.ngrid,))
+        counts = np.zeros(
+            np.array(
+                self.ngrid,
+            )
+        )
         indices = np.array(np.floor((d - self.rmin + 0.5 * self.dr) / self.dr), dtype=int)
 
         unique, val_counts = np.unique(indices, return_counts=True)
@@ -371,8 +381,8 @@ def get_rdf(
         self,
         ref_species: str | list[str],
         species: str | list[str],
-        is_average=True,
-    ):
+        is_average: bool = True,
+    ) -> tuple[np.ndarray, np.ndarray]:
         """
         Args:
             ref_species (list of species or just single specie str): the reference species.
@@ -400,8 +410,8 @@ def get_one_rdf(
         self,
         ref_species: str | list[str],
         species: str | list[str],
-        index=0,
-    ):
+        index: int = 0,
+    ) -> tuple[np.ndarray, np.ndarray]:
         """
         Get the RDF for one structure, indicated by the index of the structure
         in all structures.
@@ -439,7 +449,12 @@ def get_one_rdf(
             rdf_temp = gaussian_filter1d(rdf_temp, self.sigma)
         return self.r, rdf_temp
 
-    def get_coordination_number(self, ref_species, species, is_average=True):
+    def get_coordination_number(
+        self,
+        ref_species: SpeciesLike | list[SpeciesLike],
+        species: SpeciesLike | list[SpeciesLike],
+        is_average: bool = True,
+    ) -> tuple[np.ndarray, list]:
         """
         returns running coordination number.
 
@@ -464,7 +479,7 @@ def get_coordination_number(self, ref_species, species, is_average=True):
         return self.r, cn
 
 
-def _get_neighbor_list(structure, r) -> tuple:
+def _get_neighbor_list(structure: Structure, r: float) -> tuple:
     """
     Thin wrapper to enable parallel calculations.
 
diff --git a/pymatgen/analysis/diffusion/aimd/tests/test_clustering.py b/pymatgen/analysis/diffusion/aimd/tests/test_clustering.py
index 3db9e16..7ae59b7 100644
--- a/pymatgen/analysis/diffusion/aimd/tests/test_clustering.py
+++ b/pymatgen/analysis/diffusion/aimd/tests/test_clustering.py
@@ -14,14 +14,14 @@
 
 
 class KmeansTest(unittest.TestCase):
-    def test_cluster(self):
+    def test_cluster(self) -> None:
         data = np.random.uniform(size=(10, 5))
-        data = list(data)
+        _data = list(data)
         d2 = np.random.uniform(size=(10, 5)) + ([5] * 5)
-        data.extend(list(d2))
+        _data.extend(list(d2))
         d2 = np.random.uniform(size=(10, 5)) + ([-5] * 5)
-        data.extend(list(d2))
-        data = np.array(data)
+        _data.extend(list(d2))
+        data = np.array(_data)
 
         k = Kmeans()
         clusters = []
@@ -38,15 +38,15 @@ def test_cluster(self):
 
 
 class KmeansPBCTest(unittest.TestCase):
-    def test_cluster(self):
+    def test_cluster(self) -> None:
         lattice = Lattice.cubic(4)
 
-        pts = []
+        _pts = []
         initial = [[0, 0, 0], [0.5, 0.5, 0.5], [0.25, 0.25, 0.25], [0.5, 0, 0]]
         for c in initial:
             for _i in range(100):
-                pts.append(np.array(c) + np.random.randn(3) * 0.01 + np.random.randint(3))
-        pts = np.array(pts)
+                _pts.append(np.array(c) + np.random.randn(3) * 0.01 + np.random.randint(3))
+        pts = np.array(_pts)
         k = KmeansPBC(lattice)
         centroids, labels, ss = k.cluster(pts, 4)
         for c1 in centroids:
diff --git a/pymatgen/analysis/diffusion/aimd/tests/test_pathway.py b/pymatgen/analysis/diffusion/aimd/tests/test_pathway.py
index 8afc88b..e813f78 100644
--- a/pymatgen/analysis/diffusion/aimd/tests/test_pathway.py
+++ b/pymatgen/analysis/diffusion/aimd/tests/test_pathway.py
@@ -20,7 +20,7 @@
 
 
 class ProbabilityDensityTest(unittest.TestCase):
-    def test_probability(self):
+    def test_probability(self) -> None:
         traj_file = os.path.join(tests_dir, "cNa3PS4_trajectories.npy")
         struc_file = os.path.join(tests_dir, "cNa3PS4.cif")
 
@@ -36,7 +36,7 @@ def test_probability(self):
         self.assertAlmostEqual(pda.Pr.min(), 0.0, 12)
         self.assertAlmostEqual(Pr_tot, 1.0, 12)
 
-    def test_probability_classmethod(self):
+    def test_probability_classmethod(self) -> None:
         file = os.path.join(tests_dir, "cNa3PS4_pda.json")
         data = json.load(open(file))
         diff_analyzer = DiffusionAnalyzer.from_dict(data)
@@ -50,7 +50,7 @@ def test_probability_classmethod(self):
         self.assertAlmostEqual(pda.Pr.min(), 0.0, 12)
         self.assertAlmostEqual(Pr_tot, 1.0, 12)
 
-    def test_generate_stable_sites(self):
+    def test_generate_stable_sites(self) -> None:
         file = os.path.join(tests_dir, "cNa3PS4_pda.json")
         data = json.load(open(file))
         diff_analyzer = DiffusionAnalyzer.from_dict(data)
@@ -59,6 +59,7 @@ def test_generate_stable_sites(self):
         pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diffusion_analyzer=diff_analyzer, interval=0.1)
         pda.generate_stable_sites(p_ratio=0.25, d_cutoff=1.5)
 
+        assert pda.stable_sites is not None
         assert len(pda.stable_sites) == 50
         self.assertAlmostEqual(pda.stable_sites[1][2], 0.24113475177304966, 8)
         self.assertAlmostEqual(pda.stable_sites[7][1], 0.5193661971830985, 8)
@@ -69,7 +70,7 @@ def test_generate_stable_sites(self):
         assert s.composition["X"] == 50
         self.assertAlmostEqual(s[177].frac_coords[2], 0.57446809)
 
-    def test_to_chgcar(self):
+    def test_to_chgcar(self) -> None:
         file = os.path.join(tests_dir, "cNa3PS4_pda.json")
         data = json.load(open(file))
         diff_analyzer = DiffusionAnalyzer.from_dict(data)
@@ -83,7 +84,7 @@ def test_to_chgcar(self):
 
 
 class SiteOccupancyTest(unittest.TestCase):
-    def test_site_occupancy(self):
+    def test_site_occupancy(self) -> None:
         traj_file = os.path.join(tests_dir, "cNa3PS4_trajectories.npy")
         struc_file = os.path.join(tests_dir, "cNa3PS4.cif")
 
@@ -95,12 +96,12 @@ def test_site_occupancy(self):
         # SiteOccupancyAnalyzer object
         socc = SiteOccupancyAnalyzer(structure, coords_ref, trajectories, species=("Li", "Na"))
         site_occ = socc.site_occ
-        self.assertAlmostEqual(np.sum(site_occ), len(coords_ref), 12)
+        self.assertAlmostEqual(int(np.sum(site_occ)), len(coords_ref), 12)
         self.assertAlmostEqual(site_occ[11], 0.98, 12)
         self.assertAlmostEqual(site_occ[15], 0.875, 12)
         assert len(coords_ref) == 48
 
-    def test_site_occupancy_classmethod(self):
+    def test_site_occupancy_classmethod(self) -> None:
         file = os.path.join(tests_dir, "cNa3PS4_pda.json")
         data = json.load(open(file))
         diff_analyzer = DiffusionAnalyzer.from_dict(data)
@@ -111,7 +112,7 @@ def test_site_occupancy_classmethod(self):
         # SiteOccupancyAnalyzer object
         socc = SiteOccupancyAnalyzer.from_diffusion_analyzer(coords_ref, diffusion_analyzer=diff_analyzer)
         site_occ = socc.site_occ
-        self.assertAlmostEqual(np.sum(site_occ), len(coords_ref), 12)
+        self.assertAlmostEqual(int(np.sum(site_occ)), len(coords_ref), 12)
         self.assertAlmostEqual(site_occ[1], 0.98, 12)
         self.assertAlmostEqual(site_occ[26], 0.97, 12)
         assert len(coords_ref) == 48
diff --git a/pymatgen/analysis/diffusion/aimd/tests/test_rdf.py b/pymatgen/analysis/diffusion/aimd/tests/test_rdf.py
index eb0b30d..adba238 100644
--- a/pymatgen/analysis/diffusion/aimd/tests/test_rdf.py
+++ b/pymatgen/analysis/diffusion/aimd/tests/test_rdf.py
@@ -13,7 +13,7 @@
 
 
 class RDFTest(unittest.TestCase):
-    def test_rdf(self):
+    def test_rdf(self) -> None:
         # Parse the DiffusionAnalyzer object from json file directly
         obj = loadfn(os.path.join(tests_dir, "cNa3PS4_pda.json"))
 
@@ -31,7 +31,7 @@ def test_rdf(self):
 
         self.assertAlmostEqual(r[np.argmax(s_na_rdf)], 2.9000, 4)
 
-    def test_rdf_coordination_number(self):
+    def test_rdf_coordination_number(self) -> None:
         # create a simple cubic lattice
         coords = np.array([[0.5, 0.5, 0.5]])
         atom_list = ["S"]
diff --git a/pymatgen/analysis/diffusion/aimd/tests/test_van_hove.py b/pymatgen/analysis/diffusion/aimd/tests/test_van_hove.py
index 1ab3126..3dee26c 100644
--- a/pymatgen/analysis/diffusion/aimd/tests/test_van_hove.py
+++ b/pymatgen/analysis/diffusion/aimd/tests/test_van_hove.py
@@ -23,7 +23,7 @@
 
 
 class VanHoveTest(unittest.TestCase):
-    def test_van_hove(self):
+    def test_van_hove(self) -> None:
         # Parse the DiffusionAnalyzer object from json file directly
         obj = loadfn(os.path.join(tests_dir, "cNa3PS4_pda.json"))
 
@@ -48,7 +48,7 @@ def test_van_hove(self):
 
 
 class RDFTest(unittest.TestCase):
-    def test_rdf(self):
+    def test_rdf(self) -> None:
         # Parse the DiffusionAnalyzer object from json file directly
         obj = loadfn(os.path.join(tests_dir, "cNa3PS4_pda.json"))
 
@@ -118,7 +118,7 @@ def test_rdf(self):
         self.assertAlmostEqual(obj_1.rho, obj_2.rho)
         self.assertAlmostEqual(obj_1.rdf[0], obj_2.rdf[0])
 
-    def test_rdf_coordination_number(self):
+    def test_rdf_coordination_number(self) -> None:
         # create a simple cubic lattice
         coords = np.array([[0.5, 0.5, 0.5]])
         atom_list = ["S"]
@@ -129,7 +129,7 @@ def test_rdf_coordination_number(self):
         )
         assert rdf.coordination_number[100] == 6.0
 
-    def test_rdf_two_species_coordination_number(self):
+    def test_rdf_two_species_coordination_number(self) -> None:
         # create a structure with interpenetrating simple cubic lattice
         coords = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])
         atom_list = ["S", "Zn"]
@@ -145,25 +145,25 @@ def test_rdf_two_species_coordination_number(self):
         )
         assert rdf.coordination_number[100] == 8.0
 
-    def setUp(self):
+    def setUp(self) -> None:
         coords = np.array([[0.5, 0.5, 0.5]])
         atom_list = ["S"]
         lattice = Lattice.from_parameters(a=1.0, b=1.0, c=1.0, alpha=90, beta=90, gamma=90)
         self.structure = Structure(lattice, atom_list, coords)
 
-    def test_raises_valueerror_if_ngrid_is_less_than_2(self):
+    def test_raises_valueerror_if_ngrid_is_less_than_2(self) -> None:
         with pytest.raises(ValueError):
             RadialDistributionFunction.from_species(structures=[self.structure], ngrid=1)
 
-    def test_raises_ValueError_if_sigma_is_not_positive(self):
+    def test_raises_ValueError_if_sigma_is_not_positive(self) -> None:
         with pytest.raises(ValueError):
             RadialDistributionFunction.from_species(structures=[self.structure], sigma=0)
 
-    def test_raises_ValueError_if_species_not_in_structure(self):
+    def test_raises_ValueError_if_species_not_in_structure(self) -> None:
         with pytest.raises(ValueError):
             RadialDistributionFunction.from_species(structures=[self.structure], species=["Cl"])
 
-    def test_raises_ValueError_if_reference_species_not_in_structure(self):
+    def test_raises_ValueError_if_reference_species_not_in_structure(self) -> None:
         with pytest.raises(ValueError):
             RadialDistributionFunction.from_species(
                 structures=[self.structure], species=["S"], reference_species=["Cl"]
@@ -171,7 +171,7 @@ def test_raises_ValueError_if_reference_species_not_in_structure(self):
 
 
 class EvolutionAnalyzerTest(unittest.TestCase):
-    def test_get_df(self):
+    def test_get_df(self) -> None:
         # Parse the DiffusionAnalyzer object from json file directly
         obj = loadfn(os.path.join(tests_dir, "cNa3PS4_pda.json"))
 
diff --git a/pymatgen/analysis/diffusion/aimd/van_hove.py b/pymatgen/analysis/diffusion/aimd/van_hove.py
index 204bd4e..0880c05 100644
--- a/pymatgen/analysis/diffusion/aimd/van_hove.py
+++ b/pymatgen/analysis/diffusion/aimd/van_hove.py
@@ -9,6 +9,7 @@
 from collections import Counter
 from typing import TYPE_CHECKING, Callable
 
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
@@ -19,8 +20,13 @@
 from .rdf import RadialDistributionFunction
 
 if TYPE_CHECKING:
+    from collections.abc import Sequence
+
+    from matplotlib.axes import Axes
+
     from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
     from pymatgen.core import Structure
+    from pymatgen.util.typing import SpeciesLike
 
 __author__ = "Iek-Heng Chu"
 __version__ = "1.0"
@@ -48,10 +54,10 @@ def __init__(
         step_skip: int = 50,
         sigma: float = 0.1,
         cell_range: int = 1,
-        species: tuple | list = ("Li", "Na"),
-        reference_species: tuple | list | None = None,
+        species: Sequence[SpeciesLike] = ("Li", "Na"),
+        reference_species: Sequence | None = None,
         indices: list | None = None,
-    ):
+    ) -> None:
         """
         Initiation.
 
@@ -198,7 +204,7 @@ def __init__(
         # time interval (in ps) in gsrt and gdrt.
         self.timeskip = self.obj.time_step * self.obj.step_skip * step_skip / 1000.0
 
-    def get_3d_plot(self, figsize: tuple = (12, 8), mode: str = "distinct"):
+    def get_3d_plot(self, figsize: tuple = (12, 8), mode: str = "distinct") -> Axes:
         """
         Plot 3D self-part or distinct-part of van Hove function, which is
         specified by the input argument 'type'.
@@ -240,7 +246,8 @@ def get_3d_plot(self, figsize: tuple = (12, 8), mode: str = "distinct"):
         plt.pcolor(X, Y, grt, cmap="jet", vmin=grt.min(), vmax=vmax)
         ax.set_xlabel("Time (ps)", size=labelsize)
         ax.set_ylabel(r"$r$ ($\AA$)", size=labelsize)
-        ax.axis([x.min(), x.max(), y.min(), y.max()])
+        ax.set_xlim(x.min(), x.max())
+        ax.set_ylim(y.min(), y.max())
 
         cbar = plt.colorbar(ticks=cb_ticks)
         cbar.set_label(label=cb_label, size=labelsize)
@@ -253,7 +260,7 @@ def get_1d_plot(
         mode: str = "distinct",
         times: list | None = None,
         colors: list | None = None,
-    ):
+    ) -> Axes:
         """
         Plot the van Hove function at given r or t.
 
@@ -307,7 +314,7 @@ def get_1d_plot(
 class EvolutionAnalyzer:
     """Analyze the evolution of structures during AIMD simulations."""
 
-    def __init__(self, structures: list, rmax: float = 10, step: int = 1, time_step: int = 2):
+    def __init__(self, structures: list, rmax: float = 10, step: int = 1, time_step: int = 2) -> None:
         """
         Initialization the EvolutionAnalyzer from MD simulations. From the
         structures obtained from MD simulations, we can analyze the structure
@@ -336,7 +343,7 @@ def __init__(self, structures: list, rmax: float = 10, step: int = 1, time_step:
         self.time_step = time_step
 
     @staticmethod
-    def get_pairs(structure: Structure):
+    def get_pairs(structure: Structure) -> list:
         """
         Get all element pairs in a structure.
 
@@ -352,7 +359,7 @@ def get_pairs(structure: Structure):
         return list(pairs)
 
     @staticmethod
-    def rdf(structure: Structure, pair: tuple, ngrid: int = 101, rmax: float = 10):
+    def rdf(structure: Structure, pair: tuple, ngrid: int = 101, rmax: float = 10) -> np.ndarray:
         """
         Process rdf from a given structure and pair.
 
@@ -383,7 +390,7 @@ def atom_dist(
         ngrid: int = 101,
         window: float = 1,
         direction: str = "c",
-    ):
+    ) -> np.ndarray:
         """
         Get atomic distribution for a given specie.
 
@@ -419,7 +426,7 @@ def atom_dist(
 
         return np.array(density)
 
-    def get_df(self, func: Callable, save_csv: str | None = None, **kwargs):
+    def get_df(self, func: Callable, save_csv: str | None = None, **kwargs) -> pd.DataFrame:
         """
         Get the data frame for a given pair. This step would be very slow if
         there are hundreds or more structures to parse.
@@ -460,7 +467,7 @@ def get_df(self, func: Callable, save_csv: str | None = None, **kwargs):
         return df
 
     @staticmethod
-    def get_min_dist(df: pd.DataFrame, tol: float = 1e-10):
+    def get_min_dist(df: pd.DataFrame, tol: float = 1e-10) -> float:
         """
         Get the shortest pair distance from the given DataFrame.
 
@@ -484,8 +491,8 @@ def plot_evolution_from_data(
         df: pd.DataFrame,
         x_label: str | None = None,
         cb_label: str | None = None,
-        cmap=plt.cm.plasma,  # pylint: disable=E1101
-    ):
+        cmap: mpl.colors.Colormap = plt.cm.plasma,  # type: ignore[attr-defined] # pylint: disable=E1101
+    ) -> Axes:
         """
         Plot the evolution with time for a given DataFrame. It can be RDF,
         atomic distribution or other characterization data we might
@@ -523,11 +530,11 @@ def plot_evolution_from_data(
             rasterized=True,
         )
         ax.set_ylim(ax.get_ylim()[::-1])
-        ax.collections[0].colorbar.set_label(cb_label, fontsize=30)
+        ax.collections[0].colorbar.set_label(cb_label, fontsize=30) if cb_label and ax.collections[0].colorbar else None
 
         plt.xticks(rotation="horizontal")
 
-        ax.set_xlabel(x_label, fontsize=30)
+        ax.set_xlabel(x_label, fontsize=30) if x_label else None
         ax.set_ylabel("Time (ps)", fontsize=30)
 
         plt.yticks(rotation="horizontal")
@@ -537,9 +544,9 @@ def plot_evolution_from_data(
     def plot_rdf_evolution(
         self,
         pair: tuple,
-        cmap=plt.cm.plasma,  # pylint: disable=E1101
+        cmap: mpl.colors.Colormap = plt.cm.plasma,  # type: ignore[attr-defined] # pylint: disable=E1101
         df: pd.DataFrame = None,
-    ):
+    ) -> Axes:
         """
         Plot the RDF evolution with time for a given pair.
 
@@ -562,9 +569,9 @@ def plot_atomic_evolution(
         self,
         specie: str,
         direction: str = "c",
-        cmap=plt.cm.Blues,  # pylint: disable=E1101
+        cmap: mpl.colors.Colormap = plt.cm.Blues,  # type: ignore[attr-defined] # pylint: disable=E1101
         df: pd.DataFrame = None,
-    ):
+    ) -> Axes:
         """
         Plot the atomic distribution evolution with time for a given species.
 
diff --git a/pymatgen/analysis/diffusion/analyzer.py b/pymatgen/analysis/diffusion/analyzer.py
index 6a731c8..66b7ede 100644
--- a/pymatgen/analysis/diffusion/analyzer.py
+++ b/pymatgen/analysis/diffusion/analyzer.py
@@ -20,6 +20,7 @@
 import warnings
 from typing import TYPE_CHECKING, Literal
 
+import matplotlib.pyplot as plt
 import numpy as np
 import scipy.constants as const
 from monty.json import MSONable
@@ -32,7 +33,11 @@
 from pymatgen.util.coord import pbc_diff
 
 if TYPE_CHECKING:
-    from collections.abc import Sequence
+    from collections.abc import Generator, Sequence
+
+    from matplotlib.axes import Axes
+
+    from pymatgen.util.typing import PathLike, SpeciesLike
 
 __author__ = "Will Richards, Shyue Ping Ong"
 __version__ = "0.2"
@@ -132,20 +137,20 @@ class DiffusionAnalyzer(MSONable):
 
     def __init__(
         self,
-        structure,
-        displacements,
-        specie,
-        temperature,
-        time_step,
-        step_skip,
-        smoothed="max",
-        min_obs=30,
-        avg_nsteps=1000,
-        lattices=None,
-        c_ranges=None,
-        c_range_include_edge=False,
-        structures=None,
-    ):
+        structure: Structure,
+        displacements: np.ndarray,
+        specie: SpeciesLike,
+        temperature: float,
+        time_step: int,
+        step_skip: int,
+        smoothed: bool | str = "max",
+        min_obs: int = 30,
+        avg_nsteps: int = 1000,
+        lattices: np.ndarray | None = None,
+        c_ranges: Sequence | None = None,
+        c_range_include_edge: bool = False,
+        structures: Sequence[Structure] | None = None,
+    ) -> None:
         """
         This constructor is meant to be used with pre-processed data.
         Other convenient constructors are provided as class methods (see
@@ -222,12 +227,9 @@ def __init__(
         self.min_obs = min_obs
         self.smoothed = smoothed
         self.avg_nsteps = avg_nsteps
-        self.lattices = lattices
-
-        if lattices is None:
-            self.lattices = np.array([structure.lattice.matrix.tolist()])
+        self.lattices = lattices if lattices is not None else np.array([structure.lattice.matrix.tolist()])
 
-        indices = []
+        indices: list = []
         framework_indices = []
         for i, site in enumerate(structure):
             if site.specie.symbol == specie:
@@ -297,7 +299,7 @@ def __init__(
                 msd_components[i] = np.average(dcomponents[indices] ** 2, axis=(0, 1))
 
                 # Get regional msd
-                if c_ranges:
+                if c_ranges and structures:
                     if not c_range_include_edge:
                         for index in indices:
                             if any(
@@ -383,7 +385,9 @@ def __init__(
             self.indices = indices
             self.framework_indices = framework_indices
 
-    def get_drift_corrected_structures(self, start=None, stop=None, step=None):
+    def get_drift_corrected_structures(
+        self, start: int | None = None, stop: int | None = None, step: int | None = None
+    ) -> Generator:
         """
         Returns an iterator for the drift-corrected structures. Use of
         iterator is to reduce memory usage as # of structures in MD can be
@@ -410,7 +414,7 @@ def get_drift_corrected_structures(self, start=None, stop=None, step=None):
                 coords_are_cartesian=True,
             )
 
-    def get_summary_dict(self, include_msd_t=False, include_mscd_t=False):
+    def get_summary_dict(self, include_msd_t: bool = False, include_mscd_t: bool = False) -> dict:
         """
         Provides a summary of diffusion information.
 
@@ -450,7 +454,7 @@ def get_summary_dict(self, include_msd_t=False, include_mscd_t=False):
             d["mscd"] = self.mscd.tolist()
         return d
 
-    def get_framework_rms_plot(self, granularity=200, matching_s=None):
+    def get_framework_rms_plot(self, granularity: int = 200, matching_s: Structure = None) -> Axes:
         """
         Get the plot of rms framework displacement vs time. Useful for checking
         for melting, especially if framework atoms can move via paddle-wheel
@@ -481,22 +485,22 @@ def get_framework_rms_plot(self, granularity=200, matching_s=None):
             comparator=OrderDisorderElementComparator(),
             allow_subset=True,
         )
-        rms = []
+        _rms = []
         for s in self.get_drift_corrected_structures(step=step):
             s.remove_species([self.specie])
             d = sm.get_rms_dist(f, s)
             if d:
-                rms.append(d)
+                _rms.append(d)
             else:
-                rms.append((1, 1))
-        max_dt = (len(rms) - 1) * step * self.step_skip * self.time_step
+                _rms.append((1, 1))
+        max_dt = (len(_rms) - 1) * step * self.step_skip * self.time_step
         if max_dt > 100000:
-            plot_dt = np.linspace(0, max_dt / 1000, len(rms))
+            plot_dt = np.linspace(0, max_dt / 1000, len(_rms))
             unit = "ps"
         else:
-            plot_dt = np.linspace(0, max_dt, len(rms))
+            plot_dt = np.linspace(0, max_dt, len(_rms))
             unit = "fs"
-        rms = np.array(rms)
+        rms = np.array(_rms)
         ax.plot(plot_dt, rms[:, 0], label="RMS")
         ax.plot(plot_dt, rms[:, 1], label="max")
         ax.legend(loc="best")
@@ -504,7 +508,7 @@ def get_framework_rms_plot(self, granularity=200, matching_s=None):
         ax.set_ylabel("normalized distance")
         return ax
 
-    def get_msd_plot(self, mode="specie"):
+    def get_msd_plot(self, mode: str = "specie") -> Axes:
         """
         Get the plot of the smoothed msd vs time graph. Useful for
         checking convergence. This can be written to an image file.
@@ -517,6 +521,7 @@ def get_msd_plot(self, mode="specie"):
         from pymatgen.util.plotting import pretty_plot
 
         ax = pretty_plot(12, 8)
+        plot_dt: np.ndarray
         if np.max(self.dt) > 100000:
             plot_dt = self.dt / 1000
             unit = "ps"
@@ -553,7 +558,7 @@ def get_msd_plot(self, mode="specie"):
             ax.set_ylabel("MSD ($\\AA^2$)")
         return ax
 
-    def plot_msd(self, mode="default"):
+    def plot_msd(self, mode: str = "default") -> None:
         """
         Plot the smoothed msd vs time graph. Useful for checking convergence.
 
@@ -564,9 +569,10 @@ def plot_msd(self, mode="default"):
                 displacement by specie), or "mscd" (overall mean square charge
                 displacement for diffusing specie).
         """
-        self.get_msd_plot(mode=mode).show()
+        self.get_msd_plot(mode=mode)
+        plt.show()
 
-    def export_msdt(self, filename):
+    def export_msdt(self, filename: str) -> None:
         """
         Writes MSD data to a csv file that can be easily plotted in other
         software.
@@ -590,15 +596,15 @@ def export_msdt(self, filename):
     @classmethod
     def from_structures(
         cls,
-        structures,
-        specie,
-        temperature,
-        time_step,
-        step_skip,
-        initial_disp=None,
-        initial_structure=None,
+        structures: Sequence[Structure],
+        specie: SpeciesLike,
+        temperature: float,
+        time_step: int,
+        step_skip: int,
+        initial_disp: np.ndarray = None,
+        initial_structure: Structure = None,
         **kwargs,
-    ):
+    ) -> DiffusionAnalyzer:
         r"""
         Convenient constructor that takes in a list of Structure objects to
         perform diffusion analysis.
@@ -627,29 +633,29 @@ def from_structures(
             **kwargs: kwargs supported by the :class:`DiffusionAnalyzer`_.
                 Examples include smoothed, min_obs, avg_nsteps.
         """
-        p, lattices = [], []
+        _p, _lattices = [], []
         structure = structures[0]
         for s in structures:
-            p.append(np.array(s.frac_coords)[:, None])
-            lattices.append(s.lattice.matrix)
+            _p.append(np.array(s.frac_coords)[:, None])
+            _lattices.append(s.lattice.matrix)
         if initial_structure is not None:
-            p.insert(0, np.array(initial_structure.frac_coords)[:, None])
-            lattices.insert(0, initial_structure.lattice.matrix)
+            _p.insert(0, np.array(initial_structure.frac_coords)[:, None])
+            _lattices.insert(0, initial_structure.lattice.matrix)
         else:
-            p.insert(0, p[0])
-            lattices.insert(0, lattices[0])
+            _p.insert(0, _p[0])
+            _lattices.insert(0, _lattices[0])
 
-        p = np.concatenate(p, axis=1)
+        p = np.concatenate(_p, axis=1)
         dp = p[:, 1:] - p[:, :-1]
         dp = dp - np.round(dp)
         f_disp = np.cumsum(dp, axis=1)
         c_disp = []
         for i in f_disp:
-            c_disp.append([np.dot(d, m) for d, m in zip(i, lattices[1:])])
+            c_disp.append([np.dot(d, m) for d, m in zip(i, _lattices[1:])])
         disp = np.array(c_disp)
 
         # If is NVT-AIMD, clear lattice data.
-        lattices = np.array([lattices[0]]) if np.array_equal(lattices[0], lattices[-1]) else np.array(lattices)
+        lattices = np.array([_lattices[0]]) if np.array_equal(_lattices[0], _lattices[-1]) else np.array(_lattices)
         if initial_disp is not None:
             disp += initial_disp[:, None, :]
 
@@ -665,7 +671,14 @@ def from_structures(
         )
 
     @classmethod
-    def from_vaspruns(cls, vaspruns, specie, initial_disp=None, initial_structure=None, **kwargs):
+    def from_vaspruns(
+        cls,
+        vaspruns: Sequence[Vasprun],
+        specie: SpeciesLike,
+        initial_disp: np.ndarray = None,
+        initial_structure: Structure = None,
+        **kwargs,
+    ) -> DiffusionAnalyzer:
         r"""
         Convenient constructor that takes in a list of Vasprun objects to
         perform diffusion analysis.
@@ -690,7 +703,7 @@ def from_vaspruns(cls, vaspruns, specie, initial_disp=None, initial_structure=No
                 Examples include smoothed, min_obs, avg_nsteps.
         """
 
-        def get_structures(vaspruns):
+        def get_structures(vaspruns: Sequence[Vasprun]) -> Generator:
             step_skip = vaspruns[0].ionic_step_skip or 1
             final_structure = vaspruns[0].initial_structure
             temperature = vaspruns[0].parameters["TEEND"]
@@ -725,14 +738,14 @@ def get_structures(vaspruns):
     @classmethod
     def from_files(
         cls,
-        filepaths,
-        specie,
-        step_skip=10,
-        ncores=None,
-        initial_disp=None,
-        initial_structure=None,
+        filepaths: Sequence[PathLike],
+        specie: SpeciesLike,
+        step_skip: int = 10,
+        ncores: int | None = None,
+        initial_disp: np.ndarray = None,
+        initial_structure: Structure = None,
         **kwargs,
-    ):
+    ) -> DiffusionAnalyzer:
         r"""
         Convenient constructor that takes in a list of vasprun.xml paths to
         perform diffusion analysis.
@@ -772,14 +785,14 @@ def from_files(
             with multiprocessing.Pool(ncores) as p:
                 vaspruns = p.imap(_get_vasprun, [(fp, step_skip) for fp in filepaths])
                 return cls.from_vaspruns(
-                    vaspruns,
+                    list(vaspruns),
                     specie=specie,
                     initial_disp=initial_disp,
                     initial_structure=initial_structure,
                     **kwargs,
                 )
 
-        def vr(filepaths):
+        def vr(filepaths: Sequence[PathLike]) -> Generator[Vasprun, None, None]:
             offset = 0
             for p in filepaths:
                 v = Vasprun(p, ionic_step_offset=offset, ionic_step_skip=step_skip)
@@ -788,14 +801,14 @@ def vr(filepaths):
                 offset = (-(v.nionic_steps - offset)) % step_skip
 
         return cls.from_vaspruns(
-            vr(filepaths),
+            list(vr(filepaths)),
             specie=specie,
             initial_disp=initial_disp,
             initial_structure=initial_structure,
             **kwargs,
         )
 
-    def as_dict(self):
+    def as_dict(self) -> dict:
         """Returns: MSONable dict."""
         return {
             "@module": self.__class__.__module__,
@@ -813,7 +826,7 @@ def as_dict(self):
         }
 
     @classmethod
-    def from_dict(cls, d):
+    def from_dict(cls, d: dict) -> DiffusionAnalyzer:
         """
         Args:
             d (dict): Dict representation.
@@ -835,7 +848,7 @@ def from_dict(cls, d):
         )
 
 
-def get_conversion_factor(structure, species, temperature):
+def get_conversion_factor(structure: Structure, species: SpeciesLike, temperature: float) -> float:
     """
     Conversion factor to convert between cm^2/s diffusivity measurements and
     mS/cm conductivity measurements based on number of atoms of diffusing
@@ -861,16 +874,16 @@ def get_conversion_factor(structure, species, temperature):
     return 1000 * n / (vol * const.N_A) * z**2 * (const.N_A * const.e) ** 2 / (const.R * temperature)
 
 
-def _get_vasprun(args):
+def _get_vasprun(args: tuple[str, int]) -> Vasprun:
     """Internal method to support multiprocessing."""
     return Vasprun(args[0], ionic_step_skip=args[1], parse_dos=False, parse_eigen=False)
 
 
 def fit_arrhenius(
-    temps: Sequence[float],
-    diffusivities: Sequence[float],
+    temps: Sequence[float] | np.ndarray,
+    diffusivities: Sequence[float] | np.ndarray,
     mode: Literal["linear", "exp"] = "linear",
-    diffusivity_errors: Sequence[float] | None = None,
+    diffusivity_errors: Sequence[float] | np.ndarray | None = None,
 ) -> tuple[float, float, float | None]:
     """
     Returns Ea, c, standard error of Ea from the Arrhenius fit.
@@ -901,7 +914,7 @@ def fit_arrhenius(
         return -w[0] * const.k / const.e, np.exp(w[1]), std_Ea
     if mode == "exp":
 
-        def arrhenius(t, Ea, c):
+        def arrhenius(t: np.ndarray, Ea: float, c: float) -> np.ndarray:
             return c * np.exp(-Ea / (const.k / const.e * t))
 
         guess = fit_arrhenius(temps, diffusivities, mode="linear")[0:2]  # Use linear fit to get initial guess
@@ -911,7 +924,7 @@ def arrhenius(t, Ea, c):
     return None
 
 
-def get_diffusivity_from_msd(msd, dt, smoothed="max"):
+def get_diffusivity_from_msd(msd: np.ndarray, dt: np.ndarray, smoothed: bool | str = "max") -> tuple[float, float]:
     """
     Returns diffusivity and standard deviation of diffusivity.
 
@@ -940,7 +953,7 @@ def get_diffusivity_from_msd(msd, dt, smoothed="max"):
                smoothing.
     """
 
-    def weighted_lstsq(a, b):
+    def weighted_lstsq(a: np.ndarray, b: np.ndarray) -> tuple:
         if smoothed == "max":
             # For max smoothing, we need to weight by variance.
             w_root = (1 / dt) ** 0.5
@@ -967,12 +980,14 @@ def weighted_lstsq(a, b):
     # Pre-compute the denominator since we will use it later.
     # We divide dt by 1000 to avoid overflow errors in some systems (
     # e.g., win). This is subsequently corrected where denom is used.
-    denom = (n * np.sum((dt / 1000) ** 2) - np.sum(dt / 1000) ** 2) * (n - 2)
+    denom = (n * np.sum((np.array(dt) / 1000) ** 2) - np.sum(np.array(dt) / 1000) ** 2) * (n - 2)
     diffusivity_std_dev = np.sqrt(n * res[0] / denom) / 60 / 1000
     return diffusivity, diffusivity_std_dev
 
 
-def get_extrapolated_diffusivity(temps, diffusivities, new_temp, mode: Literal["linear", "exp"] = "linear"):
+def get_extrapolated_diffusivity(
+    temps: Sequence[float], diffusivities: Sequence[float], new_temp: float, mode: Literal["linear", "exp"] = "linear"
+) -> float:
     """
     Returns (Arrhenius) extrapolated diffusivity at new_temp.
 
@@ -990,7 +1005,13 @@ def get_extrapolated_diffusivity(temps, diffusivities, new_temp, mode: Literal["
     return c * np.exp(-Ea / (const.k / const.e * new_temp))
 
 
-def get_extrapolated_conductivity(temps, diffusivities, new_temp, structure, species):
+def get_extrapolated_conductivity(
+    temps: Sequence[float],
+    diffusivities: Sequence[float],
+    new_temp: float,
+    structure: Structure,
+    species: SpeciesLike,
+) -> float:
     """
     Returns extrapolated mS/cm conductivity.
 
@@ -1011,13 +1032,13 @@ def get_extrapolated_conductivity(temps, diffusivities, new_temp, structure, spe
 
 
 def get_arrhenius_plot(
-    temps,
-    diffusivities,
-    diffusivity_errors=None,
+    temps: Sequence[float] | np.ndarray,
+    diffusivities: Sequence[float] | np.ndarray,
+    diffusivity_errors: Sequence[float] | np.ndarray | None = None,
     mode: Literal["linear", "exp"] = "linear",
     unit: Literal["eV", "meV"] = "meV",
     **kwargs,
-):
+) -> Axes:
     r"""
     Returns an Arrhenius plot.
 
diff --git a/pymatgen/analysis/diffusion/neb/full_path_mapper.py b/pymatgen/analysis/diffusion/neb/full_path_mapper.py
index 736281a..a61373c 100644
--- a/pymatgen/analysis/diffusion/neb/full_path_mapper.py
+++ b/pymatgen/analysis/diffusion/neb/full_path_mapper.py
@@ -29,12 +29,15 @@
 from pymatgen.symmetry.structure import SymmetrizedStructure
 
 if TYPE_CHECKING:
+    from collections.abc import Generator
+
     from pymatgen.entries.computed_entries import ComputedEntry, ComputedStructureEntry
+    from pymatgen.util.typing import SpeciesLike
 
 logger = logging.getLogger(__name__)
 
 
-def generic_groupby(list_in: list, comp: Callable = operator.eq):
+def generic_groupby(list_in: list, comp: Callable = operator.eq) -> list:
     """
     Group a list of unsortable objects.
 
@@ -75,9 +78,9 @@ def __init__(
         self,
         structure: Structure,
         m_graph: StructureGraph,
-        symprec=0.1,
-        vac_mode=False,
-    ):
+        symprec: float = 0.1,
+        vac_mode: bool = False,
+    ) -> None:
         """
         Construct the MigrationGraph object using a potential_field will
         all mobile sites occupied. A potential_field graph is generated by
@@ -137,7 +140,7 @@ def symm_structure(self) -> SymmetrizedStructure:
         return sym_struct
 
     @property
-    def unique_hops(self):
+    def unique_hops(self) -> dict:
         """The unique hops dictionary keyed by the hop label."""
         # reversed so that the first instance represents the group of distinct hops
         ihop_data = list(reversed(list(self.m_graph.graph.edges(data=True))))
@@ -183,9 +186,7 @@ def with_local_env_strategy(
         return cls(structure=structure, m_graph=migration_graph, **kwargs)
 
     @classmethod
-    def with_distance(
-        cls, structure: Structure, migrating_specie: str, max_distance: float, **kwargs
-    ) -> MigrationGraph:
+    def with_distance(cls, structure: Structure, migrating_specie: SpeciesLike, max_distance: float, **kwargs):  # noqa: ANN206
         """
         Using a specific nn strategy to get the connectivity graph between all the migrating ion sites.
 
@@ -282,7 +283,7 @@ def get_structure_from_entries(
             res.append(struct)
         return res
 
-    def _get_pos_and_migration_hop(self, u, v, w):
+    def _get_pos_and_migration_hop(self, u: int, v: int, w: int) -> None:
         """
         insert a single MigrationHop object on a graph edge
         Args:
@@ -305,11 +306,11 @@ def _get_pos_and_migration_hop(self, u, v, w):
 
         edge["hop"] = MigrationHop(i_site, e_site, self.symm_structure, symprec=self.symprec)
 
-    def _populate_edges_with_migration_hops(self):
+    def _populate_edges_with_migration_hops(self) -> None:
         """Populate the edges with the data for the Migration Paths."""
         list(starmap(self._get_pos_and_migration_hop, self.m_graph.graph.edges))
 
-    def _group_and_label_hops(self):
+    def _group_and_label_hops(self) -> dict:
         """Group the MigrationHop objects together and label all the symmetrically equlivaelnt hops with the same label."""
         hops = list(nx.get_edge_attributes(self.m_graph.graph, "hop").items())
         labs = generic_groupby(hops, comp=lambda x, y: x[1] == y[1])
@@ -322,7 +323,7 @@ def add_data_to_similar_edges(
         target_label: int | str,
         data: dict,
         m_hop: MigrationHop | None = None,
-    ):
+    ) -> None:
         """
         Insert data to all edges with the same label
         Args:
@@ -346,7 +347,7 @@ def add_data_to_similar_edges(
                             continue
                         d[k] = d[k][::-1]  # flip the data in the array
 
-    def assign_cost_to_graph(self, cost_keys=None):
+    def assign_cost_to_graph(self, cost_keys: list | None = None) -> None:
         """
         Read the data dict on each add and populate a cost key
         Args:
@@ -359,7 +360,7 @@ def assign_cost_to_graph(self, cost_keys=None):
             cost_val = np.prod([v[ik] for ik in cost_keys])
             self.add_data_to_similar_edges(k, {"cost": cost_val})
 
-    def get_path(self, max_val=100000, flip_hops=True):
+    def get_path(self, max_val: float = 100000, flip_hops: bool = True) -> Generator:
         """
         obtain a pathway through the material using hops that are in the current graph
         Basic idea:
@@ -446,7 +447,7 @@ def get_summary_dict(self, added_keys: list[str] | None = None) -> dict:
         if added_keys is not None:
             keys += added_keys
 
-        def get_keys(d):
+        def get_keys(d: dict) -> dict:
             return {k_: d[k_] for k_ in keys if k_ in d}
 
         for u, v, d in self.m_graph.graph.edges(data=True):
@@ -482,7 +483,7 @@ def __init__(
         potential_field: VolumetricData,
         potential_data_key: str,
         **kwargs,
-    ):
+    ) -> None:
         """
         Construct the MigrationGraph object using a VolumetricData object.
         The graph is constructed using the structure, and cost values are assigned based on charge density analysis.
@@ -500,7 +501,7 @@ def __init__(
         super().__init__(structure=structure, m_graph=m_graph, **kwargs)
         self._setup_grids()
 
-    def _setup_grids(self):
+    def _setup_grids(self) -> None:
         """Populate the internal variables used for defining the grid points in the charge density analysis."""
         # set up the grid
         aa = np.linspace(0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False)
@@ -522,7 +523,7 @@ def _setup_grids(self):
         self._fcoords = np.vstack([AA.flatten(), BB.flatten(), CC.flatten()]).T
         self._images = np.vstack([IMA.flatten(), IMB.flatten(), IMC.flatten()]).T
 
-    def _dist_mat(self, pos_frac):
+    def _dist_mat(self, pos_frac: np.ndarray) -> np.ndarray:
         # return a matrix that contains the distances to pos_frac
         aa = np.linspace(0, 1, len(self.potential_field.get_axis_grid(0)), endpoint=False)
         bb = np.linspace(0, 1, len(self.potential_field.get_axis_grid(1)), endpoint=False)
@@ -535,7 +536,7 @@ def _dist_mat(self, pos_frac):
         )
         return dist_from_pos.reshape(AA.shape)
 
-    def _get_pathfinder_from_hop(self, migration_hop: MigrationHop, n_images=20):
+    def _get_pathfinder_from_hop(self, migration_hop: MigrationHop, n_images: int = 20) -> NEBPathfinder:
         # get migration pathfinder objects which contains the paths
         ipos = migration_hop.isite.frac_coords
         epos = migration_hop.esite.frac_coords
@@ -560,7 +561,13 @@ def _get_pathfinder_from_hop(self, migration_hop: MigrationHop, n_images=20):
             mid_struct=mid_struct,
         )
 
-    def _get_avg_chg_at_max(self, migration_hop, radius=None, chg_along_path=False, output_positions=False):
+    def _get_avg_chg_at_max(  # noqa: ANN202
+        self,
+        migration_hop: MigrationHop,
+        radius: float | None = None,
+        chg_along_path: bool = False,
+        output_positions: bool = False,
+    ):
         """Obtain the maximum average charge along the path
         Args:
             migration_hop (MigrationHop): MigrationPath object that represents a given hop
@@ -599,7 +606,7 @@ def _get_avg_chg_at_max(self, migration_hop, radius=None, chg_along_path=False,
             return max(avg_chg), avg_chg
         return max(avg_chg)
 
-    def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None):
+    def _get_chg_between_sites_tube(self, migration_hop: MigrationHop, mask_file_seedname: str | None = None) -> float:
         """
         Calculate the amount of charge that a migrating ion has to move through in order to complete a hop
         Args:
@@ -656,7 +663,7 @@ def _get_chg_between_sites_tube(self, migration_hop, mask_file_seedname=None):
             / self.potential_field.structure.volume
         )
 
-    def populate_edges_with_chg_density_info(self, tube_radius=1):
+    def populate_edges_with_chg_density_info(self, tube_radius: float = 1) -> None:
         """
         Args:
             tube_radius: Tube radius.
@@ -683,7 +690,7 @@ def populate_edges_with_chg_density_info(self, tube_radius=1):
             )
             self.add_data_to_similar_edges(k, {"max_avg_chg": max_chg})
 
-    def get_least_chg_path(self):
+    def get_least_chg_path(self) -> list:
         """
         obtain an intercolating pathway through the material that has the least amount of charge
         Returns:
@@ -701,7 +708,7 @@ def get_least_chg_path(self):
                 min_path = path
         return min_path
 
-    def get_summary_dict(self, add_keys: list[str] | None = None):
+    def get_summary_dict(self, add_keys: list[str] | None = None) -> dict:
         """Dictionary format, for saving to database."""
         a_keys = ["max_avg_chg", "chg_total"]
         if add_keys is not None:
@@ -710,7 +717,7 @@ def get_summary_dict(self, add_keys: list[str] | None = None):
 
 
 # Utility functions
-def get_only_sites_from_structure(structure: Structure, migrating_specie: str) -> Structure:
+def get_only_sites_from_structure(structure: Structure, migrating_specie: SpeciesLike) -> Structure:
     """
     Get a copy of the structure with only the migrating sites.
 
@@ -729,7 +736,7 @@ def get_only_sites_from_structure(structure: Structure, migrating_specie: str) -
     return Structure.from_sites(migrating_ion_sites)
 
 
-def _shift_grid(vv):
+def _shift_grid(vv: np.ndarray) -> np.ndarray:
     """
     Move the grid points by half a step so that they sit in the center
     Args:
@@ -844,7 +851,7 @@ def order_path(hop_list: list[dict], start_u: int | str) -> list[dict]:
 # Utility Functions for comparing UC and SC hops
 
 
-def almost(a, b):
+def almost(a, b) -> bool:  # noqa: ANN001
     """Return true if the values are almost equal."""
     SMALL_VAL = 1e-4
     try:
@@ -855,7 +862,7 @@ def almost(a, b):
         raise NotImplementedError
 
 
-def check_uc_hop(sc_hop, uc_hop):
+def check_uc_hop(sc_hop: MigrationHop, uc_hop: MigrationHop) -> tuple | None:
     """
     See if hop in the 2X2X2 supercell and a unit cell hop
     are equivalent under lattice translation.
@@ -901,7 +908,7 @@ def check_uc_hop(sc_hop, uc_hop):
     return None
 
 
-def map_hop_sc2uc(sc_hop: MigrationHop, mg: MigrationGraph):
+def map_hop_sc2uc(sc_hop: MigrationHop, mg: MigrationGraph) -> dict:
     """
     Map a given hop in the SC onto the UC.
 
diff --git a/pymatgen/analysis/diffusion/neb/io.py b/pymatgen/analysis/diffusion/neb/io.py
index 5550892..26b4205 100644
--- a/pymatgen/analysis/diffusion/neb/io.py
+++ b/pymatgen/analysis/diffusion/neb/io.py
@@ -16,7 +16,7 @@
 class MVLCINEBEndPointSet(MITRelaxSet):
     """Class for writing NEB end points relaxation inputs."""
 
-    def __init__(self, structure, **kwargs):
+    def __init__(self, structure: Structure, **kwargs) -> None:
         r"""
         Args:
             structure: Structure
@@ -48,7 +48,7 @@ class MVLCINEBSet(MITNEBSet):
     http://theory.cm.utexas.edu/vtsttools/.
     """
 
-    def __init__(self, structures, **kwargs):
+    def __init__(self, structures: list[Structure], **kwargs) -> None:
         r"""
         Args:
             structures: Input structures.
@@ -82,7 +82,7 @@ def __init__(self, structures, **kwargs):
         super().__init__(structures, **kwargs)
 
 
-def get_endpoints_from_index(structure, site_indices):
+def get_endpoints_from_index(structure: Structure, site_indices: list[int]) -> list[Structure]:
     """
     This class reads in one perfect structure and the two endpoint structures
     are generated using site_indices.
@@ -121,7 +121,7 @@ def get_endpoints_from_index(structure, site_indices):
     return [s_0, s_1]
 
 
-def get_endpoint_dist(ep_0, ep_1):
+def get_endpoint_dist(ep_0: Structure, ep_1: Structure) -> list[float]:
     """
     Calculate a list of site distances between two endpoints, assuming periodic
     boundary conditions.
diff --git a/pymatgen/analysis/diffusion/neb/pathfinder.py b/pymatgen/analysis/diffusion/neb/pathfinder.py
index c6a64ab..27028dc 100644
--- a/pymatgen/analysis/diffusion/neb/pathfinder.py
+++ b/pymatgen/analysis/diffusion/neb/pathfinder.py
@@ -25,7 +25,11 @@
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 
 if TYPE_CHECKING:
+    from collections.abc import Sequence
+
+    from pymatgen.io.vasp.outputs import Chgcar
     from pymatgen.symmetry.structure import SymmetrizedStructure
+    from pymatgen.util.typing import PathLike, SpeciesLike
 
 logger = logging.getLogger(__name__)
 
@@ -39,7 +43,7 @@ class IDPPSolver:
     Smidstrup et al., J. Chem. Phys. 140, 214106 (2014).
     """
 
-    def __init__(self, structures):
+    def __init__(self, structures: list[Structure]) -> None:
         """
         Initialization.
 
@@ -50,7 +54,7 @@ def __init__(self, structures):
         latt = structures[0].lattice
         natoms = structures[0].num_sites
         nimages = len(structures) - 2
-        target_dists = []
+        _target_dists = []
 
         # Initial guess of the path (in Cartesian coordinates) used in the IDPP
         # algo.
@@ -64,9 +68,9 @@ def __init__(self, structures):
                 structures[-1].distance_matrix - structures[0].distance_matrix
             )
 
-            target_dists.append(dist)
+            _target_dists.append(dist)
 
-        target_dists = np.array(target_dists)
+        target_dists = np.array(_target_dists)
 
         # Set of translational vector matrices (anti-symmetric) for the images.
         translations = np.zeros((nimages, natoms, natoms, 3), dtype=np.float64)
@@ -98,14 +102,14 @@ def __init__(self, structures):
 
     def run(
         self,
-        maxiter=1000,
-        tol=1e-5,
-        gtol=1e-3,
-        step_size=0.05,
-        max_disp=0.05,
-        spring_const=5.0,
-        species=None,
-    ):
+        maxiter: int = 1000,
+        tol: float = 1e-5,
+        gtol: float = 1e-3,
+        step_size: float = 0.05,
+        max_disp: float = 0.05,
+        spring_const: float = 5.0,
+        species: list[SpeciesLike] | None = None,
+    ) -> list[Structure]:
         """
         Perform iterative minimization of the set of objective functions in an
         NEB-like manner. In each iteration, the total force matrix for each
@@ -192,11 +196,11 @@ def run(
     @classmethod
     def from_endpoints(
         cls,
-        endpoints,
+        endpoints: list[Structure],
         nimages: int = 5,
         sort_tol: float = 1.0,
         interpolate_lattices: bool = False,
-    ):
+    ) -> IDPPSolver:
         """
         A class method that starts with end-point structures instead. The
         initial guess for the IDPP algo is then constructed using linear
@@ -234,7 +238,7 @@ def from_endpoints(
 
         return IDPPSolver(images)
 
-    def _get_funcs_and_forces(self, x):
+    def _get_funcs_and_forces(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
         """
         Calculate the set of objective functions as well as their gradients,
         i.e. "effective true forces".
@@ -264,7 +268,7 @@ def _get_funcs_and_forces(self, x):
         return 0.5 * np.array(funcs), -2 * np.array(funcs_prime)
 
     @staticmethod
-    def get_unit_vector(vec):
+    def get_unit_vector(vec: np.ndarray) -> np.ndarray:
         """
         Calculate the unit vector of a vector.
 
@@ -273,7 +277,7 @@ def get_unit_vector(vec):
         """
         return vec / np.sqrt(np.sum(vec**2))
 
-    def _get_total_forces(self, x, true_forces, spring_const):
+    def _get_total_forces(self, x: np.ndarray, true_forces: np.ndarray, spring_const: float) -> np.ndarray:
         """
         Calculate the total force on each image structure, which is equal to
         the spring force along the tangent + true force perpendicular to the
@@ -312,7 +316,7 @@ def __init__(
         symm_structure: SymmetrizedStructure,
         host_symm_struct: SymmetrizedStructure = None,
         symprec: float = 0.001,
-    ):
+    ) -> None:
         """
         Args:
             isite: Initial site
@@ -324,8 +328,8 @@ def __init__(
         """
         self.isite = isite
         self.esite = esite
-        self.iindex = None
-        self.eindex = None
+        self.iindex = -1
+        self.eindex = -1
         self.symm_structure = symm_structure
         self.host_symm_struct = host_symm_struct
         self.symprec = symprec
@@ -343,7 +347,7 @@ def __init__(
                 self.eindex = i
 
         # if no index was identified then loop over each site until something is found
-        if self.iindex is None:
+        if self.iindex is -1:
             for i, sites in enumerate(self.symm_structure.equivalent_sites):
                 for itr_site in sites:
                     if sg.are_symmetrically_equivalent([isite], [itr_site]):
@@ -352,7 +356,7 @@ def __init__(
                 else:
                     continue
                 break
-        if self.eindex is None:
+        if self.eindex is -1:
             for i, sites in enumerate(self.symm_structure.equivalent_sites):
                 for itr_site in sites:
                     if sg.are_symmetrically_equivalent([esite], [itr_site]):
@@ -362,12 +366,12 @@ def __init__(
                     continue
                 break
 
-        if self.iindex is None:
+        if self.iindex is -1:
             raise RuntimeError(f"No symmetrically equivalent site was found for {isite}")
-        if self.eindex is None:
+        if self.eindex is -1:
             raise RuntimeError(f"No symmetrically equivalent site was found for {esite}")
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         ifc = self.isite.frac_coords
         efc = self.esite.frac_coords
         wyk_symbols = self.symm_structure.wyckoff_symbols
@@ -381,20 +385,20 @@ def __repr__(self):
         )
 
     @property
-    def length(self):
+    def length(self) -> float:
         """
         Returns:
             (float) Length of migration path.
         """
-        return np.linalg.norm(self.isite.coords - self.esite.coords)
+        return float(np.linalg.norm(self.isite.coords - self.esite.coords))
 
-    def __hash__(self):
+    def __hash__(self) -> int:
         return self.iindex + self.eindex
 
-    def __str__(self):
+    def __str__(self) -> str:
         return self.__repr__()
 
-    def __eq__(self, other):
+    def __eq__(self, other) -> bool:  # noqa: ANN001
         if self.symm_structure != other.symm_structure:
             return False
 
@@ -407,7 +411,9 @@ def __eq__(self, other):
             self.symprec,
         )
 
-    def get_structures(self, nimages=5, vac_mode=True, idpp=False, **idpp_kwargs):
+    def get_structures(
+        self, nimages: int = 5, vac_mode: bool = True, idpp: bool = False, **idpp_kwargs
+    ) -> list[Structure]:
         r"""
         Generate structures for NEB calculation.
 
@@ -448,7 +454,7 @@ def get_structures(self, nimages=5, vac_mode=True, idpp=False, **idpp_kwargs):
 
         return structures
 
-    def _split_migrating_and_other_sites(self, vac_mode):
+    def _split_migrating_and_other_sites(self, vac_mode: bool) -> tuple[list[Site], list[Site]]:
         migrating_specie_sites = []
         other_sites = []
         for site in self.symm_structure.sites:
@@ -510,7 +516,7 @@ def get_sc_structures(
         )
         return start_struct, end_struct, base_sc
 
-    def write_path(self, fname, **kwargs):
+    def write_path(self, fname: PathLike, **kwargs) -> None:
         r"""
         Write the path to a file for easy viewing.
 
@@ -535,12 +541,12 @@ class DistinctPathFinder:
 
     def __init__(
         self,
-        structure,
-        migrating_specie,
-        max_path_length=None,
-        symprec=0.1,
-        perc_mode=">1d",
-    ):
+        structure: Structure,
+        migrating_specie: SpeciesLike,
+        max_path_length: float | None = None,
+        symprec: float = 0.1,
+        perc_mode: str = ">1d",
+    ) -> None:
         """
         Args:
             structure: Input structure that contains all sites.
@@ -598,7 +604,7 @@ def __init__(
         else:
             self.max_path_length = max_path_length
 
-    def get_paths(self):
+    def get_paths(self) -> list[MigrationHop]:
         """
         Returns:
             [MigrationHop] All distinct migration paths.
@@ -614,7 +620,7 @@ def get_paths(self):
 
         return sorted(paths, key=lambda p: p.length)
 
-    def write_all_paths(self, fname, nimages=5, **kwargs):
+    def write_all_paths(self, fname: PathLike, nimages: int = 5, **kwargs) -> None:
         r"""
         Write a file containing all paths, using hydrogen as a placeholder for
         the images. H is chosen as it is the smallest atom. This is extremely
@@ -650,7 +656,15 @@ class NEBPathfinder:
         Ceder, The Journal of Chemical Physics 145 (7), 074112
     """
 
-    def __init__(self, start_struct, end_struct, relax_sites, v, n_images=20, mid_struct=None):
+    def __init__(
+        self,
+        start_struct: Structure,
+        end_struct: Structure,
+        relax_sites: list[int],
+        v: np.ndarray,
+        n_images: int = 20,
+        mid_struct: Structure = None,
+    ) -> None:
         """
         Args:
             start_struct: Starting structure
@@ -668,10 +682,9 @@ def __init__(self, start_struct, end_struct, relax_sites, v, n_images=20, mid_st
         self.__relax_sites = relax_sites
         self.__v = v
         self.__n_images = n_images
-        self.__images = None
         self.interpolate()
 
-    def interpolate(self):
+    def interpolate(self) -> None:
         """
         Finds a set of n_images from self.s1 to self.s2, where all sites except
         for the ones given in relax_sites, the interpolation is linear (as in
@@ -717,14 +730,14 @@ def interpolate(self):
         self.__images = images
 
     @property
-    def images(self):
+    def images(self) -> list[Structure]:
         """
         Returns a list of structures interpolating between the start and
         endpoint structures.
         """
         return self.__images
 
-    def plot_images(self, outfile):
+    def plot_images(self, outfile: PathLike) -> None:
         """
         Generates a POSCAR with the calculated diffusion path with respect to the first
         endpoint.
@@ -732,6 +745,7 @@ def plot_images(self, outfile):
         Args:
             outfile: Output file for the POSCAR.
         """
+        assert self.__images is not None, "Images have not been calculated yet. Run interpolate() first."
         sum_struct = self.__images[0].sites
         for image in self.__images:
             for site_i in self.__relax_sites:
@@ -750,17 +764,17 @@ def plot_images(self, outfile):
 
     @staticmethod
     def string_relax(
-        start,
-        end,
-        V,
-        n_images=25,
-        dr=None,
-        h=3.0,
-        k=0.17,
-        min_iter=100,
-        max_iter=10000,
-        max_tol=5e-6,
-    ):
+        start: np.ndarray,
+        end: np.ndarray,
+        V: np.ndarray,
+        n_images: int = 25,
+        dr: np.ndarray | Sequence | None = None,
+        h: float = 3.0,
+        k: float = 0.17,
+        min_iter: int = 100,
+        max_iter: int = 10000,
+        max_tol: float = 5e-6,
+    ) -> np.ndarray:
         """
         Implements path relaxation via the elastic band method. In general, the
         method is to define a path by a set of points (images) connected with
@@ -883,7 +897,7 @@ def string_relax(
         return s
 
     @staticmethod
-    def __f2d(frac_coords, v):
+    def __f2d(frac_coords: np.ndarray, v: np.ndarray) -> np.ndarray:
         """
         Converts fractional coordinates to discrete coordinates with respect to
         the grid size of v.
@@ -892,7 +906,7 @@ def __f2d(frac_coords, v):
         return (np.array(frac_coords) * np.array(v.shape)).astype(int)
 
     @staticmethod
-    def __d2f(disc_coords, v):
+    def __d2f(disc_coords: np.ndarray, v: np.ndarray) -> np.ndarray:
         """
         Converts a point given in discrete coordinates with respect to the
         grid in v to fractional coordinates.
@@ -913,7 +927,7 @@ class StaticPotential:
     function to normalize the potential from 0 to 1 (recommended).
     """
 
-    def __init__(self, struct, pot):
+    def __init__(self, struct: Structure, pot: np.ndarray) -> None:
         """
         :param struct: atomic structure of the potential
         :param pot: volumentric data to be used as a potential
@@ -921,16 +935,16 @@ def __init__(self, struct, pot):
         self.__v = pot
         self.__s = struct
 
-    def get_v(self):
+    def get_v(self) -> np.ndarray:
         """Returns the potential."""
         return self.__v
 
-    def normalize(self):
+    def normalize(self) -> None:
         """Sets the potential range 0 to 1."""
         self.__v = self.__v - np.amin(self.__v)
         self.__v = self.__v / np.amax(self.__v)
 
-    def rescale_field(self, new_dim):
+    def rescale_field(self, new_dim: tuple) -> None:
         """
         Changes the discretization of the potential field by linear
         interpolation. This is necessary if the potential field
@@ -956,7 +970,7 @@ def rescale_field(self, new_dim):
         )
         self.__v = v_ngrid
 
-    def gaussian_smear(self, r):
+    def gaussian_smear(self, r: float) -> None:
         """
         Applies an isotropic Gaussian smear of width (standard deviation) r to
         the potential field. This is necessary to avoid finding paths through
@@ -1007,7 +1021,7 @@ def gaussian_smear(self, r):
 class ChgcarPotential(StaticPotential):
     """Implements a potential field based on the charge density output from VASP."""
 
-    def __init__(self, chgcar, smear=False, normalize=True):
+    def __init__(self, chgcar: Chgcar, smear: bool = False, normalize: bool = True) -> None:
         """
         :param chgcar: Chgcar object based on a VASP run of the structure of
             interest (Chgcar.from_file("CHGCAR"))
@@ -1032,7 +1046,7 @@ class FreeVolumePotential(StaticPotential):
     is lower at points farther away from any atoms in the structure.
     """
 
-    def __init__(self, struct, dim, smear=False, normalize=True):
+    def __init__(self, struct: Structure, dim: tuple, smear: bool = False, normalize: bool = True) -> None:
         """
         :param struct: Unit cell on which to base the potential
         :param dim: Grid size for the potential
@@ -1050,7 +1064,7 @@ def __init__(self, struct, dim, smear=False, normalize=True):
             self.normalize()
 
     @staticmethod
-    def __add_gaussians(s, dim, r=1.5):
+    def __add_gaussians(s: Structure, dim: Sequence, r: float = 1.5) -> np.ndarray:
         gauss_dist = np.zeros(dim)
         for a_d in np.arange(0, dim[0], 1):
             for b_d in np.arange(0, dim[1], 1):
@@ -1065,7 +1079,9 @@ def __add_gaussians(s, dim, r=1.5):
 class MixedPotential(StaticPotential):
     """Implements a potential that is a weighted sum of some other potentials."""
 
-    def __init__(self, potentials, coefficients, smear=False, normalize=True):
+    def __init__(
+        self, potentials: list[StaticPotential], coefficients: list[float], smear: bool = False, normalize: bool = True
+    ) -> None:
         """
         Args:
             potentials: List of objects extending the StaticPotential superclass
diff --git a/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py b/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py
index 4bafd9e..122787c 100644
--- a/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py
+++ b/pymatgen/analysis/diffusion/neb/periodic_dijkstra.py
@@ -18,6 +18,8 @@
 import numpy as np
 
 if TYPE_CHECKING:
+    from collections.abc import Generator
+
     from networkx.classes.graph import Graph
 
     from pymatgen.analysis.graphs import StructureGraph
@@ -36,7 +38,7 @@ def _get_adjacency_with_images(G: Graph) -> dict:
         dict: Nested dictionary with [start][end][edge_key][data_field]
     """
 
-    def copy_dict(d):
+    def copy_dict(d: dict) -> dict:
         # recursively copies the dictionary to resolve the fact that
         # two entries in the dictionary can point to the same mutable object
         # eg. changing p_graph[v][u][0]["to_jimage"] also changes
@@ -67,7 +69,7 @@ def periodic_dijkstra(
     weight: str = "weight",
     max_image: int = 2,
     target_reached: Callable = lambda idx, jimage: False,
-):
+) -> tuple | dict:
     """
     Find the lowest cost pathway from a source point in the periodic graph.
     Since the search can move many cells away without finding the target
@@ -127,7 +129,7 @@ def periodic_dijkstra_on_sgraph(
     weight: str = "weight",
     max_image: int = 1,
     target_reached: Callable = lambda idx, jimage: False,
-):
+) -> tuple:
     """
     Find the lowest cost pathway from a source point in the periodic graph.
     Since the search can move many cells away without finding the target
@@ -157,7 +159,7 @@ def periodic_dijkstra_on_sgraph(
     return best_ans, path_parent
 
 
-def get_optimal_pathway_rev(path_parent: dict, leaf_node: tuple):
+def get_optimal_pathway_rev(path_parent: dict, leaf_node: tuple) -> Generator:
     """Follow a leaf node all the way up to source."""
     cur = leaf_node
     while cur in path_parent:
diff --git a/pymatgen/analysis/diffusion/neb/tests/test_full_path_mapper.py b/pymatgen/analysis/diffusion/neb/tests/test_full_path_mapper.py
index d0c5d38..afe6a72 100644
--- a/pymatgen/analysis/diffusion/neb/tests/test_full_path_mapper.py
+++ b/pymatgen/analysis/diffusion/neb/tests/test_full_path_mapper.py
@@ -27,25 +27,25 @@
 
 
 class MigrationGraphSimpleTest(unittest.TestCase):
-    def setUp(self):
+    def setUp(self) -> None:
         struct = Structure.from_file(f"{dir_path}/full_path_files/MnO2_full_Li.vasp")
         self.fpm = MigrationGraph.with_distance(structure=struct, migrating_specie="Li", max_distance=4)
 
-    def test_get_pos_and_migration_hop(self):
+    def test_get_pos_and_migration_hop(self) -> None:
         """
         Make sure that we can populate the graph with MigrationHop Objects
         """
         self.fpm._get_pos_and_migration_hop(0, 1, 1)
         self.assertAlmostEqual(self.fpm.m_graph.graph[0][1][1]["hop"].length, 3.571248, 4)
 
-    def test_get_summary_dict(self):
+    def test_get_summary_dict(self) -> None:
         summary_dict = self.fpm.get_summary_dict()
         assert "hop_label" in summary_dict["hops"][0]
         assert "hop_label" in summary_dict["unique_hops"][0]
 
 
 class MigrationGraphFromEntriesTest(unittest.TestCase):
-    def setUp(self):
+    def setUp(self) -> None:
         self.test_ents_MOF = loadfn(f"{dir_path}/full_path_files/Mn6O5F7_cat_migration.json")
         self.aeccar_MOF = Chgcar.from_file(f"{dir_path}/full_path_files/AECCAR_Mn6O5F7.vasp")
         self.li_ent = loadfn(f"{dir_path}/full_path_files/li_ent.json")["li_ent"]
@@ -55,7 +55,7 @@ def setUp(self):
             migrating_ion_entry=self.li_ent,
         )[0]
 
-    def test_m_graph_from_entries_failed(self):
+    def test_m_graph_from_entries_failed(self) -> None:
         # only base
         s_list = MigrationGraph.get_structure_from_entries(
             entries=[self.test_ents_MOF["ent_base"]],
@@ -69,14 +69,14 @@ def test_m_graph_from_entries_failed(self):
         )
         assert len(s_list) == 0
 
-    def test_m_graph_construction(self):
+    def test_m_graph_construction(self) -> None:
         assert self.full_struct.composition["Li"] == 8
         mg = MigrationGraph.with_distance(self.full_struct, migrating_specie="Li", max_distance=4.0)
         assert len(mg.m_graph.structure) == 8
 
 
 class MigrationGraphComplexTest(unittest.TestCase):
-    def setUp(self):
+    def setUp(self) -> None:
         struct = Structure.from_file(f"{dir_path}/full_path_files/MnO2_full_Li.vasp")
         self.fpm_li = MigrationGraph.with_distance(structure=struct, migrating_specie="Li", max_distance=4)
 
@@ -85,7 +85,7 @@ def setUp(self):
         struct = Structure.from_file(f"{dir_path}/full_path_files/Mg_2atom.vasp")
         self.fpm_mg = MigrationGraph.with_distance(structure=struct, migrating_specie="Mg", max_distance=2)
 
-    def test_group_and_label_hops(self):
+    def test_group_and_label_hops(self) -> None:
         """
         Check that the set of end points in a group of similarly labeled hops are all
         the same.
@@ -107,7 +107,7 @@ def test_group_and_label_hops(self):
             for end_point_labels in sub_set:
                 assert sorted(end_point_labels) == sorted(sub_set[0])
 
-    def test_unique_hops_dict(self):
+    def test_unique_hops_dict(self) -> None:
         """
         Check that the unique hops are inequivalent
         """
@@ -117,7 +117,7 @@ def test_unique_hops_dict(self):
         for migration_hop in all_pairs:
             assert migration_hop[0]["hop"] != migration_hop[1]["hop"]
 
-    def test_add_data_to_similar_edges(self):
+    def test_add_data_to_similar_edges(self) -> None:
         # passing normal data
         self.fpm_li.add_data_to_similar_edges(0, {"key0": "data"})
         for _u, _v, d in self.fpm_li.m_graph.graph.edges(data=True):
@@ -142,7 +142,7 @@ def test_add_data_to_similar_edges(self):
             if d["hop_label"] == 2:
                 assert d["key2"] == [3, 2, 1]
 
-    def test_assign_cost_to_graph(self):
+    def test_assign_cost_to_graph(self) -> None:
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
         for _u, _v, d in self.fpm_li.m_graph.graph.edges(data=True):
             self.assertAlmostEqual(d["cost"], d["hop_distance"], 4)
@@ -151,7 +151,7 @@ def test_assign_cost_to_graph(self):
         for _u, _v, d in self.fpm_li.m_graph.graph.edges(data=True):
             self.assertAlmostEqual(d["cost"], d["hop_distance"] ** 2, 4)
 
-    def test_periodic_dijkstra(self):
+    def test_periodic_dijkstra(self) -> None:
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
 
         # test the connection graph
@@ -165,7 +165,7 @@ def test_periodic_dijkstra(self):
                     opposite_connections = [d2_["to_jimage"] for k2_, d2_ in conn_dict[v][u].items()]
                     assert neg_image in opposite_connections
 
-    def test_get_path(self):
+    def test_get_path(self) -> None:
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
         paths = [*self.fpm_li.get_path(flip_hops=False)]
         p_strings = {"->".join(map(str, get_hop_site_sequence(ipath, start_u=u))) for u, ipath in paths}
@@ -182,12 +182,12 @@ def test_get_path(self):
         p_strings = {"->".join(map(str, get_hop_site_sequence(ipath, start_u=u))) for u, ipath in paths}
         assert "1->0->1" in p_strings
 
-    def test_get_key_in_path(self):
+    def test_get_key_in_path(self) -> None:
         self.fpm_li.assign_cost_to_graph()  # use 'hop_distance'
         paths = [*self.fpm_li.get_path(flip_hops=False)]
         hop_seq_info = [get_hop_site_sequence(ipath, start_u=u, key="hop_distance") for u, ipath in paths]
 
-        hop_distances = {}
+        hop_distances: dict = {}
         for u, ipath in paths:
             hop_distances[u] = []
             for hop in ipath:
@@ -198,13 +198,13 @@ def test_get_key_in_path(self):
             assert distances == hop_seq_info[u][1]
             print(distances, hop_seq_info[u][1])
 
-    def test_not_matching_first(self):
+    def test_not_matching_first(self) -> None:
         structure = Structure.from_file(f"{dir_path}/pathfinder_files/Li6MnO4.json")
         fpm_lmo = MigrationGraph.with_distance(structure, "Li", max_distance=4)
         for _u, _v, d in fpm_lmo.m_graph.graph.edges(data=True):
             assert d["hop"].eindex in {0, 1}
 
-    def test_order_path(self):
+    def test_order_path(self) -> None:
         # add list data to migration graph - to test if list data is flipped
         for n, hop_d in self.fpm_li.unique_hops.items():
             data = {"data": [hop_d["iindex"], hop_d["eindex"]]}
@@ -227,7 +227,7 @@ def test_order_path(self):
 
 
 class ChargeBarrierGraphTest(unittest.TestCase):
-    def setUp(self):
+    def setUp(self) -> None:
         self.full_sites_MOF = loadfn(f"{dir_path}/full_path_files/LixMn6O5F7_full_sites.json")
         self.aeccar_MOF = Chgcar.from_file(f"{dir_path}/full_path_files/AECCAR_Mn6O5F7.vasp")
         self.cbg = ChargeBarrierGraph.with_distance(
@@ -239,7 +239,7 @@ def setUp(self):
         )
         self.cbg._tube_radius = 10000
 
-    def test_integration(self):
+    def test_integration(self) -> None:
         """
         Sanity check: for a long enough diagonally hop, if we turn the radius of the tube way up, it should cover the entire unit cell
         """
@@ -273,7 +273,7 @@ def test_integration(self):
             3,
         )
 
-    def test_populate_edges_with_chg_density_info(self):
+    def test_populate_edges_with_chg_density_info(self) -> None:
         """
         Test that all of the sites with similar lengths have similar charge densities,
         this will not always be true, but it valid in this Mn6O5F7
@@ -289,7 +289,7 @@ def test_populate_edges_with_chg_density_info(self):
             if 1.05 > length / prv[0] > 0.95:
                 self.assertAlmostEqual(chg, prv[1], 3)
 
-    def test_get_summary_dict(self):
+    def test_get_summary_dict(self) -> None:
         summary_dict = self.cbg.get_summary_dict()
         assert "chg_total", summary_dict["hops"][0]  # noqa: PLW0129
         assert "chg_total", summary_dict["unique_hops"][0]  # noqa: PLW0129
diff --git a/pymatgen/analysis/diffusion/neb/tests/test_io.py b/pymatgen/analysis/diffusion/neb/tests/test_io.py
index 4318c24..c07b5a6 100644
--- a/pymatgen/analysis/diffusion/neb/tests/test_io.py
+++ b/pymatgen/analysis/diffusion/neb/tests/test_io.py
@@ -2,6 +2,7 @@
 
 import os
 import unittest
+from typing import TYPE_CHECKING
 
 from pymatgen.analysis.diffusion.neb.io import (
     MVLCINEBEndPointSet,
@@ -11,12 +12,15 @@
 )
 from pymatgen.core import Structure
 
+if TYPE_CHECKING:
+    from pymatgen.util.typing import PathLike
+
 __author__ = "hat003"
 
 test_dir = os.path.join(os.path.abspath(os.path.dirname(__file__)))
 
 
-def get_path(path_str, dirname="./"):
+def get_path(path_str: PathLike, dirname: PathLike = "./") -> str:
     cwd = os.path.abspath(os.path.dirname(__file__))
     return os.path.join(cwd, dirname, path_str)
 
@@ -24,7 +28,7 @@ def get_path(path_str, dirname="./"):
 class MVLCINEBEndPointSetTest(unittest.TestCase):
     endpoint = Structure.from_file(get_path("POSCAR0", dirname="io_files"))
 
-    def test_incar(self):
+    def test_incar(self) -> None:
         m = MVLCINEBEndPointSet(self.endpoint)
         incar_string = m.incar.get_str(sort_keys=True, pretty=True)
         incar_expect = """ALGO     =  Fast
@@ -51,7 +55,7 @@ def test_incar(self):
 SIGMA    =  0.05"""
         assert incar_string == incar_expect
 
-    def test_incar_user_setting(self):
+    def test_incar_user_setting(self) -> None:
         user_incar_settings = {
             "ALGO": "Normal",
             "EDIFFG": -0.05,
@@ -91,7 +95,7 @@ def test_incar_user_setting(self):
 class MVLCINEBSetTest(unittest.TestCase):
     structures = [Structure.from_file(get_path("POSCAR" + str(i), dirname="io_files")) for i in range(3)]
 
-    def test_incar(self):
+    def test_incar(self) -> None:
         m = MVLCINEBSet(self.structures)
 
         incar_string = m.incar.get_str(sort_keys=True)
@@ -125,7 +129,7 @@ def test_incar(self):
 SPRING = -5"""
         assert incar_string.strip() == incar_expect.strip()
 
-    def test_incar_user_setting(self):
+    def test_incar_user_setting(self) -> None:
         user_incar_settings = {"IOPT": 3, "EDIFFG": -0.05, "NPAR": 4, "ISIF": 3}
         m = MVLCINEBSet(self.structures, user_incar_settings=user_incar_settings)
         incar_string = m.incar.get_str(sort_keys=True, pretty=True)
@@ -168,7 +172,7 @@ class UtilityTest(unittest.TestCase):
 
     structure = Structure.from_file(get_path("POSCAR", dirname="io_files"))
 
-    def test_get_endpoints_from_index(self):
+    def test_get_endpoints_from_index(self) -> None:
         endpoints = get_endpoints_from_index(structure=self.structure, site_indices=[0, 1])
         ep_0 = endpoints[0].as_dict()
         ep_1 = endpoints[1].as_dict()
@@ -178,7 +182,7 @@ def test_get_endpoints_from_index(self):
         assert ep_0 == ep_0_expect
         assert ep_1 == ep_1_expect
 
-    def test_get_endpoint_dist(self):
+    def test_get_endpoint_dist(self) -> None:
         ep0 = Structure.from_file(get_path("POSCAR_ep0", dirname="io_files"))
         ep1 = Structure.from_file(get_path("POSCAR_ep1", dirname="io_files"))
         distances = get_endpoint_dist(ep0, ep1)
diff --git a/pymatgen/analysis/diffusion/neb/tests/test_pathfinder.py b/pymatgen/analysis/diffusion/neb/tests/test_pathfinder.py
index 3eeeb0b..03baab6 100644
--- a/pymatgen/analysis/diffusion/neb/tests/test_pathfinder.py
+++ b/pymatgen/analysis/diffusion/neb/tests/test_pathfinder.py
@@ -5,6 +5,7 @@
 import glob
 import os
 import unittest
+from typing import TYPE_CHECKING
 
 import numpy as np
 
@@ -15,12 +16,15 @@
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from pymatgen.util.testing import PymatgenTest
 
+if TYPE_CHECKING:
+    from pymatgen.util.typing import PathLike
+
 __author__ = "Iek-Heng Chu"
 __version__ = "1.0"
 __date__ = "March 14, 2017"
 
 
-def get_path(path_str, dirname="./"):
+def get_path(path_str: PathLike, dirname: PathLike = "./") -> str:
     cwd = os.path.abspath(os.path.dirname(__file__))
     return os.path.join(cwd, dirname, path_str)
 
@@ -29,7 +33,7 @@ class IDPPSolverTest(unittest.TestCase):
     init_struct = Structure.from_file(get_path("CONTCAR-0", dirname="pathfinder_files"))
     final_struct = Structure.from_file(get_path("CONTCAR-1", dirname="pathfinder_files"))
 
-    def test_idpp_from_ep(self):
+    def test_idpp_from_ep(self) -> None:
         obj = IDPPSolver.from_endpoints([self.init_struct, self.final_struct], nimages=3, sort_tol=1.0)
         new_path = obj.run(
             maxiter=5000,
@@ -49,7 +53,7 @@ def test_idpp_from_ep(self):
         assert np.allclose(new_path[3][22].frac_coords, np.array([0.28422885, 0.62568764, 0.98975444]))
         assert np.allclose(new_path[4][47].frac_coords, np.array([0.59767531, 0.12640952, 0.37745006]))
 
-    def test_idpp_from_ep_diff_latt(self):
+    def test_idpp_from_ep_diff_latt(self) -> None:
         # This is the same test as test_idpp_from_ep where we want to interpolate
         # different lattices.
 
@@ -137,7 +141,7 @@ def test_idpp_from_ep_diff_latt(self):
             ),
         )
 
-    def test_idpp(self):
+    def test_idpp(self) -> None:
         images = self.init_struct.interpolate(self.final_struct, nimages=4, autosort_tol=1.0)
         obj = IDPPSolver(images)
         new_path = obj.run(
@@ -158,7 +162,7 @@ def test_idpp(self):
 
 
 class DistinctPathFinderTest(PymatgenTest):
-    def test_get_paths(self):
+    def test_get_paths(self) -> None:
         s = self.get_structure("LiFePO4")
         # Only one path in LiFePO4 with 4 A.
         p = DistinctPathFinder(s, "Li", max_path_length=4)
@@ -203,7 +207,7 @@ def test_get_paths(self):
         for f in glob.glob("pathfindertest_*.cif"):
             os.remove(f)
 
-    def test_max_path_length(self):
+    def test_max_path_length(self) -> None:
         s = Structure.from_file(get_path("LYPS.cif", dirname="pathfinder_files"))
         dp1 = DistinctPathFinder(s, "Li", perc_mode="1d")
         self.assertAlmostEqual(dp1.max_path_length, 4.11375354207, 7)
@@ -212,7 +216,7 @@ def test_max_path_length(self):
 
 
 class MigrationHopTest(PymatgenTest):
-    def setUp(self):
+    def setUp(self) -> None:
         self.lifepo = self.get_structure("LiFePO4")
         m_graph = MigrationGraph.with_distance(self.lifepo, max_distance=4.0, migrating_specie="Li")
         gen = iter(m_graph.m_graph.graph.edges(data=True))
@@ -223,11 +227,11 @@ def setUp(self):
         symm_structure = a.get_symmetrized_structure()
         self.m_hop = MigrationHop(i_site, e_site, symm_structure)
 
-    def test_msonable(self):
+    def test_msonable(self) -> None:
         hop_dict = self.m_hop.as_dict()
         assert isinstance(hop_dict, dict)
 
-    def test_get_start_end_structs_from_hop(self):
+    def test_get_start_end_structs_from_hop(self) -> None:
         dist_ref = self.m_hop.length
         base = self.lifepo.copy()
         base.remove_species(["Li"])
@@ -242,7 +246,7 @@ def test_get_start_end_structs_from_hop(self):
         end_site = next(filter(lambda x: x.species_string == "Li", end.sites))
         self.assertAlmostEqual(start_site.distance(end_site), dist_ref, 3)
 
-    def test_get_start_end_structs_from_hop_vac(self):
+    def test_get_start_end_structs_from_hop_vac(self) -> None:
         dist_ref = self.m_hop.length
         start, end, b_sc = get_start_end_structures(
             self.m_hop.isite,
@@ -251,21 +255,21 @@ def test_get_start_end_structs_from_hop_vac(self):
             sc_mat=[[2, 1, 0], [-1, 1, 0], [0, 0, 2]],
             vac_mode=True,
         )
-        s1 = set()
-        s2 = set()
+        _s1 = set()
+        _s2 = set()
         for itr_start, start_site in enumerate(start.sites):
             for itr_end, end_site in enumerate(end.sites):
                 dist_ = start_site.distance(end_site)
                 if dist_ < 1e-5:
-                    s1.add(itr_start)
-                    s2.add(itr_end)
-        s1 = sorted(s1, reverse=True)
-        s2 = sorted(s2, reverse=True)
+                    _s1.add(itr_start)
+                    _s2.add(itr_end)
+        s1 = sorted(_s1, reverse=True)
+        s2 = sorted(_s2, reverse=True)
         start.remove_sites(s1)
         end.remove_sites(s2)
         self.assertAlmostEqual(start.sites[0].distance(end.sites[0]), dist_ref, 3)
 
-    def test_get_sc_structures(self):
+    def test_get_sc_structures(self) -> None:
         dist_ref = self.m_hop.length
         start, end, b_sc = self.m_hop.get_sc_structures(vac_mode=False)
         start_site = next(filter(lambda x: x.species_string == "Li", start.sites))
@@ -274,7 +278,7 @@ def test_get_sc_structures(self):
         assert b_sc.composition == Composition("Fe24 P24 O96")
         self.assertAlmostEqual(start_site.distance(end_site), dist_ref, 3)
 
-    def test_get_sc_structures_vacmode(self):
+    def test_get_sc_structures_vacmode(self) -> None:
         start, end, b_sc = self.m_hop.get_sc_structures(vac_mode=True)
         assert start.composition == end.composition == Composition("Li23 Fe24 P24 O96")
         assert b_sc.composition == Composition("Li24 Fe24 P24 O96")
diff --git a/pymatgen/analysis/diffusion/tests/test_analyzer.py b/pymatgen/analysis/diffusion/tests/test_analyzer.py
index 131e292..4985826 100644
--- a/pymatgen/analysis/diffusion/tests/test_analyzer.py
+++ b/pymatgen/analysis/diffusion/tests/test_analyzer.py
@@ -26,12 +26,12 @@
 
 
 class FuncTest(PymatgenTest):
-    def test_get_conversion_factor(self):
+    def test_get_conversion_factor(self) -> None:
         s = PymatgenTest.get_structure("LiFePO4")
         # large tolerance because scipy constants changed between 0.16.1 and 0.17
         self.assertAlmostEqual(41370704.343540139, get_conversion_factor(s, "Li", 600), delta=20)
 
-    def test_fit_arrhenius(self):
+    def test_fit_arrhenius(self) -> None:
         Ea = 0.5
         k = const.k / const.e
         c = 12
@@ -65,7 +65,7 @@ def test_fit_arrhenius(self):
 
 
 class DiffusionAnalyzerTest(PymatgenTest):
-    def test_init(self):
+    def test_init(self) -> None:
         # Diffusion vasprun.xmls are rather large. We are only going to use a
         # very small preprocessed run for testing. Note that the results are
         # unreliable for short runs.
@@ -236,17 +236,17 @@ def test_init(self):
 
             d.export_msdt("test.csv")
             with open("test.csv") as f:
-                data = []
+                _data = []
                 for row in csv.reader(f):
                     if row:
-                        data.append(row)
-            data.pop(0)
-            data = np.array(data, dtype=np.float64)
-            assert data[:, 1] == pytest.approx(d.msd)
-            assert data[:, -1] == pytest.approx(d.mscd)
+                        _data.append(row)
+            _data.pop(0)
+            data = np.array(_data, dtype=np.float64)
+            assert [row[1] for row in data] == pytest.approx(d.msd)
+            assert [row[-1] for row in data] == pytest.approx(d.mscd)
             os.remove("test.csv")
 
-    def test_init_npt(self):
+    def test_init_npt(self) -> None:
         # Diffusion vasprun.xmls are rather large. We are only going to use a
         # very small preprocessed run for testing. Note that the results are
         # unreliable for short runs.
@@ -452,17 +452,17 @@ def test_init_npt(self):
 
             d.export_msdt("test.csv")
             with open("test.csv") as f:
-                data = []
+                _data = []
                 for row in csv.reader(f):
                     if row:
-                        data.append(row)
-            data.pop(0)
-            data = np.array(data, dtype=np.float64)
-            assert data[:, 1] == pytest.approx(d.msd)
-            assert data[:, -1] == pytest.approx(d.mscd)
+                        _data.append(row)
+            _data.pop(0)
+            data = np.array(_data, dtype=np.float64)
+            assert [row[1] for row in data] == pytest.approx(d.msd)
+            assert [row[-1] for row in data] == pytest.approx(d.mscd)
             os.remove("test.csv")
 
-    def test_from_structure_NPT(self):
+    def test_from_structure_NPT(self) -> None:
         coords1 = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]])
         coords2 = np.array([[0.0, 0.0, 0.0], [0.6, 0.6, 0.6]])
         coords3 = np.array([[0.0, 0.0, 0.0], [0.7, 0.7, 0.7]])
@@ -477,7 +477,7 @@ def test_from_structure_NPT(self):
             structures,
             specie="Li",
             temperature=500.0,
-            time_step=2.0,
+            time_step=2,
             step_skip=1,
             smoothed=None,
         )
diff --git a/pymatgen/analysis/diffusion/tests/test_pathfinder.py b/pymatgen/analysis/diffusion/tests/test_pathfinder.py
index 94c994e..5d8b59c 100644
--- a/pymatgen/analysis/diffusion/tests/test_pathfinder.py
+++ b/pymatgen/analysis/diffusion/tests/test_pathfinder.py
@@ -10,7 +10,7 @@
 
 
 class PathfinderTest(PymatgenTest):
-    def test_mhop_msonable(self):
+    def test_mhop_msonable(self) -> None:
         file_path = os.path.join(module_dir, "migration_graph_spinel_MgMn2O4.json")
         spinel_mg = loadfn(file_path)
         hop = spinel_mg.unique_hops[0]["hop"]
diff --git a/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py b/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py
index 629aafa..55b4a93 100644
--- a/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py
+++ b/pymatgen/analysis/diffusion/utils/edge_data_from_sc.py
@@ -116,7 +116,7 @@ def get_uc_pos(
     return p0, p1, p2
 
 
-def _get_first_close_site(frac_coord, structure, stol=0.1):
+def _get_first_close_site(frac_coord: np.ndarray, structure: Structure, stol: float = 0.1) -> np.ndarray:
     for site in structure.sites:
         dist, image = structure.lattice.get_distance_and_image(frac_coord, site.frac_coords)
         if dist < stol:
@@ -124,7 +124,7 @@ def _get_first_close_site(frac_coord, structure, stol=0.1):
     return frac_coord
 
 
-def mh_eq(mh1, mh2):
+def mh_eq(mh1: MigrationHop, mh2: MigrationHop) -> bool:
     """
     Allow for symmetric matching of MigrationPath objects with variable precession.
 
diff --git a/pymatgen/analysis/diffusion/utils/maggma.py b/pymatgen/analysis/diffusion/utils/maggma.py
index a49b305..9966e59 100644
--- a/pymatgen/analysis/diffusion/utils/maggma.py
+++ b/pymatgen/analysis/diffusion/utils/maggma.py
@@ -32,7 +32,7 @@ def get_entries_from_dbs(
     material_store: MongoStore,
     migrating_ion: str,
     material_id: str,
-):
+) -> tuple:
     """
     Get the entries needed to construct a migration from a database that
     contains topotactically matched structures.
diff --git a/pymatgen/analysis/diffusion/utils/parse_entries.py b/pymatgen/analysis/diffusion/utils/parse_entries.py
index e289ec7..c4e8464 100644
--- a/pymatgen/analysis/diffusion/utils/parse_entries.py
+++ b/pymatgen/analysis/diffusion/utils/parse_entries.py
@@ -96,7 +96,7 @@ def process_entries(
 
     results = []
 
-    def _meta_stable_sites(base_ent, inserted_ent):
+    def _meta_stable_sites(base_ent: ComputedStructureEntry, inserted_ent: ComputedStructureEntry) -> list:
         mapped_struct = get_inserted_on_base(
             base_ent,
             inserted_ent,
@@ -135,7 +135,7 @@ def _meta_stable_sites(base_ent, inserted_ent):
     return sorted(results, key=lambda x: x["inserted"].composition[working_ion], reverse=True)
 
 
-def get_matched_structure_mapping(base: Structure, inserted: Structure, sm: StructureMatcher):
+def get_matched_structure_mapping(base: Structure, inserted: Structure, sm: StructureMatcher) -> tuple | None:
     """
     Get the mapping from the inserted structure onto the base structure,
     assuming that the inserted structure sans the working ion is some kind
@@ -166,7 +166,7 @@ def get_inserted_on_base(
     inserted_ent: ComputedStructureEntry,
     migrating_ion_entry: ComputedEntry,
     sm: StructureMatcher,
-) -> Structure | None:
+) -> Structure:
     """
     For a structured-matched pair of base and inserted entries, map all of the Li
     positions in the inserted entry to positions in the base entry and return a new
@@ -185,7 +185,7 @@ def get_inserted_on_base(
     """
     mapped_result = get_matched_structure_mapping(base_ent.structure, inserted_ent.structure, sm)
     if mapped_result is None:
-        return None
+        raise ValueError("Cannot obtain inverse mapping, consider lowering tolerances in StructureMatcher")
 
     sc_m, total_t = mapped_result
     insertion_energy = get_insertion_energy(base_ent, inserted_ent, migrating_ion_entry)
diff --git a/pymatgen/analysis/diffusion/utils/supercells.py b/pymatgen/analysis/diffusion/utils/supercells.py
index 8b49e79..2ab33e7 100644
--- a/pymatgen/analysis/diffusion/utils/supercells.py
+++ b/pymatgen/analysis/diffusion/utils/supercells.py
@@ -106,7 +106,7 @@ def get_start_end_structures(
         initial structure, final structure, empty structure all in the supercell
     """
 
-    def remove_site_at_pos(structure: Structure, site: PeriodicSite, tol: float):
+    def remove_site_at_pos(structure: Structure, site: PeriodicSite, tol: float) -> Structure:
         new_struct_sites = []
         for isite in structure:
             if not vac_mode or (isite.distance(site) <= tol):
diff --git a/pymatgen/analysis/diffusion/utils/tests/test_edge_data_from_sc.py b/pymatgen/analysis/diffusion/utils/tests/test_edge_data_from_sc.py
index 1c8f4a0..0fa0153 100644
--- a/pymatgen/analysis/diffusion/utils/tests/test_edge_data_from_sc.py
+++ b/pymatgen/analysis/diffusion/utils/tests/test_edge_data_from_sc.py
@@ -30,7 +30,7 @@
 mg_Mg = MigrationGraph.with_distance(structure=mg_uc_full_sites, migrating_specie="Mg", max_distance=4)
 
 
-def test_add_edge_data_from_sc():
+def test_add_edge_data_from_sc() -> None:
     errors = []
 
     test_key = "test_key"
@@ -70,7 +70,7 @@ def test_add_edge_data_from_sc():
     assert not errors, "errors occurred:\n" + "\n".join(errors)
 
 
-def test_get_uc_pos():
+def test_get_uc_pos() -> None:
     errors = []
 
     # set up parameters to initiate get_uc_pos
@@ -114,7 +114,7 @@ def test_get_uc_pos():
     assert not errors, "errors occurred:\n" + "\n".join(errors)
 
 
-def test_get_unique_hop_host():
+def test_get_unique_hop_host() -> None:
     results = get_unique_hop(
         mg_Mg,
         mg_input_struct_i,
@@ -125,7 +125,7 @@ def test_get_unique_hop_host():
     assert results[0] == 2
 
 
-def test_get_unique_host_nonhost():
+def test_get_unique_host_nonhost() -> None:
     with pytest.raises(RuntimeError) as exc_info:
         get_unique_hop(
             mg_Mg,
diff --git a/pymatgen/analysis/diffusion/utils/tests/test_maggma.py b/pymatgen/analysis/diffusion/utils/tests/test_maggma.py
index fa47b8b..1a3aec6 100644
--- a/pymatgen/analysis/diffusion/utils/tests/test_maggma.py
+++ b/pymatgen/analysis/diffusion/utils/tests/test_maggma.py
@@ -17,14 +17,14 @@
 
 
 @pytest.fixture()
-def maggma_stores():
+def maggma_stores() -> dict[str, JSONStore]:
     return {
         "sgroups": JSONStore(f"{dir_path}/maggma_sgroup_store.json", key="group_id"),
         "materials": JSONStore(f"{dir_path}/maggma_materials_store.json", key="material_id"),
     }
 
 
-def test(maggma_stores):
+def test(maggma_stores: dict) -> None:
     base_ents, inserted_ents = get_entries_from_dbs(
         maggma_stores["sgroups"],
         maggma_stores["materials"],
@@ -33,7 +33,7 @@ def test(maggma_stores):
     )
 
     # check that the entries have been created
-    def has_mg(ent):
+    def has_mg(ent) -> bool:  # noqa: ANN001
         return "Mg" in ent.composition.as_dict()
 
     assert all(map(has_mg, inserted_ents))
diff --git a/pymatgen/analysis/diffusion/utils/tests/test_parse_entries.py b/pymatgen/analysis/diffusion/utils/tests/test_parse_entries.py
index 15ee5f1..9a11431 100644
--- a/pymatgen/analysis/diffusion/utils/tests/test_parse_entries.py
+++ b/pymatgen/analysis/diffusion/utils/tests/test_parse_entries.py
@@ -27,7 +27,7 @@
 
 
 class ParseEntriesTest(unittest.TestCase):
-    def setUp(self):
+    def setUp(self) -> None:
         d = loadfn(f"{dir_path}/parse_entry_test_vars.json")
         struct_uc = d["struct_uc"]
         self.li_ent = d["li_ent"]
@@ -53,7 +53,7 @@ def setUp(self):
         self.struct_inserted_1Li2 = get_inserted_on_base(self.base, self.inserted_1Li2, self.li_ent, self.sm)
         self.struct_inserted_2Li = get_inserted_on_base(self.base, self.inserted_2Li, self.li_ent, self.sm)
 
-    def _is_valid_inserted_ent(self, mapped_struct):
+    def _is_valid_inserted_ent(self, mapped_struct: Structure) -> bool:
         res = True
         for isite in mapped_struct.sites:
             if isite.species_string == "Li":
@@ -64,7 +64,7 @@ def _is_valid_inserted_ent(self, mapped_struct):
                     return False
         return res
 
-    def test_get_inserted_on_base(self):
+    def test_get_inserted_on_base(self) -> None:
         mapped_struct = get_inserted_on_base(self.base, self.inserted_1Li1, self.li_ent, self.sm)
         assert self._is_valid_inserted_ent(mapped_struct)
         assert mapped_struct[0].properties["insertion_energy"] == 5.0
@@ -75,7 +75,7 @@ def test_get_inserted_on_base(self):
         assert self._is_valid_inserted_ent(mapped_struct)
         assert mapped_struct[0].properties["insertion_energy"] == 4.0
 
-    def test_process_ents(self):
+    def test_process_ents(self) -> None:
         base_2_ent = ComputedStructureEntry(
             structure=self.base.structure * [[1, 1, 0], [1, -1, 0], [0, 0, 2]],
             energy=self.base.energy * 4,
@@ -90,16 +90,17 @@ def test_process_ents(self):
                 if i_insert_site.species_string == "Li":
                     assert i_insert_site.properties["insertion_energy"] == 4
 
-    def test_filter_and_merge(self):
+    def test_filter_and_merge(self) -> None:
         combined_struct = Structure.from_sites(
             self.struct_inserted_1Li1.sites + self.struct_inserted_1Li2.sites + self.struct_inserted_2Li.sites
         )
         filtered_struct = _filter_and_merge(combined_struct)
+        assert filtered_struct is not None
         for i_insert_site in filtered_struct:
             if i_insert_site.species_string == "Li":
                 assert i_insert_site.properties["insertion_energy"] in {4.5, 5.5}
 
-    def test_get_insertion_energy(self):
+    def test_get_insertion_energy(self) -> None:
         insert_energy = get_insertion_energy(self.base, self.inserted_1Li1, self.li_ent)
         basex2_ = ComputedStructureEntry(structure=self.base.structure * [1, 1, 2], energy=self.base.energy * 2)
         insert_energyx2 = get_insertion_energy(basex2_, self.inserted_1Li1, self.li_ent)
@@ -108,7 +109,7 @@ def test_get_insertion_energy(self):
         insert_energy = get_insertion_energy(self.base, self.inserted_2Li, self.li_ent)
         self.assertAlmostEqual(insert_energy, 4)
 
-    def test_get_all_sym_sites(self):
+    def test_get_all_sym_sites(self) -> None:
         struct11 = get_sym_migration_ion_sites(self.base.structure, self.inserted_1Li1.structure, migrating_ion="Li")
         assert struct11.composition["Li"] == 4
         struct12 = get_sym_migration_ion_sites(self.base.structure, self.inserted_1Li2.structure, migrating_ion="Li")
diff --git a/pyproject.toml b/pyproject.toml
index 60416b6..c91fb04 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -87,6 +87,7 @@ include = ["pymatgen.*"]
 target-version = "py39"
 line-length = 120
 lint.select = [
+    "ANN", # flake8-annotations
     "B", # flake8-bugbear
     "C4", # flake8-comprehensions
     "D", # pydocstyle
@@ -118,6 +119,9 @@ lint.select = [
     "YTT", # flake8-2020
 ]
 lint.ignore = [
+    "ANN003", # Missing type annotation for kwargs
+    "ANN101", # Missing type annotation for self in method
+    "ANN102", # Missing type annotation for cls in classmethod
     "B023", # Function definition does not bind loop variable
     "B028", # No explicit stacklevel keyword argument found
     "B904", # Within an except clause, raise exceptions with ...
diff --git a/tasks.py b/tasks.py
index 9667e7c..1b660d3 100644
--- a/tasks.py
+++ b/tasks.py
@@ -23,7 +23,7 @@
 
 
 @task
-def make_doc(ctx):
+def make_doc(ctx) -> None:
     with cd("docs_rst"):
         ctx.run("cp ../CHANGES.rst change_log.rst")
         ctx.run(
@@ -64,7 +64,7 @@ def make_doc(ctx):
 
 
 @task
-def set_ver(ctx):
+def set_ver(ctx) -> None:
 
     lines = []
     with open("pyproject.toml") as f:
@@ -80,7 +80,7 @@ def set_ver(ctx):
 
 
 @task
-def update_doc(ctx):
+def update_doc(ctx) -> None:
     make_doc(ctx)
     with cd("docs"):
         ctx.run("git add .")
@@ -89,7 +89,7 @@ def update_doc(ctx):
 
 
 @task
-def publish(ctx):
+def publish(ctx) -> None:
     ctx.run("rm dist/*.*", warn=True)
     ctx.run("python -m build", warn=True)
     ctx.run("twine upload --skip-existing dist/*.whl", warn=True)
@@ -97,7 +97,7 @@ def publish(ctx):
 
 
 @task
-def release_github(ctx):
+def release_github(ctx) -> None:
     payload = {
         "tag_name": "v" + NEW_VER,
         "target_commitish": "master",
@@ -115,12 +115,12 @@ def release_github(ctx):
 
 
 @task
-def test(ctx):
+def test(ctx) -> None:
     ctx.run("pytest pymatgen")
 
 
 @task
-def release(ctx):
+def release(ctx) -> None:
     set_ver(ctx)
     # test(ctx)
     update_doc(ctx)