Add property chemical_system_set to Composition and `SiteCollecti…

…on` (#3864)

* add property Composition.chemical_system_set -> set[str]

* Add test for Composition.elements and chemical_system_set properties

test_elements
test_chemical_system_set

* add properties chemical_system
chemical_system_set to SiteCollection

* test chemical_system and chemical_system_set in TestIStructure and TestIMolecule

* simplify count increments

* rename unreadable variables

pymatgen/analysis/chemenv/utils/coordination_geometry_utils.py

@@ -370,14 +370,19 @@ def solid_angle(center, 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]))
+    cross_products = [np.cross(r[i + 1], r[i]) for i in range(len(r) - 1)]
+    cross_products.append(np.cross(r[1], r[0]))
     phi = 0.0
-    for idx in range(len(n) - 1):
+    for idx in range(len(cross_products) - 1):
         try:
-            value = math.acos(-np.dot(n[idx], n[idx + 1]) / (np.linalg.norm(n[idx]) * np.linalg.norm(n[idx + 1])))
+            value = math.acos(
+                -np.dot(cross_products[idx], cross_products[idx + 1])
+                / (np.linalg.norm(cross_products[idx]) * np.linalg.norm(cross_products[idx + 1]))
+            )
         except ValueError:
-            cos = -np.dot(n[idx], n[idx + 1]) / (np.linalg.norm(n[idx]) * np.linalg.norm(n[idx + 1]))
+            cos = -np.dot(cross_products[idx], cross_products[idx + 1]) / (
+                np.linalg.norm(cross_products[idx]) * np.linalg.norm(cross_products[idx + 1])
+            )
             if 0.999999999999 < cos < 1.000000000001:
                 value = math.acos(1.0)
             elif -0.999999999999 > cos > -1.000000000001:

pymatgen/analysis/local_env.py

@@ -102,13 +102,13 @@ def _get_ionic_radii(self):
         vnn = VoronoiNN()
 
         def nearest_key(sorted_vals: list[int], skey: int) -> int:
-            n = bisect_left(sorted_vals, skey)
-            if n == len(sorted_vals):
+            idx = bisect_left(sorted_vals, skey)
+            if idx == len(sorted_vals):
                 return sorted_vals[-1]
-            if n == 0:
+            if idx == 0:
                 return sorted_vals[0]
-            before = sorted_vals[n - 1]
-            after = sorted_vals[n]
+            before = sorted_vals[idx - 1]
+            after = sorted_vals[idx]
             if after - skey < skey - before:
                 return after
             return before
@@ -133,7 +133,7 @@ def nearest_key(sorted_vals: list[int], skey: int) -> int:
                 oxi_state = nearest_key(tab_oxi_states, oxi_state)
                 radius = _ion_radii[el][str(oxi_state)][str(coord_no)]
             except KeyError:
-                new_coord_no = coord_no + 1 if vnn.get_cn(self._structure, idx) - coord_no > 0 else coord_no - 1
+                new_coord_no = coord_no + (1 if vnn.get_cn(self._structure, idx) - coord_no > 0 else -1)
                 try:
                     radius = _ion_radii[el][str(oxi_state)][str(new_coord_no)]
                     coord_no = new_coord_no
@@ -144,7 +144,7 @@ def nearest_key(sorted_vals: list[int], skey: int) -> int:
                     for val in tab_coords:
                         if val > coord_no:
                             break
-                        idx = idx + 1
+                        idx += 1
                     if idx == len(tab_coords):
                         key = str(tab_coords[-1])
                         radius = _ion_radii[el][str(oxi_state)][key]
