# 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 COOPs/COHPs or COBIs automatically."""
from __future__ import annotations
import warnings
from collections import Counter
from itertools import combinations_with_replacement, permutations
from pathlib import Path
from typing import Literal
import numpy as np
from monty.json import MSONable
from pymatgen.analysis.lobster_env import LobsterNeighbors
from pymatgen.core.structure import Structure
from pymatgen.electronic_structure.cohp import CompleteCohp
from pymatgen.electronic_structure.core import Spin
from pymatgen.io.lobster import (
Charge,
Icohplist,
MadelungEnergies,
)
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from scipy.integrate import trapezoid
from lobsterpy.utils import get_file_paths
POSCAR_WARNING = (
"Falling back to POSCAR, translations between individual atoms may differ from LOBSTER outputs. "
"Please note that translations in the LOBSTER outputs are consistent with CONTCAR "
"(also with POSCAR.lobster.vasp or POSCAR.vasp : written by LOBSTER >=v5)."
)
[docs]
class Analysis(MSONable):
"""
Class to automatically analyze COHP/COOP/COBI populations from Lobster.
Can be initialized using either file paths or pymatgen objects.
Pymatgen objects will be preferred in case both are supplied.
: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 charge_obj: pymatgen lobster.io.charge object
:param completecoxx_obj: pymatgen.electronic_structure.cohp.CompleteCohp object
:param icoxxlist_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 charge_obj is provided, Valences will be used (see pymatgen BVAnalyzer).
Otherwise, Mulliken charges from CHARGE.lobster are used by default.
: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,
structure: Structure | None = None,
icoxxlist_obj: Icohplist | None = None,
completecoxx_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: Literal["Mulliken", "Loewdin", "Valences"] = "Mulliken",
which_bonds: Literal["cation-anion", "all"] = "cation-anion",
):
"""
Initialize automatic bonding analysis.
Can be initialized using pymatgen objects.
: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 charge_obj: pymatgen lobster.io.charge object (Optional)
:param completecoxx_obj: pymatgen.electronic_structure.cohp.CompleteCohp object
:param icoxxlist_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 charge_obj is provided, Valences will be used
(see pymatgen BVAnalyzer). Otherwise, Mulliken charges from CHARGE.lobster are used by default.
: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.structure = structure
self.which_bonds = which_bonds
self.cutoff_icohp = cutoff_icohp
self.orbital_cutoff = orbital_cutoff
self.orbital_resolved = orbital_resolved
self.are_cobis = are_cobis
self.are_coops = are_coops
self.noise_cutoff = noise_cutoff
self.summed_spins = summed_spins
self.completecoxx_obj = completecoxx_obj
self.icoxxlist_obj = icoxxlist_obj
self.charge_obj = charge_obj
self.madelung_obj = madelung_obj
self.type_charge = type_charge.capitalize()
if self.type_charge not in {"Mulliken", "Loewdin", "Valences"}:
raise ValueError(f"type_charge must be 'Mulliken', 'Loewdin', or 'Valences'. Got '{type_charge}'.")
if self.type_charge == "Loewdin":
warnings.warn("Support for Loewdin charges is currently experimental. Use with caution!")
if self.type_charge == "Valences":
warnings.warn(
"Using Valences for chemical environment analysis. It is recommended to use "
" 'Mulliken' or 'Loewdin' charges."
)
if self.which_bonds not in ["cation-anion", "all"]:
raise ValueError("only cation-anion or all bonds analysis is supported so far.")
self.setup_env()
self.get_information_all_bonds(summed_spins=self.summed_spins)
self.orbital_resolved = orbital_resolved
self.set_condensed_bonding_analysis()
self.set_summary_dicts()
[docs]
def setup_env(self):
"""
Set up the light structure environments based on COHPs using this method.
Returns:
None
"""
sga = SpacegroupAnalyzer(self.structure)
symmetry_dataset = sga.get_symmetry_dataset()
self.list_equivalent_sites = symmetry_dataset.equivalent_atoms
self.seq_equivalent_sites = list(set(self.list_equivalent_sites))
self.spg = symmetry_dataset.international
lob_neigh_kwargs = {
"perc_strength_icohp": self.cutoff_icohp,
"valences_from_charges": self.type_charge != "Valences",
"adapt_extremum_to_add_cond": True,
"are_cobis": self.are_cobis,
"are_coops": self.are_coops,
"noise_cutoff": self.noise_cutoff,
"which_charge": self.type_charge,
"structure": self.structure,
"icoxxlist_obj": self.icoxxlist_obj,
"charge_obj": self.charge_obj,
"additional_condition": 1 if self.which_bonds == "cation-anion" else 0,
}
try:
self.chemenv = LobsterNeighbors(**lob_neigh_kwargs)
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"
):
if self.which_bonds == "cation-anion":
raise ValueError(
"No cations detected. Consider analyzing all bonds instead of only cation-anion bonds."
)
else:
raise err
# determine cations and anions
# try:
self.lse = self.chemenv.get_light_structure_environment(
only_cation_environments=(self.which_bonds == "cation-anion"),
on_error="warn",
)
[docs]
@classmethod
def from_files(
cls,
structure_path: str | Path,
icoxxlist_path: str | Path,
coxxcar_path: str | Path,
charge_path: str | Path | None = None,
madelung_path: str | Path | None = None,
**kwargs,
) -> Analysis:
"""
Create Analysis from explicit file paths.
:param structure_path: path to structure (e.g., `CONTCAR` (preferred), `POSCAR.lobster` or `POSCAR`)
:param icoxxlist_path: path to `ICOHPLIST.lobster` or `ICOBILIST.lobster` or `ICOOPLIST.lobster`.
:param coxxcar_path: path to `COHPCAR.lobster` or `COBICAR.lobster` or `COOPCAR.lobster` .
:param charge_path: path to `CHARGE.lobster`.
:param madelung_path: path to `MadelungEnergies.lobster`.
:param kwargs: Additional arguments passed to Analysis constructor
Returns:
Analysis object
"""
are_cobis = kwargs.get("are_cobis")
are_coops = kwargs.get("are_coops")
# infer file type
if are_cobis is None or are_coops is None:
fname = str(icoxxlist_path).lower()
if "icobilist" in fname:
kwargs["are_cobis"] = True if are_cobis is None else are_cobis
kwargs["are_coops"] = False if are_coops is None else are_coops
elif "icooplist" in fname:
kwargs["are_coops"] = True if are_coops is None else are_coops
kwargs["are_cobis"] = False if are_cobis is None else are_cobis
else:
kwargs["are_cobis"] = False if are_cobis is None else are_cobis
kwargs["are_coops"] = False if are_coops is None else are_coops
if Path(structure_path).name in ("POSCAR", "POSCAR.gz"):
warnings.warn(POSCAR_WARNING, UserWarning)
structure = Structure.from_file(structure_path)
icoxxlist_obj = Icohplist(
filename=icoxxlist_path, are_cobis=kwargs.get("are_cobis"), are_coops=kwargs.get("are_coops")
)
completecoxx_obj = CompleteCohp.from_file(
fmt="LOBSTER",
filename=coxxcar_path,
structure_file=structure_path,
are_cobis=kwargs.get("are_cobis"),
are_coops=kwargs.get("are_coops"),
)
charge_obj = Charge(filename=charge_path) if charge_path else None
madelung_obj = MadelungEnergies(filename=madelung_path) if madelung_path else None
return cls(
structure=structure,
icoxxlist_obj=icoxxlist_obj,
completecoxx_obj=completecoxx_obj,
charge_obj=charge_obj,
madelung_obj=madelung_obj,
**kwargs,
)
[docs]
@classmethod
def from_directory(
cls,
path_to_lobster_calc: str | Path,
read_madelung_energies: bool = False,
analyze_coops: bool = False,
analyze_cobis: bool = False,
type_charge: Literal["Mulliken", "Loewdin", "Valences"] = "Mulliken",
**kwargs,
) -> Analysis:
"""
Create Analysis from a directory containing LOBSTER calculation files.
:param path_to_lobster_calc: Path to directory with LOBSTER calculation files
:param read_madelung_energies: Whether to read Madelung energies from file
:param analyze_coops: Whether to analyze COOPs instead of COHPs
:param analyze_cobis: Whether to analyze COBIs instead of COHPs
:param type_charge: If no charge file is provided, Valences will be used (see pymatgen BVAnalyzer).
Otherwise, Mulliken charges from CHARGE.lobster are used by default.
:param kwargs: Additional arguments passed to Analysis constructor
Returns:
Analysis object
"""
if analyze_coops:
requested_files = ["structure", "icooplist", "coopcar"]
elif analyze_cobis:
requested_files = ["structure", "icobilist", "cobicar"]
else:
requested_files = ["structure", "icohplist", "cohpcar"]
if type_charge != "Valences":
requested_files.append("charge")
if read_madelung_energies:
requested_files.append("madelung")
file_paths = get_file_paths(
path_to_lobster_calc=path_to_lobster_calc,
requested_files=requested_files,
)
structure_path = file_paths.get("structure")
icoxxlist_path = file_paths.get("icohplist")
cooxxar_path = file_paths.get("cohpcar")
charge_path = file_paths.get("charge")
madelung_path = file_paths.get("madelung")
kwargs["type_charge"] = type_charge
if analyze_coops:
icoxxlist_path = file_paths.get("icooplist")
cooxxar_path = file_paths.get("coopcar")
elif analyze_cobis:
icoxxlist_path = file_paths.get("icobilist")
cooxxar_path = file_paths.get("cobicar")
return cls.from_files(
structure_path=structure_path,
icoxxlist_path=icoxxlist_path,
coxxcar_path=cooxxar_path,
charge_path=charge_path if type_charge != "Valences" else None,
madelung_path=madelung_path if read_madelung_energies else None,
**kwargs,
)
[docs]
def as_dict(self) -> dict:
"""
Serialize the Analysis object to a dictionary.
Returns:
Dictionary representation of the Analysis object
"""
return {
"structure": self.structure.as_dict() if self.structure else None,
"completecoxx_obj": self.completecoxx_obj.as_dict() if self.completecoxx_obj else None,
"icoxxlist_obj": self.icoxxlist_obj.as_dict() if self.icoxxlist_obj else None,
"charge_obj": self.charge_obj.as_dict() if self.charge_obj else None,
"madelung_obj": self.madelung_obj.as_dict() if self.madelung_obj else None,
"are_cobis": self.are_cobis,
"are_coops": self.are_coops,
"cutoff_icohp": self.cutoff_icohp,
"noise_cutoff": self.noise_cutoff,
"orbital_cutoff": self.orbital_cutoff,
"orbital_resolved": self.orbital_resolved,
"summed_spins": self.summed_spins,
"type_charge": self.type_charge,
"which_bonds": self.which_bonds,
"start": self.start,
"@module": self.__class__.__module__,
"@class": self.__class__.__name__,
}
[docs]
@classmethod
def from_dict(cls, d: dict) -> Analysis:
"""
Reconstruct an Analysis object from a dictionary.
This method properly handles deserialization of nested pymatgen objects
that have @module and @class keys.
Args:
d: Dictionary representation of the Analysis object
Returns:
Analysis object
"""
d = d.copy()
# Remove monty metadata keys
d.pop("@module", None)
d.pop("@class", None)
d.pop("@version", None)
# Manually deserialize nested MSONable objects that have @module and @class
if d.get("structure") and isinstance(d["structure"], dict) and "@class" in d["structure"]:
d["structure"] = Structure.from_dict(d["structure"])
if d.get("completecoxx_obj") and isinstance(d["completecoxx_obj"], dict) and "@class" in d["completecoxx_obj"]:
d["completecoxx_obj"] = CompleteCohp.from_dict(d["completecoxx_obj"])
if d.get("icoxxlist_obj") and isinstance(d["icoxxlist_obj"], dict) and "@class" in d["icoxxlist_obj"]:
d["icoxxlist_obj"] = Icohplist.from_dict(d["icoxxlist_obj"])
if d.get("charge_obj") and isinstance(d["charge_obj"], dict) and "@class" in d["charge_obj"]:
d["charge_obj"] = Charge.from_dict(d["charge_obj"])
if d.get("madelung_obj") and isinstance(d["madelung_obj"], dict) and "@class" in d["madelung_obj"]:
d["madelung_obj"] = MadelungEnergies.from_dict(d["madelung_obj"])
return cls(**d)
[docs]
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 index, 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[index][i] for i in indices]
label_data.update({search_item: filtered_bond_label_list})
return label_data
@property
def charges(self) -> list[float]:
"""Charges used for chemical environment analysis.
List of charges for each site in the structure.
"""
return self.chemenv.valences
@property
def completecoxx(self) -> CompleteCohp:
"""
Pymatgen CompleteCohp object.
Depending on the type of files read during Analyse class initialization,
it can contain COHP, COOP or COBI data.
"""
return self.chemenv.completecohp
@property
def icoxxlist(self) -> Icohplist:
"""
Pymatgen Icohplist object.
Depending on the type of files read during Analyse class initialization,
it can contain ICOHP, ICOOP or ICOBI data.
"""
return self.chemenv.ICOHP
@property
def lobsterneighbors(self) -> LobsterNeighbors:
"""Pymatgen LobsterNeighbors object."""
return self.chemenv
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_contri_perc > 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_contri_perc > 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
[docs]
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,
)
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])
else:
pos.append(0)
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])
else:
neg.append(0)
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]))]
[docs]
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.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:
if self.type_charge == "Mulliken":
madelung_energy = self.madelung_obj.madelungenergies_mulliken
elif self.type_charge == "Loewdin":
madelung_energy = self.madelung_obj.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,
}
[docs]
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))
