diff --git a/.github/workflows/lint.yml b/.github/workflows/lint.yml
index cacd81dcf36..ec9da68e969 100644
--- a/.github/workflows/lint.yml
+++ b/.github/workflows/lint.yml
@@ -22,7 +22,7 @@ jobs:
 
       - name: Install dependencies
         run: |
-          pip install --upgrade ruff mypy
+          pip install --upgrade ruff mypy pyright
 
       - name: ruff
         run: |
@@ -31,7 +31,7 @@ jobs:
           ruff format --check .
 
       - name: mypy
-        run: |
-          mypy --version
-          rm -rf .mypy_cache
-          mypy ${{ github.event.repository.name }}
+        run: mypy ${{ github.event.repository.name }}
+
+      - name: pyright
+        run: pyright
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index a6745467041..f0b810db1d0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,13 +2,13 @@ exclude: ^(docs|tests/files|cmd_line|tasks.py)
 
 ci:
   autoupdate_schedule: monthly
-  skip: [mypy]
+  skip: [mypy, pyright]
   autofix_commit_msg: pre-commit auto-fixes
   autoupdate_commit_msg: pre-commit autoupdate
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.3.7
+    rev: v0.4.1
     hooks:
       - id: ruff
         args: [--fix, --unsafe-fixes]
@@ -35,7 +35,7 @@ repos:
         additional_dependencies: [tomli] # needed to read pyproject.toml below py3.11
 
   - repo: https://github.com/MarcoGorelli/cython-lint
-    rev: v0.16.0
+    rev: v0.16.2
     hooks:
       - id: cython-lint
         args: [--no-pycodestyle]
@@ -62,3 +62,8 @@ repos:
     hooks:
       - id: nbstripout
         args: [--drop-empty-cells, --keep-output]
+
+  - repo: https://github.com/RobertCraigie/pyright-python
+    rev: v1.1.359
+    hooks:
+      - id: pyright
diff --git a/dev_scripts/update_pt_data.py b/dev_scripts/update_pt_data.py
index b6a59fdd249..30f7f3c8b15 100644
--- a/dev_scripts/update_pt_data.py
+++ b/dev_scripts/update_pt_data.py
@@ -234,7 +234,7 @@ def gen_iupac_ordering():
 def add_electron_affinities():
     """Update the periodic table data file with electron affinities."""
 
-    req = requests.get("https://wikipedia.org/wiki/Electron_affinity_(data_page)", timeout=60)
+    req = requests.get("https://wikipedia.org/wiki/Electron_affinity_(data_page)", timeout=600)
     soup = BeautifulSoup(req.text, "html.parser")
     table = None
     for table in soup.find_all("table"):
