# 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
[docs]
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
[docs]
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)
[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 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_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,
)
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]))]
[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.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,
}
[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))
[docs]
@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,
fp_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,
fp_type=orb.name,
)
tani_orb = round(
dos_vasp.get_dos_fp_similarity(fp_lobster_orb, fp_vasp_orb, metric="tanimoto"),
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,
fp_type="summed_pdos",
)
fp_vasp = dos_vasp.get_dos_fp(
min_e=min_e,
max_e=max_e,
n_bins=n_bins,
normalize=True,
fp_type="summed_pdos",
)
tanimoto_summed = round(dos_vasp.get_dos_fp_similarity(fp_lobster, fp_vasp, metric="tanimoto"), 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
