New class to handle NcICOBILIST.lobster files

pymatgen/io/lobster/__init__.py

@@ -18,6 +18,7 @@
     Icohplist,
     Lobsterout,
     MadelungEnergies,
+    Ncicobilist,
     SitePotential,
     Wavefunction,
 )

pymatgen/io/lobster/outputs.py

@@ -48,7 +48,8 @@
 
 
 class Cohpcar:
-    """Class to read COHPCAR/COOPCAR files generated by LOBSTER.
+    """
+    Class to read COHPCAR/COOPCAR/COBICAR files generated by LOBSTER.
 
     Attributes:
         cohp_data (dict[str, Dict[str, Any]]): A dictionary containing the COHP data of the form:
@@ -76,7 +77,7 @@ def __init__(self, are_coops: bool = False, are_cobis: bool = False, filename: s
         Args:
             are_coops: Determines if the file is a list of COHPs or COOPs.
               Default is False for COHPs.
-            are_cobis: Determines if the file is a list of COHPs or COOPs.
+            are_cobis: Determines if the file is a list of COHPs or COBIs.
               Default is False for COHPs.
 
             filename: Name of the COHPCAR file. If it is None, the default
@@ -408,6 +409,113 @@ def icohpcollection(self):
         return self._icohpcollection
 
 
+class Ncicobilist:
+    """
+    Class to read NcICOBILIST (multi-center ICOBI) files generated by LOBSTER.
+
+     Attributes:
+         is_spin_polarized (bool): Boolean to indicate if the calculation is spin polarized.
+         Ncicobilist (dict): Dict containing the listfile data of the form:
+           {bond: "number_of_atoms": number of atoms involved in the multi-center interaction,
+                  "ncicobi": {Spin.up: Nc-ICOBI(Ef) spin up, Spin.down: ...}},
+                  "interaction_type": type of the multi-center interaction
+    """
+
+    def __init__(self, filename: str | None = None):
+        """
+        Args:
+            filename: Name of the NcICOBILIST file.
+        """
+
+        if filename is None:
+            warnings.warn("Please consider using the newest LOBSTER version (4.1.0+). See http://www.cohp.de/.")
+            filename = "NcICOBILIST.lobster"
+
+        # LOBSTER list files have an extra trailing blank line
+        # and we don't need the header.
+        with zopen(filename, "rt") as f:
+            data = f.read().split("\n")[1:-1]
+        if len(data) == 0:
+            raise OSError("NcICOBILIST file contains no data.")
+
+        # If the calculation is spin polarized, the line in the middle
+        # of the file will be another header line.
+        if "spin" in data[len(data) // 2]:
+            # TODO: adapt this for orbitalwise stuff
+            self.is_spin_polarized = True
+        else:
+            self.is_spin_polarized = False
+
+        # check if orbitalwise NcICOBILIST
+        # include case when there is only one NcICOBI
+        for entry in data:  # NcICOBIs orbitalwise and non-orbitalwise can be mixed
+            if len(data) > 2 and "s]" in str(entry.split()[3:]):
+                self.orbitalwise = True
+                warnings.warn(
+                    "This is an orbitalwise NcICOBILIST.lobster file. Currently, the orbitalwise information is not "
+                    "read!"
+                )
+                break  # condition has only to be met once
+
+        if self.orbitalwise:
+            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.append(line)
+        else:
+            data_without_orbitals = data
+
+        if "spin" in data_without_orbitals[len(data_without_orbitals) // 2]:
+            # TODO: adapt this for orbitalwise stuff
+            num_bonds = len(data_without_orbitals) // 2
+            if num_bonds == 0:
+                raise OSError("NcICOBILIST file contains no data.")
+        else:
+            num_bonds = len(data_without_orbitals)
+
+        self.list_labels = []
+        self.list_numofatoms = []
+        self.list_ncicobi = []
+        self.list_interactiontype = []
+        self.list_num = []
+
+        for bond in range(num_bonds):
+            line = data_without_orbitals[bond].split()
+            ncicobi = {}
+
+            label = f"{line[0]}"
+            numofatoms = str(line[1])
+            ncicobi[Spin.up] = float(line[2])
+            interactiontype = str(line[3:]).replace("'", "").replace(" ", "")
+            num = 1
+
+            if self.is_spin_polarized:
+                ncicobi[Spin.down] = float(data_without_orbitals[bond + num_bonds + 1].split()[2])
+
+            self.list_labels.append(label)
+            self.list_numofatoms.append(numofatoms)
+            self.list_ncicobi.append(ncicobi)
+            self.list_interactiontype.append(interactiontype)
+            self.list_num.append(num)
+
+        # TODO: add functions to get orbital resolved NcICOBIs
+
+    @property
+    def ncicobilist(self) -> dict[Any, dict[str, Any]]:
+        """
+        Returns: ncicobilist.
+        """
+        ncicobilist = {}
+        for key, _entry in enumerate(self.list_labels):
+            ncicobilist[str(key + 1)] = {
+                "number_of_atoms": int(self.list_numofatoms[key]),
+                "ncicobi": self.list_ncicobi[key],
+                "interaction_type": self.list_interactiontype[key],
+            }
+
+        return ncicobilist
+
+
 class Doscar:
     """
     Class to deal with Lobster's projected DOS and local projected DOS.

test_files/cohp/NcICOBILIST.lobster

@@ -0,0 +1,50 @@
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 1  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X20 0 0]
+      2              3                      0.00009  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]
+      2              3                      0.00116  X22[5s]->Xs42[3s]->X31[5s]
+      2              3                      0.00001  X22[5s]->Xs42[3s]->X31[4d_xy]
+      2              3                     -0.00003  X22[5s]->Xs42[3s]->X31[4d_yz]
+      2              3                     -0.00000  X22[5s]->Xs42[3s]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3s]->X31[4d_xz]
+      2              3                      0.00001  X22[5s]->Xs42[3s]->X31[4d_x^2-y^2]
+      2              3                      0.00117  X22[5s]->Xs42[3p_y]->X31[5s]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_y]->X31[4d_xy]
+      2              3                     -0.00003  X22[5s]->Xs42[3p_y]->X31[4d_yz]
+      2              3                      0.00001  X22[5s]->Xs42[3p_y]->X31[4d_z^2]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_y]->X31[4d_xz]
+      2              3                      0.00003  X22[5s]->Xs42[3p_y]->X31[4d_x^2-y^2]
+      2              3                     -0.00071  X22[5s]->Xs42[3p_z]->X31[5s]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xy]
+      2              3                      0.00004  X22[5s]->Xs42[3p_z]->X31[4d_yz]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_z]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xz]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_z]->X31[4d_x^2-y^2]
+      2              3                     -0.00149  X22[5s]->Xs42[3p_x]->X31[5s]
+      2              3                      0.00000  X22[4d_x^2-y^2]->Xs42[3p]->X31[4d_x^2-y^2]
+      2              3                     -0.00007  X22[4d_x^2-y^2]->Xs42[3p_y]->X31[5s]
+      2              3                      0.00016  X22[4d_x^2-y^2]->Xs42[3p_x]->X31[5s]
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 2  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X2[0 0 0]
+      2              3                      0.00009  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]
+      2              3                      0.00116  X22[5s]->Xs42[3s]->X31[5s]
+      2              3                      0.00001  X22[5s]->Xs42[3s]->X31[4d_xy]
+      2              3                     -0.00003  X22[5s]->Xs42[3s]->X31[4d_yz]
+      2              3                     -0.00000  X22[5s]->Xs42[3s]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3s]->X31[4d_xz]
+      2              3                      0.00001  X22[5s]->Xs42[3s]->X31[4d_x^2-y^2]
+      2              3                      0.00117  X22[5s]->Xs42[3p_y]->X31[5s]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_y]->X31[4d_xy]
+      2              3                     -0.00003  X22[5s]->Xs42[3p_y]->X31[4d_yz]
+      2              3                      0.00001  X22[5s]->Xs42[3p_y]->X31[4d_z^2]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_y]->X31[4d_xz]
+      2              3                      0.00003  X22[5s]->Xs42[3p_y]->X31[4d_x^2-y^2]
+      2              3                     -0.00071  X22[5s]->Xs42[3p_z]->X31[5s]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xy]
+      2              3                      0.00004  X22[5s]->Xs42[3p_z]->X31[4d_yz]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_z]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xz]
+      2              3                     -0.00002  X22[5s]->Xs42[3p_z]->X31[4d_x^2-y^2]
+      2              3                     -0.00149  X22[5s]->Xs42[3p_x]->X31[5s]
+      2              3                      0.00000  X22[4d_x^2-y^2]->Xs42[3p]->X31[4d_x^2-y^2]
+      2              3                     -0.00007  X22[4d_x^2-y^2]->Xs42[3p_y]->X31[5s]
+      2              3                      0.00016  X22[4d_x^2-y^2]->Xs42[3p_x]->X31[5s]

test_files/cohp/NcICOBILIST.lobster.nospin

@@ -0,0 +1,25 @@
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 1  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X20 0 0]
+      2              3                      0.00018  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]
+      2              3                      0.00232  X22[5s]->Xs42[3s]->X31[5s]
+      2              3                      0.00002  X22[5s]->Xs42[3s]->X31[4d_xy]
+      2              3                     -0.00006  X22[5s]->Xs42[3s]->X31[4d_yz]
+      2              3                     -0.00000  X22[5s]->Xs42[3s]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3s]->X31[4d_xz]
+      2              3                      0.00002  X22[5s]->Xs42[3s]->X31[4d_x^2-y^2]
+      2              3                      0.00234  X22[5s]->Xs42[3p_y]->X31[5s]
+      2              3                     -0.00004  X22[5s]->Xs42[3p_y]->X31[4d_xy]
+      2              3                     -0.00006  X22[5s]->Xs42[3p_y]->X31[4d_yz]
+      2              3                      0.00002  X22[5s]->Xs42[3p_y]->X31[4d_z^2]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_y]->X31[4d_xz]
+      2              3                      0.00006  X22[5s]->Xs42[3p_y]->X31[4d_x^2-y^2]
+      2              3                     -0.00142  X22[5s]->Xs42[3p_z]->X31[5s]
+      2              3                     -0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xy]
+      2              3                      0.00008  X22[5s]->Xs42[3p_z]->X31[4d_yz]
+      2              3                     -0.00004  X22[5s]->Xs42[3p_z]->X31[4d_z^2]
+      2              3                      0.00000  X22[5s]->Xs42[3p_z]->X31[4d_xz]
+      2              3                     -0.00004  X22[5s]->Xs42[3p_z]->X31[4d_x^2-y^2]
+      2              3                     -0.00298  X22[5s]->Xs42[3p_x]->X31[5s]
+      2              3                      0.00000  X22[4d_x^2-y^2]->Xs42[3p]->X31[4d_x^2-y^2]
+      2              3                     -0.00014  X22[4d_x^2-y^2]->Xs42[3p_y]->X31[5s]
+      2              3                      0.00032  X22[4d_x^2-y^2]->Xs42[3p_x]->X31[5s]

test_files/cohp/NcICOBILIST.lobster.nospin.withoutorbitals

@@ -0,0 +1,3 @@
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 1  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X20 0 0]
+      2              3                      0.00018  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]

test_files/cohp/NcICOBILIST.lobster.withoutorbitals

@@ -0,0 +1,6 @@
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 1  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X20 0 0]
+      2              3                      0.00009  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]
+  COBI#   No. of atoms  Nc-ICOBI (at) eF for spin 2  Atoms for Nc-ICOBI
+      1              2                      0.00000  X1[0 0 0]->X2[0 0 0]
+      2              3                      0.00009  X22[0 0 0]->Xs42[0 0 0]->X31[0 0 0]

tests/io/lobster/test_inputs.py

@@ -23,6 +23,7 @@
     Lobsterin,
     Lobsterout,
     MadelungEnergies,
+    Ncicobilist,
     SitePotential,
     Wavefunction,
 )
@@ -670,6 +671,36 @@ def test_values(self):
         assert self.icobi_orbitalwise_spinpolarized.icohplist["2"]["orbitals"]["2s-6s"]["icohp"][Spin.up] == 0.0247
 
 
+class TestNcicobilist(unittest.TestCase):
+    def setUp(self):
+        self.ncicobi = Ncicobilist(filename=f"{TEST_FILES_DIR}/cohp/NcICOBILIST.lobster")
+        self.ncicobigz = Ncicobilist(filename=f"{TEST_FILES_DIR}/cohp/NcICOBILIST.lobster.gz")
+        self.ncicobinospin = Ncicobilist(filename=f"{TEST_FILES_DIR}/cohp/NcICOBILIST.lobster.nospin")
+        self.ncicobinospinwo = Ncicobilist(filename=f"{TEST_FILES_DIR}/cohp/NcICOBILIST.lobster.nospin.withoutorbitals")
+        self.ncicobiwo = Ncicobilist(filename=f"{TEST_FILES_DIR}/cohp/NcICOBILIST.lobster.withoutorbitals")
+
+    def test_ncicobilist(self):
+        assert self.ncicobi.is_spin_polarized
+        assert not self.ncicobinospin.is_spin_polarized
+        assert self.ncicobiwo.is_spin_polarized
+        assert not self.ncicobinospinwo.is_spin_polarized
+        assert self.ncicobi.orbitalwise
+        assert self.ncicobinospin.orbitalwise
+        assert not self.ncicobiwo.orbitalwise
+        assert not self.ncicobinospinwo.orbitalwise
+        assert len(self.ncicobi.ncicobilist) == 2
+        assert self.ncicobi.ncicobilist["2"]["number_of_atoms"] == 3
+        assert self.ncicobi.ncicobilist["2"]["ncicobi"][Spin.up] == approx(0.00009)
+        assert self.ncicobi.ncicobilist["2"]["ncicobi"][Spin.down] == approx(0.00009)
+        assert self.ncicobi.ncicobilist["2"]["interaction_type"] == "[X22[0,0,0]->Xs42[0,0,0]->X31[0,0,0]]"
+        assert self.ncicobi.ncicobilist["2"]["ncicobi"][Spin.up] == self.ncicobiwo.ncicobilist["2"]["ncicobi"][Spin.up]
+        assert self.ncicobi.ncicobilist["2"]["ncicobi"][Spin.up] == self.ncicobigz.ncicobilist["2"]["ncicobi"][Spin.up]
+        assert self.ncicobi.ncicobilist["2"]["interaction_type"] == self.ncicobigz.ncicobilist["2"]["interaction_type"]
+        assert sum(list(self.ncicobi.ncicobilist["2"]["ncicobi"].values())) == approx(
+            self.ncicobinospin.ncicobilist["2"]["ncicobi"][Spin.up]
+        )
+
+
 class TestDoscar(unittest.TestCase):
     def setUp(self):
         # first for spin polarized version