diff --git a/src/pymatgen/electronic_structure/bandstructure.py b/src/pymatgen/electronic_structure/bandstructure.py index f844674be2b..73378506807 100644 --- a/src/pymatgen/electronic_structure/bandstructure.py +++ b/src/pymatgen/electronic_structure/bandstructure.py @@ -1,4 +1,4 @@ -"""This module provides classes to define everything related to band structures.""" +"""This module provides classes to define things related to band structures.""" from __future__ import annotations @@ -7,7 +7,7 @@ import re import warnings from collections import defaultdict -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, overload import numpy as np from monty.json import MSONable @@ -20,6 +20,7 @@ if TYPE_CHECKING: from typing import Any + from numpy.typing import NDArray from typing_extensions import Self __author__ = "Geoffroy Hautier, Shyue Ping Ong, Michael Kocher" @@ -32,14 +33,13 @@ class Kpoint(MSONable): - """Store kpoint objects. A kpoint is defined with a lattice and frac - or Cartesian coordinates syntax similar than the site object in - pymatgen.core.structure. + """A kpoint defined with a lattice and frac or Cartesian coordinates, + similar to the Site object in pymatgen.core.structure. """ def __init__( self, - coords: np.ndarray, + coords: NDArray, lattice: Lattice, to_unit_cell: bool = False, coords_are_cartesian: bool = False, @@ -47,15 +47,14 @@ def __init__( ) -> None: """ Args: - coords: coordinate of the kpoint as a numpy array - lattice: A pymatgen.core.Lattice object representing - the reciprocal lattice of the kpoint - to_unit_cell: Translates fractional coordinate to the basic unit + coords (NDArray): Coordinate of the Kpoint. + lattice (Lattice): The reciprocal lattice of the kpoint. + to_unit_cell (bool): Translate fractional coordinate to the basic unit cell, i.e., all fractional coordinates satisfy 0 <= a < 1. Defaults to False. - coords_are_cartesian: Boolean indicating if the coordinates given are - in Cartesian or fractional coordinates (by default fractional) - label: the label of the kpoint if any (None by default). + coords_are_cartesian (bool): Whether the coordinates given are + in Cartesian (True) or fractional coordinates (by default fractional). + label (str): The label of the Kpoint if any (None by default). """ self._lattice = lattice self._frac_coords = lattice.get_fractional_coords(coords) if coords_are_cartesian else coords @@ -67,11 +66,24 @@ def __init__( self._cart_coords = lattice.get_cartesian_coords(self._frac_coords) + def __str__(self) -> str: + """String with fractional, Cartesian coordinates and label.""" + return f"{self.frac_coords} {self.cart_coords} {self.label}" + + def __eq__(self, other: object) -> bool: + """Whether two Kpoints are equal.""" + if not isinstance(other, type(self)): + return NotImplemented + + return ( + np.allclose(self.frac_coords, other.frac_coords) + and self.lattice == other.lattice + and self.label == other.label + ) + @property def lattice(self) -> Lattice: - """The lattice associated with the kpoint. It's a - pymatgen.core.Lattice object. - """ + """The lattice associated with the kpoint, as a Lattice object.""" return self._lattice @property @@ -85,13 +97,13 @@ def label(self, label: str | None) -> None: self._label = label @property - def frac_coords(self) -> np.ndarray: - """The fractional coordinates of the kpoint as a numpy array.""" + def frac_coords(self) -> NDArray: + """The fractional coordinates of the kpoint as a NumPy array.""" return np.copy(self._frac_coords) @property - def cart_coords(self) -> np.ndarray: - """The Cartesian coordinates of the kpoint as a numpy array.""" + def cart_coords(self) -> NDArray: + """The Cartesian coordinates of the kpoint as a NumPy array.""" return np.copy(self._cart_coords) @property @@ -109,22 +121,8 @@ def c(self) -> float: """Fractional c coordinate of the kpoint.""" return self._frac_coords[2] - def __str__(self) -> str: - """Get a string with fractional, Cartesian coordinates and label.""" - return f"{self.frac_coords} {self.cart_coords} {self.label}" - - def __eq__(self, other: object) -> bool: - """Check if two kpoints are equal.""" - if not isinstance(other, Kpoint): - return NotImplemented - return ( - np.allclose(self.frac_coords, other.frac_coords) - and self.lattice == other.lattice - and self.label == other.label - ) - def as_dict(self) -> dict[str, Any]: - """JSON-serializable dict representation of a kpoint.""" + """JSON-serializable dict representation of the kpoint.""" return { "lattice": self.lattice.as_dict(), "fcoords": self.frac_coords.tolist(), @@ -136,7 +134,7 @@ def as_dict(self) -> dict[str, Any]: @classmethod def from_dict(cls, dct: dict) -> Self: - """Create from dict. + """Create from a dict. Args: dct (dict): A dict with all data for a kpoint object. @@ -149,59 +147,57 @@ def from_dict(cls, dct: dict) -> Self: class BandStructure: - """This is the most generic band structure data possible - it's defined by a list of kpoints + energies for each of them. + """Generic band structure data, defined by a list of Kpoints + and corresponding energies for each of them. Attributes: - kpoints (list): The list of kpoints (as Kpoint objects) in the band structure. + kpoints (list[Kpoint]): Kpoints in the band structure. lattice_rec (Lattice): The reciprocal lattice of the band structure. - efermi (float): The Fermi energy. - is_spin_polarized (bool): True if the band structure is spin-polarized. - bands (dict): The energy eigenvalues as a {spin: array}. Note that the use of an - array is necessary for computational as well as memory efficiency due to the large - amount of numerical data. The indices of the array are [band_index, kpoint_index]. - nb_bands (int): Returns the number of bands in the band structure. - structure (Structure): Returns the structure. - projections (dict): The projections as a {spin: array}. Note that the use of an - array is necessary for computational as well as memory efficiency due to the large - amount of numerical data. The indices of the array are [band_index, kpoint_index, - orbital_index, ion_index]. + efermi (float): The Fermi level. + is_spin_polarized (bool): Whether the band structure is spin-polarized. + bands (dict[Spin, NDArray]): The energy eigenvalues. Note that the use of an + array is necessary for computational and memory efficiency due to the large + amount of numerical data. The indices of the array are (band_index, kpoint_index). + nb_bands (int): The number of bands in the band structure. + structure (Structure): The structure. + projections (dict[Spin, NDArray]): The projections. Note that the use of an + array is necessary for computational and memory efficiency due to the large + amount of numerical data. The indices of the array are (band_index, kpoint_index, + orbital_index, ion_index). """ def __init__( self, - kpoints: np.ndarray, - eigenvals: dict[Spin, np.ndarray], + kpoints: NDArray, + eigenvals: dict[Spin, NDArray], lattice: Lattice, efermi: float, - labels_dict=None, + labels_dict: dict[str, Kpoint] | None = None, coords_are_cartesian: bool = False, structure: Structure | None = None, - projections: dict[Spin, np.ndarray] | None = None, + projections: dict[Spin, NDArray] | None = None, ) -> None: """ Args: - kpoints: list of kpoint as numpy arrays, in frac_coords of the - given lattice by default - eigenvals: dict of energies for spin up and spin down - {Spin.up:[][],Spin.down:[][]}, the first index of the array - [][] refers to the band and the second to the index of the - kpoint. The kpoints are ordered according to the order of the - kpoints array. If the band structure is not spin polarized, we - only store one data set under Spin.up - lattice: The reciprocal lattice as a pymatgen Lattice object. - Pymatgen uses the physics convention of reciprocal lattice vectors - WITH a 2*pi coefficient - efermi (float): Fermi energy - labels_dict: (dict) of {} this links a kpoint (in frac coords or - Cartesian coordinates depending on the coords) to a label. - coords_are_cartesian: Whether coordinates are cartesian. - structure: The crystal structure (as a pymatgen Structure object) + kpoints (NDArray): Kpoint as NumPy array, in frac_coords of the + given lattice by default. + eigenvals (dict): Energies for spin up and spin down as + {Spin.up:[][], Spin.down:[][]}, the first index of the array + [][] refers to the band and the second to the index of the kpoint. + The kpoints are ordered according to the kpoints array. + If the band structure is not spin polarized, we + only store one data set under Spin.up. + lattice (Lattice): The reciprocal lattice. Pymatgen uses the physics + convention of reciprocal lattice vectors with a 2*pi coefficient. + efermi (float): The Fermi level. + labels_dict (dict[str, Kpoint]): Dict mapping label to Kpoint. + coords_are_cartesian (bool): Whether coordinates are cartesian. + structure (Structure): The crystal structure associated with the band structure. This is needed if we - provide projections to the band structure - projections: dict of orbital projections as {spin: array}. The - indices of the array are [band_index, kpoint_index, orbital_index, - ion_index].If the band structure is not spin polarized, we only + provide projections to the band structure. + projections (dict[Spin, NDArray]): Orbital projections. The + indices of the array are (band_index, kpoint_index, orbital_index, + ion_index). If the band structure is not spin polarized, we only store one data set under Spin.up. """ self.efermi = efermi @@ -215,58 +211,58 @@ def __init__( if labels_dict is None: labels_dict = {} - if len(self.projections) != 0 and self.structure is None: + if self.projections and self.structure is None: raise RuntimeError("if projections are provided a structure object is also required") - for k in kpoints: - # let see if this kpoint has been assigned a label + for kpt in kpoints: + # Check if this Kpoint has a label label = None for c in labels_dict: - if np.linalg.norm(k - np.array(labels_dict[c])) < 0.0001: + if np.linalg.norm(kpt - np.array(labels_dict[c])) < 0.0001: label = c self.labels_dict[label] = Kpoint( - k, + kpt, lattice, label=label, coords_are_cartesian=coords_are_cartesian, ) - self.kpoints.append(Kpoint(k, lattice, label=label, coords_are_cartesian=coords_are_cartesian)) + self.kpoints.append(Kpoint(kpt, lattice, label=label, coords_are_cartesian=coords_are_cartesian)) self.bands = {spin: np.array(v) for spin, v in eigenvals.items()} self.nb_bands = len(eigenvals[Spin.up]) self.is_spin_polarized = len(self.bands) == 2 - def get_projection_on_elements(self): - """Get a dictionary of projections on elements. + def get_projection_on_elements(self) -> dict[Spin, NDArray]: + """Get projections on elements. Returns: - a dictionary in the {Spin.up:[][{Element: [values]}], - Spin.down:[][{Element: [values]}]} format - if there is no projections in the band structure - returns an empty dict + dict[Spin, NDArray]: Dict in {Spin.up:[][{Element: [values]}], + Spin.down: [][{Element: [values]}]} format. + If there is no projections in the band structure, return {}. """ - result = {} - for spin, v in self.projections.items(): - result[spin] = [[defaultdict(float) for i in range(len(self.kpoints))] for j in range(self.nb_bands)] + assert self.structure is not None + result: dict[Spin, NDArray] = {} + for spin, val in self.projections.items(): + result[spin] = [[defaultdict(float) for _ in range(len(self.kpoints))] for _ in range(self.nb_bands)] for i, j, k in itertools.product( range(self.nb_bands), range(len(self.kpoints)), range(len(self.structure)), ): - result[spin][i][j][str(self.structure[k].specie)] += np.sum(v[i, j, :, k]) + result[spin][i][j][str(self.structure[k].specie)] += np.sum(val[i, j, :, k]) return result def get_projections_on_elements_and_orbitals(self, el_orb_spec: dict[str, list[str]]): - """Get a dictionary of projections on elements and specific orbitals. + """Get projections on elements and specific orbitals. Args: - el_orb_spec (dict[str, list[str]]): A dictionary of elements and orbitals which - to project onto. Format is {Element: [orbitals]}, e.g. {'Cu':['d','s']}. + el_orb_spec (dict[str, list[str]]): Elements and orbitals to project onto. + Format is {Element: [orbitals]}, e.g. {"Cu": ["d", "s"]}. Returns: - A dictionary of projections on elements in the - {Spin.up:[][{Element:{orb:values}}], - Spin.down:[][{Element:{orb:values}}]} format - if there is no projections in the band structure returns an empty dict. + dict[str, list[str]: Projections on elements in the + {Spin.up: [][{Element: {orb: values}}], + Spin.down: [][{Element: {orb: values}}]} format. + If there is no projections in the band structure, return {}. """ if self.structure is None: raise ValueError("Structure is required for this method") @@ -274,8 +270,8 @@ def get_projections_on_elements_and_orbitals(self, el_orb_spec: dict[str, list[s species_orb_spec = {get_el_sp(el): orbs for el, orbs in el_orb_spec.items()} for spin, v in self.projections.items(): result[spin] = [ - [{str(e): defaultdict(float) for e in species_orb_spec} for i in range(len(self.kpoints))] - for j in range(self.nb_bands) + [{str(e): defaultdict(float) for e in species_orb_spec} for _ in range(len(self.kpoints))] + for _ in range(self.nb_bands) ] for i, j, k in itertools.product( @@ -290,12 +286,12 @@ def get_projections_on_elements_and_orbitals(self, el_orb_spec: dict[str, list[s result[spin][i][j][str(sp)][o] += v[i][j][orb_i][k] return result - def is_metal(self, efermi_tol=1e-4) -> bool: - """Check if the band structure indicates a metal by looking if the fermi - level crosses a band. + def is_metal(self, efermi_tol: float = 1e-4) -> bool: + """Check if the band structure indicates a metal, + by looking at if the fermi level crosses a band. Returns: - bool: True if a metal. + bool: True if is metal. """ for vals in self.bands.values(): for idx in range(self.nb_bands): @@ -303,26 +299,26 @@ def is_metal(self, efermi_tol=1e-4) -> bool: return True return False - def get_vbm(self): - """Get data about the VBM. + def get_vbm(self) -> dict[str, Any]: + """Get data about the valence band maximum (VBM). Returns: - dict: With keys "band_index", "kpoint_index", "kpoint", "energy" - - "band_index": A dict with spin keys pointing to a list of the - indices of the band containing the VBM (please note that you - can have several bands sharing the VBM) {Spin.up:[], - Spin.down:[]} - - "kpoint_index": The list of indices in self.kpoints for the - kpoint VBM. Please note that there can be several - kpoint_indices relating to the same kpoint (e.g., Gamma can - occur at different spots in the band structure line plot) - - "kpoint": The kpoint (as a kpoint object) - - "energy": The energy of the VBM - - "projections": The projections along sites and orbitals of the - VBM if any projection data is available (else it is an empty - dictionary). The format is similar to the projections field in - BandStructure: {spin:{'Orbital': [proj]}} where the array - [proj] is ordered according to the sites in structure + dict of keys "band_index", "kpoint_index", "kpoint", "energy": + - "band_index" (dict): A dict with spin keys pointing to a list of the + indices of the band containing the VBM (please note that you + can have several bands sharing the VBM) {Spin.up:[], + Spin.down:[]}. + - "kpoint_index": The list of indices in self.kpoints for the + kpoint VBM. Please note that there can be several + kpoint_indices relating to the same kpoint (e.g., Gamma can + occur at different spots in the band structure line plot). + - "kpoint" (Kpoint): The kpoint. + - "energy" (float): The energy of the VBM. + - "projections": The projections along sites and orbitals of the + VBM if any projection data is available (else it is an empty + dictionary). The format is similar to the projections field in + BandStructure: {spin:{'Orbital': [proj]}} where the array + [proj] is ordered according to the sites in structure. """ if self.is_metal(): return { @@ -332,33 +328,36 @@ def get_vbm(self): "energy": None, "projections": {}, } + max_tmp = -float("inf") index = kpoint_vbm = None for value in self.bands.values(): - for i, j in zip(*np.where(value < self.efermi)): - if value[i, j] > max_tmp: - max_tmp = float(value[i, j]) + for idx, j in zip(*np.where(value < self.efermi)): + if value[idx, j] > max_tmp: + max_tmp = float(value[idx, j]) index = j kpoint_vbm = self.kpoints[j] list_ind_kpts = [] - if kpoint_vbm.label is not None: - for i, kpt in enumerate(self.kpoints): + if kpoint_vbm is not None and kpoint_vbm.label is not None: + for idx, kpt in enumerate(self.kpoints): if kpt.label == kpoint_vbm.label: - list_ind_kpts.append(i) + list_ind_kpts.append(idx) else: list_ind_kpts.append(index) - # get all other bands sharing the vbm + + # Get all other bands sharing the VBM list_ind_band = defaultdict(list) for spin in self.bands: - for i in range(self.nb_bands): - if math.fabs(self.bands[spin][i][index] - max_tmp) < 0.001: - list_ind_band[spin].append(i) + for idx in range(self.nb_bands): + if math.fabs(self.bands[spin][idx][index] - max_tmp) < 0.001: + list_ind_band[spin].append(idx) proj = {} for spin, value in self.projections.items(): if len(list_ind_band[spin]) == 0: continue proj[spin] = value[list_ind_band[spin][0]][list_ind_kpts[0]] + return { "band_index": list_ind_band, "kpoint_index": list_ind_kpts, @@ -367,25 +366,25 @@ def get_vbm(self): "projections": proj, } - def get_cbm(self): - """Get data about the CBM. + def get_cbm(self) -> dict[str, Any]: + """Get data about the conduction band minimum (CBM). Returns: - dict[str, Any]: with keys band_index, kpoint_index, kpoint, energy. - - "band_index": A dict with spin keys pointing to a list of the + dict of keys "band_index", "kpoint_index", "kpoint", "energy": + - "band_index" (dict): A dict with spin keys pointing to a list of the indices of the band containing the CBM (please note that you - can have several bands sharing the CBM) {Spin.up:[], Spin.down:[]} + can have several bands sharing the CBM) {Spin.up:[], Spin.down:[]}. - "kpoint_index": The list of indices in self.kpoints for the kpoint CBM. Please note that there can be several kpoint_indices relating to the same kpoint (e.g., Gamma can - occur at different spots in the band structure line plot) - - "kpoint": The kpoint (as a kpoint object) - - "energy": The energy of the CBM + occur at different spots in the band structure line plot). + - "kpoint" (Kpoint): The kpoint. + - "energy" (float): The energy of the CBM. - "projections": The projections along sites and orbitals of the CBM if any projection data is available (else it is an empty dictionary). The format is similar to the projections field in BandStructure: {spin:{'Orbital': [proj]}} where the array - [proj] is ordered according to the sites in structure + [proj] is ordered according to the sites in structure. """ if self.is_metal(): return { @@ -395,30 +394,30 @@ def get_cbm(self): "energy": None, "projections": {}, } - max_tmp = float("inf") + max_tmp = float("inf") index = kpoint_cbm = None for value in self.bands.values(): - for i, j in zip(*np.where(value >= self.efermi)): - if value[i, j] < max_tmp: - max_tmp = float(value[i, j]) + for idx, j in zip(*np.where(value >= self.efermi)): + if value[idx, j] < max_tmp: + max_tmp = float(value[idx, j]) index = j kpoint_cbm = self.kpoints[j] list_index_kpoints = [] - if kpoint_cbm.label is not None: - for i, kpt in enumerate(self.kpoints): + if kpoint_cbm is not None and kpoint_cbm.label is not None: + for idx, kpt in enumerate(self.kpoints): if kpt.label == kpoint_cbm.label: - list_index_kpoints.append(i) + list_index_kpoints.append(idx) else: list_index_kpoints.append(index) - # get all other bands sharing the cbm + # Get all other bands sharing the CBM list_index_band = defaultdict(list) for spin in self.bands: - for i in range(self.nb_bands): - if math.fabs(self.bands[spin][i][index] - max_tmp) < 0.001: - list_index_band[spin].append(i) + for idx in range(self.nb_bands): + if math.fabs(self.bands[spin][idx][index] - max_tmp) < 0.001: + list_index_band[spin].append(idx) proj = {} for spin, value in self.projections.items(): if len(list_index_band[spin]) == 0: @@ -433,22 +432,25 @@ def get_cbm(self): "projections": proj, } - def get_band_gap(self): - r"""Get band gap data. + def get_band_gap(self) -> dict[str, Any]: + r"""Get band gap. Returns: - A dict {"energy","direct","transition"}: - "energy": band gap energy - "direct": A boolean telling if the gap is direct or not - "transition": kpoint labels of the transition (e.g., "\\Gamma-X") + dict of keys "energy", "direct", "transition": + "energy" (float): Band gap energy. + "direct" (bool): Whether the gap is direct. + "transition" (str): Kpoint labels of the transition (e.g., "\\Gamma-X"). """ if self.is_metal(): return {"energy": 0.0, "direct": False, "transition": None} + cbm = self.get_cbm() vbm = self.get_vbm() - result = {"direct": False, "energy": 0.0, "transition": None} - - result["energy"] = cbm["energy"] - vbm["energy"] + result = { + "direct": False, + "transition": None, + "energy": cbm["energy"] - vbm["energy"], + } if (cbm["kpoint"].label is not None and cbm["kpoint"].label == vbm["kpoint"].label) or np.linalg.norm( cbm["kpoint"].cart_coords - vbm["kpoint"].cart_coords @@ -464,16 +466,16 @@ def get_band_gap(self): return result - def get_direct_band_gap_dict(self): - """Get a dictionary of information about the direct - band gap. + def get_direct_band_gap_dict(self) -> dict[Spin, dict[str, Any]]: + """Get information about the direct band gap. Returns: - a dictionary of the band gaps indexed by spin - along with their band indices and k-point index + dict[Spin, dict[str, Any]]: The band gaps indexed by spin + along with their band indices and kpoint index. """ if self.is_metal(): raise ValueError("get_direct_band_gap_dict should only be used with non-metals") + direct_gap_dict = {} for spin, v in self.bands.items(): above = v[np.all(v > self.efermi, axis=1)] @@ -493,35 +495,42 @@ def get_direct_band_gap_dict(self): } return direct_gap_dict - def get_direct_band_gap(self): + def get_direct_band_gap(self) -> float: """Get the direct band gap. Returns: - the value of the direct band gap + float: The direct band gap value. """ if self.is_metal(): return 0.0 + dg = self.get_direct_band_gap_dict() return min(v["value"] for v in dg.values()) - def get_sym_eq_kpoints(self, kpoint, cartesian=False, tol: float = 1e-2): - """Get a list of unique symmetrically equivalent k-points. + def get_sym_eq_kpoints( + self, + kpoint: NDArray, + cartesian: bool = False, + tol: float = 1e-2, + ) -> NDArray: + """Get unique symmetrically equivalent Kpoints. Args: - kpoint (1x3 array): coordinate of the k-point - cartesian (bool): kpoint is in Cartesian or fractional coordinates - tol (float): tolerance below which coordinates are considered equal + kpoint (1x3 array): Coordinate of the Kpoint. + cartesian (bool): Whether kpoint is in Cartesian or fractional coordinates. + tol (float): Tolerance below which coordinates are considered equal. Returns: - list[1x3 array] | None: if structure is not available returns None + (1x3 NDArray) | None: None if structure is not available. """ if not self.structure: return None + sg = SpacegroupAnalyzer(self.structure) symm_ops = sg.get_point_group_operations(cartesian=cartesian) points = np.dot(kpoint, [m.rotation_matrix for m in symm_ops]) rm_list = [] - # identify and remove duplicates from the list of equivalent k-points: + # Identify and remove duplicates from equivalent k-points for i in range(len(points) - 1): for j in range(i + 1, len(points)): if np.allclose(pbc_diff(points[i], points[j]), [0, 0, 0], tol): @@ -529,25 +538,28 @@ def get_sym_eq_kpoints(self, kpoint, cartesian=False, tol: float = 1e-2): break return np.delete(points, rm_list, axis=0) - def get_kpoint_degeneracy(self, kpoint, cartesian=False, tol: float = 1e-2): - """Get degeneracy of a given k-point based on structure symmetry. + def get_kpoint_degeneracy( + self, + kpoint: NDArray, + cartesian: bool = False, + tol: float = 1e-2, + ) -> NDArray | None: + """Get degeneracy of a given kpoint based on structure symmetry. Args: - kpoint (1x3 array): coordinate of the k-point - cartesian (bool): kpoint is in Cartesian or fractional coordinates - tol (float): tolerance below which coordinates are considered equal. + kpoint (1x3 NDArray): Coordinate of the k-point. + cartesian (bool): Whether kpoint is in Cartesian or fractional coordinates. + tol (float): Tolerance below which coordinates are considered equal. Returns: - int | None: degeneracy or None if structure is not available + int | None: Degeneracy, or None if structure is not available. """ all_kpts = self.get_sym_eq_kpoints(kpoint, cartesian, tol=tol) - if all_kpts is not None: - return len(all_kpts) - return None + return len(all_kpts) if all_kpts is not None else None - def as_dict(self): + def as_dict(self) -> dict[str, Any]: """JSON-serializable dict representation of BandStructure.""" - dct = { + dct: dict[str, Any] = { "@module": type(self).__module__, "@class": type(self).__name__, "lattice_rec": self.lattice_rec.as_dict(), @@ -555,7 +567,7 @@ def as_dict(self): "kpoints": [], } # kpoints are not kpoint objects dicts but are frac coords (this makes - # the dict smaller and avoids the repetition of the lattice + # the dict smaller and avoids the repetition of the lattice). for k in self.kpoints: dct["kpoints"].append(k.as_dict()["fcoords"]) @@ -579,39 +591,40 @@ def as_dict(self): dct["labels_dict"] = {} dct["is_spin_polarized"] = self.is_spin_polarized - # MongoDB does not accept keys starting with $. Add a blank space to fix the problem + # MongoDB does not accept keys starting with "$", add a space to fix this. for c, label in self.labels_dict.items(): - mongo_key = c if not c.startswith("$") else f" {c}" + mongo_key = f" {c}" if c.startswith("$") else c dct["labels_dict"][mongo_key] = label.as_dict()["fcoords"] dct["projections"] = {} - if len(self.projections) != 0: + if len(self.projections) != 0 and self.structure is not None: dct["structure"] = self.structure.as_dict() dct["projections"] = {str(int(spin)): np.array(v).tolist() for spin, v in self.projections.items()} return dct @classmethod - def from_dict(cls, dct: dict) -> Self: - """Create from dict. + def from_dict(cls, dct: dict[str, Any]) -> Self: + """Create from a dict. Args: - dct: A dict with all data for a band structure object. + dct: A dict with all data for a BandStructure. Returns: - A BandStructure object + A BandStructure object. """ # Strip the label to recover initial string - # (see trick used in as_dict to handle $ chars) + # (see trick used in as_dict to handle "$"" chars) labels_dict = {k.strip(): v for k, v in dct["labels_dict"].items()} - projections = {} - structure = None + if isinstance(next(iter(dct["bands"].values())), dict): eigenvals = {Spin(int(k)): np.array(dct["bands"][k]["data"]) for k in dct["bands"]} else: eigenvals = {Spin(int(k)): dct["bands"][k] for k in dct["bands"]} + structure = None if "structure" in dct: structure = Structure.from_dict(dct["structure"]) + projections = {} try: if dct.get("projections"): if isinstance(dct["projections"]["1"][0][0], dict): @@ -638,15 +651,16 @@ def from_dict(cls, dct: dict) -> Self: return cls.from_old_dict(dct) @classmethod - def from_old_dict(cls, dct) -> Self: + def from_old_dict(cls, dct: dict[str, Any]) -> Self: """ Args: - dct (dict): A dict with all data for a band structure symmetry line object. + dct (dict): A dict with all data for a BandStructure object. Returns: - A BandStructureSymmLine object + A BandStructure object. """ - # Strip the label to recover initial string (see trick used in as_dict to handle $ chars) + # Strip the label to recover initial string + # (see trick used in as_dict to handle "$" chars) labels_dict = {k.strip(): v for k, v in dct["labels_dict"].items()} projections: dict = {} structure = None @@ -687,37 +701,35 @@ class BandStructureSymmLine(BandStructure, MSONable): def __init__( self, - kpoints, - eigenvals, - lattice, - efermi, - labels_dict, - coords_are_cartesian=False, - structure=None, - projections=None, + kpoints: NDArray, + eigenvals: dict[Spin, list], + lattice: Lattice, + efermi: float, + labels_dict: dict[str, Kpoint], + coords_are_cartesian: bool = False, + structure: Structure | None = None, + projections: dict[Spin, NDArray] | None = None, ) -> None: """ Args: - kpoints: list of kpoint as numpy arrays, in frac_coords of the + kpoints (NDArray): Array of kpoint, in frac_coords of the given lattice by default - eigenvals: dict of energies for spin up and spin down + eigenvals (dict[Spin, list]): Energies for spin up and spin down {Spin.up:[][],Spin.down:[][]}, the first index of the array [][] refers to the band and the second to the index of the kpoint. The kpoints are ordered according to the order of the kpoints array. If the band structure is not spin polarized, we only store one data set under Spin.up. - lattice: The reciprocal lattice. - Pymatgen uses the physics convention of reciprocal lattice vectors - WITH a 2*pi coefficient - efermi: fermi energy - labels_dict: (dict) of {} this link a kpoint (in frac coords or - Cartesian coordinates depending on the coords). - coords_are_cartesian: Whether coordinates are cartesian. - structure: The crystal structure (as a pymatgen Structure object) - associated with the band structure. This is needed if we - provide projections to the band structure. - projections: dict of orbital projections as {spin: array}. The - indices of the array are [band_index, kpoint_index, orbital_index, + lattice (Lattice): The reciprocal lattice. Pymatgen uses the physics + convention of reciprocal lattice vectors with a 2*pi coefficient. + efermi (float): The Fermi level. + labels_dict (dict[str, Kpoint]): Dict mapping label to Kpoint. + coords_are_cartesian (bool): Whether coordinates are cartesian. + structure (Structure): The crystal structure associated with the + band structure. This is needed if we provide projections to + the band structure. + projections (dict[Spin, NDArray]): Orbital projections as {spin: array}. + The indices of the array are [band_index, kpoint_index, orbital_index, ion_index].If the band structure is not spin polarized, we only store one data set under Spin.up. """ @@ -735,7 +747,7 @@ def __init__( self.branches = [] one_group: list = [] branches_tmp = [] - # get labels and distance for each kpoint + # Get labels and distance for each kpoint previous_kpoint = self.kpoints[0] previous_distance = 0.0 @@ -770,44 +782,40 @@ def __init__( if len(self.bands) == 2: self.is_spin_polarized = True - def get_equivalent_kpoints(self, index): - """Get the list of kpoint indices equivalent (meaning they are the - same frac coords) to the given one. + def get_equivalent_kpoints(self, index: int) -> list[int]: + """Get kpoint indices equivalent (having the same coords) to the given one. Args: - index: the kpoint index + index (int): The kpoint index Returns: - a list of equivalent indices + list[int]: Equivalent indices. - TODO: now it uses the label we might want to use coordinates instead - (in case there was a mislabel) + TODO: now it uses the label, we might want to use coordinates + instead in case there was a mislabel. """ - # if the kpoint has no label it can't have a repetition along the band - # structure line object - + # If the kpoint has no label it can't have a repetition + # along the BandStructureSymmLine object if self.kpoints[index].label is None: return [index] list_index_kpoints = [] - for i, kpt in enumerate(self.kpoints): + for idx, kpt in enumerate(self.kpoints): if kpt.label == self.kpoints[index].label: - list_index_kpoints.append(i) + list_index_kpoints.append(idx) return list_index_kpoints - def get_branch(self, index): - r"""Get in what branch(es) is the kpoint. There can be several - branches. + def get_branch(self, index: int) -> list[dict[str, Any]]: + """Get what branch(es) is the kpoint. It takes into account the + fact that one kpoint (e.g., Gamma) can be in several branches. Args: - index: the kpoint index + index (int): The kpoint index. Returns: - A list of dictionaries [{"name","start_index","end_index","index"}] - indicating all branches in which the k_point is. It takes into - account the fact that one kpoint (e.g., \\Gamma) can be in several - branches + A list of dicts [{"name", "start_index", "end_index", "index"}] + indicating all branches in which the k_point is. """ to_return = [] for idx in self.get_equivalent_kpoints(index): @@ -823,19 +831,19 @@ def get_branch(self, index): ) return to_return - def apply_scissor(self, new_band_gap): + def apply_scissor(self, new_band_gap: float) -> Self: """Apply a scissor operator (shift of the CBM) to fit the given band gap. If it's a metal, we look for the band crossing the Fermi level and shift this one up. This will not work all the time for metals! Args: - new_band_gap: the band gap the scissor band structure need to have. + new_band_gap (float): The band gap the scissor band structure need to have. Returns: - BandStructureSymmLine: with the applied scissor shift + BandStructureSymmLine: With the applied scissor shift. """ if self.is_metal(): - # moves then the highest index band crossing the Fermi level find this band... + # Move then the highest index band crossing the Fermi level find this band... max_index = -1000 # spin_index = None for idx in range(self.nb_bands): @@ -867,6 +875,7 @@ def apply_scissor(self, new_band_gap): for v in range(len(old_dict["bands"][spin][k])): if k >= max_index: old_dict["bands"][spin][k][v] = old_dict["bands"][spin][k][v] + shift + else: shift = new_band_gap - self.get_band_gap()["energy"] old_dict = self.as_dict() @@ -876,9 +885,10 @@ def apply_scissor(self, new_band_gap): if old_dict["bands"][spin][k][v] >= old_dict["cbm"]["energy"]: old_dict["bands"][spin][k][v] = old_dict["bands"][spin][k][v] + shift old_dict["efermi"] = old_dict["efermi"] + shift + return self.from_dict(old_dict) - def as_dict(self): + def as_dict(self) -> dict[str, Any]: """JSON-serializable dict representation of BandStructureSymmLine.""" dct = super().as_dict() dct["branches"] = self.branches @@ -886,11 +896,11 @@ def as_dict(self): class LobsterBandStructureSymmLine(BandStructureSymmLine): - """Lobster subclass of BandStructure with customized functions.""" + """LOBSTER subclass of BandStructure with customized functions.""" - def as_dict(self): + def as_dict(self) -> dict[str, Any]: """JSON-serializable dict representation of BandStructureSymmLine.""" - dct = { + dct: dict[str, Any] = { "@module": type(self).__module__, "@class": type(self).__name__, "lattice_rec": self.lattice_rec.as_dict(), @@ -921,27 +931,28 @@ def as_dict(self): dct["band_gap"] = self.get_band_gap() dct["labels_dict"] = {} dct["is_spin_polarized"] = self.is_spin_polarized - # MongoDB does not accept keys starting with $. Add a blank space to fix the problem + + # MongoDB does not accept keys starting with "$", add a space to fix this. for c, label in self.labels_dict.items(): - mongo_key = c if not c.startswith("$") else " " + c + mongo_key = f" {c}" if c.startswith("$") else c dct["labels_dict"][mongo_key] = label.as_dict()["fcoords"] - if len(self.projections) != 0: + if len(self.projections) != 0 and self.structure is not None: dct["structure"] = self.structure.as_dict() dct["projections"] = {str(int(spin)): np.array(v).tolist() for spin, v in self.projections.items()} return dct @classmethod - def from_dict(cls, dct: dict) -> Self: + def from_dict(cls, dct: dict[str, Any]) -> Self: """ Args: - dct (dict): A dict with all data for a band structure symmetry line - object. + dct (dict): All data for a LobsterBandStructureSymmLine object. Returns: - A BandStructureSymmLine object + A LobsterBandStructureSymmLine object. """ try: - # Strip the label to recover initial string (see trick used in as_dict to handle $ chars) + # Strip the label to recover initial string + # (see trick used in as_dict to handle "$" chars) labels_dict = {k.strip(): v for k, v in dct["labels_dict"].items()} projections = {} structure = None @@ -970,16 +981,16 @@ def from_dict(cls, dct: dict) -> Self: return cls.from_old_dict(dct) @classmethod - def from_old_dict(cls, dct) -> Self: + def from_old_dict(cls, dct: dict[str, Any]) -> Self: """ Args: - dct (dict): A dict with all data for a band structure symmetry line - object. + dct (dict): All data for a LobsterBandStructureSymmLine object. Returns: - A BandStructureSymmLine object + A LobsterBandStructureSymmLine object """ - # Strip the label to recover initial string (see trick used in as_dict to handle $ chars) + # Strip the label to recover initial string + # (see trick used in as_dict to handle "$" chars) labels_dict = {k.strip(): v for k, v in dct["labels_dict"].items()} projections: dict = {} structure = None @@ -1005,19 +1016,18 @@ def from_old_dict(cls, dct) -> Self: projections=projections, ) - def get_projection_on_elements(self): - """Get a dictionary of projections on elements. It sums over all available orbitals + def get_projection_on_elements(self) -> dict[Spin, list]: + """Get projections on elements. It sums over all available orbitals for each element. Returns: - a dictionary in the {Spin.up:[][{Element:values}], - Spin.down:[][{Element:values}]} format - if there is no projections in the band structure - returns an empty dict + dict[Spin, list]: dict in the {Spin.up:[][{Element:values}], + Spin.down:[][{Element:values}]} format. + If there is no projections in the band structure, return {}. """ - result = {} + result: dict[Spin, list] = {} for spin, v in self.projections.items(): - result[spin] = [[defaultdict(float) for i in range(len(self.kpoints))] for j in range(self.nb_bands)] + result[spin] = [[defaultdict(float) for _ in range(len(self.kpoints))] for _ in range(self.nb_bands)] for i, j in itertools.product(range(self.nb_bands), range(len(self.kpoints))): for key, item in v[i][j].items(): for item2 in item.values(): @@ -1025,13 +1035,17 @@ def get_projection_on_elements(self): result[spin][i][j][specie] += item2 return result - def get_projections_on_elements_and_orbitals(self, el_orb_spec): - """Return a dictionary of projections on elements and specific orbitals. + def get_projections_on_elements_and_orbitals( + self, + el_orb_spec: dict[Element, list], # type: ignore[override] + ) -> dict[Spin, list]: + """Get projections on elements and specific orbitals. Args: - el_orb_spec: A dictionary of Elements and Orbitals for which we want - to have projections on. It is given as: {Element:[orbitals]}, - e.g. {'Si':['3s','3p']} or {'Si':['3s','3p_x', '3p_y', '3p_z']} depending on input files + el_orb_spec (dict): Elements and Orbitals for which we want + to project on. It is given as {Element: [orbitals]}, + e.g. {"Si": ["3s", "3p"]} or {"Si": ["3s", "3p_x", "3p_y", "3p_z']} + depending on input files. Returns: A dictionary of projections on elements in the @@ -1040,12 +1054,12 @@ def get_projections_on_elements_and_orbitals(self, el_orb_spec): if there is no projections in the band structure returns an empty dict. """ - result = {} + result: dict[Spin, list] = {} el_orb_spec = {get_el_sp(el): orbs for el, orbs in el_orb_spec.items()} for spin, v in self.projections.items(): result[spin] = [ - [{str(e): defaultdict(float) for e in el_orb_spec} for i in range(len(self.kpoints))] - for j in range(self.nb_bands) + [{str(e): defaultdict(float) for e in el_orb_spec} for _ in range(len(self.kpoints))] + for _ in range(self.nb_bands) ] for i, j in itertools.product(range(self.nb_bands), range(len(self.kpoints))): @@ -1057,39 +1071,52 @@ def get_projections_on_elements_and_orbitals(self, el_orb_spec): return result -def get_reconstructed_band_structure(list_bs, efermi=None): - """Take a list of band structures and reconstructs - one band structure object from all of them. +@overload +def get_reconstructed_band_structure( # type: ignore[overload-overlap] + list_bs: list[BandStructure], + efermi: float | None = None, +) -> BandStructure: + pass - This is typically very useful when you split non self consistent - band structure runs in several independent jobs and want to merge back - the results + +@overload +def get_reconstructed_band_structure( + list_bs: list[BandStructureSymmLine], + efermi: float | None = None, +) -> BandStructureSymmLine: + pass + + +def get_reconstructed_band_structure( + list_bs: list[BandStructure] | list[BandStructureSymmLine], + efermi: float | None = None, +) -> BandStructure | BandStructureSymmLine: + """Merge multiple BandStructure(SymmLine) objects to a single one. + + This is typically useful when you split non self-consistent band + structure runs to several independent jobs and want to merge the results. Args: - list_bs: A list of BandStructure or BandStructureSymmLine objects. - efermi: The Fermi energy of the reconstructed band structure. If - None is assigned an average of all the Fermi energy in each + list_bs (list): BandStructure or BandStructureSymmLine objects. + efermi (float): The Fermi level of the reconstructed band structure. + If None, an average of all the Fermi levels in each object in the list_bs is used. Returns: A BandStructure or BandStructureSymmLine object (depending on - the type of the list_bs objects) + the type of the objects in list_bs). """ if efermi is None: efermi = sum(b.efermi for b in list_bs) / len(list_bs) - kpoints = [] - labels_dict = {} rec_lattice = list_bs[0].lattice_rec nb_bands = min(list_bs[i].nb_bands for i in range(len(list_bs))) - kpoints = np.concatenate([[k.frac_coords for k in bs.kpoints] for bs in list_bs]) + kpoints = np.concatenate([[kpt.frac_coords for kpt in bs.kpoints] for bs in list_bs]) dicts = [bs.labels_dict for bs in list_bs] - labels_dict = {k: v.frac_coords for d in dicts for k, v in d.items()} - - eigenvals = {} - eigenvals[Spin.up] = np.concatenate([bs.bands[Spin.up][:nb_bands] for bs in list_bs], axis=1) + labels_dict = {key: val.frac_coords for dct in dicts for key, val in dct.items()} + eigenvals = {Spin.up: np.concatenate([bs.bands[Spin.up][:nb_bands] for bs in list_bs], axis=1)} if list_bs[0].is_spin_polarized: eigenvals[Spin.down] = np.concatenate([bs.bands[Spin.down][:nb_bands] for bs in list_bs], axis=1) diff --git a/src/pymatgen/electronic_structure/boltztrap2.py b/src/pymatgen/electronic_structure/boltztrap2.py index d3ed99d5167..13d7891629e 100644 --- a/src/pymatgen/electronic_structure/boltztrap2.py +++ b/src/pymatgen/electronic_structure/boltztrap2.py @@ -1,29 +1,27 @@ -"""BoltzTraP2 is a python software interpolating band structures and -computing materials properties from dft band structure using Boltzmann -semi-classical transport theory. -This module provides a pymatgen interface to BoltzTraP2. +"""This module provides an interface to BoltzTraP2. Some of the code is written following the examples provided in BoltzTraP2. -BoltzTraP2 has been developed by Georg Madsen, Jesús Carrete, Matthieu J. Verstraete. +BoltzTraP2 is a Python software interpolating band structures and +computing materials properties from DFT band structure using Boltzmann +semi-classical transport theory, developed by Georg Madsen, Jesús Carrete, +Matthieu J. Verstraete. https://gitlab.com/sousaw/BoltzTraP2 https://www.sciencedirect.com/science/article/pii/S0010465518301632 -References are: - +References: Georg K.H.Madsen, Jesús Carrete, Matthieu J.Verstraete BoltzTraP2, a program for interpolating band structures and calculating semi-classical transport coefficients - Computer Physics Communications 231, 140-145, 2018 + Computer Physics Communications 231, 140-145, 2018. Madsen, G. K. H., and Singh, D. J. (2006). BoltzTraP. A code for calculating band-structure dependent quantities. - Computer Physics Communications, 175, 67-71 + Computer Physics Communications, 175, 67-71. Todo: -- DONE: spin polarized bands -- read first derivative of the eigenvalues from vasprun.xml (mommat) -- handle magnetic moments (magmom) +- Read first derivative of the eigenvalues from vasprun.xml (mommat) +- Handle magnetic moments (MAGMOM) """ from __future__ import annotations diff --git a/src/pymatgen/electronic_structure/cohp.py b/src/pymatgen/electronic_structure/cohp.py index b641d2a3992..acf4f8f45b5 100644 --- a/src/pymatgen/electronic_structure/cohp.py +++ b/src/pymatgen/electronic_structure/cohp.py @@ -1,6 +1,8 @@ -"""This module defines classes to represent crystal orbital Hamilton -populations (COHP) and integrated COHP (ICOHP), but can also be used -for crystal orbital overlap populations (COOP) or crystal orbital bond indices (COBIs). +"""This module defines classes to represent: + - Crystal orbital Hamilton population (COHP) and integrated COHP (ICOHP). + - Crystal orbital overlap population (COOP). + - Crystal orbital bond index (COBI). + If you use this module, please cite: J. George, G. Petretto, A. Naik, M. Esters, A. J. Jackson, R. Nelson, R. Dronskowski, G.-M. Rignanese, G. Hautier, "Automated Bonding Analysis with Crystal Orbital Hamilton Populations", @@ -29,10 +31,14 @@ from pymatgen.util.num import round_to_sigfigs if TYPE_CHECKING: - from typing import Any + from collections.abc import Sequence + from typing import Any, Literal + from numpy.typing import NDArray from typing_extensions import Self + from pymatgen.util.typing import PathLike, SpinLike, Vector3D + __author__ = "Marco Esters, Janine George" __copyright__ = "Copyright 2017, The Materials Project" __version__ = "0.2" @@ -50,17 +56,24 @@ class Cohp(MSONable): """Basic COHP object.""" def __init__( - self, efermi, energies, cohp, are_coops=False, are_cobis=False, are_multi_center_cobis=False, icohp=None + self, + efermi: float, + energies: Sequence[float], + cohp: dict[Spin, NDArray], + are_coops: bool = False, + are_cobis: bool = False, + are_multi_center_cobis: bool = False, + icohp: dict[Spin, NDArray] | None = None, ) -> None: """ Args: - are_coops: Indicates whether this object describes COOPs. - are_cobis: Indicates whether this object describes COBIs. - are_multi_center_cobis: Indicates whether this object describes multi-center COBIs - efermi: Fermi energy. - energies: A sequence of energies. - cohp ({Spin: np.array}): representing the COHP for each spin. - icohp ({Spin: np.array}): representing the ICOHP for each spin. + efermi (float): The Fermi level. + energies (Sequence[float]): Energies. + cohp ({Spin: NDArrary}): The COHP for each spin. + are_coops (bool): Whether this object describes COOPs. + are_cobis (bool): Whether this object describes COBIs. + are_multi_center_cobis (bool): Whether this object describes multi-center COBIs. + icohp ({Spin: NDArrary}): The ICOHP for each spin. """ self.are_coops = are_coops self.are_cobis = are_cobis @@ -71,7 +84,7 @@ def __init__( self.icohp = icohp def __repr__(self) -> str: - """Get a string that can be easily plotted (e.g. using gnuplot).""" + """A string that can be easily plotted (e.g. using gnuplot).""" if self.are_coops: cohp_str = "COOP" elif self.are_cobis or self.are_multi_center_cobis: @@ -97,7 +110,7 @@ def __repr__(self) -> str: str_arr.append(format_data.format(*(d[idx] for d in data))) return "\n".join(str_arr) - def as_dict(self): + def as_dict(self) -> dict[str, Any]: """JSON-serializable dict representation of COHP.""" dct = { "@module": type(self).__module__, @@ -113,20 +126,22 @@ def as_dict(self): dct["ICOHP"] = {str(spin): pops.tolist() for spin, pops in self.icohp.items()} return dct - def get_cohp(self, spin=None, integrated=False): + def get_cohp( + self, + spin: SpinLike | None = None, + integrated: bool = False, + ) -> dict[Spin, NDArray] | None: """Get the COHP or ICOHP for a particular spin. Args: - spin: Spin. Can be parsed as spin object, integer (-1/1) - or str ("up"/"down") - integrated: Return COHP (False) or ICOHP (True) + spin (SpinLike): Selected spin. If is None and both + spins are present, both will be returned. + integrated: Return ICOHP (True) or COHP (False). Returns: - Returns the CHOP or ICOHP for the input spin. If Spin is - None and both spins are present, both spins will be returned - as a dictionary. + dict: The COHP or ICOHP for the selected spin. """ - populations = self.cohp if not integrated else self.icohp + populations = self.icohp if integrated else self.cohp if populations is None: return None @@ -138,118 +153,127 @@ def get_cohp(self, spin=None, integrated=False): spin = Spin({"up": 1, "down": -1}[spin.lower()]) return {spin: populations[spin]} - def get_icohp(self, spin=None): - """Convenient alternative to get the ICOHP for a particular spin.""" + def get_icohp( + self, + spin: SpinLike | None = None, + ) -> dict[Spin, NDArray] | None: + """Convenient wrapper to get the ICOHP for a particular spin.""" return self.get_cohp(spin=spin, integrated=True) - def get_interpolated_value(self, energy, integrated=False): - """Get the COHP for a particular energy. + def get_interpolated_value( + self, + energy: float, + integrated: bool = False, + ) -> dict[Spin, float]: + """Get the interpolated COHP for a particular energy. Args: - energy: Energy to return the COHP value for. - integrated: Return COHP (False) or ICOHP (True) + energy (float): Energy to get the COHP value for. + integrated (bool): Return ICOHP (True) or COHP (False). """ - inter = {} + inters = {} for spin in self.cohp: if not integrated: - inter[spin] = get_linear_interpolated_value(self.energies, self.cohp[spin], energy) + inters[spin] = get_linear_interpolated_value(self.energies, self.cohp[spin], energy) elif self.icohp is not None: - inter[spin] = get_linear_interpolated_value(self.energies, self.icohp[spin], energy) + inters[spin] = get_linear_interpolated_value(self.energies, self.icohp[spin], energy) else: raise ValueError("ICOHP is empty.") - return inter + return inters + + def has_antibnd_states_below_efermi( + self, + spin: SpinLike | None = None, + limit: float = 0.01, + ) -> dict[Spin, bool] | None: + """Get dict of antibonding states below the Fermi level for the spin. - def has_antibnd_states_below_efermi(self, spin=None, limit=0.01): - """Get dict indicating if there are antibonding states below the Fermi level depending on the spin - spin: Spin - limit: -COHP smaller -limit will be considered. + Args: + spin (SpinLike): Selected spin. + limit (float): Only COHP higher than this value will be considered. """ populations = self.cohp n_energies_below_efermi = len([energy for energy in self.energies if energy <= self.efermi]) if populations is None: return None + + dict_to_return = {} if spin is None: - dict_to_return = {} for sp, cohp_vals in populations.items(): - if (max(cohp_vals[:n_energies_below_efermi])) > limit: - dict_to_return[sp] = True - else: - dict_to_return[sp] = False + # NOTE: Casting to bool is necessary, otherwise ended up + # getting "bool_" instead of "bool" from NumPy + dict_to_return[sp] = bool((max(cohp_vals[:n_energies_below_efermi])) > limit) + else: - dict_to_return = {} if isinstance(spin, int): spin = Spin(spin) elif isinstance(spin, str): spin = Spin({"up": 1, "down": -1}[spin.lower()]) - if (max(populations[spin][:n_energies_below_efermi])) > limit: - dict_to_return[spin] = True - else: - dict_to_return[spin] = False + dict_to_return[spin] = bool((max(populations[spin][:n_energies_below_efermi])) > limit) return dict_to_return @classmethod def from_dict(cls, dct: dict[str, Any]) -> Self: - """Get a COHP object from a dict representation of the COHP.""" + """Generate Cohp from a dict representation.""" icohp = {Spin(int(key)): np.array(val) for key, val in dct["ICOHP"].items()} if "ICOHP" in dct else None - are_cobis = dct.get("are_cobis", False) - are_multi_center_cobis = dct.get("are_multi_center_cobis", False) + return cls( dct["efermi"], dct["energies"], {Spin(int(key)): np.array(val) for key, val in dct["COHP"].items()}, icohp=icohp, are_coops=dct["are_coops"], - are_cobis=are_cobis, - are_multi_center_cobis=are_multi_center_cobis, + are_cobis=dct.get("are_cobis", False), + are_multi_center_cobis=dct.get("are_multi_center_cobis", False), ) class CompleteCohp(Cohp): - """A wrapper class that defines an average COHP, and individual COHPs. + """A wrapper that defines an average COHP, and individual COHPs. Attributes: - are_coops (bool): Indicates whether the object is consisting of COOPs. - are_cobis (bool): Indicates whether the object is consisting of COBIs. - efermi (float): Fermi energy. + are_coops (bool): Whether the object is consisting of COOPs. + are_cobis (bool): Whether the object is consisting of COBIs. + efermi (float): The Fermi level. energies (Sequence[float]): Sequence of energies. - structure (pymatgen.Structure): Structure associated with the COHPs. + structure (Structure): Structure associated with the COHPs. cohp (Sequence[float]): The average COHP. icohp (Sequence[float]): The average ICOHP. - all_cohps (dict[str, Sequence[float]]): A dict of COHPs for individual bonds of the form {label: COHP}. + all_cohps (dict[str, Sequence[float]]): COHPs for individual bonds of the form {label: COHP}. orb_res_cohp (dict[str, Dict[str, Sequence[float]]]): Orbital-resolved COHPs. """ def __init__( self, - structure, - avg_cohp, - cohp_dict, - bonds=None, - are_coops=False, - are_cobis=False, - are_multi_center_cobis=False, - orb_res_cohp=None, + structure: Structure, + avg_cohp: Cohp, + cohp_dict: dict[str, Cohp], + bonds: dict[str, Any] | None = None, + are_coops: bool = False, + are_cobis: bool = False, + are_multi_center_cobis: bool = False, + orb_res_cohp: dict[str, dict] | None = None, ) -> None: """ Args: - structure: Structure associated with this COHP. - avg_cohp: The average cohp as a COHP object. - cohp_dict: A dict of COHP objects for individual bonds of the form - {label: COHP} - bonds: A dict containing information on the bonds of the form - {label: {key: val}}. The key-val pair can be any information - the user wants to put in, but typically contains the sites, - the bond length, and the number of bonds. If nothing is + structure (Structure): Structure associated with this COHP. + avg_cohp (Cohp): The average COHP. + cohp_dict (dict[str, Cohp]): COHP for individual bonds of the form + {label: COHP}. + bonds (dict[str, Any]): Information on the bonds of the form + {label: {key: val}}. The value can be any information, + but typically contains the sites, the bond length, + and the number of bonds. If nothing is supplied, it will default to an empty dict. - are_coops: indicates whether the Cohp objects are COOPs. + are_coops (bool): Whether the Cohp objects are COOPs. Defaults to False for COHPs. - are_cobis: indicates whether the Cohp objects are COBIs. + are_cobis (bool): Whether the Cohp objects are COBIs. Defaults to False for COHPs. - are_multi_center_cobis: indicates whether the Cohp objects are multi-center COBIs. + are_multi_center_cobis (bool): Whether the Cohp objects are multi-center COBIs. Defaults to False for COHPs. - orb_res_cohp: Orbital-resolved COHPs. + orb_res_cohp (dict): Orbital-resolved COHPs. """ if ( (are_coops and are_cobis) @@ -257,6 +281,7 @@ def __init__( or (are_cobis and are_multi_center_cobis) ): raise ValueError("You cannot have info about COOPs, COBIs and/or multi-center COBIS in the same file.") + super().__init__( avg_cohp.efermi, avg_cohp.energies, @@ -276,12 +301,15 @@ def __init__( def __str__(self) -> str: if self.are_coops: - return f"Complete COOPs for {self.structure}" - if self.are_cobis: - return f"Complete COBIs for {self.structure}" - return f"Complete COHPs for {self.structure}" + header = "COOPs" + elif self.are_cobis: + header = "COBIs" + else: + header = "COHPs" + + return f"Complete {header} for {self.structure}" - def as_dict(self): + def as_dict(self) -> dict[str, Any]: """JSON-serializable dict representation of CompleteCohp.""" dct = { "@module": type(self).__module__, @@ -299,16 +327,14 @@ def as_dict(self): dct["ICOHP"] = {"average": {str(spin): pops.tolist() for spin, pops in self.icohp.items()}} for label in self.all_cohps: - dct["COHP"].update({label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].cohp.items()}}) - if self.all_cohps[label].icohp is not None: + dct["COHP"] |= {label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].cohp.items()}} + icohp = self.all_cohps[label].icohp + if icohp is not None: if "ICOHP" not in dct: - dct["ICOHP"] = { - label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].icohp.items()} - } + dct["ICOHP"] = {label: {str(spin): pops.tolist() for spin, pops in icohp.items()}} else: - dct["ICOHP"].update( - {label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].icohp.items()}} - ) + dct["ICOHP"] |= {label: {str(spin): pops.tolist() for spin, pops in icohp.items()}} + if False in [bond_dict == {} for bond_dict in self.bonds.values()]: dct["bonds"] = { bond: { @@ -317,44 +343,55 @@ def as_dict(self): } for bond in self.bonds } + if self.orb_res_cohp: - orb_dict = {} + orb_dict: dict[str, Any] = {} for label in self.orb_res_cohp: orb_dict[label] = {} for orbs in self.orb_res_cohp[label]: - cohp = {str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["COHP"].items()} - orb_dict[label][orbs] = {"COHP": cohp} - icohp = {str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["ICOHP"].items()} - orb_dict[label][orbs]["ICOHP"] = icohp - orbitals = [[orb[0], orb[1].name] for orb in self.orb_res_cohp[label][orbs]["orbitals"]] - orb_dict[label][orbs]["orbitals"] = orbitals + orb_dict[label][orbs] = { + "COHP": { + str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["COHP"].items() + }, + "ICOHP": { + str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["ICOHP"].items() + }, + "orbitals": [[orb[0], orb[1].name] for orb in self.orb_res_cohp[label][orbs]["orbitals"]], + } + dct["orb_res_cohp"] = orb_dict return dct - def get_cohp_by_label(self, label, summed_spin_channels=False): - """Get specific COHP object. + def get_cohp_by_label( + self, + label: str, + summed_spin_channels: bool = False, + ) -> Cohp: + """Get specific Cohp by the label, to simplify plotting. Args: - label: string (for newer Lobster versions: a number) - summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true + label (str): Label for the interaction. + summed_spin_channels (bool): Sum the spin channels and return the sum as Spin.up. Returns: - Returns the COHP object to simplify plotting + The Cohp. """ if label.lower() == "average": - divided_cohp = self.cohp - divided_icohp = self.icohp - + divided_cohp: dict[Spin, Any] | None = self.cohp + divided_icohp: dict[Spin, Any] | None = self.icohp else: divided_cohp = self.all_cohps[label].get_cohp(spin=None, integrated=False) divided_icohp = self.all_cohps[label].get_icohp(spin=None) + assert divided_cohp is not None + if summed_spin_channels and Spin.down in self.cohp: - final_cohp = {} - final_icohp = {} - final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) - final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) + assert divided_icohp is not None + final_cohp: dict[Spin, Any] = {Spin.up: np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0)} + final_icohp: dict[Spin, Any] | None = { + Spin.up: np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) + } else: final_cohp = divided_cohp final_icohp = divided_icohp @@ -368,46 +405,50 @@ def get_cohp_by_label(self, label, summed_spin_channels=False): icohp=final_icohp, ) - def get_summed_cohp_by_label_list(self, label_list, divisor=1, summed_spin_channels=False): - """Get a COHP object that includes a summed COHP divided by divisor. + def get_summed_cohp_by_label_list( + self, + label_list: list[str], + divisor: float = 1, + summed_spin_channels: bool = False, + ) -> Cohp: + """Get a Cohp object that includes a summed COHP divided by divisor. Args: - label_list: list of labels for the COHP that should be included in the summed cohp - divisor: float/int, the summed cohp will be divided by this divisor - summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true + label_list (list[str]): Labels for the COHP to include. + divisor (float): The summed COHP will be divided by this divisor. + summed_spin_channels (bool): Sum the spin channels and return the sum in Spin.up. Returns: - Returns a COHP object including a summed COHP + A Cohp object for the summed COHP. """ - # check if cohps are spinpolarized or not + # Check if COHPs are spin polarized first_cohpobject = self.get_cohp_by_label(label_list[0]) summed_cohp = first_cohpobject.cohp.copy() + assert first_cohpobject.icohp is not None summed_icohp = first_cohpobject.icohp.copy() for label in label_list[1:]: - cohp_here = self.get_cohp_by_label(label) - summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp_here.cohp[Spin.up]], axis=0) + cohp = self.get_cohp_by_label(label) + icohp = cohp.icohp + assert icohp is not None + summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp.cohp[Spin.up]], axis=0) if Spin.down in summed_cohp: - summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp_here.cohp[Spin.down]], axis=0) + summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp.cohp[Spin.down]], axis=0) - summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp_here.icohp[Spin.up]], axis=0) + summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], icohp[Spin.up]], axis=0) if Spin.down in summed_icohp: - summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.icohp[Spin.down]], axis=0) + summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], icohp[Spin.down]], axis=0) - divided_cohp = {} - divided_icohp = {} - divided_cohp[Spin.up] = np.divide(summed_cohp[Spin.up], divisor) - divided_icohp[Spin.up] = np.divide(summed_icohp[Spin.up], divisor) + divided_cohp = {Spin.up: np.divide(summed_cohp[Spin.up], divisor)} + divided_icohp = {Spin.up: np.divide(summed_icohp[Spin.up], divisor)} if Spin.down in summed_cohp: divided_cohp[Spin.down] = np.divide(summed_cohp[Spin.down], divisor) divided_icohp[Spin.down] = np.divide(summed_icohp[Spin.down], divisor) if summed_spin_channels and Spin.down in summed_cohp: - final_cohp = {} - final_icohp = {} - final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) - final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) + final_cohp = {Spin.up: np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0)} + final_icohp = {Spin.up: np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0)} else: final_cohp = divided_cohp final_icohp = divided_icohp @@ -422,50 +463,56 @@ def get_summed_cohp_by_label_list(self, label_list, divisor=1, summed_spin_chann ) def get_summed_cohp_by_label_and_orbital_list( - self, label_list, orbital_list, divisor=1, summed_spin_channels=False - ): - """Get a COHP object that includes a summed COHP divided by divisor. + self, + label_list: list[str], + orbital_list: list[str], + divisor: float = 1, + summed_spin_channels: bool = False, + ) -> Cohp: + """Get a Cohp object that includes a summed COHP divided by divisor. Args: - label_list: list of labels for the COHP that should be included in the summed cohp - orbital_list: list of orbitals for the COHPs that should be included in the summed cohp (same order as - label_list) - divisor: float/int, the summed cohp will be divided by this divisor - summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true + label_list (list[str]): Labels for the COHP that should be included. + orbital_list (list[str]): Orbitals for the COHPs that should be included + (same order as label_list). + divisor (float): The summed COHP will be divided by this divisor. + summed_spin_channels (bool): Sum the spin channels and return the sum in Spin.up. Returns: - Returns a COHP object including a summed COHP + A Cohp object including the summed COHP. """ - # check length of label_list and orbital_list: + # Check length of label_list and orbital_list if not len(label_list) == len(orbital_list): raise ValueError("label_list and orbital_list don't have the same length!") - # check if cohps are spinpolarized or not + + # Check if COHPs are spin polarized first_cohpobject = self.get_orbital_resolved_cohp(label_list[0], orbital_list[0]) + assert first_cohpobject is not None + assert first_cohpobject.icohp is not None summed_cohp = first_cohpobject.cohp.copy() summed_icohp = first_cohpobject.icohp.copy() - for ilabel, label in enumerate(label_list[1:], start=1): - cohp_here = self.get_orbital_resolved_cohp(label, orbital_list[ilabel]) - summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp_here.cohp.copy()[Spin.up]], axis=0) + + for idx, label in enumerate(label_list[1:], start=1): + cohp = self.get_orbital_resolved_cohp(label, orbital_list[idx]) + assert cohp is not None + assert cohp.icohp is not None + summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp.cohp.copy()[Spin.up]], axis=0) if Spin.down in summed_cohp: - summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp_here.cohp.copy()[Spin.down]], axis=0) - summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp_here.icohp.copy()[Spin.up]], axis=0) + summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp.cohp.copy()[Spin.down]], axis=0) + + summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp.icohp.copy()[Spin.up]], axis=0) if Spin.down in summed_icohp: - summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.icohp.copy()[Spin.down]], axis=0) + summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp.icohp.copy()[Spin.down]], axis=0) - divided_cohp = {} - divided_icohp = {} - divided_cohp[Spin.up] = np.divide(summed_cohp[Spin.up], divisor) - divided_icohp[Spin.up] = np.divide(summed_icohp[Spin.up], divisor) + divided_cohp = {Spin.up: np.divide(summed_cohp[Spin.up], divisor)} + divided_icohp = {Spin.up: np.divide(summed_icohp[Spin.up], divisor)} if Spin.down in summed_cohp: divided_cohp[Spin.down] = np.divide(summed_cohp[Spin.down], divisor) divided_icohp[Spin.down] = np.divide(summed_icohp[Spin.down], divisor) if summed_spin_channels and Spin.down in divided_cohp: - final_cohp = {} - final_icohp = {} - - final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) - final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) + final_cohp = {Spin.up: np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0)} + final_icohp = {Spin.up: np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0)} else: final_cohp = divided_cohp final_icohp = divided_icohp @@ -479,30 +526,34 @@ def get_summed_cohp_by_label_and_orbital_list( icohp=final_icohp, ) - def get_orbital_resolved_cohp(self, label, orbitals, summed_spin_channels=False): + def get_orbital_resolved_cohp( + self, + label: str, + orbitals: str | list[tuple[str, Orbital]] | tuple[tuple[str, Orbital], ...], + summed_spin_channels: bool = False, + ) -> Cohp | None: """Get orbital-resolved COHP. Args: - label: bond label (Lobster: labels as in ICOHPLIST/ICOOPLIST.lobster). - - orbitals: The orbitals as a label, or list or tuple of the form - [(n1, orbital1), (n2, orbital2)]. Orbitals can either be str, - int, or Orbital. - - summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true + label (str): Bond labels as in ICOHPLIST/ICOOPLIST.lobster. + orbitals: The orbitals as a label, or list/tuple of + [(n1, orbital1), (n2, orbital2), ...]. + Where each orbital can either be str, int, or Orbital. + summed_spin_channels (bool): Sum the spin channels and return the sum as Spin.up. Returns: - A Cohp object if CompleteCohp contains orbital-resolved cohp, - or None if it doesn't. + A Cohp object if CompleteCohp contains orbital-resolved COHP, + or None if it doesn't. Note: It currently assumes that orbitals are str if they aren't the - other valid types. This is not ideal, but the easiest way to - avoid unicode issues between python 2 and python 3. + other valid types. This is not ideal, but is the easiest way to + avoid unicode issues between Python 2 and Python 3. """ if self.orb_res_cohp is None: return None + if isinstance(orbitals, (list, tuple)): - cohp_orbs = [d["orbitals"] for d in self.orb_res_cohp[label].values()] + cohp_orbs = [val["orbitals"] for val in self.orb_res_cohp[label].values()] orbs = [] for orbital in orbitals: if isinstance(orbital[1], int): @@ -515,10 +566,12 @@ def get_orbital_resolved_cohp(self, label, orbitals, summed_spin_channels=False) raise TypeError("Orbital must be str, int, or Orbital.") orb_index = cohp_orbs.index(orbs) orb_label = list(self.orb_res_cohp[label])[orb_index] + elif isinstance(orbitals, str): orb_label = orbitals else: raise TypeError("Orbitals must be str, list, or tuple.") + try: icohp = self.orb_res_cohp[label][orb_label]["ICOHP"] except KeyError: @@ -547,9 +600,11 @@ def get_orbital_resolved_cohp(self, label, orbitals, summed_spin_channels=False) ) @classmethod - def from_dict(cls, dct: dict) -> Self: - """Get CompleteCohp object from dict representation.""" - # TODO: clean that mess up? + def from_dict(cls, dct: dict[str, Any]) -> Self: + """Get CompleteCohp object from a dict representation. + + TODO: Clean this up. + """ cohp_dict = {} efermi = dct["efermi"] energies = dct["energies"] @@ -570,6 +625,7 @@ def from_dict(cls, dct: dict) -> Self: } else: bonds = None + for label in dct["COHP"]: cohp = {Spin(int(spin)): np.array(dct["COHP"][label][spin]) for spin in dct["COHP"][label]} try: @@ -590,7 +646,7 @@ def from_dict(cls, dct: dict) -> Self: cohp_dict[label] = Cohp(efermi, energies, cohp, icohp=icohp) if "orb_res_cohp" in dct: - orb_cohp: dict[str, dict] = {} + orb_cohp: dict[str, dict[Orbital, dict[str, Any]]] = {} for label in dct["orb_res_cohp"]: orb_cohp[label] = {} for orb in dct["orb_res_cohp"][label]: @@ -612,9 +668,9 @@ def from_dict(cls, dct: dict) -> Self: "orbitals": orbitals, } # If no total COHPs are present, calculate the total - # COHPs from the single-orbital populations. Total COHPs - # may not be present when the COHP generator keyword is used - # in LOBSTER versions 2.2.0 and earlier. + # COHPs from the single-orbital populations. + # Total COHPs may not be present when the COHP generator keyword + # is used in LOBSTER versions 2.2.0 and earlier. if label not in dct["COHP"] or dct["COHP"][label] is None: cohp = { Spin.up: np.sum( @@ -648,52 +704,52 @@ def from_dict(cls, dct: dict) -> Self: else: orb_cohp = {} - are_cobis = dct.get("are_cobis", False) - + assert avg_cohp is not None return cls( structure, avg_cohp, cohp_dict, bonds=bonds, are_coops=dct["are_coops"], - are_cobis=are_cobis, + are_cobis=dct.get("are_cobis", False), are_multi_center_cobis=are_multi_center_cobis, orb_res_cohp=orb_cohp, ) @classmethod def from_file( - cls, fmt, filename=None, structure_file=None, are_coops=False, are_cobis=False, are_multi_center_cobis=False + cls, + fmt: Literal["LMTO", "LOBSTER"], + filename: PathLike | None = None, + structure_file: PathLike | None = None, + are_coops: bool = False, + are_cobis: bool = False, + are_multi_center_cobis: bool = False, ) -> Self: - """ - Creates a CompleteCohp object from an output file of a COHP - calculation. Valid formats are either LMTO (for the Stuttgart - LMTO-ASA code) or LOBSTER (for the LOBSTER code). + """Create CompleteCohp from an output file of a COHP calculation. Args: - fmt: A string for the code that was used to calculate - the COHPs so that the output file can be handled - correctly. Can take the values "LMTO" or "LOBSTER". - filename: Name of the COHP output file. Defaults to COPL - for LMTO and COHPCAR.lobster/COOPCAR.lobster for LOBSTER. - structure_file: Name of the file containing the structure. - If no file name is given, use CTRL for LMTO and POSCAR - for LOBSTER. - are_coops: Indicates whether the populations are COOPs or - COHPs. Defaults to False for COHPs. - are_cobis: Indicates whether the populations are COBIs or - COHPs. Defaults to False for COHPs. - are_multi_center_cobis: Indicates whether this file - includes information on multi-center COBIs + fmt (Literal["LMTO", "LOBSTER"]): The code used to calculate COHPs. + filename (PathLike): The COHP output file. Defaults to "COPL" + for LMTO and "COHPCAR.lobster/COOPCAR.lobster" for LOBSTER. + structure_file (PathLike): The file containing the structure. + If None, use "CTRL" for LMTO and "POSCAR" for LOBSTER. + are_coops (bool): Whether the populations are COOPs or COHPs. + Defaults to False for COHPs. + are_cobis (bool): Whether the populations are COBIs or COHPs. + Defaults to False for COHPs. + are_multi_center_cobis (bool): Whether this file + includes information on multi-center COBIs. Returns: A CompleteCohp object. """ if are_coops and are_cobis: raise ValueError("You cannot have info about COOPs and COBIs in the same file.") - fmt = fmt.upper() + + fmt = fmt.upper() # type: ignore[assignment] if fmt == "LMTO": - # LMTO COOPs and orbital-resolved COHP cannot be handled yet. + # TODO: LMTO COOPs and orbital-resolved COHP cannot be handled yet are_coops = False are_cobis = False orb_res_cohp = None @@ -701,7 +757,9 @@ def from_file( structure_file = "CTRL" if filename is None: filename = "COPL" + cohp_file: LMTOCopl | Cohpcar = LMTOCopl(filename=filename, to_eV=True) + elif fmt == "LOBSTER": if ( (are_coops and are_cobis) @@ -725,6 +783,7 @@ def from_file( are_multi_center_cobis=are_multi_center_cobis, ) orb_res_cohp = cohp_file.orb_res_cohp + else: raise ValueError(f"Unknown format {fmt}. Valid formats are LMTO and LOBSTER.") @@ -733,9 +792,8 @@ def from_file( cohp_data = cohp_file.cohp_data energies = cohp_file.energies - # Lobster shifts the energies so that the Fermi energy is at zero. + # LOBSTER shifts the energies so that the Fermi level is at zero. # Shifting should be done by the plotter object though. - spins = [Spin.up, Spin.down] if cohp_file.is_spin_polarized else [Spin.up] if fmt == "LOBSTER": energies += efermi @@ -745,6 +803,7 @@ def from_file( # COHPs from the single-orbital populations. Total COHPs # may not be present when the cohpgenerator keyword is used # in LOBSTER versions 2.2.0 and earlier. + # TODO: Test this more extensively for label in orb_res_cohp: @@ -766,16 +825,16 @@ def from_file( } if fmt == "LMTO": - # Calculate the average COHP for the LMTO file to be - # consistent with LOBSTER output. - avg_data: dict[str, dict] = {"COHP": {}, "ICOHP": {}} - for i in avg_data: + # Calculate the average COHP for the LMTO file to be consistent with LOBSTER + avg_data: dict[Literal["COHP", "ICOHP"], dict] = {"COHP": {}, "ICOHP": {}} + for dtype in avg_data: for spin in spins: - rows = np.array([v[i][spin] for v in cohp_data.values()]) + rows = np.array([v[dtype][spin] for v in cohp_data.values()]) avg = np.mean(rows, axis=0) - # LMTO COHPs have 5 significant figures - avg_data[i].update({spin: np.array([round_to_sigfigs(a, 5) for a in avg], dtype=float)}) + # LMTO COHPs have 5 significant digits + avg_data[dtype] |= {spin: np.array([round_to_sigfigs(a, 5) for a in avg], dtype=float)} avg_cohp = Cohp(efermi, energies, avg_data["COHP"], icohp=avg_data["ICOHP"]) + elif not are_multi_center_cobis: avg_cohp = Cohp( efermi, @@ -787,9 +846,9 @@ def from_file( are_multi_center_cobis=are_multi_center_cobis, ) del cohp_data["average"] + else: - # only include two-center cobis in average - # do this for both spin channels + # Only include two-center COBIs in average for both spin channels cohp = {} cohp[Spin.up] = np.array( [np.array(c["COHP"][Spin.up]) for c in cohp_file.cohp_data.values() if len(c["sites"]) <= 2] @@ -800,6 +859,7 @@ def from_file( ).mean(axis=0) except KeyError: pass + try: icohp = {} icohp[Spin.up] = np.array( @@ -813,6 +873,7 @@ def from_file( pass except KeyError: icohp = None + avg_cohp = Cohp( efermi, energies, @@ -856,38 +917,53 @@ def from_file( class IcohpValue(MSONable): - """Store information on an ICOHP or ICOOP value. + """Information for an ICOHP or ICOOP value. Attributes: - energies (ndarray): Energy values for the COHP/ICOHP/COOP/ICOOP. - densities (ndarray): Density of states values for the COHP/ICOHP/COOP/ICOOP. - energies_are_cartesian (bool): Whether the energies are cartesian or not. - are_coops (bool): Whether the object is a COOP/ICOOP or not. - are_cobis (bool): Whether the object is a COBIS/ICOBIS or not. - icohp (dict): A dictionary of the ICOHP/COHP values. The keys are Spin.up and Spin.down. + energies (NDArray): Energy values for the COHP/ICOHP/COOP/ICOOP. + densities (NDArray): Density of states for the COHP/ICOHP/COOP/ICOOP. + energies_are_cartesian (bool): Whether the energies are cartesian. + are_coops (bool): Whether the object is COOP/ICOOP. + are_cobis (bool): Whether the object is COBIS/ICOBIS. + icohp (dict): The ICOHP/COHP values, whose keys are Spin.up and Spin.down. summed_icohp (float): The summed ICOHP/COHP values. - num_bonds (int): The number of bonds used for the average COHP (relevant for Lobster versions <3.0). + num_bonds (int): The number of bonds used for the average COHP (for LOBSTER versions <3.0). """ def __init__( - self, label, atom1, atom2, length, translation, num, icohp, are_coops=False, are_cobis=False, orbitals=None + self, + label: str, + atom1: str, + atom2: str, + length: float, + translation: Vector3D, + num: int, + icohp: dict[Spin, float], + are_coops: bool = False, + are_cobis: bool = False, + orbitals: dict[str, dict[Literal["icohp", "orbitals"], Any]] | None = None, ) -> None: """ Args: - label: label for the icohp - atom1: str of atom that is contributing to the bond - atom2: str of second atom that is contributing to the bond - length: float of bond lengths - translation: translation list, e.g. [0,0,0] - num: integer describing how often the bond exists - icohp: dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} - are_coops: if True, this are COOPs - are_cobis: if True, this are COBIs - orbitals: {[str(Orbital1)-str(Orbital2)]: {"icohp":{Spin.up: icohpvalue for spin.up, Spin.down: - icohpvalue for spin.down}, "orbitals":[Orbital1, Orbital2]}}. + label (str): Label for the ICOHP. + atom1 (str): The first atom that contributes to the bond. + atom2 (str): The second atom that contributes to the bond. + length (float): Bond length. + translation (Vector3D): cell translation vector, e.g. (0, 0, 0). + num (int): The number of equivalent bonds. + icohp (dict[Spin, float]): {Spin.up: ICOHP_up, Spin.down: ICOHP_down} + are_coops (bool): Whether these are COOPs. + are_cobis (bool): Whether these are COBIs. + orbitals (dict): {[str(Orbital1)-str(Orbital2)]: { + "icohp": { + Spin.up: IcohpValue for spin.up, + Spin.down: IcohpValue for spin.down + }, + "orbitals": [Orbital1, Orbital2, ...]}. """ if are_coops and are_cobis: raise ValueError("You cannot have info about COOPs and COBIs in the same file.") + self._are_coops = are_coops self._are_cobis = are_cobis self._label = label @@ -898,136 +974,127 @@ def __init__( self._num = num self._icohp = icohp self._orbitals = orbitals - if Spin.down in self._icohp: - self._is_spin_polarized = True - else: - self._is_spin_polarized = False + self._is_spin_polarized = Spin.down in self._icohp def __str__(self) -> str: """String representation of the ICOHP/ICOOP.""" - if not self._are_coops and not self._are_cobis: - if self._is_spin_polarized: - return ( - f"ICOHP {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up) and {self._icohp[Spin.down]} eV (Spin down)" - ) - return ( - f"ICOHP {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up)" - ) - if self._are_coops and not self._are_cobis: - if self._is_spin_polarized: - return ( - f"ICOOP {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up) and {self._icohp[Spin.down]} eV (Spin down)" - ) - return ( - f"ICOOP {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up)" - ) - if self._are_cobis and not self._are_coops: - if self._is_spin_polarized: - return ( - f"ICOBI {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up) and {self._icohp[Spin.down]} eV (Spin down)" - ) + # (are_coops and are_cobis) is never True + if self._are_coops: + header = "ICOOP" + elif self._are_cobis: + header = "ICOBI" + else: + header = "ICOHP" + + if self._is_spin_polarized: return ( - f"ICOBI {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " - f"{self._icohp[Spin.up]} eV (Spin up)" + f"{header} {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " + f"{self._icohp[Spin.up]} eV (Spin up) and {self._icohp[Spin.down]} eV (Spin down)" ) - return "" + return ( + f"{header} {self._label} between {self._atom1} and {self._atom2} ({self._translation}): " + f"{self._icohp[Spin.up]} eV (Spin up)" + ) @property - def num_bonds(self): - """Tells the number of bonds for which the ICOHP value is an average. + def num_bonds(self) -> int: + """The number of bonds for which the ICOHP value is an average. Returns: - Int. + int """ return self._num @property def are_coops(self) -> bool: - """Tells if ICOOPs or not. + """Whether these are ICOOPs. Returns: - Boolean. + bool """ return self._are_coops @property def are_cobis(self) -> bool: - """Tells if ICOBIs or not. + """Whether these are ICOBIs. Returns: - Boolean. + bool """ return self._are_cobis @property def is_spin_polarized(self) -> bool: - """Tells if spin polarized calculation or not. + """Whether this is a spin polarized calculation. Returns: - Boolean. + bool """ return self._is_spin_polarized - def icohpvalue(self, spin=Spin.up): + def icohpvalue(self, spin: Spin = Spin.up) -> float: """ Args: spin: Spin.up or Spin.down. Returns: - float: corresponding to chosen spin. + float: ICOHP value corresponding to chosen spin. """ if not self.is_spin_polarized and spin == Spin.down: raise ValueError("The calculation was not performed with spin polarization") return self._icohp[spin] - def icohpvalue_orbital(self, orbitals, spin=Spin.up) -> float: + def icohpvalue_orbital( + self, + orbitals: tuple[Orbital, Orbital] | str, + spin: Spin = Spin.up, + ) -> float: """ Args: - orbitals: List of Orbitals or "str(Orbital1)-str(Orbital2)" - spin: Spin.up or Spin.down. + orbitals: tuple[Orbital, Orbital] or "str(Orbital0)-str(Orbital1)". + spin (Spin): Spin.up or Spin.down. Returns: - float: corresponding to chosen spin. + float: ICOHP value corresponding to chosen spin. """ if not self.is_spin_polarized and spin == Spin.down: raise ValueError("The calculation was not performed with spin polarization") - if isinstance(orbitals, list): + + if isinstance(orbitals, (tuple, list)): orbitals = f"{orbitals[0]}-{orbitals[1]}" + + assert self._orbitals is not None return self._orbitals[orbitals]["icohp"][spin] @property - def icohp(self): + def icohp(self) -> dict[Spin, float]: """Dict with ICOHPs for spin up and spin down. Returns: - dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down}. + dict[Spin, float]: {Spin.up: ICOHP_up, Spin.down: ICOHP_down}. """ return self._icohp @property - def summed_icohp(self): - """Sums ICOHPs of both spin channels for spin polarized compounds. + def summed_icohp(self) -> float: + """Summed ICOHPs of both spin channels if spin polarized. Returns: - float: icohp value in eV. + float: ICOHP value in eV. """ return self._icohp[Spin.down] + self._icohp[Spin.up] if self._is_spin_polarized else self._icohp[Spin.up] @property - def summed_orbital_icohp(self): - """Sums orbital-resolved ICOHPs of both spin channels for spin-polarized compounds. + def summed_orbital_icohp(self) -> dict[str, float]: + """Summed orbital-resolved ICOHPs of both spin channels if spin-polarized. Returns: - dict[str, float]: "str(Orbital1)-str(Ortibal2)" mapped to ICOHP value in eV. + dict[str, float]: "str(Orbital1)-str(Ortibal2)": ICOHP value in eV. """ orbital_icohp = {} + assert self._orbitals is not None for orb, item in self._orbitals.items(): orbital_icohp[orb] = ( item["icohp"][Spin.up] + item["icohp"][Spin.down] if self._is_spin_polarized else item["icohp"][Spin.up] @@ -1036,49 +1103,49 @@ def summed_orbital_icohp(self): class IcohpCollection(MSONable): - """Store IcohpValues. + """Collection of IcohpValues. Attributes: - are_coops (bool): Boolean to indicate if these are ICOOPs. - are_cobis (bool): Boolean to indicate if these are ICOOPs. - is_spin_polarized (bool): Boolean to indicate if the Lobster calculation was done spin polarized or not. + are_coops (bool): Whether these are ICOOPs. + are_cobis (bool): Whether these are ICOOPs. + is_spin_polarized (bool): Whether the calculation is spin polarized. """ def __init__( self, - list_labels, - list_atom1, - list_atom2, - list_length, - list_translation, - list_num, - list_icohp, - is_spin_polarized, - list_orb_icohp=None, - are_coops=False, - are_cobis=False, + list_labels: list[str], + list_atom1: list[str], + list_atom2: list[str], + list_length: list[float], + list_translation: list[Vector3D], + list_num: list[int], + list_icohp: list[dict[Spin, float]], + is_spin_polarized: bool, + list_orb_icohp: list[dict[str, dict[Literal["icohp", "orbitals"], Any]]] | None = None, + are_coops: bool = False, + are_cobis: bool = False, ) -> None: """ Args: - list_labels: list of labels for ICOHP/ICOOP values - list_atom1: list of str of atomnames e.g. "O1" - list_atom2: list of str of atomnames e.g. "O1" - list_length: list of lengths of corresponding bonds in Angstrom - list_translation: list of translation list, e.g. [0,0,0] - list_num: list of equivalent bonds, usually 1 starting from Lobster 3.0.0 - list_icohp: list of dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} - is_spin_polarized: Boolean to indicate if the Lobster calculation was done spin polarized or not Boolean to - indicate if the Lobster calculation was done spin polarized or not - list_orb_icohp: list of dict={[str(Orbital1)-str(Orbital2)]: {"icohp":{Spin.up: icohpvalue for spin.up, - Spin.down: icohpvalue for spin.down}, "orbitals":[Orbital1, Orbital2]}} - are_coops: Boolean to indicate whether ICOOPs are stored - are_cobis: Boolean to indicate whether ICOBIs are stored. + list_labels (list[str]): Labels for ICOHP/ICOOP values. + list_atom1 (list[str]): Atom names, e.g. "O1". + list_atom2 (list[str]): Atom names, e.g. "O1". + list_length (list[float]): Bond lengths in Angstrom. + list_translation (list[Vector3D]): Cell translation vectors. + list_num (list[int]): Numbers of equivalent bonds, usually 1 starting from LOBSTER 3.0.0. + list_icohp (list[dict]): Dicts as {Spin.up: ICOHP_up, Spin.down: ICOHP_down}. + is_spin_polarized (bool): Whether the calculation is spin polarized. + list_orb_icohp (list[dict]): Dicts as {[str(Orbital1)-str(Orbital2)]: { + "icohp": {Spin.up: IcohpValue for spin.up, Spin.down: IcohpValue for spin.down}, + "orbitals": [Orbital1, Orbital2]}. + are_coops (bool): Whether ICOOPs are stored. + are_cobis (bool): Whether ICOBIs are stored. """ if are_coops and are_cobis: raise ValueError("You cannot have info about COOPs and COBIs in the same file.") + self._are_coops = are_coops self._are_cobis = are_cobis - self._icohplist = {} self._is_spin_polarized = is_spin_polarized self._list_labels = list_labels self._list_atom1 = list_atom1 @@ -1089,132 +1156,153 @@ def __init__( self._list_icohp = list_icohp self._list_orb_icohp = list_orb_icohp - for ilist, listel in enumerate(list_labels): - self._icohplist[listel] = IcohpValue( - label=listel, - atom1=list_atom1[ilist], - atom2=list_atom2[ilist], - length=list_length[ilist], - translation=list_translation[ilist], - num=list_num[ilist], - icohp=list_icohp[ilist], + # TODO: DanielYang: self._icohplist name is misleading + # (not list), and confuses with self._list_icohp + self._icohplist: dict[str, IcohpValue] = {} + for idx, label in enumerate(list_labels): + self._icohplist[label] = IcohpValue( + label=label, + atom1=list_atom1[idx], + atom2=list_atom2[idx], + length=list_length[idx], + translation=list_translation[idx], + num=list_num[idx], + icohp=list_icohp[idx], are_coops=are_coops, are_cobis=are_cobis, - orbitals=None if list_orb_icohp is None else list_orb_icohp[ilist], + orbitals=None if list_orb_icohp is None else list_orb_icohp[idx], ) def __str__(self) -> str: - lst = [] - for value in self._icohplist.values(): - lst.append(str(value)) - return "\n".join(lst) + return "\n".join([str(value) for value in self._icohplist.values()]) - def get_icohp_by_label(self, label, summed_spin_channels=True, spin=Spin.up, orbitals=None) -> float: - """Get an icohp value for a certain bond as indicated by the label (bond labels starting by "1" as in - ICOHPLIST/ICOOPLIST). + def get_icohp_by_label( + self, + label: str, + summed_spin_channels: bool = True, + spin: Spin = Spin.up, + orbitals: str | tuple[Orbital, Orbital] | None = None, + ) -> float: + """Get an ICOHP value for a certain bond indicated by the label. Args: - label: label in str format (usually the bond number in Icohplist.lobster/Icooplist.lobster - summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed - spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned - orbitals: List of Orbital or "str(Orbital1)-str(Orbital2)" + label (str): The bond number in Icohplist.lobster/Icooplist.lobster, + starting from "1". + summed_spin_channels (bool): Whether the ICOHPs/ICOOPs of both + spin channels should be summed. + spin (Spin): If not summed_spin_channels, indicate + which spin channel should be returned. + orbitals: List of Orbital or "str(Orbital1)-str(Orbital2)". Returns: - float: ICOHP/ICOOP value + float: ICOHP/ICOOP value. """ - icohp_here: IcohpValue = self._icohplist[label] + icohp: IcohpValue = self._icohplist[label] + if orbitals is None: - if summed_spin_channels: - return icohp_here.summed_icohp - return icohp_here.icohpvalue(spin) + return icohp.summed_icohp if summed_spin_channels else icohp.icohpvalue(spin) - if isinstance(orbitals, list): + if isinstance(orbitals, (tuple, list)): orbitals = f"{orbitals[0]}-{orbitals[1]}" + if summed_spin_channels: - return icohp_here.summed_orbital_icohp[orbitals] + return icohp.summed_orbital_icohp[orbitals] - return icohp_here.icohpvalue_orbital(spin=spin, orbitals=orbitals) + return icohp.icohpvalue_orbital(spin=spin, orbitals=orbitals) - def get_summed_icohp_by_label_list(self, label_list, divisor=1.0, summed_spin_channels=True, spin=Spin.up) -> float: - """Get the sum of several ICOHP values that are indicated by a list of labels - (labels of the bonds are the same as in ICOHPLIST/ICOOPLIST). + def get_summed_icohp_by_label_list( + self, + label_list: list[str], + divisor: float = 1.0, + summed_spin_channels: bool = True, + spin: Spin = Spin.up, + ) -> float: + """Get the sum of ICOHP values. Args: - label_list: list of labels of the ICOHPs/ICOOPs that should be summed - divisor: is used to divide the sum - summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed - spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned + label_list (list[str]): Labels of the ICOHPs/ICOOPs that should be summed, + the same as in ICOHPLIST/ICOOPLIST. + divisor (float): Divisor used to divide the sum. + summed_spin_channels (bool): Whether the ICOHPs/ICOOPs of both + spin channels should be summed. + spin (Spin): If not summed_spin_channels, indicate + which spin channel should be returned. Returns: - float: sum of all ICOHPs/ICOOPs as indicated with label_list + float: Sum of ICOHPs selected with label_list. """ - sum_icohp = 0 + sum_icohp: float = 0 for label in label_list: - icohp_here = self._icohplist[label] - if icohp_here.num_bonds != 1: + icohp = self._icohplist[label] + if icohp.num_bonds != 1: warnings.warn("One of the ICOHP values is an average over bonds. This is currently not considered.") - if icohp_here._is_spin_polarized: - if summed_spin_channels: - sum_icohp = sum_icohp + icohp_here.summed_icohp - else: - sum_icohp = sum_icohp + icohp_here.icohpvalue(spin) + + if icohp._is_spin_polarized and summed_spin_channels: + sum_icohp = sum_icohp + icohp.summed_icohp else: - sum_icohp = sum_icohp + icohp_here.icohpvalue(spin) + sum_icohp = sum_icohp + icohp.icohpvalue(spin) + return sum_icohp / divisor - def get_icohp_dict_by_bondlengths(self, minbondlength=0.0, maxbondlength=8.0): - """Get a dict of IcohpValues corresponding to certain bond lengths. + def get_icohp_dict_by_bondlengths( + self, + minbondlength: float = 0.0, + maxbondlength: float = 8.0, + ) -> dict[str, IcohpValue]: + """Get IcohpValues within certain bond length range. Args: - minbondlength: defines the minimum of the bond lengths of the bonds - maxbondlength: defines the maximum of the bond lengths of the bonds. + minbondlength (float): The minimum bond length. + maxbondlength (float): The maximum bond length. Returns: - dict of IcohpValues, the keys correspond to the values from the initial list_labels. + dict[str, IcohpValue]: Keys are the labels from the initial list_labels. """ new_icohp_dict = {} for value in self._icohplist.values(): - if value._length >= minbondlength and value._length <= maxbondlength: + if minbondlength <= value._length <= maxbondlength: new_icohp_dict[value._label] = value return new_icohp_dict def get_icohp_dict_of_site( self, - site, - minsummedicohp=None, - maxsummedicohp=None, - minbondlength=0.0, - maxbondlength=8.0, - only_bonds_to=None, - ): - """Get a dict of IcohpValue for a certain site (indicated by integer). + site: int, + minsummedicohp: float | None = None, + maxsummedicohp: float | None = None, + minbondlength: float = 0.0, + maxbondlength: float = 8.0, + only_bonds_to: list[str] | None = None, + ) -> dict[str, IcohpValue]: + """Get IcohpValues for a certain site. Args: - site: integer describing the site of interest, order as in Icohplist.lobster/Icooplist.lobster, starts at 0 - minsummedicohp: float, minimal icohp/icoop of the bonds that are considered. It is the summed ICOHP value - from both spin channels for spin polarized cases - maxsummedicohp: float, maximal icohp/icoop of the bonds that are considered. It is the summed ICOHP value - from both spin channels for spin polarized cases - minbondlength: float, defines the minimum of the bond lengths of the bonds - maxbondlength: float, defines the maximum of the bond lengths of the bonds - only_bonds_to: list of strings describing the bonding partners that are allowed, e.g. ['O'] + site (int): The site of interest, ordered as in Icohplist.lobster/Icooplist.lobster, + starts from 0. + minsummedicohp (float): Minimal ICOHP/ICOOP of the bonds that are considered. + It is the summed ICOHP value from both spin channels for spin polarized cases + maxsummedicohp (float): Maximal ICOHP/ICOOP of the bonds that are considered. + It is the summed ICOHP value from both spin channels for spin polarized cases + minbondlength (float): The minimum bond length. + maxbondlength (float): The maximum bond length. + only_bonds_to (list[str]): The bonding partners that are allowed, e.g. ["O"]. Returns: - dict of IcohpValues, the keys correspond to the values from the initial list_labels + Dict of IcohpValues, the keys correspond to the values from the initial list_labels. """ new_icohp_dict = {} for key, value in self._icohplist.items(): atomnumber1 = int(re.split(r"(\d+)", value._atom1)[1]) - 1 atomnumber2 = int(re.split(r"(\d+)", value._atom2)[1]) - 1 if site in (atomnumber1, atomnumber2): - # manipulate order of atoms so that searched one is always atom1 + # Swap order of atoms so that searched one is always atom1 if site == atomnumber2: save = value._atom1 value._atom1 = value._atom2 value._atom2 = save second_test = True if only_bonds_to is None else re.split("(\\d+)", value._atom2)[0] in only_bonds_to - if value._length >= minbondlength and value._length <= maxbondlength and second_test: + if minbondlength <= value._length <= maxbondlength and second_test: + # TODO: DanielYang: merge the following condition blocks if minsummedicohp is not None: if value.summed_icohp >= minsummedicohp: if maxsummedicohp is not None: @@ -1230,16 +1318,21 @@ def get_icohp_dict_of_site( return new_icohp_dict - def extremum_icohpvalue(self, summed_spin_channels=True, spin=Spin.up): - """Get ICOHP/ICOOP of strongest bond. + def extremum_icohpvalue( + self, + summed_spin_channels: bool = True, + spin: Spin = Spin.up, + ) -> float: + """Get ICOHP/ICOOP of the strongest bond. Args: - summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed. - - spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned + summed_spin_channels (bool): Whether the ICOHPs/ICOOPs of both + spin channels should be summed. + spin (Spin): If not summed_spin_channels, this indicates which + spin channel should be returned. Returns: - lowest ICOHP/largest ICOOP value (i.e. ICOHP/ICOOP value of strongest bond) + Lowest ICOHP/largest ICOOP value (i.e. ICOHP/ICOOP value of strongest bond). """ extremum = -sys.float_info.max if self._are_coops or self._are_cobis else sys.float_info.max @@ -1255,60 +1348,71 @@ def extremum_icohpvalue(self, summed_spin_channels=True, spin=Spin.up): extremum = value.icohpvalue(spin) elif value.icohpvalue(spin) > extremum: extremum = value.icohpvalue(spin) + elif not self._are_coops and not self._are_cobis: if value.summed_icohp < extremum: extremum = value.summed_icohp + elif value.summed_icohp > extremum: extremum = value.summed_icohp + return extremum @property def is_spin_polarized(self) -> bool: - """Whether it is spin polarized.""" + """Whether this is spin polarized.""" return self._is_spin_polarized @property def are_coops(self) -> bool: - """Whether this is a coop.""" + """Whether this is COOP.""" return self._are_coops @property def are_cobis(self) -> bool: - """Whether this a cobi.""" + """Whether this is COBI.""" return self._are_cobis def get_integrated_cohp_in_energy_range( - cohp, label, orbital=None, energy_range=None, relative_E_Fermi=True, summed_spin_channels=True -): - """Integrate CompleteCohp objects which include data on integrated COHPs + cohp: CompleteCohp, + label: str, + orbital: str | None = None, + energy_range: float | tuple[float, float] | None = None, + relative_E_Fermi: bool = True, + summed_spin_channels: bool = True, +) -> float | dict[Spin, float]: + """Integrate CompleteCohps which include data of integrated COHPs (ICOHPs). + Args: - cohp: CompleteCohp object - label: label of the COHP data - orbital: If not None, a orbital resolved integrated COHP will be returned - energy_range: If None, returns icohp value at Fermi level. - If float, integrates from this float up to the Fermi level. - If [float,float], will integrate in between. - relative_E_Fermi: if True, energy scale with E_Fermi at 0 eV is chosen - summed_spin_channels: if True, Spin channels will be summed. + cohp (CompleteCohp): CompleteCohp object. + label (str): Label of the COHP data. + orbital (str): If not None, a orbital resolved integrated COHP will be returned. + energy_range: If None, return the ICOHP value at Fermi level. + If float, integrate from this value up to Fermi level. + If (float, float), integrate in between. + relative_E_Fermi (bool): Whether energy scale with Fermi level at 0 eV is chosen. + summed_spin_channels (bool): Whether Spin channels will be summed. Returns: - float indicating the integrated COHP if summed_spin_channels==True, otherwise dict of the following form { - Spin.up:float, Spin.down:float} + If summed_spin_channels: + float: the ICOHP. + else: + dict: {Spin.up: float, Spin.down: float} """ - summedicohp = {} if orbital is None: icohps = cohp.all_cohps[label].get_icohp(spin=None) - if summed_spin_channels and Spin.down in icohps: - summedicohp[Spin.up] = icohps[Spin.up] + icohps[Spin.down] - else: - summedicohp = icohps else: - icohps = cohp.get_orbital_resolved_cohp(label=label, orbitals=orbital).icohp - if summed_spin_channels and Spin.down in icohps: - summedicohp[Spin.up] = icohps[Spin.up] + icohps[Spin.down] - else: - summedicohp = icohps + _icohps = cohp.get_orbital_resolved_cohp(label=label, orbitals=orbital) + assert _icohps is not None + icohps = _icohps.icohp + + assert icohps is not None + summedicohp = {} + if summed_spin_channels and Spin.down in icohps: + summedicohp[Spin.up] = icohps[Spin.up] + icohps[Spin.down] + else: + summedicohp = icohps if energy_range is None: energies_corrected = cohp.energies - cohp.efermi @@ -1317,12 +1421,10 @@ def get_integrated_cohp_in_energy_range( if not summed_spin_channels and Spin.down in icohps: spl_spindown = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.down], ext=0) return {Spin.up: spl_spinup(0.0), Spin.down: spl_spindown(0.0)} - if summed_spin_channels: - return spl_spinup(0.0) - return {Spin.up: spl_spinup(0.0)} + return spl_spinup(0.0) if summed_spin_channels else {Spin.up: spl_spinup(0.0)} - # returns icohp value at the Fermi level! + # Return ICOHP value at the Fermi level if isinstance(energy_range, float): if relative_E_Fermi: energies_corrected = cohp.energies - cohp.efermi diff --git a/src/pymatgen/electronic_structure/core.py b/src/pymatgen/electronic_structure/core.py index c10ee208137..cfac056dd42 100644 --- a/src/pymatgen/electronic_structure/core.py +++ b/src/pymatgen/electronic_structure/core.py @@ -13,6 +13,7 @@ if TYPE_CHECKING: from collections.abc import Sequence + from numpy.typing import NDArray from typing_extensions import Self from pymatgen.core import Lattice @@ -77,7 +78,7 @@ def __str__(self) -> str: return str(self.name) @property - def orbital_type(self): + def orbital_type(self) -> OrbitalType: """OrbitalType of an orbital.""" return OrbitalType[self.name[0]] @@ -127,19 +128,18 @@ class Magmom(MSONable): """ def __init__( - self, moment: float | Sequence[float] | np.ndarray | Magmom, saxis: Sequence[float] = (0, 0, 1) + self, + moment: float | Sequence[float] | NDArray | Magmom, + saxis: Sequence[float] = (0, 0, 1), ) -> None: """ Args: moment: magnetic moment, supplied as float or list/np.ndarray saxis: spin axis, supplied as list/np.ndarray, parameter will be converted to unit vector (default is [0, 0, 1]). - - Returns: - Magmom object """ - # to init from another Magmom instance - if isinstance(moment, Magmom): + # Init from another Magmom instance + if isinstance(moment, type(self)): saxis = moment.saxis # type: ignore[has-type] moment = moment.moment # type: ignore[has-type] @@ -153,6 +153,63 @@ def __init__( self.saxis = saxis / np.linalg.norm(saxis) + def __getitem__(self, key): + return self.moment[key] + + def __iter__(self): + return iter(self.moment) + + def __abs__(self) -> float: + return np.linalg.norm(self.moment) + + def __eq__(self, other: object) -> bool: + """Equal if 'global' magnetic moments are the same, saxis can differ.""" + try: + other_magmom = type(self)(other) + except (TypeError, ValueError): + return NotImplemented + + return np.allclose(self.global_moment, other_magmom.global_moment) + + def __lt__(self, other: Self) -> bool: + return abs(self) < abs(other) + + def __neg__(self) -> Self: + return type(self)(-self.moment, saxis=self.saxis) + + def __hash__(self) -> int: + return hash(tuple(self.moment) + tuple(self.saxis)) + + def __float__(self) -> float: + """Get magnitude of magnetic moment with a sign with respect to + an arbitrary direction. + + Should give unsurprising output if Magmom is treated like a + scalar or if a set of Magmoms describes a collinear structure. + + Implemented this way rather than simpler abs(self) so that + moments will have a consistent sign in case of e.g. + antiferromagnetic collinear structures without additional + user intervention. + + However, should be used with caution for non-collinear + structures and might give nonsensical results except in the case + of only slightly non-collinear structures (e.g. small canting). + + This approach is also used to obtain "diff" VolumetricDensity + in pymatgen.io.vasp.outputs.VolumetricDensity when processing + Chgcars from SOC calculations. + """ + return float(self.get_00t_magmom_with_xyz_saxis()[2]) + + def __str__(self) -> str: + return str(float(self)) + + def __repr__(self) -> str: + if np.allclose(self.saxis, (0, 0, 1)): + return f"Magnetic moment {self.moment}" + return f"Magnetic moment {self.moment} (spin axis = {self.saxis})" + @classmethod def from_global_moment_and_saxis(cls, global_moment, saxis) -> Self: """Convenience method to initialize Magmom from a given global @@ -166,7 +223,7 @@ def from_global_moment_and_saxis(cls, global_moment, saxis) -> Self: global_moment: global magnetic moment saxis: desired saxis """ - magmom = Magmom(global_moment) + magmom = cls(global_moment) return cls(magmom.get_moment(saxis=saxis), saxis=saxis) @classmethod @@ -217,26 +274,26 @@ def get_moment(self, saxis=(0, 0, 1)): Returns: np.ndarray of length 3 """ - # transform back to moment with spin axis [0, 0, 1] + # Transform back to moment with spin axis [0, 0, 1] trafo_mat_inv = self._get_transformation_matrix_inv(self.saxis) moment = np.matmul(self.moment, trafo_mat_inv) - # transform to new saxis + # Transform to new saxis trafo_mat = self._get_transformation_matrix(saxis) moment = np.matmul(moment, trafo_mat) - # round small values to zero + # Round small values to zero moment[np.abs(moment) < 1e-8] = 0 return moment @property - def global_moment(self) -> np.ndarray: + def global_moment(self) -> NDArray: """The magnetic moment defined in an arbitrary global reference frame as an np.array of length 3.""" return self.get_moment() @property - def projection(self): + def projection(self) -> float: """Projects moment along spin quantization axis. Useful for obtaining collinear approximation for slightly non-collinear magmoms. @@ -245,16 +302,16 @@ def projection(self): """ return np.dot(self.moment, self.saxis) - def get_xyz_magmom_with_001_saxis(self): + def get_xyz_magmom_with_001_saxis(self) -> Self: """Get a Magmom in the default setting of saxis = [0, 0, 1] and the magnetic moment rotated as required. Returns: Magmom """ - return Magmom(self.get_moment()) + return type(self)(self.get_moment()) - def get_00t_magmom_with_xyz_saxis(self): + def get_00t_magmom_with_xyz_saxis(self) -> Self: """For internal implementation reasons, in non-collinear calculations VASP prefers the following. MAGMOM = 0 0 total_magnetic_moment @@ -276,7 +333,7 @@ def get_00t_magmom_with_xyz_saxis(self): Returns: Magmom """ - # reference direction gives sign of moment + # Reference direction gives sign of moment # entirely arbitrary, there will always be a pathological case # where a consistent sign is not possible if the magnetic moments # are aligned along the reference direction, but in practice this @@ -288,8 +345,8 @@ def get_00t_magmom_with_xyz_saxis(self): if np.dot(ref_direction, new_saxis) < 0: t = -t new_saxis = -new_saxis - return Magmom([0, 0, t], saxis=new_saxis) - return Magmom(self) + return type(self)([0, 0, t], saxis=new_saxis) + return type(self)(self) @staticmethod def have_consistent_saxis(magmoms) -> bool: @@ -339,14 +396,14 @@ def get_suggested_saxis(magmoms): Returns: np.ndarray of length 3 """ - # heuristic, will pick largest magmom as reference + # Heuristic, will pick largest magmom as reference # useful for creating collinear approximations of # e.g. slightly canted magnetic structures # for fully collinear structures, will return expected # result magmoms = [Magmom(magmom) for magmom in magmoms] - # filter only non-zero magmoms + # Filter only non-zero magmoms magmoms = [magmom for magmom in magmoms if abs(magmom)] magmoms.sort(reverse=True) if len(magmoms) > 0: @@ -367,20 +424,24 @@ def are_collinear(magmoms) -> bool: if not Magmom.have_consistent_saxis(magmoms): magmoms = Magmom.get_consistent_set_and_saxis(magmoms)[0] - # convert to numpy array for convenience + # Convert to numpy array for convenience magmoms = np.array([list(magmom) for magmom in magmoms]) magmoms = magmoms[np.any(magmoms, axis=1)] # remove zero magmoms if len(magmoms) == 0: return True - # use first moment as reference to compare against + # Use first moment as reference to compare against ref_magmom = magmoms[0] - # magnitude of cross products != 0 if non-collinear with reference + # Magnitude of cross products != 0 if non-collinear with reference num_ncl = np.count_nonzero(np.linalg.norm(np.cross(ref_magmom, magmoms), axis=1)) return num_ncl == 0 @classmethod - def from_moment_relative_to_crystal_axes(cls, moment: list[float], lattice: Lattice) -> Self: + def from_moment_relative_to_crystal_axes( + cls, + moment: list[float], + lattice: Lattice, + ) -> Self: """Obtaining a Magmom object from a magnetic moment provided relative to crystal axes. @@ -393,14 +454,14 @@ def from_moment_relative_to_crystal_axes(cls, moment: list[float], lattice: Latt Returns: Magmom """ - # get matrix representing unit lattice vectors + # Get matrix representing unit lattice vectors unit_m = lattice.matrix / np.linalg.norm(lattice.matrix, axis=1)[:, None] moment = np.matmul(list(moment), unit_m) - # round small values to zero + # Round small values to zero moment[np.abs(moment) < 1e-8] = 0 return cls(moment) - def get_moment_relative_to_crystal_axes(self, lattice): + def get_moment_relative_to_crystal_axes(self, lattice: Lattice): """If scalar magmoms, moments will be given arbitrarily along z. Used for writing moments to magCIF file. @@ -410,67 +471,9 @@ def get_moment_relative_to_crystal_axes(self, lattice): Returns: vector as list of floats """ - # get matrix representing unit lattice vectors + # Get matrix representing unit lattice vectors unit_m = lattice.matrix / np.linalg.norm(lattice.matrix, axis=1)[:, None] - # note np.matmul() requires numpy version >= 1.10 moment = np.matmul(self.global_moment, np.linalg.inv(unit_m)) - # round small values to zero + # Round small values to zero moment[np.abs(moment) < 1e-8] = 0 return moment - - def __getitem__(self, key): - return self.moment[key] - - def __iter__(self): - return iter(self.moment) - - def __abs__(self): - return np.linalg.norm(self.moment) - - def __eq__(self, other: object) -> bool: - """Equal if 'global' magnetic moments are the same, saxis can differ.""" - try: - other_magmom = Magmom(other) - except (TypeError, ValueError): - return NotImplemented - - return np.allclose(self.global_moment, other_magmom.global_moment) - - def __lt__(self, other): - return abs(self) < abs(other) - - def __neg__(self): - return Magmom(-self.moment, saxis=self.saxis) - - def __hash__(self) -> int: - return hash(tuple(self.moment) + tuple(self.saxis)) - - def __float__(self) -> float: - """Get magnitude of magnetic moment with a sign with respect to - an arbitrary direction. - - Should give unsurprising output if Magmom is treated like a - scalar or if a set of Magmoms describes a collinear structure. - - Implemented this way rather than simpler abs(self) so that - moments will have a consistent sign in case of e.g. - antiferromagnetic collinear structures without additional - user intervention. - - However, should be used with caution for non-collinear - structures and might give nonsensical results except in the case - of only slightly non-collinear structures (e.g. small canting). - - This approach is also used to obtain "diff" VolumetricDensity - in pymatgen.io.vasp.outputs.VolumetricDensity when processing - Chgcars from SOC calculations. - """ - return float(self.get_00t_magmom_with_xyz_saxis()[2]) - - def __str__(self) -> str: - return str(float(self)) - - def __repr__(self) -> str: - if np.allclose(self.saxis, (0, 0, 1)): - return f"Magnetic moment {self.moment}" - return f"Magnetic moment {self.moment} (spin axis = {self.saxis})" diff --git a/src/pymatgen/electronic_structure/dos.py b/src/pymatgen/electronic_structure/dos.py index 43605e08ade..3f1ed3feebe 100644 --- a/src/pymatgen/electronic_structure/dos.py +++ b/src/pymatgen/electronic_structure/dos.py @@ -52,6 +52,18 @@ def __init__(self, energies: ArrayLike, densities: ArrayLike, efermi: float) -> super().__init__(energies, densities, efermi) self.efermi = efermi + def __str__(self) -> str: + """Get a string which can be easily plotted (using gnuplot).""" + if Spin.down in self.densities: + str_arr = [f"#{'Energy':30s} {'DensityUp':30s} {'DensityDown':30s}"] + for idx, energy in enumerate(self.energies): + str_arr.append(f"{energy:.5f} {self.densities[Spin.up][idx]:.5f} {self.densities[Spin.down][idx]:.5f}") + else: + str_arr = [f"#{'Energy':30s} {'DensityUp':30s}"] + for idx, energy in enumerate(self.energies): + str_arr.append(f"{energy:.5f} {self.densities[Spin.up][idx]:.5f}") + return "\n".join(str_arr) + def get_interpolated_gap(self, tol: float = 0.001, abs_tol: bool = False, spin: Spin | None = None): """Expects a DOS object and finds the gap. @@ -148,18 +160,6 @@ def get_gap(self, tol: float = 0.001, abs_tol: bool = False, spin: Spin | None = cbm, vbm = self.get_cbm_vbm(tol, abs_tol, spin) return max(cbm - vbm, 0.0) - def __str__(self) -> str: - """Get a string which can be easily plotted (using gnuplot).""" - if Spin.down in self.densities: - str_arr = [f"#{'Energy':30s} {'DensityUp':30s} {'DensityDown':30s}"] - for i, energy in enumerate(self.energies): - str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f} {self.densities[Spin.down][i]:.5f}") - else: - str_arr = [f"#{'Energy':30s} {'DensityUp':30s}"] - for i, energy in enumerate(self.energies): - str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f}") - return "\n".join(str_arr) - class Dos(MSONable): """Basic DOS object. All other DOS objects are extended versions of this @@ -189,6 +189,33 @@ def __init__( vol = norm_vol or 1 self.densities = {k: np.array(d) / vol for k, d in densities.items()} + def __add__(self, other): + """Add two DOS together. Checks that energy scales are the same. + Otherwise, a ValueError is thrown. + + Args: + other: Another DOS object. + + Returns: + Sum of the two DOSs. + """ + if not all(np.equal(self.energies, other.energies)): + raise ValueError("Energies of both DOS are not compatible!") + densities = {spin: self.densities[spin] + other.densities[spin] for spin in self.densities} + return Dos(self.efermi, self.energies, densities) + + def __str__(self) -> str: + """Get a string which can be easily plotted (using gnuplot).""" + if Spin.down in self.densities: + str_arr = [f"#{'Energy':30s} {'DensityUp':30s} {'DensityDown':30s}"] + for i, energy in enumerate(self.energies): + str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f} {self.densities[Spin.down][i]:.5f}") + else: + str_arr = [f"#{'Energy':30s} {'DensityUp':30s}"] + for i, energy in enumerate(self.energies): + str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f}") + return "\n".join(str_arr) + def get_densities(self, spin: Spin | None = None): """Get the density of states for a particular spin. @@ -200,8 +227,9 @@ def get_densities(self, spin: Spin | None = None): None, the sum of all spins is returned. """ if self.densities is None: - result = None - elif spin is None: + return None + + if spin is None: if Spin.down in self.densities: result = self.densities[Spin.up] + self.densities[Spin.down] else: @@ -227,21 +255,6 @@ def get_smeared_densities(self, sigma: float): smeared_dens[spin] = gaussian_filter1d(dens, sigma / avg_diff) return smeared_dens - def __add__(self, other): - """Add two DOS together. Checks that energy scales are the same. - Otherwise, a ValueError is thrown. - - Args: - other: Another DOS object. - - Returns: - Sum of the two DOSs. - """ - if not all(np.equal(self.energies, other.energies)): - raise ValueError("Energies of both DOS are not compatible!") - densities = {spin: self.densities[spin] + other.densities[spin] for spin in self.densities} - return Dos(self.efermi, self.energies, densities) - def get_interpolated_value(self, energy: float) -> dict[Spin, float]: """Get interpolated density for a particular energy. @@ -341,18 +354,6 @@ def get_gap(self, tol: float = 0.001, abs_tol: bool = False, spin: Spin | None = cbm, vbm = self.get_cbm_vbm(tol, abs_tol, spin) return max(cbm - vbm, 0.0) - def __str__(self) -> str: - """Get a string which can be easily plotted (using gnuplot).""" - if Spin.down in self.densities: - str_arr = [f"#{'Energy':30s} {'DensityUp':30s} {'DensityDown':30s}"] - for i, energy in enumerate(self.energies): - str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f} {self.densities[Spin.down][i]:.5f}") - else: - str_arr = [f"#{'Energy':30s} {'DensityUp':30s}"] - for i, energy in enumerate(self.energies): - str_arr.append(f"{energy:.5f} {self.densities[Spin.up][i]:.5f}") - return "\n".join(str_arr) - @classmethod def from_dict(cls, dct: dict) -> Self: """Get Dos object from dict representation of Dos.""" @@ -626,6 +627,9 @@ def __init__( self.pdos = pdoss self.structure = structure + def __str__(self) -> str: + return f"Complete DOS for {self.structure}" + def get_normalized(self) -> CompleteDos: """Get a normalized version of the CompleteDos.""" if self.norm_vol is not None: @@ -775,7 +779,7 @@ def spin_polarization(self) -> float | None: n_F_down = n_F[Spin.down] if (n_F_up + n_F_down) == 0: - # only well defined for metals or half-metals + # Only well defined for metals or half-metals return float("NaN") spin_polarization = (n_F_up - n_F_down) / (n_F_up + n_F_down) @@ -1285,9 +1289,6 @@ def as_dict(self) -> dict: dct["spd_dos"] = {str(orb): dos.as_dict() for orb, dos in self.get_spd_dos().items()} return dct - def __str__(self) -> str: - return f"Complete DOS for {self.structure}" - class LobsterCompleteDos(CompleteDos): """Extended CompleteDOS for Lobster.""" @@ -1452,8 +1453,24 @@ def _get_orb_type_lobster(orb) -> OrbitalType | None: Returns: OrbitalType """ - orb_labs = ["s", "p_y", "p_z", "p_x", "d_xy", "d_yz", "d_z^2", "d_xz", "d_x^2-y^2"] - orb_labs += ["f_y(3x^2-y^2)", "f_xyz", "f_yz^2", "f_z^3", "f_xz^2", "f_z(x^2-y^2)", "f_x(x^2-3y^2)"] + orb_labs = ( + "s", + "p_y", + "p_z", + "p_x", + "d_xy", + "d_yz", + "d_z^2", + "d_xz", + "d_x^2-y^2", + "f_y(3x^2-y^2)", + "f_xyz", + "f_yz^2", + "f_z^3", + "f_xz^2", + "f_z(x^2-y^2)", + "f_x(x^2-3y^2)", + ) try: orbital = Orbital(orb_labs.index(orb[1:])) @@ -1471,7 +1488,7 @@ def _get_orb_lobster(orb): Returns: Orbital. """ - orb_labs = [ + orb_labs = ( "s", "p_y", "p_z", @@ -1488,10 +1505,10 @@ def _get_orb_lobster(orb): "f_xz^2", "f_z(x^2-y^2)", "f_x(x^2-3y^2)", - ] + ) try: return Orbital(orb_labs.index(orb[1:])) except AttributeError: print("Orb not in list") - return None + return None diff --git a/src/pymatgen/electronic_structure/plotter.py b/src/pymatgen/electronic_structure/plotter.py index f4cfa52be78..d7737e80a32 100644 --- a/src/pymatgen/electronic_structure/plotter.py +++ b/src/pymatgen/electronic_structure/plotter.py @@ -51,22 +51,27 @@ class DosPlotter: - """Plot DOSs. The interface is extremely flexible given there are many + """Plot DOS. The interface is extremely flexible given there are many different ways in which people want to view DOS. Typical usage is: - # Initializes plotter with some optional args. Defaults are usually fine + # Initialize plotter with some optional args. Defaults are usually fine plotter = PhononDosPlotter(). # Add DOS with a label plotter.add_dos("Total DOS", dos) - # Alternatively, you can add a dict of DOSes. This is the typical form + # Alternatively, you can add a dict of DOS. This is the typical form # returned by CompletePhononDos.get_element_dos(). plotter.add_dos_dict({"dos1": dos1, "dos2": dos2}) plotter.add_dos_dict(complete_dos.get_spd_dos()) """ - def __init__(self, zero_at_efermi: bool = True, stack: bool = False, sigma: float | None = None) -> None: + def __init__( + self, + zero_at_efermi: bool = True, + stack: bool = False, + sigma: float | None = None, + ) -> None: """ Args: zero_at_efermi (bool): Whether to shift all Dos to have zero energy at the diff --git a/src/pymatgen/io/lobster/__init__.py b/src/pymatgen/io/lobster/__init__.py index f5e4a5f4234..e4b1100a2d6 100644 --- a/src/pymatgen/io/lobster/__init__.py +++ b/src/pymatgen/io/lobster/__init__.py @@ -1,5 +1,5 @@ """ -This package implements modules for input and output to and from Lobster. It +This package implements modules for input and output to and from LOBSTER. It imports the key classes form both lobster.inputs and lobster_outputs to allow most classes to be simply called as pymatgen.io.lobster.Lobsterin for example, to retain backwards compatibility. diff --git a/src/pymatgen/io/lobster/outputs.py b/src/pymatgen/io/lobster/outputs.py index aacc38332b4..31ebe16b04d 100644 --- a/src/pymatgen/io/lobster/outputs.py +++ b/src/pymatgen/io/lobster/outputs.py @@ -279,6 +279,7 @@ def _get_bond_data(line: str, are_multi_center_cobis: bool = False) -> dict: indices, a tuple containing the orbitals (if orbital-resolved), and a label for the orbitals (if orbital-resolved). """ + if not are_multi_center_cobis: line_new = line.rsplit("(", 1) length = float(line_new[-1][:-1]) @@ -357,7 +358,7 @@ def __init__( filename: Name of the ICOHPLIST file. If it is None, the default file name will be chosen, depending on the value of are_coops is_spin_polarized: Boolean to indicate if the calculation is spin polarized - icohpcollection: IcohpCollection Object. + icohpcollection: IcohpCollection Object """ self._filename = filename @@ -384,14 +385,14 @@ def __init__( if len(data) == 0: raise RuntimeError("ICOHPLIST file contains no data.") - # Which Lobster version? + # Determine LOBSTER version if len(data[0].split()) == 8: version = "3.1.1" elif len(data[0].split()) == 6: version = "2.2.1" - warnings.warn("Please consider using the new Lobster version. See www.cohp.de.") + warnings.warn("Please consider using the new LOBSTER version. See www.cohp.de.") else: - raise ValueError + raise ValueError("Unsupported LOBSTER version.") # If the calculation is spin polarized, the line in the middle # of the file will be another header line. @@ -408,9 +409,9 @@ def __init__( data_orbitals = [] for line in data: if "_" not in line.split()[1]: - data_without_orbitals += [line] + data_without_orbitals.append(line) else: - data_orbitals += [line] + data_orbitals.append(line) else: data_without_orbitals = data @@ -423,49 +424,45 @@ def __init__( else: n_bonds = len(data_without_orbitals) - labels, atoms1, atoms2, lens, translations, nums, icohps = [], [], [], [], [], [], [] - - # initialize static variables - label = "" - atom1 = "" - atom2 = "" - length = None - num = None - translation = [] + labels: list[str] = [] + atoms1: list[str] = [] + atoms2: list[str] = [] + lens: list[float] = [] + translations: list[tuple[int, int, int]] = [] + nums: list[int] = [] + icohps: list[dict[Spin, float]] = [] for bond in range(n_bonds): - line = data_without_orbitals[bond].split() - icohp = {} + line_parts = data_without_orbitals[bond].split() + + label = f"{line_parts[0]}" + atom1 = str(line_parts[1]) + atom2 = str(line_parts[2]) + length = float(line_parts[3]) + + icohp: dict[Spin, float] = {} if version == "2.2.1": - label = f"{line[0]}" - atom1 = str(line[1]) - atom2 = str(line[2]) - length = float(line[3]) - icohp[Spin.up] = float(line[4]) - num = int(line[5]) - translation = [0, 0, 0] + icohp[Spin.up] = float(line_parts[4]) + num = int(line_parts[5]) + translation = (0, 0, 0) if self.is_spin_polarized: icohp[Spin.down] = float(data_without_orbitals[bond + n_bonds + 1].split()[4]) - elif version == "3.1.1": - label = f"{line[0]}" - atom1 = str(line[1]) - atom2 = str(line[2]) - length = float(line[3]) - translation = [int(line[4]), int(line[5]), int(line[6])] - icohp[Spin.up] = float(line[7]) + else: # version == "3.1.1" + translation = (int(line_parts[4]), int(line_parts[5]), int(line_parts[6])) + icohp[Spin.up] = float(line_parts[7]) num = 1 if self.is_spin_polarized: icohp[Spin.down] = float(data_without_orbitals[bond + n_bonds + 1].split()[7]) - labels += [label] - atoms1 += [atom1] - atoms2 += [atom2] - lens += [length] - translations += [translation] - nums += [num] - icohps += [icohp] + labels.append(label) + atoms1.append(atom1) + atoms2.append(atom2) + lens.append(length) + translations.append(translation) + nums.append(num) + icohps.append(icohp) list_orb_icohp: list[dict] | None = None if self.orbitalwise: @@ -475,17 +472,17 @@ def __init__( for i_data_orb in range(n_orbs): data_orb = data_orbitals[i_data_orb] icohp = {} - line = data_orb.split() - label = f"{line[0]}" + line_parts = data_orb.split() + label = f"{line_parts[0]}" orbs = re.findall(r"_(.*?)(?=\s)", data_orb) orb_label, orbitals = get_orb_from_str(orbs) - icohp[Spin.up] = float(line[7]) + icohp[Spin.up] = float(line_parts[7]) if self.is_spin_polarized: icohp[Spin.down] = float(data_orbitals[n_orbs + i_data_orb].split()[7]) if len(list_orb_icohp) < int(label): - list_orb_icohp += [{orb_label: {"icohp": icohp, "orbitals": orbitals}}] + list_orb_icohp.append({orb_label: {"icohp": icohp, "orbitals": orbitals}}) else: list_orb_icohp[int(label) - 1][orb_label] = {"icohp": icohp, "orbitals": orbitals} @@ -499,7 +496,7 @@ def __init__( list_atom1=atoms1, list_atom2=atoms2, list_length=lens, - list_translation=translations, + list_translation=translations, # type: ignore[arg-type] list_num=nums, list_icohp=icohps, is_spin_polarized=self.is_spin_polarized, @@ -542,11 +539,12 @@ def __init__( filename: PathLike | None = "NcICOBILIST.lobster", ) -> None: """ - LOBSTER < 4.1.0: no COBI/ICOBI/NcICOBI. + LOBSTER < 4.1.0: no COBI/ICOBI/NcICOBI Args: filename: Name of the NcICOBILIST file. """ + # LOBSTER list files have an extra trailing blank line # and we don't need the header. with zopen(filename, mode="rt") as file: @@ -574,7 +572,7 @@ def __init__( data_without_orbitals = [] for line in data: if "_" not in str(line.split()[3:]) and "s]" not in str(line.split()[3:]): - data_without_orbitals += [line] + data_without_orbitals.append(line) else: data_without_orbitals = data @@ -605,17 +603,19 @@ def __init__( if self.is_spin_polarized: ncicobi[Spin.down] = float(data_without_orbitals[bond + n_bonds + 1].split()[2]) - self.list_labels += [label] - self.list_n_atoms += [n_atoms] - self.list_ncicobi += [ncicobi] - self.list_interaction_type += [interaction_type] - self.list_num += [num] + self.list_labels.append(label) + self.list_n_atoms.append(n_atoms) + self.list_ncicobi.append(ncicobi) + self.list_interaction_type.append(interaction_type) + self.list_num.append(num) # TODO: add functions to get orbital resolved NcICOBIs @property def ncicobi_list(self) -> dict[Any, dict[str, Any]]: - """Returns: ncicobilist.""" + """ + Returns: ncicobilist. + """ ncicobi_list = {} for idx in range(len(self.list_labels)): ncicobi_list[str(idx + 1)] = { @@ -688,7 +688,7 @@ def _parse_doscar(self): for nd in range(1, ndos): line = file.readline().split() cdos[nd] = np.array(line) - dos += [cdos] + dos.append(cdos) doshere = np.array(dos[0]) if len(doshere[0, :]) == 5: self._is_spin_polarized = True @@ -710,7 +710,7 @@ def _parse_doscar(self): for orb_num, j in enumerate(range(1, ncol)): orb = orbitals[atom + 1][orb_num] pdos[orb][spin] = data[:, j] - pdoss += [pdos] + pdoss.append(pdos) else: tdensities[Spin.up] = doshere[:, 1] tdensities[Spin.down] = doshere[:, 2] @@ -728,7 +728,7 @@ def _parse_doscar(self): pdos[orb][spin] = data[:, j] if j % 2 == 0: orb_num += 1 - pdoss += [pdos] + pdoss.append(pdos) self._efermi = efermi self._pdos = pdoss @@ -744,32 +744,32 @@ def _parse_doscar(self): @property def completedos(self) -> LobsterCompleteDos: - """LobsterCompleteDos.""" + """LobsterCompleteDos""" return self._completedos @property def pdos(self) -> list: - """Projected DOS.""" + """Projected DOS""" return self._pdos @property def tdos(self) -> Dos: - """Total DOS.""" + """Total DOS""" return self._tdos @property def energies(self) -> np.ndarray: - """Energies.""" + """Energies""" return self._energies @property def tdensities(self) -> dict[Spin, np.ndarray]: - """Total densities as a np.ndarray.""" + """total densities as a np.ndarray""" return self._tdensities @property def itdensities(self) -> dict[Spin, np.ndarray]: - """Integrated total densities as a np.ndarray.""" + """integrated total densities as a np.ndarray""" return self._itdensities @property @@ -805,7 +805,7 @@ def __init__( atomlist: list of atoms in the structure types: list of unique species in the structure mulliken: list of Mulliken charges - loewdin: list of Loewdin charges. + loewdin: list of Loewdin charges """ self._filename = filename self.num_atoms = num_atoms @@ -821,15 +821,15 @@ def __init__( raise RuntimeError("CHARGE file contains no data.") self.num_atoms = len(data) - for atom in range(self.num_atoms): - line = data[atom].split() - self.atomlist += [line[1] + line[0]] - self.types += [line[1]] - self.mulliken += [float(line[2])] - self.loewdin += [float(line[3])] + for atom_idx in range(self.num_atoms): + line_parts = data[atom_idx].split() + self.atomlist.append(line_parts[1] + line_parts[0]) + self.types.append(line_parts[1]) + self.mulliken.append(float(line_parts[2])) + self.loewdin.append(float(line_parts[3])) def get_structure_with_charges(self, structure_filename: PathLike) -> Structure: - """Get a Structure with Mulliken and Loewdin charges as site properties. + """Get a Structure with Mulliken and Loewdin charges as site properties Args: structure_filename: filename of POSCAR @@ -926,7 +926,7 @@ def __init__(self, filename: PathLike | None, **kwargs) -> None: """ Args: filename: The lobsterout file. - **kwargs: dict to initialize Lobsterout instance. + **kwargs: dict to initialize Lobsterout instance """ self.filename = filename if kwargs: @@ -1008,7 +1008,7 @@ def __init__(self, filename: PathLike | None, **kwargs) -> None: raise ValueError("must provide either filename or kwargs to initialize Lobsterout") def get_doc(self) -> dict[str, Any]: - """Get the LobsterDict with all the information stored in lobsterout.""" + """Get a dict with all the information in lobsterout.""" return { # Check if LOBSTER starts from a projection "restart_from_projection": self.is_restart_from_projection, @@ -1039,7 +1039,7 @@ def get_doc(self) -> dict[str, Any]: } def as_dict(self) -> dict: - """MSONable dict.""" + """MSONable dict""" dct = dict(vars(self)) dct["@module"] = type(self).__module__ dct["@class"] = type(self).__name__ @@ -1090,9 +1090,9 @@ def _get_spillings(data, number_of_spins): splitrow = row.split() if len(splitrow) > 2 and splitrow[2] == "spilling:": if splitrow[1] == "charge": - charge_spilling += [np.float64(splitrow[3].replace("%", "")) / 100.0] + charge_spilling.append(np.float64(splitrow[3].replace("%", "")) / 100.0) if splitrow[1] == "total": - total_spilling += [np.float64(splitrow[3].replace("%", "")) / 100.0] + total_spilling.append(np.float64(splitrow[3].replace("%", "")) / 100.0) if len(charge_spilling) == number_of_spins and len(total_spilling) == number_of_spins: break @@ -1108,8 +1108,8 @@ def _get_elements_basistype_basisfunctions(data): basisfunctions = [] for row in data: if begin and not end: - splitrow = row.split() - if splitrow[0] not in [ + row_parts = row.split() + if row_parts[0] not in { "INFO:", "WARNING:", "setting", @@ -1118,11 +1118,11 @@ def _get_elements_basistype_basisfunctions(data): "saving", "spillings", "writing", - ]: - elements += [splitrow[0]] - basistype += [splitrow[1].replace("(", "").replace(")", "")] + }: + elements.append(row_parts[0]) + basistype.append(row_parts[1].replace("(", "").replace(")", "")) # last sign is a '' - basisfunctions += [splitrow[2:]] + basisfunctions += [row_parts[2:]] else: end = True if "setting up local basis functions..." in row: @@ -1159,7 +1159,7 @@ def _get_warning_orthonormalization(data): for row in data: splitrow = row.split() if "orthonormalized" in splitrow: - orthowarning += [" ".join(splitrow[1:])] + orthowarning.append(" ".join(splitrow[1:])) return orthowarning @staticmethod @@ -1168,7 +1168,7 @@ def _get_all_warning_lines(data): for row in data: splitrow = row.split() if len(splitrow) > 0 and splitrow[0] == "WARNING:": - ws += [" ".join(splitrow[1:])] + ws.append(" ".join(splitrow[1:])) return ws @staticmethod @@ -1177,7 +1177,7 @@ def _get_all_info_lines(data): for row in data: splitrow = row.split() if len(splitrow) > 0 and splitrow[0] == "INFO:": - infos += [" ".join(splitrow[1:])] + infos.append(" ".join(splitrow[1:])) return infos @@ -1218,10 +1218,10 @@ def __init__( "FATBAND_*" files will be read kpoints_file (PathLike): KPOINTS file for bandstructure calculation, typically "KPOINTS". vasprun_file (PathLike): Corresponding vasprun file. - Instead, the Fermi energy from the DFT run can be provided. Then, + Instead, the Fermi level from the DFT run can be provided. Then, this value should be set to None. structure (Structure): Structure object. - efermi (float): fermi energy in eV. + efermi (float): Fermi level in eV. """ warnings.warn("Make sure all relevant FATBAND files were generated and read in!") warnings.warn("Use Lobster 3.2.0 or newer for fatband calculations!") @@ -1260,7 +1260,7 @@ def __init__( filenames = "." for name in os.listdir(filenames): if fnmatch.fnmatch(name, "FATBAND_*.lobster"): - filenames_new += [os.path.join(filenames, name)] + filenames_new.append(os.path.join(filenames, name)) filenames = filenames_new if len(filenames) == 0: raise ValueError("No FATBAND files in folder or given") @@ -1268,17 +1268,17 @@ def __init__( with zopen(name, mode="rt") as file: contents = file.read().split("\n") - atom_names += [os.path.split(name)[1].split("_")[1].capitalize()] + atom_names.append(os.path.split(name)[1].split("_")[1].capitalize()) parameters = contents[0].split() - atom_type += [re.split(r"[0-9]+", parameters[3])[0].capitalize()] - orbital_names += [parameters[4]] + atom_type.append(re.split(r"[0-9]+", parameters[3])[0].capitalize()) + orbital_names.append(parameters[4]) # get atomtype orbital dict atom_orbital_dict = {} # type: dict - for iatom, atom in enumerate(atom_names): + for idx, atom in enumerate(atom_names): if atom not in atom_orbital_dict: atom_orbital_dict[atom] = [] - atom_orbital_dict[atom] += [orbital_names[iatom]] + atom_orbital_dict[atom].append(orbital_names[idx]) # test if there are the same orbitals twice or if two different formats were used or if all necessary orbitals # are there for items in atom_orbital_dict.values(): @@ -1286,9 +1286,9 @@ def __init__( raise ValueError("The are two FATBAND files for the same atom and orbital. The program will stop.") split = [] for item in items: - split += [item.split("_")[0]] + split.append(item.split("_")[0]) for number in collections.Counter(split).values(): - if number not in (1, 3, 5, 7): + if number not in {1, 3, 5, 7}: raise ValueError( "Make sure all relevant orbitals were generated and that no duplicates (2p and 2p_x) are " "present" @@ -1313,7 +1313,7 @@ def __init__( linenumbers = [] for iline, line in enumerate(contents[1 : self.nbands * 2 + 4]): if line.split()[0] == "#": - linenumbers += [iline] + linenumbers.append(iline) if ifilename == 0: self.is_spinpolarized = len(linenumbers) == 2 @@ -1363,7 +1363,7 @@ def __init__( ] ) if ifilename == 0: - kpoints_array += [KPOINT] + kpoints_array.append(KPOINT) linenumber = 0 iband = 0 @@ -1406,7 +1406,7 @@ def get_bandstructure(self) -> LobsterBandStructureSymmLine: kpoints=self.kpoints_array, eigenvals=self.eigenvals, lattice=self.lattice, - efermi=self.efermi, + efermi=self.efermi, # type: ignore[arg-type] labels_dict=self.label_dict, structure=self.structure, projections=self.p_eigenvals, @@ -1415,7 +1415,6 @@ def get_bandstructure(self) -> LobsterBandStructureSymmLine: class Bandoverlaps(MSONable): """Read in bandOverlaps.lobster files. These files are not created during every Lobster run. - Attributes: band_overlaps_dict (dict[Spin, Dict[str, Dict[str, Union[float, np.ndarray]]]]): A dictionary containing the band overlap data of the form: {spin: {"kpoint as string": {"maxDeviation": @@ -1455,7 +1454,7 @@ def __init__( def _read(self, contents: list, spin_numbers: list): """ - Will read in all contents of the file. + Will read in all contents of the file Args: contents: list of strings @@ -1475,7 +1474,7 @@ def _read(self, contents: list, spin_numbers: list): kpoint_array = [] for kpointel in kpoint: if kpointel not in {"at", "k-point", ""}: - kpoint_array += [float(kpointel)] + kpoint_array.append(float(kpointel)) elif "maxDeviation" in line: if spin not in self.band_overlaps_dict: @@ -1487,23 +1486,23 @@ def _read(self, contents: list, spin_numbers: list): if "matrices" not in self.band_overlaps_dict[spin]: self.band_overlaps_dict[spin]["matrices"] = [] maxdev = line.split(" ")[2] - self.band_overlaps_dict[spin]["max_deviations"] += [float(maxdev)] + self.band_overlaps_dict[spin]["max_deviations"].append(float(maxdev)) self.band_overlaps_dict[spin]["k_points"] += [kpoint_array] - self.max_deviation += [float(maxdev)] + self.max_deviation.append(float(maxdev)) overlaps = [] else: rows = [] for el in line.split(" "): if el != "": - rows += [float(el)] + rows.append(float(el)) overlaps += [rows] if len(overlaps) == len(rows): self.band_overlaps_dict[spin]["matrices"] += [np.matrix(overlaps)] def has_good_quality_maxDeviation(self, limit_maxDeviation: float = 0.1) -> bool: """ - Will check if the maxDeviation from the ideal bandoverlap is smaller or equal to limit_maxDeviation. + Will check if the maxDeviation from the ideal bandoverlap is smaller or equal to limit_maxDeviation Args: limit_maxDeviation: limit of the maxDeviation @@ -1584,7 +1583,7 @@ def __init__(self, filename: str = "GROSSPOP.lobster", list_dict_grosspop: list[ """ Args: filename: filename of the "GROSSPOP.lobster" file - list_dict_grosspop: List of dictionaries including all information about the gross populations. + list_dict_grosspop: List of dictionaries including all information about the gross populations """ # opens file self._filename = filename @@ -1607,10 +1606,10 @@ def __init__(self, filename: str = "GROSSPOP.lobster", list_dict_grosspop: list[ small_dict["Mulliken GP"][cleanline[0]] = float(cleanline[1]) small_dict["Loewdin GP"][cleanline[0]] = float(cleanline[2]) if "total" in cleanline[0]: - self.list_dict_grosspop += [small_dict] + self.list_dict_grosspop.append(small_dict) def get_structure_with_total_grosspop(self, structure_filename: str) -> Structure: - """Get a Structure with Mulliken and Loewdin total grosspopulations as site properties. + """Get a Structure with Mulliken and Loewdin total grosspopulations as site properties Args: structure_filename (str): filename of POSCAR @@ -1619,7 +1618,6 @@ def get_structure_with_total_grosspop(self, structure_filename: str) -> Structur Structure Object with Mulliken and Loewdin total grosspopulations as site properties. """ struct = Structure.from_file(structure_filename) - # site_properties: dict[str, Any] = {} mullikengp = [] loewdingp = [] for grosspop in self.list_dict_grosspop: @@ -1668,9 +1666,9 @@ def _parse_file(filename): splitline = line.split() if len(splitline) >= 6: points += [[float(splitline[0]), float(splitline[1]), float(splitline[2])]] - distance += [float(splitline[3])] - real += [float(splitline[4])] - imaginary += [float(splitline[5])] + distance.append(float(splitline[3])) + real.append(float(splitline[4])) + imaginary.append(float(splitline[5])) if len(real) != grid[0] * grid[1] * grid[2] or len(imaginary) != grid[0] * grid[1] * grid[2]: raise ValueError("Something went wrong while reading the file") @@ -1716,9 +1714,9 @@ def set_volumetric_data(self, grid, structure): "coordinates 0.0 0.0 0.0 coordinates 1.0 1.0 1.0 box bandlist 1 " ) - new_x += [x_here] - new_y += [y_here] - new_z += [z_here] + new_x.append(x_here) + new_y.append(y_here) + new_z.append(z_here) new_real += [self.real[runner]] new_imaginary += [self.imaginary[runner]] @@ -1885,7 +1883,7 @@ def __init__( sitepotentials_loewdin: Loewdin site potential sitepotentials_mulliken: Mulliken site potential madelungenergies_loewdin: Madelung energy based on the Loewdin approach - madelungenergies_mulliken: Madelung energy based on the Mulliken approach. + madelungenergies_mulliken: Madelung energy based on the Mulliken approach """ self._filename = filename self.ewald_splitting = [] if ewald_splitting is None else ewald_splitting @@ -1910,17 +1908,17 @@ def __init__( data = data[5:-1] self.num_atoms = len(data) - 2 for atom in range(self.num_atoms): - line = data[atom].split() - self.atomlist += [line[1] + str(line[0])] - self.types += [line[1]] - self.sitepotentials_mulliken += [float(line[2])] - self.sitepotentials_loewdin += [float(line[3])] + line_parts = data[atom].split() + self.atomlist.append(line_parts[1] + str(line_parts[0])) + self.types.append(line_parts[1]) + self.sitepotentials_mulliken.append(float(line_parts[2])) + self.sitepotentials_loewdin.append(float(line_parts[3])) self.madelungenergies_mulliken = float(data[self.num_atoms + 1].split()[3]) self.madelungenergies_loewdin = float(data[self.num_atoms + 1].split()[4]) def get_structure_with_site_potentials(self, structure_filename): - """Get a Structure with Mulliken and Loewdin charges as site properties. + """Get a Structure with Mulliken and Loewdin charges as site properties Args: structure_filename: filename of POSCAR @@ -1983,7 +1981,7 @@ def get_orb_from_str(orbs): list of tw Orbital objects """ # TODO: also useful for plotting of DOS - orb_labs = [ + orb_labs = ( "s", "p_y", "p_z", @@ -2000,7 +1998,7 @@ def get_orb_from_str(orbs): "f_xz^2", "f_z(x^2-y^2)", "f_x(x^2-3y^2)", - ] + ) orbitals = [(int(orb[0]), Orbital(orb_labs.index(orb[1:]))) for orb in orbs] orb_label = "" @@ -2055,8 +2053,9 @@ def __init__(self, e_fermi=None, filename: str = "hamiltonMatrices.lobster"): Args: filename: filename for the hamiltonMatrices file, typically "hamiltonMatrices.lobster". e_fermi: fermi level in eV for the structure only - relevant if input file contains hamilton matrices data. + relevant if input file contains hamilton matrices data """ + self._filename = filename # hamiltonMatrices with zopen(self._filename, mode="rt") as file: @@ -2101,21 +2100,21 @@ def _parse_matrix(file_data, pattern, e_fermi): for idx, line in enumerate(file_data): line = line.strip() if "Real parts" in line: - start_inxs_real += [idx + 1] + start_inxs_real.append(idx + 1) if idx == 1: # ignore the first occurrence as files start with real matrices pass else: - end_inxs_imag += [idx - 1] + end_inxs_imag.append(idx - 1) matches = re.search(pattern, file_data[idx - 1]) if matches and len(matches.groups()) == 2: k_point = matches.group(2) complex_matrices[k_point] = {} if "Imag parts" in line: - end_inxs_real += [idx - 1] - start_inxs_imag += [idx + 1] + end_inxs_real.append(idx - 1) + start_inxs_imag.append(idx + 1) # explicitly add the last line as files end with imaginary matrix if idx == len(file_data) - 1: - end_inxs_imag += [len(file_data)] + end_inxs_imag.append(len(file_data)) # extract matrix data and store diagonal elements matrix_real = [] @@ -2136,13 +2135,13 @@ def _parse_matrix(file_data, pattern, e_fermi): matches = re.search(pattern, file_data[start_inx_real - 2]) if matches and len(matches.groups()) == 2: - spin = Spin.up if matches.group(1) == "1" else Spin.down - k_point = matches.group(2) + spin = Spin.up if matches[1] == "1" else Spin.down + k_point = matches[2] complex_matrices[k_point].update({spin: comp_matrix}) elif matches and len(matches.groups()) == 1: - k_point = matches.group(1) + k_point = matches[1] complex_matrices.update({k_point: comp_matrix}) - matrix_diagonal_values += [comp_matrix.real.diagonal() - e_fermi] + matrix_diagonal_values.append(comp_matrix.real.diagonal() - e_fermi) # extract elements basis functions as list elements_basis_functions = [ diff --git a/src/pymatgen/io/vasp/outputs.py b/src/pymatgen/io/vasp/outputs.py index e5da42e9c0b..5ec3838e64a 100644 --- a/src/pymatgen/io/vasp/outputs.py +++ b/src/pymatgen/io/vasp/outputs.py @@ -1033,7 +1033,7 @@ def get_band_structure( eigenvals, lattice_new, e_fermi, - labels_dict, + labels_dict, # type: ignore[arg-type] structure=self.final_structure, projections=p_eig_vals, ) diff --git a/src/pymatgen/util/typing.py b/src/pymatgen/util/typing.py index 7f498c5e93f..52246c384df 100644 --- a/src/pymatgen/util/typing.py +++ b/src/pymatgen/util/typing.py @@ -7,9 +7,10 @@ from collections.abc import Sequence from os import PathLike as OsPathLike -from typing import TYPE_CHECKING, Any, Union +from typing import TYPE_CHECKING, Any, Literal, Union from pymatgen.core import Composition, DummySpecies, Element, Species +from pymatgen.electronic_structure.core import Spin if TYPE_CHECKING: # needed to avoid circular imports from pymatgen.analysis.cost import CostEntry # type: ignore[attr-defined] @@ -25,6 +26,9 @@ PathLike = Union[str, OsPathLike] PbcLike = tuple[bool, bool, bool] +# Things that can be cast to a Spin +SpinLike = Union[Spin, Literal[-1, 1, "up", "down"]] + # Things that can be cast to a Species-like object using get_el_sp SpeciesLike = Union[str, Element, Species, DummySpecies] diff --git a/tests/io/lobster/test_inputs.py b/tests/io/lobster/test_inputs.py index 646a2f70eb3..6675eaf70c8 100644 --- a/tests/io/lobster/test_inputs.py +++ b/tests/io/lobster/test_inputs.py @@ -444,77 +444,77 @@ def test_values(self): "length": 2.88231, "number_of_bonds": 3, "icohp": {Spin.up: -2.18042}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "2": { "length": 3.10144, "number_of_bonds": 3, "icohp": {Spin.up: -1.14347}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "3": { "length": 2.88231, "number_of_bonds": 3, "icohp": {Spin.up: -2.18042}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "4": { "length": 3.10144, "number_of_bonds": 3, "icohp": {Spin.up: -1.14348}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "5": { "length": 3.05001, "number_of_bonds": 3, "icohp": {Spin.up: -1.30006}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "6": { "length": 2.91676, "number_of_bonds": 3, "icohp": {Spin.up: -1.96843}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "7": { "length": 3.05001, "number_of_bonds": 3, "icohp": {Spin.up: -1.30006}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "8": { "length": 2.91676, "number_of_bonds": 3, "icohp": {Spin.up: -1.96843}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "9": { "length": 3.37522, "number_of_bonds": 3, "icohp": {Spin.up: -0.47531}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "10": { "length": 3.07294, "number_of_bonds": 3, "icohp": {Spin.up: -2.38796}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "11": { "length": 3.37522, "number_of_bonds": 3, "icohp": {Spin.up: -0.47531}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, } @@ -523,77 +523,77 @@ def test_values(self): "length": 2.88231, "number_of_bonds": 3, "icohp": {Spin.up: 0.14245}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "2": { "length": 3.10144, "number_of_bonds": 3, "icohp": {Spin.up: -0.04118}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "3": { "length": 2.88231, "number_of_bonds": 3, "icohp": {Spin.up: 0.14245}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "4": { "length": 3.10144, "number_of_bonds": 3, "icohp": {Spin.up: -0.04118}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "5": { "length": 3.05001, "number_of_bonds": 3, "icohp": {Spin.up: -0.03516}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "6": { "length": 2.91676, "number_of_bonds": 3, "icohp": {Spin.up: 0.10745}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "7": { "length": 3.05001, "number_of_bonds": 3, "icohp": {Spin.up: -0.03516}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "8": { "length": 2.91676, "number_of_bonds": 3, "icohp": {Spin.up: 0.10745}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "9": { "length": 3.37522, "number_of_bonds": 3, "icohp": {Spin.up: -0.12395}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "10": { "length": 3.07294, "number_of_bonds": 3, "icohp": {Spin.up: 0.24714}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "11": { "length": 3.37522, "number_of_bonds": 3, "icohp": {Spin.up: -0.12395}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, } @@ -602,14 +602,14 @@ def test_values(self): "length": 2.83189, "number_of_bonds": 2, "icohp": {Spin.up: -0.10218, Spin.down: -0.19701}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, "2": { "length": 2.45249, "number_of_bonds": 1, "icohp": {Spin.up: -0.28485, Spin.down: -0.58279}, - "translation": [0, 0, 0], + "translation": (0, 0, 0), "orbitals": None, }, }