diff --git a/pymatgen/analysis/adsorption.py b/pymatgen/analysis/adsorption.py
index 527b0304085..b3e34285316 100644
--- a/pymatgen/analysis/adsorption.py
+++ b/pymatgen/analysis/adsorption.py
@@ -680,19 +680,19 @@ def plot_slab(
         sites = list(reversed(sites))
         coords = np.array(reversed(coords))
     # Draw circles at sites and stack them accordingly
-    for n, coord in enumerate(coords):
-        radius = sites[n].species.elements[0].atomic_radius * scale
-        ax.add_patch(patches.Circle(coord[:2] - lattice_sum * (repeat // 2), radius, color="w", zorder=2 * n))
-        color = color_dict[sites[n].species.elements[0].symbol]
+    for idx, coord in enumerate(coords):
+        radius = sites[idx].species.elements[0].atomic_radius * scale
+        ax.add_patch(patches.Circle(coord[:2] - lattice_sum * (repeat // 2), radius, color="w", zorder=2 * idx))
+        color = color_dict[sites[idx].species.elements[0].symbol]
         ax.add_patch(
             patches.Circle(
                 coord[:2] - lattice_sum * (repeat // 2),
                 radius,
                 facecolor=color,
-                alpha=alphas[n],
+                alpha=alphas[idx],
                 edgecolor="k",
                 lw=0.3,
-                zorder=2 * n + 1,
+                zorder=2 * idx + 1,
             )
         )
     # Adsorption sites
@@ -714,7 +714,7 @@ def plot_slab(
         codes = [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY]
         vertices = [(np.array(vert) + corner).tolist() for vert in vertices]
         path = Path(vertices, codes)
-        patch = patches.PathPatch(path, facecolor="none", lw=2, alpha=0.5, zorder=2 * n + 2)
+        patch = patches.PathPatch(path, facecolor="none", lw=2, alpha=0.5, zorder=2 * len(coords) + 2)
         ax.add_patch(patch)
     ax.set_aspect("equal")
     center = corner + lattice_sum / 2.0
diff --git a/pymatgen/analysis/bond_dissociation.py b/pymatgen/analysis/bond_dissociation.py
index 7f6ddbbff44..b448ad06537 100644
--- a/pymatgen/analysis/bond_dissociation.py
+++ b/pymatgen/analysis/bond_dissociation.py
@@ -118,8 +118,9 @@ def fragment_and_process(self, bonds):
             bonds (list): bonds to process.
         """
         # Try to split the principle:
+        fragments: list[MoleculeGraph] = []
         try:
-            frags = self.mol_graph.split_molecule_subgraphs(bonds, allow_reverse=True)
+            fragments = self.mol_graph.split_molecule_subgraphs(bonds, allow_reverse=True)
             frag_success = True
 
         except MolGraphSplitError:
@@ -166,7 +167,7 @@ def fragment_and_process(self, bonds):
                         # We shouldn't ever encounter more than one good entry.
                         raise RuntimeError("There should only be one valid ring opening fragment! Exiting...")
             elif len(bonds) == 2:
-                raise RuntimeError("Should only be trying to break two bonds if multibreak is true! Exiting...")
+                raise RuntimeError("Should only be trying to break two bonds if multibreak=True! Exiting...")
             else:
                 raise ValueError("No reason to try and break more than two bonds at once! Exiting...")
             frag_success = False
@@ -176,19 +177,19 @@ def fragment_and_process(self, bonds):
             # As above, we begin by making sure we haven't already encountered an identical pair of fragments:
             frags_done = False
             for frag_pair in self.done_frag_pairs:
-                if frag_pair[0].isomorphic_to(frags[0]):
-                    if frag_pair[1].isomorphic_to(frags[1]):
+                if frag_pair[0].isomorphic_to(fragments[0]):
+                    if frag_pair[1].isomorphic_to(fragments[1]):
                         frags_done = True
                         break
-                elif frag_pair[1].isomorphic_to(frags[0]) and frag_pair[0].isomorphic_to(frags[1]):
+                elif frag_pair[1].isomorphic_to(fragments[0]) and frag_pair[0].isomorphic_to(fragments[1]):
                     frags_done = True
                     break
             if not frags_done:
                 # If we haven't, we save this pair and search for the relevant fragment entries:
-                self.done_frag_pairs += [frags]
+                self.done_frag_pairs += [fragments]
                 n_entries_for_this_frag_pair = 0
-                frag1_entries = self.search_fragment_entries(frags[0])
-                frag2_entries = self.search_fragment_entries(frags[1])
+                frag1_entries = self.search_fragment_entries(fragments[0])
+                frag2_entries = self.search_fragment_entries(fragments[1])
                 frag1_charges_found = []
                 frag2_charges_found = []
                 # We then check for our expected charges of each fragment:
@@ -200,14 +201,14 @@ def fragment_and_process(self, bonds):
                         frag2_charges_found += [frag2["initial_molecule"]["charge"]]
                 # If we're missing some of either, tell the user:
                 if len(frag1_charges_found) < len(self.expected_charges):
-                    bb = BabelMolAdaptor(frags[0].molecule)
+                    bb = BabelMolAdaptor(fragments[0].molecule)
                     pb_mol = bb.pybel_mol
                     smiles = pb_mol.write("smi").split()[0]
                     for charge in self.expected_charges:
                         if charge not in frag1_charges_found:
                             warnings.warn(f"Missing {charge=} for fragment {smiles}")
                 if len(frag2_charges_found) < len(self.expected_charges):
-                    bb = BabelMolAdaptor(frags[1].molecule)
+                    bb = BabelMolAdaptor(fragments[1].molecule)
                     pb_mol = bb.pybel_mol
                     smiles = pb_mol.write("smi").split()[0]
                     for charge in self.expected_charges:
diff --git a/pymatgen/analysis/bond_valence.py b/pymatgen/analysis/bond_valence.py
index 36099fb645e..b721595c77e 100644
--- a/pymatgen/analysis/bond_valence.py
+++ b/pymatgen/analysis/bond_valence.py
@@ -268,6 +268,7 @@ def get_valences(self, structure: Structure):
                 valences.append(vals)
 
         # make variables needed for recursion
+        attrib = []
         if structure.is_ordered:
             n_sites = np.array(list(map(len, equi_sites)))
             valence_min = np.array(list(map(min, valences)))
@@ -324,7 +325,6 @@ def _recurse(assigned=None):
         else:
             n_sites = np.array([len(sites) for sites in equi_sites])
             tmp = []
-            attrib = []
             for idx, n_site in enumerate(n_sites):
                 for _ in valences[idx]:
                     tmp.append(n_site)
diff --git a/pymatgen/analysis/chemenv/coordination_environments/chemenv_strategies.py b/pymatgen/analysis/chemenv/coordination_environments/chemenv_strategies.py
index 63a32c88460..b5c2c2f61db 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/chemenv_strategies.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/chemenv_strategies.py
@@ -9,7 +9,7 @@
 
 import abc
 import os
-from typing import TYPE_CHECKING, ClassVar
+from typing import TYPE_CHECKING
 
 import numpy as np
 from monty.json import MSONable
@@ -31,6 +31,8 @@
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 
 if TYPE_CHECKING:
+    from typing import ClassVar
+
     from typing_extensions import Self
 
 __author__ = "David Waroquiers"
@@ -256,6 +258,7 @@ def equivalent_site_index_and_transform(self, psite):
             Equivalent site in the unit cell, translations and symmetry transformation.
         """
         # Get the index of the site in the unit cell of which the PeriodicSite psite is a replica.
+        isite = 0
         try:
             isite = self.structure_environments.structure.index(psite)
         except ValueError:
@@ -283,6 +286,8 @@ def equivalent_site_index_and_transform(self, psite):
         # that gets back the site to the unit cell (Translation III)
         # TODO: check that these tolerances are needed, now that the structures are refined before analyzing envs
         tolerances = [1e-8, 1e-7, 1e-6, 1e-5, 1e-4]
+        d_this_site2 = (0, 0, 0)
+        sym_trafo = None
         for tolerance in tolerances:
             for sym_op in self.symops:
                 new_site = PeriodicSite(
diff --git a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometries.py b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometries.py
index de6d88630f1..2f0ecc3e23b 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometries.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometries.py
@@ -808,6 +808,9 @@ def edges(self, sites, permutation=None, input="sites"):
             coords = [site.coords for site in sites]
         elif input == "coords":
             coords = sites
+        else:
+            raise RuntimeError("Invalid input for edges.")
+
         # if permutation is None:
         #     coords = [site.coords for site in sites]
         # else:
diff --git a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
index 42bc3c2f305..34a529294ea 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/coordination_geometry_finder.py
@@ -775,6 +775,7 @@ def compute_structure_environments(
 
             to_add_from_hints = []
             nb_sets_info = {}
+            cn = 0
 
             for cn, nb_sets in struct_envs.neighbors_sets[isite].items():
                 if cn not in all_cns:
@@ -1534,8 +1535,7 @@ def coordination_geometry_symmetry_measures_separation_plane(
         algos = []
         perfect2local_maps = []
         local2perfect_maps = []
-        if testing:
-            separation_permutations = []
+        separation_permutations = []
         nplanes = 0
         for npoints in range(
             separation_plane_algo.minimum_number_of_points,
@@ -1793,9 +1793,10 @@ def _cg_csm_separation_plane(
         plane_found = False
         permutations = []
         permutations_symmetry_measures = []
-        if testing:
-            separation_permutations = []
+        separation_permutations = []
         dist_tolerances = dist_tolerances or DIST_TOLERANCES
+        algo = ""
+
         for dist_tolerance in dist_tolerances:
             algo = "NOT_FOUND"
             separation = local_plane.indices_separate(self.local_geometry._coords, dist_tolerance)
diff --git a/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py b/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
index 74afc331db8..9d2bd88e7d2 100644
--- a/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
+++ b/pymatgen/analysis/chemenv/coordination_environments/structure_environments.py
@@ -2230,6 +2230,7 @@ def __str__(self):
         if len(self.coord_geoms) == 0:
             out += " => No coordination in it <=\n"
             return out
+        mp_symbol = ""
         for key in self.coord_geoms:
             mp_symbol = key
             break
diff --git a/pymatgen/analysis/chemenv/utils/chemenv_config.py b/pymatgen/analysis/chemenv/utils/chemenv_config.py
index 99eeb8d1a62..667531261a3 100644
--- a/pymatgen/analysis/chemenv/utils/chemenv_config.py
+++ b/pymatgen/analysis/chemenv/utils/chemenv_config.py
@@ -78,12 +78,11 @@ def setup(self):
                 break
             elif test == "S":
                 config_file = self.save()
+                print(f"Configuration has been saved to file {config_file!r}")
                 break
             else:
                 print(" ... wrong key, try again ...")
             print()
-        if test == "S":
-            print(f"Configuration has been saved to file {config_file!r}")
 
     @property
     def has_materials_project_access(self):
diff --git a/pymatgen/analysis/chemenv/utils/scripts_utils.py b/pymatgen/analysis/chemenv/utils/scripts_utils.py
index 5b17eb0df15..31af2809ac3 100644
--- a/pymatgen/analysis/chemenv/utils/scripts_utils.py
+++ b/pymatgen/analysis/chemenv/utils/scripts_utils.py
@@ -78,6 +78,8 @@ def draw_cg(
         show_distorted:
         faces_color_override:
     """
+    csm_suffix = ""
+    perf_radius = 0
     if show_perfect:
         if csm_info is None:
             raise ValueError("Not possible to show perfect environment without csm_info")
@@ -219,6 +221,7 @@ def compute_environments(chemenv_configuration):
             for key_character, qq in questions.items():
                 print(f" - <{key_character}> for a structure from {string_sources[qq]['string']}")
             test = input(" ... ")
+            source_type = ""
             if test == "q":
                 break
             if test not in list(questions):
@@ -231,16 +234,22 @@ def compute_environments(chemenv_configuration):
                     continue
             else:
                 source_type = questions[test]
+
         else:
             found = False
             source_type = next(iter(questions.values()))
+
+        input_source = ""
         if found and len(questions) > 1:
-            input_source = test
+            input_source = test  # type: ignore[reportPossiblyUnboundVariable]
+
+        structure = None
         if source_type == "cif":
             if not found:
                 input_source = input("Enter path to CIF file : ")
             parser = CifParser(input_source)
             structure = parser.parse_structures(primitive=True)[0]
+
         elif source_type == "mp":
             if not found:
                 input_source = input('Enter materials project id (e.g. "mp-1902") : ')
@@ -248,6 +257,7 @@ def compute_environments(chemenv_configuration):
 
             with MPRester() as mpr:
                 structure = mpr.get_structure_by_material_id(input_source)
+
         lgf.setup_structure(structure)
         print(f"Computing environments for {structure.reduced_formula} ... ")
         se = lgf.compute_structure_environments(maximum_distance_factor=max_dist_factor)
@@ -258,7 +268,7 @@ def compute_environments(chemenv_configuration):
                 '("y" or "n", "d" with details, "g" to see the grid) : '
             )
             strategy = default_strategy
-            if test in ["y", "d", "g"]:
+            if test in {"y", "d", "g"}:
                 strategy.set_structure_environments(se)
                 for equiv_list in se.equivalent_sites:
                     site = equiv_list[0]
@@ -277,6 +287,7 @@ def compute_environments(chemenv_configuration):
                     comp = site.species
                     # ce = strategy.get_site_coordination_environment(site)
                     reduced_formula = comp.get_reduced_formula_and_factor()[0]
+                    the_cg = None
                     if strategy.uniquely_determines_coordination_environments:
                         ce = ces[0]
                         if ce is None:
@@ -349,6 +360,9 @@ def compute_environments(chemenv_configuration):
                     vis = StructureVis(show_polyhedron=False, show_unit_cell=True)
                     vis.show_help = False
                     first_time = False
+                else:
+                    vis = None  # TODO: following code logic seems buggy
+
                 vis.set_structure(se.structure)
                 strategy.set_structure_environments(se)
                 for site in se.structure:
diff --git a/pymatgen/analysis/chempot_diagram.py b/pymatgen/analysis/chempot_diagram.py
index 7a1469147aa..abfa6bf5e59 100644
--- a/pymatgen/analysis/chempot_diagram.py
+++ b/pymatgen/analysis/chempot_diagram.py
@@ -268,6 +268,8 @@ def _get_2d_plot(self, elements: list[Element], label_stable: bool | None, eleme
 
         if element_padding is not None and element_padding > 0:
             new_lims = self._get_new_limits_from_padding(domains, elem_indices, element_padding, self.default_min_limit)
+        else:
+            new_lims = []
 
         for formula, pts in domains.items():
             formula_elems = set(Composition(formula).elements)
@@ -326,6 +328,8 @@ def _get_3d_plot(
 
         if element_padding and element_padding > 0:
             new_lims = self._get_new_limits_from_padding(domains, elem_indices, element_padding, self.default_min_limit)
+        else:
+            new_lims = []
 
         for formula, pts in domains.items():
             entry = self.entry_dict[formula]
diff --git a/pymatgen/analysis/diffraction/tem.py b/pymatgen/analysis/diffraction/tem.py
index bab041edf66..22fc48292e8 100644
--- a/pymatgen/analysis/diffraction/tem.py
+++ b/pymatgen/analysis/diffraction/tem.py
@@ -474,12 +474,16 @@ def get_positions(self, structure: Structure, points: list) -> dict[tuple[int, i
         points = self.zone_axis_filter(points)
         # first is the max_d, min_r
         first_point_dict = self.get_first_point(structure, points)
+        first_point = (0, 0, 0)
+        first_d = 0.0
         for point, v in first_point_dict.items():
             first_point = point
             first_d = v
         spacings = self.get_interplanar_spacings(structure, points)
         # second is the first non-parallel-to-first-point vector when sorted.
         # note 000 is "parallel" to every plane vector.
+        second_point = (0, 0, 0)
+        second_d = 0.0
         for plane in sorted(spacings):
             second_point, second_d = plane, spacings[plane]
             if not self.is_parallel(structure, first_point, second_point):
diff --git a/pymatgen/analysis/elasticity/elastic.py b/pymatgen/analysis/elasticity/elastic.py
index 8db995325e7..9a6549d3010 100644
--- a/pymatgen/analysis/elasticity/elastic.py
+++ b/pymatgen/analysis/elasticity/elastic.py
@@ -979,7 +979,7 @@ def get_strain_state_dict(strains, stresses, eq_stress=None, tol: float = 1e-10,
         if add_eq:
             # add zero strain state
             mstrains = np.vstack([mstrains, np.zeros(6)])
-            mstresses = np.vstack([mstresses, veq_stress])
+            mstresses = np.vstack([mstresses, veq_stress])  # type: ignore[reportPossiblyUnboundVariable]
         # sort strains/stresses by strain values
         if sort:
             mstresses = mstresses[mstrains[:, ind[0]].argsort()]
diff --git a/pymatgen/analysis/functional_groups.py b/pymatgen/analysis/functional_groups.py
index b54a9c8fb06..f21f3ea894b 100644
--- a/pymatgen/analysis/functional_groups.py
+++ b/pymatgen/analysis/functional_groups.py
@@ -279,6 +279,7 @@ def get_basic_functional_groups(self, func_groups=None):
                         num_deviants += 1
 
                 if num_deviants <= 1:
+                    ring_group = []
                     for node in ring:
                         ring_group = copy.deepcopy(ring)
                         neighbors = self.molgraph.graph[node]
diff --git a/pymatgen/analysis/graphs.py b/pymatgen/analysis/graphs.py
index 974a02572b6..061b13b5567 100644
--- a/pymatgen/analysis/graphs.py
+++ b/pymatgen/analysis/graphs.py
@@ -11,7 +11,7 @@
 from itertools import combinations
 from operator import itemgetter
 from shutil import which
-from typing import TYPE_CHECKING, Any, Callable, cast
+from typing import TYPE_CHECKING, cast
 
 import networkx as nx
 import networkx.algorithms.isomorphism as iso
@@ -35,6 +35,7 @@
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
+    from typing import Any, Callable
 
     from igraph import Graph
     from numpy.typing import ArrayLike
@@ -571,9 +572,8 @@ def alter_edge(
                 f"Edge between {from_index} and {to_index} cannot be altered; no edge exists between those sites."
             )
 
-        if to_jimage is None:
-            edge_index = 0
-        else:
+        edge_index = 0
+        if to_jimage is not None:
             for idx, properties in existing_edges.items():
                 if properties["to_jimage"] == to_jimage:
                     edge_index = idx
@@ -606,6 +606,7 @@ def break_edge(
         if to_jimage is None:
             raise ValueError("Image must be supplied, to avoid ambiguity.")
 
+        edge_index = 0
         if existing_edges:
             for idx, props in existing_edges.items():
                 if props["to_jimage"] == to_jimage:
@@ -1489,6 +1490,8 @@ def get_subgraphs_as_molecules(self, use_weights: bool = False) -> list[Molecule
         # without adding extra logic
         if getattr(self, "_supercell_sg", None) is None:
             self._supercell_sg = supercell_sg = self * (3, 3, 3)
+        else:
+            raise RuntimeError("Supercell spacegroup is not None.")
 
         # make undirected to find connected subgraphs
         supercell_sg.graph = nx.Graph(supercell_sg.graph)
diff --git a/pymatgen/analysis/local_env.py b/pymatgen/analysis/local_env.py
index aaa2b8336fc..643721995ad 100644
--- a/pymatgen/analysis/local_env.py
+++ b/pymatgen/analysis/local_env.py
@@ -1631,13 +1631,15 @@ def get_nn_info(self, structure: Structure, n: int):
         siw = []
 
         for bond in bonds:
-            capture_bond = False
             if bond.site1 == structure[n]:
                 site = bond.site2
                 capture_bond = True
             elif bond.site2 == structure[n]:
                 site = bond.site1
                 capture_bond = True
+            else:
+                site = None
+                capture_bond = False
 
             if capture_bond:
                 index = structure.index(site)
@@ -2927,6 +2929,7 @@ def get_order_parameters(
         # norms = [[[] for j in range(nneigh)] for t in self._types]
 
         # First, coordination number and distance-based OPs.
+        typ = ""
         for idx, typ in enumerate(self._types):
             if typ == "cn":
                 ops[idx] = n_neighbors / self._params[idx]["norm"]
@@ -2978,9 +2981,9 @@ def get_order_parameters(
         # (Peters, J. Chem. Phys., 131, 244103, 2009;
         #  Zimmermann et al., J. Am. Chem. Soc., under revision, 2015).
         if self._geomops:
-            gaussthetak: list[float] = [0 for t in self._types]  # not used by all OPs
-            qsp_theta = [[[] for j in range(n_neighbors)] for t in self._types]  # type: ignore
-            norms = [[[] for j in range(n_neighbors)] for t in self._types]  # type: ignore
+            gaussthetak: list[float] = [0 for _t in self._types]  # not used by all OPs
+            qsp_theta = [[[] for _j in range(n_neighbors)] for _t in self._types]  # type: ignore
+            norms = [[[] for _j in range(n_neighbors)] for _t in self._types]  # type: ignore
             ipi = 1 / pi
             piover2 = pi / 2.0
             onethird = 1 / 3
@@ -2988,6 +2991,7 @@ def get_order_parameters(
             for j in range(n_neighbors):  # Neighbor j is put to the North pole.
                 zaxis = rij_norm[j]
                 kc = 0
+                idx = 0
                 for k in range(n_neighbors):  # From neighbor k, we construct
                     if j != k:  # the prime meridian.
                         for idx in range(len(self._types)):
@@ -3001,42 +3005,46 @@ def get_order_parameters(
                         else:
                             xaxis = xaxis / np.linalg.norm(xaxis)
                             flag_xaxis = False
+
                         if self._comp_azi:
                             flag_yaxis = True
                             yaxis = np.cross(zaxis, xaxis)
                             if np.linalg.norm(yaxis) > very_small:
                                 yaxis = yaxis / np.linalg.norm(yaxis)
                                 flag_yaxis = False
+                        else:
+                            yaxis = None
+                            flag_yaxis = False
 
                         # Contributions of j-i-k angles, where i represents the
                         # central atom and j and k two of the neighbors.
                         for idx, typ in enumerate(self._types):
-                            if typ in ["bent", "sq_pyr_legacy"]:
+                            if typ in {"bent", "sq_pyr_legacy"}:
                                 tmp = self._params[idx]["IGW_TA"] * (thetak * ipi - self._params[idx]["TA"])
                                 qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp)
                                 norms[idx][j][kc] += 1
-                            elif typ in ["tri_plan", "tri_plan_max", "tet", "tet_max"]:
+                            elif typ in {"tri_plan", "tri_plan_max", "tet", "tet_max"}:
                                 tmp = self._params[idx]["IGW_TA"] * (thetak * ipi - self._params[idx]["TA"])
                                 gaussthetak[idx] = exp(-0.5 * tmp * tmp)
                                 if typ in ["tri_plan_max", "tet_max"]:
                                     qsp_theta[idx][j][kc] += gaussthetak[idx]
                                     norms[idx][j][kc] += 1
-                            elif typ in ["T", "tri_pyr", "sq_pyr", "pent_pyr", "hex_pyr"]:
+                            elif typ in {"T", "tri_pyr", "sq_pyr", "pent_pyr", "hex_pyr"}:
                                 tmp = self._params[idx]["IGW_EP"] * (thetak * ipi - 0.5)
                                 qsp_theta[idx][j][kc] += exp(-0.5 * tmp * tmp)
                                 norms[idx][j][kc] += 1
-                            elif typ in [
+                            elif typ in {
                                 "sq_plan",
                                 "oct",
                                 "oct_legacy",
                                 "cuboct",
                                 "cuboct_max",
-                            ]:
+                            }:
                                 if thetak >= self._params[idx]["min_SPP"]:
                                     tmp = self._params[idx]["IGW_SPP"] * (thetak * ipi - 1.0)
                                     qsp_theta[idx][j][kc] += self._params[idx]["w_SPP"] * exp(-0.5 * tmp * tmp)
                                     norms[idx][j][kc] += self._params[idx]["w_SPP"]
-                            elif typ in [
+                            elif typ in {
                                 "see_saw_rect",
                                 "tri_bipyr",
                                 "sq_bipyr",
@@ -3045,7 +3053,7 @@ def get_order_parameters(
                                 "oct_max",
                                 "sq_plan_max",
                                 "hex_plan_max",
-                            ]:
+                            }:
                                 if thetak < self._params[idx]["min_SPP"]:
                                     tmp = (
                                         self._params[idx]["IGW_EP"] * (thetak * ipi - 0.5)
@@ -3078,8 +3086,10 @@ def get_order_parameters(
                                 thetam = acos(tmp)
                                 x_two_axis_tmp = gramschmidt(rij_norm[m], zaxis)
                                 norm = np.linalg.norm(x_two_axis_tmp)
+                                phi2 = 0.0
                                 if norm < very_small:
                                     flag_xtwoaxis = True
+                                    phi = 0.0
                                 else:
                                     xtwoaxis = x_two_axis_tmp / norm
                                     phi = acos(max(-1.0, min(np.inner(xtwoaxis, xaxis), 1.0)))
@@ -3092,7 +3102,7 @@ def get_order_parameters(
                                 # South pole contributions of m.
                                 if (
                                     typ
-                                    in [
+                                    in {
                                         "tri_bipyr",
                                         "sq_bipyr",
                                         "pent_bipyr",
@@ -3101,7 +3111,7 @@ def get_order_parameters(
                                         "sq_plan_max",
                                         "hex_plan_max",
                                         "see_saw_rect",
-                                    ]
+                                    }
                                     and thetam >= self._params[idx]["min_SPP"]
                                 ):
                                     tmp = self._params[idx]["IGW_SPP"] * (thetam * ipi - 1.0)
@@ -3112,12 +3122,12 @@ def get_order_parameters(
                                 # angles between plane j-i-k and i-m vector.
                                 if not flag_xaxis and not flag_xtwoaxis:
                                     for idx, typ in enumerate(self._types):
-                                        if typ in [
+                                        if typ in {
                                             "tri_plan",
                                             "tri_plan_max",
                                             "tet",
                                             "tet_max",
-                                        ]:
+                                        }:
                                             tmp = self._params[idx]["IGW_TA"] * (thetam * ipi - self._params[idx]["TA"])
                                             tmp2 = (
                                                 cos(self._params[idx]["fac_AA"] * phi)
@@ -3133,13 +3143,13 @@ def get_order_parameters(
                                             tmp4 = 1 if typ == "pent_plan_max" else gaussthetak[idx]
                                             qsp_theta[idx][j][kc] += tmp4 * exp(-0.5 * tmp2 * tmp2) * tmp3 * tmp3
                                             norms[idx][j][kc] += 1
-                                        elif typ in [
+                                        elif typ in {
                                             "T",
                                             "tri_pyr",
                                             "sq_pyr",
                                             "pent_pyr",
                                             "hex_pyr",
-                                        ]:
+                                        }:
                                             tmp = (
                                                 cos(self._params[idx]["fac_AA"] * phi)
                                                 ** self._params[idx]["exp_cos_AA"]
@@ -3163,7 +3173,7 @@ def get_order_parameters(
                                                         tmp * self._params[idx][6] * self._params[idx][7]
                                                     )
                                                 norms[idx][j][kc] += 1
-                                        elif typ in [
+                                        elif typ in {
                                             "tri_bipyr",
                                             "sq_bipyr",
                                             "pent_bipyr",
@@ -3171,7 +3181,7 @@ def get_order_parameters(
                                             "oct_max",
                                             "sq_plan_max",
                                             "hex_plan_max",
-                                        ]:
+                                        }:
                                             if (
                                                 thetam < self._params[idx]["min_SPP"]
                                                 and thetak < self._params[idx]["min_SPP"]
@@ -3265,7 +3275,7 @@ def get_order_parameters(
 
             # Normalize Peters-style OPs.
             for idx, typ in enumerate(self._types):
-                if typ in [
+                if typ in {
                     "tri_plan",
                     "tet",
                     "bent",
@@ -3274,13 +3284,13 @@ def get_order_parameters(
                     "oct_legacy",
                     "cuboct",
                     "pent_plan",
-                ]:
+                }:
                     ops[idx] = tmp_norm = 0.0
                     for j in range(n_neighbors):
                         ops[idx] += sum(qsp_theta[idx][j])
                         tmp_norm += float(sum(norms[idx][j]))
                     ops[idx] = ops[idx] / tmp_norm if tmp_norm > 1.0e-12 else None  # type: ignore
-                elif typ in [
+                elif typ in {
                     "T",
                     "tri_pyr",
                     "see_saw_rect",
@@ -3299,7 +3309,7 @@ def get_order_parameters(
                     "cuboct_max",
                     "hex_plan_max",
                     "sq_face_cap_trig_pris",
-                ]:
+                }:
                     ops[idx] = None  # type: ignore[call-overload]
                     if n_neighbors > 1:
                         for j in range(n_neighbors):
@@ -3353,7 +3363,7 @@ def get_order_parameters(
             dhalf = max(distjk_unique) / 2 if len(distjk_unique) > 0 else 0
 
             for idx, typ in enumerate(self._types):
-                if typ in ("reg_tri", "sq"):
+                if typ in {"reg_tri", "sq"}:
                     if n_neighbors < 3:
                         ops[idx] = None  # type: ignore[call-overload]
                     else:
@@ -3361,9 +3371,11 @@ def get_order_parameters(
                         if typ == "reg_tri":
                             a = 2 * asin(b / (2 * sqrt(h * h + (b / (2 * cos(3 * pi / 18))) ** 2)))  # type: ignore
                             nmax = 3
-                        elif typ == "sq":
+
+                        else:
                             a = 2 * asin(b / (2 * sqrt(h * h + dhalf * dhalf)))  # type: ignore
                             nmax = 4
+
                         for j in range(min([n_neighbors, nmax])):
                             ops[idx] = ops[idx] * exp(-0.5 * ((aijs[j] - a) * self._params[idx][0]) ** 2)
 
diff --git a/pymatgen/analysis/magnetism/heisenberg.py b/pymatgen/analysis/magnetism/heisenberg.py
index 8d5922aa4db..868c00f4fa1 100644
--- a/pymatgen/analysis/magnetism/heisenberg.py
+++ b/pymatgen/analysis/magnetism/heisenberg.py
@@ -95,6 +95,7 @@ def __init__(self, ordered_structures, energies, cutoff=0, tol: float = 0.02):
         # Check how many commensurate graphs we found
         if len(self.sgraphs) < 2:
             raise SystemExit("We need at least 2 unique orderings.")
+
         # Set attributes
         self._get_nn_dict()
         self._get_exchange_df()
@@ -211,6 +212,7 @@ def _get_nn_dict(self):
                 dist = round(cs[-1], 2)  # i_j distance
 
                 j = cs[2]  # j index
+                j_key = None
                 for key, value in unique_site_ids.items():
                     if j in key:
                         j_key = value
@@ -273,6 +275,7 @@ def _get_exchange_df(self):
 
             # Loop over all sites in each graph and compute |S_i . S_j|
             # for n+1 unique graphs to compute n exchange params
+            order = ""
             for _graph in sgraphs:
                 sgraph = sgraphs_copy.pop(0)
                 ex_row = pd.DataFrame(np.zeros((1, n_nn_j + 1)), index=[sgraph_index], columns=columns)
@@ -281,6 +284,7 @@ def _get_exchange_df(self):
                     # s_i_sign = np.sign(sgraph.structure.site_properties['magmom'][i])
                     s_i = sgraph.structure.site_properties["magmom"][idx]
 
+                    i_index = None
                     for k, v in unique_site_ids.items():
                         if idx in k:
                             i_index = v
@@ -297,6 +301,7 @@ def _get_exchange_df(self):
                         # s_j_sign = np.sign(sgraph.structure.site_properties['magmom'][j_site])
                         s_j = sgraph.structure.site_properties["magmom"][j_site]
 
+                        j_index = None
                         for k, v in unique_site_ids.items():
                             if j_site in k:
                                 j_index = v
@@ -308,6 +313,7 @@ def _get_exchange_df(self):
                             order = "-nnn"
                         elif abs(dist - dists["nnnn"]) <= tol:
                             order = "-nnnn"
+
                         j_ij = f"{i_index}-{j_index}{order}"
                         j_ji = f"{j_index}-{i_index}{order}"
 
@@ -393,6 +399,8 @@ def get_low_energy_orderings(self):
         fm_struct, afm_struct = None, None
         mag_min = np.inf
         mag_max = 0.001
+        fm_e = 0
+        afm_e = 0
         fm_e_min = 0
         afm_e_min = 0
 
@@ -588,15 +596,16 @@ def _get_j_exc(self, i, j, dist):
             j_exc (float): Exchange parameter in meV
         """
         # Get unique site identifiers
+        i_index = 0
+        j_index = 0
         for k, v in self.unique_site_ids.items():
             if i in k:
                 i_index = v
             if j in k:
                 j_index = v
 
-        order = ""
-
         # Determine order of interaction
+        order = ""
         if abs(dist - self.dists["nn"]) <= self.tol:
             order = "-nn"
         elif abs(dist - self.dists["nnn"]) <= self.tol:
@@ -958,15 +967,16 @@ def _get_j_exc(self, i, j, dist):
             j_exc (float): Exchange parameter in meV
         """
         # Get unique site identifiers
+        i_index = 0
+        j_index = 0
         for k in self.unique_site_ids:
             if i in k:
                 i_index = self.unique_site_ids[k]
             if j in k:
                 j_index = self.unique_site_ids[k]
 
-        order = ""
-
         # Determine order of interaction
+        order = ""
         if abs(dist - self.dists["nn"]) <= self.tol:
             order = "-nn"
         elif abs(dist - self.dists["nnn"]) <= self.tol:
diff --git a/pymatgen/analysis/molecule_matcher.py b/pymatgen/analysis/molecule_matcher.py
index 88b3115b8c2..ab99882b565 100644
--- a/pymatgen/analysis/molecule_matcher.py
+++ b/pymatgen/analysis/molecule_matcher.py
@@ -32,7 +32,7 @@
 
     from pymatgen.io.babel import BabelMolAdaptor
 except ImportError:
-    openbabel = None
+    openbabel = BabelMolAdaptor = None  # type: ignore[misc]
 
 if TYPE_CHECKING:
     from typing_extensions import Self
@@ -894,6 +894,8 @@ def match(self, mol: Molecule, ignore_warning: bool = False) -> tuple[np.ndarray
         rmsd = np.inf
 
         # Generate all permutation grouped/sorted by the elements
+        p_inds = []
+        U = np.empty(0)
         for p_inds_test in self.permutations(p_atoms):
             p_centroid_test = p_centroid[p_inds_test]
             U_test = self.kabsch(p_centroid_test, q_centroid)
@@ -991,6 +993,8 @@ def match(self, p: Molecule):
         rmsd = np.inf
 
         # Generate all permutation grouped/sorted by the elements
+        inds = []
+        U = np.empty(0)
         for p_inds_test in self.permutations(p_atoms, p_centroid, p_weights, q_atoms, q_centroid, q_weights):
             p_centroid_test = p_centroid[p_inds_test]
             U_test = self.kabsch(p_centroid_test, q_centroid)
@@ -1262,6 +1266,7 @@ def permutations(self, p: Molecule):
 
         # starting matches (only based on element)
         partial_matches = [[j] for j in range(self.N) if p_atoms[j] == q_atoms[0]]
+        matches: list = []
 
         for idx in range(1, self.N):
             # extending the target fragment with then next atom
diff --git a/pymatgen/analysis/phase_diagram.py b/pymatgen/analysis/phase_diagram.py
index 397cd6cc90d..20623224b9b 100644
--- a/pymatgen/analysis/phase_diagram.py
+++ b/pymatgen/analysis/phase_diagram.py
@@ -1258,7 +1258,9 @@ def getmu_vertices_stability_phase(self, target_comp, dep_elt, tol_en=1e-2):
                                 break
                         if not already_in:
                             all_coords.append(res)
-        return all_coords
+
+                return all_coords
+        return None
 
     def get_chempot_range_stability_phase(self, target_comp, open_elt):
         """
@@ -2939,6 +2941,11 @@ def _create_plotly_element_annotations(self):
                 font_dict = {"color": "#000000", "size": 24.0}
                 opacity = 1.0
 
+            else:
+                clean_formula = ""
+                font_dict = {}
+                opacity = 0
+
             offset = 0.03 if self._dim == 2 else 0.06
 
             if x < 0.4:
@@ -3835,6 +3842,9 @@ def order_phase_diagram(lines, stable_entries, unstable_entries, ordering):
     xleft = 1000.0
     xright = -1000.0
 
+    nameup = ""
+    nameleft = ""
+    nameright = ""
     for coord in stable_entries:
         if coord[0] > xright:
             xright = coord[0]
diff --git a/pymatgen/analysis/piezo_sensitivity.py b/pymatgen/analysis/piezo_sensitivity.py
index 5cb94c29863..0ba017a1240 100644
--- a/pymatgen/analysis/piezo_sensitivity.py
+++ b/pymatgen/analysis/piezo_sensitivity.py
@@ -15,7 +15,7 @@
     from phonopy import Phonopy
     from phonopy.harmonic import dynmat_to_fc as dyntofc
 except ImportError:
-    Phonopy = None
+    Phonopy = dyntofc = None
 
 
 if TYPE_CHECKING:
diff --git a/pymatgen/analysis/quasiharmonic.py b/pymatgen/analysis/quasiharmonic.py
index 63fea3bcc04..e743b7e7f63 100644
--- a/pymatgen/analysis/quasiharmonic.py
+++ b/pymatgen/analysis/quasiharmonic.py
@@ -133,6 +133,7 @@ def optimize_gibbs_free_energy(self):
         )
 
         for temp in temperatures:
+            G_opt = V_opt = None
             try:
                 G_opt, V_opt = self.optimizer(temp)
             except Exception:
diff --git a/pymatgen/analysis/structure_analyzer.py b/pymatgen/analysis/structure_analyzer.py
index 7984561366b..2557c25e997 100644
--- a/pymatgen/analysis/structure_analyzer.py
+++ b/pymatgen/analysis/structure_analyzer.py
@@ -447,6 +447,9 @@ def parse_oxide(self) -> tuple[str, int]:
                 for species, occu in site.species.items():
                     elem_map[species.element] += occu
             comp = Composition(elem_map)
+        else:
+            raise TypeError("Invalid type for element.")
+
         if Element("O") not in comp or comp.is_element:
             return "None", 0
 
@@ -466,6 +469,7 @@ def parse_oxide(self) -> tuple[str, int]:
         is_superoxide = False
         is_peroxide = False
         is_ozonide = False
+        bond_atoms = []
         if np.any(dist_matrix < relative_cutoff * 1.35):
             bond_atoms = np.where(dist_matrix < relative_cutoff * 1.35)[0]
             is_superoxide = True
@@ -475,10 +479,9 @@ def parse_oxide(self) -> tuple[str, int]:
         if is_superoxide and len(bond_atoms) > len(set(bond_atoms)):
             is_superoxide = False
             is_ozonide = True
-        try:
-            n_bonds = len(set(bond_atoms))
-        except UnboundLocalError:
-            n_bonds = 0
+
+        n_bonds = len(set(bond_atoms))
+
         if is_ozonide:
             str_oxide = "ozonide"
         elif is_superoxide:
diff --git a/pymatgen/analysis/structure_matcher.py b/pymatgen/analysis/structure_matcher.py
index 59b3a5a206b..8cd2364d6f6 100644
--- a/pymatgen/analysis/structure_matcher.py
+++ b/pymatgen/analysis/structure_matcher.py
@@ -1006,6 +1006,7 @@ def get_best_electronegativity_anonymous_mapping(self, struct1: Structure, struc
 
         if matches:
             min_X_diff = np.inf
+            best = None
             for match in matches:
                 X_diff = 0
                 for key, val in match[0].items():
diff --git a/pymatgen/analysis/surface_analysis.py b/pymatgen/analysis/surface_analysis.py
index 9535116bbc0..51d526c2a5f 100644
--- a/pymatgen/analysis/surface_analysis.py
+++ b/pymatgen/analysis/surface_analysis.py
@@ -214,6 +214,7 @@ def surface_energy(self, ucell_entry, ref_entries=None):
         # Next, we add the contribution to the bulk energy from
         # the variable element (the element without a ref_entry),
         # as a function of the other elements
+        ref_el = None
         for ref_el in ucell_comp.as_dict():
             if str(ref_el) not in ref_entries_dict:
                 break
@@ -955,6 +956,8 @@ def chempot_vs_gamma(
                 stable_u_range_dict = self.stable_u_range_dict(
                     chempot_range, ref_delu, no_doped=no_doped, delu_dict=delu_dict, miller_index=hkl
                 )
+            else:
+                stable_u_range_dict = {}
 
             already_labelled = []
             label = ""
@@ -1021,6 +1024,7 @@ def monolayer_vs_BE(self, plot_eads=False):
             Plot: Plot of binding energy vs monolayer for all facets.
         """
         ax = pretty_plot(width=8, height=7)
+        ads_entry = clean_entry = None
         for hkl in self.all_slab_entries:
             ml_be_dict = {}
             for clean_entry in self.all_slab_entries[hkl]:
diff --git a/pymatgen/analysis/topological/spillage.py b/pymatgen/analysis/topological/spillage.py
index ed1b9300bb3..1481e4ada79 100644
--- a/pymatgen/analysis/topological/spillage.py
+++ b/pymatgen/analysis/topological/spillage.py
@@ -67,13 +67,16 @@ def overlap_so_spinpol(self):
         # so_occs = np.array([np.array(so.band_energy)[:, :, 2]])
         so_nkpts = len(so_k)
 
-        nelec_list = []
+        n_elec_list: list[list[int]] = []
+        cup = cdn = 0
         for nk1 in range(1, n_kpts_noso + 1):  # no spin orbit kpoints loop
-            knoso = noso_kvecs[nk1 - 1, :]
+            k_no_spin_orbit = noso_kvecs[nk1 - 1, :]
             for nk2 in range(1, so_nkpts + 1):  # spin orbit
                 kso = so_kvecs[nk2 - 1, :]
                 if (
-                    self.isclose(kso[0], knoso[0]) and self.isclose(kso[1], knoso[1]) and self.isclose(kso[2], knoso[2])
+                    self.isclose(kso[0], k_no_spin_orbit[0])
+                    and self.isclose(kso[1], k_no_spin_orbit[1])
+                    and self.isclose(kso[2], k_no_spin_orbit[2])
                 ):  # do kpoints match?
                     for c, e in enumerate(noso_occs[0, nk1 - 1, :]):
                         if e < 0.5:
@@ -84,14 +87,14 @@ def overlap_so_spinpol(self):
                             cdn = c
                             break
 
-                    nelec_list.append([cup, cdn, cup + cdn])
-        n_arr = np.array(nelec_list)
+                    n_elec_list.append([cup, cdn, cup + cdn])
+        n_arr = np.array(n_elec_list)
 
         n_up = int(round(np.mean(n_arr[:, 0])))
         n_dn = int(round(np.mean(n_arr[:, 1])))
         n_tot = int(round(np.mean(n_arr[:, 2])))
 
-        nelec = int(n_tot)
+        n_elec = int(n_tot)
 
         # noso_homo_up = np.max(noso_bands[0, :, n_up - 1])
         # noso_lumo_up = np.min(noso_bands[0, :, n_up])
@@ -99,13 +102,13 @@ def overlap_so_spinpol(self):
         # noso_homo_dn = np.max(noso_bands[1, :, n_dn - 1])
         # noso_lumo_dn = np.min(noso_bands[1, :, n_dn])
 
-        so_homo = np.max(so_bands[0, :, nelec - 1])
-        so_lumo = np.min(so_bands[0, :, nelec])
+        so_homo = np.max(so_bands[0, :, n_elec - 1])
+        so_lumo = np.min(so_bands[0, :, n_elec])
 
         # noso_direct_up = np.min(noso_bands[0, :, n_up] - noso_bands[0, :, n_up - 1])
         # noso_direct_dn = np.min(noso_bands[1, :, n_dn] - noso_bands[1, :, n_dn - 1])
 
-        so_direct = np.min(so_bands[0, :, nelec] - so_bands[0, :, nelec - 1])
+        so_direct = np.min(so_bands[0, :, n_elec] - so_bands[0, :, n_elec - 1])
 
         noso_direct = 1000000.0
         noso_homo = -10000000.0
@@ -127,11 +130,13 @@ def overlap_so_spinpol(self):
 
         nelec_tot = 0.0
         for nk1 in range(1, n_kpts_noso + 1):  # no spin orbit kpoints loop
-            knoso = noso_kvecs[nk1 - 1, :]
+            k_no_spin_orbit = noso_kvecs[nk1 - 1, :]
             for nk2 in range(1, so_nkpts + 1):  # spin orbit
                 kso = so_kvecs[nk2 - 1, :]
                 if (
-                    self.isclose(kso[0], knoso[0]) and self.isclose(kso[1], knoso[1]) and self.isclose(kso[2], knoso[2])
+                    self.isclose(kso[0], k_no_spin_orbit[0])
+                    and self.isclose(kso[1], k_no_spin_orbit[1])
+                    and self.isclose(kso[2], k_no_spin_orbit[2])
                 ):  # do kpoints match?
                     # changes section 2
                     nelec_up = n_arr[nk1 - 1, 0]
@@ -193,7 +198,7 @@ def overlap_so_spinpol(self):
                                 nk1,
                                 nk2,
                                 kso,
-                                knoso,
+                                k_no_spin_orbit,
                                 np.real(gamma_k[-1]),
                                 "!!!!!!!!!!",
                             )
diff --git a/pymatgen/analysis/transition_state.py b/pymatgen/analysis/transition_state.py
index 56be52650a9..dab2503fd24 100644
--- a/pymatgen/analysis/transition_state.py
+++ b/pymatgen/analysis/transition_state.py
@@ -305,6 +305,10 @@ def combine_neb_plots(neb_analyses, arranged_neb_analyses=False, reverse_plot=Fa
         a NEBAnalysis object
     """
     x = StructureMatcher()
+    neb1_structures = []
+    neb1_energies = []
+    neb1_forces = []
+    neb1_r = []
     for neb_index, neb in enumerate(neb_analyses):
         if neb_index == 0:
             neb1 = neb
diff --git a/pymatgen/analysis/wulff.py b/pymatgen/analysis/wulff.py
index 10c0c13362b..2f558f26ec7 100644
--- a/pymatgen/analysis/wulff.py
+++ b/pymatgen/analysis/wulff.py
@@ -320,7 +320,7 @@ def get_line_in_facet(self, facet):
         """Returns the sorted pts in a facet used to draw a line."""
         lines = list(facet.outer_lines)
         pt = []
-        prev = None
+        prev = line = None
         while len(lines) > 0:
             if prev is None:
                 line = lines.pop(0)
diff --git a/pymatgen/analysis/xas/spectrum.py b/pymatgen/analysis/xas/spectrum.py
index cb9c68d77dc..5af607f3160 100644
--- a/pymatgen/analysis/xas/spectrum.py
+++ b/pymatgen/analysis/xas/spectrum.py
@@ -4,7 +4,7 @@
 
 import math
 import warnings
-from typing import Literal
+from typing import TYPE_CHECKING
 
 import numpy as np
 from scipy.interpolate import interp1d
@@ -13,6 +13,9 @@
 from pymatgen.core.spectrum import Spectrum
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 
+if TYPE_CHECKING:
+    from typing import Literal
+
 __author__ = "Chen Zheng, Yiming Chen"
 __copyright__ = "Copyright 2012, The Materials Project"
 __version__ = "3.0"
@@ -234,6 +237,7 @@ def site_weighted_spectrum(xas_list: list[XAS], num_samples: int = 500) -> XAS:
     maxes, mines = [], []
     fs = []
     multiplicities = []
+    xas = None
 
     for xas in xas_list:
         multiplicity = len(ss.find_equivalent_sites(ss[xas.absorbing_index]))
diff --git a/pymatgen/command_line/bader_caller.py b/pymatgen/command_line/bader_caller.py
index 18c5413a437..c2cd5d49840 100644
--- a/pymatgen/command_line/bader_caller.py
+++ b/pymatgen/command_line/bader_caller.py
@@ -533,8 +533,9 @@ def _get_filepath(filename: str, msg: str = "") -> str | None:
         return paths[0]
 
     chgcar_path = _get_filepath("CHGCAR", "Could not find CHGCAR!")
-    if chgcar_path is not None:
-        chgcar = Chgcar.from_file(chgcar_path)
+    if chgcar_path is None:
+        raise FileNotFoundError("Could not find CHGCAR!")
+    chgcar = Chgcar.from_file(chgcar_path)
 
     aeccar0_path = _get_filepath("AECCAR0")
     if not aeccar0_path:
diff --git a/pymatgen/command_line/chargemol_caller.py b/pymatgen/command_line/chargemol_caller.py
index 7f2fc1bcd91..7fafc7e00da 100644
--- a/pymatgen/command_line/chargemol_caller.py
+++ b/pymatgen/command_line/chargemol_caller.py
@@ -58,6 +58,7 @@
 
 if TYPE_CHECKING:
     from pathlib import Path
+    from typing import Literal
 
     from pymatgen.core import Structure
 
@@ -265,15 +266,15 @@ def get_charge_transfer(self, atom_index, charge_type="ddec"):
         Returns:
             float: charge transferred at atom_index
         """
-        if charge_type.lower() not in ["ddec", "cm5"]:
-            raise ValueError(f"Invalid {charge_type=}")
         if charge_type.lower() == "ddec":
-            charge_transfer = -self.ddec_charges[atom_index]
-        elif charge_type.lower() == "cm5":
-            charge_transfer = -self.cm5_charges[atom_index]
-        return charge_transfer
+            return -self.ddec_charges[atom_index]
 
-    def get_charge(self, atom_index, nelect=None, charge_type="ddec"):
+        if charge_type.lower() == "cm5":
+            return -self.cm5_charges[atom_index]
+
+        raise ValueError(f"Invalid {charge_type=}")
+
+    def get_charge(self, atom_index, nelect=None, charge_type: Literal["ddec", "cm5"] = "ddec"):
         """Convenience method to get the charge on a particular atom using the same
         sign convention as the BaderAnalysis. Note that this is *not* the partial
         atomic charge. This value is nelect (e.g. ZVAL from the POTCAR) + the
@@ -305,7 +306,7 @@ def get_charge(self, atom_index, nelect=None, charge_type="ddec"):
             charge = None
         return charge
 
-    def get_partial_charge(self, atom_index, charge_type="ddec"):
+    def get_partial_charge(self, atom_index, charge_type: Literal["ddec", "cm5"] = "ddec"):
         """Convenience method to get the partial atomic charge on a particular atom.
         This is the value printed in the Chargemol analysis.
 
@@ -313,13 +314,13 @@ def get_partial_charge(self, atom_index, charge_type="ddec"):
             atom_index (int): Index of atom to get charge for.
             charge_type (str): Type of charge to use ("ddec" or "cm5").
         """
-        if charge_type.lower() not in ["ddec", "cm5"]:
-            raise ValueError(f"Invalid charge_type: {charge_type}")
         if charge_type.lower() == "ddec":
-            partial_charge = self.ddec_charges[atom_index]
-        elif charge_type.lower() == "cm5":
-            partial_charge = self.cm5_charges[atom_index]
-        return partial_charge
+            return self.ddec_charges[atom_index]
+
+        if charge_type.lower() == "cm5":
+            return self.cm5_charges[atom_index]
+
+        raise ValueError(f"Invalid charge_type: {charge_type}")
 
     def get_bond_order(self, index_from, index_to):
         """Convenience method to get the bond order between two atoms.
@@ -440,13 +441,15 @@ def _get_bond_order_info(filename):
         # Get where relevant info for each atom starts
         bond_order_info = {}
 
-        with open(filename) as r:
+        with open(filename, encoding="utf-8") as r:
+            start_idx = 0
             for line in r:
                 split = line.strip().split()
                 if "Printing BOs" in line:
                     start_idx = int(split[5]) - 1
                     start_el = Element(split[7])
                     bond_order_info[start_idx] = {"element": start_el, "bonded_to": []}
+
                 elif "Bonded to the" in line:
                     direction = tuple(int(i.split(")")[0].split(",")[0]) for i in split[4:7])
                     end_idx = int(split[12]) - 1
@@ -461,6 +464,7 @@ def _get_bond_order_info(filename):
                         "spin_polarization": spin_bo,
                     }
                     bond_order_info[start_idx]["bonded_to"].append(bonded_to)
+
                 elif "The sum of bond orders for this atom" in line:
                     bond_order_info[start_idx]["bond_order_sum"] = float(split[-1])
 
diff --git a/pymatgen/command_line/critic2_caller.py b/pymatgen/command_line/critic2_caller.py
index 569abd29eea..bc12eeabdc4 100644
--- a/pymatgen/command_line/critic2_caller.py
+++ b/pymatgen/command_line/critic2_caller.py
@@ -104,6 +104,7 @@ def __init__(self, input_script: str):
             _stdout, _stderr = rs.communicate()
         stdout = _stdout.decode()
 
+        stderr = ""
         if _stderr:
             stderr = _stderr.decode()
             warnings.warn(stderr)
@@ -791,6 +792,7 @@ def _parse_stdout(self, stdout):
         unique_critical_points = []
 
         # parse unique critical points
+        start_i = end_i = unique_idx = None
         for idx, line in enumerate(stdout):
             if "mult  name            f             |grad|           lap" in line:
                 start_i = idx + 1
diff --git a/pymatgen/command_line/vampire_caller.py b/pymatgen/command_line/vampire_caller.py
index 6f83e25a16f..d7b3ad24d6a 100644
--- a/pymatgen/command_line/vampire_caller.py
+++ b/pymatgen/command_line/vampire_caller.py
@@ -205,8 +205,10 @@ def _create_mat(self):
 
                 if magmoms[site] > 0:
                     spin = 1
-                if magmoms[site] < 0:
+                elif magmoms[site] < 0:
                     spin = -1
+                else:
+                    spin = 0
 
                 atom = structure[i].species.reduced_formula
 
diff --git a/pymatgen/core/interface.py b/pymatgen/core/interface.py
index ba472c4159f..a6862c6d0fb 100644
--- a/pymatgen/core/interface.py
+++ b/pymatgen/core/interface.py
@@ -590,6 +590,7 @@ def gb_from_parameters(
                 else:
                     mu, mv = ratio
                 trans_cry1 = np.array([[1, 0, 0], [-0.5, np.sqrt(3.0) / 2.0, 0], [0, 0, np.sqrt(mu / mv)]])
+
             elif lat_type.lower() == "r":
                 if ratio is None:
                     c2_a2_ratio = 1.0
@@ -603,6 +604,7 @@ def gb_from_parameters(
                         [0, -1 * np.sqrt(3.0) / 3.0, 1.0 / 3 * np.sqrt(c2_a2_ratio)],
                     ]
                 )
+
             else:
                 if lat_type.lower() == "t":
                     if ratio is None:
@@ -610,12 +612,19 @@ def gb_from_parameters(
                     else:
                         mu, mv = ratio
                     lam = mv
+
                 elif lat_type.lower() == "o":
                     new_ratio = [1 if v is None else v for v in ratio]
                     mu, lam, mv = new_ratio
+
+                else:
+                    raise RuntimeError("Invalid lattice type.")
+
                 trans_cry1 = np.array([[1, 0, 0], [0, np.sqrt(lam / mv), 0], [0, 0, np.sqrt(mu / mv)]])
+
         else:
             trans_cry1 = trans_cry
+
         grain_matrix = np.dot(t2, trans_cry1)
         plane_init = np.cross(grain_matrix[0], grain_matrix[1])
         if lat_type.lower() != "c":
@@ -977,6 +986,7 @@ def get_trans_mat(
             index = reduce(gcd, surface)
             surface = [int(round(x / index)) for x in surface]
 
+        lam = None
         if lat_type.lower() == "h":
             # set the value for u,v,w,mu,mv,m,n,d,x
             # check the reference for the meaning of these parameters
@@ -1033,6 +1043,7 @@ def get_trans_mat(
             com_fac = reduce(gcd, all_list)
             sigma = F / com_fac
             r_matrix = (np.array(r_list) / com_fac / sigma).reshape(3, 3)
+
         elif lat_type.lower() == "r":
             # set the value for u,v,w,mu,mv,m,n,d
             # check the reference for the meaning of these parameters
@@ -1107,8 +1118,10 @@ def get_trans_mat(
             com_fac = reduce(gcd, all_list)
             sigma = F / com_fac
             r_matrix = (np.array(r_list) / com_fac / sigma).reshape(3, 3)
+
         else:
             u, v, w = r_axis
+            mu = mv = None
             if lat_type.lower() == "c":
                 mu = 1
                 lam = 1
@@ -1216,6 +1229,7 @@ def get_trans_mat(
         # set one vector of the basis to the rotation axis direction, and
         # obtain the corresponding transform matrix
         eye = np.eye(3, dtype=int)
+        hh = kk = ll = None
         for hh in range(3):
             if abs(r_axis[hh]) != 0:
                 eye[hh] = np.array(r_axis)
@@ -1232,14 +1246,14 @@ def get_trans_mat(
         scale[hh, hh] = 1
         scale[kk, kk] = least_mul
         scale[ll, ll] = sigma / least_mul
+        n_final = None
         for idx in range(least_mul):
             check_int = idx * new_rot[:, kk] + (sigma / least_mul) * new_rot[:, ll]
             if all(np.round(x, 5).is_integer() for x in list(check_int)):
                 n_final = idx
                 break
-        try:
-            n_final  # noqa: B018
-        except NameError:
+
+        if n_final is None:
             raise RuntimeError("Something is wrong. Check if this GB exists or not")
         scale[kk, ll] = n_final
         # each row of mat_csl is the CSL lattice vector
@@ -2031,6 +2045,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals
                 warnings.warn("Too large matrix. Suggest to use quick_gen=False")
             return t_matrix
 
+        c_index = 0
         for ii, jj in enumerate(miller):
             if jj == 0:
                 ab_vector.append(csl[ii])
@@ -2068,6 +2083,8 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals
         if normal:
             c_cross = np.cross(np.matmul(t_matrix[2], trans), np.matmul(surface, ctrans))
             normal_init = np.linalg.norm(c_cross) < 1e-8
+        else:
+            normal_init = False
 
         jj = np.arange(0, max_j + 1)
         combination = []
@@ -2151,6 +2168,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals
                     logger.info("Found perpendicular c vector")
 
         # find the best a, b vectors with their formed area smallest and average norm of a,b smallest.
+        ab_norm = None
         for ii in combinations(ab_vector, 2):
             area_temp = np.linalg.norm(np.cross(np.matmul(ii[0], trans), np.matmul(ii[1], trans)))
             if abs(area_temp - 0) > 1.0e-8:
@@ -2160,6 +2178,7 @@ def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=Fals
                     ab_norm = ab_norm_temp
                     t_matrix[0] = ii[0]
                     t_matrix[1] = ii[1]
+
                 elif area_temp < area or (abs(area - area_temp) < 1.0e-8 and ab_norm_temp < ab_norm):
                     t_matrix[0] = ii[0]
                     t_matrix[1] = ii[1]
diff --git a/pymatgen/core/structure.py b/pymatgen/core/structure.py
index 8fd45120a22..c96ea89f7b4 100644
--- a/pymatgen/core/structure.py
+++ b/pymatgen/core/structure.py
@@ -751,7 +751,7 @@ def _calculate(self, calculator: str | Calculator, verbose: bool = False) -> Cal
 
     def _relax(
         self,
-        calculator: str | Calculator,
+        calculator: Literal["M3GNet", "gfn2-xtb"] | Calculator,
         relax_cell: bool = True,
         optimizer: str | Optimizer = "FIRE",
         steps: int = 500,
@@ -805,6 +805,8 @@ def is_ase_optimizer(key):
             if optimizer not in valid_keys:
                 raise ValueError(f"Unknown {optimizer=}, must be one of {valid_keys}")
             opt_class = getattr(optimize, optimizer)
+        else:
+            opt_class = optimizer
 
         # Get Atoms object
         adaptor = AseAtomsAdaptor()
@@ -843,11 +845,12 @@ def is_ase_optimizer(key):
 
         if return_trajectory:
             if run_uip:
-                traj_observer()  # save properties of the Atoms during the relaxation
+                # Save properties of the Atoms during the relaxation
+                traj_observer()  # type: ignore[reportPossiblyUnboundVariable]
             else:
                 traj_file = opt_kwargs["trajectory"]
                 traj_observer = read(traj_file, index=":")
-            return system, traj_observer
+            return system, traj_observer  # type: ignore[reportPossiblyUnboundVariable]
 
         return system
 
@@ -1672,7 +1675,7 @@ def get_neighbor_list(
     def get_symmetric_neighbor_list(
         self,
         r: float,
-        sg: str,
+        sg: str | None,
         unique: bool = False,
         numerical_tol: float = 1e-8,
         exclude_self: bool = True,
@@ -1716,6 +1719,7 @@ def get_symmetric_neighbor_list(
 
         if sg is None:
             ops = SpaceGroup(self.get_space_group_info()[0]).symmetry_ops
+
         else:
             try:  # first assume sg is int
                 sgp = SpaceGroup.from_int_number(int(sg))
@@ -1723,13 +1727,13 @@ def get_symmetric_neighbor_list(
                 sgp = SpaceGroup(sg)
             ops = sgp.symmetry_ops
 
-        lattice = self.lattice
+            lattice = self.lattice
 
-        if not sgp.is_compatible(lattice):
-            raise ValueError(
-                f"Supplied lattice with parameters {lattice.parameters} is incompatible with "
-                f"supplied spacegroup {sgp.symbol}!"
-            )
+            if not sgp.is_compatible(lattice):
+                raise ValueError(
+                    f"Supplied lattice with parameters {lattice.parameters} is incompatible with "
+                    f"supplied spacegroup {sgp.symbol}!"
+                )
 
         # get a list of neighbors up to distance r
         bonds = self.get_neighbor_list(r)
@@ -2055,7 +2059,7 @@ def get_all_neighbors_old(self, r, include_index=False, include_image=False, inc
                 for i in indices[within_r]:
                     item = []
                     if include_site:
-                        item.append(nnsite)
+                        item.append(nnsite)  # type: ignore[reportPossiblyUnboundVariable]
                     item.append(d[i])
                     if include_index:
                         item.append(j)
@@ -2303,7 +2307,7 @@ def interpolate(
 
         for x in images:
             if interpolate_lattices:
-                l_a = np.dot(np.identity(3) + x * lvec, lstart).T
+                l_a = np.dot(np.identity(3) + x * lvec, lstart).T  # type: ignore[reportPossiblyUnboundVariable]
                 lattice = Lattice(l_a)
             else:
                 lattice = self.lattice
diff --git a/pymatgen/core/surface.py b/pymatgen/core/surface.py
index 03af6ad463d..44e41fdb788 100644
--- a/pymatgen/core/surface.py
+++ b/pymatgen/core/surface.py
@@ -420,6 +420,7 @@ def get_symmetric_site(
 
         # Each operation on a site will return an equivalent site.
         # We want to find the site on the other side of the slab.
+        site_other = None
         for op in ops:
             slab = self.copy()
             site_other = op.operate(point)
@@ -437,6 +438,9 @@ def get_symmetric_site(
             slab.remove_sites([len(slab) - 1])
             slab.remove_sites([len(slab) - 1])
 
+        if site_other is None:
+            raise RuntimeError("Failed to get symmetric site.")
+
         return site_other
 
     def get_orthogonal_c_slab(self) -> Slab:
@@ -686,6 +690,8 @@ def get_equi_sites(slab: Slab, sites: list[int]) -> list[int]:
                 list[int]: Indices of the equivalent sites.
             """
             equi_sites = []
+            eq_indices = []
+            eq_sites: list = []
 
             for pt in sites:
                 # Get the index of the original site
@@ -830,6 +836,8 @@ def get_slab_regions(
     # If slab is noncontiguous
     if frac_coords:
         # Locate the lowest site within the upper Slab
+        last_fcoords = []
+        last_indices = []
         while frac_coords:
             last_fcoords = copy.copy(frac_coords)
             last_indices = copy.copy(indices)
@@ -1690,12 +1698,16 @@ def get_d(slab: Slab) -> float:
     # Sort all sites by z-coordinates
     sorted_sites = sorted(slab, key=lambda site: site.frac_coords[2])
 
+    distance = None
     for site, next_site in zip(sorted_sites, sorted_sites[1:]):
         if not isclose(site.frac_coords[2], next_site.frac_coords[2], abs_tol=1e-6):
             # DEBUG (@DanielYang59): code will break if no distinguishable layers found
             distance = next_site.frac_coords[2] - site.frac_coords[2]
             break
 
+    if distance is None:
+        raise RuntimeError("Cannot identify any layer.")
+
     return slab.lattice.get_cartesian_coords([0, 0, distance])[2]
 
 
@@ -1959,6 +1971,9 @@ def get_symmetrically_equivalent_miller_indices(
     # Convert to hkl if hkil, because in_coord_list only handles tuples of 3
     if len(miller_index) >= 3:
         _miller_index: tuple[int, int, int] = (miller_index[0], miller_index[1], miller_index[-1])
+    else:
+        _miller_index = (miller_index[0], miller_index[1], miller_index[2])
+
     max_idx = max(np.abs(miller_index))
     idx_range = list(range(-max_idx, max_idx + 1))
     idx_range.reverse()
diff --git a/pymatgen/core/tensors.py b/pymatgen/core/tensors.py
index 684cc15520e..68916e30d9e 100644
--- a/pymatgen/core/tensors.py
+++ b/pymatgen/core/tensors.py
@@ -584,6 +584,7 @@ def populate(
             Tensor: Populated tensor
         """
         guess = Tensor(np.zeros(self.shape))
+        mask = None
         if precond:
             # Generate the guess from populated
             sops = SpacegroupAnalyzer(structure).get_symmetry_operations()
diff --git a/pymatgen/core/trajectory.py b/pymatgen/core/trajectory.py
index 966781c7654..abb429087f2 100644
--- a/pymatgen/core/trajectory.py
+++ b/pymatgen/core/trajectory.py
@@ -322,10 +322,14 @@ def __getitem__(self, frames: int | slice | list[int]) -> Molecule | Structure |
                 raise IndexError(f"Frame index {frames} out of range.")
 
             if self.lattice is None:
+                charge = 0
                 if self.charge is not None:
                     charge = int(self.charge)
+
+                spin = None
                 if self.spin_multiplicity is not None:
                     spin = int(self.spin_multiplicity)
+
                 return Molecule(
                     self.species,
                     self.coords[frames],
@@ -544,11 +548,15 @@ def from_file(cls, filename: str | Path, constant_lattice: bool = True, **kwargs
         """
         filename = str(Path(filename).expanduser().resolve())
         is_mol = False
+        molecules = []
+        structures = []
 
         if fnmatch(filename, "*XDATCAR*"):
             structures = Xdatcar(filename).structures
+
         elif fnmatch(filename, "vasprun*.xml*"):
             structures = Vasprun(filename).structures
+
         elif fnmatch(filename, "*.traj"):
             try:
                 from ase.io.trajectory import Trajectory as AseTrajectory
@@ -571,6 +579,7 @@ def from_file(cls, filename: str | Path, constant_lattice: bool = True, **kwargs
 
         if is_mol:
             return cls.from_molecules(molecules, **kwargs)
+
         return cls.from_structures(structures, constant_lattice=constant_lattice, **kwargs)
 
     @staticmethod
diff --git a/pymatgen/electronic_structure/boltztrap.py b/pymatgen/electronic_structure/boltztrap.py
index d2582cb4689..3c16678b7bb 100644
--- a/pymatgen/electronic_structure/boltztrap.py
+++ b/pymatgen/electronic_structure/boltztrap.py
@@ -316,7 +316,7 @@ def write_struct(self, output_file) -> None:
 
         with open(output_file, mode="w") as file:
             if self._symprec is not None:
-                file.write(f"{self._bs.structure.formula} {sym.get_space_group_symbol()}\n")
+                file.write(f"{self._bs.structure.formula} {sym.get_space_group_symbol()}\n")  # type: ignore[reportPossiblyUnboundVariable]
             elif self._symprec is None:
                 file.write(f"{self._bs.structure.formula} symmetries disabled\n")
 
@@ -329,12 +329,12 @@ def write_struct(self, output_file) -> None:
             )
 
             if self._symprec is not None:
-                ops = sym.get_symmetry_dataset()["rotations"]
+                ops = sym.get_symmetry_dataset()["rotations"]  # type: ignore[reportPossiblyUnboundVariable]
             elif self._symprec is None:
                 ops = [np.eye(3)]
-            file.write(f"{len(ops)}\n")
+            file.write(f"{len(ops)}\n")  # type: ignore[reportPossiblyUnboundVariable]
 
-            for op in ops:
+            for op in ops:  # type: ignore[reportPossiblyUnboundVariable]
                 for row in op:
                     file.write(f"{' '.join(map(str, row))}\n")
 
@@ -1291,6 +1291,7 @@ def get_average_eff_mass(self, output="eigs", doping_levels=True):
                             pass
         else:
             result = {t: [] for t in self._seebeck}
+            cond_inv = None
             for temp in result:
                 for i in range(len(self.mu_steps)):
                     try:
@@ -1359,7 +1360,9 @@ def get_seebeck_eff_mass(self, output="average", temp=300, doping_levels=False,
                         )
         return sbk_mass
 
-    def get_complexity_factor(self, output="average", temp=300, doping_levels=False, Lambda=0.5):
+    def get_complexity_factor(
+        self, output: Literal["average", "tensor"] = "average", temp=300, doping_levels=False, Lambda=0.5
+    ):
         """Fermi surface complexity factor respect to calculated as explained in Ref.
         Gibbs, Z. M. et al., Effective mass and fermi surface complexity factor
         from ab initio band structure calculations.
@@ -1394,27 +1397,29 @@ def get_complexity_factor(self, output="average", temp=300, doping_levels=False,
 
                 if output == "average":
                     cmplx_fact[dt] = [(m_s / abs(m_c)) ** 1.5 for m_s, m_c in zip(sbk_mass, cond_mass)]
-                elif output == "tensor":
+
+                else:
                     cmplx_fact[dt] = []
                     for i, sm in enumerate(sbk_mass):
                         cmplx_fact[dt].append([])
                         for j in range(3):
                             cmplx_fact[dt][-1].append((sm[j] / abs(cond_mass[i][j][j])) ** 1.5)
 
-        else:
-            sbk_mass = self.get_seebeck_eff_mass(output, temp, doping_levels=False, Lambda=Lambda)
-            cond_mass = self.get_average_eff_mass(output=output, doping_levels=False)[temp]
-
-            if output == "average":
-                cmplx_fact = [(m_s / abs(m_c)) ** 1.5 for m_s, m_c in zip(sbk_mass, cond_mass)]
-            elif output == "tensor":
-                cmplx_fact = []
-                for i, sm in enumerate(sbk_mass):
-                    cmplx_fact.append([])
-                    for j in range(3):
-                        cmplx_fact[-1].append((sm[j] / abs(cond_mass[i][j][j])) ** 1.5)
+            return cmplx_fact
+
+        sbk_mass = self.get_seebeck_eff_mass(output, temp, doping_levels=False, Lambda=Lambda)
+        cond_mass = self.get_average_eff_mass(output=output, doping_levels=False)[temp]
 
-        return cmplx_fact
+        if output == "average":
+            return [(m_s / abs(m_c)) ** 1.5 for m_s, m_c in zip(sbk_mass, cond_mass)]
+
+        cmplx_fact_list: list = []
+        for i, sm in enumerate(sbk_mass):
+            cmplx_fact_list.append([])
+            for j in range(3):
+                cmplx_fact_list[-1].append((sm[j] / abs(cond_mass[i][j][j])) ** 1.5)
+
+        return cmplx_fact_list
 
     def get_extreme(
         self,
@@ -1592,12 +1597,12 @@ def get_complete_dos(self, structure: Structure, analyzer_for_second_spin=None):
         """
         pdoss: dict[PeriodicSite, dict[Orbital, dict[Spin, ArrayLike]]] = {}
         spin_1 = next(iter(self.dos.densities))
+        spin_2 = next(iter(analyzer_for_second_spin.dos.densities))
 
         if analyzer_for_second_spin:
             if not np.all(self.dos.energies == analyzer_for_second_spin.dos.energies):
                 raise BoltztrapError("Dos merging error: energies of the two dos are different")
 
-            spin_2 = next(iter(analyzer_for_second_spin.dos.densities))
             if spin_1 == spin_2:
                 raise BoltztrapError("Dos merging error: spin component are the same")
 
@@ -1611,6 +1616,7 @@ def get_complete_dos(self, structure: Structure, analyzer_for_second_spin=None):
                 pdoss[structure[idx]][Orbital[o]][spin_1] = self._dos_partial[s][o]
                 if analyzer_for_second_spin:
                     pdoss[structure[idx]][Orbital[o]][spin_2] = analyzer_for_second_spin._dos_partial[s][o]
+
         if analyzer_for_second_spin:
             total_dos = Dos(
                 self.dos.efermi,
@@ -2155,8 +2161,8 @@ def read_cube_file(filename):
     Returns:
         Energy data.
     """
-    with open(filename) as file:
-        n_atoms = 0
+    with open(filename, encoding="utf-8") as file:
+        n_atoms = n1 = n2 = n3 = None
         for idx, line in enumerate(file):
             line = line.rstrip("\n")
             if idx == 0 and "CUBE" not in line:
@@ -2182,9 +2188,13 @@ def read_cube_file(filename):
         n_lines_data = len(energy_data)
         last_line = np.genfromtxt(filename, skip_header=n_lines_data + n_atoms + 6)
         energy_data = np.append(energy_data.flatten(), last_line).reshape(n1, n2, n3)
+
     elif "boltztrap_BZ.cube" in filename:
         energy_data = np.loadtxt(filename, skiprows=n_atoms + 6).reshape(n1, n2, n3)
 
+    else:
+        raise RuntimeError("Cannot read cube file.")
+
     energy_data /= Energy(1, "eV").to("Ry")
 
     return energy_data
diff --git a/pymatgen/electronic_structure/boltztrap2.py b/pymatgen/electronic_structure/boltztrap2.py
index 1b78950cf74..9de856fec01 100644
--- a/pymatgen/electronic_structure/boltztrap2.py
+++ b/pymatgen/electronic_structure/boltztrap2.py
@@ -560,6 +560,7 @@ def get_dos(self, partial_dos=False, npts_mu=10000, T=None, progress=False):
             vvband_ud = [self.vvband]
             spins = [Spin.up]
 
+        energies = []
         for spin, eb, vvb in zip(spins, eband_ud, vvband_ud):
             energies, densities, _vvdos, _cdos = BL.BTPDOS(eb, vvb, npts=npts_mu, erange=enr)
 
@@ -591,14 +592,17 @@ def get_partial_doses(self, tdos, eband_ud, spins, enr, npts_mu, T, progress):
         pdoss = {}
         if progress:
             n_iter = np.prod(np.sum([np.array(i.shape)[2:] for i in self.data.proj.values()]))
-            t = tqdm(total=n_iter * 2)
+            bar = tqdm(total=n_iter * 2)
+        else:
+            bar = None
+
         for spin, eb in zip(spins, eband_ud):
             for idx, site in enumerate(self.data.structure):
                 if site not in pdoss:
                     pdoss[site] = {}
                 for iorb, orb in enumerate(Orbital):
                     if progress:
-                        t.update()
+                        bar.update()
                     if iorb == self.data.proj[spin].shape[-1]:
                         break
 
@@ -1059,6 +1063,7 @@ def plot_props(
         if temps is None:
             temps = self.bzt_transP.temp_r.tolist()
 
+        doping_all = []
         if isinstance(self.bzt_transP.doping, np.ndarray):
             doping_all = self.bzt_transP.doping.tolist()
             if doping is None:
diff --git a/pymatgen/electronic_structure/cohp.py b/pymatgen/electronic_structure/cohp.py
index 4eddad85a6c..21840772dbd 100644
--- a/pymatgen/electronic_structure/cohp.py
+++ b/pymatgen/electronic_structure/cohp.py
@@ -555,6 +555,8 @@ def from_dict(cls, dct: dict) -> Self:
         are_cobis = dct.get("are_cobis", False)
         are_multi_center_cobis = dct.get("are_multi_center_cobis", False)
         are_coops = dct["are_coops"]
+        avg_cohp = None
+
         if "bonds" in dct:
             bonds = {
                 bond: {
diff --git a/pymatgen/electronic_structure/plotter.py b/pymatgen/electronic_structure/plotter.py
index 17416024aed..2be7a76248f 100644
--- a/pymatgen/electronic_structure/plotter.py
+++ b/pymatgen/electronic_structure/plotter.py
@@ -9,7 +9,7 @@
 import typing
 import warnings
 from collections import Counter
-from typing import TYPE_CHECKING, Literal, cast, no_type_check
+from typing import TYPE_CHECKING, cast, no_type_check
 
 import matplotlib.lines as mlines
 import matplotlib.pyplot as plt
@@ -34,11 +34,13 @@
 
 if TYPE_CHECKING:
     from collections.abc import Sequence
+    from typing import Literal
 
     from numpy.typing import ArrayLike
 
     from pymatgen.electronic_structure.dos import CompleteDos, Dos
 
+logger = logging.getLogger(__name__)
 
 __author__ = "Shyue Ping Ong, Geoffroy Hautier, Anubhav Jain"
 __copyright__ = "Copyright 2012, The Materials Project"
@@ -47,8 +49,6 @@
 __email__ = "shyuep@gmail.com"
 __date__ = "May 1, 2012"
 
-logger = logging.getLogger(__name__)
-
 
 class DosPlotter:
     """Class for plotting phonon DOSs. The interface is extremely flexible given there are many
@@ -441,6 +441,9 @@ def bs_plot_data(self, zero_to_efermi=True, bs=None, bs_ref=None, split_branches
         if not bs_is_metal:
             vbm = bs.get_vbm()
             cbm = bs.get_cbm()
+        else:
+            vbm = {}
+            cbm = {}
 
         zero_energy = 0.0
         if zero_to_efermi:
@@ -604,6 +607,8 @@ def get_plot(
 
         handles = []
         vbm_min, cbm_max = [], []
+        data = []
+        one_is_metal = False
 
         colors = next(iter(plt.rcParams["axes.prop_cycle"].by_key().values()))
         for ibs, bs in enumerate(self._bs):
@@ -617,9 +622,7 @@ def get_plot(
                 data = self.bs_plot_data(zero_to_efermi, bs, bs_ref, split_branches=False)
 
             # remember if one bs is a metal for setting the ylim later
-            one_is_metal = False
-            if not one_is_metal and data["is_metal"]:
-                one_is_metal = data["is_metal"]
+            one_is_metal = data["is_metal"]
 
             # remember all the cbm and vbm for setting the ylim later
             if not data["is_metal"]:
@@ -1709,6 +1712,7 @@ def _Orbitals_SumOrbitals(cls, dictio, sum_morbs):
         if len(dictio) == 0:
             raise KeyError("The 'dictio' is empty. We cannot do anything.")
 
+        orb = None
         for elt in dictio:
             if Element.is_valid_symbol(elt):
                 if isinstance(dictio[elt], list):
@@ -2318,8 +2322,7 @@ def get_plot(
         fig = plt.figure(figsize=self.fig_size)
         fig.patch.set_facecolor("white")
         bs_ax = plt.subplot(gs[0])
-        if dos:
-            dos_ax = plt.subplot(gs[1])
+        dos_ax = plt.subplot(gs[1]) if dos else None
 
         # set basic axes limits for the plot
         bs_ax.set_xlim(0, x_distances_list[-1][-1])
@@ -2354,8 +2357,7 @@ def get_plot(
                     band_energies[spin].append([e - bs.efermi for e in band])  # type: ignore
 
         # renormalize the DOS energies to Fermi level
-        if dos:
-            dos_energies = [e - dos.efermi for e in dos.energies]
+        dos_energies = [e - dos.efermi for e in dos.energies] if dos else []
 
         # get the projection data to set colors for the band structure
         colordata = self._get_colordata(bs, elements, bs_projection)
@@ -2534,6 +2536,7 @@ def _get_colordata(bs, elements, bs_projection):
             Dictionary representation of color data.
         """
         contribs = {}
+        projections = None
         if bs_projection and bs_projection.lower() == "elements":
             projections = bs.get_projection_on_elements()
 
@@ -3037,7 +3040,7 @@ def plot_zt_mu(
         plt.tight_layout()
         return ax
 
-    def plot_seebeck_temp(self, doping="all", output="average"):
+    def plot_seebeck_temp(self, doping="all", output: Literal["average", "eigs"] = "average"):
         """Plot the Seebeck coefficient in function of temperature for different
         doping levels.
 
@@ -3050,10 +3053,7 @@ def plot_seebeck_temp(self, doping="all", output="average"):
         Returns:
             a matplotlib object
         """
-        if output == "average":
-            sbk = self._bz.get_seebeck(output="average")
-        elif output == "eigs":
-            sbk = self._bz.get_seebeck(output="eigs")
+        sbk = self._bz.get_seebeck(output="average") if output == "average" else self._bz.get_seebeck(output="eigs")
 
         ax = pretty_plot(22, 14)
         tlist = sorted(sbk["n"])
@@ -3089,7 +3089,9 @@ def plot_seebeck_temp(self, doping="all", output="average"):
 
         return ax
 
-    def plot_conductivity_temp(self, doping="all", output="average", relaxation_time=1e-14):
+    def plot_conductivity_temp(
+        self, doping="all", output: Literal["average", "eigs"] = "average", relaxation_time=1e-14
+    ):
         """Plot the conductivity in function of temperature for different doping levels.
 
         Args:
@@ -3104,7 +3106,7 @@ def plot_conductivity_temp(self, doping="all", output="average", relaxation_time
         """
         if output == "average":
             cond = self._bz.get_conductivity(relaxation_time=relaxation_time, output="average")
-        elif output == "eigs":
+        else:
             cond = self._bz.get_conductivity(relaxation_time=relaxation_time, output="eigs")
 
         ax = pretty_plot(22, 14)
@@ -3142,7 +3144,9 @@ def plot_conductivity_temp(self, doping="all", output="average", relaxation_time
 
         return ax
 
-    def plot_power_factor_temp(self, doping="all", output="average", relaxation_time=1e-14):
+    def plot_power_factor_temp(
+        self, doping="all", output: Literal["average", "eigs"] = "average", relaxation_time=1e-14
+    ):
         """Plot the Power Factor in function of temperature for different doping levels.
 
         Args:
@@ -3157,7 +3161,7 @@ def plot_power_factor_temp(self, doping="all", output="average", relaxation_time
         """
         if output == "average":
             pow_factor = self._bz.get_power_factor(relaxation_time=relaxation_time, output="average")
-        elif output == "eigs":
+        else:
             pow_factor = self._bz.get_power_factor(relaxation_time=relaxation_time, output="eigs")
 
         ax = pretty_plot(22, 14)
@@ -3210,7 +3214,7 @@ def plot_zt_temp(self, doping="all", output: Literal["average", "eigs"] = "avera
         Returns:
             a matplotlib object
         """
-        if output not in ("average", "eigs"):
+        if output not in {"average", "eigs"}:
             raise ValueError(f"{output=} must be 'average' or 'eigs'")
         zt = self._bz.get_zt(relaxation_time=relaxation_time, output=output)
 
@@ -3262,7 +3266,7 @@ def plot_eff_mass_temp(self, doping="all", output: Literal["average", "eigs"] =
         """
         if output == "average":
             eff_mass = self._bz.get_average_eff_mass(output="average")
-        elif output == "eigs":
+        else:
             eff_mass = self._bz.get_average_eff_mass(output="eigs")
 
         ax_main = pretty_plot(22, 14)
@@ -3290,7 +3294,7 @@ def plot_eff_mass_temp(self, doping="all", output: Literal["average", "eigs"] =
         plt.tight_layout()
         return ax_main
 
-    def plot_seebeck_dop(self, temps="all", output="average"):
+    def plot_seebeck_dop(self, temps="all", output: Literal["average", "eigs"] = "average"):
         """Plot the Seebeck in function of doping levels for different temperatures.
 
         Args:
@@ -3302,10 +3306,7 @@ def plot_seebeck_dop(self, temps="all", output="average"):
         Returns:
             a matplotlib object
         """
-        if output == "average":
-            sbk = self._bz.get_seebeck(output="average")
-        elif output == "eigs":
-            sbk = self._bz.get_seebeck(output="eigs")
+        sbk = self._bz.get_seebeck(output="average") if output == "average" else self._bz.get_seebeck(output="eigs")
 
         tlist = sorted(sbk["n"]) if temps == "all" else temps
         ax = pretty_plot(22, 14)
@@ -3334,7 +3335,7 @@ def plot_seebeck_dop(self, temps="all", output="average"):
 
         return ax
 
-    def plot_conductivity_dop(self, temps="all", output="average", relaxation_time=1e-14):
+    def plot_conductivity_dop(self, temps="all", output: Literal["average", "eigs"] = "average", relaxation_time=1e-14):
         """Plot the conductivity in function of doping levels for different
         temperatures.
 
@@ -3350,7 +3351,7 @@ def plot_conductivity_dop(self, temps="all", output="average", relaxation_time=1
         """
         if output == "average":
             cond = self._bz.get_conductivity(relaxation_time=relaxation_time, output="average")
-        elif output == "eigs":
+        else:
             cond = self._bz.get_conductivity(relaxation_time=relaxation_time, output="eigs")
 
         tlist = sorted(cond["n"]) if temps == "all" else temps
@@ -3384,7 +3385,7 @@ def plot_conductivity_dop(self, temps="all", output="average", relaxation_time=1
 
         return ax
 
-    def plot_power_factor_dop(self, temps="all", output="average", relaxation_time=1e-14):
+    def plot_power_factor_dop(self, temps="all", output: Literal["average", "eigs"] = "average", relaxation_time=1e-14):
         """Plot the Power Factor in function of doping levels for different temperatures.
 
         Args:
@@ -3399,7 +3400,7 @@ def plot_power_factor_dop(self, temps="all", output="average", relaxation_time=1
         """
         if output == "average":
             pow_factor = self._bz.get_power_factor(relaxation_time=relaxation_time, output="average")
-        elif output == "eigs":
+        else:
             pow_factor = self._bz.get_power_factor(relaxation_time=relaxation_time, output="eigs")
 
         tlist = sorted(pow_factor["n"]) if temps == "all" else temps
@@ -3432,7 +3433,7 @@ def plot_power_factor_dop(self, temps="all", output="average", relaxation_time=1
 
         return ax
 
-    def plot_zt_dop(self, temps="all", output="average", relaxation_time=1e-14):
+    def plot_zt_dop(self, temps="all", output: Literal["average", "eigs"] = "average", relaxation_time=1e-14):
         """Plot the figure of merit zT in function of doping levels for different
         temperatures.
 
@@ -3448,7 +3449,7 @@ def plot_zt_dop(self, temps="all", output="average", relaxation_time=1e-14):
         """
         if output == "average":
             zt = self._bz.get_zt(relaxation_time=relaxation_time, output="average")
-        elif output == "eigs":
+        else:
             zt = self._bz.get_zt(relaxation_time=relaxation_time, output="eigs")
 
         tlist = sorted(zt["n"]) if temps == "all" else temps
@@ -3483,7 +3484,7 @@ def plot_zt_dop(self, temps="all", output="average", relaxation_time=1e-14):
 
         return ax
 
-    def plot_eff_mass_dop(self, temps="all", output="average"):
+    def plot_eff_mass_dop(self, temps="all", output: Literal["average", "eigs"] = "average"):
         """Plot the average effective mass in function of doping levels
         for different temperatures.
 
@@ -3499,7 +3500,7 @@ def plot_eff_mass_dop(self, temps="all", output="average"):
         """
         if output == "average":
             em = self._bz.get_average_eff_mass(output="average")
-        elif output == "eigs":
+        else:
             em = self._bz.get_average_eff_mass(output="eigs")
 
         tlist = sorted(em["n"]) if temps == "all" else temps
@@ -3714,7 +3715,9 @@ def get_plot(
 
         allpts = []
         keys = list(self._cohps)
-        for i, key in enumerate(keys):
+        idx = key = None
+
+        for idx, key in enumerate(keys):
             energies = self._cohps[key]["energies"]
             populations = self._cohps[key]["COHP"] if not integrated else self._cohps[key]["ICOHP"]
             for spin in [Spin.up, Spin.down]:
@@ -3730,13 +3733,13 @@ def get_plot(
                         ax.plot(
                             x,
                             y,
-                            color=colors[i % ncolors],
+                            color=colors[idx % ncolors],
                             linestyle="-",
                             label=str(key),
                             linewidth=3,
                         )
                     else:
-                        ax.plot(x, y, color=colors[i % ncolors], linestyle="--", linewidth=3)
+                        ax.plot(x, y, color=colors[idx % ncolors], linestyle="--", linewidth=3)
 
         if xlim:
             ax.set_xlim(xlim)
@@ -3761,7 +3764,7 @@ def get_plot(
                 ax.plot(
                     [self._cohps[key]["efermi"], self._cohps[key]["efermi"]],
                     ylim,
-                    color=colors[i % ncolors],
+                    color=colors[idx % ncolors],
                     linestyle="--",
                     linewidth=2,
                 )
@@ -3773,7 +3776,7 @@ def get_plot(
                 ax.plot(
                     xlim,
                     [self._cohps[key]["efermi"], self._cohps[key]["efermi"]],
-                    color=colors[i % ncolors],
+                    color=colors[idx % ncolors],
                     linestyle="--",
                     linewidth=2,
                 )
@@ -3905,8 +3908,7 @@ def plot_fermi_surface(
     if transparency_factor is None:
         transparency_factor = [1] * n_surfaces
 
-    if mlab_figure:
-        fig = mlab_figure
+    fig = mlab_figure if mlab_figure else None
 
     if kpoints_dict is None:
         kpoints_dict = {}
diff --git a/pymatgen/ext/cod.py b/pymatgen/ext/cod.py
index 6c25f7d1d80..190b9940488 100644
--- a/pymatgen/ext/cod.py
+++ b/pymatgen/ext/cod.py
@@ -57,7 +57,7 @@ def query(self, sql: str) -> str:
         return response.decode("utf-8")
 
     @requires(which("mysql"), "mysql must be installed to use this query.")
-    def get_cod_ids(self, formula):
+    def get_cod_ids(self, formula) -> list[int]:
         """Queries the COD for all cod ids associated with a formula. Requires
         mysql executable to be in the path.
 
@@ -80,17 +80,18 @@ def get_cod_ids(self, formula):
                 cod_ids.append(int(match.group(1)))
         return cod_ids
 
-    def get_structure_by_id(self, cod_id, **kwargs):
+    def get_structure_by_id(self, cod_id: int, timeout: int = 600, **kwargs) -> Structure:
         """Queries the COD for a structure by id.
 
         Args:
             cod_id (int): COD id.
+            timeout (int): Timeout for the request in seconds. Default = 600.
             kwargs: All kwargs supported by Structure.from_str.
 
         Returns:
             A Structure.
         """
-        response = requests.get(f"http://{self.url}/cod/{cod_id}.cif", timeout=60)
+        response = requests.get(f"http://{self.url}/cod/{cod_id}.cif", timeout=timeout)
         return Structure.from_str(response.text, fmt="cif", **kwargs)
 
     @requires(which("mysql"), "mysql must be installed to use this query.")
@@ -112,7 +113,7 @@ def get_structure_by_formula(self, formula: str, **kwargs) -> list[dict[str, str
         for line in text:
             if line.strip():
                 cod_id, sg = line.split("\t")
-                response = requests.get(f"http://www.crystallography.net/cod/{cod_id.strip()}.cif", timeout=60)
+                response = requests.get(f"http://www.crystallography.net/cod/{cod_id.strip()}.cif", timeout=600)
                 try:
                     struct = Structure.from_str(response.text, fmt="cif", **kwargs)
                     structures.append({"structure": struct, "cod_id": int(cod_id), "sg": sg})
diff --git a/pymatgen/ext/matproj_legacy.py b/pymatgen/ext/matproj_legacy.py
index 68ad6d2ee5f..f720c5984cc 100644
--- a/pymatgen/ext/matproj_legacy.py
+++ b/pymatgen/ext/matproj_legacy.py
@@ -1570,7 +1570,7 @@ def _check_get_download_info_url_by_task_id(self, prefix, task_ids) -> list[str]
 
     @staticmethod
     def _check_nomad_exist(url) -> bool:
-        response = requests.get(url=url, timeout=60)
+        response = requests.get(url=url, timeout=600)
         if response.status_code != 200:
             return False
         content = json.loads(response.text)
diff --git a/pymatgen/io/vasp/help.py b/pymatgen/io/vasp/help.py
index 376e114d0a8..bd6b3cbad49 100644
--- a/pymatgen/io/vasp/help.py
+++ b/pymatgen/io/vasp/help.py
@@ -48,7 +48,7 @@ def get_help(cls, tag: str, fmt: str = "text") -> str:
             Help text.
         """
         tag = tag.upper()
-        response = requests.get(f"https://www.vasp.at/wiki/index.php/{tag}", verify=False, timeout=60)
+        response = requests.get(f"https://www.vasp.at/wiki/index.php/{tag}", verify=False, timeout=600)
         soup = BeautifulSoup(response.text)
         main_doc = soup.find(id="mw-content-text")
         if fmt == "text":
@@ -67,7 +67,7 @@ def get_incar_tags(cls) -> list[str]:
             "https://www.vasp.at/wiki/index.php/Category:INCAR",
             "https://www.vasp.at/wiki/index.php?title=Category:INCAR&pagefrom=ML+FF+LCONF+DISCARD#mw-pages",
         ]:
-            response = requests.get(page, verify=False, timeout=60)
+            response = requests.get(page, verify=False, timeout=600)
             soup = BeautifulSoup(response.text)
             for div in soup.findAll("div", {"class": "mw-category-group"}):
                 children = div.findChildren("li")
diff --git a/tasks.py b/tasks.py
index f9c9f99f9e9..a5ad2e06a1e 100644
--- a/tasks.py
+++ b/tasks.py
@@ -178,7 +178,7 @@ def update_changelog(ctx: Context, version: str | None = None, dry_run: bool = F
             pr_number = re_match.group(1)
             contributor, pr_name = re_match.group(2).split("/", 1)
             response = requests.get(
-                f"https://api.github.com/repos/materialsproject/pymatgen/pulls/{pr_number}", timeout=60
+                f"https://api.github.com/repos/materialsproject/pymatgen/pulls/{pr_number}", timeout=600
             )
             lines += [f"* PR #{pr_number} from @{contributor} {pr_name}"]
             json_resp = response.json()
diff --git a/tests/ext/test_cod.py b/tests/ext/test_cod.py
index 6616205909e..9a3164d03d0 100644
--- a/tests/ext/test_cod.py
+++ b/tests/ext/test_cod.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import os
 from shutil import which
 from unittest import TestCase
 
@@ -8,25 +9,29 @@
 
 from pymatgen.ext.cod import COD
 
+if "CI" in os.environ:  # test is slow and flaky, skip in CI. see
+    # https://github.com/materialsproject/pymatgen/pull/3777#issuecomment-2071217785
+    pytest.skip(allow_module_level=True)
+
 try:
-    website_down = requests.get("https://www.crystallography.net", timeout=60).status_code != 200
+    website_down = requests.get("https://www.crystallography.net", timeout=600).status_code != 200
 except requests.exceptions.ConnectionError:
     website_down = True
 
 
-@pytest.mark.skipif(website_down, reason="www.crystallography.net is down.")
+@pytest.mark.skipif(website_down, reason="www.crystallography.net is down")
 class TestCOD(TestCase):
-    @pytest.mark.skipif(not which("mysql"), reason="No mysql.")
+    @pytest.mark.skipif(not which("mysql"), reason="No mysql")
     def test_get_cod_ids(self):
         ids = COD().get_cod_ids("Li2O")
         assert len(ids) > 15
 
-    @pytest.mark.skipif(not which("mysql"), reason="No mysql.")
+    @pytest.mark.skipif(not which("mysql"), reason="No mysql")
     def test_get_structure_by_formula(self):
         data = COD().get_structure_by_formula("Li2O")
         assert len(data) > 15
         assert data[0]["structure"].reduced_formula == "Li2O"
 
     def test_get_structure_by_id(self):
-        struct = COD().get_structure_by_id(2002926)
+        struct = COD().get_structure_by_id(2_002_926)
         assert struct.formula == "Be8 H64 N16 F32"
diff --git a/tests/ext/test_matproj.py b/tests/ext/test_matproj.py
index 2d4961518c9..ec5ee928620 100644
--- a/tests/ext/test_matproj.py
+++ b/tests/ext/test_matproj.py
@@ -26,7 +26,7 @@
 from pymatgen.util.testing import TEST_FILES_DIR, PymatgenTest
 
 try:
-    skip_mprester_tests = requests.get("https://materialsproject.org", timeout=60).status_code != 200
+    skip_mprester_tests = requests.get("https://materialsproject.org", timeout=600).status_code != 200
 
 except (ModuleNotFoundError, ImportError, requests.exceptions.ConnectionError):
     # Skip all MPRester tests if some downstream problem on the website, mp-api or whatever.
@@ -76,7 +76,7 @@ def test_get_data(self):
             "total_magnetization",
         }
         mp_id = "mp-1143"
-        vals = requests.get(f"http://legacy.materialsproject.org/materials/{mp_id}/json/", timeout=60)
+        vals = requests.get(f"http://legacy.materialsproject.org/materials/{mp_id}/json/", timeout=600)
         expected_vals = vals.json()
 
         for prop in props:
@@ -570,7 +570,7 @@ def test_get_all_materials_ids_doc(self):
     #         "total_magnetization",
     #     }
     #     mp_id = "mp-1143"
-    #     vals = requests.get(f"http://legacy.materialsproject.org/materials/{mp_id}/json/", timeout=60)
+    #     vals = requests.get(f"http://legacy.materialsproject.org/materials/{mp_id}/json/", timeout=600)
     #     expected_vals = vals.json()
     #
     #     for prop in props:
diff --git a/tests/ext/test_optimade.py b/tests/ext/test_optimade.py
index 564694b48f0..ea351acb8c7 100644
--- a/tests/ext/test_optimade.py
+++ b/tests/ext/test_optimade.py
@@ -13,12 +13,12 @@
     website_down = True
 
 try:
-    optimade_providers_down = requests.get("https://providers.optimade.org", timeout=60).status_code not in (200, 403)
+    optimade_providers_down = requests.get("https://providers.optimade.org", timeout=600).status_code not in (200, 403)
 except requests.exceptions.ConnectionError:
     optimade_providers_down = True
 
 try:
-    mc3d_down = requests.get(OptimadeRester.aliases["mcloud.mc3d"] + "/v1/info", timeout=60).status_code not in (
+    mc3d_down = requests.get(OptimadeRester.aliases["mcloud.mc3d"] + "/v1/info", timeout=600).status_code not in (
         200,
         403,
         301,
@@ -27,7 +27,7 @@
     mc3d_down = True
 
 try:
-    mc2d_down = requests.get(OptimadeRester.aliases["mcloud.mc2d"] + "/v1/info", timeout=60).status_code not in (
+    mc2d_down = requests.get(OptimadeRester.aliases["mcloud.mc2d"] + "/v1/info", timeout=600).status_code not in (
         200,
         403,
         301,