@@ -2887,7 +2887,7 @@ def get_order_parameters(
                         if k > j:
                             distjk_unique.append(distjk[j][kk])
                         rjknorm[j].append(rjk[j][kk] / distjk[j][kk])
-                        kk = kk + 1
+                        kk += 1
         # Initialize OP list and, then, calculate OPs.
         ops = [0.0 for t in self._types]
         # norms = [[[] for j in range(nneigh)] for t in self._types]

pymatgen/analysis/magnetism/jahnteller.py

@@ -452,8 +452,7 @@ def _estimate_spin_state(
 
     @staticmethod
     def mu_so(species: str | Species, motif: Literal["oct", "tet"], spin_state: Literal["high", "low"]) -> float | None:
-        """Calculate the spin-only magnetic moment for a
-        given species. Only supports transition metals.
+        """Calculate the spin-only magnetic moment for a given species. Only supports transition metals.
 
         Args:
             species: Species
@@ -466,8 +465,8 @@ def mu_so(species: str | Species, motif: Literal["oct", "tet"], spin_state: Lite
         """
         try:
             sp = get_el_sp(species)
-            n = sp.get_crystal_field_spin(coordination=motif, spin_config=spin_state)
+            unpaired_spins = sp.get_crystal_field_spin(coordination=motif, spin_config=spin_state)
             # calculation spin-only magnetic moment for this number of unpaired spins
-            return np.sqrt(n * (n + 2))
+            return np.sqrt(unpaired_spins * (unpaired_spins + 2))
         except AttributeError:
             return None

pymatgen/analysis/phase_diagram.py

@@ -2421,7 +2421,7 @@ class (pymatgen.analysis.chempot_diagram).
             comp = entry.composition
             is_x = comp.get_atomic_fraction(el0) < 0.01
             is_y = comp.get_atomic_fraction(el1) < 0.01
-            n = len(coords)
+            n_coords = len(coords)
             if not (is_x and is_y):
                 if is_x:
                     coords = sorted(coords, key=lambda c: c[1])
@@ -2439,11 +2439,11 @@ class (pymatgen.analysis.chempot_diagram).
                         plt.plot(x, y, "k")
                         center_x += coords[idx][0]
                         center_y += min(ylim)
-                xy = (center_x / (n + 2), center_y / (n + 2))
+                xy = (center_x / (n_coords + 2), center_y / (n_coords + 2))
             else:
                 center_x = sum(coord[0] for coord in coords) + xlim[0]
                 center_y = sum(coord[1] for coord in coords) + ylim[0]
-                xy = (center_x / (n + 1), center_y / (n + 1))
+                xy = (center_x / (n_coords + 1), center_y / (n_coords + 1))
 
             ax.annotate(
                 latexify(entry.name),

pymatgen/analysis/piezo_sensitivity.py

@@ -499,7 +499,7 @@ def get_stable_FCM(self, fcm, fcmasum=10):
                 if eigs[eig_sort[idx]] > 1e-6:
                     unstable_modes = 1
             if unstable_modes == 1:
-                count = count + 1
+                count += 1
                 continue
             check = 1
 

pymatgen/analysis/structure_analyzer.py

@@ -347,11 +347,13 @@ def solid_angle(center, coords):
     origin = np.array(center)
     radii = [np.array(c) - origin for c in coords]
     radii.append(radii[0])
-    n = [np.cross(radii[i + 1], radii[i]) for i in range(len(radii) - 1)]
-    n.append(np.cross(radii[1], radii[0]))
+    cross_products = [np.cross(radii[i + 1], radii[i]) for i in range(len(radii) - 1)]
+    cross_products.append(np.cross(radii[1], radii[0]))
     vals = []
-    for i in range(len(n) - 1):
-        v = -np.dot(n[i], n[i + 1]) / (np.linalg.norm(n[i]) * np.linalg.norm(n[i + 1]))
+    for i in range(len(cross_products) - 1):
+        v = -np.dot(cross_products[i], cross_products[i + 1]) / (
+            np.linalg.norm(cross_products[i]) * np.linalg.norm(cross_products[i + 1])
+        )
         vals.append(acos(np.clip(v, -1, 1)))
     phi = sum(vals)
     return phi + (3 - len(radii)) * pi

pymatgen/analysis/surface_analysis.py

@@ -156,10 +156,12 @@ def gibbs_binding_energy(self, eads=False):
                 (True) or the binding energy (False) which is just
                 adsorption energy normalized by number of adsorbates.
         """
-        n = self.get_unit_primitive_area
+        surface_area = self.get_unit_primitive_area
         n_ads = self.Nads_in_slab
 
-        BE = (self.energy - n * self.clean_entry.energy) / n_ads - sum(ads.energy_per_atom for ads in self.adsorbates)
+        BE = (self.energy - surface_area * self.clean_entry.energy) / n_ads - sum(
+            ads.energy_per_atom for ads in self.adsorbates
+        )
         return BE * n_ads if eads else BE
 
     def surface_energy(self, ucell_entry, ref_entries=None):

pymatgen/core/composition.py

@@ -462,6 +462,20 @@ def elements(self) -> list[Element | Species | DummySpecies]:
         """List of elements in Composition."""
         return list(self)
 
+    @property
+    def chemical_system_set(self) -> set[str]:
+        """The set of elements in the Composition. E.g. {"O", "Si"} for SiO2."""
+        return {el.symbol for el in self.elements}
+
+    @property
+    def chemical_system(self) -> str:
+        """The chemical system of a Composition, for example "O-Si" for
+        SiO2. Chemical system is a string of a list of elements
+        sorted alphabetically and joined by dashes, by convention for use
+        in database keys.
+        """
+        return "-".join(sorted(el.symbol for el in self.elements))
+
     def __str__(self) -> str:
         return " ".join(f"{key}{formula_double_format(val, ignore_ones=False)}" for key, val in self.as_dict().items())
 
@@ -599,15 +613,6 @@ def anonymized_formula(self) -> str:
             anon += f"{elem}{amt_str}"
         return anon
 
-    @property
-    def chemical_system(self) -> str:
-        """The chemical system of a Composition, for example "O-Si" for
-        SiO2. Chemical system is a string of a list of elements
-        sorted alphabetically and joined by dashes, by convention for use
-        in database keys.
-        """
-        return "-".join(sorted(el.symbol for el in self.elements))
-
     @property
     def valid(self) -> bool:
         """True if Composition contains valid elements or species and

pymatgen/core/periodic_table.py

@@ -1231,7 +1231,7 @@ def get_crystal_field_spin(
             spin_config ("low" | "high"): Whether the species is in a high or low spin state
 
         Returns:
-            Crystal field spin in Bohr magneton.
+            float: Crystal field spin in Bohr magneton.
 
         Raises:
             AttributeError if species is not a valid transition metal or has

pymatgen/core/structure.py

@@ -420,6 +420,18 @@ def composition(self) -> Composition:
                 elem_map[species] += occu
         return Composition(elem_map)
 
+    @property
+    def chemical_system(self) -> str:
+        """The chemical system of the structure."""
+        return self.composition.chemical_system
+
+    @property
+    def chemical_system_set(self) -> set[str]:
+        """The set of chemical systems in the structure. E.g. {"Al", "Ga", "In", "N"} for
+        a AlGaInN quaternary.
+        """
+        return self.composition.chemical_system_set
+
     @property
     def charge(self) -> float:
         """The net charge of the structure based on oxidation states. If

pymatgen/io/exciting/inputs.py

@@ -249,7 +249,7 @@ def write_etree(self, celltype, cartesian=False, bandstr=False, symprec: float =
                     coord = f"{coord2[0]:16.8f} {coord2[1]:16.8f} {coord2[2]:16.8f}"
 
                 # write atomic positions
-                index = index + 1
+                index += 1
                 _ = ET.SubElement(species, "atom", coord=coord)
 
         # write bandstructure if needed

pymatgen/io/lobster/inputs.py

@@ -240,7 +240,7 @@ def _get_nbands(self, structure: Structure):
         no_basis_functions = 0
         for basis in basis_functions:
             if "s" in basis:
-                no_basis_functions = no_basis_functions + 1
+                no_basis_functions += 1
             elif "p" in basis:
                 no_basis_functions = no_basis_functions + 3
             elif "d" in basis:

pymatgen/io/lobster/lobsterenv.py

@@ -1312,25 +1312,24 @@ def from_Lobster(
                 ce_dict = None
 
             if list_neighisite[isite] is not None:
-                for ineighsite, neighsite in enumerate(list_neighsite[isite]):
-                    diff = neighsite.frac_coords - structure[list_neighisite[isite][ineighsite]].frac_coords
-                    rounddiff = np.round(diff)
-                    if not np.allclose(diff, rounddiff):
+                for idx_neigh_site, neigh_site in enumerate(list_neighsite[isite]):
+                    diff = neigh_site.frac_coords - structure[list_neighisite[isite][idx_neigh_site]].frac_coords
+                    round_diff = np.round(diff)
+                    if not np.allclose(diff, round_diff):
                         raise ValueError(
                             "Weird, differences between one site in a periodic image cell is not integer ..."
                         )
-                    nb_image_cell = np.array(rounddiff, int)
+                    nb_image_cell = np.array(round_diff, int)
 
                     all_nbs_sites_indices_here.append(counter)
 
-                    all_nbs_sites.append(
-                        {
-                            "site": neighsite,
-                            "index": list_neighisite[isite][ineighsite],
-                            "image_cell": nb_image_cell,
-                        }
-                    )
-                    counter = counter + 1
+                    neighbor = {
+                        "site": neigh_site,
+                        "index": list_neighisite[isite][idx_neigh_site],
+                        "image_cell": nb_image_cell,
+                    }
+                    all_nbs_sites.append(neighbor)
+                    counter += 1
 
                 all_nbs_sites_indices.append(all_nbs_sites_indices_here)
             else:

pymatgen/io/lobster/outputs.py

@@ -154,7 +154,7 @@ def __init__(
                 cohp = {spin: data[2 * (bond + s * (num_bonds + 1)) + 3] for s, spin in enumerate(spins)}
                 icohp = {spin: data[2 * (bond + s * (num_bonds + 1)) + 4] for s, spin in enumerate(spins)}
                 if orbs is None:
-                    bond_num = bond_num + 1
+                    bond_num += 1
                     label = str(bond_num)
                     cohp_data[label] = {
                         "COHP": cohp,
@@ -207,7 +207,7 @@ def __init__(
 
                 icohp = {spin: data[2 * (bond + s * (num_bonds)) + 2] for s, spin in enumerate(spins)}
                 if orbs is None:
-                    bond_num = bond_num + 1
+                    bond_num += 1
                     label = str(bond_num)
                     cohp_data[label] = {
                         "COHP": cohp,
@@ -704,11 +704,10 @@ def _parse_doscar(self):
                 pdos = defaultdict(dict)
                 data = dos[atom + 1]
                 _, ncol = data.shape
-                orbnumber = 0
-                for j in range(1, ncol):
-                    orb = orbitals[atom + 1][orbnumber]
+
+                for orb_num, j in enumerate(range(1, ncol)):
+                    orb = orbitals[atom + 1][orb_num]
                     pdos[orb][spin] = data[:, j]
-                    orbnumber = orbnumber + 1
                 pdoss += [pdos]
         else:
             tdensities[Spin.up] = doshere[:, 1]
@@ -720,13 +719,13 @@ def _parse_doscar(self):
                 pdos = defaultdict(dict)
                 data = dos[atom + 1]
                 _, ncol = data.shape
-                orbnumber = 0
+                orb_num = 0
                 for j in range(1, ncol):
                     spin = Spin.down if j % 2 == 0 else Spin.up
-                    orb = orbitals[atom + 1][orbnumber]
+                    orb = orbitals[atom + 1][orb_num]
                     pdos[orb][spin] = data[:, j]
                     if j % 2 == 0:
-                        orbnumber = orbnumber + 1
+                        orb_num += 1
                 pdoss += [pdos]
 
         self._efermi = efermi

pymatgen/io/xtb/outputs.py

@@ -55,8 +55,8 @@ def _parse_crest_output(self):
 
         # Get CREST command
         crest_cmd = None
-        with open(output_filepath, encoding="utf-8") as xtbout_file:
-            for line in xtbout_file:
+        with open(output_filepath, encoding="utf-8") as xtb_out_file:
+            for line in xtb_out_file:
                 if "> crest" in line:
                     crest_cmd = line.strip()[8:]
                     break
@@ -88,11 +88,11 @@ def _parse_crest_output(self):
             chg = int(str_chg) if "-" in str_chg else int(str_chg[-1])
 
         # Check for proper termination
-        with open(output_filepath, "rb+") as xtbout_file:
-            xtbout_file.seek(-2, 2)
-            while xtbout_file.read(1) != b"\n":
-                xtbout_file.seek(-2, 1)
-            end_bstring = xtbout_file.read()
+        with open(output_filepath, "rb+") as xtb_out_file:
+            xtb_out_file.seek(-2, 2)
+            while xtb_out_file.read(1) != b"\n":
+                xtb_out_file.seek(-2, 1)
+            end_bstring = xtb_out_file.read()
             if b"CREST terminated normally." in end_bstring:
                 self.properly_terminated = True
 
@@ -111,8 +111,8 @@ def _parse_crest_output(self):
             )
             conformer_degeneracies = []
             energies = []
-            with open(output_filepath, encoding="utf-8") as xtbout_file:
-                for line in xtbout_file:
+            with open(output_filepath, encoding="utf-8") as xtb_out_file:
+                for line in xtb_out_file:
                     conformer_match = conformer_pattern.match(line)
                     rotamer_match = rotamer_pattern.match(line)
                     if conformer_match:
@@ -148,19 +148,19 @@ def _parse_crest_output(self):
                 print(f"{final_rotamer_filename} not found, no rotamer list processed")
 
             # Get lowest energy conformer from 'crest_best.xyz'
-            crestbest_path = os.path.join(self.path, "crest_best.xyz")
+            crest_best_path = os.path.join(self.path, "crest_best.xyz")
             try:
-                lowest_e_struct = Molecule.from_file(crestbest_path)
+                lowest_e_struct = Molecule.from_file(crest_best_path)
                 lowest_e_struct.set_charge_and_spin(charge=chg)
                 self.lowest_energy_structure = lowest_e_struct
             except FileNotFoundError:
-                print(f"{crestbest_path} not found")
+                print(f"{crest_best_path} not found")
 
         else:
-            crestbest_path = os.path.join(self.path, "crest_best.xyz")
+            crest_best_path = os.path.join(self.path, "crest_best.xyz")
             try:
-                lowest_e_struct = Molecule.from_file(crestbest_path)
+                lowest_e_struct = Molecule.from_file(crest_best_path)
                 lowest_e_struct.set_charge_and_spin(charge=chg)
                 self.lowest_energy_structure = lowest_e_struct
             except FileNotFoundError:
-                print(f"{crestbest_path} not found")
+                print(f"{crest_best_path} not found")