diff --git a/pymatgen/io/lobster/lobsterenv.py b/pymatgen/io/lobster/lobsterenv.py index f24ec340f54..d3b201286e2 100644 --- a/pymatgen/io/lobster/lobsterenv.py +++ b/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 @@ -1121,6 +1131,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, @@ -1134,12 +1157,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: @@ -1154,7 +1179,7 @@ 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 @@ -1162,7 +1187,8 @@ def _get_limit_from_extremum( 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 @@ -1179,13 +1205,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 = [] @@ -1197,7 +1227,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 @@ -1210,13 +1241,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): diff --git a/pymatgen/io/lobster/tests/test_lobsterenv.py b/pymatgen/io/lobster/tests/test_lobsterenv.py index 8da3414d128..43654869693 100644 --- a/pymatgen/io/lobster/tests/test_lobsterenv.py +++ b/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( diff --git a/test_files/cohp/environments/CHARGE.lobster.mp_632319.gz b/test_files/cohp/environments/CHARGE.lobster.mp_632319.gz new file mode 100755 index 00000000000..2e8b3834b52 Binary files /dev/null and b/test_files/cohp/environments/CHARGE.lobster.mp_632319.gz differ diff --git a/test_files/cohp/environments/ICOBILIST.lobster.mp_470.gz b/test_files/cohp/environments/ICOBILIST.lobster.mp_470.gz new file mode 100755 index 00000000000..8f6a116d214 Binary files /dev/null and b/test_files/cohp/environments/ICOBILIST.lobster.mp_470.gz differ diff --git a/test_files/cohp/environments/ICOHPLIST.lobster.mp_632319.gz b/test_files/cohp/environments/ICOHPLIST.lobster.mp_632319.gz new file mode 100755 index 00000000000..6d77ecdffd9 Binary files /dev/null and b/test_files/cohp/environments/ICOHPLIST.lobster.mp_632319.gz differ diff --git a/test_files/cohp/environments/POSCAR.mp_470.gz b/test_files/cohp/environments/POSCAR.mp_470.gz new file mode 100755 index 00000000000..7eb8d52eff5 Binary files /dev/null and b/test_files/cohp/environments/POSCAR.mp_470.gz differ diff --git a/test_files/cohp/environments/POSCAR.mp_632319.gz b/test_files/cohp/environments/POSCAR.mp_632319.gz new file mode 100755 index 00000000000..7bbc6a50c77 Binary files /dev/null and b/test_files/cohp/environments/POSCAR.mp_632319.gz differ