From 163f0699939be314ddb6b92645a16bdbf65975d9 Mon Sep 17 00:00:00 2001 From: naik-aakash Date: Sat, 20 Jul 2024 18:49:29 +0530 Subject: [PATCH 1/3] move to src layout --- pyproject.toml | 5 +- {lobsterpy => src/lobsterpy}/__init__.py | 0 {lobsterpy => src/lobsterpy}/cli.py | 0 {lobsterpy => src/lobsterpy}/cohp/__init__.py | 0 {lobsterpy => src/lobsterpy}/cohp/analyze.py | 3464 ++++++++--------- {lobsterpy => src/lobsterpy}/cohp/describe.py | 1714 ++++---- .../lobsterpy}/featurize/__init__.py | 0 .../lobsterpy}/featurize/batch.py | 0 .../lobsterpy}/featurize/core.py | 0 .../lobsterpy}/featurize/utils.py | 0 .../lobsterpy}/plotting/__init__.py | 0 .../lobsterpy}/plotting/layout_dicts.py | 282 +- .../plotting/lobsterpy_base.mplstyle | 0 .../lobsterpy}/structuregraph/__init__.py | 0 .../lobsterpy}/structuregraph/graph.py | 0 15 files changed, 2733 insertions(+), 2732 deletions(-) rename {lobsterpy => src/lobsterpy}/__init__.py (100%) rename {lobsterpy => src/lobsterpy}/cli.py (100%) rename {lobsterpy => src/lobsterpy}/cohp/__init__.py (100%) rename {lobsterpy => src/lobsterpy}/cohp/analyze.py (97%) rename {lobsterpy => src/lobsterpy}/cohp/describe.py (97%) rename {lobsterpy => src/lobsterpy}/featurize/__init__.py (100%) rename {lobsterpy => src/lobsterpy}/featurize/batch.py (100%) rename {lobsterpy => src/lobsterpy}/featurize/core.py (100%) rename {lobsterpy => src/lobsterpy}/featurize/utils.py (100%) rename {lobsterpy => src/lobsterpy}/plotting/__init__.py (100%) rename {lobsterpy => src/lobsterpy}/plotting/layout_dicts.py (96%) rename {lobsterpy => src/lobsterpy}/plotting/lobsterpy_base.mplstyle (100%) rename {lobsterpy => src/lobsterpy}/structuregraph/__init__.py (100%) rename {lobsterpy => src/lobsterpy}/structuregraph/graph.py (100%) diff --git a/pyproject.toml b/pyproject.toml index b62bc95b..d8e84f50 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -6,8 +6,9 @@ requires = [ build-backend = "setuptools.build_meta" [tool.setuptools.packages.find] +where = ["src"] include = ["lobsterpy*"] -exclude = ["docs*", "tests*", "paper*"] +#exclude = ["docs*", "tests*", "paper*"] [project] name = "lobsterpy" @@ -165,7 +166,7 @@ exclude_lines = [ ] [tool.coverage.run] -include = ["lobsterpy/*"] +include = ["src/"] omit = [ # omit anything in test directory anywhere "tests/*", diff --git a/lobsterpy/__init__.py b/src/lobsterpy/__init__.py similarity index 100% rename from lobsterpy/__init__.py rename to src/lobsterpy/__init__.py diff --git a/lobsterpy/cli.py b/src/lobsterpy/cli.py similarity index 100% rename from lobsterpy/cli.py rename to src/lobsterpy/cli.py diff --git a/lobsterpy/cohp/__init__.py b/src/lobsterpy/cohp/__init__.py similarity index 100% rename from lobsterpy/cohp/__init__.py rename to src/lobsterpy/cohp/__init__.py diff --git a/lobsterpy/cohp/analyze.py b/src/lobsterpy/cohp/analyze.py similarity index 97% rename from lobsterpy/cohp/analyze.py rename to src/lobsterpy/cohp/analyze.py index 38c6a260..ee9be0e8 100644 --- a/lobsterpy/cohp/analyze.py +++ b/src/lobsterpy/cohp/analyze.py @@ -1,1732 +1,1732 @@ -# Copyright (c) lobsterpy development team -# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License - -"""This module defines classes to analyze the COHPs automatically.""" - -from __future__ import annotations - -import warnings -from collections import Counter -from itertools import combinations_with_replacement, permutations -from pathlib import Path - -import numpy as np -from pymatgen.analysis.bond_valence import BVAnalyzer -from pymatgen.core.structure import Structure -from pymatgen.electronic_structure.cohp import CompleteCohp -from pymatgen.electronic_structure.core import Spin -from pymatgen.electronic_structure.dos import LobsterCompleteDos -from pymatgen.io.lobster import ( - Bandoverlaps, - Charge, - Doscar, - Icohplist, - Lobsterin, - Lobsterout, - MadelungEnergies, -) -from pymatgen.io.lobster.lobsterenv import LobsterNeighbors -from pymatgen.io.vasp.outputs import Vasprun -from pymatgen.symmetry.analyzer import SpacegroupAnalyzer - - -class Analysis: - """ - Class to automatically analyze COHP/COOP/COBI populations from Lobster. - - :param are_cobis: bool indicating if file contains COBI/ICOBI data - :param are_coops: bool indicating if file contains COOP/ICOOP data - :param cutoff_icohp: Cutoff in percentage for evaluating neighbors based on ICOHP values. - cutoff_icohp*max_icohp limits the number of considered neighbours for evaluating environments. - :param path_to_cohpcar: path to `COHPCAR.lobster` or `COBICAR.lobster` or `COOPCAR.lobster` . - :param path_to_charge: path to `CHARGE.lobster`. - :param path_to_icohplist: path to `ICOHPLIST.lobster` or `ICOBILIST.lobster` or `ICOOPLIST.lobster`. - :param path_to_poscar: path to structure (e.g., `POSCAR` or `POSCAR.lobster`) - :param path_to_madelung: path to `MadelungEnergies.lobster`. - :param charge_obj: pymatgen lobster.io.charge object - :param completecohp_obj: pymatgen.electronic_structure.cohp.CompleteCohp object - :param icohplist_obj: pymatgen lobster.io.Icohplist object - :param madelung_obj: pymatgen lobster.io.MadelungEnergies object - :param noise_cutoff: Sets the lower limit tolerance for ICOHPs or ICOOPs or ICOBIs considered - in analysis. - :param orbital_cutoff: Sets the minimum percentage for the orbital contribution considered to be - relevant in orbital resolved analysis. (Affects only when orbital_resolved argument is set to True) - Set it to 0 to get results of all orbitals in the detected relevant bonds. Default is to 0.05 i.e. - only analyzes if orbital contribution is 5 % or more. - :param orbital_resolved: bool indicating whether orbital wise analysis is performed - :param type_charge: If no path_to_charge is given, Valences will be used. Otherwise, Mulliken charges. - Löwdin charges cannot be selected at the moment. - :param which_bonds: Selects kinds of bonds that are analyzed. `cation-anion` is the default. - Alternatively, `all` bonds can also be selected. Support to other kinds of bonds will be - added soon. - :param summed_spins: if True, COHP `Spin.up` and `Spin.down` populations will be summed - :param start: sets the lower limit of energy for evaluation of bonding and antibonding - percentages below efermi. Defaults to None (i.e., all populations below efermi are included) - - Attributes: - - condensed_bonding_analysis: dict including a summary of the most important bonding properties - - final_dict_bonds: dict including information on ICOHPs per bond type - - final_dict_ions: dict including information on environments of cations - - chemenv: pymatgen.io.lobster.lobsterenv.LobsterNeighbors object - - lse: LightStructureEnvironment from pymatgen - - anion_types: Set of Element objects from pymatgen - - list_equivalent_sites: list of site indices of sites that indicate which sites are equivalent - e.g., [0 1 2 2 2] where site 0, 1, 2 indicate sites that are independent from each other - - seq_cohps: list of cohps - - seq_coord_ions: list of co-ordination environment strings for each cation - - seq_equivalent_sites: seq of inequivalent sites - - seq_ineq_ions: seq of inequivalent cations/sites in the structure - - seq_infos_bonds (list): information on cation anion bonds (lists - of pymatgen.io.lobster.lobsterenv.ICOHPNeighborsInfo) - - spg: space group information - - structure: Structure object - - """ - - def __init__( - self, - path_to_poscar: str | Path | None, - path_to_icohplist: str | Path | None, - path_to_cohpcar: str | Path | None, - path_to_charge: str | Path | None = None, - path_to_madelung: str | Path | None = None, - icohplist_obj: Icohplist | None = None, - completecohp_obj: CompleteCohp | None = None, - charge_obj: Charge | None = None, - madelung_obj: MadelungEnergies | None = None, - are_cobis: bool = False, - are_coops: bool = False, - cutoff_icohp: float = 0.1, - noise_cutoff: float = 0.1, - orbital_cutoff: float = 0.05, - orbital_resolved: bool = False, - start: float | None = None, - summed_spins: bool = True, - type_charge: str | None = None, - which_bonds: str = "cation-anion", - ): - """ - Initialize automatic bonding analysis. - - :param are_cobis: bool indicating if file contains COBI/ICOBI data - :param are_coops: bool indicating if file contains COOP/ICOOP data - :param cutoff_icohp: Cutoff in percentage for evaluating neighbors based on ICOHP values. - cutoff_icohp*max_icohp limits the number of considered neighbours for evaluating environments. - :param path_to_cohpcar: path to `COHPCAR.lobster` or `COBICAR.lobster` or `COOPCAR.lobster` . - :param path_to_charge: path to `CHARGE.lobster`. - :param path_to_icohplist: path to `ICOHPLIST.lobster` or `ICOBILIST.lobster` or `ICOOPLIST.lobster`. - :param path_to_poscar: path to structure (e.g., `POSCAR` or `POSCAR.lobster`) - :param path_to_madelung: path to `MadelungEnergies.lobster`. - :param charge_obj: pymatgen lobster.io.charge object (Optional) - :param completecohp_obj: pymatgen.electronic_structure.cohp.CompleteCohp object - :param icohplist_obj: pymatgen lobster.io.Icohplist object - :param madelung_obj: pymatgen lobster.io.MadelungEnergies object - :param noise_cutoff: Sets the lower limit tolerance for ICOHPs or ICOOPs or ICOBIs considered - in analysis. - :param orbital_cutoff: Sets the minimum percentage for the orbital contribution considered to be - relevant in orbital resolved analysis. (Affects only when orbital_resolved argument is set to True) - Set it to 0 to get results of all orbitals in the detected relevant bonds. Default is to 0.05 i.e. - only analyzes if orbital contribution is 5 % or more. - :param orbital_resolved: bool indicating whether orbital wise analysis is performed - :param type_charge: If no path_to_charge is given, Valences will be used. Otherwise, Mulliken charges. - Löwdin charges cannot be selected at the moment. - :param which_bonds: Selects kinds of bonds that are analyzed. `cation-anion` is the default. - Alternatively, `all` bonds can also be selected. Support to other kinds of bonds will be - added soon. - :param summed_spins: if True, COHP `Spin.up` and `Spin.down` populations will be summed - :param start: sets the lower limit of energy for evaluation of bonding and antibonding - percentages below efermi. Defaults to None (i.e., all populations below efermi are included) - - """ - self.start = start - self.completecohp_obj = completecohp_obj - self.icohplist_obj = icohplist_obj - # checks to ensure LobsterEnv inputs are not duplicated in case users provide both path and obj - if self.completecohp_obj is not None and self.icohplist_obj is not None: - self.path_to_poscar = None - self.path_to_cohpcar = None - self.path_to_icohplist = None - else: - self.path_to_poscar = path_to_poscar - self.path_to_icohplist = path_to_icohplist - self.path_to_cohpcar = path_to_cohpcar - self.which_bonds = which_bonds - self.cutoff_icohp = cutoff_icohp - self.orbital_cutoff = orbital_cutoff - self.path_to_charge = path_to_charge - self.charge_obj = charge_obj - self.path_to_madelung = path_to_madelung - self.madelung_obj = madelung_obj - self.are_cobis = are_cobis - self.are_coops = are_coops - self.noise_cutoff = noise_cutoff - self.setup_env() - self.get_information_all_bonds(summed_spins=summed_spins) - self.orbital_resolved = orbital_resolved - - # This determines how cations and anions - if path_to_charge is None and charge_obj is None: - self.type_charge = "Valences" - else: - if type_charge is None: - self.type_charge = "Mulliken" - elif type_charge == "Mulliken": - self.type_charge = "Mulliken" - elif type_charge == "Löwdin": - raise ValueError("Only Mulliken charges can be used here at the moment. Implementation will follow.") - else: - self.type_charge = "Valences" - print("type_charge cannot be read! Please use Mulliken/Löwdin. Now, we will use valences") - - self.set_condensed_bonding_analysis() - self.set_summary_dicts() - self.path_to_madelung = path_to_madelung - - def setup_env(self): - """ - Set up the light structure environments based on COHPs using this method. - - Returns: - None - - """ - self.structure = ( - Structure.from_file(self.path_to_poscar) if self.path_to_poscar else self.completecohp_obj.structure - ) - sga = SpacegroupAnalyzer(structure=self.structure) - symmetry_dataset = sga.get_symmetry_dataset() - equivalent_sites = symmetry_dataset["equivalent_atoms"] - self.list_equivalent_sites = equivalent_sites - self.seq_equivalent_sites = list(set(equivalent_sites)) - self.spg = symmetry_dataset["international"] - - if self.which_bonds == "cation-anion": - try: - self.chemenv = LobsterNeighbors( - filename_icohp=self.path_to_icohplist, - obj_icohp=self.icohplist_obj, - structure=self.structure, - additional_condition=1, - perc_strength_icohp=self.cutoff_icohp, - filename_charge=self.path_to_charge, - obj_charge=self.charge_obj, - valences=None, - valences_from_charges=True, - adapt_extremum_to_add_cond=True, - are_cobis=self.are_cobis, - are_coops=self.are_coops, - noise_cutoff=self.noise_cutoff, - ) - except ValueError as err: - if ( - str(err) == "min() arg is an empty sequence" - or str(err) == "All valences are equal to 0, additional_conditions 1, 3, 5 and 6 will not work" - ): - raise ValueError( - "Consider switching to an analysis of all bonds and not only cation-anion bonds." - " It looks like no cations are detected." - ) - raise err - elif self.which_bonds == "all": - # raise ValueError("only cation anion bonds implemented so far") - self.chemenv = LobsterNeighbors( - filename_icohp=self.path_to_icohplist, - obj_icohp=self.icohplist_obj, - structure=self.structure, - additional_condition=0, - perc_strength_icohp=self.cutoff_icohp, - filename_charge=self.path_to_charge, - obj_charge=self.charge_obj, - valences_from_charges=True, - adapt_extremum_to_add_cond=True, - are_cobis=self.are_cobis, - are_coops=self.are_coops, - noise_cutoff=self.noise_cutoff, - ) - - else: - raise ValueError("only cation anion and all bonds implemented so far") - - # determine cations and anions - try: - if self.which_bonds == "cation-anion": - self.lse = self.chemenv.get_light_structure_environment(only_cation_environments=True) - elif self.which_bonds == "all": - self.lse = self.chemenv.get_light_structure_environment(only_cation_environments=False) - except ValueError: - - class Lse: - """Test class when error was raised.""" - - def __init__(self, chemenv, valences=None): - """ - Test class when error was raised. - - :param chemenv: LobsterNeighbors object - :param valences: list of valences - - """ - if valences is None: - self.coordination_environments = [] - for coord in chemenv: - if len(coord) > 0: - self.coordination_environments.append([{"ce_symbol": str(len(coord))}]) - else: - self.coordination_environments.append([{"ce_symbol": None}]) - else: - self.coordination_environments = [] - - for val, coord in zip(valences, chemenv): - if val >= 0.0 and len(coord) > 0: - self.coordination_environments.append([{"ce_symbol": str(len(coord))}]) - else: - self.coordination_environments.append([{"ce_symbol": None}]) - - if self.which_bonds == "all": - self.lse = Lse(self.chemenv.list_coords) - elif self.which_bonds == "cation-anion": - # make a new list - self.lse = Lse(self.chemenv.list_coords, self.chemenv.valences) - - def get_information_all_bonds(self, summed_spins: bool = True): - """ - Gather all information on the bonds within the compound with this method. - - Returns: - None - - """ - if self.which_bonds == "cation-anion": - # this will only analyze cation anion bonds which simplifies the analysis - self.seq_ineq_ions = [] - self.seq_coord_ions = [] - self.seq_infos_bonds = [] - self.seq_labels_cohps = [] - self.seq_cohps = [] - # only_bonds_to - - self.anion_types = self.chemenv.anion_types - for ice, ce in enumerate(self.lse.coordination_environments): - # only look at inequivalent sites (use of symmetry to speed everything up!)! - # only look at those cations that have cation-anion bonds - if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None: - self.seq_ineq_ions.append(ice) - - self.seq_coord_ions.append(ce[0]["ce_symbol"]) - self.seq_infos_bonds.append(self.chemenv.get_info_icohps_to_neighbors([ice])) - - aniontype_labels = [] - aniontype_cohps = [] - - # go through all anions in the structure! - for anion in self.anion_types: - # get labels and summed cohp objects - labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors( - path_to_cohpcar=self.path_to_cohpcar, - obj_cohpcar=self.completecohp_obj, - isites=[ice], - summed_spin_channels=summed_spins, - per_bond=False, - only_bonds_to=[str(anion)], - ) - - aniontype_labels.append(labels) - aniontype_cohps.append(summedcohps) - - self.seq_labels_cohps.append(aniontype_labels) - self.seq_cohps.append(aniontype_cohps) - - elif self.which_bonds == "all": - # this will only analyze all bonds - - self.seq_ineq_ions = [] - self.seq_coord_ions = [] - self.seq_infos_bonds = [] - self.seq_labels_cohps = [] - self.seq_cohps = [] - # only_bonds_to - self.elements = self.structure.composition.elements - # self.anion_types = self.chemenv.get_anion_types() - for ice, ce in enumerate(self.lse.coordination_environments): - # only look at inequivalent sites (use of symmetry to speed everything up!)! - # only look at those cations that have cation-anion bonds - if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None: - self.seq_ineq_ions.append(ice) - self.seq_coord_ions.append(ce[0]["ce_symbol"]) - self.seq_infos_bonds.append(self.chemenv.get_info_icohps_to_neighbors([ice])) - - type_labels = [] - type_cohps = [] - - for element in self.elements: - # get labels and summed cohp objects - labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors( - path_to_cohpcar=self.path_to_cohpcar, - obj_cohpcar=self.completecohp_obj, - isites=[ice], - onlycation_isites=False, - summed_spin_channels=summed_spins, - per_bond=False, - only_bonds_to=[str(element)], - ) - - type_labels.append(labels) - type_cohps.append(summedcohps) - - self.seq_labels_cohps.append(type_labels) - self.seq_cohps.append(type_cohps) - - def get_site_bond_resolved_labels(self): - """ - Return relevant bond labels for each symmetrically independent site. - - Returns: - dict with bond labels for each site, e.g. - {'Na1: Na-Cl': ['21', '23', '24', '27', '28', '30']} - - """ - bonds = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore - labels = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore - for inx, bond_info in enumerate(self.seq_infos_bonds): - for ixx, val in enumerate(bond_info.atoms): - label_srt = sorted(val.copy()) - bonds[inx].append( - self.structure.sites[bond_info.central_isites[0]].species_string - + str(bond_info.central_isites[0] + 1) - + ": " - + label_srt[0].strip("0123456789") - + "-" - + label_srt[1].strip("0123456789") - ) - labels[inx].append(bond_info.labels[ixx]) - - label_data = {} - for indx, atom_pairs in enumerate(bonds): - searched_atom_pairs = set(atom_pairs) - for search_item in searched_atom_pairs: - indices = [i for i, pair in enumerate(atom_pairs) if pair == search_item] - filtered_bond_label_list = [labels[indx][i] for i in indices] - label_data.update({search_item: filtered_bond_label_list}) - - return label_data - - def _get_orbital_resolved_data( - self, - nameion: str, - iion: int, - labels: list[str], - bond_resolved_labels: dict[str, list[str]], - type_pop: str, - ): - """ - Retrieve orbital-wise analysis data. - - :param nameion: name of symmetrically relevant cation or anion - :param iion: index of symmetrically relevant cation or anion - :param labels: list of bond label names - :param bond_resolved_labels: dict of bond labels from ICOHPLIST resolved for each bond - :param type_pop: population type analyzed. e.g. COHP or COOP or COBI - - Returns: - dict consisting of relevant orbitals (contribution > 5 % to overall ICOHP or ICOBI or ICOOP), - bonding and antibonding percentages with bond label names as keys. - """ - orb_resolved_bond_info = {} - for label in labels: - if label is not None: - bond_resolved_label_key = nameion + str(iion + 1) + ":" + label.split("x")[-1] - bond_labels = bond_resolved_labels[bond_resolved_label_key] - orb_combinations = self._get_orb_combinations() - grouped_orb_pairs = self._group_orb_pairs(bond_label=bond_labels[0], orb_combinations=orb_combinations) - # available_orbitals = list(self.chemenv.completecohp.orb_res_cohp[bond_labels[0]].keys()) - # initialize empty list to store orb paris for bonding, - # antibonding integrals and percentages - bndg_orb_pair_list = [] - bndg_orb_integral_list = [] - bndg_orb_perc_list = [] - bndg_orb_icohp_list = [] - antibndg_orb_integral_list = [] - antibndg_orb_perc_list = [] - antibndg_orb_pair_list = [] - antibndg_orb_icohp_list = [] - - # get total summed cohps using label list - cohp_summed = self.chemenv.completecohp.get_summed_cohp_by_label_list(label_list=bond_labels) - if type_pop.lower() == "cohp": - ( - antibndg_tot, - per_anti_tot, - bndg_tot, - per_bndg_tot, - ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed, start=self.start) - else: - ( - bndg_tot, - per_bndg_tot, - antibndg_tot, - per_anti_tot, - ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed, start=self.start) - - orb_bonding_dict_data = {} # type: ignore - # For each orbital collect the contributions of summed bonding - # and antibonding interactions separately - for orb in grouped_orb_pairs: - mapped_bond_labels = self._get_bond_orb_label_mapped_list( - bond_labels=bond_labels, orb_pair=grouped_orb_pairs, root_orb_pair=orb - ) - # for orb in available_orbitals: - cohp_summed_orb = self.chemenv.completecohp.get_summed_cohp_by_label_and_orbital_list( - label_list=mapped_bond_labels, orbital_list=grouped_orb_pairs[orb] * len(bond_labels) - ) - - if type_pop.lower() == "cohp": - ( - antibndg_orb, - per_anti_orb, - bndg_orb, - per_bndg_orb, - ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed_orb, start=self.start) - else: - ( - bndg_orb, - per_bndg_orb, - antibndg_orb, - per_anti_orb, - ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed_orb, start=self.start) - - # replace nan values with zero (tackle numerical integration issues) - bndg_orb = bndg_orb if not np.isnan(bndg_orb) else 0 - per_bndg_orb = per_bndg_orb if not np.isnan(per_bndg_orb) else 0 - bndg_tot = bndg_tot if not np.isnan(bndg_tot) else 0 - per_bndg_tot = per_bndg_tot if not np.isnan(per_bndg_tot) else 0 - # skip collecting orb contributions if no summed bonding contribution exists - if bndg_tot > 0: - orb_icohps_bndg = [] - for bond_label in bond_labels: - for sub_orb in grouped_orb_pairs[orb]: - orb_icohp_bn = self.chemenv.Icohpcollection.get_icohp_by_label( - label=bond_label, orbitals=sub_orb - ) - orb_icohps_bndg.append(orb_icohp_bn) - bndg_orb_pair_list.append(orb) - bndg_orb_icohp_list.append(orb_icohps_bndg) - bndg_orb_integral_list.append(bndg_orb) - bndg_orb_perc_list.append(per_bndg_orb) - - # replace nan values with zero (tackle numerical integration issues) - antibndg_orb = antibndg_orb if not np.isnan(antibndg_orb) else 0 - per_anti_orb = per_anti_orb if not np.isnan(per_anti_orb) else 0 - antibndg_tot = antibndg_tot if not np.isnan(antibndg_tot) else 0 - per_anti_tot = per_anti_tot if not np.isnan(per_anti_tot) else 0 - # skip collecting orb contributions if no summed antibonding contribution exists - if antibndg_tot > 0: - orb_icohps_anti = [] - for bond_label in bond_labels: - for sub_orb in grouped_orb_pairs[orb]: - orb_icohp_an = self.chemenv.Icohpcollection.get_icohp_by_label( - label=bond_label, orbitals=sub_orb - ) - orb_icohps_anti.append(orb_icohp_an) - - antibndg_orb_pair_list.append(orb) - antibndg_orb_icohp_list.append(orb_icohps_anti) - antibndg_orb_integral_list.append(antibndg_orb) - antibndg_orb_perc_list.append(per_anti_orb) - - # Populate the dictionary with relevant orbitals for bonding interactions - for inx, bndg_orb_pair in enumerate(bndg_orb_pair_list): - bndg_contri_perc = round(bndg_orb_integral_list[inx] / sum(bndg_orb_integral_list), 2) - # filter out very small bonding interactions ( self.orbital_cutoff: - if bndg_orb_pair in orb_bonding_dict_data: - orb_bonding_dict_data[bndg_orb_pair].update( - { - "orb_contribution_perc_bonding": bndg_contri_perc, - "bonding": { - "integral": bndg_orb_integral_list[inx], - "perc": bndg_orb_perc_list[inx], - }, - } - ) - else: - orb_bonding_dict_data[bndg_orb_pair] = { - f"I{type_pop}_mean": round(np.mean(bndg_orb_icohp_list[inx]), 4), - f"I{type_pop}_sum": round(np.sum(bndg_orb_icohp_list[inx]), 4), - "orb_contribution_perc_bonding": round( - bndg_orb_integral_list[inx] / sum(bndg_orb_integral_list), - 2, - ), - "bonding": { - "integral": bndg_orb_integral_list[inx], - "perc": bndg_orb_perc_list[inx], - }, - "relevant_sub_orbitals": grouped_orb_pairs[bndg_orb_pair], - } - - # Populate the dictionary with relevant orbitals for antibonding interactions - for inx, antibndg_orb_pair in enumerate(antibndg_orb_pair_list): - antibndg_contri_perc = round( - antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), - 2, - ) - # filter out very small antibonding interactions ( self.orbital_cutoff: - if antibndg_orb_pair in orb_bonding_dict_data: - orb_bonding_dict_data[antibndg_orb_pair].update( - { - "orb_contribution_perc_antibonding": round( - antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), - 2, - ), - "antibonding": { - "integral": antibndg_orb_integral_list[inx], - "perc": antibndg_orb_perc_list[inx], - }, - } - ) - else: - orb_bonding_dict_data[antibndg_orb_pair] = { - f"I{type_pop}_mean": round(np.mean(antibndg_orb_icohp_list[inx]), 4), - f"I{type_pop}_sum": round(np.sum(antibndg_orb_icohp_list[inx]), 4), - "orb_contribution_perc_antibonding": round( - antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), - 2, - ), - "antibonding": { - "integral": antibndg_orb_integral_list[inx], - "perc": antibndg_orb_perc_list[inx], - }, - "relevant_sub_orbitals": grouped_orb_pairs[antibndg_orb_pair], - } - - orb_bonding_dict_data["relevant_bonds"] = bond_labels # type: ignore - - orb_resolved_bond_info[bond_resolved_label_key] = orb_bonding_dict_data - - return orb_resolved_bond_info - - def _get_bond_resolved_data_stats(self, orb_resolved_bond_data: dict): - """ - Retrieve the maximum bonding and anti-bonding orbital contributions. - - :param orb_resolved_bond_data: A dictionary with orbital names as keys and corresponding bonding data - - Returns: - dict with orbital data stats the site for relevant orbitals, e.g. - {'orbital_summary_stats': {'max_bonding_contribution': {'3p-3p': 0.41}, - 'max_antibonding_contribution': {'3s-3p': 0.39}}} - - """ - # get max orbital bonding and contribution for the site - orb_pairs_bndg = [] - orb_pairs_antibndg = [] - orb_contri_bndg = [] - orb_contri_antibndg = [] - orbital_summary_stats = {"orbital_summary_stats": {}} # type: ignore - if orb_resolved_bond_data: - for orb_pair, data in orb_resolved_bond_data.items(): - if "orb_contribution_perc_bonding" in data: - orb_pairs_bndg.append(orb_pair) - orb_contri_bndg.append(data["orb_contribution_perc_bonding"]) - - if "orb_contribution_perc_antibonding" in data: - orb_pairs_antibndg.append(orb_pair) - orb_contri_antibndg.append(data["orb_contribution_perc_antibonding"]) - - if orb_contri_bndg: - max_orb_contri_bndg = max(orb_contri_bndg) - max_orb_contri_bndg_inxs = [ - inx for inx, orb_contri in enumerate(orb_contri_bndg) if orb_contri == max_orb_contri_bndg - ] - max_orb_contri_bndg_dict = {} - for inx in max_orb_contri_bndg_inxs: - max_orb_contri_bndg_dict[orb_pairs_bndg[inx]] = orb_contri_bndg[inx] - orbital_summary_stats["orbital_summary_stats"]["max_bonding_contribution"] = max_orb_contri_bndg_dict - if orb_contri_antibndg: - max_orb_contri_antibndg = max(orb_contri_antibndg) - max_antibndg_contri_inxs = [ - inx - for inx, orb_anti_per in enumerate(orb_contri_antibndg) - if orb_anti_per == max_orb_contri_antibndg - ] - max_antibndg_contri_dict = {} - for inx in max_antibndg_contri_inxs: - max_antibndg_contri_dict[orb_pairs_antibndg[inx]] = orb_contri_antibndg[inx] - orbital_summary_stats["orbital_summary_stats"]["max_antibonding_contribution"] = ( - max_antibndg_contri_dict - ) - - return orbital_summary_stats - - def get_site_orbital_resolved_labels(self): - """ - Return relevant orbitals and bond labels for each symmetrically independent site. - - Returns: - dict with bond labels for each site for relevant orbitals, e.g. - {'Na1: Na-Cl': {'3p-3s': {'bond_labels': ['21', '23', '24', '27', '28', '30'], - 'relevant_sub_orbitals': ['3py-3s', '3pz-3s', '3px-3s']}} - """ - site_bond_labels = self.get_site_bond_resolved_labels() - orb_plot_data = {atom_pair: {} for atom_pair in site_bond_labels} - if self.orbital_resolved: - for site_index, cba_data in self.condensed_bonding_analysis["sites"].items(): - for atom in cba_data["bonds"]: - for orb_pair in cba_data["bonds"][atom]["orbital_data"]: - if orb_pair not in ("orbital_summary_stats", "relevant_bonds"): - atom_pair = [cba_data["ion"], atom] - atom_pair.sort() - key = ( - self.structure.sites[site_index].species_string - + str(site_index + 1) - + ": " - + "-".join(atom_pair) - ) - label_list = site_bond_labels[key] - relevant_sub_orbitals = cba_data["bonds"][atom]["orbital_data"][orb_pair][ - "relevant_sub_orbitals" - ] - orb_plot_data[key].update( - {orb_pair: {"bond_labels": label_list, "relevant_sub_orbitals": relevant_sub_orbitals}} - ) - else: - print("Please set orbital_resolved to True when instantiating Analysis object, to get this data") - - return orb_plot_data - - @staticmethod - def _get_strenghts_for_each_bond(pairs: list[list[str]], strengths: list[float], nameion: str | None = None): - """ - Return a dictionary of bond strengths. - - :param pairs: list of list including labels for the atoms, e.g., [['O3', 'Cu1'], ['O3', 'Cu1']] - :param strengths: list that gives the icohp strengths as a float, [-1.86287, -1.86288] - :param nameion: string including the name of the cation in the list, e.g Cu1 - - Returns: - dict including inormation on icohps for each bond type, e.g. - {'Yb-Sb': [-1.59769, -2.14723, -1.7925, -1.60773, -1.80149, -2.14335]} - - - """ - dict_strenghts = {} # type: ignore - - for pair, strength in zip(pairs, strengths): - if nameion is not None: - new = [ - LobsterNeighbors._split_string(pair[0])[0], - LobsterNeighbors._split_string(pair[1])[0], - ] - new = Analysis._sort_name(new, nameion) - string_here = new[0] + "-" + new[1] - else: - new = sorted( - [ - LobsterNeighbors._split_string(pair[0])[0], - LobsterNeighbors._split_string(pair[1])[0], - ] - ) - string_here = new[0] + "-" + new[1] - - if string_here not in dict_strenghts: - dict_strenghts[string_here] = [] - dict_strenghts[string_here].append(strength) - return dict_strenghts - - @staticmethod - def _sort_name(pair: list[str], nameion: str | None = None): - """ - Place the cation first in a list of name strings. - - :param pair: ["O","Cu"] - :param nameion: "Cu" - - Returns: - will return list of str, e.g. ["Cu", "O"] - - """ - if nameion is not None: - new = [] - if pair[0] == nameion: - new.append(pair[0]) - new.append(pair[1]) - - elif pair[1] == nameion: - new.append(pair[1]) - new.append(pair[0]) - - return new - - @staticmethod - def _sort_orbital_atom_pair( - atom_pair: list[str], - label: str, - complete_cohp: CompleteCohp, - orb_pair: str, - ): - """ - Place the cation first in a list of name strings and add the associated orbital name alongside the atom name. - - :param atom_pair: list of atom pair with cation first eg., ["Cl","Na"] - :param label: LOBSTER relevant bond label eg ., "3" - :param complete_cohp: pymatgen CompleteCohp object - :param orb_pair: relevant orbital pair eg., "2p-3s" - - Returns: - will return list of str, e.g. ["Na(2p)", "Cl(3s)"] - - """ - orb_atom = {} # type: ignore - orb_pair_list = orb_pair.split("-") - # get orbital associated to the atom and store in a dict - for _inx, (site, site_orb) in enumerate(zip(complete_cohp.bonds[label]["sites"], orb_pair_list)): - if site.species_string in orb_atom: # check necessary for bonds between same atoms - orb_atom[site.species_string].append(site_orb) - else: - orb_atom[site.species_string] = [site_orb] - - orb_atom_list = [] - # add orbital name next to atom_pair - for inx, atom in enumerate(atom_pair): - # check to ensure getting 2nd orbital if bond is between same atomic species - if inx == 1 and len(orb_atom.get(atom)) > 1: # type: ignore - atom_with_orb_name = f"{atom}({orb_atom.get(atom)[1]})" # type: ignore - else: - atom_with_orb_name = f"{atom}({orb_atom.get(atom)[0]})" # type: ignore - orb_atom_list.append(atom_with_orb_name) - - return orb_atom_list - - def _get_antibdg_states(self, cohps, labels: list[str], nameion: str | None = None, limit=0.01): - """ - Return a dictionary containing information on anti-bonding states. - - e.g., similar to: {'Cu-O': True, 'Cu-F': True} - - :param cohps: list of pymatgen.electronic_structure.cohp.Cohp objects - :param labels: ['2 x Cu-O', '4 x Cu-F'] - :param nameion: string of the cation name, e.g. "Cu" - :param limit: limit to detect antibonding states - - Returns: - dict including in formation on whether antibonding interactions exist, - e.g., {'Cu-O': True, 'Cu-F': True} - - - """ - dict_antibd = {} - for label, cohp in zip(labels, cohps): - if label is not None: - if nameion is not None: - new = label.split(" ")[2].split("-") - sorted_new = self._sort_name(new, nameion) - new_label = sorted_new[0] + "-" + sorted_new[1] - else: - new = label.split(" ")[2].split("-") - sorted_new = sorted(new.copy()) - new_label = sorted_new[0] + "-" + sorted_new[1] - - antbd = cohp.has_antibnd_states_below_efermi(limit=limit) - if Spin.down in antbd: - dict_antibd[new_label] = antbd[Spin.up] or antbd[Spin.down] - else: - dict_antibd[new_label] = antbd[Spin.up] - - return dict_antibd - - def _integrate_antbdstates_below_efermi_for_set_cohps(self, labels: list[str], cohps, nameion: str): - """ - Return a dictionary containing information on antibonding states. - - .. warning:: NEEDS MORE TESTS - - It is important to note that only the energy range that has been computed can be considered - (i.e., this might not be all) - - e.g. output: {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995}, - 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}} - - :param cohps: list of pymatgen.electronic_structure.cohp.Cohp objects - :param labels: ['2 x Cu-O', '4 x Cu-F'] - :param nameion: string of the cation name, e.g. "Cu" - - Returns: - dict including in formation on whether antibonding interactions exist, - e.g., {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995}, - 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}}} - """ - dict_bd_antibd = {} - for label, cohp in zip(labels, cohps): - if label is not None: - new = label.split(" ")[2].split("-") - sorted_new = self._sort_name(new, nameion) - new_label = sorted_new[0] + "-" + sorted_new[1] - if not self.are_cobis and not self.are_coops: - ( - integral, - perc, - integral2, - perc2, - ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start) - else: - ( - integral2, - perc2, - integral, - perc, - ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start) - - if integral == 0 and integral2 != 0.0: - dict_bd_antibd[new_label] = { - "bonding": {"integral": integral2, "perc": perc2}, - "antibonding": {"integral": integral, "perc": 0.0}, - } - elif integral2 == 0.0 and integral != 0.0: - dict_bd_antibd[new_label] = { - "bonding": {"integral": integral2, "perc": 0.0}, - "antibonding": {"integral": integral, "perc": perc}, - } - elif integral == 0.0 and integral2 == 0.0: - dict_bd_antibd[new_label] = { - "bonding": {"integral": integral2, "perc": 0.0}, - "antibonding": {"integral": integral, "perc": 0.0}, - } - else: - dict_bd_antibd[new_label] = { - "bonding": {"integral": integral2, "perc": perc2}, - "antibonding": {"integral": integral, "perc": perc}, - } - - return dict_bd_antibd - - def _integrate_antbdstates_below_efermi(self, cohp, start: float | None): - """ - Integrate the cohp data to compute bonding and anti-bonding contribution below efermi. - - .. warning:: NEEDS MORE TESTS - - This integrates the whole COHP curve that has been computed. - The energy range is very important. - At present the energy range considered is dependent on COHPstartEnergy - set during lobster runs. The bonding / antibonding integral values are sensitive to this parameter. - If COHPstartEnergy value does not cover entire range of VASP calculations then - absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values. - - :param cohp: cohp object - :param start: integration start energy in eV , eg start = -15 - - Returns: - absolute value of antibonding, percentage value of antibonding, - absolute value of bonding and percentage value of bonding interactions - """ - warnings.warn( - "The bonding, antibonding integral/percent values are numerical estimate." - " These values are sensitive to COHPstartEnergy parameter." - " If COHPstartEnergy value does not cover entire range of VASP calculations then" - " absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values.", - stacklevel=2, - ) - - from scipy.integrate import trapezoid - - def integrate_positive(y, x): - """ - Integrate only bonding interactions of COHPs. - - :param y: COHP values - :param x: Energy values - - Returns: - integrated value of bonding interactions - """ - y = np.asanyarray(y) - x = np.asanyarray(x) - - bonding = trapezoid(y, x) - - return np.round(bonding, 2) - - def integrate_negative(y, x): - """ - Integrate only anti-bonding interactions of COHPs. - - :param y: COHP values - :param x: Energy values - - Returns: - integrated value of anti-bonding interactions - """ - y = np.asanyarray(y) - x = np.asanyarray(x) - antibonding = trapezoid(y, x) - - return np.round(antibonding, 2) - - # will integrate spin.up and spin.down only below efermi - energies_corrected = cohp.energies - cohp.efermi - summedcohp = cohp.cohp[Spin.up] + cohp.cohp[Spin.down] if Spin.down in cohp.cohp else cohp.cohp[Spin.up] - - cohp_bf = [] - en_bf = [] - - for i, en in enumerate(energies_corrected): - if (start is None) and en <= 0: - en_bf.append(en) - cohp_bf.append(-1 * summedcohp[i]) - if (start is not None) and 0 >= en >= start: - en_bf.append(en) - cohp_bf.append(-1 * summedcohp[i]) - - # Separate the bonding and antibonding COHP values in separate lists - pos = [] - en_pos = [] - neg = [] - en_neg = [] - - for i, scohp in enumerate(cohp_bf): - if scohp >= 0: - pos.append(scohp) - en_pos.append(energies_corrected[i]) - - for i, scohp in enumerate(cohp_bf): - if scohp <= 0: - neg.append(-1 * scohp) - en_neg.append(energies_corrected[i]) - - antibonding = integrate_negative(y=neg, x=en_neg) - - bonding = integrate_positive(y=pos, x=en_pos) - - return ( - antibonding, - np.round(abs(antibonding) / (abs(bonding) + abs(antibonding)), 5), - bonding, - np.round(abs(bonding) / (abs(bonding) + abs(antibonding)), 5), - ) - - def _get_pop_type(self): - """ - Return the type of the input population file. - - Returns: - A String of analysed population can be COOP/COBI/COHP - """ - if self.are_cobis: - type_pop = "COBI" - elif self.are_coops: - type_pop = "COOP" - else: - type_pop = "COHP" - - return type_pop - - @staticmethod - def _get_bond_dict( - bond_strength_dict: dict, - small_antbd_dict: dict, - nameion: str | None = None, - large_antbd_dict: dict | None = None, - type_pop: str | None = None, - ): - """ - Return a bond_dict that contains information for each site. - - :param bond_strength_dict: dict with bond names as key and lists of bond strengths as items - :param small_antbd_dict: dict including if there are antibonding interactions, {'Yb-Sb': False} - :param nameion: name of the cation, e.g. Yb - :param large_antbd_dict: will be implemented later - :param type_pop: population type analyzed. eg. COHP - - Returns: - Eg., if type_pop == 'COHP', will return - dict including information on the anion (as label) and the ICOHPs in the item of the dict - ICOHP_mean refers to the mean ICOHP in eV - ICOHP_sum refers to the sum of the ICOHPs in eV - has_antibdg_states_below_Efermi is True if there are antibonding interactions below Efermi - "number_of_bonds" will count the numbers of bonds to the cation - - Example: - {'Sb': {'ICOHP_mean': '-1.85', 'ICOHP_sum': '-11.09', - 'has_antibdg_states_below_Efermi': False, 'number_of_bonds': 6}} - - """ - bond_dict = {} - - for key, item in bond_strength_dict.items(): - if nameion is not None: - a = key.split("-")[0] - b = key.split("-")[1] - if a == nameion: - key_here = b - elif b == nameion: - key_here = a - - if large_antbd_dict is None: - bond_dict[key_here] = { - f"I{type_pop}_mean": str(round(np.mean(item), 2)), - f"I{type_pop}_sum": str(round(np.sum(item), 2)), - "has_antibdg_states_below_Efermi": small_antbd_dict[key], - "number_of_bonds": len(item), - } - else: - bond_dict[key_here] = { - f"I{type_pop}_mean": str(round(np.mean(item), 2)), - f"I{type_pop}_sum": str(round(np.sum(item), 2)), - "has_antibdg_states_below_Efermi": small_antbd_dict[key], - "number_of_bonds": len(item), - "perc_antibdg_states_below_Efermi": large_antbd_dict[key], - } - - return bond_dict - - @staticmethod - def _get_orb_combinations(): - """ - Get a list of unique 2-length permutations and combinations. - - Generates 2-length permutations and combinations with replacement from the set ['s', 'p', 'd', 'f'] - - Returns: - list[tuple[str, str]]: A list of tuples, each containing two elements. - """ - list_orbs_comb: list[tuple[str, str]] = [] # type: ignore - - # Add all 2-length permutations - for perm in permutations(["s", "p", "d", "f"], 2): - list_orbs_comb.append(perm) # noqa : PERF402 - - # Add 2-length combinations with replacement, ensuring no duplicates - for comb in combinations_with_replacement(["s", "p", "d", "f"], 2): - if comb not in list_orbs_comb: - list_orbs_comb.append(comb) - - return list_orbs_comb - - def _group_orb_pairs(self, bond_label: str, orb_combinations: list[tuple[str, str]]) -> dict[str, list[str]]: - """ - Group orbital pairs based on the provided bond label. - - :param bond_label: The bond label to filter the orbitals. - :param orb_combinations: A list of tuples containing orbital combinations. - - Returns: - dict[str, List[str]]: A dictionary where the keys are top level orbital pairs and the values are lists of - sub orbitals associated to the bond label. - """ - orb_pair = {} # type: ignore - bond_label_int = int(bond_label) - 1 # convert label to int to access orbital data from icohpcollection - for sub_orb, data in self.chemenv.Icohpcollection._list_orb_icohp[bond_label_int].items(): - for orb1, orb2 in orb_combinations: - if orb1 in data["orbitals"][0][1].name and orb2 in data["orbitals"][1][1].name: - root_orb_pair = f"{data['orbitals'][0][0]}{orb1}-{data['orbitals'][1][0]}{orb2}" - if root_orb_pair not in orb_pair: - orb_pair[root_orb_pair] = [] - orb_pair[root_orb_pair].append(sub_orb) - return orb_pair - - @staticmethod - def _get_bond_orb_label_mapped_list(orb_pair: dict[str, list[str]], bond_labels: list[str], root_orb_pair: str): - """ - Get a lists of bond labels mapped to the corresponding orbital pair. - - :param orb_pair: A dictionary containing orbital pairs as keys and lists of - sub orbitals as values. - :param bond_labels: A list of bond labels. - :param root_orb_pair: The root key in orb_pair use to map bond labels list. - - Returns: - list: A list where the items of bond_labels are repeated based on the - length of orb_pair[root_orb_pair]. - """ - return [item for item in bond_labels for _ in range(len(orb_pair[root_orb_pair]))] - - def set_condensed_bonding_analysis(self): - """ - Condense the bonding analysis into a summary dictionary. - - Returns: - None - - """ - self.condensed_bonding_analysis = {} - # which icohps are considered - if self.which_bonds == "cation-anion": - limit_icohps = self.chemenv._get_limit_from_extremum( - self.chemenv.Icohpcollection, - self.cutoff_icohp, - adapt_extremum_to_add_cond=True, - additional_condition=1, - ) - elif self.which_bonds == "all": - limit_icohps = self.chemenv._get_limit_from_extremum( - self.chemenv.Icohpcollection, - self.cutoff_icohp, - adapt_extremum_to_add_cond=True, - additional_condition=0, - ) - # formula of the compound - formula = str(self.structure.composition.reduced_formula) - # set population type - type_pop = self._get_pop_type() - # how many inequivalent cations are in the structure - if self.which_bonds == "cation-anion": - number_considered_ions = len(self.seq_ineq_ions) - elif self.which_bonds == "all": - number_considered_ions = len(self.seq_ineq_ions) - - # what was the maximum bond lengths that was considered - max_bond_lengths = max(self.chemenv.Icohpcollection._list_length) - - # what are the charges for the cations in the structure - charge_list = self.chemenv.valences - - # dictionary including bonding information for each site - site_dict = {} - if self.which_bonds == "cation-anion": - for ication, ce, cation_anion_infos, labels, cohps in zip( - self.seq_ineq_ions, - self.seq_coord_ions, - self.seq_infos_bonds, - self.seq_labels_cohps, - self.seq_cohps, - ): - namecation = str(self.structure[ication].specie) - - # This will compute the mean strengths of ICOHPs - mean_icohps = self._get_strenghts_for_each_bond( - pairs=cation_anion_infos[4], - strengths=cation_anion_infos[1], - nameion=namecation, - ) - # pairs, strengths, nameion - # will collect if there are antibonding states present - antbdg = self._get_antibdg_states(cohps, labels, namecation) - dict_antibonding = self._integrate_antbdstates_below_efermi_for_set_cohps( - labels, cohps, nameion=namecation - ) - bond_dict = self._get_bond_dict(mean_icohps, antbdg, namecation, type_pop=type_pop) - bond_resolved_labels = self.get_site_bond_resolved_labels() - - for cation_name, icohp_data in bond_dict.items(): - for atom_pair, bonding_data in dict_antibonding.items(): - if namecation == atom_pair.split("-")[0] and cation_name == atom_pair.split("-")[1]: - icohp_data["bonding"] = bonding_data["bonding"] - icohp_data["antibonding"] = bonding_data["antibonding"] - if self.orbital_resolved: - # get orb resolved data to be added - orb_resolved_bond_info = self._get_orbital_resolved_data( - nameion=namecation, - iion=ication, - labels=labels, - bond_resolved_labels=bond_resolved_labels, - type_pop=type_pop, - ) - # match the dict key in bond_dict and get corresponding orbital data - for ion_atom_pair_orb in orb_resolved_bond_info: - orb_data_atom_pair = ion_atom_pair_orb.split(": ")[-1] - atom_pair_here = atom_pair.split("-") - atom_pair_here.sort() - if ( - orb_data_atom_pair == "-".join(atom_pair_here) - and (namecation + str(ication + 1) + ":") in ion_atom_pair_orb - ): - icohp_data["orbital_data"] = orb_resolved_bond_info[ion_atom_pair_orb] - - orb_data_stats = self._get_bond_resolved_data_stats( - orb_resolved_bond_data=orb_resolved_bond_info[ion_atom_pair_orb], - ) - - icohp_data["orbital_data"].update(orb_data_stats) - - site_dict[ication] = { - "env": ce, - "bonds": bond_dict, - "ion": namecation, - "charge": charge_list[ication], - "relevant_bonds": cation_anion_infos[3], - } - elif self.which_bonds == "all": - for iion, ce, bond_infos, labels, cohps in zip( - self.seq_ineq_ions, - self.seq_coord_ions, - self.seq_infos_bonds, - self.seq_labels_cohps, - self.seq_cohps, - ): - nameion = str(self.structure[iion].specie) - - # This will compute the mean strengths of ICOHPs - mean_icohps = self._get_strenghts_for_each_bond( - pairs=bond_infos[4], strengths=bond_infos[1], nameion=None - ) - # pairs, strengths, nameion - # will collect if there are antibonding states present - antbdg = self._get_antibdg_states(cohps, labels, nameion=None) - - dict_antibonding = self._integrate_antbdstates_below_efermi_for_set_cohps(labels, cohps, nameion) - - bond_dict = self._get_bond_dict(mean_icohps, antbdg, nameion=nameion, type_pop=type_pop) - bond_resolved_labels = self.get_site_bond_resolved_labels() - - for cation_name, icohp_data in bond_dict.items(): - for atom_pair, bonding_data in dict_antibonding.items(): - if nameion == atom_pair.split("-")[0] and cation_name == atom_pair.split("-")[1]: - icohp_data["bonding"] = bonding_data["bonding"] - icohp_data["antibonding"] = bonding_data["antibonding"] - if self.orbital_resolved: - # get orb resolved data to be added - orb_resolved_bond_info = self._get_orbital_resolved_data( - nameion=nameion, - iion=iion, - labels=labels, - bond_resolved_labels=bond_resolved_labels, - type_pop=type_pop, - ) - # match the dict key in bond_dict and get corresponding orbital data - for ion_atom_pair_orb in orb_resolved_bond_info: - orb_data_atom_pair = ion_atom_pair_orb.split(": ")[-1] - atom_pair_here = atom_pair.split("-") - atom_pair_here.sort() - if ( - orb_data_atom_pair == "-".join(atom_pair_here) - and (nameion + str(iion + 1) + ":") in ion_atom_pair_orb - ): - icohp_data["orbital_data"] = orb_resolved_bond_info[ion_atom_pair_orb] - - orb_data_stats = self._get_bond_resolved_data_stats( - orb_resolved_bond_data=orb_resolved_bond_info[ion_atom_pair_orb], - ) - - icohp_data["orbital_data"].update(orb_data_stats) - - site_dict[iion] = { - "env": ce, - "bonds": bond_dict, - "ion": nameion, - "charge": charge_list[iion], - "relevant_bonds": bond_infos[3], - } - - if self.path_to_madelung is None and self.madelung_obj is None: - if self.which_bonds == "cation-anion": - # This sets the dictionary including the most important information on the compound - self.condensed_bonding_analysis = { - "formula": formula, - "max_considered_bond_length": max_bond_lengths, - f"limit_i{type_pop.lower()}": limit_icohps, - "number_of_considered_ions": number_considered_ions, - "sites": site_dict, - "type_charges": self.type_charge, - } - elif self.which_bonds == "all": - self.condensed_bonding_analysis = { - "formula": formula, - "max_considered_bond_length": max_bond_lengths, - f"limit_i{type_pop.lower()}": limit_icohps, - "number_of_considered_ions": number_considered_ions, - "sites": site_dict, - "type_charges": self.type_charge, - } - else: - madelung = MadelungEnergies(self.path_to_madelung) if self.path_to_madelung else self.madelung_obj - if self.type_charge == "Mulliken": - madelung_energy = madelung.madelungenergies_mulliken - elif self.type_charge == "Löwdin": - madelung_energy = madelung.madelungenergies_loewdin - else: - madelung_energy = None - # This sets the dictionary including the most important information on the compound - if self.which_bonds == "cation-anion": - self.condensed_bonding_analysis = { - "formula": formula, - "max_considered_bond_length": max_bond_lengths, - f"limit_i{type_pop.lower()}": limit_icohps, - "number_of_considered_ions": number_considered_ions, - "sites": site_dict, - "type_charges": self.type_charge, - "madelung_energy": madelung_energy, - } - elif self.which_bonds == "all": - self.condensed_bonding_analysis = { - "formula": formula, - "max_considered_bond_length": max_bond_lengths, - f"limit_i{type_pop.lower()}": limit_icohps, - "number_of_considered_ions": number_considered_ions, - "sites": site_dict, - "type_charges": self.type_charge, - "madelung_energy": madelung_energy, - } - - def set_summary_dicts(self): - """ - Set summary dict that can be used for correlations. - - bond_dict that includes information on each bond - - "has_antbd" tells if there are antbonding states - "ICOHP_mean" shows the mean of all ICOHPs in EV - - {'Yb-Sb': { 'has_antbdg': False, 'ICOHP_mean': -1.7448}, - 'Mn-Sb': { 'has_antbdg': True, 'ICOHP_mean': -1.525}} - - a cation dict that includes all different coordination environments and counts for them - {'Na': {'T:4': 4, 'A:2': 4}, 'Si': {'T:6': 4, 'PP:6': 4}} - - Returns: - None - - """ - relevant_ion_ids = [isite for isite in self.list_equivalent_sites if isite in self.seq_ineq_ions] - # set population type - type_pop = self._get_pop_type() - - final_dict_bonds = {} - for key in relevant_ion_ids: - item = self.condensed_bonding_analysis["sites"][key] - for type, properties in item["bonds"].items(): - label_list = [item["ion"], str(type)] - new_label = sorted(label_list.copy()) - label = str(new_label[0]) + "-" + str(new_label[1]) - - if label not in final_dict_bonds: - final_dict_bonds[label] = { - "number_of_bonds": int(properties["number_of_bonds"]), - f"I{type_pop}_sum": float(properties[f"I{type_pop}_sum"]), - "has_antbdg": properties["has_antibdg_states_below_Efermi"], - } - else: - final_dict_bonds[label]["number_of_bonds"] += int(properties["number_of_bonds"]) - final_dict_bonds[label][f"I{type_pop}_sum"] += float(properties[f"I{type_pop}_sum"]) - final_dict_bonds[label]["has_antbdg"] = ( - final_dict_bonds[label]["has_antbdg"] or properties["has_antibdg_states_below_Efermi"] - ) - self.final_dict_bonds = {} - for key, item in final_dict_bonds.items(): - self.final_dict_bonds[key] = {} - self.final_dict_bonds[key][f"I{type_pop}_mean"] = item[f"I{type_pop}_sum"] / (item["number_of_bonds"]) - self.final_dict_bonds[key]["has_antbdg"] = item["has_antbdg"] - - # rework, add all environments! - final_dict_ions = {} - for key in relevant_ion_ids: - if self.condensed_bonding_analysis["sites"][key]["ion"] not in final_dict_ions: - final_dict_ions[self.condensed_bonding_analysis["sites"][key]["ion"]] = [ - self.condensed_bonding_analysis["sites"][key]["env"] - ] - else: - final_dict_ions[self.condensed_bonding_analysis["sites"][key]["ion"]].append( - self.condensed_bonding_analysis["sites"][key]["env"] - ) - - self.final_dict_ions = {} - for key, item in final_dict_ions.items(): - self.final_dict_ions[key] = dict(Counter(item)) - - @staticmethod - def get_lobster_calc_quality_summary( - path_to_poscar: str | None = None, - path_to_lobsterout: str | None = None, - path_to_lobsterin: str | None = None, - path_to_potcar: str | None = None, - potcar_symbols: list | None = None, - path_to_charge: str | None = None, - path_to_bandoverlaps: str | None = None, - path_to_doscar: str | None = None, - path_to_vasprun: str | None = None, - structure_obj: Structure | None = None, - lobsterin_obj: Lobsterin | None = None, - lobsterout_obj: Lobsterout | None = None, - charge_obj: Charge | None = None, - bandoverlaps_obj: Bandoverlaps | None = None, - lobster_completedos_obj: LobsterCompleteDos | None = None, - vasprun_obj: Vasprun | None = None, - dos_comparison: bool = False, - e_range: list = [-5, 0], - n_bins: int | None = None, - bva_comp: bool = False, - ) -> dict: - """ - Analyze LOBSTER calculation quality. - - :param path_to_poscar: path to structure file - :param path_to_lobsterout: path to lobsterout file - :param path_to_lobsterin: path to lobsterin file - :param path_to_potcar: path to VASP potcar file - :param potcar_symbols: list of potcar symbols from potcar file (can be used if no potcar available) - :param path_to_charge: path to CHARGE.lobster file - :param path_to_bandoverlaps: path to bandOverlaps.lobster file - :param path_to_doscar: path to DOSCAR.lobster or DOSCAR.LSO.lobster file - :param path_to_vasprun: path to vasprun.xml file - :param structure_obj: pymatgen pymatgen.core.structure.Structure object - :param lobsterin_obj: pymatgen.lobster.io.Lobsterin object - :param lobsterout_obj: pymatgen lobster.io.Lobsterout object - :param charge_obj: pymatgen lobster.io.Charge object - :param bandoverlaps_obj: pymatgen lobster.io.BandOverlaps object - :param lobster_completedos_obj: pymatgen.electronic_structure.dos.LobsterCompleteDos object - :param vasprun_obj: pymatgen vasp.io.Vasprun object - :param dos_comparison: will compare DOS from VASP and LOBSTER and return tanimoto index - :param e_range: energy range for DOS comparisons - :param n_bins: number of bins to discretize DOS for comparisons - :param bva_comp: Compares LOBSTER charge signs with Bond valence charge signs - - Returns: - A dict of summary of LOBSTER calculation quality by analyzing basis set used, - charge spilling from lobsterout/ PDOS comparisons of VASP and LOBSTER / - BVA charge comparisons - - """ - quality_dict = {} - - if path_to_potcar and not potcar_symbols and not path_to_vasprun and not vasprun_obj: - potcar_names = Lobsterin._get_potcar_symbols(POTCAR_input=path_to_potcar) - elif not path_to_potcar and not path_to_vasprun and not vasprun_obj and potcar_symbols: - potcar_names = potcar_symbols - elif path_to_vasprun and not vasprun_obj: - vasprun = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False) - potcar_names = [potcar.split(" ")[1] for potcar in vasprun.potcar_symbols] - elif vasprun_obj and not path_to_vasprun: - potcar_names = [potcar.split(" ")[1] for potcar in vasprun_obj.potcar_symbols] - else: - raise ValueError( - "Please provide either path_to_potcar or list of " - "potcar_symbols or path to vasprun.xml or vasprun object. " - "Crucial to identify basis used for projections" - ) - - if path_to_poscar: - struct = Structure.from_file(path_to_poscar) - elif structure_obj: - struct = structure_obj - else: - raise ValueError("Please provide path_to_poscar or structure_obj") - - ref_bases = Lobsterin.get_all_possible_basis_functions(structure=struct, potcar_symbols=potcar_names) - - if path_to_lobsterin: - lobs_in = Lobsterin.from_file(path_to_lobsterin) - elif lobsterin_obj: - lobs_in = lobsterin_obj - else: - raise ValueError("Please provide path_to_lobsterin or lobsterin_obj") - - calc_basis = [] - for basis in lobs_in["basisfunctions"]: - basis_sep = basis.split()[1:] - basis_comb = " ".join(basis_sep) - calc_basis.append(basis_comb) - - if calc_basis == list(ref_bases[0].values()): - quality_dict["minimal_basis"] = True # type: ignore - else: - quality_dict["minimal_basis"] = False # type: ignore - warnings.warn( - "Consider rerunning the calc with the minimum basis as well. Choosing is " - "larger basis set is recommended if you see a significant improvement of " - "the charge spilling and material has non-zero band gap.", - stacklevel=2, - ) - - if path_to_lobsterout: - lob_out = Lobsterout(path_to_lobsterout) - elif lobsterout_obj: - lob_out = lobsterout_obj - else: - raise ValueError("Please provide path_to_lobsterout or lobsterout_obj") - - quality_dict["charge_spilling"] = { - "abs_charge_spilling": round((sum(lob_out.charge_spilling) / 2) * 100, 4), - "abs_total_spilling": round((sum(lob_out.total_spilling) / 2) * 100, 4), - } # type: ignore - - if path_to_bandoverlaps is not None and not bandoverlaps_obj: - band_overlaps = Bandoverlaps(filename=path_to_bandoverlaps) if Path(path_to_bandoverlaps).exists() else None - elif path_to_bandoverlaps is None and bandoverlaps_obj: - band_overlaps = bandoverlaps_obj - else: - band_overlaps = None - - if band_overlaps is not None: - for line in lob_out.warning_lines: - if "k-points could not be orthonormalized" in line: - total_kpoints = int(line.split(" ")[2]) - - # store actual number of devations above pymatgen default limit of 0.1 - dev_val = [] - for dev in band_overlaps.max_deviation: - if dev > 0.1: - dev_val.append(dev) - - quality_dict["band_overlaps_analysis"] = { # type: ignore - "file_exists": True, - "limit_maxDeviation": 0.1, - "has_good_quality_maxDeviation": band_overlaps.has_good_quality_maxDeviation(limit_maxDeviation=0.1), - "max_deviation": round(max(band_overlaps.max_deviation), 4), - "percent_kpoints_abv_limit": round((len(dev_val) / total_kpoints) * 100, 4), - } - - else: - quality_dict["band_overlaps_analysis"] = { # type: ignore - "file_exists": False, - "limit_maxDeviation": None, - "has_good_quality_maxDeviation": True, - "max_deviation": None, - "percent_kpoints_abv_limit": None, - } - - if bva_comp: - try: - bond_valence = BVAnalyzer() - - bva_oxi = [] - if path_to_charge and not charge_obj: - lobs_charge = Charge(filename=path_to_charge) - elif not path_to_charge and charge_obj: - lobs_charge = charge_obj - else: - raise Exception("BVA comparison is requested, thus please provide path_to_charge or charge_obj") - for i in bond_valence.get_valences(structure=struct): - if i >= 0: - bva_oxi.append("POS") - else: - bva_oxi.append("NEG") - - mull_oxi = [] - for i in lobs_charge.Mulliken: - if i >= 0: - mull_oxi.append("POS") - else: - mull_oxi.append("NEG") - - loew_oxi = [] - for i in lobs_charge.Loewdin: - if i >= 0: - loew_oxi.append("POS") - else: - loew_oxi.append("NEG") - - quality_dict["charge_comparisons"] = {} # type: ignore - if mull_oxi == bva_oxi: - quality_dict["charge_comparisons"]["bva_mulliken_agree"] = True # type: ignore - else: - quality_dict["charge_comparisons"]["bva_mulliken_agree"] = False # type: ignore - - if mull_oxi == bva_oxi: - quality_dict["charge_comparisons"]["bva_loewdin_agree"] = True # type: ignore - else: - quality_dict["charge_comparisons"]["bva_loewdin_agree"] = False # type: ignore - - except ValueError: - quality_dict["charge_comparisons"] = {} # type: ignore - warnings.warn( - "Oxidation states from BVA analyzer cannot be determined. " - "Thus BVA charge comparison will be skipped", - stacklevel=2, - ) - if dos_comparison: - if "LSO" not in str(path_to_doscar).split("."): - warnings.warn( - "Consider using DOSCAR.LSO.lobster, as non LSO DOS from LOBSTER can have negative DOS values", - stacklevel=2, - ) - if path_to_doscar: - doscar_lobster = Doscar( - doscar=path_to_doscar, - structure_file=path_to_poscar, - structure=structure_obj, - ) - - dos_lobster = doscar_lobster.completedos - elif lobster_completedos_obj: - dos_lobster = lobster_completedos_obj - else: - raise ValueError( - "Dos comparison is requested, so please provide either path_to_doscar or lobster_completedos_obj" - ) - - if path_to_vasprun: - vasprun = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False) - elif vasprun_obj: - vasprun = vasprun_obj - else: - raise ValueError( - "Dos comparison is requested, so please provide either path to vasprun.xml or vasprun_obj" - ) - dos_vasp = vasprun.complete_dos - - quality_dict["dos_comparisons"] = {} # type: ignore - - for orb in dos_lobster.get_spd_dos(): - if e_range[0] >= min(dos_vasp.energies) and e_range[0] >= min(dos_lobster.energies): - min_e = e_range[0] - else: - warnings.warn( - "Minimum energy range requested for DOS comparisons is not available " - "in VASP or LOBSTER calculation. Thus, setting min_e to -5 eV", - stacklevel=2, - ) - min_e = -5 - - if e_range[-1] <= max(dos_vasp.energies) and e_range[-1] <= max(dos_lobster.energies): - max_e = e_range[-1] - else: - warnings.warn( - "Maximum energy range requested for DOS comparisons is not available " - "in VASP or LOBSTER calculation. Thus, setting max_e to 0 eV", - stacklevel=2, - ) - max_e = 0 - - if np.diff(dos_vasp.energies)[0] >= 0.1 or np.diff(dos_lobster.energies)[0] >= 0.1: - warnings.warn( - "Input DOS files have very few points in the energy interval and thus " - "comparisons will not be reliable. Please rerun the calculations with " - "higher number of DOS points. Set NEDOS and COHPSteps tags to >= 2000 in VASP and LOBSTER " - "calculations, respectively.", - stacklevel=2, - ) - - if not n_bins: - n_bins = 56 - - fp_lobster_orb = dos_lobster.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, - type=orb.name, - ) - fp_vasp_orb = dos_vasp.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, - type=orb.name, - ) - - tani_orb = round( - dos_vasp.get_dos_fp_similarity(fp_lobster_orb, fp_vasp_orb, tanimoto=True), - 4, - ) - quality_dict["dos_comparisons"][f"tanimoto_orb_{orb.name}"] = tani_orb # type: ignore - - fp_lobster = dos_lobster.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, - type="summed_pdos", - ) - fp_vasp = dos_vasp.get_dos_fp( - min_e=min_e, - max_e=max_e, - n_bins=n_bins, - normalize=True, - type="summed_pdos", - ) - - tanimoto_summed = round(dos_vasp.get_dos_fp_similarity(fp_lobster, fp_vasp, tanimoto=True), 4) - quality_dict["dos_comparisons"]["tanimoto_summed"] = tanimoto_summed # type: ignore - quality_dict["dos_comparisons"]["e_range"] = [min_e, max_e] # type: ignore - quality_dict["dos_comparisons"]["n_bins"] = n_bins # type: ignore - - return quality_dict +# Copyright (c) lobsterpy development team +# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License + +"""This module defines classes to analyze the COHPs automatically.""" + +from __future__ import annotations + +import warnings +from collections import Counter +from itertools import combinations_with_replacement, permutations +from pathlib import Path + +import numpy as np +from pymatgen.analysis.bond_valence import BVAnalyzer +from pymatgen.core.structure import Structure +from pymatgen.electronic_structure.cohp import CompleteCohp +from pymatgen.electronic_structure.core import Spin +from pymatgen.electronic_structure.dos import LobsterCompleteDos +from pymatgen.io.lobster import ( + Bandoverlaps, + Charge, + Doscar, + Icohplist, + Lobsterin, + Lobsterout, + MadelungEnergies, +) +from pymatgen.io.lobster.lobsterenv import LobsterNeighbors +from pymatgen.io.vasp.outputs import Vasprun +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + + +class Analysis: + """ + Class to automatically analyze COHP/COOP/COBI populations from Lobster. + + :param are_cobis: bool indicating if file contains COBI/ICOBI data + :param are_coops: bool indicating if file contains COOP/ICOOP data + :param cutoff_icohp: Cutoff in percentage for evaluating neighbors based on ICOHP values. + cutoff_icohp*max_icohp limits the number of considered neighbours for evaluating environments. + :param path_to_cohpcar: path to `COHPCAR.lobster` or `COBICAR.lobster` or `COOPCAR.lobster` . + :param path_to_charge: path to `CHARGE.lobster`. + :param path_to_icohplist: path to `ICOHPLIST.lobster` or `ICOBILIST.lobster` or `ICOOPLIST.lobster`. + :param path_to_poscar: path to structure (e.g., `POSCAR` or `POSCAR.lobster`) + :param path_to_madelung: path to `MadelungEnergies.lobster`. + :param charge_obj: pymatgen lobster.io.charge object + :param completecohp_obj: pymatgen.electronic_structure.cohp.CompleteCohp object + :param icohplist_obj: pymatgen lobster.io.Icohplist object + :param madelung_obj: pymatgen lobster.io.MadelungEnergies object + :param noise_cutoff: Sets the lower limit tolerance for ICOHPs or ICOOPs or ICOBIs considered + in analysis. + :param orbital_cutoff: Sets the minimum percentage for the orbital contribution considered to be + relevant in orbital resolved analysis. (Affects only when orbital_resolved argument is set to True) + Set it to 0 to get results of all orbitals in the detected relevant bonds. Default is to 0.05 i.e. + only analyzes if orbital contribution is 5 % or more. + :param orbital_resolved: bool indicating whether orbital wise analysis is performed + :param type_charge: If no path_to_charge is given, Valences will be used. Otherwise, Mulliken charges. + Löwdin charges cannot be selected at the moment. + :param which_bonds: Selects kinds of bonds that are analyzed. `cation-anion` is the default. + Alternatively, `all` bonds can also be selected. Support to other kinds of bonds will be + added soon. + :param summed_spins: if True, COHP `Spin.up` and `Spin.down` populations will be summed + :param start: sets the lower limit of energy for evaluation of bonding and antibonding + percentages below efermi. Defaults to None (i.e., all populations below efermi are included) + + Attributes: + - condensed_bonding_analysis: dict including a summary of the most important bonding properties + - final_dict_bonds: dict including information on ICOHPs per bond type + - final_dict_ions: dict including information on environments of cations + - chemenv: pymatgen.io.lobster.lobsterenv.LobsterNeighbors object + - lse: LightStructureEnvironment from pymatgen + - anion_types: Set of Element objects from pymatgen + - list_equivalent_sites: list of site indices of sites that indicate which sites are equivalent + e.g., [0 1 2 2 2] where site 0, 1, 2 indicate sites that are independent from each other + - seq_cohps: list of cohps + - seq_coord_ions: list of co-ordination environment strings for each cation + - seq_equivalent_sites: seq of inequivalent sites + - seq_ineq_ions: seq of inequivalent cations/sites in the structure + - seq_infos_bonds (list): information on cation anion bonds (lists + of pymatgen.io.lobster.lobsterenv.ICOHPNeighborsInfo) + - spg: space group information + - structure: Structure object + + """ + + def __init__( + self, + path_to_poscar: str | Path | None, + path_to_icohplist: str | Path | None, + path_to_cohpcar: str | Path | None, + path_to_charge: str | Path | None = None, + path_to_madelung: str | Path | None = None, + icohplist_obj: Icohplist | None = None, + completecohp_obj: CompleteCohp | None = None, + charge_obj: Charge | None = None, + madelung_obj: MadelungEnergies | None = None, + are_cobis: bool = False, + are_coops: bool = False, + cutoff_icohp: float = 0.1, + noise_cutoff: float = 0.1, + orbital_cutoff: float = 0.05, + orbital_resolved: bool = False, + start: float | None = None, + summed_spins: bool = True, + type_charge: str | None = None, + which_bonds: str = "cation-anion", + ): + """ + Initialize automatic bonding analysis. + + :param are_cobis: bool indicating if file contains COBI/ICOBI data + :param are_coops: bool indicating if file contains COOP/ICOOP data + :param cutoff_icohp: Cutoff in percentage for evaluating neighbors based on ICOHP values. + cutoff_icohp*max_icohp limits the number of considered neighbours for evaluating environments. + :param path_to_cohpcar: path to `COHPCAR.lobster` or `COBICAR.lobster` or `COOPCAR.lobster` . + :param path_to_charge: path to `CHARGE.lobster`. + :param path_to_icohplist: path to `ICOHPLIST.lobster` or `ICOBILIST.lobster` or `ICOOPLIST.lobster`. + :param path_to_poscar: path to structure (e.g., `POSCAR` or `POSCAR.lobster`) + :param path_to_madelung: path to `MadelungEnergies.lobster`. + :param charge_obj: pymatgen lobster.io.charge object (Optional) + :param completecohp_obj: pymatgen.electronic_structure.cohp.CompleteCohp object + :param icohplist_obj: pymatgen lobster.io.Icohplist object + :param madelung_obj: pymatgen lobster.io.MadelungEnergies object + :param noise_cutoff: Sets the lower limit tolerance for ICOHPs or ICOOPs or ICOBIs considered + in analysis. + :param orbital_cutoff: Sets the minimum percentage for the orbital contribution considered to be + relevant in orbital resolved analysis. (Affects only when orbital_resolved argument is set to True) + Set it to 0 to get results of all orbitals in the detected relevant bonds. Default is to 0.05 i.e. + only analyzes if orbital contribution is 5 % or more. + :param orbital_resolved: bool indicating whether orbital wise analysis is performed + :param type_charge: If no path_to_charge is given, Valences will be used. Otherwise, Mulliken charges. + Löwdin charges cannot be selected at the moment. + :param which_bonds: Selects kinds of bonds that are analyzed. `cation-anion` is the default. + Alternatively, `all` bonds can also be selected. Support to other kinds of bonds will be + added soon. + :param summed_spins: if True, COHP `Spin.up` and `Spin.down` populations will be summed + :param start: sets the lower limit of energy for evaluation of bonding and antibonding + percentages below efermi. Defaults to None (i.e., all populations below efermi are included) + + """ + self.start = start + self.completecohp_obj = completecohp_obj + self.icohplist_obj = icohplist_obj + # checks to ensure LobsterEnv inputs are not duplicated in case users provide both path and obj + if self.completecohp_obj is not None and self.icohplist_obj is not None: + self.path_to_poscar = None + self.path_to_cohpcar = None + self.path_to_icohplist = None + else: + self.path_to_poscar = path_to_poscar + self.path_to_icohplist = path_to_icohplist + self.path_to_cohpcar = path_to_cohpcar + self.which_bonds = which_bonds + self.cutoff_icohp = cutoff_icohp + self.orbital_cutoff = orbital_cutoff + self.path_to_charge = path_to_charge + self.charge_obj = charge_obj + self.path_to_madelung = path_to_madelung + self.madelung_obj = madelung_obj + self.are_cobis = are_cobis + self.are_coops = are_coops + self.noise_cutoff = noise_cutoff + self.setup_env() + self.get_information_all_bonds(summed_spins=summed_spins) + self.orbital_resolved = orbital_resolved + + # This determines how cations and anions + if path_to_charge is None and charge_obj is None: + self.type_charge = "Valences" + else: + if type_charge is None: + self.type_charge = "Mulliken" + elif type_charge == "Mulliken": + self.type_charge = "Mulliken" + elif type_charge == "Löwdin": + raise ValueError("Only Mulliken charges can be used here at the moment. Implementation will follow.") + else: + self.type_charge = "Valences" + print("type_charge cannot be read! Please use Mulliken/Löwdin. Now, we will use valences") + + self.set_condensed_bonding_analysis() + self.set_summary_dicts() + self.path_to_madelung = path_to_madelung + + def setup_env(self): + """ + Set up the light structure environments based on COHPs using this method. + + Returns: + None + + """ + self.structure = ( + Structure.from_file(self.path_to_poscar) if self.path_to_poscar else self.completecohp_obj.structure + ) + sga = SpacegroupAnalyzer(structure=self.structure) + symmetry_dataset = sga.get_symmetry_dataset() + equivalent_sites = symmetry_dataset["equivalent_atoms"] + self.list_equivalent_sites = equivalent_sites + self.seq_equivalent_sites = list(set(equivalent_sites)) + self.spg = symmetry_dataset["international"] + + if self.which_bonds == "cation-anion": + try: + self.chemenv = LobsterNeighbors( + filename_icohp=self.path_to_icohplist, + obj_icohp=self.icohplist_obj, + structure=self.structure, + additional_condition=1, + perc_strength_icohp=self.cutoff_icohp, + filename_charge=self.path_to_charge, + obj_charge=self.charge_obj, + valences=None, + valences_from_charges=True, + adapt_extremum_to_add_cond=True, + are_cobis=self.are_cobis, + are_coops=self.are_coops, + noise_cutoff=self.noise_cutoff, + ) + except ValueError as err: + if ( + str(err) == "min() arg is an empty sequence" + or str(err) == "All valences are equal to 0, additional_conditions 1, 3, 5 and 6 will not work" + ): + raise ValueError( + "Consider switching to an analysis of all bonds and not only cation-anion bonds." + " It looks like no cations are detected." + ) + raise err + elif self.which_bonds == "all": + # raise ValueError("only cation anion bonds implemented so far") + self.chemenv = LobsterNeighbors( + filename_icohp=self.path_to_icohplist, + obj_icohp=self.icohplist_obj, + structure=self.structure, + additional_condition=0, + perc_strength_icohp=self.cutoff_icohp, + filename_charge=self.path_to_charge, + obj_charge=self.charge_obj, + valences_from_charges=True, + adapt_extremum_to_add_cond=True, + are_cobis=self.are_cobis, + are_coops=self.are_coops, + noise_cutoff=self.noise_cutoff, + ) + + else: + raise ValueError("only cation anion and all bonds implemented so far") + + # determine cations and anions + try: + if self.which_bonds == "cation-anion": + self.lse = self.chemenv.get_light_structure_environment(only_cation_environments=True) + elif self.which_bonds == "all": + self.lse = self.chemenv.get_light_structure_environment(only_cation_environments=False) + except ValueError: + + class Lse: + """Test class when error was raised.""" + + def __init__(self, chemenv, valences=None): + """ + Test class when error was raised. + + :param chemenv: LobsterNeighbors object + :param valences: list of valences + + """ + if valences is None: + self.coordination_environments = [] + for coord in chemenv: + if len(coord) > 0: + self.coordination_environments.append([{"ce_symbol": str(len(coord))}]) + else: + self.coordination_environments.append([{"ce_symbol": None}]) + else: + self.coordination_environments = [] + + for val, coord in zip(valences, chemenv): + if val >= 0.0 and len(coord) > 0: + self.coordination_environments.append([{"ce_symbol": str(len(coord))}]) + else: + self.coordination_environments.append([{"ce_symbol": None}]) + + if self.which_bonds == "all": + self.lse = Lse(self.chemenv.list_coords) + elif self.which_bonds == "cation-anion": + # make a new list + self.lse = Lse(self.chemenv.list_coords, self.chemenv.valences) + + def get_information_all_bonds(self, summed_spins: bool = True): + """ + Gather all information on the bonds within the compound with this method. + + Returns: + None + + """ + if self.which_bonds == "cation-anion": + # this will only analyze cation anion bonds which simplifies the analysis + self.seq_ineq_ions = [] + self.seq_coord_ions = [] + self.seq_infos_bonds = [] + self.seq_labels_cohps = [] + self.seq_cohps = [] + # only_bonds_to + + self.anion_types = self.chemenv.anion_types + for ice, ce in enumerate(self.lse.coordination_environments): + # only look at inequivalent sites (use of symmetry to speed everything up!)! + # only look at those cations that have cation-anion bonds + if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None: + self.seq_ineq_ions.append(ice) + + self.seq_coord_ions.append(ce[0]["ce_symbol"]) + self.seq_infos_bonds.append(self.chemenv.get_info_icohps_to_neighbors([ice])) + + aniontype_labels = [] + aniontype_cohps = [] + + # go through all anions in the structure! + for anion in self.anion_types: + # get labels and summed cohp objects + labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors( + path_to_cohpcar=self.path_to_cohpcar, + obj_cohpcar=self.completecohp_obj, + isites=[ice], + summed_spin_channels=summed_spins, + per_bond=False, + only_bonds_to=[str(anion)], + ) + + aniontype_labels.append(labels) + aniontype_cohps.append(summedcohps) + + self.seq_labels_cohps.append(aniontype_labels) + self.seq_cohps.append(aniontype_cohps) + + elif self.which_bonds == "all": + # this will only analyze all bonds + + self.seq_ineq_ions = [] + self.seq_coord_ions = [] + self.seq_infos_bonds = [] + self.seq_labels_cohps = [] + self.seq_cohps = [] + # only_bonds_to + self.elements = self.structure.composition.elements + # self.anion_types = self.chemenv.get_anion_types() + for ice, ce in enumerate(self.lse.coordination_environments): + # only look at inequivalent sites (use of symmetry to speed everything up!)! + # only look at those cations that have cation-anion bonds + if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None: + self.seq_ineq_ions.append(ice) + self.seq_coord_ions.append(ce[0]["ce_symbol"]) + self.seq_infos_bonds.append(self.chemenv.get_info_icohps_to_neighbors([ice])) + + type_labels = [] + type_cohps = [] + + for element in self.elements: + # get labels and summed cohp objects + labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors( + path_to_cohpcar=self.path_to_cohpcar, + obj_cohpcar=self.completecohp_obj, + isites=[ice], + onlycation_isites=False, + summed_spin_channels=summed_spins, + per_bond=False, + only_bonds_to=[str(element)], + ) + + type_labels.append(labels) + type_cohps.append(summedcohps) + + self.seq_labels_cohps.append(type_labels) + self.seq_cohps.append(type_cohps) + + def get_site_bond_resolved_labels(self): + """ + Return relevant bond labels for each symmetrically independent site. + + Returns: + dict with bond labels for each site, e.g. + {'Na1: Na-Cl': ['21', '23', '24', '27', '28', '30']} + + """ + bonds = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore + labels = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore + for inx, bond_info in enumerate(self.seq_infos_bonds): + for ixx, val in enumerate(bond_info.atoms): + label_srt = sorted(val.copy()) + bonds[inx].append( + self.structure.sites[bond_info.central_isites[0]].species_string + + str(bond_info.central_isites[0] + 1) + + ": " + + label_srt[0].strip("0123456789") + + "-" + + label_srt[1].strip("0123456789") + ) + labels[inx].append(bond_info.labels[ixx]) + + label_data = {} + for indx, atom_pairs in enumerate(bonds): + searched_atom_pairs = set(atom_pairs) + for search_item in searched_atom_pairs: + indices = [i for i, pair in enumerate(atom_pairs) if pair == search_item] + filtered_bond_label_list = [labels[indx][i] for i in indices] + label_data.update({search_item: filtered_bond_label_list}) + + return label_data + + def _get_orbital_resolved_data( + self, + nameion: str, + iion: int, + labels: list[str], + bond_resolved_labels: dict[str, list[str]], + type_pop: str, + ): + """ + Retrieve orbital-wise analysis data. + + :param nameion: name of symmetrically relevant cation or anion + :param iion: index of symmetrically relevant cation or anion + :param labels: list of bond label names + :param bond_resolved_labels: dict of bond labels from ICOHPLIST resolved for each bond + :param type_pop: population type analyzed. e.g. COHP or COOP or COBI + + Returns: + dict consisting of relevant orbitals (contribution > 5 % to overall ICOHP or ICOBI or ICOOP), + bonding and antibonding percentages with bond label names as keys. + """ + orb_resolved_bond_info = {} + for label in labels: + if label is not None: + bond_resolved_label_key = nameion + str(iion + 1) + ":" + label.split("x")[-1] + bond_labels = bond_resolved_labels[bond_resolved_label_key] + orb_combinations = self._get_orb_combinations() + grouped_orb_pairs = self._group_orb_pairs(bond_label=bond_labels[0], orb_combinations=orb_combinations) + # available_orbitals = list(self.chemenv.completecohp.orb_res_cohp[bond_labels[0]].keys()) + # initialize empty list to store orb paris for bonding, + # antibonding integrals and percentages + bndg_orb_pair_list = [] + bndg_orb_integral_list = [] + bndg_orb_perc_list = [] + bndg_orb_icohp_list = [] + antibndg_orb_integral_list = [] + antibndg_orb_perc_list = [] + antibndg_orb_pair_list = [] + antibndg_orb_icohp_list = [] + + # get total summed cohps using label list + cohp_summed = self.chemenv.completecohp.get_summed_cohp_by_label_list(label_list=bond_labels) + if type_pop.lower() == "cohp": + ( + antibndg_tot, + per_anti_tot, + bndg_tot, + per_bndg_tot, + ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed, start=self.start) + else: + ( + bndg_tot, + per_bndg_tot, + antibndg_tot, + per_anti_tot, + ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed, start=self.start) + + orb_bonding_dict_data = {} # type: ignore + # For each orbital collect the contributions of summed bonding + # and antibonding interactions separately + for orb in grouped_orb_pairs: + mapped_bond_labels = self._get_bond_orb_label_mapped_list( + bond_labels=bond_labels, orb_pair=grouped_orb_pairs, root_orb_pair=orb + ) + # for orb in available_orbitals: + cohp_summed_orb = self.chemenv.completecohp.get_summed_cohp_by_label_and_orbital_list( + label_list=mapped_bond_labels, orbital_list=grouped_orb_pairs[orb] * len(bond_labels) + ) + + if type_pop.lower() == "cohp": + ( + antibndg_orb, + per_anti_orb, + bndg_orb, + per_bndg_orb, + ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed_orb, start=self.start) + else: + ( + bndg_orb, + per_bndg_orb, + antibndg_orb, + per_anti_orb, + ) = self._integrate_antbdstates_below_efermi(cohp=cohp_summed_orb, start=self.start) + + # replace nan values with zero (tackle numerical integration issues) + bndg_orb = bndg_orb if not np.isnan(bndg_orb) else 0 + per_bndg_orb = per_bndg_orb if not np.isnan(per_bndg_orb) else 0 + bndg_tot = bndg_tot if not np.isnan(bndg_tot) else 0 + per_bndg_tot = per_bndg_tot if not np.isnan(per_bndg_tot) else 0 + # skip collecting orb contributions if no summed bonding contribution exists + if bndg_tot > 0: + orb_icohps_bndg = [] + for bond_label in bond_labels: + for sub_orb in grouped_orb_pairs[orb]: + orb_icohp_bn = self.chemenv.Icohpcollection.get_icohp_by_label( + label=bond_label, orbitals=sub_orb + ) + orb_icohps_bndg.append(orb_icohp_bn) + bndg_orb_pair_list.append(orb) + bndg_orb_icohp_list.append(orb_icohps_bndg) + bndg_orb_integral_list.append(bndg_orb) + bndg_orb_perc_list.append(per_bndg_orb) + + # replace nan values with zero (tackle numerical integration issues) + antibndg_orb = antibndg_orb if not np.isnan(antibndg_orb) else 0 + per_anti_orb = per_anti_orb if not np.isnan(per_anti_orb) else 0 + antibndg_tot = antibndg_tot if not np.isnan(antibndg_tot) else 0 + per_anti_tot = per_anti_tot if not np.isnan(per_anti_tot) else 0 + # skip collecting orb contributions if no summed antibonding contribution exists + if antibndg_tot > 0: + orb_icohps_anti = [] + for bond_label in bond_labels: + for sub_orb in grouped_orb_pairs[orb]: + orb_icohp_an = self.chemenv.Icohpcollection.get_icohp_by_label( + label=bond_label, orbitals=sub_orb + ) + orb_icohps_anti.append(orb_icohp_an) + + antibndg_orb_pair_list.append(orb) + antibndg_orb_icohp_list.append(orb_icohps_anti) + antibndg_orb_integral_list.append(antibndg_orb) + antibndg_orb_perc_list.append(per_anti_orb) + + # Populate the dictionary with relevant orbitals for bonding interactions + for inx, bndg_orb_pair in enumerate(bndg_orb_pair_list): + bndg_contri_perc = round(bndg_orb_integral_list[inx] / sum(bndg_orb_integral_list), 2) + # filter out very small bonding interactions ( self.orbital_cutoff: + if bndg_orb_pair in orb_bonding_dict_data: + orb_bonding_dict_data[bndg_orb_pair].update( + { + "orb_contribution_perc_bonding": bndg_contri_perc, + "bonding": { + "integral": bndg_orb_integral_list[inx], + "perc": bndg_orb_perc_list[inx], + }, + } + ) + else: + orb_bonding_dict_data[bndg_orb_pair] = { + f"I{type_pop}_mean": round(np.mean(bndg_orb_icohp_list[inx]), 4), + f"I{type_pop}_sum": round(np.sum(bndg_orb_icohp_list[inx]), 4), + "orb_contribution_perc_bonding": round( + bndg_orb_integral_list[inx] / sum(bndg_orb_integral_list), + 2, + ), + "bonding": { + "integral": bndg_orb_integral_list[inx], + "perc": bndg_orb_perc_list[inx], + }, + "relevant_sub_orbitals": grouped_orb_pairs[bndg_orb_pair], + } + + # Populate the dictionary with relevant orbitals for antibonding interactions + for inx, antibndg_orb_pair in enumerate(antibndg_orb_pair_list): + antibndg_contri_perc = round( + antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), + 2, + ) + # filter out very small antibonding interactions ( self.orbital_cutoff: + if antibndg_orb_pair in orb_bonding_dict_data: + orb_bonding_dict_data[antibndg_orb_pair].update( + { + "orb_contribution_perc_antibonding": round( + antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), + 2, + ), + "antibonding": { + "integral": antibndg_orb_integral_list[inx], + "perc": antibndg_orb_perc_list[inx], + }, + } + ) + else: + orb_bonding_dict_data[antibndg_orb_pair] = { + f"I{type_pop}_mean": round(np.mean(antibndg_orb_icohp_list[inx]), 4), + f"I{type_pop}_sum": round(np.sum(antibndg_orb_icohp_list[inx]), 4), + "orb_contribution_perc_antibonding": round( + antibndg_orb_integral_list[inx] / sum(antibndg_orb_integral_list), + 2, + ), + "antibonding": { + "integral": antibndg_orb_integral_list[inx], + "perc": antibndg_orb_perc_list[inx], + }, + "relevant_sub_orbitals": grouped_orb_pairs[antibndg_orb_pair], + } + + orb_bonding_dict_data["relevant_bonds"] = bond_labels # type: ignore + + orb_resolved_bond_info[bond_resolved_label_key] = orb_bonding_dict_data + + return orb_resolved_bond_info + + def _get_bond_resolved_data_stats(self, orb_resolved_bond_data: dict): + """ + Retrieve the maximum bonding and anti-bonding orbital contributions. + + :param orb_resolved_bond_data: A dictionary with orbital names as keys and corresponding bonding data + + Returns: + dict with orbital data stats the site for relevant orbitals, e.g. + {'orbital_summary_stats': {'max_bonding_contribution': {'3p-3p': 0.41}, + 'max_antibonding_contribution': {'3s-3p': 0.39}}} + + """ + # get max orbital bonding and contribution for the site + orb_pairs_bndg = [] + orb_pairs_antibndg = [] + orb_contri_bndg = [] + orb_contri_antibndg = [] + orbital_summary_stats = {"orbital_summary_stats": {}} # type: ignore + if orb_resolved_bond_data: + for orb_pair, data in orb_resolved_bond_data.items(): + if "orb_contribution_perc_bonding" in data: + orb_pairs_bndg.append(orb_pair) + orb_contri_bndg.append(data["orb_contribution_perc_bonding"]) + + if "orb_contribution_perc_antibonding" in data: + orb_pairs_antibndg.append(orb_pair) + orb_contri_antibndg.append(data["orb_contribution_perc_antibonding"]) + + if orb_contri_bndg: + max_orb_contri_bndg = max(orb_contri_bndg) + max_orb_contri_bndg_inxs = [ + inx for inx, orb_contri in enumerate(orb_contri_bndg) if orb_contri == max_orb_contri_bndg + ] + max_orb_contri_bndg_dict = {} + for inx in max_orb_contri_bndg_inxs: + max_orb_contri_bndg_dict[orb_pairs_bndg[inx]] = orb_contri_bndg[inx] + orbital_summary_stats["orbital_summary_stats"]["max_bonding_contribution"] = max_orb_contri_bndg_dict + if orb_contri_antibndg: + max_orb_contri_antibndg = max(orb_contri_antibndg) + max_antibndg_contri_inxs = [ + inx + for inx, orb_anti_per in enumerate(orb_contri_antibndg) + if orb_anti_per == max_orb_contri_antibndg + ] + max_antibndg_contri_dict = {} + for inx in max_antibndg_contri_inxs: + max_antibndg_contri_dict[orb_pairs_antibndg[inx]] = orb_contri_antibndg[inx] + orbital_summary_stats["orbital_summary_stats"]["max_antibonding_contribution"] = ( + max_antibndg_contri_dict + ) + + return orbital_summary_stats + + def get_site_orbital_resolved_labels(self): + """ + Return relevant orbitals and bond labels for each symmetrically independent site. + + Returns: + dict with bond labels for each site for relevant orbitals, e.g. + {'Na1: Na-Cl': {'3p-3s': {'bond_labels': ['21', '23', '24', '27', '28', '30'], + 'relevant_sub_orbitals': ['3py-3s', '3pz-3s', '3px-3s']}} + """ + site_bond_labels = self.get_site_bond_resolved_labels() + orb_plot_data = {atom_pair: {} for atom_pair in site_bond_labels} + if self.orbital_resolved: + for site_index, cba_data in self.condensed_bonding_analysis["sites"].items(): + for atom in cba_data["bonds"]: + for orb_pair in cba_data["bonds"][atom]["orbital_data"]: + if orb_pair not in ("orbital_summary_stats", "relevant_bonds"): + atom_pair = [cba_data["ion"], atom] + atom_pair.sort() + key = ( + self.structure.sites[site_index].species_string + + str(site_index + 1) + + ": " + + "-".join(atom_pair) + ) + label_list = site_bond_labels[key] + relevant_sub_orbitals = cba_data["bonds"][atom]["orbital_data"][orb_pair][ + "relevant_sub_orbitals" + ] + orb_plot_data[key].update( + {orb_pair: {"bond_labels": label_list, "relevant_sub_orbitals": relevant_sub_orbitals}} + ) + else: + print("Please set orbital_resolved to True when instantiating Analysis object, to get this data") + + return orb_plot_data + + @staticmethod + def _get_strenghts_for_each_bond(pairs: list[list[str]], strengths: list[float], nameion: str | None = None): + """ + Return a dictionary of bond strengths. + + :param pairs: list of list including labels for the atoms, e.g., [['O3', 'Cu1'], ['O3', 'Cu1']] + :param strengths: list that gives the icohp strengths as a float, [-1.86287, -1.86288] + :param nameion: string including the name of the cation in the list, e.g Cu1 + + Returns: + dict including inormation on icohps for each bond type, e.g. + {'Yb-Sb': [-1.59769, -2.14723, -1.7925, -1.60773, -1.80149, -2.14335]} + + + """ + dict_strenghts = {} # type: ignore + + for pair, strength in zip(pairs, strengths): + if nameion is not None: + new = [ + LobsterNeighbors._split_string(pair[0])[0], + LobsterNeighbors._split_string(pair[1])[0], + ] + new = Analysis._sort_name(new, nameion) + string_here = new[0] + "-" + new[1] + else: + new = sorted( + [ + LobsterNeighbors._split_string(pair[0])[0], + LobsterNeighbors._split_string(pair[1])[0], + ] + ) + string_here = new[0] + "-" + new[1] + + if string_here not in dict_strenghts: + dict_strenghts[string_here] = [] + dict_strenghts[string_here].append(strength) + return dict_strenghts + + @staticmethod + def _sort_name(pair: list[str], nameion: str | None = None): + """ + Place the cation first in a list of name strings. + + :param pair: ["O","Cu"] + :param nameion: "Cu" + + Returns: + will return list of str, e.g. ["Cu", "O"] + + """ + if nameion is not None: + new = [] + if pair[0] == nameion: + new.append(pair[0]) + new.append(pair[1]) + + elif pair[1] == nameion: + new.append(pair[1]) + new.append(pair[0]) + + return new + + @staticmethod + def _sort_orbital_atom_pair( + atom_pair: list[str], + label: str, + complete_cohp: CompleteCohp, + orb_pair: str, + ): + """ + Place the cation first in a list of name strings and add the associated orbital name alongside the atom name. + + :param atom_pair: list of atom pair with cation first eg., ["Cl","Na"] + :param label: LOBSTER relevant bond label eg ., "3" + :param complete_cohp: pymatgen CompleteCohp object + :param orb_pair: relevant orbital pair eg., "2p-3s" + + Returns: + will return list of str, e.g. ["Na(2p)", "Cl(3s)"] + + """ + orb_atom = {} # type: ignore + orb_pair_list = orb_pair.split("-") + # get orbital associated to the atom and store in a dict + for _inx, (site, site_orb) in enumerate(zip(complete_cohp.bonds[label]["sites"], orb_pair_list)): + if site.species_string in orb_atom: # check necessary for bonds between same atoms + orb_atom[site.species_string].append(site_orb) + else: + orb_atom[site.species_string] = [site_orb] + + orb_atom_list = [] + # add orbital name next to atom_pair + for inx, atom in enumerate(atom_pair): + # check to ensure getting 2nd orbital if bond is between same atomic species + if inx == 1 and len(orb_atom.get(atom)) > 1: # type: ignore + atom_with_orb_name = f"{atom}({orb_atom.get(atom)[1]})" # type: ignore + else: + atom_with_orb_name = f"{atom}({orb_atom.get(atom)[0]})" # type: ignore + orb_atom_list.append(atom_with_orb_name) + + return orb_atom_list + + def _get_antibdg_states(self, cohps, labels: list[str], nameion: str | None = None, limit=0.01): + """ + Return a dictionary containing information on anti-bonding states. + + e.g., similar to: {'Cu-O': True, 'Cu-F': True} + + :param cohps: list of pymatgen.electronic_structure.cohp.Cohp objects + :param labels: ['2 x Cu-O', '4 x Cu-F'] + :param nameion: string of the cation name, e.g. "Cu" + :param limit: limit to detect antibonding states + + Returns: + dict including in formation on whether antibonding interactions exist, + e.g., {'Cu-O': True, 'Cu-F': True} + + + """ + dict_antibd = {} + for label, cohp in zip(labels, cohps): + if label is not None: + if nameion is not None: + new = label.split(" ")[2].split("-") + sorted_new = self._sort_name(new, nameion) + new_label = sorted_new[0] + "-" + sorted_new[1] + else: + new = label.split(" ")[2].split("-") + sorted_new = sorted(new.copy()) + new_label = sorted_new[0] + "-" + sorted_new[1] + + antbd = cohp.has_antibnd_states_below_efermi(limit=limit) + if Spin.down in antbd: + dict_antibd[new_label] = antbd[Spin.up] or antbd[Spin.down] + else: + dict_antibd[new_label] = antbd[Spin.up] + + return dict_antibd + + def _integrate_antbdstates_below_efermi_for_set_cohps(self, labels: list[str], cohps, nameion: str): + """ + Return a dictionary containing information on antibonding states. + + .. warning:: NEEDS MORE TESTS + + It is important to note that only the energy range that has been computed can be considered + (i.e., this might not be all) + + e.g. output: {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995}, + 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}} + + :param cohps: list of pymatgen.electronic_structure.cohp.Cohp objects + :param labels: ['2 x Cu-O', '4 x Cu-F'] + :param nameion: string of the cation name, e.g. "Cu" + + Returns: + dict including in formation on whether antibonding interactions exist, + e.g., {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995}, + 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}}} + """ + dict_bd_antibd = {} + for label, cohp in zip(labels, cohps): + if label is not None: + new = label.split(" ")[2].split("-") + sorted_new = self._sort_name(new, nameion) + new_label = sorted_new[0] + "-" + sorted_new[1] + if not self.are_cobis and not self.are_coops: + ( + integral, + perc, + integral2, + perc2, + ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start) + else: + ( + integral2, + perc2, + integral, + perc, + ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start) + + if integral == 0 and integral2 != 0.0: + dict_bd_antibd[new_label] = { + "bonding": {"integral": integral2, "perc": perc2}, + "antibonding": {"integral": integral, "perc": 0.0}, + } + elif integral2 == 0.0 and integral != 0.0: + dict_bd_antibd[new_label] = { + "bonding": {"integral": integral2, "perc": 0.0}, + "antibonding": {"integral": integral, "perc": perc}, + } + elif integral == 0.0 and integral2 == 0.0: + dict_bd_antibd[new_label] = { + "bonding": {"integral": integral2, "perc": 0.0}, + "antibonding": {"integral": integral, "perc": 0.0}, + } + else: + dict_bd_antibd[new_label] = { + "bonding": {"integral": integral2, "perc": perc2}, + "antibonding": {"integral": integral, "perc": perc}, + } + + return dict_bd_antibd + + def _integrate_antbdstates_below_efermi(self, cohp, start: float | None): + """ + Integrate the cohp data to compute bonding and anti-bonding contribution below efermi. + + .. warning:: NEEDS MORE TESTS + + This integrates the whole COHP curve that has been computed. + The energy range is very important. + At present the energy range considered is dependent on COHPstartEnergy + set during lobster runs. The bonding / antibonding integral values are sensitive to this parameter. + If COHPstartEnergy value does not cover entire range of VASP calculations then + absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values. + + :param cohp: cohp object + :param start: integration start energy in eV , eg start = -15 + + Returns: + absolute value of antibonding, percentage value of antibonding, + absolute value of bonding and percentage value of bonding interactions + """ + warnings.warn( + "The bonding, antibonding integral/percent values are numerical estimate." + " These values are sensitive to COHPstartEnergy parameter." + " If COHPstartEnergy value does not cover entire range of VASP calculations then" + " absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values.", + stacklevel=2, + ) + + from scipy.integrate import trapezoid + + def integrate_positive(y, x): + """ + Integrate only bonding interactions of COHPs. + + :param y: COHP values + :param x: Energy values + + Returns: + integrated value of bonding interactions + """ + y = np.asanyarray(y) + x = np.asanyarray(x) + + bonding = trapezoid(y, x) + + return np.round(bonding, 2) + + def integrate_negative(y, x): + """ + Integrate only anti-bonding interactions of COHPs. + + :param y: COHP values + :param x: Energy values + + Returns: + integrated value of anti-bonding interactions + """ + y = np.asanyarray(y) + x = np.asanyarray(x) + antibonding = trapezoid(y, x) + + return np.round(antibonding, 2) + + # will integrate spin.up and spin.down only below efermi + energies_corrected = cohp.energies - cohp.efermi + summedcohp = cohp.cohp[Spin.up] + cohp.cohp[Spin.down] if Spin.down in cohp.cohp else cohp.cohp[Spin.up] + + cohp_bf = [] + en_bf = [] + + for i, en in enumerate(energies_corrected): + if (start is None) and en <= 0: + en_bf.append(en) + cohp_bf.append(-1 * summedcohp[i]) + if (start is not None) and 0 >= en >= start: + en_bf.append(en) + cohp_bf.append(-1 * summedcohp[i]) + + # Separate the bonding and antibonding COHP values in separate lists + pos = [] + en_pos = [] + neg = [] + en_neg = [] + + for i, scohp in enumerate(cohp_bf): + if scohp >= 0: + pos.append(scohp) + en_pos.append(energies_corrected[i]) + + for i, scohp in enumerate(cohp_bf): + if scohp <= 0: + neg.append(-1 * scohp) + en_neg.append(energies_corrected[i]) + + antibonding = integrate_negative(y=neg, x=en_neg) + + bonding = integrate_positive(y=pos, x=en_pos) + + return ( + antibonding, + np.round(abs(antibonding) / (abs(bonding) + abs(antibonding)), 5), + bonding, + np.round(abs(bonding) / (abs(bonding) + abs(antibonding)), 5), + ) + + def _get_pop_type(self): + """ + Return the type of the input population file. + + Returns: + A String of analysed population can be COOP/COBI/COHP + """ + if self.are_cobis: + type_pop = "COBI" + elif self.are_coops: + type_pop = "COOP" + else: + type_pop = "COHP" + + return type_pop + + @staticmethod + def _get_bond_dict( + bond_strength_dict: dict, + small_antbd_dict: dict, + nameion: str | None = None, + large_antbd_dict: dict | None = None, + type_pop: str | None = None, + ): + """ + Return a bond_dict that contains information for each site. + + :param bond_strength_dict: dict with bond names as key and lists of bond strengths as items + :param small_antbd_dict: dict including if there are antibonding interactions, {'Yb-Sb': False} + :param nameion: name of the cation, e.g. Yb + :param large_antbd_dict: will be implemented later + :param type_pop: population type analyzed. eg. COHP + + Returns: + Eg., if type_pop == 'COHP', will return + dict including information on the anion (as label) and the ICOHPs in the item of the dict + ICOHP_mean refers to the mean ICOHP in eV + ICOHP_sum refers to the sum of the ICOHPs in eV + has_antibdg_states_below_Efermi is True if there are antibonding interactions below Efermi + "number_of_bonds" will count the numbers of bonds to the cation + + Example: + {'Sb': {'ICOHP_mean': '-1.85', 'ICOHP_sum': '-11.09', + 'has_antibdg_states_below_Efermi': False, 'number_of_bonds': 6}} + + """ + bond_dict = {} + + for key, item in bond_strength_dict.items(): + if nameion is not None: + a = key.split("-")[0] + b = key.split("-")[1] + if a == nameion: + key_here = b + elif b == nameion: + key_here = a + + if large_antbd_dict is None: + bond_dict[key_here] = { + f"I{type_pop}_mean": str(round(np.mean(item), 2)), + f"I{type_pop}_sum": str(round(np.sum(item), 2)), + "has_antibdg_states_below_Efermi": small_antbd_dict[key], + "number_of_bonds": len(item), + } + else: + bond_dict[key_here] = { + f"I{type_pop}_mean": str(round(np.mean(item), 2)), + f"I{type_pop}_sum": str(round(np.sum(item), 2)), + "has_antibdg_states_below_Efermi": small_antbd_dict[key], + "number_of_bonds": len(item), + "perc_antibdg_states_below_Efermi": large_antbd_dict[key], + } + + return bond_dict + + @staticmethod + def _get_orb_combinations(): + """ + Get a list of unique 2-length permutations and combinations. + + Generates 2-length permutations and combinations with replacement from the set ['s', 'p', 'd', 'f'] + + Returns: + list[tuple[str, str]]: A list of tuples, each containing two elements. + """ + list_orbs_comb: list[tuple[str, str]] = [] # type: ignore + + # Add all 2-length permutations + for perm in permutations(["s", "p", "d", "f"], 2): + list_orbs_comb.append(perm) # noqa : PERF402 + + # Add 2-length combinations with replacement, ensuring no duplicates + for comb in combinations_with_replacement(["s", "p", "d", "f"], 2): + if comb not in list_orbs_comb: + list_orbs_comb.append(comb) + + return list_orbs_comb + + def _group_orb_pairs(self, bond_label: str, orb_combinations: list[tuple[str, str]]) -> dict[str, list[str]]: + """ + Group orbital pairs based on the provided bond label. + + :param bond_label: The bond label to filter the orbitals. + :param orb_combinations: A list of tuples containing orbital combinations. + + Returns: + dict[str, List[str]]: A dictionary where the keys are top level orbital pairs and the values are lists of + sub orbitals associated to the bond label. + """ + orb_pair = {} # type: ignore + bond_label_int = int(bond_label) - 1 # convert label to int to access orbital data from icohpcollection + for sub_orb, data in self.chemenv.Icohpcollection._list_orb_icohp[bond_label_int].items(): + for orb1, orb2 in orb_combinations: + if orb1 in data["orbitals"][0][1].name and orb2 in data["orbitals"][1][1].name: + root_orb_pair = f"{data['orbitals'][0][0]}{orb1}-{data['orbitals'][1][0]}{orb2}" + if root_orb_pair not in orb_pair: + orb_pair[root_orb_pair] = [] + orb_pair[root_orb_pair].append(sub_orb) + return orb_pair + + @staticmethod + def _get_bond_orb_label_mapped_list(orb_pair: dict[str, list[str]], bond_labels: list[str], root_orb_pair: str): + """ + Get a lists of bond labels mapped to the corresponding orbital pair. + + :param orb_pair: A dictionary containing orbital pairs as keys and lists of + sub orbitals as values. + :param bond_labels: A list of bond labels. + :param root_orb_pair: The root key in orb_pair use to map bond labels list. + + Returns: + list: A list where the items of bond_labels are repeated based on the + length of orb_pair[root_orb_pair]. + """ + return [item for item in bond_labels for _ in range(len(orb_pair[root_orb_pair]))] + + def set_condensed_bonding_analysis(self): + """ + Condense the bonding analysis into a summary dictionary. + + Returns: + None + + """ + self.condensed_bonding_analysis = {} + # which icohps are considered + if self.which_bonds == "cation-anion": + limit_icohps = self.chemenv._get_limit_from_extremum( + self.chemenv.Icohpcollection, + self.cutoff_icohp, + adapt_extremum_to_add_cond=True, + additional_condition=1, + ) + elif self.which_bonds == "all": + limit_icohps = self.chemenv._get_limit_from_extremum( + self.chemenv.Icohpcollection, + self.cutoff_icohp, + adapt_extremum_to_add_cond=True, + additional_condition=0, + ) + # formula of the compound + formula = str(self.structure.composition.reduced_formula) + # set population type + type_pop = self._get_pop_type() + # how many inequivalent cations are in the structure + if self.which_bonds == "cation-anion": + number_considered_ions = len(self.seq_ineq_ions) + elif self.which_bonds == "all": + number_considered_ions = len(self.seq_ineq_ions) + + # what was the maximum bond lengths that was considered + max_bond_lengths = max(self.chemenv.Icohpcollection._list_length) + + # what are the charges for the cations in the structure + charge_list = self.chemenv.valences + + # dictionary including bonding information for each site + site_dict = {} + if self.which_bonds == "cation-anion": + for ication, ce, cation_anion_infos, labels, cohps in zip( + self.seq_ineq_ions, + self.seq_coord_ions, + self.seq_infos_bonds, + self.seq_labels_cohps, + self.seq_cohps, + ): + namecation = str(self.structure[ication].specie) + + # This will compute the mean strengths of ICOHPs + mean_icohps = self._get_strenghts_for_each_bond( + pairs=cation_anion_infos[4], + strengths=cation_anion_infos[1], + nameion=namecation, + ) + # pairs, strengths, nameion + # will collect if there are antibonding states present + antbdg = self._get_antibdg_states(cohps, labels, namecation) + dict_antibonding = self._integrate_antbdstates_below_efermi_for_set_cohps( + labels, cohps, nameion=namecation + ) + bond_dict = self._get_bond_dict(mean_icohps, antbdg, namecation, type_pop=type_pop) + bond_resolved_labels = self.get_site_bond_resolved_labels() + + for cation_name, icohp_data in bond_dict.items(): + for atom_pair, bonding_data in dict_antibonding.items(): + if namecation == atom_pair.split("-")[0] and cation_name == atom_pair.split("-")[1]: + icohp_data["bonding"] = bonding_data["bonding"] + icohp_data["antibonding"] = bonding_data["antibonding"] + if self.orbital_resolved: + # get orb resolved data to be added + orb_resolved_bond_info = self._get_orbital_resolved_data( + nameion=namecation, + iion=ication, + labels=labels, + bond_resolved_labels=bond_resolved_labels, + type_pop=type_pop, + ) + # match the dict key in bond_dict and get corresponding orbital data + for ion_atom_pair_orb in orb_resolved_bond_info: + orb_data_atom_pair = ion_atom_pair_orb.split(": ")[-1] + atom_pair_here = atom_pair.split("-") + atom_pair_here.sort() + if ( + orb_data_atom_pair == "-".join(atom_pair_here) + and (namecation + str(ication + 1) + ":") in ion_atom_pair_orb + ): + icohp_data["orbital_data"] = orb_resolved_bond_info[ion_atom_pair_orb] + + orb_data_stats = self._get_bond_resolved_data_stats( + orb_resolved_bond_data=orb_resolved_bond_info[ion_atom_pair_orb], + ) + + icohp_data["orbital_data"].update(orb_data_stats) + + site_dict[ication] = { + "env": ce, + "bonds": bond_dict, + "ion": namecation, + "charge": charge_list[ication], + "relevant_bonds": cation_anion_infos[3], + } + elif self.which_bonds == "all": + for iion, ce, bond_infos, labels, cohps in zip( + self.seq_ineq_ions, + self.seq_coord_ions, + self.seq_infos_bonds, + self.seq_labels_cohps, + self.seq_cohps, + ): + nameion = str(self.structure[iion].specie) + + # This will compute the mean strengths of ICOHPs + mean_icohps = self._get_strenghts_for_each_bond( + pairs=bond_infos[4], strengths=bond_infos[1], nameion=None + ) + # pairs, strengths, nameion + # will collect if there are antibonding states present + antbdg = self._get_antibdg_states(cohps, labels, nameion=None) + + dict_antibonding = self._integrate_antbdstates_below_efermi_for_set_cohps(labels, cohps, nameion) + + bond_dict = self._get_bond_dict(mean_icohps, antbdg, nameion=nameion, type_pop=type_pop) + bond_resolved_labels = self.get_site_bond_resolved_labels() + + for cation_name, icohp_data in bond_dict.items(): + for atom_pair, bonding_data in dict_antibonding.items(): + if nameion == atom_pair.split("-")[0] and cation_name == atom_pair.split("-")[1]: + icohp_data["bonding"] = bonding_data["bonding"] + icohp_data["antibonding"] = bonding_data["antibonding"] + if self.orbital_resolved: + # get orb resolved data to be added + orb_resolved_bond_info = self._get_orbital_resolved_data( + nameion=nameion, + iion=iion, + labels=labels, + bond_resolved_labels=bond_resolved_labels, + type_pop=type_pop, + ) + # match the dict key in bond_dict and get corresponding orbital data + for ion_atom_pair_orb in orb_resolved_bond_info: + orb_data_atom_pair = ion_atom_pair_orb.split(": ")[-1] + atom_pair_here = atom_pair.split("-") + atom_pair_here.sort() + if ( + orb_data_atom_pair == "-".join(atom_pair_here) + and (nameion + str(iion + 1) + ":") in ion_atom_pair_orb + ): + icohp_data["orbital_data"] = orb_resolved_bond_info[ion_atom_pair_orb] + + orb_data_stats = self._get_bond_resolved_data_stats( + orb_resolved_bond_data=orb_resolved_bond_info[ion_atom_pair_orb], + ) + + icohp_data["orbital_data"].update(orb_data_stats) + + site_dict[iion] = { + "env": ce, + "bonds": bond_dict, + "ion": nameion, + "charge": charge_list[iion], + "relevant_bonds": bond_infos[3], + } + + if self.path_to_madelung is None and self.madelung_obj is None: + if self.which_bonds == "cation-anion": + # This sets the dictionary including the most important information on the compound + self.condensed_bonding_analysis = { + "formula": formula, + "max_considered_bond_length": max_bond_lengths, + f"limit_i{type_pop.lower()}": limit_icohps, + "number_of_considered_ions": number_considered_ions, + "sites": site_dict, + "type_charges": self.type_charge, + } + elif self.which_bonds == "all": + self.condensed_bonding_analysis = { + "formula": formula, + "max_considered_bond_length": max_bond_lengths, + f"limit_i{type_pop.lower()}": limit_icohps, + "number_of_considered_ions": number_considered_ions, + "sites": site_dict, + "type_charges": self.type_charge, + } + else: + madelung = MadelungEnergies(self.path_to_madelung) if self.path_to_madelung else self.madelung_obj + if self.type_charge == "Mulliken": + madelung_energy = madelung.madelungenergies_mulliken + elif self.type_charge == "Löwdin": + madelung_energy = madelung.madelungenergies_loewdin + else: + madelung_energy = None + # This sets the dictionary including the most important information on the compound + if self.which_bonds == "cation-anion": + self.condensed_bonding_analysis = { + "formula": formula, + "max_considered_bond_length": max_bond_lengths, + f"limit_i{type_pop.lower()}": limit_icohps, + "number_of_considered_ions": number_considered_ions, + "sites": site_dict, + "type_charges": self.type_charge, + "madelung_energy": madelung_energy, + } + elif self.which_bonds == "all": + self.condensed_bonding_analysis = { + "formula": formula, + "max_considered_bond_length": max_bond_lengths, + f"limit_i{type_pop.lower()}": limit_icohps, + "number_of_considered_ions": number_considered_ions, + "sites": site_dict, + "type_charges": self.type_charge, + "madelung_energy": madelung_energy, + } + + def set_summary_dicts(self): + """ + Set summary dict that can be used for correlations. + + bond_dict that includes information on each bond + + "has_antbd" tells if there are antbonding states + "ICOHP_mean" shows the mean of all ICOHPs in EV + + {'Yb-Sb': { 'has_antbdg': False, 'ICOHP_mean': -1.7448}, + 'Mn-Sb': { 'has_antbdg': True, 'ICOHP_mean': -1.525}} + + a cation dict that includes all different coordination environments and counts for them + {'Na': {'T:4': 4, 'A:2': 4}, 'Si': {'T:6': 4, 'PP:6': 4}} + + Returns: + None + + """ + relevant_ion_ids = [isite for isite in self.list_equivalent_sites if isite in self.seq_ineq_ions] + # set population type + type_pop = self._get_pop_type() + + final_dict_bonds = {} + for key in relevant_ion_ids: + item = self.condensed_bonding_analysis["sites"][key] + for type, properties in item["bonds"].items(): + label_list = [item["ion"], str(type)] + new_label = sorted(label_list.copy()) + label = str(new_label[0]) + "-" + str(new_label[1]) + + if label not in final_dict_bonds: + final_dict_bonds[label] = { + "number_of_bonds": int(properties["number_of_bonds"]), + f"I{type_pop}_sum": float(properties[f"I{type_pop}_sum"]), + "has_antbdg": properties["has_antibdg_states_below_Efermi"], + } + else: + final_dict_bonds[label]["number_of_bonds"] += int(properties["number_of_bonds"]) + final_dict_bonds[label][f"I{type_pop}_sum"] += float(properties[f"I{type_pop}_sum"]) + final_dict_bonds[label]["has_antbdg"] = ( + final_dict_bonds[label]["has_antbdg"] or properties["has_antibdg_states_below_Efermi"] + ) + self.final_dict_bonds = {} + for key, item in final_dict_bonds.items(): + self.final_dict_bonds[key] = {} + self.final_dict_bonds[key][f"I{type_pop}_mean"] = item[f"I{type_pop}_sum"] / (item["number_of_bonds"]) + self.final_dict_bonds[key]["has_antbdg"] = item["has_antbdg"] + + # rework, add all environments! + final_dict_ions = {} + for key in relevant_ion_ids: + if self.condensed_bonding_analysis["sites"][key]["ion"] not in final_dict_ions: + final_dict_ions[self.condensed_bonding_analysis["sites"][key]["ion"]] = [ + self.condensed_bonding_analysis["sites"][key]["env"] + ] + else: + final_dict_ions[self.condensed_bonding_analysis["sites"][key]["ion"]].append( + self.condensed_bonding_analysis["sites"][key]["env"] + ) + + self.final_dict_ions = {} + for key, item in final_dict_ions.items(): + self.final_dict_ions[key] = dict(Counter(item)) + + @staticmethod + def get_lobster_calc_quality_summary( + path_to_poscar: str | None = None, + path_to_lobsterout: str | None = None, + path_to_lobsterin: str | None = None, + path_to_potcar: str | None = None, + potcar_symbols: list | None = None, + path_to_charge: str | None = None, + path_to_bandoverlaps: str | None = None, + path_to_doscar: str | None = None, + path_to_vasprun: str | None = None, + structure_obj: Structure | None = None, + lobsterin_obj: Lobsterin | None = None, + lobsterout_obj: Lobsterout | None = None, + charge_obj: Charge | None = None, + bandoverlaps_obj: Bandoverlaps | None = None, + lobster_completedos_obj: LobsterCompleteDos | None = None, + vasprun_obj: Vasprun | None = None, + dos_comparison: bool = False, + e_range: list = [-5, 0], + n_bins: int | None = None, + bva_comp: bool = False, + ) -> dict: + """ + Analyze LOBSTER calculation quality. + + :param path_to_poscar: path to structure file + :param path_to_lobsterout: path to lobsterout file + :param path_to_lobsterin: path to lobsterin file + :param path_to_potcar: path to VASP potcar file + :param potcar_symbols: list of potcar symbols from potcar file (can be used if no potcar available) + :param path_to_charge: path to CHARGE.lobster file + :param path_to_bandoverlaps: path to bandOverlaps.lobster file + :param path_to_doscar: path to DOSCAR.lobster or DOSCAR.LSO.lobster file + :param path_to_vasprun: path to vasprun.xml file + :param structure_obj: pymatgen pymatgen.core.structure.Structure object + :param lobsterin_obj: pymatgen.lobster.io.Lobsterin object + :param lobsterout_obj: pymatgen lobster.io.Lobsterout object + :param charge_obj: pymatgen lobster.io.Charge object + :param bandoverlaps_obj: pymatgen lobster.io.BandOverlaps object + :param lobster_completedos_obj: pymatgen.electronic_structure.dos.LobsterCompleteDos object + :param vasprun_obj: pymatgen vasp.io.Vasprun object + :param dos_comparison: will compare DOS from VASP and LOBSTER and return tanimoto index + :param e_range: energy range for DOS comparisons + :param n_bins: number of bins to discretize DOS for comparisons + :param bva_comp: Compares LOBSTER charge signs with Bond valence charge signs + + Returns: + A dict of summary of LOBSTER calculation quality by analyzing basis set used, + charge spilling from lobsterout/ PDOS comparisons of VASP and LOBSTER / + BVA charge comparisons + + """ + quality_dict = {} + + if path_to_potcar and not potcar_symbols and not path_to_vasprun and not vasprun_obj: + potcar_names = Lobsterin._get_potcar_symbols(POTCAR_input=path_to_potcar) + elif not path_to_potcar and not path_to_vasprun and not vasprun_obj and potcar_symbols: + potcar_names = potcar_symbols + elif path_to_vasprun and not vasprun_obj: + vasprun = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False) + potcar_names = [potcar.split(" ")[1] for potcar in vasprun.potcar_symbols] + elif vasprun_obj and not path_to_vasprun: + potcar_names = [potcar.split(" ")[1] for potcar in vasprun_obj.potcar_symbols] + else: + raise ValueError( + "Please provide either path_to_potcar or list of " + "potcar_symbols or path to vasprun.xml or vasprun object. " + "Crucial to identify basis used for projections" + ) + + if path_to_poscar: + struct = Structure.from_file(path_to_poscar) + elif structure_obj: + struct = structure_obj + else: + raise ValueError("Please provide path_to_poscar or structure_obj") + + ref_bases = Lobsterin.get_all_possible_basis_functions(structure=struct, potcar_symbols=potcar_names) + + if path_to_lobsterin: + lobs_in = Lobsterin.from_file(path_to_lobsterin) + elif lobsterin_obj: + lobs_in = lobsterin_obj + else: + raise ValueError("Please provide path_to_lobsterin or lobsterin_obj") + + calc_basis = [] + for basis in lobs_in["basisfunctions"]: + basis_sep = basis.split()[1:] + basis_comb = " ".join(basis_sep) + calc_basis.append(basis_comb) + + if calc_basis == list(ref_bases[0].values()): + quality_dict["minimal_basis"] = True # type: ignore + else: + quality_dict["minimal_basis"] = False # type: ignore + warnings.warn( + "Consider rerunning the calc with the minimum basis as well. Choosing is " + "larger basis set is recommended if you see a significant improvement of " + "the charge spilling and material has non-zero band gap.", + stacklevel=2, + ) + + if path_to_lobsterout: + lob_out = Lobsterout(path_to_lobsterout) + elif lobsterout_obj: + lob_out = lobsterout_obj + else: + raise ValueError("Please provide path_to_lobsterout or lobsterout_obj") + + quality_dict["charge_spilling"] = { + "abs_charge_spilling": round((sum(lob_out.charge_spilling) / 2) * 100, 4), + "abs_total_spilling": round((sum(lob_out.total_spilling) / 2) * 100, 4), + } # type: ignore + + if path_to_bandoverlaps is not None and not bandoverlaps_obj: + band_overlaps = Bandoverlaps(filename=path_to_bandoverlaps) if Path(path_to_bandoverlaps).exists() else None + elif path_to_bandoverlaps is None and bandoverlaps_obj: + band_overlaps = bandoverlaps_obj + else: + band_overlaps = None + + if band_overlaps is not None: + for line in lob_out.warning_lines: + if "k-points could not be orthonormalized" in line: + total_kpoints = int(line.split(" ")[2]) + + # store actual number of devations above pymatgen default limit of 0.1 + dev_val = [] + for dev in band_overlaps.max_deviation: + if dev > 0.1: + dev_val.append(dev) + + quality_dict["band_overlaps_analysis"] = { # type: ignore + "file_exists": True, + "limit_maxDeviation": 0.1, + "has_good_quality_maxDeviation": band_overlaps.has_good_quality_maxDeviation(limit_maxDeviation=0.1), + "max_deviation": round(max(band_overlaps.max_deviation), 4), + "percent_kpoints_abv_limit": round((len(dev_val) / total_kpoints) * 100, 4), + } + + else: + quality_dict["band_overlaps_analysis"] = { # type: ignore + "file_exists": False, + "limit_maxDeviation": None, + "has_good_quality_maxDeviation": True, + "max_deviation": None, + "percent_kpoints_abv_limit": None, + } + + if bva_comp: + try: + bond_valence = BVAnalyzer() + + bva_oxi = [] + if path_to_charge and not charge_obj: + lobs_charge = Charge(filename=path_to_charge) + elif not path_to_charge and charge_obj: + lobs_charge = charge_obj + else: + raise Exception("BVA comparison is requested, thus please provide path_to_charge or charge_obj") + for i in bond_valence.get_valences(structure=struct): + if i >= 0: + bva_oxi.append("POS") + else: + bva_oxi.append("NEG") + + mull_oxi = [] + for i in lobs_charge.Mulliken: + if i >= 0: + mull_oxi.append("POS") + else: + mull_oxi.append("NEG") + + loew_oxi = [] + for i in lobs_charge.Loewdin: + if i >= 0: + loew_oxi.append("POS") + else: + loew_oxi.append("NEG") + + quality_dict["charge_comparisons"] = {} # type: ignore + if mull_oxi == bva_oxi: + quality_dict["charge_comparisons"]["bva_mulliken_agree"] = True # type: ignore + else: + quality_dict["charge_comparisons"]["bva_mulliken_agree"] = False # type: ignore + + if mull_oxi == bva_oxi: + quality_dict["charge_comparisons"]["bva_loewdin_agree"] = True # type: ignore + else: + quality_dict["charge_comparisons"]["bva_loewdin_agree"] = False # type: ignore + + except ValueError: + quality_dict["charge_comparisons"] = {} # type: ignore + warnings.warn( + "Oxidation states from BVA analyzer cannot be determined. " + "Thus BVA charge comparison will be skipped", + stacklevel=2, + ) + if dos_comparison: + if "LSO" not in str(path_to_doscar).split("."): + warnings.warn( + "Consider using DOSCAR.LSO.lobster, as non LSO DOS from LOBSTER can have negative DOS values", + stacklevel=2, + ) + if path_to_doscar: + doscar_lobster = Doscar( + doscar=path_to_doscar, + structure_file=path_to_poscar, + structure=structure_obj, + ) + + dos_lobster = doscar_lobster.completedos + elif lobster_completedos_obj: + dos_lobster = lobster_completedos_obj + else: + raise ValueError( + "Dos comparison is requested, so please provide either path_to_doscar or lobster_completedos_obj" + ) + + if path_to_vasprun: + vasprun = Vasprun(path_to_vasprun, parse_potcar_file=False, parse_eigen=False) + elif vasprun_obj: + vasprun = vasprun_obj + else: + raise ValueError( + "Dos comparison is requested, so please provide either path to vasprun.xml or vasprun_obj" + ) + dos_vasp = vasprun.complete_dos + + quality_dict["dos_comparisons"] = {} # type: ignore + + for orb in dos_lobster.get_spd_dos(): + if e_range[0] >= min(dos_vasp.energies) and e_range[0] >= min(dos_lobster.energies): + min_e = e_range[0] + else: + warnings.warn( + "Minimum energy range requested for DOS comparisons is not available " + "in VASP or LOBSTER calculation. Thus, setting min_e to -5 eV", + stacklevel=2, + ) + min_e = -5 + + if e_range[-1] <= max(dos_vasp.energies) and e_range[-1] <= max(dos_lobster.energies): + max_e = e_range[-1] + else: + warnings.warn( + "Maximum energy range requested for DOS comparisons is not available " + "in VASP or LOBSTER calculation. Thus, setting max_e to 0 eV", + stacklevel=2, + ) + max_e = 0 + + if np.diff(dos_vasp.energies)[0] >= 0.1 or np.diff(dos_lobster.energies)[0] >= 0.1: + warnings.warn( + "Input DOS files have very few points in the energy interval and thus " + "comparisons will not be reliable. Please rerun the calculations with " + "higher number of DOS points. Set NEDOS and COHPSteps tags to >= 2000 in VASP and LOBSTER " + "calculations, respectively.", + stacklevel=2, + ) + + if not n_bins: + n_bins = 56 + + fp_lobster_orb = dos_lobster.get_dos_fp( + min_e=min_e, + max_e=max_e, + n_bins=n_bins, + normalize=True, + type=orb.name, + ) + fp_vasp_orb = dos_vasp.get_dos_fp( + min_e=min_e, + max_e=max_e, + n_bins=n_bins, + normalize=True, + type=orb.name, + ) + + tani_orb = round( + dos_vasp.get_dos_fp_similarity(fp_lobster_orb, fp_vasp_orb, tanimoto=True), + 4, + ) + quality_dict["dos_comparisons"][f"tanimoto_orb_{orb.name}"] = tani_orb # type: ignore + + fp_lobster = dos_lobster.get_dos_fp( + min_e=min_e, + max_e=max_e, + n_bins=n_bins, + normalize=True, + type="summed_pdos", + ) + fp_vasp = dos_vasp.get_dos_fp( + min_e=min_e, + max_e=max_e, + n_bins=n_bins, + normalize=True, + type="summed_pdos", + ) + + tanimoto_summed = round(dos_vasp.get_dos_fp_similarity(fp_lobster, fp_vasp, tanimoto=True), 4) + quality_dict["dos_comparisons"]["tanimoto_summed"] = tanimoto_summed # type: ignore + quality_dict["dos_comparisons"]["e_range"] = [min_e, max_e] # type: ignore + quality_dict["dos_comparisons"]["n_bins"] = n_bins # type: ignore + + return quality_dict diff --git a/lobsterpy/cohp/describe.py b/src/lobsterpy/cohp/describe.py similarity index 97% rename from lobsterpy/cohp/describe.py rename to src/lobsterpy/cohp/describe.py index 07fbca10..41f22a2f 100644 --- a/lobsterpy/cohp/describe.py +++ b/src/lobsterpy/cohp/describe.py @@ -1,857 +1,857 @@ -# Copyright (c) lobsterpy development team -# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License - -"""This module defines classes to describe the COHPs automatically.""" - -from __future__ import annotations - -import warnings -from pathlib import Path - -from lobsterpy.plotting import InteractiveCohpPlotter, PlainCohpPlotter - - -class Description: - """ - Base class that will write generate a text description for all relevant bonds. - - It analyses all relevant coordination environments in the system based on electronic structure theory. - - """ - - def __init__(self, analysis_object): - """ - Generate a text description for all relevant bonds. - - :param analysis_object: Analysis object from lobsterpy.analysis - """ - self.analysis_object = analysis_object - self.set_description() - - def set_description(self): - """ - Set the descriptions of the structures using the cation names, starting with numbers at 1. - - Uses the cation names from the lobster files. - - Returns: - None - - """ - self.condensed_bonding_analysis = self.analysis_object.condensed_bonding_analysis - # set type of population analyzed - type_pop = self.analysis_object._get_pop_type() - # set units for populations - units = " eV" if type_pop == "COHP" else "" - if self.analysis_object.which_bonds == "cation-anion": - relevant_cations = ", ".join( - [ - str(site.specie) + str(isite + 1) - for isite, site in enumerate(self.analysis_object.structure) - if isite in self.analysis_object.seq_ineq_ions - ] - ) - self.text = [] - self.text.append( - "The compound " - + str(self.condensed_bonding_analysis["formula"]) - + " has " - + str(self.condensed_bonding_analysis["number_of_considered_ions"]) - + " symmetry-independent cation(s) with relevant cation-anion interactions: " - + str(relevant_cations) - + "." - ) - - for key, item in self.condensed_bonding_analysis["sites"].items(): - # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi), - bond_info = [] - orb_info = [] - for type, properties in item["bonds"].items(): - if not properties["has_antibdg_states_below_Efermi"]: - bond_info.append( - str(properties["number_of_bonds"]) - + " " - + item["ion"] - + "-" - + str(type) - + f" (mean I{type_pop}: " - "" - + properties[f"I{type_pop}_mean"] - + f"{units}, 0.0 percent antibonding interaction below EFermi)" - ) - if self.analysis_object.orbital_resolved: - text_orbital = self._generate_orbital_resolved_analysis_text( - orbital_resolved_data=properties, - type_pop=type_pop, - atom_name=str(type), - ion=item["ion"], - ) - orb_info.extend(text_orbital) - else: - bond_info.append( - str(properties["number_of_bonds"]) - + " " - + item["ion"] - + "-" - + str(type) - + f" (mean I{type_pop}: " - "" - + properties[f"I{type_pop}_mean"] - + f"{units}, " - + str(round(properties["antibonding"]["perc"] * 100, 3)) - + " percent antibonding interaction below EFermi)" - ) - if self.analysis_object.orbital_resolved: - text_orbital = self._generate_orbital_resolved_analysis_text( - orbital_resolved_data=properties, - type_pop=type_pop, - atom_name=str(type), - ion=item["ion"], - ) - orb_info.extend(text_orbital) - - bonds = ",".join(bond_info[0:-1]) + ", and " + bond_info[-1] if len(bond_info) > 1 else bond_info[0] - - if len(orb_info) > 1: - orb_bonds = "".join(orb_info).replace(".In", ". In") - else: - orb_bonds = orb_info[0] if orb_info else "" - if item["env"] == "O:6": - self.text.append( - str(item["ion"]) - + str(key + 1) - + " has an " - + str(self._coordination_environment_to_text(item["env"])) - + " coordination environment. It has " - + str(bonds) - + " bonds." - ) - if orb_bonds: - self.text.append(orb_bonds) - else: - self.text.append( - str(item["ion"]) - + str(key + 1) - + " has a " - + str(self._coordination_environment_to_text(item["env"])) - + " coordination environment. It has " - + str(bonds) - + " bonds." - ) - if orb_bonds: - self.text.append(orb_bonds) - - elif self.analysis_object.which_bonds == "all": - relevant_ions = ", ".join( - [ - str(site.specie) + str(isite + 1) - for isite, site in enumerate(self.analysis_object.structure) - if isite in self.analysis_object.seq_ineq_ions - ] - ) - self.text = [] - self.text.append( - "The compound " - + str(self.condensed_bonding_analysis["formula"]) - + " has " - + str(self.condensed_bonding_analysis["number_of_considered_ions"]) - + " symmetry-independent atoms(s) with relevant bonds: " - + str(relevant_ions) - + "." - ) - - for key, item in self.condensed_bonding_analysis["sites"].items(): - # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi), - bond_info = [] - orb_info = [] - for type, properties in item["bonds"].items(): - if not properties["has_antibdg_states_below_Efermi"]: - bond_info.append( - str(properties["number_of_bonds"]) - + " " - + item["ion"] - + "-" - + str(type) - + f" (mean I{type_pop}: " - "" - + properties[f"I{type_pop}_mean"] - + f"{units}, 0.0 percent antibonding interaction below EFermi)" - ) - if self.analysis_object.orbital_resolved: - text_orbital = self._generate_orbital_resolved_analysis_text( - orbital_resolved_data=properties, - type_pop=type_pop, - atom_name=str(type), - ion=item["ion"], - ) - orb_info.extend(text_orbital) - else: - bond_info.append( - str(properties["number_of_bonds"]) - + " " - + item["ion"] - + "-" - + str(type) - + f" (mean I{type_pop}: " - "" - + properties[f"I{type_pop}_mean"] - + f"{units}, " - + str(round(properties["antibonding"]["perc"] * 100, 3)) - + " percent antibonding interaction below EFermi)" - ) - - if self.analysis_object.orbital_resolved: - text_orbital = self._generate_orbital_resolved_analysis_text( - orbital_resolved_data=properties, - type_pop=type_pop, - atom_name=str(type), - ion=item["ion"], - ) - orb_info.extend(text_orbital) - - bonds = ",".join(bond_info[0:-1]) + ", and " + bond_info[-1] if len(bond_info) > 1 else bond_info[0] - - if len(orb_info) > 1: - orb_bonds = "".join(orb_info).replace(".In", ". In") - else: - orb_bonds = orb_info[0] if orb_info else "" - if item["env"] == "O:6": - self.text.append( - str(item["ion"]) - + str(key + 1) - + " has an " - + str(self._coordination_environment_to_text(item["env"])) - + " coordination environment. It has " - + str(bonds) - + " bonds." - ) - if orb_bonds: - self.text.append(orb_bonds) - else: - self.text.append( - str(item["ion"]) - + str(key + 1) - + " has a " - + str(self._coordination_environment_to_text(item["env"])) - + " coordination environment. It has " - + str(bonds) - + " bonds." - ) - if orb_bonds: - self.text.append(orb_bonds) - - if ( - "madelung_energy" in self.analysis_object.condensed_bonding_analysis - and self.analysis_object.condensed_bonding_analysis["madelung_energy"] is not None - ): - self.text.append( - "The Madelung energy of this crystal structure per unit cell is: " - + str(self.analysis_object.condensed_bonding_analysis["madelung_energy"]) - + " eV." - ) - - def _generate_orbital_resolved_analysis_text( - self, - orbital_resolved_data: dict, - ion: str, - atom_name: str, - type_pop: str, - ): - """ - Generate text from orbital-resolved analysis data of the most relevant COHP, COOP, or COBI. - - :param orbital_resolved_data: dict of orbital data from condensed bonding analysis object - :param ion: name of ion at the site - :param atom_name: name of atomic speice to which ion is bonded - :param type_pop: population type analysed could be "COHP" or "COOP" or "COBI" - - Returns: - A python list with text describing the orbital which contributes - the most to the bonding and antibonding in the bond at site - """ - orb_info = [] - if orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: - orb_names = [] - orb_contri = [] - # get atom-pair list with ion placed first - atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion) - if "max_bonding_contribution" in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: - for orb, data in orbital_resolved_data["orbital_data"]["orbital_summary_stats"][ - "max_bonding_contribution" - ].items(): - atom_pair_with_orb_name = self.analysis_object._sort_orbital_atom_pair( - atom_pair=atom_pair, - complete_cohp=self.analysis_object.chemenv.completecohp, - label=orbital_resolved_data["orbital_data"]["relevant_bonds"][0], - orb_pair=orb, - ) - orb_names.append("-".join(atom_pair_with_orb_name)) - orb_contri.append( - str( - round( - data * 100, - 3, - ) - ) - ) - orb_names_anti = [] - orb_antibonding = [] - if "max_antibonding_contribution" in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: - for orb, data in orbital_resolved_data["orbital_data"]["orbital_summary_stats"][ - "max_antibonding_contribution" - ].items(): - atom_pair_with_orb_name = self.analysis_object._sort_orbital_atom_pair( - atom_pair=atom_pair, - complete_cohp=self.analysis_object.chemenv.completecohp, - label=orbital_resolved_data["orbital_data"]["relevant_bonds"][0], - orb_pair=orb, - ) - orb_names_anti.append("-".join(atom_pair_with_orb_name)) - orb_antibonding.append( - str( - round( - data * 100, - 3, - ) - ) - ) - if len(orb_contri) > 1: - orb_name_contri = "" - for inx, name in enumerate(orb_names): - if len(orb_contri) == 2 and inx + 1 != len(orb_contri): - orb_name_contri += f"{name} " - elif 2 < len(orb_contri) != inx + 1: - orb_name_contri += f"{name}, " - else: - orb_name_contri += f"and {name}" - - orb_name_contri += " orbitals, contributing " - for inx, contribution in enumerate(orb_contri): - if len(orb_contri) == 2 and inx + 1 != len(orb_contri): - orb_name_contri += f"{contribution} " - elif 2 < len(orb_contri) != inx + 1: - orb_name_contri += f"{contribution}, " - else: - orb_name_contri += f"and {contribution} percent, respectively" - num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) - bonds = "bonds" if num_bonds > 1 else "bond" - orb_info.append( - f"In the {num_bonds} " - + "-".join(atom_pair) - + f" {bonds}, relative to the summed I{type_pop}s, " - + "the maximum bonding contribution is from " - + orb_name_contri - ) - elif not orb_contri: - num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) - bonds = "bonds" if num_bonds > 1 else "bond" - orb_info.append( - f"In the {num_bonds} " - + "-".join(atom_pair) - + f" {bonds}, relative to the summed I{type_pop}s, " - + f"no orbital has a bonding contribution greater than " - f"{self.analysis_object.orbital_cutoff*100} percent" - ) - else: - num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) - bonds = "bonds" if num_bonds > 1 else "bond" - orb_info.append( - f"In the {num_bonds} " - + "-".join(atom_pair) - + f" {bonds}, relative to the summed I{type_pop}s, " - + "the maximum bonding contribution is from the " - + f"{orb_names[0]}" - + f" orbital, contributing {orb_contri[0]} percent" - ) - - if len(orb_antibonding) > 1: - orb_anti = "" - for inx, name in enumerate(orb_names_anti): - if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti): - orb_anti += f"{name} " - elif 2 < len(orb_antibonding) != inx + 1: - orb_anti += f"{name}, " - else: - orb_anti += f"and {name}" - - orb_anti += " orbitals, contributing " - for inx, contribution in enumerate(orb_antibonding): - if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti): - orb_anti += f"{contribution} " - elif 2 < len(orb_antibonding) != inx + 1: - orb_anti += f"{contribution}, " - else: - orb_anti += f"and {contribution} percent, respectively." - orb_info.append(f", whereas the maximum antibonding contribution is from {orb_anti}") - elif not orb_antibonding: - orb_info.append(", whereas no significant antibonding contribution is found in this bond.") - else: - orb_info.append( - f", whereas the maximum antibonding contribution is from the " - f"{orb_names_anti[0]} orbital, contributing {orb_antibonding[0]} percent." - ) - else: - # get atom-pair list with ion placed first - atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion) - percentage_cutoff = round(self.analysis_object.orbital_cutoff * 100, 2) - orb_info.append( - f"No individual orbital interactions detected above {percentage_cutoff} percent" - f" with summed I{type_pop} as reference for the " + "-".join(atom_pair) + " bond." - ) - - return orb_info - - def plot_cohps( - self, - xlim: list[float] | None = None, - ylim: list[float] | None = [-4, 2], - integrated: bool = False, - title: str = "", - save: bool = False, - filename: str | None = None, - sigma: float | None = None, - hide: bool = False, - ): - """ - Automatically generate plots of the most relevant COHPs, COOPs, or COBIs. - - :param save: will save the plot to a file - :param filename: name of the file to save the plot. - :param ylim: energy scale that is shown in plot (eV) - :param xlim: energy range for COHPs in eV - :param integrated: if True, integrated COHPs will be shown - :param sigma: Standard deviation of Gaussian broadening applied to - population data. If None, no broadening will be added. - :param title: sets the title of figure generated - :param hide: if True, the plot will not be shown. - - Returns: - A matplotlib object. - - """ - seq_cohps = self.analysis_object.seq_cohps - if self.analysis_object.which_bonds == "cation-anion": - seq_ineq_cations = self.analysis_object.seq_ineq_ions - elif self.analysis_object.which_bonds == "all": - seq_ineq_cations = self.analysis_object.seq_ineq_ions - seq_labels = self.analysis_object.seq_labels_cohps - structure = self.analysis_object.structure - - if len(seq_ineq_cations) >= 20: - warnings.warn( - "We will switch of displaying all plots " - "as there are more than 20 inequivalent ions. " - "We will instead save them in files called " - "'automatic-analysis-*.png'.", - stacklevel=2, - ) - hide = True - save = True - if filename is None: - filename = "./automatic_analysis.png" - - for iplot, (ication, labels, cohps) in enumerate(zip(seq_ineq_cations, seq_labels, seq_cohps)): - namecation = str(structure[ication].specie) - - cp = PlainCohpPlotter( - are_coops=self.analysis_object.are_coops, - are_cobis=self.analysis_object.are_cobis, - ) - for label, cohp in zip(labels, cohps): - if label is not None: - cp.add_cohp(namecation + str(ication + 1) + ": " + label, cohp) - - plot = cp.get_plot(integrated=integrated, sigma=sigma) - plot.ylim(ylim) - if xlim is not None: - plot.xlim(xlim) - - plot.title(title) - if save: - if len(seq_ineq_cations) > 1: - if isinstance(filename, str): - filename = Path(filename) # type: ignore - filename_new = ( - filename.parent / f"{filename.stem}-{iplot}{filename.suffix}" # type: ignore - ) - else: - filename_new = filename - plot.savefig(filename_new) - if hide: - plot.close() - if not hide: - plot.show() - - def plot_interactive_cohps( - self, - ylim: list[float] | None = None, - xlim: list[float] | None = None, - save_as_html: bool = False, - filename: str | None = None, - integrated: bool = False, - title: str = "", - sigma: float | None = None, - label_resolved: bool = False, - orbital_resolved: bool = False, - hide: bool = False, - ): - """ - Automatically generate interactive plots of the most relevant COHPs, COBIs or COOPs. - - :param save_as_html: will save the plot to a html file - :param filename: name of the file to save the plot. - :param ylim: energy scale that is shown in plot (eV) - :param xlim: energy range for COHPs in eV - :param integrated: if True, integrated COHPs will be shown - :param sigma: Standard deviation of Gaussian broadening applied to - population data. If None, no broadening will be added. - :param label_resolved: if true, relevant cohp curves will be further resolved based on band labels - :param orbital_resolved: if true, relevant orbital interactions in cohp curves will be added to figure - :param title: Title of the interactive plot - :param hide: if True, the plot will not be shown. - - Returns: - A plotly.graph_objects.Figure object. - """ - cba_cohp_plot_data = {} # Initialize dict to store plot data - set_cohps = self.analysis_object.seq_cohps - set_labels_cohps = self.analysis_object.seq_labels_cohps - set_inequivalent_cations = self.analysis_object.seq_ineq_ions - structure = self.analysis_object.structure - - for _iplot, (ication, labels, cohps) in enumerate(zip(set_inequivalent_cations, set_labels_cohps, set_cohps)): - label_str = f"{structure[ication].specie!s}{ication + 1!s}: " - for label, cohp in zip(labels, cohps): - if label is not None: - cba_cohp_plot_data[label_str + label] = cohp - - ip = InteractiveCohpPlotter( - are_coops=self.analysis_object.are_coops, - are_cobis=self.analysis_object.are_cobis, - ) - if label_resolved or orbital_resolved: - ip.add_all_relevant_cohps( - analyse=self.analysis_object, - label_resolved=label_resolved, - orbital_resolved=orbital_resolved, - ) - else: - ip.add_cohps_from_plot_data(plot_data_dict=cba_cohp_plot_data) - - plot = ip.get_plot(integrated=integrated, xlim=xlim, ylim=ylim, sigma=sigma) - - plot.update_layout(title_text=title) - if save_as_html: - plot.write_html(filename, include_mathjax="cdn") - if not hide: - return plot.show() - - return plot - - @staticmethod - def _coordination_environment_to_text(ce: str): - """ - Convert a coordination environment string into a text description of the environment. - - :param ce: output from ChemEnv package (e.g., "O:6") - - Returns: - A text description of coordination environment - """ - if ce == "S:1": - return "single (CN=1)" - if ce == "L:2": - return "linear (CN=2)" - if ce == "A:2": - return "angular (CN=2)" - if ce == "TL:3": - return "trigonal planar (CN=3)" - if ce == "TY:3": - return "triangular non-coplanar (CN=3)" - if ce == "TS:3": - return "t-shaped (CN=3)" - if ce == "T:4": - return "tetrahedral (CN=4)" - if ce == "S:4": - return "square planar (CN=4)" - if ce == "SY:4": - return "square non-coplanar (CN=4)" - if ce == "SS:4": - return "see-saw like (CN=4)" - if ce == "PP:5": - return "pentagonal (CN=5)" - if ce == "S:5": - return "square pyramidal (CN=5)" - if ce == "T:5": - return "trigonal bipyramidal (CN=5)" - if ce == "O:6": - return "octahedral (CN=6)" - if ce == "T:6": - return "trigonal prismatic (CN=6)" - if ce == "PP:6": - return "pentagonal pyramidal (CN=6)" - if ce == "PB:7": - return "pentagonal bipyramidal (CN=7)" - if ce == "ST:7": - return "square-face capped trigonal prismatic (CN=7)" - if ce == "ET:7": - return "end-trigonal-face capped trigonal prismatic (CN=7)" - if ce == "FO:7": - return "face-capped octahedron (CN=7)" - if ce == "C:8": - return "cubic (CN=8)" - if ce == "SA:8": - return "square antiprismatic (CN=8)" - if ce == "SBT:8": - return "square-face bicapped trigonal prismatic (CN=8)" - if ce == "TBT:8": - return "triangular-face bicapped trigonal prismatic (CN=8)" - if ce == "DD:8": - return "dodecahedronal (with triangular faces) (CN=8)" - if ce == "DDPN:8": - return "dodecahedronal (with triangular faces - p2345 plane normalized) (CN=8)" - if ce == "HB:8": - return "hexagonal bipyramidal (CN=8)" - if ce == "BO_1:8": - return "bicapped octahedral (opposed cap faces) (CN=8)" - if ce == "BO_2:8": - return "bicapped octahedral (cap faces with one atom in common) (CN=8)" - if ce == "BO_3:8": - return "bicapped octahedral (cap faces with one edge in common) (CN=8)" - if ce == "TC:9": - return "triangular cupola (CN=9)" - if ce == "TT_1:9": - return "Tricapped triangular prismatic (three square - face caps) (CN=9)" - if ce == "TT_2:9": - return "Tricapped triangular prismatic (two square - face caps and one triangular - face cap) (CN=9)" - if ce == "TT_3:9": - return "Tricapped triangular prism (one square - face cap and two triangular - face caps) (CN=9)" - if ce == "HD:9": - return "Heptagonal dipyramidal (CN=9)" - if ce == "TI:9": - return "tridiminished icosohedral (CN=9)" - if ce == "SMA:9": - return "Square-face monocapped antiprism (CN=9)" - if ce == "SS:9": - return "Square-face capped square prismatic (CN=9)" - if ce == "TO_1:9": - return "Tricapped octahedral (all 3 cap faces share one atom) (CN=9)" - if ce == "TO_2:9": - return "Tricapped octahedral (cap faces are aligned) (CN=9)" - if ce == "TO_3:9": - return "Tricapped octahedron (all 3 cap faces are sharing one edge of a face) (CN=9)" - if ce == "PP:10": - return "Pentagonal prismatic (CN=10)" - if ce == "PA:10": - return "Pentagonal antiprismatic (CN=10)" - if ce == "SBSA:10": - return "Square-face bicapped square antiprismatic (CN=10)" - if ce == "MI:10": - return "Metabidiminished icosahedral (CN=10)" - if ce == "S:10": - return "sphenocoronal (CN=10)" - if ce == "H:10": - return "Hexadecahedral (CN=10)" - if ce == "BS_1:10": - return "Bicapped square prismatic (opposite faces) (CN=10)" - if ce == "BS_1:10": - return "Bicapped square prismatic (opposite faces) (CN=10)" - if ce == "BS_2:10": - return "Bicapped square prism(adjacent faces) (CN=10)" - if ce == "TBSA:10": - return "Trigonal-face bicapped square antiprismatic (CN=10)" - if ce == "PCPA:11": - return "Pentagonal - face capped pentagonal antiprismatic (CN=11)" - if ce == "H:11": - return "Hendecahedral (CN=11)" - if ce == "SH:11": - return "Sphenoid hendecahedral (CN=11)" - if ce == "CO:11": - return "Cs - octahedral (CN=11)" - if ce == "DI:11": - return "Diminished icosahedral (CN=12)" - if ce == "I:12": - return "Icosahedral (CN=12)" - if ce == "PBP: 12": - return "Pentagonal - face bicapped pentagonal prismatic (CN=12)" - if ce == "TT:12": - return "Truncated tetrahedral (CN=12)" - if ce == "C:12": - return "Cuboctahedral (CN=12)" - if ce == "AC:12": - return "Anticuboctahedral (CN=12)" - if ce == "SC:12": - return "Square cupola (CN=12)" - if ce == "S:12": - return "Sphenomegacorona (CN=12)" - if ce == "HP:12": - return "Hexagonal prismatic (CN=12)" - if ce == "HA:12": - return "Hexagonal antiprismatic (CN=12)" - if ce == "SH:13": - return "Square-face capped hexagonal prismatic (CN=13)" - if ce == "1": - return "1-fold" - if ce == "2": - return "2-fold" - if ce == "3": - return "3-fold" - if ce == "4": - return "4-fold" - if ce == "5": - return "5-fold" - if ce == "6": - return "6-fold" - if ce == "7": - return "7-fold" - if ce == "8": - return "8-fold" - if ce == "9": - return "9-fold" - if ce == "10": - return "10-fold" - if ce == "11": - return "11-fold" - if ce == "12": - return "12-fold" - if ce == "13": - return "13-fold" - if ce == "14": - return "14-fold" - if ce == "15": - return "15-fold" - if ce == "16": - return "16-fold" - if ce == "17": - return "17-fold" - if ce == "18": - return "18-fold" - if ce == "19": - return "19-fold" - if ce == "20": - return "20-fold" - if ce == "21": - return "21-fold" - if ce == "22": - return "22-fold" - if ce == "23": - return "23-fold" - if ce == "24": - return "24-fold" - if ce == "25": - return "25-fold" - if ce == "26": - return "26-fold" - if ce == "27": - return "27-fold" - if ce == "28": - return "28-fold" - if ce == "29": - return "29-fold" - if ce == "30": - return "30-fold" - return ce - - def write_description(self): - """Print the description of the COHPs or COBIs or COOPs to the screen.""" - for textpart in self.text: - print(textpart) - - @staticmethod - def get_calc_quality_description(quality_dict): - """ - Generate a text description of the LOBSTER calculation quality. - - :param quality_dict: python dictionary from lobsterpy.analysis.get_lobster_calc_quality_summary - """ - text_des = [] - - for key, val in quality_dict.items(): - if key == "minimal_basis": - if val: - text_des.append("The LOBSTER calculation used minimal basis.") - if not val: - text_des.append( - "Consider rerunning the calculation with the minimum basis as well. Choosing a " - "larger basis set is only recommended if you see a significant improvement of " - "the charge spilling." - ) - - elif key == "charge_spilling": - text_des.append( - "The absolute and total charge spilling for the calculation is {} and {} %, " - "respectively.".format( - quality_dict[key]["abs_charge_spilling"], - quality_dict[key]["abs_total_spilling"], - ) - ) - elif key == "band_overlaps_analysis": - if quality_dict[key]["file_exists"]: - if quality_dict[key]["has_good_quality_maxDeviation"]: - text_des.append( - "The bandOverlaps.lobster file is generated during the LOBSTER run. This " - "indicates that the projected wave function is not completely orthonormalized; " - "however, the maximal deviation values observed compared to the identity matrix " - "is below the threshold of 0.1." - ) - else: - text_des.append( - "The bandOverlaps.lobster file is generated during the LOBSTER run. This " - "indicates that the projected wave function is not completely orthonormalized. " - "The maximal deviation value from the identity matrix is {}, and there are " - "{} percent k-points above the deviation threshold of 0.1. Please check the " - "results of other quality checks like dos comparisons, charges, " - "charge spillings before using the results for further " - "analysis.".format( - quality_dict[key]["max_deviation"], - quality_dict[key]["percent_kpoints_abv_limit"], - ) - ) - else: - text_des.append( - "The projected wave function is completely orthonormalized as no " - "bandOverlaps.lobster file is generated during the LOBSTER run." - ) - - elif key == "charge_comparisons": - if val: - for charge in ["mulliken", "loewdin"]: - if val[f"bva_{charge}_agree"]: - text_des.append( - f"The atomic charge signs from {charge.capitalize()} population analysis " - f"agree with the bond valence analysis." - ) - if not val[f"bva_{charge}_agree"]: - text_des.append( - f"The atomic charge signs from {charge.capitalize()} population analysis " - f"do not agree with the bond valence analysis." - ) - else: - text_des.append( - "Oxidation states from BVA analyzer cannot be determined. " - "Thus BVA charge comparison is not conducted." - ) - - elif key == "dos_comparisons": - comp_types = [] - tani_index = [] - for orb in val: - if orb.split("_")[-1] in ["s", "p", "d", "f", "summed"]: - comp_types.append(orb.split("_")[-1]) - tani_index.append(str(val[orb])) - text_des.append( - "The Tanimoto index from DOS comparisons in the energy range between {}, {} eV " - "for {} orbitals are: {}.".format( - val["e_range"][0], - val["e_range"][1], - ", ".join(comp_types), - ", ".join(tani_index), - ) - ) - - return text_des - - @staticmethod - def write_calc_quality_description(calc_quality_text): - """Print the calculation quality description to the screen.""" - print(" ".join(calc_quality_text)) +# Copyright (c) lobsterpy development team +# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License + +"""This module defines classes to describe the COHPs automatically.""" + +from __future__ import annotations + +import warnings +from pathlib import Path + +from lobsterpy.plotting import InteractiveCohpPlotter, PlainCohpPlotter + + +class Description: + """ + Base class that will write generate a text description for all relevant bonds. + + It analyses all relevant coordination environments in the system based on electronic structure theory. + + """ + + def __init__(self, analysis_object): + """ + Generate a text description for all relevant bonds. + + :param analysis_object: Analysis object from lobsterpy.analysis + """ + self.analysis_object = analysis_object + self.set_description() + + def set_description(self): + """ + Set the descriptions of the structures using the cation names, starting with numbers at 1. + + Uses the cation names from the lobster files. + + Returns: + None + + """ + self.condensed_bonding_analysis = self.analysis_object.condensed_bonding_analysis + # set type of population analyzed + type_pop = self.analysis_object._get_pop_type() + # set units for populations + units = " eV" if type_pop == "COHP" else "" + if self.analysis_object.which_bonds == "cation-anion": + relevant_cations = ", ".join( + [ + str(site.specie) + str(isite + 1) + for isite, site in enumerate(self.analysis_object.structure) + if isite in self.analysis_object.seq_ineq_ions + ] + ) + self.text = [] + self.text.append( + "The compound " + + str(self.condensed_bonding_analysis["formula"]) + + " has " + + str(self.condensed_bonding_analysis["number_of_considered_ions"]) + + " symmetry-independent cation(s) with relevant cation-anion interactions: " + + str(relevant_cations) + + "." + ) + + for key, item in self.condensed_bonding_analysis["sites"].items(): + # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi), + bond_info = [] + orb_info = [] + for type, properties in item["bonds"].items(): + if not properties["has_antibdg_states_below_Efermi"]: + bond_info.append( + str(properties["number_of_bonds"]) + + " " + + item["ion"] + + "-" + + str(type) + + f" (mean I{type_pop}: " + "" + + properties[f"I{type_pop}_mean"] + + f"{units}, 0.0 percent antibonding interaction below EFermi)" + ) + if self.analysis_object.orbital_resolved: + text_orbital = self._generate_orbital_resolved_analysis_text( + orbital_resolved_data=properties, + type_pop=type_pop, + atom_name=str(type), + ion=item["ion"], + ) + orb_info.extend(text_orbital) + else: + bond_info.append( + str(properties["number_of_bonds"]) + + " " + + item["ion"] + + "-" + + str(type) + + f" (mean I{type_pop}: " + "" + + properties[f"I{type_pop}_mean"] + + f"{units}, " + + str(round(properties["antibonding"]["perc"] * 100, 3)) + + " percent antibonding interaction below EFermi)" + ) + if self.analysis_object.orbital_resolved: + text_orbital = self._generate_orbital_resolved_analysis_text( + orbital_resolved_data=properties, + type_pop=type_pop, + atom_name=str(type), + ion=item["ion"], + ) + orb_info.extend(text_orbital) + + bonds = ",".join(bond_info[0:-1]) + ", and " + bond_info[-1] if len(bond_info) > 1 else bond_info[0] + + if len(orb_info) > 1: + orb_bonds = "".join(orb_info).replace(".In", ". In") + else: + orb_bonds = orb_info[0] if orb_info else "" + if item["env"] == "O:6": + self.text.append( + str(item["ion"]) + + str(key + 1) + + " has an " + + str(self._coordination_environment_to_text(item["env"])) + + " coordination environment. It has " + + str(bonds) + + " bonds." + ) + if orb_bonds: + self.text.append(orb_bonds) + else: + self.text.append( + str(item["ion"]) + + str(key + 1) + + " has a " + + str(self._coordination_environment_to_text(item["env"])) + + " coordination environment. It has " + + str(bonds) + + " bonds." + ) + if orb_bonds: + self.text.append(orb_bonds) + + elif self.analysis_object.which_bonds == "all": + relevant_ions = ", ".join( + [ + str(site.specie) + str(isite + 1) + for isite, site in enumerate(self.analysis_object.structure) + if isite in self.analysis_object.seq_ineq_ions + ] + ) + self.text = [] + self.text.append( + "The compound " + + str(self.condensed_bonding_analysis["formula"]) + + " has " + + str(self.condensed_bonding_analysis["number_of_considered_ions"]) + + " symmetry-independent atoms(s) with relevant bonds: " + + str(relevant_ions) + + "." + ) + + for key, item in self.condensed_bonding_analysis["sites"].items(): + # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi), + bond_info = [] + orb_info = [] + for type, properties in item["bonds"].items(): + if not properties["has_antibdg_states_below_Efermi"]: + bond_info.append( + str(properties["number_of_bonds"]) + + " " + + item["ion"] + + "-" + + str(type) + + f" (mean I{type_pop}: " + "" + + properties[f"I{type_pop}_mean"] + + f"{units}, 0.0 percent antibonding interaction below EFermi)" + ) + if self.analysis_object.orbital_resolved: + text_orbital = self._generate_orbital_resolved_analysis_text( + orbital_resolved_data=properties, + type_pop=type_pop, + atom_name=str(type), + ion=item["ion"], + ) + orb_info.extend(text_orbital) + else: + bond_info.append( + str(properties["number_of_bonds"]) + + " " + + item["ion"] + + "-" + + str(type) + + f" (mean I{type_pop}: " + "" + + properties[f"I{type_pop}_mean"] + + f"{units}, " + + str(round(properties["antibonding"]["perc"] * 100, 3)) + + " percent antibonding interaction below EFermi)" + ) + + if self.analysis_object.orbital_resolved: + text_orbital = self._generate_orbital_resolved_analysis_text( + orbital_resolved_data=properties, + type_pop=type_pop, + atom_name=str(type), + ion=item["ion"], + ) + orb_info.extend(text_orbital) + + bonds = ",".join(bond_info[0:-1]) + ", and " + bond_info[-1] if len(bond_info) > 1 else bond_info[0] + + if len(orb_info) > 1: + orb_bonds = "".join(orb_info).replace(".In", ". In") + else: + orb_bonds = orb_info[0] if orb_info else "" + if item["env"] == "O:6": + self.text.append( + str(item["ion"]) + + str(key + 1) + + " has an " + + str(self._coordination_environment_to_text(item["env"])) + + " coordination environment. It has " + + str(bonds) + + " bonds." + ) + if orb_bonds: + self.text.append(orb_bonds) + else: + self.text.append( + str(item["ion"]) + + str(key + 1) + + " has a " + + str(self._coordination_environment_to_text(item["env"])) + + " coordination environment. It has " + + str(bonds) + + " bonds." + ) + if orb_bonds: + self.text.append(orb_bonds) + + if ( + "madelung_energy" in self.analysis_object.condensed_bonding_analysis + and self.analysis_object.condensed_bonding_analysis["madelung_energy"] is not None + ): + self.text.append( + "The Madelung energy of this crystal structure per unit cell is: " + + str(self.analysis_object.condensed_bonding_analysis["madelung_energy"]) + + " eV." + ) + + def _generate_orbital_resolved_analysis_text( + self, + orbital_resolved_data: dict, + ion: str, + atom_name: str, + type_pop: str, + ): + """ + Generate text from orbital-resolved analysis data of the most relevant COHP, COOP, or COBI. + + :param orbital_resolved_data: dict of orbital data from condensed bonding analysis object + :param ion: name of ion at the site + :param atom_name: name of atomic speice to which ion is bonded + :param type_pop: population type analysed could be "COHP" or "COOP" or "COBI" + + Returns: + A python list with text describing the orbital which contributes + the most to the bonding and antibonding in the bond at site + """ + orb_info = [] + if orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: + orb_names = [] + orb_contri = [] + # get atom-pair list with ion placed first + atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion) + if "max_bonding_contribution" in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: + for orb, data in orbital_resolved_data["orbital_data"]["orbital_summary_stats"][ + "max_bonding_contribution" + ].items(): + atom_pair_with_orb_name = self.analysis_object._sort_orbital_atom_pair( + atom_pair=atom_pair, + complete_cohp=self.analysis_object.chemenv.completecohp, + label=orbital_resolved_data["orbital_data"]["relevant_bonds"][0], + orb_pair=orb, + ) + orb_names.append("-".join(atom_pair_with_orb_name)) + orb_contri.append( + str( + round( + data * 100, + 3, + ) + ) + ) + orb_names_anti = [] + orb_antibonding = [] + if "max_antibonding_contribution" in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]: + for orb, data in orbital_resolved_data["orbital_data"]["orbital_summary_stats"][ + "max_antibonding_contribution" + ].items(): + atom_pair_with_orb_name = self.analysis_object._sort_orbital_atom_pair( + atom_pair=atom_pair, + complete_cohp=self.analysis_object.chemenv.completecohp, + label=orbital_resolved_data["orbital_data"]["relevant_bonds"][0], + orb_pair=orb, + ) + orb_names_anti.append("-".join(atom_pair_with_orb_name)) + orb_antibonding.append( + str( + round( + data * 100, + 3, + ) + ) + ) + if len(orb_contri) > 1: + orb_name_contri = "" + for inx, name in enumerate(orb_names): + if len(orb_contri) == 2 and inx + 1 != len(orb_contri): + orb_name_contri += f"{name} " + elif 2 < len(orb_contri) != inx + 1: + orb_name_contri += f"{name}, " + else: + orb_name_contri += f"and {name}" + + orb_name_contri += " orbitals, contributing " + for inx, contribution in enumerate(orb_contri): + if len(orb_contri) == 2 and inx + 1 != len(orb_contri): + orb_name_contri += f"{contribution} " + elif 2 < len(orb_contri) != inx + 1: + orb_name_contri += f"{contribution}, " + else: + orb_name_contri += f"and {contribution} percent, respectively" + num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) + bonds = "bonds" if num_bonds > 1 else "bond" + orb_info.append( + f"In the {num_bonds} " + + "-".join(atom_pair) + + f" {bonds}, relative to the summed I{type_pop}s, " + + "the maximum bonding contribution is from " + + orb_name_contri + ) + elif not orb_contri: + num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) + bonds = "bonds" if num_bonds > 1 else "bond" + orb_info.append( + f"In the {num_bonds} " + + "-".join(atom_pair) + + f" {bonds}, relative to the summed I{type_pop}s, " + + f"no orbital has a bonding contribution greater than " + f"{self.analysis_object.orbital_cutoff*100} percent" + ) + else: + num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"]) + bonds = "bonds" if num_bonds > 1 else "bond" + orb_info.append( + f"In the {num_bonds} " + + "-".join(atom_pair) + + f" {bonds}, relative to the summed I{type_pop}s, " + + "the maximum bonding contribution is from the " + + f"{orb_names[0]}" + + f" orbital, contributing {orb_contri[0]} percent" + ) + + if len(orb_antibonding) > 1: + orb_anti = "" + for inx, name in enumerate(orb_names_anti): + if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti): + orb_anti += f"{name} " + elif 2 < len(orb_antibonding) != inx + 1: + orb_anti += f"{name}, " + else: + orb_anti += f"and {name}" + + orb_anti += " orbitals, contributing " + for inx, contribution in enumerate(orb_antibonding): + if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti): + orb_anti += f"{contribution} " + elif 2 < len(orb_antibonding) != inx + 1: + orb_anti += f"{contribution}, " + else: + orb_anti += f"and {contribution} percent, respectively." + orb_info.append(f", whereas the maximum antibonding contribution is from {orb_anti}") + elif not orb_antibonding: + orb_info.append(", whereas no significant antibonding contribution is found in this bond.") + else: + orb_info.append( + f", whereas the maximum antibonding contribution is from the " + f"{orb_names_anti[0]} orbital, contributing {orb_antibonding[0]} percent." + ) + else: + # get atom-pair list with ion placed first + atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion) + percentage_cutoff = round(self.analysis_object.orbital_cutoff * 100, 2) + orb_info.append( + f"No individual orbital interactions detected above {percentage_cutoff} percent" + f" with summed I{type_pop} as reference for the " + "-".join(atom_pair) + " bond." + ) + + return orb_info + + def plot_cohps( + self, + xlim: list[float] | None = None, + ylim: list[float] | None = [-4, 2], + integrated: bool = False, + title: str = "", + save: bool = False, + filename: str | None = None, + sigma: float | None = None, + hide: bool = False, + ): + """ + Automatically generate plots of the most relevant COHPs, COOPs, or COBIs. + + :param save: will save the plot to a file + :param filename: name of the file to save the plot. + :param ylim: energy scale that is shown in plot (eV) + :param xlim: energy range for COHPs in eV + :param integrated: if True, integrated COHPs will be shown + :param sigma: Standard deviation of Gaussian broadening applied to + population data. If None, no broadening will be added. + :param title: sets the title of figure generated + :param hide: if True, the plot will not be shown. + + Returns: + A matplotlib object. + + """ + seq_cohps = self.analysis_object.seq_cohps + if self.analysis_object.which_bonds == "cation-anion": + seq_ineq_cations = self.analysis_object.seq_ineq_ions + elif self.analysis_object.which_bonds == "all": + seq_ineq_cations = self.analysis_object.seq_ineq_ions + seq_labels = self.analysis_object.seq_labels_cohps + structure = self.analysis_object.structure + + if len(seq_ineq_cations) >= 20: + warnings.warn( + "We will switch of displaying all plots " + "as there are more than 20 inequivalent ions. " + "We will instead save them in files called " + "'automatic-analysis-*.png'.", + stacklevel=2, + ) + hide = True + save = True + if filename is None: + filename = "./automatic_analysis.png" + + for iplot, (ication, labels, cohps) in enumerate(zip(seq_ineq_cations, seq_labels, seq_cohps)): + namecation = str(structure[ication].specie) + + cp = PlainCohpPlotter( + are_coops=self.analysis_object.are_coops, + are_cobis=self.analysis_object.are_cobis, + ) + for label, cohp in zip(labels, cohps): + if label is not None: + cp.add_cohp(namecation + str(ication + 1) + ": " + label, cohp) + + plot = cp.get_plot(integrated=integrated, sigma=sigma) + plot.ylim(ylim) + if xlim is not None: + plot.xlim(xlim) + + plot.title(title) + if save: + if len(seq_ineq_cations) > 1: + if isinstance(filename, str): + filename = Path(filename) # type: ignore + filename_new = ( + filename.parent / f"{filename.stem}-{iplot}{filename.suffix}" # type: ignore + ) + else: + filename_new = filename + plot.savefig(filename_new) + if hide: + plot.close() + if not hide: + plot.show() + + def plot_interactive_cohps( + self, + ylim: list[float] | None = None, + xlim: list[float] | None = None, + save_as_html: bool = False, + filename: str | None = None, + integrated: bool = False, + title: str = "", + sigma: float | None = None, + label_resolved: bool = False, + orbital_resolved: bool = False, + hide: bool = False, + ): + """ + Automatically generate interactive plots of the most relevant COHPs, COBIs or COOPs. + + :param save_as_html: will save the plot to a html file + :param filename: name of the file to save the plot. + :param ylim: energy scale that is shown in plot (eV) + :param xlim: energy range for COHPs in eV + :param integrated: if True, integrated COHPs will be shown + :param sigma: Standard deviation of Gaussian broadening applied to + population data. If None, no broadening will be added. + :param label_resolved: if true, relevant cohp curves will be further resolved based on band labels + :param orbital_resolved: if true, relevant orbital interactions in cohp curves will be added to figure + :param title: Title of the interactive plot + :param hide: if True, the plot will not be shown. + + Returns: + A plotly.graph_objects.Figure object. + """ + cba_cohp_plot_data = {} # Initialize dict to store plot data + set_cohps = self.analysis_object.seq_cohps + set_labels_cohps = self.analysis_object.seq_labels_cohps + set_inequivalent_cations = self.analysis_object.seq_ineq_ions + structure = self.analysis_object.structure + + for _iplot, (ication, labels, cohps) in enumerate(zip(set_inequivalent_cations, set_labels_cohps, set_cohps)): + label_str = f"{structure[ication].specie!s}{ication + 1!s}: " + for label, cohp in zip(labels, cohps): + if label is not None: + cba_cohp_plot_data[label_str + label] = cohp + + ip = InteractiveCohpPlotter( + are_coops=self.analysis_object.are_coops, + are_cobis=self.analysis_object.are_cobis, + ) + if label_resolved or orbital_resolved: + ip.add_all_relevant_cohps( + analyse=self.analysis_object, + label_resolved=label_resolved, + orbital_resolved=orbital_resolved, + ) + else: + ip.add_cohps_from_plot_data(plot_data_dict=cba_cohp_plot_data) + + plot = ip.get_plot(integrated=integrated, xlim=xlim, ylim=ylim, sigma=sigma) + + plot.update_layout(title_text=title) + if save_as_html: + plot.write_html(filename, include_mathjax="cdn") + if not hide: + return plot.show() + + return plot + + @staticmethod + def _coordination_environment_to_text(ce: str): + """ + Convert a coordination environment string into a text description of the environment. + + :param ce: output from ChemEnv package (e.g., "O:6") + + Returns: + A text description of coordination environment + """ + if ce == "S:1": + return "single (CN=1)" + if ce == "L:2": + return "linear (CN=2)" + if ce == "A:2": + return "angular (CN=2)" + if ce == "TL:3": + return "trigonal planar (CN=3)" + if ce == "TY:3": + return "triangular non-coplanar (CN=3)" + if ce == "TS:3": + return "t-shaped (CN=3)" + if ce == "T:4": + return "tetrahedral (CN=4)" + if ce == "S:4": + return "square planar (CN=4)" + if ce == "SY:4": + return "square non-coplanar (CN=4)" + if ce == "SS:4": + return "see-saw like (CN=4)" + if ce == "PP:5": + return "pentagonal (CN=5)" + if ce == "S:5": + return "square pyramidal (CN=5)" + if ce == "T:5": + return "trigonal bipyramidal (CN=5)" + if ce == "O:6": + return "octahedral (CN=6)" + if ce == "T:6": + return "trigonal prismatic (CN=6)" + if ce == "PP:6": + return "pentagonal pyramidal (CN=6)" + if ce == "PB:7": + return "pentagonal bipyramidal (CN=7)" + if ce == "ST:7": + return "square-face capped trigonal prismatic (CN=7)" + if ce == "ET:7": + return "end-trigonal-face capped trigonal prismatic (CN=7)" + if ce == "FO:7": + return "face-capped octahedron (CN=7)" + if ce == "C:8": + return "cubic (CN=8)" + if ce == "SA:8": + return "square antiprismatic (CN=8)" + if ce == "SBT:8": + return "square-face bicapped trigonal prismatic (CN=8)" + if ce == "TBT:8": + return "triangular-face bicapped trigonal prismatic (CN=8)" + if ce == "DD:8": + return "dodecahedronal (with triangular faces) (CN=8)" + if ce == "DDPN:8": + return "dodecahedronal (with triangular faces - p2345 plane normalized) (CN=8)" + if ce == "HB:8": + return "hexagonal bipyramidal (CN=8)" + if ce == "BO_1:8": + return "bicapped octahedral (opposed cap faces) (CN=8)" + if ce == "BO_2:8": + return "bicapped octahedral (cap faces with one atom in common) (CN=8)" + if ce == "BO_3:8": + return "bicapped octahedral (cap faces with one edge in common) (CN=8)" + if ce == "TC:9": + return "triangular cupola (CN=9)" + if ce == "TT_1:9": + return "Tricapped triangular prismatic (three square - face caps) (CN=9)" + if ce == "TT_2:9": + return "Tricapped triangular prismatic (two square - face caps and one triangular - face cap) (CN=9)" + if ce == "TT_3:9": + return "Tricapped triangular prism (one square - face cap and two triangular - face caps) (CN=9)" + if ce == "HD:9": + return "Heptagonal dipyramidal (CN=9)" + if ce == "TI:9": + return "tridiminished icosohedral (CN=9)" + if ce == "SMA:9": + return "Square-face monocapped antiprism (CN=9)" + if ce == "SS:9": + return "Square-face capped square prismatic (CN=9)" + if ce == "TO_1:9": + return "Tricapped octahedral (all 3 cap faces share one atom) (CN=9)" + if ce == "TO_2:9": + return "Tricapped octahedral (cap faces are aligned) (CN=9)" + if ce == "TO_3:9": + return "Tricapped octahedron (all 3 cap faces are sharing one edge of a face) (CN=9)" + if ce == "PP:10": + return "Pentagonal prismatic (CN=10)" + if ce == "PA:10": + return "Pentagonal antiprismatic (CN=10)" + if ce == "SBSA:10": + return "Square-face bicapped square antiprismatic (CN=10)" + if ce == "MI:10": + return "Metabidiminished icosahedral (CN=10)" + if ce == "S:10": + return "sphenocoronal (CN=10)" + if ce == "H:10": + return "Hexadecahedral (CN=10)" + if ce == "BS_1:10": + return "Bicapped square prismatic (opposite faces) (CN=10)" + if ce == "BS_1:10": + return "Bicapped square prismatic (opposite faces) (CN=10)" + if ce == "BS_2:10": + return "Bicapped square prism(adjacent faces) (CN=10)" + if ce == "TBSA:10": + return "Trigonal-face bicapped square antiprismatic (CN=10)" + if ce == "PCPA:11": + return "Pentagonal - face capped pentagonal antiprismatic (CN=11)" + if ce == "H:11": + return "Hendecahedral (CN=11)" + if ce == "SH:11": + return "Sphenoid hendecahedral (CN=11)" + if ce == "CO:11": + return "Cs - octahedral (CN=11)" + if ce == "DI:11": + return "Diminished icosahedral (CN=12)" + if ce == "I:12": + return "Icosahedral (CN=12)" + if ce == "PBP: 12": + return "Pentagonal - face bicapped pentagonal prismatic (CN=12)" + if ce == "TT:12": + return "Truncated tetrahedral (CN=12)" + if ce == "C:12": + return "Cuboctahedral (CN=12)" + if ce == "AC:12": + return "Anticuboctahedral (CN=12)" + if ce == "SC:12": + return "Square cupola (CN=12)" + if ce == "S:12": + return "Sphenomegacorona (CN=12)" + if ce == "HP:12": + return "Hexagonal prismatic (CN=12)" + if ce == "HA:12": + return "Hexagonal antiprismatic (CN=12)" + if ce == "SH:13": + return "Square-face capped hexagonal prismatic (CN=13)" + if ce == "1": + return "1-fold" + if ce == "2": + return "2-fold" + if ce == "3": + return "3-fold" + if ce == "4": + return "4-fold" + if ce == "5": + return "5-fold" + if ce == "6": + return "6-fold" + if ce == "7": + return "7-fold" + if ce == "8": + return "8-fold" + if ce == "9": + return "9-fold" + if ce == "10": + return "10-fold" + if ce == "11": + return "11-fold" + if ce == "12": + return "12-fold" + if ce == "13": + return "13-fold" + if ce == "14": + return "14-fold" + if ce == "15": + return "15-fold" + if ce == "16": + return "16-fold" + if ce == "17": + return "17-fold" + if ce == "18": + return "18-fold" + if ce == "19": + return "19-fold" + if ce == "20": + return "20-fold" + if ce == "21": + return "21-fold" + if ce == "22": + return "22-fold" + if ce == "23": + return "23-fold" + if ce == "24": + return "24-fold" + if ce == "25": + return "25-fold" + if ce == "26": + return "26-fold" + if ce == "27": + return "27-fold" + if ce == "28": + return "28-fold" + if ce == "29": + return "29-fold" + if ce == "30": + return "30-fold" + return ce + + def write_description(self): + """Print the description of the COHPs or COBIs or COOPs to the screen.""" + for textpart in self.text: + print(textpart) + + @staticmethod + def get_calc_quality_description(quality_dict): + """ + Generate a text description of the LOBSTER calculation quality. + + :param quality_dict: python dictionary from lobsterpy.analysis.get_lobster_calc_quality_summary + """ + text_des = [] + + for key, val in quality_dict.items(): + if key == "minimal_basis": + if val: + text_des.append("The LOBSTER calculation used minimal basis.") + if not val: + text_des.append( + "Consider rerunning the calculation with the minimum basis as well. Choosing a " + "larger basis set is only recommended if you see a significant improvement of " + "the charge spilling." + ) + + elif key == "charge_spilling": + text_des.append( + "The absolute and total charge spilling for the calculation is {} and {} %, " + "respectively.".format( + quality_dict[key]["abs_charge_spilling"], + quality_dict[key]["abs_total_spilling"], + ) + ) + elif key == "band_overlaps_analysis": + if quality_dict[key]["file_exists"]: + if quality_dict[key]["has_good_quality_maxDeviation"]: + text_des.append( + "The bandOverlaps.lobster file is generated during the LOBSTER run. This " + "indicates that the projected wave function is not completely orthonormalized; " + "however, the maximal deviation values observed compared to the identity matrix " + "is below the threshold of 0.1." + ) + else: + text_des.append( + "The bandOverlaps.lobster file is generated during the LOBSTER run. This " + "indicates that the projected wave function is not completely orthonormalized. " + "The maximal deviation value from the identity matrix is {}, and there are " + "{} percent k-points above the deviation threshold of 0.1. Please check the " + "results of other quality checks like dos comparisons, charges, " + "charge spillings before using the results for further " + "analysis.".format( + quality_dict[key]["max_deviation"], + quality_dict[key]["percent_kpoints_abv_limit"], + ) + ) + else: + text_des.append( + "The projected wave function is completely orthonormalized as no " + "bandOverlaps.lobster file is generated during the LOBSTER run." + ) + + elif key == "charge_comparisons": + if val: + for charge in ["mulliken", "loewdin"]: + if val[f"bva_{charge}_agree"]: + text_des.append( + f"The atomic charge signs from {charge.capitalize()} population analysis " + f"agree with the bond valence analysis." + ) + if not val[f"bva_{charge}_agree"]: + text_des.append( + f"The atomic charge signs from {charge.capitalize()} population analysis " + f"do not agree with the bond valence analysis." + ) + else: + text_des.append( + "Oxidation states from BVA analyzer cannot be determined. " + "Thus BVA charge comparison is not conducted." + ) + + elif key == "dos_comparisons": + comp_types = [] + tani_index = [] + for orb in val: + if orb.split("_")[-1] in ["s", "p", "d", "f", "summed"]: + comp_types.append(orb.split("_")[-1]) + tani_index.append(str(val[orb])) + text_des.append( + "The Tanimoto index from DOS comparisons in the energy range between {}, {} eV " + "for {} orbitals are: {}.".format( + val["e_range"][0], + val["e_range"][1], + ", ".join(comp_types), + ", ".join(tani_index), + ) + ) + + return text_des + + @staticmethod + def write_calc_quality_description(calc_quality_text): + """Print the calculation quality description to the screen.""" + print(" ".join(calc_quality_text)) diff --git a/lobsterpy/featurize/__init__.py b/src/lobsterpy/featurize/__init__.py similarity index 100% rename from lobsterpy/featurize/__init__.py rename to src/lobsterpy/featurize/__init__.py diff --git a/lobsterpy/featurize/batch.py b/src/lobsterpy/featurize/batch.py similarity index 100% rename from lobsterpy/featurize/batch.py rename to src/lobsterpy/featurize/batch.py diff --git a/lobsterpy/featurize/core.py b/src/lobsterpy/featurize/core.py similarity index 100% rename from lobsterpy/featurize/core.py rename to src/lobsterpy/featurize/core.py diff --git a/lobsterpy/featurize/utils.py b/src/lobsterpy/featurize/utils.py similarity index 100% rename from lobsterpy/featurize/utils.py rename to src/lobsterpy/featurize/utils.py diff --git a/lobsterpy/plotting/__init__.py b/src/lobsterpy/plotting/__init__.py similarity index 100% rename from lobsterpy/plotting/__init__.py rename to src/lobsterpy/plotting/__init__.py diff --git a/lobsterpy/plotting/layout_dicts.py b/src/lobsterpy/plotting/layout_dicts.py similarity index 96% rename from lobsterpy/plotting/layout_dicts.py rename to src/lobsterpy/plotting/layout_dicts.py index 106d6a25..197670c6 100644 --- a/lobsterpy/plotting/layout_dicts.py +++ b/src/lobsterpy/plotting/layout_dicts.py @@ -1,141 +1,141 @@ -""" -Dictionaries for plotly figure layouts. - -Layout inspired by -https://github.com/materialsproject/crystaltoolkit/blob/main/crystal_toolkit/components/bandstructure.py - -Contains dicts for - 1. general figure layout - 2. axes ((I)COHP / energy) - 3. traces (Spin.up / Spin.down) - 4. legend - -""" - -layout_dict = { - "titlefont": {"size": 18}, - "showlegend": True, - "title_x": 0.1, - "title_y": 0.9, - "hovermode": "closest", - "paper_bgcolor": "rgba(0,0,0,0)", - "plot_bgcolor": "rgba(0,0,0,0)", -} - -""" - General layout of Plotly figure. - - :param dict titlefont: Line style dictionary. - :param bool showlegend: Legend hide or show. - :param float title_x: x axis title size. - :param float title_y: y axis title size. - :param str hovermode: hover behaviour. - :param str paper_bgcolor: background color. - :param str plot_bgcolor: plot background color. -""" - -cohp_axis_style_dict = { - "titlefont": {"size": 20}, - "tickfont": {"size": 16}, - "ticks": "inside", - "tickwidth": 2, - "showgrid": False, - "showline": True, - "zeroline": True, - "zerolinecolor": "black", - "zerolinewidth": 3, - "mirror": True, - "linewidth": 2, - "linecolor": "black", -} - -""" - COXX axis style. - - :param dict titlefont: axis title font style. - :param dict tickfont: axis tick font style. - :param str ticks: ticks style. - :param float tickwidth: width of ticks. - :param bool showgrid: show or hide axis grid. - :param bool showline: show or hide axis. - :param bool zeroline: show or hide zero line. - :param str zerolinecolor: color of zero line. - :param float zerolinewidth: width of zero line. - :param bool mirror: mirror axis. - :param float linewidth: width of axis line. - :param str linecolor: color of axis line. -""" - -energy_axis_style_dict = { - "titlefont": {"size": 20}, - "tickfont": {"size": 16}, - "ticks": "inside", - "tickwidth": 2, - "showgrid": False, - "showline": True, - "zeroline": True, - "zerolinecolor": "rgba(0,0,0,0.5)", - "zerolinewidth": 3, - "mirror": True, - "linewidth": 2, - "linecolor": "black", -} - -""" - Energy axis style. - - :param dict titlefont: axis title font style. - :param dict tickfont: axis tick font style. - :param str ticks: ticks style. - :param float tickwidth: width of ticks. - :param bool showgrid: show or hide axis grid. - :param bool showline: show or hide axis. - :param bool zeroline: show or hide zero line. - :param str zerolinecolor: color of zero line. - :param float zerolinewidth: width of zero line. - :param bool mirror: mirror axis. - :param float linewidth: width of axis line. - :param str linecolor: color of axis line. -""" - -spin_up_trace_style_dict = { - "line": {"width": 3}, - "mode": "lines", - "visible": "legendonly", -} - -""" - Line style for Spin.up traces. - - :param dict line: Line style dictionary. - :param str mode: Plotly mode. - :param str visible: Trace visibility. -""" - -spin_down_trace_style_dict = { - "line": {"width": 3, "dash": "dash"}, - "mode": "lines", - "visible": "legendonly", -} - -""" - Line style for Spin.down traces. - - :param dict line: Line style dictionary. - :param str mode: Plotly mode. - :param str visible: Trace visibility. -""" - -legend_style_dict = { - "bordercolor": "black", - "borderwidth": 2, - "traceorder": "normal", -} - -""" - Legend style. - - :param str bordercolor: border color of legend. - :param int borderwidth: width of border. - :param str traceorder: order of trace. -""" +""" +Dictionaries for plotly figure layouts. + +Layout inspired by +https://github.com/materialsproject/crystaltoolkit/blob/main/crystal_toolkit/components/bandstructure.py + +Contains dicts for + 1. general figure layout + 2. axes ((I)COHP / energy) + 3. traces (Spin.up / Spin.down) + 4. legend + +""" + +layout_dict = { + "titlefont": {"size": 18}, + "showlegend": True, + "title_x": 0.1, + "title_y": 0.9, + "hovermode": "closest", + "paper_bgcolor": "rgba(0,0,0,0)", + "plot_bgcolor": "rgba(0,0,0,0)", +} + +""" + General layout of Plotly figure. + + :param dict titlefont: Line style dictionary. + :param bool showlegend: Legend hide or show. + :param float title_x: x axis title size. + :param float title_y: y axis title size. + :param str hovermode: hover behaviour. + :param str paper_bgcolor: background color. + :param str plot_bgcolor: plot background color. +""" + +cohp_axis_style_dict = { + "titlefont": {"size": 20}, + "tickfont": {"size": 16}, + "ticks": "inside", + "tickwidth": 2, + "showgrid": False, + "showline": True, + "zeroline": True, + "zerolinecolor": "black", + "zerolinewidth": 3, + "mirror": True, + "linewidth": 2, + "linecolor": "black", +} + +""" + COXX axis style. + + :param dict titlefont: axis title font style. + :param dict tickfont: axis tick font style. + :param str ticks: ticks style. + :param float tickwidth: width of ticks. + :param bool showgrid: show or hide axis grid. + :param bool showline: show or hide axis. + :param bool zeroline: show or hide zero line. + :param str zerolinecolor: color of zero line. + :param float zerolinewidth: width of zero line. + :param bool mirror: mirror axis. + :param float linewidth: width of axis line. + :param str linecolor: color of axis line. +""" + +energy_axis_style_dict = { + "titlefont": {"size": 20}, + "tickfont": {"size": 16}, + "ticks": "inside", + "tickwidth": 2, + "showgrid": False, + "showline": True, + "zeroline": True, + "zerolinecolor": "rgba(0,0,0,0.5)", + "zerolinewidth": 3, + "mirror": True, + "linewidth": 2, + "linecolor": "black", +} + +""" + Energy axis style. + + :param dict titlefont: axis title font style. + :param dict tickfont: axis tick font style. + :param str ticks: ticks style. + :param float tickwidth: width of ticks. + :param bool showgrid: show or hide axis grid. + :param bool showline: show or hide axis. + :param bool zeroline: show or hide zero line. + :param str zerolinecolor: color of zero line. + :param float zerolinewidth: width of zero line. + :param bool mirror: mirror axis. + :param float linewidth: width of axis line. + :param str linecolor: color of axis line. +""" + +spin_up_trace_style_dict = { + "line": {"width": 3}, + "mode": "lines", + "visible": "legendonly", +} + +""" + Line style for Spin.up traces. + + :param dict line: Line style dictionary. + :param str mode: Plotly mode. + :param str visible: Trace visibility. +""" + +spin_down_trace_style_dict = { + "line": {"width": 3, "dash": "dash"}, + "mode": "lines", + "visible": "legendonly", +} + +""" + Line style for Spin.down traces. + + :param dict line: Line style dictionary. + :param str mode: Plotly mode. + :param str visible: Trace visibility. +""" + +legend_style_dict = { + "bordercolor": "black", + "borderwidth": 2, + "traceorder": "normal", +} + +""" + Legend style. + + :param str bordercolor: border color of legend. + :param int borderwidth: width of border. + :param str traceorder: order of trace. +""" diff --git a/lobsterpy/plotting/lobsterpy_base.mplstyle b/src/lobsterpy/plotting/lobsterpy_base.mplstyle similarity index 100% rename from lobsterpy/plotting/lobsterpy_base.mplstyle rename to src/lobsterpy/plotting/lobsterpy_base.mplstyle diff --git a/lobsterpy/structuregraph/__init__.py b/src/lobsterpy/structuregraph/__init__.py similarity index 100% rename from lobsterpy/structuregraph/__init__.py rename to src/lobsterpy/structuregraph/__init__.py diff --git a/lobsterpy/structuregraph/graph.py b/src/lobsterpy/structuregraph/graph.py similarity index 100% rename from lobsterpy/structuregraph/graph.py rename to src/lobsterpy/structuregraph/graph.py From 072c4e732499c25ed69f70556b49868d4ebd68d8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sat, 20 Jul 2024 13:24:37 +0000 Subject: [PATCH 2/3] pre-commit auto-fixes --- tests/cli/test_cli.py | 3 +-- tests/cohp/test_analyze.py | 3 +-- tests/conftest.py | 5 ++--- tests/featurize/test_batch.py | 3 +-- tests/featurize/test_core.py | 1 - tests/featurize/test_utils.py | 3 +-- tests/plotting/test_plotting.py | 7 +++---- tests/structuregraph/test_graph.py | 1 - 8 files changed, 9 insertions(+), 17 deletions(-) diff --git a/tests/cli/test_cli.py b/tests/cli/test_cli.py index e41bb021..420478fe 100644 --- a/tests/cli/test_cli.py +++ b/tests/cli/test_cli.py @@ -14,9 +14,8 @@ import matplotlib.style import numpy as np import pytest -from matplotlib.figure import Figure - from lobsterpy.cli import get_parser, run +from matplotlib.figure import Figure CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/cohp/test_analyze.py b/tests/cohp/test_analyze.py index d24afc1d..3180aaea 100644 --- a/tests/cohp/test_analyze.py +++ b/tests/cohp/test_analyze.py @@ -5,12 +5,11 @@ from pathlib import Path import pytest +from lobsterpy.cohp.analyze import Analysis from pymatgen.core import Structure from pymatgen.io.lobster import Bandoverlaps, Charge, Doscar, Lobsterin, Lobsterout from pymatgen.io.vasp import Vasprun -from lobsterpy.cohp.analyze import Analysis - CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/conftest.py b/tests/conftest.py index 60228fe4..ce950666 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,11 +5,10 @@ from pathlib import Path import pytest -from pymatgen.electronic_structure.cohp import CompleteCohp -from pymatgen.io.lobster import Charge, Doscar, Icohplist - from lobsterpy.cohp.analyze import Analysis from lobsterpy.cohp.describe import Description +from pymatgen.electronic_structure.cohp import CompleteCohp +from pymatgen.io.lobster import Charge, Doscar, Icohplist TestDir = Path(__file__).absolute().parent diff --git a/tests/featurize/test_batch.py b/tests/featurize/test_batch.py index c6573516..c5a30f2b 100644 --- a/tests/featurize/test_batch.py +++ b/tests/featurize/test_batch.py @@ -2,8 +2,6 @@ import pandas as pd import pytest -from pymatgen.analysis.graphs import StructureGraph - from lobsterpy.featurize.batch import ( BatchCoxxFingerprint, BatchDosFeaturizer, @@ -11,6 +9,7 @@ BatchSummaryFeaturizer, ) from lobsterpy.featurize.core import CoxxFingerprint, DosFingerprint +from pymatgen.analysis.graphs import StructureGraph CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/featurize/test_core.py b/tests/featurize/test_core.py index a0363aa4..44658df2 100644 --- a/tests/featurize/test_core.py +++ b/tests/featurize/test_core.py @@ -3,7 +3,6 @@ import numpy as np import pandas as pd import pytest - from lobsterpy.featurize.core import ( DosFingerprint, FeaturizeCharges, diff --git a/tests/featurize/test_utils.py b/tests/featurize/test_utils.py index 2240be94..ccd769d4 100644 --- a/tests/featurize/test_utils.py +++ b/tests/featurize/test_utils.py @@ -2,9 +2,8 @@ import shutil from pathlib import Path -from pymatgen.core import Structure - from lobsterpy.featurize import get_file_paths, get_structure_path +from pymatgen.core import Structure CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/plotting/test_plotting.py b/tests/plotting/test_plotting.py index 9c2cddc2..4f54251d 100644 --- a/tests/plotting/test_plotting.py +++ b/tests/plotting/test_plotting.py @@ -4,10 +4,6 @@ import numpy as np import pytest -from plotly.io import read_json -from pymatgen.electronic_structure.cohp import Cohp -from pymatgen.electronic_structure.core import Spin - from lobsterpy.cohp.describe import Description from lobsterpy.plotting import ( IcohpDistancePlotter, @@ -15,6 +11,9 @@ PlainCohpPlotter, PlainDosPlotter, ) +from plotly.io import read_json +from pymatgen.electronic_structure.cohp import Cohp +from pymatgen.electronic_structure.core import Spin CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/structuregraph/test_graph.py b/tests/structuregraph/test_graph.py index 69f96ce2..616080cb 100644 --- a/tests/structuregraph/test_graph.py +++ b/tests/structuregraph/test_graph.py @@ -1,7 +1,6 @@ from pathlib import Path import pytest - from lobsterpy.structuregraph.graph import LobsterGraph CurrentDir = Path(__file__).absolute().parent From c63ad3c376387142cfab9dc4aedded5437ef8409 Mon Sep 17 00:00:00 2001 From: naik-aakash Date: Sat, 20 Jul 2024 18:59:07 +0530 Subject: [PATCH 3/3] fix linting errros --- tests/cli/test_cli.py | 3 +-- tests/cohp/test_analyze.py | 3 +-- tests/conftest.py | 5 ++--- tests/featurize/test_batch.py | 3 +-- tests/featurize/test_core.py | 1 - tests/featurize/test_utils.py | 3 +-- tests/plotting/test_plotting.py | 7 +++---- tests/structuregraph/test_graph.py | 1 - 8 files changed, 9 insertions(+), 17 deletions(-) diff --git a/tests/cli/test_cli.py b/tests/cli/test_cli.py index e41bb021..420478fe 100644 --- a/tests/cli/test_cli.py +++ b/tests/cli/test_cli.py @@ -14,9 +14,8 @@ import matplotlib.style import numpy as np import pytest -from matplotlib.figure import Figure - from lobsterpy.cli import get_parser, run +from matplotlib.figure import Figure CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/cohp/test_analyze.py b/tests/cohp/test_analyze.py index d24afc1d..3180aaea 100644 --- a/tests/cohp/test_analyze.py +++ b/tests/cohp/test_analyze.py @@ -5,12 +5,11 @@ from pathlib import Path import pytest +from lobsterpy.cohp.analyze import Analysis from pymatgen.core import Structure from pymatgen.io.lobster import Bandoverlaps, Charge, Doscar, Lobsterin, Lobsterout from pymatgen.io.vasp import Vasprun -from lobsterpy.cohp.analyze import Analysis - CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/conftest.py b/tests/conftest.py index 60228fe4..ce950666 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,11 +5,10 @@ from pathlib import Path import pytest -from pymatgen.electronic_structure.cohp import CompleteCohp -from pymatgen.io.lobster import Charge, Doscar, Icohplist - from lobsterpy.cohp.analyze import Analysis from lobsterpy.cohp.describe import Description +from pymatgen.electronic_structure.cohp import CompleteCohp +from pymatgen.io.lobster import Charge, Doscar, Icohplist TestDir = Path(__file__).absolute().parent diff --git a/tests/featurize/test_batch.py b/tests/featurize/test_batch.py index c6573516..c5a30f2b 100644 --- a/tests/featurize/test_batch.py +++ b/tests/featurize/test_batch.py @@ -2,8 +2,6 @@ import pandas as pd import pytest -from pymatgen.analysis.graphs import StructureGraph - from lobsterpy.featurize.batch import ( BatchCoxxFingerprint, BatchDosFeaturizer, @@ -11,6 +9,7 @@ BatchSummaryFeaturizer, ) from lobsterpy.featurize.core import CoxxFingerprint, DosFingerprint +from pymatgen.analysis.graphs import StructureGraph CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/featurize/test_core.py b/tests/featurize/test_core.py index a0363aa4..44658df2 100644 --- a/tests/featurize/test_core.py +++ b/tests/featurize/test_core.py @@ -3,7 +3,6 @@ import numpy as np import pandas as pd import pytest - from lobsterpy.featurize.core import ( DosFingerprint, FeaturizeCharges, diff --git a/tests/featurize/test_utils.py b/tests/featurize/test_utils.py index 2240be94..ccd769d4 100644 --- a/tests/featurize/test_utils.py +++ b/tests/featurize/test_utils.py @@ -2,9 +2,8 @@ import shutil from pathlib import Path -from pymatgen.core import Structure - from lobsterpy.featurize import get_file_paths, get_structure_path +from pymatgen.core import Structure CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/plotting/test_plotting.py b/tests/plotting/test_plotting.py index 9c2cddc2..4f54251d 100644 --- a/tests/plotting/test_plotting.py +++ b/tests/plotting/test_plotting.py @@ -4,10 +4,6 @@ import numpy as np import pytest -from plotly.io import read_json -from pymatgen.electronic_structure.cohp import Cohp -from pymatgen.electronic_structure.core import Spin - from lobsterpy.cohp.describe import Description from lobsterpy.plotting import ( IcohpDistancePlotter, @@ -15,6 +11,9 @@ PlainCohpPlotter, PlainDosPlotter, ) +from plotly.io import read_json +from pymatgen.electronic_structure.cohp import Cohp +from pymatgen.electronic_structure.core import Spin CurrentDir = Path(__file__).absolute().parent TestDir = CurrentDir / "../" diff --git a/tests/structuregraph/test_graph.py b/tests/structuregraph/test_graph.py index 69f96ce2..616080cb 100644 --- a/tests/structuregraph/test_graph.py +++ b/tests/structuregraph/test_graph.py @@ -1,7 +1,6 @@ from pathlib import Path import pytest - from lobsterpy.structuregraph.graph import LobsterGraph CurrentDir = Path(__file__).absolute().parent