Extend lobsterenv for coop/cobi (#3050)

* add method to adapt extremum values

Author: naik-aakash

pymatgen/io/lobster/lobsterenv.py

@@ -58,11 +58,13 @@ def __init__(
         structure: Structure,
         filename_ICOHP: str = "ICOHPLIST.lobster",
         are_coops: bool = False,
+        are_cobis: bool = False,
         valences: list[int | float] | None = None,
         limits: tuple[float, float] | None = None,
         additional_condition: int = 0,
         only_bonds_to: list[str] | None = None,
         perc_strength_ICOHP: float = 0.15,
+        noise_cutoff: float = 0.1,
         valences_from_charges: bool = False,
         filename_CHARGE: str | None = None,
         which_charge: str = "Mulliken",
@@ -76,12 +78,13 @@ def __init__(
         """
 
         Args:
-            filename_ICOHP: (str) Path to ICOHPLIST.lobster
-            structure: (Structure): typically constructed by Structure.from_file("POSCAR")
+            filename_ICOHP: (str) Path to ICOHPLIST.lobster or ICOOPLIST.lobster or ICOBILIST.lobster
+            structure: (Structure) typically constructed by Structure.from_file("POSCAR")
             are_coops: (bool) if True, the file is a ICOOPLIST.lobster and not a ICOHPLIST.lobster; only tested for
                 ICOHPLIST.lobster so far
+            are_cobis: (bool) if True, the file is a ICOBILIST.lobster and not a ICOHPLIST.lobster
             valences: (list[int | float]): gives valence/charge for each element
-            limits (tuple[float, float] | None): limit to decide which ICOHPs should be considered
+            limits (tuple[float, float] | None): limit to decide which ICOHPs (ICOOP or ICOBI) should be considered
             additional_condition (int): Additional condition that decides which kind of bonds will be considered
                 NO_ADDITIONAL_CONDITION = 0
                 ONLY_ANION_CATION_BONDS = 1
@@ -93,7 +96,8 @@ def __init__(
             only_bonds_to: (list[str]) will only consider bonds to certain elements (e.g. ["O"] for oxygen)
             perc_strength_ICOHP: if no limits are given, this will decide which icohps will still be considered (
             relative to
-            the strongest ICOHP)
+            the strongest ICOHP (ICOOP or ICOBI)
+            noise_cutoff: if provided hardcodes the lower limit of icohps considered
             valences_from_charges: if True and path to CHARGE.lobster is provided, will use Lobster charges (
             Mulliken) instead of valences
             filename_CHARGE: (str) Path to Charge.lobster
@@ -108,16 +112,18 @@ def __init__(
             id_blist_sg2: (str) Identity of data in filename_blist_sg2,
                 e.g., "icoop" or "icobi".
         """
-        self.ICOHP = Icohplist(are_coops=are_coops, filename=filename_ICOHP)
+        self.ICOHP = Icohplist(are_coops=are_coops, are_cobis=are_cobis, filename=filename_ICOHP)
         self.Icohpcollection = self.ICOHP.icohpcollection
         self.structure = structure
         self.limits = limits
         self.only_bonds_to = only_bonds_to
         self.adapt_extremum_to_add_cond = adapt_extremum_to_add_cond
         self.are_coops = are_coops
+        self.are_cobis = are_cobis
         self.add_additional_data_sg = add_additional_data_sg
         self.filename_blist_sg1 = filename_blist_sg1
         self.filename_blist_sg2 = filename_blist_sg2
+        self.noise_cutoff = noise_cutoff
 
         allowed_arguments = ["icoop", "icobi"]
         if id_blist_sg1.lower() not in allowed_arguments or id_blist_sg2.lower() not in allowed_arguments:
@@ -152,9 +158,6 @@ def __init__(
                 are_cobis=are_cobis_id2,
             )
 
-        if are_coops:
-            raise ValueError("Algorithm only works correctly for ICOHPLIST.lobster")
-
         # will check if the additional condition is correctly delivered
         if additional_condition not in range(7):
             raise ValueError(f"Unexpected {additional_condition=}, must be one of {list(range(7))}")
@@ -364,12 +367,12 @@ def get_info_icohps_to_neighbors(self, isites=None, onlycation_isites=True):
         """
         This method returns information on the icohps of neighbors for certain sites as identified by their site id.
         This is useful for plotting the relevant cohps of a site in the structure.
+        (could be ICOOPLIST.lobster or ICOHPLIST.lobster or ICOBILIST.lobster)
 
         Args:
             isites: list of site ids. If isite==None, all isites will be used to add the icohps of the neighbors
             onlycation_isites: if True and if isite==None, it will only analyse the sites of the cations
 
-
         Returns: ICOHPNeighborsInfo
         """
         if self.valences is None and onlycation_isites:
@@ -415,7 +418,8 @@ def plot_cohps_of_neighbors(
         integrated=False,
     ):
         """
-        Will plot summed cohps (please be careful in the spin polarized case (plots might overlap (exactly!)).
+        Will plot summed cohps or cobis or coops
+        (please be careful in the spin polarized case (plots might overlap (exactly!)).
 
         Args:
             isites: list of site ids, if isite==[], all isites will be used to add the icohps of the neighbors
@@ -429,12 +433,12 @@ def plot_cohps_of_neighbors(
             integrated: bool, if true will show integrated cohp instead of cohp
 
         Returns:
-            plt of the cohps
+            plt of the cohps or coops or cobis
         """
         # include COHPPlotter and plot a sum of these COHPs
         # might include option to add Spin channels
         # implement only_bonds_to
-        cp = CohpPlotter()
+        cp = CohpPlotter(are_cobis=self.are_cobis, are_coops=self.are_coops)
 
         plotlabel, summed_cohp = self.get_info_cohps_to_neighbors(
             path_to_COHPCAR,
@@ -465,19 +469,19 @@ def get_info_cohps_to_neighbors(
         summed_spin_channels=False,
     ):
         """
-        Return info about the cohps as a summed cohp object and a label
+        Return info about the cohps (coops or cobis) as a summed cohp object and a label
          from all sites mentioned in isites with neighbors.
 
         Args:
-            path_to_COHPCAR: str, path to COHPCAR
+            path_to_COHPCAR: str, path to COHPCAR or COOPCAR or COBICAR
             isites: list of int that indicate the number of the site
             only_bonds_to: list of str, e.g. ["O"] to only show cohps of anything to oxygen
             onlycation_isites: if isites=None, only cation sites will be returned
             per_bond: will normalize per bond
             summed_spin_channels: will sum all spin channels
 
-        Returns: label for cohp (str), CompleteCohp object which describes all cohps of the sites as given by isites
-        and the other parameters
+        Returns: label for cohp (str), CompleteCohp object which describes all cohps (coops or cobis) of the sites
+        as given by isites and the other parameters
         """
         # TODO: add options for orbital-resolved cohps
         (
@@ -497,7 +501,13 @@ def get_info_cohps_to_neighbors(
             self.structure.to(filename=path, fmt="POSCAR")
 
             if not hasattr(self, "completecohp"):
-                self.completecohp = CompleteCohp.from_file(fmt="LOBSTER", filename=path_to_COHPCAR, structure_file=path)
+                self.completecohp = CompleteCohp.from_file(
+                    fmt="LOBSTER",
+                    filename=path_to_COHPCAR,
+                    structure_file=path,
+                    are_coops=self.are_coops,
+                    are_cobis=self.are_cobis,
+                )
 
         # will check that the number of bonds in ICOHPLIST and COHPCAR are identical
         # further checks could be implemented
@@ -1117,6 +1127,19 @@ def _determine_unit_cell(site):
 
         return unitcell
 
+    def _adapt_extremum_to_add_cond(self, list_icohps, percentage):
+        """
+        Convinicence method for returning the extremum of the given icohps or icoops or icobis list
+
+        Args:
+            list_icohps: can be a list of icohps or icobis or icobis
+
+        Returns: min value of input list of icohps / max value of input list of icobis or icobis
+        """
+
+        which_extr = min if not self.are_coops and not self.are_cobis else max
+        return which_extr(list_icohps) * percentage
+
     def _get_limit_from_extremum(
         self,
         icohpcollection,
@@ -1130,12 +1153,14 @@ def _get_limit_from_extremum(
 
         Args:
             icohpcollection: icohpcollection object
-            percentage: will determine which ICOHPs will be considered (only 0.15 from the maximum value)
+            percentage: will determine which ICOHPs or ICOOP or ICOBI will be considered
+            (only 0.15 from the maximum value)
             adapt_extremum_to_add_cond: should the extrumum be adapted to the additional condition
             additional_condition: additional condition to determine which bonds are relevant
-        Returns: [-float('inf'), min(max_icohp*0.15,-0.1)]
+
+        Returns: [-inf, min(strongest_icohp*0.15,-noise_cutoff)] / [max(strongest_icohp*0.15, noise_cutoff),inf]
         """
-        # TODO: make it work for COOPs/COBIs
+
         if not adapt_extremum_to_add_cond or additional_condition == 0:
             extremum_based = icohpcollection.extremum_icohpvalue(summed_spin_channels=True) * percentage
         elif additional_condition == 1:
@@ -1150,15 +1175,16 @@ def _get_limit_from_extremum(
                 if (val1 < 0.0 < val2) or (val2 < 0.0 < val1):
                     list_icohps.append(value.summed_icohp)
 
-            extremum_based = min(list_icohps) * percentage
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
 
         elif additional_condition == 2:
             # NO_ELEMENT_TO_SAME_ELEMENT_BONDS
             list_icohps = []
             for value in icohpcollection._icohplist.values():
                 if value._atom1.rstrip("0123456789") != value._atom2.rstrip("0123456789"):
                     list_icohps.append(value.summed_icohp)
-            extremum_based = min(list_icohps) * percentage
+
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
 
         elif additional_condition == 3:
             # ONLY_ANION_CATION_BONDS_AND_NO_ELEMENT_TO_SAME_ELEMENT_BONDS = 3
@@ -1175,13 +1201,17 @@ def _get_limit_from_extremum(
                     and value._atom1.rstrip("0123456789") != value._atom2.rstrip("0123456789")
                 ):
                     list_icohps.append(value.summed_icohp)
-            extremum_based = min(list_icohps) * percentage
+
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
+
         elif additional_condition == 4:
             list_icohps = []
             for value in icohpcollection._icohplist.values():
                 if value._atom1.rstrip("0123456789") == "O" or value._atom2.rstrip("0123456789") == "O":
                     list_icohps.append(value.summed_icohp)
-            extremum_based = min(list_icohps) * percentage
+
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
+
         elif additional_condition == 5:
             # DO_NOT_CONSIDER_ANION_CATION_BONDS=5
             list_icohps = []
@@ -1193,7 +1223,8 @@ def _get_limit_from_extremum(
 
                 if (val1 > 0.0 and val2 > 0.0) or (val1 < 0.0 and val2 < 0.0):
                     list_icohps.append(value.summed_icohp)
-            extremum_based = min(list_icohps) * percentage
+
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
 
         elif additional_condition == 6:
             # ONLY_CATION_CATION_BONDS=6
@@ -1206,13 +1237,17 @@ def _get_limit_from_extremum(
 
                 if val1 > 0.0 and val2 > 0.0:
                     list_icohps.append(value.summed_icohp)
-            extremum_based = min(list_icohps) * percentage
 
-        # if not self.are_coops:
-        max_here = min(extremum_based, -0.1)
-        return -float("inf"), max_here
-        # else:
-        #    return extremum_based, 100000
+            extremum_based = self._adapt_extremum_to_add_cond(list_icohps, percentage)
+
+        if not self.are_coops and not self.are_cobis:
+            max_here = min(extremum_based, -self.noise_cutoff) if self.noise_cutoff is not None else extremum_based
+            return -float("inf"), max_here
+        if self.are_coops or self.are_cobis:
+            min_here = max(extremum_based, self.noise_cutoff) if self.noise_cutoff is not None else extremum_based
+            return min_here, float("inf")
+
+        return None
 
 
 class LobsterLightStructureEnvironments(LightStructureEnvironments):

pymatgen/io/lobster/tests/test_lobsterenv.py

@@ -125,6 +125,29 @@ def setUp(self):
             structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_353.gz")),
             additional_condition=6,
         )
+        # coop / cobi
+        self.chemenvlobster1_coop_NaCl = LobsterNeighbors(
+            are_coops=True,
+            filename_ICOHP=os.path.join(test_dir_env, "ICOOPLIST.lobster.NaCl.gz"),
+            structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.NaCl.gz")),
+            additional_condition=1,
+            noise_cutoff=None,
+        )
+
+        self.chemenvlobster1_cobi_NaCl = LobsterNeighbors(
+            are_coops=True,
+            filename_ICOHP=os.path.join(test_dir_env, "ICOBILIST.lobster.NaCl.gz"),
+            structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.NaCl.gz")),
+            additional_condition=1,
+            noise_cutoff=None,
+        )
+
+        self.chemenvlobster1_cobi_mp470 = LobsterNeighbors(
+            are_coops=True,
+            filename_ICOHP=os.path.join(test_dir_env, "ICOBILIST.lobster.mp_470.gz"),
+            structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_470.gz")),
+            additional_condition=1,
+        )
 
         # TODO: use charge instead of valence
         self.chemenvlobster1_charges = LobsterNeighbors(
@@ -135,6 +158,26 @@ def setUp(self):
             filename_CHARGE=os.path.join(test_dir_env, "CHARGE.lobster.mp-353.gz"),
             additional_condition=1,
         )
+        self.chemenvlobster1_charges_noisecutoff = LobsterNeighbors(
+            are_coops=False,
+            filename_ICOHP=os.path.join(test_dir_env, "ICOHPLIST.lobster.mp_632319.gz"),
+            structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_632319.gz")),
+            valences_from_charges=True,
+            filename_CHARGE=os.path.join(test_dir_env, "CHARGE.lobster.mp_632319.gz"),
+            additional_condition=1,
+            perc_strength_ICOHP=0.05,
+            noise_cutoff=0.1,
+        )
+        self.chemenvlobster1_charges_wo_noisecutoff = LobsterNeighbors(
+            are_coops=False,
+            filename_ICOHP=os.path.join(test_dir_env, "ICOHPLIST.lobster.mp_632319.gz"),
+            structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_632319.gz")),
+            valences_from_charges=True,
+            filename_CHARGE=os.path.join(test_dir_env, "CHARGE.lobster.mp_632319.gz"),
+            additional_condition=1,
+            perc_strength_ICOHP=0.05,
+            noise_cutoff=None,
+        )
         self.chemenvlobster1_charges_loewdin = LobsterNeighbors(
             are_coops=False,
             filename_ICOHP=os.path.join(test_dir_env, "ICOHPLIST.lobster.mp_353.gz"),
@@ -218,17 +261,6 @@ def setUp(self):
             adapt_extremum_to_add_cond=True,
         )
 
-    def test_use_of_coop(self):
-        with pytest.raises(ValueError, match="Algorithm only works correctly for ICOHPLIST.lobster"):
-            _ = LobsterNeighbors(
-                are_coops=True,
-                filename_ICOHP=os.path.join(test_dir_env, "ICOHPLIST.lobster.mp_353.gz"),
-                structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_353.gz")),
-                valences_from_charges=True,
-                filename_CHARGE=os.path.join(test_dir_env, "CHARGE.lobster.mp-353.gz"),
-                additional_condition=1,
-            )
-
     def test_cation_anion_mode_without_ions(self):
         with pytest.raises(
             ValueError, match="Valences cannot be assigned, additional_conditions 1, 3, 5 and 6 will not work"
@@ -332,6 +364,24 @@ def test_get_nn_info(self):
             )
             == 2
         )
+        assert (
+            len(
+                self.chemenvlobster1_charges_noisecutoff.get_nn(
+                    structure=self.chemenvlobster1_charges_noisecutoff.structure,
+                    n=1,
+                )
+            )
+            == 0
+        )
+        assert (
+            len(
+                self.chemenvlobster1_charges_wo_noisecutoff.get_nn(
+                    structure=self.chemenvlobster1_charges_wo_noisecutoff.structure,
+                    n=1,
+                )
+            )
+            == 8
+        )
         # NO_ELEMENT_TO_SAME_ELEMENT_BONDS = 2
         assert (
             len(
@@ -453,6 +503,36 @@ def test_get_nn_info(self):
             == 2
         )
 
+        assert (
+            len(
+                self.chemenvlobster1_coop_NaCl.get_nn(
+                    structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.NaCl.gz")),
+                    n=0,
+                )
+            )
+            == 6
+        )
+
+        assert (
+            len(
+                self.chemenvlobster1_cobi_NaCl.get_nn(
+                    structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.NaCl.gz")),
+                    n=0,
+                )
+            )
+            == 6
+        )
+
+        assert (
+            len(
+                self.chemenvlobster1_cobi_mp470.get_nn(
+                    structure=Structure.from_file(os.path.join(test_dir_env, "POSCAR.mp_470.gz")),
+                    n=3,
+                )
+            )
+            == 3
+        )
+
         # NO_ELEMENT_TO_SAME_ELEMENT_BONDS = 2
         assert (
             len(