diff --git a/pymatgen/io/lobster/outputs.py b/pymatgen/io/lobster/outputs.py index cd4425b4a93..807be33d568 100644 --- a/pymatgen/io/lobster/outputs.py +++ b/pymatgen/io/lobster/outputs.py @@ -20,6 +20,7 @@ import numpy as np from monty.io import zopen +from monty.json import MSONable from pymatgen.core.structure import Structure from pymatgen.electronic_structure.bandstructure import LobsterBandStructureSymmLine @@ -225,7 +226,7 @@ def _get_bond_data(line: str) -> dict: } -class Icohplist: +class Icohplist(MSONable): """ Class to read ICOHPLIST/ICOOPLIST files generated by LOBSTER. @@ -241,7 +242,15 @@ class Icohplist: IcohpCollection (IcohpCollection): IcohpCollection Object. """ - def __init__(self, are_coops: bool = False, are_cobis: bool = False, filename: str | None = None): + def __init__( + self, + are_coops: bool = False, + are_cobis: bool = False, + filename: str | None = None, + is_spin_polarized: bool = False, + orbitalwise: bool = False, + icohpcollection=None, + ): """ Args: are_coops: Determines if the file is a list of ICOOPs. @@ -249,8 +258,15 @@ def __init__(self, are_coops: bool = False, are_cobis: bool = False, filename: s are_cobis: Determines if the file is a list of ICOBIs. Defaults to False for ICOHPs. filename: Name of the ICOHPLIST file. If it is None, the default - file name will be chosen, depending on the value of are_coops. + 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 + """ + self._filename = filename + self.is_spin_polarized = is_spin_polarized + self.orbitalwise = orbitalwise + self._icohpcollection = icohpcollection 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 @@ -265,129 +281,130 @@ def __init__(self, are_coops: bool = False, are_cobis: bool = False, filename: s # LOBSTER list files have an extra trailing blank line # and we don't need the header. - with zopen(filename, mode="rt") as file: - data = file.read().split("\n")[1:-1] - if len(data) == 0: - raise OSError("ICOHPLIST file contains no data.") - - # Which 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.") - else: - raise ValueError - - # If the calculation is spin polarized, the line in the middle - # of the file will be another header line. - # TODO: adapt this for orbital-wise stuff - self.is_spin_polarized = "distance" in data[len(data) // 2] - - # check if orbital-wise ICOHPLIST - # include case when there is only one ICOHP!!! - self.orbitalwise = len(data) > 2 and "_" in data[1].split()[1] - - if self.orbitalwise: - data_without_orbitals = [] - data_orbitals = [] - for line in data: - if "_" not in line.split()[1]: - data_without_orbitals.append(line) - else: - data_orbitals.append(line) + if self._icohpcollection is None: + with zopen(filename, mode="rt") as file: + data = file.read().split("\n")[1:-1] + if len(data) == 0: + raise OSError("ICOHPLIST file contains no data.") - else: - data_without_orbitals = data + # Which 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.") + else: + raise ValueError - if "distance" in data_without_orbitals[len(data_without_orbitals) // 2]: + # If the calculation is spin polarized, the line in the middle + # of the file will be another header line. # TODO: adapt this for orbital-wise stuff - num_bonds = len(data_without_orbitals) // 2 - if num_bonds == 0: - raise OSError("ICOHPLIST file contains no data.") - else: - num_bonds = len(data_without_orbitals) + self.is_spin_polarized = "distance" in data[len(data) // 2] - list_labels = [] - list_atom1 = [] - list_atom2 = [] - list_length = [] - list_translation = [] - list_num = [] - list_icohp = [] + # check if orbital-wise ICOHPLIST + # include case when there is only one ICOHP!!! + self.orbitalwise = len(data) > 2 and "_" in data[1].split()[1] - for bond in range(num_bonds): - line = data_without_orbitals[bond].split() - icohp = {} - 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] - if self.is_spin_polarized: - icohp[Spin.down] = float(data_without_orbitals[bond + num_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]) - num = 1 - - if self.is_spin_polarized: - icohp[Spin.down] = float(data_without_orbitals[bond + num_bonds + 1].split()[7]) - - list_labels.append(label) - list_atom1.append(atom1) - list_atom2.append(atom2) - list_length.append(length) - list_translation.append(translation) - list_num.append(num) - list_icohp.append(icohp) - - list_orb_icohp: list[dict] | None = None - if self.orbitalwise: - list_orb_icohp = [] - num_orbs = len(data_orbitals) // 2 if self.is_spin_polarized else len(data_orbitals) - - for i_data_orb in range(num_orbs): - data_orb = data_orbitals[i_data_orb] - icohp = {} - line = data_orb.split() - label = f"{line[0]}" - orbs = re.findall(r"_(.*?)(?=\s)", data_orb) - orb_label, orbitals = get_orb_from_str(orbs) - icohp[Spin.up] = float(line[7]) + if self.orbitalwise: + data_without_orbitals = [] + data_orbitals = [] + for line in data: + if "_" not in line.split()[1]: + data_without_orbitals.append(line) + else: + data_orbitals.append(line) - if self.is_spin_polarized: - icohp[Spin.down] = float(data_orbitals[num_orbs + i_data_orb].split()[7]) + else: + data_without_orbitals = data - if len(list_orb_icohp) < int(label): - list_orb_icohp.append({orb_label: {"icohp": icohp, "orbitals": orbitals}}) - else: - list_orb_icohp[int(label) - 1][orb_label] = {"icohp": icohp, "orbitals": orbitals} - - # to avoid circular dependencies - from pymatgen.electronic_structure.cohp import IcohpCollection - - self._icohpcollection = IcohpCollection( - are_coops=are_coops, - are_cobis=are_cobis, - list_labels=list_labels, - list_atom1=list_atom1, - list_atom2=list_atom2, - list_length=list_length, - list_translation=list_translation, - list_num=list_num, - list_icohp=list_icohp, - is_spin_polarized=self.is_spin_polarized, - list_orb_icohp=list_orb_icohp, - ) + if "distance" in data_without_orbitals[len(data_without_orbitals) // 2]: + # TODO: adapt this for orbital-wise stuff + num_bonds = len(data_without_orbitals) // 2 + if num_bonds == 0: + raise OSError("ICOHPLIST file contains no data.") + else: + num_bonds = len(data_without_orbitals) + + list_labels = [] + list_atom1 = [] + list_atom2 = [] + list_length = [] + list_translation = [] + list_num = [] + list_icohp = [] + + for bond in range(num_bonds): + line = data_without_orbitals[bond].split() + icohp = {} + 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] + if self.is_spin_polarized: + icohp[Spin.down] = float(data_without_orbitals[bond + num_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]) + num = 1 + + if self.is_spin_polarized: + icohp[Spin.down] = float(data_without_orbitals[bond + num_bonds + 1].split()[7]) + + list_labels.append(label) + list_atom1.append(atom1) + list_atom2.append(atom2) + list_length.append(length) + list_translation.append(translation) + list_num.append(num) + list_icohp.append(icohp) + + list_orb_icohp: list[dict] | None = None + if self.orbitalwise: + list_orb_icohp = [] + num_orbs = len(data_orbitals) // 2 if self.is_spin_polarized else len(data_orbitals) + + for i_data_orb in range(num_orbs): + data_orb = data_orbitals[i_data_orb] + icohp = {} + line = data_orb.split() + label = f"{line[0]}" + orbs = re.findall(r"_(.*?)(?=\s)", data_orb) + orb_label, orbitals = get_orb_from_str(orbs) + icohp[Spin.up] = float(line[7]) + + if self.is_spin_polarized: + icohp[Spin.down] = float(data_orbitals[num_orbs + i_data_orb].split()[7]) + + if len(list_orb_icohp) < int(label): + list_orb_icohp.append({orb_label: {"icohp": icohp, "orbitals": orbitals}}) + else: + list_orb_icohp[int(label) - 1][orb_label] = {"icohp": icohp, "orbitals": orbitals} + + # to avoid circular dependencies + from pymatgen.electronic_structure.cohp import IcohpCollection + + self._icohpcollection = IcohpCollection( + are_coops=are_coops, + are_cobis=are_cobis, + list_labels=list_labels, + list_atom1=list_atom1, + list_atom2=list_atom2, + list_length=list_length, + list_translation=list_translation, + list_num=list_num, + list_icohp=list_icohp, + is_spin_polarized=self.is_spin_polarized, + list_orb_icohp=list_orb_icohp, + ) @property def icohplist(self) -> dict[Any, dict[str, Any]]: @@ -662,39 +679,56 @@ def is_spin_polarized(self) -> bool: return self._is_spin_polarized -class Charge: +class Charge(MSONable): """ Class to read CHARGE files generated by LOBSTER. Attributes: atomlist (list[str]): List of atoms in CHARGE.lobster. types (list[str]): List of types of atoms in CHARGE.lobster. - Mulliken (list[float]): List of Mulliken charges of atoms in CHARGE.lobster. - Loewdin (list[float]): List of Loewdin charges of atoms in CHARGE.Loewdin. + mulliken (list[float]): List of Mulliken charges of atoms in CHARGE.lobster. + loewdin (list[float]): List of Loewdin charges of atoms in CHARGE.Loewdin. num_atoms (int): Number of atoms in CHARGE.lobster. """ - def __init__(self, filename: str = "CHARGE.lobster"): + def __init__( + self, + filename: str = "CHARGE.lobster", + num_atoms: int | None = None, + atomlist: list[str] | None = None, + types: list[str] | None = None, + mulliken: list[float] | None = None, + loewdin: list[float] | None = None, + ): """ Args: filename: filename for the CHARGE file, typically "CHARGE.lobster". + num_atoms: number of atoms in the structure + 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 """ - with zopen(filename, mode="rt") as file: - data = file.read().split("\n")[3:-3] - if len(data) == 0: - raise OSError("CHARGES file contains no data.") - - self.num_atoms = len(data) - self.atomlist: list[str] = [] - self.types: list[str] = [] - self.Mulliken: list[float] = [] - self.Loewdin: list[float] = [] - for atom in range(self.num_atoms): - line = data[atom].split() - self.atomlist.append(line[1] + line[0]) - self.types.append(line[1]) - self.Mulliken.append(float(line[2])) - self.Loewdin.append(float(line[3])) + self._filename = filename + self.num_atoms = num_atoms + self.types = [] if types is None else types + self.atomlist = [] if atomlist is None else atomlist + self.mulliken = [] if mulliken is None else mulliken + self.loewdin = [] if loewdin is None else loewdin + + if self.num_atoms is None: + with zopen(filename, mode="rt") as file: + data = file.read().split("\n")[3:-3] + if len(data) == 0: + raise OSError("CHARGES file contains no data.") + + self.num_atoms = len(data) + for atom in range(self.num_atoms): + line = data[atom].split() + self.atomlist.append(line[1] + line[0]) + self.types.append(line[1]) + self.mulliken.append(float(line[2])) + self.loewdin.append(float(line[3])) def get_structure_with_charges(self, structure_filename): """ @@ -707,11 +741,21 @@ def get_structure_with_charges(self, structure_filename): Structure Object with Mulliken and Loewdin charges as site properties. """ struct = Structure.from_file(structure_filename) - Mulliken = self.Mulliken - Loewdin = self.Loewdin - site_properties = {"Mulliken Charges": Mulliken, "Loewdin Charges": Loewdin} + mulliken = self.mulliken + loewdin = self.loewdin + site_properties = {"Mulliken Charges": mulliken, "Loewdin Charges": loewdin} return struct.copy(site_properties=site_properties) + @property + def Mulliken(self): + warnings.warn("`Mulliken` attribute is deprecated. Use `mulliken` instead.", DeprecationWarning, stacklevel=2) + return self.mulliken + + @property + def Loewdin(self): + warnings.warn("`Loewdin` attribute is deprecated. Use `loewdin` instead.", DeprecationWarning, stacklevel=2) + return self.loewdin + class Lobsterout: """ @@ -1212,28 +1256,43 @@ def get_bandstructure(self): ) -class Bandoverlaps: +class Bandoverlaps(MSONable): """ Class to read in bandOverlaps.lobster files. These files are not created during every Lobster run. Attributes: - bandoverlapsdict (dict[Spin, Dict[str, Dict[str, Union[float, np.ndarray]]]]): A dictionary + 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": float that describes the max deviation, "matrix": 2D array of the size number of bands times number of bands including the overlap matrices with}}}. - maxDeviation (list[float]): A list of floats describing the maximal deviation for each problematic kpoint. + max_deviation (list[float]): A list of floats describing the maximal deviation for each problematic kpoint. """ - def __init__(self, filename: str = "bandOverlaps.lobster"): + def __init__( + self, + filename: str = "bandOverlaps.lobster", + band_overlaps_dict: dict[Any, dict] | None = None, # Any is spin number 1 or -1 + max_deviation: list[float] | None = None, + ): """ Args: filename: filename of the "bandOverlaps.lobster" file. + band_overlaps_dict: A dictionary containing the band overlap data of the form: {spin: {"kpoint as string": + {"max_deviation":float that describes the max deviation, "matrix": 2D array of the size number of bands + times number of bands including the overlap matrices with}}}. + max_deviation (list[float]): A list of floats describing the maximal deviation for each problematic kpoint. """ - with zopen(filename, mode="rt") as file: - contents = file.read().split("\n") + self._filename = filename + self.band_overlaps_dict = {} if band_overlaps_dict is None else band_overlaps_dict + self.max_deviation = [] if max_deviation is None else max_deviation - spin_numbers = [0, 1] if contents[0].split()[-1] == "0" else [1, 2] + if not self.band_overlaps_dict: + with zopen(filename, mode="rt") as file: + contents = file.read().split("\n") + + spin_numbers = [0, 1] if contents[0].split()[-1] == "0" else [1, 2] - self._read(contents, spin_numbers) + self._filename = filename + self._read(contents, spin_numbers) def _read(self, contents: list, spin_numbers: list): """ @@ -1243,8 +1302,6 @@ def _read(self, contents: list, spin_numbers: list): contents: list of strings spin_numbers: list of spin numbers depending on `Lobster` version. """ - self.bandoverlapsdict: dict[Any, dict] = {} # Any is spin number 1 or -1 - self.max_deviation = [] # This has to be done like this because there can be different numbers of problematic k-points per spin for line in contents: if f"Overlap Matrix (abs) of the orthonormalized projected bands for spin {spin_numbers[0]}" in line: @@ -1259,21 +1316,21 @@ def _read(self, contents: list, spin_numbers: list): kpoint_array.append(str(kpointel)) elif "maxDeviation" in line: - if spin not in self.bandoverlapsdict: - self.bandoverlapsdict[spin] = {} - if " ".join(kpoint_array) not in self.bandoverlapsdict[spin]: - self.bandoverlapsdict[spin][" ".join(kpoint_array)] = {} + if spin not in self.band_overlaps_dict: + self.band_overlaps_dict[spin] = {} + if " ".join(kpoint_array) not in self.band_overlaps_dict[spin]: + self.band_overlaps_dict[spin][" ".join(kpoint_array)] = {} maxdev = line.split(" ")[2] - self.bandoverlapsdict[spin][" ".join(kpoint_array)]["maxDeviation"] = float(maxdev) + self.band_overlaps_dict[spin][" ".join(kpoint_array)]["maxDeviation"] = float(maxdev) self.max_deviation.append(float(maxdev)) - self.bandoverlapsdict[spin][" ".join(kpoint_array)]["matrix"] = [] + self.band_overlaps_dict[spin][" ".join(kpoint_array)]["matrix"] = [] else: overlaps = [] for el in line.split(" "): if el not in [""]: overlaps.append(float(el)) - self.bandoverlapsdict[spin][" ".join(kpoint_array)]["matrix"].append(overlaps) + self.band_overlaps_dict[spin][" ".join(kpoint_array)]["matrix"].append(overlaps) def has_good_quality_maxDeviation(self, limit_maxDeviation: float = 0.1) -> bool: """ @@ -1307,7 +1364,7 @@ def has_good_quality_check_occupied_bands( Returns: Boolean that will give you information about the quality of the projection """ - for matrix in self.bandoverlapsdict[Spin.up].values(): + for matrix in self.band_overlaps_dict[Spin.up].values(): for iband1, band1 in enumerate(matrix["matrix"]): for iband2, band2 in enumerate(band1): if iband1 < number_occ_bands_spin_up and iband2 < number_occ_bands_spin_up: @@ -1318,7 +1375,7 @@ def has_good_quality_check_occupied_bands( return False if spin_polarized: - for matrix in self.bandoverlapsdict[Spin.down].values(): + for matrix in self.band_overlaps_dict[Spin.down].values(): for iband1, band1 in enumerate(matrix["matrix"]): for iband2, band2 in enumerate(band1): if number_occ_bands_spin_down is not None: @@ -1332,8 +1389,17 @@ def has_good_quality_check_occupied_bands( ValueError("number_occ_bands_spin_down has to be specified") return True + @property + def bandoverlapsdict(self): + warnings.warn( + "`bandoverlapsdict` attribute is deprecated. Use `band_overlaps_dict` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.band_overlaps_dict + -class Grosspop: +class Grosspop(MSONable): """ Class to read in GROSSPOP.lobster files. @@ -1348,31 +1414,34 @@ class Grosspop: The 0th entry of the list refers to the first atom in GROSSPOP.lobster and so on. """ - def __init__(self, filename: str = "GROSSPOP.lobster"): + def __init__(self, filename: str = "GROSSPOP.lobster", list_dict_grosspop: list[dict] | None = None): """ Args: - filename: filename of the "GROSSPOP.lobster" file. + filename: filename of the "GROSSPOP.lobster" file + list_dict_grosspop: List of dictionaries including all information about the gross populations """ # opens file - with zopen(filename, mode="rt") as file: - contents = file.read().split("\n") + self._filename = filename + self.list_dict_grosspop = [] if list_dict_grosspop is None else list_dict_grosspop - self.list_dict_grosspop = [] - # transfers content of file to list of dict - for line in contents[3:]: - cleanline = [i for i in line.split(" ") if i != ""] - if len(cleanline) == 5: - small_dict = {} - small_dict["element"] = cleanline[1] - small_dict["Mulliken GP"] = {} - small_dict["Loewdin GP"] = {} - small_dict["Mulliken GP"][cleanline[2]] = float(cleanline[3]) - small_dict["Loewdin GP"][cleanline[2]] = float(cleanline[4]) - elif len(cleanline) > 0: - 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.append(small_dict) + if not self.list_dict_grosspop: + with zopen(filename, mode="rt") as file: + contents = file.read().split("\n") + # transfers content of file to list of dict + for line in contents[3:]: + cleanline = [i for i in line.split(" ") if i != ""] + if len(cleanline) == 5: + small_dict = {} + small_dict["element"] = cleanline[1] + small_dict["Mulliken GP"] = {} + small_dict["Loewdin GP"] = {} + small_dict["Mulliken GP"][cleanline[2]] = float(cleanline[3]) + small_dict["Loewdin GP"][cleanline[2]] = float(cleanline[4]) + elif len(cleanline) > 0: + 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.append(small_dict) def get_structure_with_total_grosspop(self, structure_filename: str) -> Structure: """ @@ -1385,7 +1454,7 @@ 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] = {} + # site_properties: dict[str, Any] = {} mullikengp = [] loewdingp = [] for grosspop in self.list_dict_grosspop: @@ -1567,32 +1636,63 @@ def write_file(self, filename="WAVECAR.vasp", part="real"): # madleung and sitepotential classes -class MadelungEnergies: +class MadelungEnergies(MSONable): """ Class to read MadelungEnergies.lobster files generated by LOBSTER. Attributes: - madelungenergies_Mulliken (float): Float that gives the Madelung energy based on the Mulliken approach. - madelungenergies_Loewdin (float): Float that gives the Madelung energy based on the Loewdin approach. + madelungenergies_mulliken (float): Float that gives the Madelung energy based on the Mulliken approach. + madelungenergies_loewdin (float): Float that gives the Madelung energy based on the Loewdin approach. ewald_splitting (float): Ewald splitting parameter to compute SitePotentials. """ - def __init__(self, filename: str = "MadelungEnergies.lobster"): + def __init__( + self, + filename: str = "MadelungEnergies.lobster", + ewald_splitting: float | None = None, + madelungenergies_mulliken: float | None = None, + madelungenergies_loewdin: float | None = None, + ): """ Args: filename: filename of the "MadelungEnergies.lobster" file. """ - with zopen(filename, mode="rt") as file: - data = file.read().split("\n")[5] - if len(data) == 0: - raise OSError("MadelungEnergies file contains no data.") - line = data.split() - self.ewald_splitting = float(line[0]) - self.madelungenergies_Mulliken = float(line[1]) - self.madelungenergies_Loewdin = float(line[2]) + self._filename = filename + self.ewald_splitting = None if ewald_splitting is None else ewald_splitting + self.madelungenergies_loewdin = None if madelungenergies_loewdin is None else madelungenergies_loewdin + self.madelungenergies_mulliken = None if madelungenergies_mulliken is None else madelungenergies_mulliken + + if self.ewald_splitting is None: + with zopen(filename, mode="rt") as file: + data = file.read().split("\n")[5] + if len(data) == 0: + raise OSError("MadelungEnergies file contains no data.") + line = data.split() + self._filename = filename + self.ewald_splitting = float(line[0]) + self.madelungenergies_mulliken = float(line[1]) + self.madelungenergies_loewdin = float(line[2]) + + @property + def madelungenergies_Loewdin(self): + warnings.warn( + "`madelungenergies_Loewdin` attribute is deprecated. Use `madelungenergies_loewdin` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.madelungenergies_loewdin + @property + def madelungenergies_Mulliken(self): + warnings.warn( + "`madelungenergies_Mulliken` attribute is deprecated. Use `madelungenergies_mulliken` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.madelungenergies_mulliken -class SitePotential: + +class SitePotential(MSONable): """ Class to read SitePotentials.lobster files generated by LOBSTER. @@ -1600,41 +1700,68 @@ class SitePotential: atomlist (list[str]): List of atoms in SitePotentials.lobster. types (list[str]): List of types of atoms in SitePotentials.lobster. num_atoms (int): Number of atoms in SitePotentials.lobster. - sitepotentials_Mulliken (list[float]): List of Mulliken potentials of sites in SitePotentials.lobster. - sitepotentials_Loewdin (list[float]): List of Loewdin potentials of sites in SitePotentials.lobster. - madelung_Mulliken (float): Float that gives the Madelung energy based on the Mulliken approach. - madelung_Loewdin (float): Float that gives the Madelung energy based on the Loewdin approach. + sitepotentials_mulliken (list[float]): List of Mulliken potentials of sites in SitePotentials.lobster. + sitepotentials_loewdin (list[float]): List of Loewdin potentials of sites in SitePotentials.lobster. + madelungenergies_mulliken (float): Float that gives the Madelung energy based on the Mulliken approach. + madelungenergies_loewdin (float): Float that gives the Madelung energy based on the Loewdin approach. ewald_splitting (float): Ewald Splitting parameter to compute SitePotentials. """ - def __init__(self, filename: str = "SitePotentials.lobster"): + def __init__( + self, + filename: str = "SitePotentials.lobster", + ewald_splitting: float | None = None, + num_atoms: int | None = None, + atomlist: list[str] | None = None, + types: list[str] | None = None, + sitepotentials_loewdin: list[float] | None = None, + sitepotentials_mulliken: list[float] | None = None, + madelungenergies_mulliken: float | None = None, + madelungenergies_loewdin: float | None = None, + ): """ Args: - filename: filename for the SitePotentials file, typically "SitePotentials.lobster". + filename: filename for the SitePotentials file, typically "SitePotentials.lobster" + ewald_splitting: ewald splitting parameter used for computing madelung energies + num_atoms: number of atoms in the structure + atomlist: list of atoms in the structure + types: list of unique atom types in the structure + 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 """ - # site_potentials - with zopen(filename, mode="rt") as file: - data = file.read().split("\n") - if len(data) == 0: - raise OSError("SitePotentials file contains no data.") - - self.ewald_splitting = float(data[0].split()[9]) - - data = data[5:-1] - self.num_atoms = len(data) - 2 - self.atomlist: list[str] = [] - self.types: list[str] = [] - self.sitepotentials_Mulliken: list[float] = [] - self.sitepotentials_Loewdin: list[float] = [] - for atom in range(self.num_atoms): - line = data[atom].split() - self.atomlist.append(line[1] + str(line[0])) - self.types.append(line[1]) - self.sitepotentials_Mulliken.append(float(line[2])) - self.sitepotentials_Loewdin.append(float(line[3])) - - self.madelungenergies_Mulliken = float(data[self.num_atoms + 1].split()[3]) - self.madelungenergies_Loewdin = float(data[self.num_atoms + 1].split()[4]) + self._filename = filename + self.ewald_splitting = [] if ewald_splitting is None else ewald_splitting + self.num_atoms = num_atoms + self.types = [] if types is None else types + self.atomlist = [] if atomlist is None else atomlist + self.sitepotentials_loewdin = [] if sitepotentials_loewdin is None else sitepotentials_loewdin + self.sitepotentials_mulliken = [] if sitepotentials_mulliken is None else sitepotentials_mulliken + self.madelungenergies_loewdin = [] if madelungenergies_loewdin is None else madelungenergies_loewdin + self.madelungenergies_mulliken = [] if madelungenergies_mulliken is None else madelungenergies_mulliken + + if self.num_atoms is None: + # site_potentials + with zopen(filename, mode="rt") as file: + data = file.read().split("\n") + if len(data) == 0: + raise OSError("SitePotentials file contains no data.") + + self._filename = filename + self.ewald_splitting = float(data[0].split()[9]) + + data = data[5:-1] + self.num_atoms = len(data) - 2 + for atom in range(self.num_atoms): + line = data[atom].split() + self.atomlist.append(line[1] + str(line[0])) + self.types.append(line[1]) + self.sitepotentials_mulliken.append(float(line[2])) + self.sitepotentials_loewdin.append(float(line[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): """ @@ -1647,14 +1774,50 @@ def get_structure_with_site_potentials(self, structure_filename): Structure Object with Mulliken and Loewdin charges as site properties. """ struct = Structure.from_file(structure_filename) - Mulliken = self.sitepotentials_Mulliken - Loewdin = self.sitepotentials_Loewdin + mulliken = self.sitepotentials_mulliken + loewdin = self.sitepotentials_loewdin site_properties = { - "Mulliken Site Potentials (eV)": Mulliken, - "Loewdin Site Potentials (eV)": Loewdin, + "Mulliken Site Potentials (eV)": mulliken, + "Loewdin Site Potentials (eV)": loewdin, } return struct.copy(site_properties=site_properties) + @property + def sitepotentials_Mulliken(self): + warnings.warn( + "`sitepotentials_Mulliken` attribute is deprecated. Use `sitepotentials_mulliken` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.sitepotentials_mulliken + + @property + def sitepotentials_Loewdin(self): + warnings.warn( + "`sitepotentials_Loewdin` attribute is deprecated. Use `sitepotentials_loewdin` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.sitepotentials_loewdin + + @property + def madelungenergies_Mulliken(self): + warnings.warn( + "`madelungenergies_Mulliken` attribute is deprecated. Use `madelungenergies_mulliken` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.madelungenergies_mulliken + + @property + def madelungenergies_Loewdin(self): + warnings.warn( + "`madelungenergies_Loewdin` attribute is deprecated. Use `madelungenergies_loewdin` instead.", + DeprecationWarning, + stacklevel=2, + ) + return self.madelungenergies_loewdin + def get_orb_from_str(orbs): """ diff --git a/tests/io/lobster/test_inputs.py b/tests/io/lobster/test_inputs.py index eb445fe42ed..cb5790c2b0a 100644 --- a/tests/io/lobster/test_inputs.py +++ b/tests/io/lobster/test_inputs.py @@ -11,6 +11,7 @@ from pytest import approx from pymatgen.core.structure import Structure +from pymatgen.electronic_structure.cohp import IcohpCollection from pymatgen.electronic_structure.core import Orbital, Spin from pymatgen.io.lobster import ( Bandoverlaps, @@ -549,6 +550,16 @@ def test_values(self): assert self.icobi.icohpcollection.extremum_icohpvalue() == 0.58649 assert self.icobi_orbitalwise_spinpolarized.icohplist["2"]["orbitals"]["2s-6s"]["icohp"][Spin.up] == 0.0247 + def test_msonable(self): + dict_data = self.icobi_orbitalwise_spinpolarized.as_dict() + icohplist_from_dict = Icohplist.from_dict(dict_data) + all_attributes = vars(self.icobi_orbitalwise_spinpolarized) + for attr_name, attr_value in all_attributes.items(): + if isinstance(attr_value, IcohpCollection): + assert getattr(icohplist_from_dict, attr_name).as_dict() == attr_value.as_dict() + else: + assert getattr(icohplist_from_dict, attr_name) == attr_value + class TestNciCobiList(unittest.TestCase): def setUp(self): @@ -815,6 +826,13 @@ def test_get_structure_with_charges(self): s2 = Structure.from_dict(structure_dict2) assert s2 == self.charge2.get_structure_with_charges(f"{TEST_FILES_DIR}/POSCAR.MnO") + def test_msonable(self): + dict_data = self.charge2.as_dict() + charge_from_dict = Charge.from_dict(dict_data) + all_attributes = vars(self.charge2) + for attr_name, attr_value in all_attributes.items(): + assert getattr(charge_from_dict, attr_name) == attr_value + class TestLobsterout(PymatgenTest): def setUp(self): @@ -2062,6 +2080,13 @@ def test_has_good_quality(self): number_occ_bands_spin_up=1, limit_deviation=0.1 ) + def test_msonable(self): + dict_data = self.bandoverlaps2_new.as_dict() + bandoverlaps_from_dict = Bandoverlaps.from_dict(dict_data) + all_attributes = vars(self.bandoverlaps2_new) + for attr_name, attr_value in all_attributes.items(): + assert getattr(bandoverlaps_from_dict, attr_name) == attr_value + class TestGrosspop(unittest.TestCase): def setUp(self): @@ -2179,6 +2204,13 @@ def test_structure_with_grosspop(self): new_structure = self.grosspop1.get_structure_with_total_grosspop(f"{TEST_FILES_DIR}/cohp/POSCAR.SiO2") assert_allclose(new_structure.frac_coords, Structure.from_dict(struct_dict).frac_coords) + def test_msonable(self): + dict_data = self.grosspop1.as_dict() + grosspop_from_dict = Grosspop.from_dict(dict_data) + all_attributes = vars(self.grosspop1) + for attr_name, attr_value in all_attributes.items(): + assert getattr(grosspop_from_dict, attr_name) == attr_value + class TestUtils(PymatgenTest): def test_get_all_possible_basis_combinations(self): @@ -2330,6 +2362,13 @@ def test_get_structure(self): assert structure.site_properties["Loewdin Site Potentials (eV)"] == [-8.77, -17.08, 9.57, 9.57, 8.45] assert structure.site_properties["Mulliken Site Potentials (eV)"] == [-11.38, -19.62, 11.18, 11.18, 10.09] + def test_msonable(self): + dict_data = self.sitepotential.as_dict() + sitepotential_from_dict = SitePotential.from_dict(dict_data) + all_attributes = vars(self.sitepotential) + for attr_name, attr_value in all_attributes.items(): + assert getattr(sitepotential_from_dict, attr_name) == attr_value + class TestMadelungEnergies(PymatgenTest): def setUp(self) -> None: @@ -2340,6 +2379,13 @@ def test_attributes(self): assert self.madelungenergies.madelungenergies_Mulliken == approx(-40.02) assert self.madelungenergies.ewald_splitting == approx(3.14) + def test_msonable(self): + dict_data = self.madelungenergies.as_dict() + madelung_from_dict = MadelungEnergies.from_dict(dict_data) + all_attributes = vars(self.madelungenergies) + for attr_name, attr_value in all_attributes.items(): + assert getattr(madelung_from_dict, attr_name) == attr_value + class TestLobsterMatrices(PymatgenTest): def setUp(self) -> None: