Test EntrySet.ground_states and CIF writing in NEBSet.write_input (…

…#3732)

* use PbcLike typing in coord.py

* rename variables for readability

* fix structure matching and reduction algo doc strings

* breaking: fix typo in method name has_antiband_states_below_efermi

* refactor ground_states and uniquelines using set comprehension

* add test_ground_states

* check CIF parsable and correct len in TestMITNEBSet.test_write_input

* doc str fixes

pymatgen/analysis/adsorption.py

@@ -672,22 +672,22 @@ def plot_slab(
     alphas = alphas.clip(min=0)
     corner = [0, 0, slab.lattice.get_fractional_coords(coords[-1])[-1]]
     corner = slab.lattice.get_cartesian_coords(corner)[:2]
-    verts = orig_cell[:2, :2]
-    lattice_sum = verts[0] + verts[1]
+    vertices = orig_cell[:2, :2]
+    lattice_sum = vertices[0] + vertices[1]
     # inverse coords, sites, alphas, to show other side of slab
     if inverse:
         alphas = np.array(reversed(alphas))
         sites = list(reversed(sites))
         coords = np.array(reversed(coords))
     # Draw circles at sites and stack them accordingly
     for n, coord in enumerate(coords):
-        r = sites[n].species.elements[0].atomic_radius * scale
-        ax.add_patch(patches.Circle(coord[:2] - lattice_sum * (repeat // 2), r, color="w", zorder=2 * n))
+        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]
         ax.add_patch(
             patches.Circle(
                 coord[:2] - lattice_sum * (repeat // 2),
-                r,
+                radius,
                 facecolor=color,
                 alpha=alphas[n],
                 edgecolor="k",
@@ -708,12 +708,12 @@ def plot_slab(
         ax.plot(*zip(*ads_sites), color="k", marker="x", markersize=10, mew=1, linestyle="", zorder=10000)
     # Draw unit cell
     if draw_unit_cell:
-        verts = np.insert(verts, 1, lattice_sum, axis=0).tolist()
-        verts += [[0.0, 0.0]]
-        verts = [[0.0, 0.0], *verts]
+        vertices = np.insert(vertices, 1, lattice_sum, axis=0).tolist()
+        vertices += [[0.0, 0.0]]
+        vertices = [[0.0, 0.0], *vertices]
         codes = [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY]
-        verts = [(np.array(vert) + corner).tolist() for vert in verts]
-        path = Path(verts, codes)
+        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)
         ax.add_patch(patch)
     ax.set_aspect("equal")

pymatgen/analysis/bond_valence.py

@@ -36,8 +36,7 @@
 
 
 def calculate_bv_sum(site, nn_list, scale_factor=1.0):
-    """
-    Calculates the BV sum of a site.
+    """Calculates the BV sum of a site.
 
     Args:
         site (PeriodicSite): The central site to calculate the bond valence
@@ -63,8 +62,7 @@ def calculate_bv_sum(site, nn_list, scale_factor=1.0):
 
 
 def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):
-    """
-    Calculates the BV sum of a site for unordered structures.
+    """Calculates the BV sum of a site for unordered structures.
 
     Args:
         site (PeriodicSite): The central site to calculate the bond valence
@@ -82,7 +80,7 @@ def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):
     # site "site" is obtained as :
     # \sum_{nn} \sum_j^N \sum_k^{N_{nn}} f_{site}_j f_{nn_i}_k vij_full
     # where vij_full is the valence bond of the fully occupied bond
-    bvsum = 0
+    bv_sum = 0
     for specie1, occu1 in site.species.items():
         el1 = Element(specie1.symbol)
         for nn in nn_list:
@@ -95,8 +93,8 @@ def calculate_bv_sum_unordered(site, nn_list, scale_factor=1):
                     c2 = BV_PARAMS[el2]["c"]
                     R = r1 + r2 - r1 * r2 * (sqrt(c1) - sqrt(c2)) ** 2 / (c1 * r1 + c2 * r2)
                     vij = exp((R - nn.nn_distance * scale_factor) / 0.31)
-                    bvsum += occu1 * occu2 * vij * (1 if el1.X < el2.X else -1)
-    return bvsum
+                    bv_sum += occu1 * occu2 * vij * (1 if el1.X < el2.X else -1)
+    return bv_sum
 
 
 class BVAnalyzer:

pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py

@@ -342,13 +342,13 @@ def rectangle_surface_intersection(
     xmax = min(x2, bounds_lower[1])
 
     def diff(x):
-        flwx = f_lower(x)
-        fupx = f_upper(x)
-        minup = np.min([fupx, y2 * np.ones_like(fupx)], axis=0)
-        maxlw = np.max([flwx, y1 * np.ones_like(flwx)], axis=0)
-        zeros = np.zeros_like(fupx)
-        upper = np.where(y2 >= flwx, np.where(y1 <= fupx, minup, zeros), zeros)
-        lower = np.where(y1 <= fupx, np.where(y2 >= flwx, maxlw, zeros), zeros)
+        f_low_x = f_lower(x)
+        f_up_x = f_upper(x)
+        min_up = np.min([f_up_x, y2 * np.ones_like(f_up_x)], axis=0)
+        max_lw = np.max([f_low_x, y1 * np.ones_like(f_low_x)], axis=0)
+        zeros = np.zeros_like(f_up_x)
+        upper = np.where(y2 >= f_low_x, np.where(y1 <= f_up_x, min_up, zeros), zeros)
+        lower = np.where(y1 <= f_up_x, np.where(y2 >= f_low_x, max_lw, zeros), zeros)
         return upper - lower
 
     return quad(diff, xmin, xmax)
@@ -359,16 +359,14 @@ def solid_angle(center, coords):
     Helper method to calculate the solid angle of a set of coords from the center.
 
     Args:
-        center:
-            Center to measure solid angle from.
-        coords:
-            List of coords to determine solid angle.
+        center: Center to measure solid angle from.
+        coords: List of coords to determine solid angle.
 
     Returns:
         The solid angle.
     """
-    o = np.array(center)
-    r = [np.array(c) - o for c in coords]
+    origin = np.array(center)
+    r = [np.array(c) - origin for c in coords]
     r.append(r[0])
     n = [np.cross(r[i + 1], r[i]) for i in range(len(r) - 1)]
     n.append(np.cross(r[1], r[0]))

pymatgen/analysis/ferroelectricity/polarization.py

@@ -125,8 +125,8 @@ def get_nearest_site(struct: Structure, coords: Sequence[float], site: PeriodicS
         Closest site and distance.
     """
     index = struct.index(site)
-    r = r or np.linalg.norm(np.sum(struct.lattice.matrix, axis=0))
-    ns = struct.get_sites_in_sphere(coords, r, include_index=True)
+    radius = r or np.linalg.norm(np.sum(struct.lattice.matrix, axis=0))
+    ns = struct.get_sites_in_sphere(coords, radius, include_index=True)
     # Get sites with identical index to site
     ns = [n for n in ns if n[2] == index]
     # Sort by distance to coords

pymatgen/analysis/local_env.py

@@ -4088,19 +4088,21 @@ def _semicircle_integral(dist_bins, idx):
         Returns:
             float: integral of portion of unit semicircle
         """
-        r = 1
+        radius = 1
 
         x1 = dist_bins[idx]
         x2 = dist_bins[idx + 1]
 
         if dist_bins[idx] == 1:
-            area1 = 0.25 * math.pi * r**2
+            area1 = 0.25 * math.pi * radius**2
         else:
-            area1 = 0.5 * ((x1 * math.sqrt(r**2 - x1**2)) + (r**2 * math.atan(x1 / math.sqrt(r**2 - x1**2))))
+            area1 = 0.5 * (
+                (x1 * math.sqrt(radius**2 - x1**2)) + (radius**2 * math.atan(x1 / math.sqrt(radius**2 - x1**2)))
+            )
 
-        area2 = 0.5 * ((x2 * math.sqrt(r**2 - x2**2)) + (r**2 * math.atan(x2 / math.sqrt(r**2 - x2**2))))
+        area2 = 0.5 * ((x2 * math.sqrt(radius**2 - x2**2)) + (radius**2 * math.atan(x2 / math.sqrt(radius**2 - x2**2))))
 
-        return (area1 - area2) / (0.25 * math.pi * r**2)
+        return (area1 - area2) / (0.25 * math.pi * radius**2)
 
     @staticmethod
     def transform_to_length(nn_data, length):

pymatgen/analysis/phase_diagram.py

@@ -618,10 +618,10 @@ def _get_facet_and_simplex(self, comp: Composition) -> tuple[Simplex, Simplex]:
         Args:
             comp (Composition): A composition
         """
-        c = self.pd_coords(comp)
-        for f, s in zip(self.facets, self.simplexes):
-            if s.in_simplex(c, PhaseDiagram.numerical_tol / 10):
-                return f, s
+        coord = self.pd_coords(comp)
+        for facet, simplex in zip(self.facets, self.simplexes):
+            if simplex.in_simplex(coord, PhaseDiagram.numerical_tol / 10):
+                return facet, simplex
 
         raise RuntimeError(f"No facet found for {comp = }")
 
@@ -632,10 +632,12 @@ def _get_all_facets_and_simplexes(self, comp):
         Args:
             comp (Composition): A composition
         """
-        c = self.pd_coords(comp)
+        coords = self.pd_coords(comp)
 
         all_facets = [
-            f for f, s in zip(self.facets, self.simplexes) if s.in_simplex(c, PhaseDiagram.numerical_tol / 10)
+            facet
+            for facet, simplex in zip(self.facets, self.simplexes)
+            if simplex.in_simplex(coords, PhaseDiagram.numerical_tol / 10)
         ]
 
         if not all_facets:
@@ -1973,8 +1975,8 @@ def fmt(fl):
 
                         for c, entry in zip(coeffs[:-1], face_entries):
                             if c > tol:
-                                r = entry.composition.reduced_composition
-                                products.append(f"{fmt(c / r.num_atoms * factor)} {r.reduced_formula}")
+                                redu_comp = entry.composition.reduced_composition
+                                products.append(f"{fmt(c / redu_comp.num_atoms * factor)} {redu_comp.reduced_formula}")
                                 product_entries.append((c, entry))
                                 energy += c * entry.energy_per_atom
 
@@ -3756,11 +3758,7 @@ def uniquelines(q):
         setoflines:
             A set of tuple of lines. E.g., ((1,2), (1,3), (2,3), ....)
     """
-    setoflines = set()
-    for facets in q:
-        for line in itertools.combinations(facets, 2):
-            setoflines.add(tuple(sorted(line)))
-    return setoflines
+    return {tuple(sorted(line)) for facets in q for line in itertools.combinations(facets, 2)}
 
 
 def triangular_coord(coord):

pymatgen/analysis/structure_matcher.py

@@ -969,9 +969,10 @@ def get_rms_anonymous(self, struct1, struct2):
             struct2 (Structure): 2nd structure
 
         Returns:
-            tuple: [min_rms, min_mapping]: min_rms is the minimum rms distance, and min_mapping is the
-                corresponding minimal species mapping that would map
-                struct1 to struct2. (None, None) is returned if the minimax_rms exceeds the threshold.
+            tuple[float, float] | tuple[None, None]: 1st element is min_rms, 2nd is min_mapping.
+                min_rms is the minimum RMS distance, and min_mapping is the corresponding
+                minimal species mapping that would map struct1 to struct2. (None, None) is
+                returned if the minimax_rms exceeds the threshold.
         """
         struct1, struct2 = self._process_species([struct1, struct2])
         struct1, struct2, fu, s1_supercell = self._preprocess(struct1, struct2)

pymatgen/analysis/surface_analysis.py

@@ -1655,13 +1655,13 @@ def solve_equilibrium_point(self, analyzer1, analyzer2, delu_dict=None, delu_def
         # Now calculate r
         delta_gamma = wulff1.weighted_surface_energy - wulff2.weighted_surface_energy
         delta_E = self.bulk_gform(analyzer1.ucell_entry) - self.bulk_gform(analyzer2.ucell_entry)
-        r = (-3 * delta_gamma) / (delta_E)
+        radius = (-3 * delta_gamma) / (delta_E)
 
-        return r / 10 if units == "nanometers" else r
+        return radius / 10 if units == "nanometers" else radius
 
     def wulff_gform_and_r(
         self,
-        wulffshape,
+        wulff_shape,
         bulk_entry,
         r,
         from_sphere_area=False,
@@ -1674,7 +1674,7 @@ def wulff_gform_and_r(
         Calculates the formation energy of the particle with arbitrary radius r.
 
         Args:
-            wulffshape (WulffShape): Initial, unscaled WulffShape
+            wulff_shape (WulffShape): Initial unscaled WulffShape
             bulk_entry (ComputedStructureEntry): Entry of the corresponding bulk.
             r (float (Ang)): Arbitrary effective radius of the WulffShape
             from_sphere_area (bool): There are two ways to calculate the bulk
@@ -1690,8 +1690,8 @@ def wulff_gform_and_r(
             particle formation energy (float in keV), effective radius
         """
         # Set up
-        miller_se_dict = wulffshape.miller_energy_dict
-        new_wulff = self.scaled_wulff(wulffshape, r)
+        miller_se_dict = wulff_shape.miller_energy_dict
+        new_wulff = self.scaled_wulff(wulff_shape, r)
         new_wulff_area = new_wulff.miller_area_dict
 
         # calculate surface energy of the particle
@@ -1708,7 +1708,7 @@ def wulff_gform_and_r(
             # By approximating the particle as a perfect sphere
             w_vol = (4 / 3) * np.pi * r**3
             sphere_sa = 4 * np.pi * r**2
-            tot_wulff_se = wulffshape.weighted_surface_energy * sphere_sa
+            tot_wulff_se = wulff_shape.weighted_surface_energy * sphere_sa
             Ebulk = self.bulk_gform(bulk_entry) * w_vol
             new_r = r
 
@@ -1735,7 +1735,7 @@ def bulk_gform(bulk_entry):
         """
         return bulk_entry.energy / bulk_entry.structure.volume
 
-    def scaled_wulff(self, wulffshape, r):
+    def scaled_wulff(self, wulff_shape, r):
         """
         Scales the Wulff shape with an effective radius r. Note that the resulting
             Wulff does not necessarily have the same effective radius as the one
@@ -1744,22 +1744,22 @@ def scaled_wulff(self, wulffshape, r):
             multiplied by the given effective radius.
 
         Args:
-            wulffshape (WulffShape): Initial, unscaled WulffShape
+            wulff_shape (WulffShape): Initial, unscaled WulffShape
             r (float): Arbitrary effective radius of the WulffShape
 
         Returns:
             WulffShape (scaled by r)
         """
         # get the scaling ratio for the energies
-        r_ratio = r / wulffshape.effective_radius
-        miller_list = list(wulffshape.miller_energy_dict)
+        r_ratio = r / wulff_shape.effective_radius
+        miller_list = list(wulff_shape.miller_energy_dict)
         # Normalize the magnitude of the facet normal vectors
         # of the Wulff shape by the minimum surface energy.
-        se_list = np.array(list(wulffshape.miller_energy_dict.values()))
+        se_list = np.array(list(wulff_shape.miller_energy_dict.values()))
         # Scale the magnitudes by r_ratio
         scaled_se = se_list * r_ratio
 
-        return WulffShape(wulffshape.lattice, miller_list, scaled_se, symprec=self.symprec)
+        return WulffShape(wulff_shape.lattice, miller_list, scaled_se, symprec=self.symprec)
 
     def plot_one_stability_map(
         self,
@@ -1800,22 +1800,22 @@ def plot_one_stability_map(
         """
         plt = plt or pretty_plot(width=8, height=7)
 
-        wulffshape = analyzer.wulff_from_chempot(delu_dict=delu_dict, delu_default=delu_default, symprec=self.symprec)
+        wulff_shape = analyzer.wulff_from_chempot(delu_dict=delu_dict, delu_default=delu_default, symprec=self.symprec)
 
         gform_list, r_list = [], []
-        for r in np.linspace(1e-6, max_r, increments):
-            gform, r = self.wulff_gform_and_r(
-                wulffshape,
+        for radius in np.linspace(1e-6, max_r, increments):
+            gform, radius = self.wulff_gform_and_r(
+                wulff_shape,
                 analyzer.ucell_entry,
-                r,
+                radius,
                 from_sphere_area=from_sphere_area,
                 r_units=r_units,
                 e_units=e_units,
                 normalize=normalize,
                 scale_per_atom=scale_per_atom,
             )
             gform_list.append(gform)
-            r_list.append(r)
+            r_list.append(radius)
 
         ru = "nm" if r_units == "nanometers" else r"\AA"
         plt.xlabel(rf"Particle radius (${ru}$)")

pymatgen/core/lattice.py

@@ -1332,14 +1332,13 @@ def get_points_in_sphere(
             return self.get_points_in_sphere_py(frac_points=frac_points, center=center, r=r, zip_results=zip_results)
         else:
             frac_points = np.ascontiguousarray(frac_points, dtype=float)
-            lattice_matrix = np.ascontiguousarray(self.matrix, dtype=float)
+            latt_matrix = np.ascontiguousarray(self.matrix, dtype=float)
             cart_coords = np.ascontiguousarray(self.get_cartesian_coords(frac_points), dtype=float)
             pbc = np.ascontiguousarray(self.pbc, dtype=int)
-            r = float(r)
             center_coords = np.ascontiguousarray([center], dtype=float)
 
             _, indices, images, distances = find_points_in_spheres(
-                all_coords=cart_coords, center_coords=center_coords, r=r, pbc=pbc, lattice=lattice_matrix, tol=1e-8
+                all_coords=cart_coords, center_coords=center_coords, r=float(r), pbc=pbc, lattice=latt_matrix, tol=1e-8
             )
             if len(indices) < 1:
                 return [] if zip_results else [()] * 4