Use strict=True with zip to ensure length equality (#4011)

* global replacement of strict=True, very prone to error and need careful check

* revert changes on non-zip functions

* revert to strict=F for failing tests

* revert a missing one

* remove a call-overload ignore tag

* fix another unit test

* fix typos

---------

Co-authored-by: Janosh Riebesell <janosh.riebesell@gmail.com>

dev_scripts/update_pt_data.py

@@ -219,7 +219,7 @@ def gen_iupac_ordering():
 
     order = sum((list(product(x, y)) for x, y in order), [])  # noqa: RUF017
     iupac_ordering_dict = dict(
-        zip([Element.from_row_and_group(row, group) for group, row in order], range(len(order)), strict=False)
+        zip([Element.from_row_and_group(row, group) for group, row in order], range(len(order)), strict=True)
     )
 
     # first clean periodic table of any IUPAC ordering

src/pymatgen/analysis/adsorption.py

@@ -178,7 +178,7 @@ def find_surface_sites_by_height(self, slab: Slab, height=0.9, xy_tol=0.05):
         surf_sites = [slab.sites[n] for n in np.where(mask)[0]]
         if xy_tol:
             # sort surface sites by height
-            surf_sites = [s for (h, s) in zip(m_projs[mask], surf_sites, strict=False)]
+            surf_sites = [s for (h, s) in zip(m_projs[mask], surf_sites, strict=True)]
             surf_sites.reverse()
             unique_sites: list = []
             unique_perp_fracs: list = []
@@ -235,7 +235,7 @@ def find_adsorption_sites(
 
         Args:
             distance (float): distance from the coordinating ensemble
-                of atoms along the miller index for the site (i. e.
+                of atoms along the miller index for the site (i.e.
                 the distance from the slab itself)
             put_inside (bool): whether to put the site inside the cell
             symm_reduce (float): symm reduction threshold
@@ -268,7 +268,7 @@ def find_adsorption_sites(
             for v in dt.simplices:
                 if -1 not in v:
                     dots = []
-                    for i_corner, i_opp in zip(range(3), ((1, 2), (0, 2), (0, 1)), strict=False):
+                    for i_corner, i_opp in zip(range(3), ((1, 2), (0, 2), (0, 1)), strict=True):
                         corner, opp = v[i_corner], [v[o] for o in i_opp]
                         vecs = [mesh[d].coords - mesh[corner].coords for d in opp]
                         vecs = [vec / np.linalg.norm(vec) for vec in vecs]
@@ -701,7 +701,7 @@ def plot_slab(
         ads_sites = asf.find_adsorption_sites()["all"]
         symm_op = get_rot(orig_slab)
         ads_sites = [symm_op.operate(ads_site)[:2].tolist() for ads_site in ads_sites]
-        ax.plot(*zip(*ads_sites, strict=False), color="k", marker="x", markersize=10, mew=1, linestyle="", zorder=10000)
+        ax.plot(*zip(*ads_sites, strict=True), color="k", marker="x", markersize=10, mew=1, linestyle="", zorder=10000)
     # Draw unit cell
     if draw_unit_cell:
         vertices = np.insert(vertices, 1, lattice_sum, axis=0).tolist()

src/pymatgen/analysis/bond_valence.py

@@ -403,7 +403,7 @@ def _recurse(assigned=None):
         if self._best_vset:
             if structure.is_ordered:
                 assigned = {}
-                for val, sites in zip(self._best_vset, equi_sites, strict=False):
+                for val, sites in zip(self._best_vset, equi_sites, strict=True):
                     for site in sites:
                         assigned[site] = val
 
@@ -414,7 +414,7 @@ def _recurse(assigned=None):
                 new_best_vset.append([])
             for ival, val in enumerate(self._best_vset):
                 new_best_vset[attrib[ival]].append(val)
-            for val, sites in zip(new_best_vset, equi_sites, strict=False):
+            for val, sites in zip(new_best_vset, equi_sites, strict=True):
                 for site in sites:
                     assigned[site] = val
 

src/pymatgen/analysis/chempot_diagram.py

@@ -213,7 +213,7 @@ def _get_domains(self) -> dict[str, np.ndarray]:
 
         domains: dict[str, list] = {entry.reduced_formula: [] for entry in entries}
 
-        for intersection, facet in zip(hs_int.intersections, hs_int.dual_facets, strict=False):
+        for intersection, facet in zip(hs_int.intersections, hs_int.dual_facets, strict=True):
             for v in facet:
                 if v < len(entries):
                     this_entry = entries[v]
@@ -578,7 +578,7 @@ def get_chempot_axis_title(element) -> str:
             return f"<br> μ<sub>{element}</sub> - μ<sub>{element}</sub><sup>o</sup> (eV)"
 
         axes_layout = {}
-        for ax, el in zip(axes, elements, strict=False):
+        for ax, el in zip(axes, elements, strict=True):
             layout = plotly_layouts[layout_name].copy()
             layout["title"] = get_chempot_axis_title(el)
             axes_layout[ax] = layout

src/pymatgen/analysis/diffraction/core.py

@@ -105,7 +105,7 @@ def get_plot(
         xrd = self.get_pattern(structure, two_theta_range=two_theta_range)
         imax = max(xrd.y)
 
-        for two_theta, i, hkls in zip(xrd.x, xrd.y, xrd.hkls, strict=False):
+        for two_theta, i, hkls in zip(xrd.x, xrd.y, xrd.hkls, strict=True):
             if two_theta_range[0] <= two_theta <= two_theta_range[1]:
                 hkl_tuples = [hkl["hkl"] for hkl in hkls]
                 label = ", ".join(map(str, hkl_tuples))  # 'full' label
@@ -188,7 +188,7 @@ def plot_structures(self, structures, fontsize=6, **kwargs):
         n_rows = len(structures)
         fig, axes = plt.subplots(nrows=n_rows, ncols=1, sharex=True, squeeze=False)
 
-        for i, (ax, structure) in enumerate(zip(axes.ravel(), structures, strict=False)):
+        for i, (ax, structure) in enumerate(zip(axes.ravel(), structures, strict=True)):
             self.get_plot(structure, fontsize=fontsize, ax=ax, with_labels=i == n_rows - 1, **kwargs)
             spg_symbol, spg_number = structure.get_space_group_info()
             ax.set_title(f"{structure.formula} {spg_symbol} ({spg_number}) ")

src/pymatgen/analysis/diffraction/tem.py

@@ -139,7 +139,7 @@ def get_interplanar_spacings(
         if (0, 0, 0) in points_filtered:
             points_filtered.remove((0, 0, 0))
         interplanar_spacings_val = np.array([structure.lattice.d_hkl(x) for x in points_filtered])
-        return dict(zip(points_filtered, interplanar_spacings_val, strict=False))
+        return dict(zip(points_filtered, interplanar_spacings_val, strict=True))
 
     def bragg_angles(self, interplanar_spacings: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, float]:
         """Get the Bragg angles for every hkl point passed in (where n = 1).
@@ -153,7 +153,7 @@ def bragg_angles(self, interplanar_spacings: dict[Tuple3Ints, float]) -> dict[Tu
         plane = list(interplanar_spacings)
         interplanar_spacings_val = np.array(list(interplanar_spacings.values()))
         bragg_angles_val = np.arcsin(self.wavelength_rel() / (2 * interplanar_spacings_val))
-        return dict(zip(plane, bragg_angles_val, strict=False))
+        return dict(zip(plane, bragg_angles_val, strict=True))
 
     def get_s2(self, bragg_angles: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, float]:
         """
@@ -169,7 +169,7 @@ def get_s2(self, bragg_angles: dict[Tuple3Ints, float]) -> dict[Tuple3Ints, floa
         plane = list(bragg_angles)
         bragg_angles_val = np.array(list(bragg_angles.values()))
         s2_val = (np.sin(bragg_angles_val) / self.wavelength_rel()) ** 2
-        return dict(zip(plane, s2_val, strict=False))
+        return dict(zip(plane, s2_val, strict=True))
 
     def x_ray_factors(
         self, structure: Structure, bragg_angles: dict[Tuple3Ints, float]
@@ -269,7 +269,7 @@ def cell_intensity(self, structure: Structure, bragg_angles: dict[Tuple3Ints, fl
         csf = self.cell_scattering_factors(structure, bragg_angles)
         csf_val = np.array(list(csf.values()))
         cell_intensity_val = (csf_val * csf_val.conjugate()).real
-        return dict(zip(bragg_angles, cell_intensity_val, strict=False))
+        return dict(zip(bragg_angles, cell_intensity_val, strict=True))
 
     def get_pattern(
         self,

src/pymatgen/analysis/elasticity/elastic.py

@@ -137,7 +137,7 @@ class ElasticTensor(NthOrderElasticTensor):
     def __new__(cls, input_array, tol: float = 1e-4) -> Self:
         """
         Create an ElasticTensor object. The constructor throws an error if the shape of
-        the input_matrix argument is not 3x3x3x3, i. e. in true tensor notation. Issues a
+        the input_matrix argument is not 3x3x3x3, i.e. in true tensor notation. Issues a
         warning if the input_matrix argument does not satisfy standard symmetries. Note
         that the constructor uses __new__ rather than __init__ according to the standard
         method of subclassing numpy ndarrays.
@@ -564,7 +564,7 @@ def from_diff_fit(cls, strains, stresses, eq_stress=None, tol: float = 1e-10, or
     @property
     def order(self) -> int:
         """
-        Order of the elastic tensor expansion, i. e. the order of the
+        Order of the elastic tensor expansion, i.e. the order of the
         highest included set of elastic constants.
         """
         return self[-1].order
@@ -619,7 +619,7 @@ def get_tgt(self, temperature: float | None = None, structure: Structure = None,
         points = quad["points"]
         weights = quad["weights"]
         num, denom, c = np.zeros((3, 3)), 0, 1
-        for p, w in zip(points, weights, strict=False):
+        for p, w in zip(points, weights, strict=True):
             gk = ElasticTensor(self[0]).green_kristoffel(p)
             _rho_wsquareds, us = np.linalg.eigh(gk)
             us = [u / np.linalg.norm(u) for u in np.transpose(us)]
@@ -856,7 +856,7 @@ def diff_fit(strains, stresses, eq_stress=None, order=2, tol: float = 1e-10):
         stresses (nx3x3 array-like): Array of 3x3 stresses
             to use in fitting ECs. These should be PK2 stresses.
         eq_stress (3x3 array-like): stress corresponding to
-            equilibrium strain (i. e. "0" strain state).
+            equilibrium strain (i.e. "0" strain state).
             If not specified, function will try to find
             the state in the list of provided stresses
             and strains. If not found, defaults to 0.
@@ -882,7 +882,7 @@ def diff_fit(strains, stresses, eq_stress=None, order=2, tol: float = 1e-10):
     for _ord in range(1, order):
         cvec, carr = get_symbol_list(_ord + 1)
         svec = np.ravel(dei_dsi[_ord - 1].T)
-        cmap = dict(zip(cvec, np.dot(m[_ord - 1], svec), strict=False))
+        cmap = dict(zip(cvec, np.dot(m[_ord - 1], svec), strict=True))
         c_list.append(v_subs(carr, cmap))
     return [Tensor.from_voigt(c) for c in c_list]
 
@@ -916,7 +916,7 @@ def find_eq_stress(strains, stresses, tol: float = 1e-10):
 
 def get_strain_state_dict(strains, stresses, eq_stress=None, tol: float = 1e-10, add_eq=True, sort=True):
     """Create a dictionary of voigt notation stress-strain sets
-    keyed by "strain state", i. e. a tuple corresponding to
+    keyed by "strain state", i.e. a tuple corresponding to
     the non-zero entries in ratios to the lowest nonzero value,
     e.g. [0, 0.1, 0, 0.2, 0, 0] -> (0,1,0,2,0,0)
     This allows strains to be collected in stencils as to
@@ -974,7 +974,7 @@ def generate_pseudo(strain_states, order=3):
 
     Args:
         strain_states (6xN array like): a list of Voigt-notation strain-states,
-            i. e. perturbed indices of the strain as a function of the smallest
+            i.e. perturbed indices of the strain as a function of the smallest
             strain e.g. (0, 1, 0, 0, 1, 0)
         order (int): order of pseudo-inverse to calculate
 
@@ -1010,7 +1010,7 @@ def generate_pseudo(strain_states, order=3):
 def get_symbol_list(rank, dim=6):
     """Get a symbolic representation of the Voigt-notation
     tensor that places identical symbols for entries related
-    by index transposition, i. e. C_1121 = C_1211 etc.
+    by index transposition, i.e. C_1121 = C_1211 etc.
 
     Args:
         dim (int): dimension of matrix/tensor, e.g. 6 for

src/pymatgen/analysis/elasticity/strain.py

@@ -43,7 +43,7 @@ class Deformation(SquareTensor):
     def __new__(cls, deformation_gradient) -> Self:
         """
         Create a Deformation object. Note that the constructor uses __new__ rather than
-        __init__ according to the standard method of subclassing numpy ndarrays.
+        __init__ according to the standard method of subclassing numpy arrays.
 
         Args:
             deformation_gradient (3x3 array-like): the 3x3 array-like
@@ -56,11 +56,11 @@ def is_independent(self, tol: float = 1e-8):
         """Check to determine whether the deformation is independent."""
         return len(self.get_perturbed_indices(tol)) == 1
 
-    def get_perturbed_indices(self, tol: float = 1e-8):
+    def get_perturbed_indices(self, tol: float = 1e-8) -> list[tuple[int, int]]:
         """Get indices of perturbed elements of the deformation gradient,
-        i. e. those that differ from the identity.
+        i.e. those that differ from the identity.
         """
-        return list(zip(*np.where(abs(self - np.eye(3)) > tol), strict=False))
+        return list(zip(*np.where(abs(self - np.eye(3)) > tol), strict=True))
 
     @property
     def green_lagrange_strain(self):

src/pymatgen/analysis/eos.py

@@ -431,7 +431,7 @@ def get_rms(x, y):
             return np.sqrt(np.sum((np.array(x) - np.array(y)) ** 2) / len(x))
 
         # list of (energy, volume) tuples
-        e_v = list(zip(self.energies, self.volumes, strict=False))
+        e_v = list(zip(self.energies, self.volumes, strict=True))
         n_data = len(e_v)
         # minimum number of data points used for fitting
         n_data_min = max(n_data - 2 * min_ndata_factor, min_poly_order + 1)

src/pymatgen/analysis/ewald.py

@@ -332,7 +332,7 @@ def _calc_recip(self):
         s_reals = np.sum(oxi_states[None, :] * np.cos(grs), 1)
         s_imags = np.sum(oxi_states[None, :] * np.sin(grs), 1)
 
-        for g, g2, gr, exp_val, s_real, s_imag in zip(gs, g2s, grs, exp_vals, s_reals, s_imags, strict=False):
+        for g, g2, gr, exp_val, s_real, s_imag in zip(gs, g2s, grs, exp_vals, s_reals, s_imags, strict=True):
             # Uses the identity sin(x)+cos(x) = 2**0.5 sin(x + pi/4)
             m = np.sin((gr[None, :] + math.pi / 4) - gr[:, None])
             m *= exp_val / g2

src/pymatgen/analysis/ferroelectricity/polarization.py

@@ -307,7 +307,7 @@ def get_same_branch_polarization_data(self, convert_to_muC_per_cm2=True, all_in_
             sites.append(new_site[0])
 
         adjust_pol = []
-        for site, struct in zip(sites, d_structs, strict=False):
+        for site, struct in zip(sites, d_structs, strict=True):
             adjust_pol.append(np.multiply(site.frac_coords, np.array(struct.lattice.lengths)).ravel())
         return np.array(adjust_pol)
 

src/pymatgen/analysis/interface_reactions.py

@@ -186,7 +186,7 @@ def get_kinks(self) -> list[tuple[int, float, float, Reaction, float]]:
 
         index_kink = range(1, len(critical_comp) + 1)
 
-        return list(zip(index_kink, x_kink, energy_kink, react_kink, energy_per_rxt_formula, strict=False))
+        return list(zip(index_kink, x_kink, energy_kink, react_kink, energy_per_rxt_formula, strict=True))
 
     def plot(self, backend: Literal["plotly", "matplotlib"] = "plotly") -> Figure | plt.Figure:
         """
@@ -326,7 +326,7 @@ def _get_elem_amt_in_rxn(self, rxn: Reaction) -> float:
 
     def _get_plotly_figure(self) -> Figure:
         """Get a Plotly figure of reaction kinks diagram."""
-        kinks = map(list, zip(*self.get_kinks(), strict=False))
+        kinks = map(list, zip(*self.get_kinks(), strict=True))
         _, x, energy, reactions, _ = kinks
 
         lines = Scatter(
@@ -348,7 +348,7 @@ def _get_plotly_figure(self) -> Figure:
 
         labels = [
             f"{htmlify(str(r))} <br>\u0394E<sub>rxn</sub> = {round(e, 3)} eV/atom"
-            for r, e in zip(reactions, energy, strict=False)
+            for r, e in zip(reactions, energy, strict=True)
         ]
 
         markers = Scatter(
@@ -392,13 +392,13 @@ def _get_matplotlib_figure(self) -> plt.Figure:
         ax = pretty_plot(8, 5)
         plt.xlim([-0.05, 1.05])  # plot boundary is 5% wider on each side
 
-        kinks = list(zip(*self.get_kinks(), strict=False))
+        kinks = list(zip(*self.get_kinks(), strict=True))
         _, x, energy, reactions, _ = kinks
 
         plt.plot(x, energy, "o-", markersize=8, c="navy", zorder=1)
         plt.scatter(self.minimum[0], self.minimum[1], marker="*", c="red", s=400, zorder=2)
 
-        for x_coord, y_coord, rxn in zip(x, energy, reactions, strict=False):
+        for x_coord, y_coord, rxn in zip(x, energy, reactions, strict=True):
             products = ", ".join(
                 [latexify(p.reduced_formula) for p in rxn.products if not np.isclose(rxn.get_coeff(p), 0)]
             )
@@ -438,7 +438,7 @@ def _get_xaxis_title(self, latex: bool = True) -> str:
     def _get_plotly_annotations(x: list[float], y: list[float], reactions: list[Reaction]):
         """Get dictionary of annotations for the Plotly figure layout."""
         annotations = []
-        for x_coord, y_coord, rxn in zip(x, y, reactions, strict=False):
+        for x_coord, y_coord, rxn in zip(x, y, reactions, strict=True):
             products = ", ".join(
                 [htmlify(p.reduced_formula) for p in rxn.products if not np.isclose(rxn.get_coeff(p), 0)]
             )

src/pymatgen/analysis/interfaces/coherent_interfaces.py

@@ -143,7 +143,7 @@ def _find_terminations(self):
         self._terminations = {
             (film_label, sub_label): (film_shift, sub_shift)
             for (film_label, film_shift), (sub_label, sub_shift) in product(
-                zip(film_terminations, film_shifts, strict=False), zip(sub_terminations, sub_shifts, strict=False)
+                zip(film_terminations, film_shifts, strict=True), zip(sub_terminations, sub_shifts, strict=True)
             )
         }
         self.terminations = list(self._terminations)

src/pymatgen/analysis/interfaces/zsl.py

@@ -172,7 +172,7 @@ def get_equiv_transformations(self, transformation_sets, film_vectors, substrate
             for (f_trans, s_trans), (f, s) in zip(
                 product(film_transformations, substrate_transformations),
                 product(films, substrates),
-                strict=False,
+                strict=True,
             ):
                 if is_same_vectors(
                     f,