Orbital-resolved icohplist

pymatgen/electronic_structure/cohp.py

@@ -358,7 +358,6 @@ def get_cohp_by_label(self, label, summed_spin_channels=False):
         Returns:
             Returns the COHP object to simplify plotting
         """
-        # TODO: possibly add option to only get one spin channel as Spin.down in the future!
         if label.lower() == "average":
             divided_cohp = self.cohp
             divided_icohp = self.icohp
@@ -834,7 +833,9 @@ class IcohpValue(MSONable):
 
     """
 
-    def __init__(self, label, atom1, atom2, length, translation, num, icohp, are_coops=False, are_cobis=False):
+    def __init__(
+        self, label, atom1, atom2, length, translation, num, icohp, are_coops=False, are_cobis=False, orbitals=None
+    ):
         """
         Args:
             label: label for the icohp
@@ -846,6 +847,9 @@ def __init__(self, label, atom1, atom2, length, translation, num, icohp, are_coo
             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]}}
+
         """
         if are_coops and are_cobis:
             raise ValueError("You cannot have info about COOPs and COBIs in the same file.")
@@ -858,6 +862,7 @@ def __init__(self, label, atom1, atom2, length, translation, num, icohp, are_coo
         self._translation = translation
         self._num = num
         self._icohp = icohp
+        self._orbitals = orbitals
         if Spin.down in self._icohp:
             self._is_spin_polarized = True
         else:
@@ -1004,6 +1009,20 @@ def icohpvalue(self, spin=Spin.up):
 
         return self._icohp[spin]
 
+    def icohpvalue_orbital(self, orbitals, spin=Spin.up):
+        """
+        Args:
+            orbitals: List of Orbitals or "str(Orbital1)-str(Orbital2)"
+            spin: Spin.up or Spin.down
+        Returns:
+            icohpvalue (float) 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):
+            orbitals = str(orbitals[0]) + "-" + str(orbitals[1])
+        return self._orbitals[orbitals]["icohp"][spin]
+
     @property
     def icohp(self):
         """
@@ -1016,23 +1035,38 @@ def icohp(self):
     @property
     def summed_icohp(self):
         """
-        Adds ICOHPs of both spin channels for spin polarized compounds
+        Sums ICOHPs of both spin channels for spin polarized compounds
         Returns:
              icohp value in eV
         """
         sum_icohp = self._icohp[Spin.down] + self._icohp[Spin.up] if self._is_spin_polarized else self._icohp[Spin.up]
         return sum_icohp
 
+    @property
+    def summed_orbital_icohp(self):
+        """
+        Sums orbitals-resolved ICOHPs of both spin channels for spin-plarized compounds
+        Returns:
+            {"str(Orbital1)-str(Ortibal2)": icohp value in eV}
+
+        """
+        orbital_icohp = {}
+        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]
+            )
+        return orbital_icohp
+
 
 class IcohpCollection(MSONable):
     """
     Class to store IcohpValues
 
     .. attribute:: are_coops
-        Boolean to indicate if this are ICOOPs
+        Boolean to indicate if these are ICOOPs
 
     .. attribute:: are_cobis
-        Boolean to indicate if this are ICOOPs
+        Boolean to indicate if these are ICOOPs
 
     .. attribute:: is_spin_polarized
         Boolean to indicate if the Lobster calculation was done spin polarized or not
@@ -1049,22 +1083,25 @@ def __init__(
         list_num,
         list_icohp,
         is_spin_polarized,
+        list_orb_icohp=None,
         are_coops=False,
         are_cobis=False,
     ):
         """
         Args:
-            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
-            are_coops: Boolean to indicate whether ICOOPs are stored
-            are_cobis: Boolean to indicate whether ICOBIs are stored
             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
         """
         if are_coops and are_cobis:
             raise ValueError("You cannot have info about COOPs and COBIs in the same file.")
@@ -1079,6 +1116,7 @@ def __init__(
         self._list_translation = list_translation
         self._list_num = list_num
         self._list_icohp = list_icohp
+        self._list_orb_icohp = list_orb_icohp
 
         for ilist, listel in enumerate(list_labels):
             self._icohplist[listel] = IcohpValue(
@@ -1091,6 +1129,7 @@ def __init__(
                 icohp=list_icohp[ilist],
                 are_coops=are_coops,
                 are_cobis=are_cobis,
+                orbitals=None if list_orb_icohp is None else list_orb_icohp[ilist],
             )
 
     def __str__(self):
@@ -1099,7 +1138,7 @@ def __str__(self):
             joinstr.append(str(value))
         return "\n".join(joinstr)
 
-    def get_icohp_by_label(self, label, summed_spin_channels=True, spin=Spin.up):
+    def get_icohp_by_label(self, label, summed_spin_channels=True, spin=Spin.up, orbitals=None):
         """
         get an icohp value for a certain bond as indicated by the label (bond labels starting by "1" as in
         ICOHPLIST/ICOOPLIST)
@@ -1108,16 +1147,22 @@ def get_icohp_by_label(self, label, summed_spin_channels=True, spin=Spin.up):
             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)"
         Returns:
             float describing ICOHP/ICOOP value
         """
         icohp_here = self._icohplist[label]
-        if icohp_here._is_spin_polarized:
+        if orbitals is None:
             if summed_spin_channels:
                 return icohp_here.summed_icohp
             return icohp_here.icohpvalue(spin)
-        return icohp_here.icohpvalue(spin)
+
+        if isinstance(orbitals, list):
+            orbitals = str(orbitals[0]) + "-" + str(orbitals[1])
+        if summed_spin_channels:
+            return icohp_here.summed_orbital_icohp[orbitals]
+
+        return icohp_here.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):
         """

pymatgen/electronic_structure/tests/test_cohp.py

@@ -208,7 +208,45 @@ def setUp(self):
             list_num=list_num,
             list_icohp=list_icoop,
         )
-
+        self.icohpcollection_orbitalwise = IcohpCollection.from_dict(
+            {
+                "@module": "pymatgen.electronic_structure.cohp",
+                "@class": "IcohpCollection",
+                "@version": None,
+                "list_labels": ["1", "2"],
+                "list_atom1": ["O5", "O5"],
+                "list_atom2": ["Ta2", "Ta2"],
+                "list_length": [1.99474, 1.99474],
+                "list_translation": [[0, 0, -1], [0, 0, 0]],
+                "list_num": [1, 1],
+                "list_icohp": [{Spin.up: 0.29324, Spin.down: 0.29324}, {Spin.up: 0.29324, Spin.down: 0.29324}],
+                "is_spin_polarized": True,
+                "list_orb_icohp": [
+                    {
+                        "2s-6s": {
+                            "icohp": {Spin.up: 0.0247, Spin.down: 0.0247},
+                            "orbitals": [[2, Orbital.s], [6, Orbital.s]],
+                        },
+                        "2s-5py": {
+                            "icohp": {Spin.up: 8e-05, Spin.down: 8e-05},
+                            "orbitals": [[2, Orbital.s], [5, Orbital.py]],
+                        },
+                    },
+                    {
+                        "2s-6s": {
+                            "icohp": {Spin.up: 0.0247, Spin.down: 0.0247},
+                            "orbitals": [[2, Orbital.s], [6, Orbital.s]],
+                        },
+                        "2s-5py": {
+                            "icohp": {Spin.up: 0.5, Spin.down: 0},
+                            "orbitals": [[2, Orbital.s], [5, Orbital.py]],
+                        },
+                    },
+                ],
+                "are_coops": False,
+                "are_cobis": True,
+            }
+        )
         # with spin polarization:
         list_atom2_sp = ["Fe7", "Fe9"]
         list_labels_sp = ["1", "2"]
@@ -278,6 +316,27 @@ def test_get_icohp_by_label(self):
         assert self.icohpcollection_Fe.get_icohp_by_label("1", summed_spin_channels=False, spin=Spin.down) == -0.19701
         assert self.icohpcollection_Fe.get_icohp_by_label("2", summed_spin_channels=False, spin=Spin.down) == -0.58279
 
+        # orbitalwise
+        assert self.icohpcollection_orbitalwise.get_icohp_by_label("1", orbitals="2s-6s") == 0.0494
+        assert (
+            self.icohpcollection_orbitalwise.get_icohp_by_label(
+                "1", orbitals="2s-6s", spin=Spin.up, summed_spin_channels=False
+            )
+            == 0.0247
+        )
+        assert (
+            self.icohpcollection_orbitalwise.get_icohp_by_label(
+                "1", orbitals="2s-6s", spin=Spin.down, summed_spin_channels=False
+            )
+            == 0.0247
+        )
+        assert (
+            self.icohpcollection_orbitalwise.get_icohp_by_label(
+                "2", orbitals="2s-5py", spin=Spin.up, summed_spin_channels=False
+            )
+            == 0.5
+        )
+
     def test_get_summed_icohp_by_label_list(self):
         # without spin polarization
         assert self.icohpcollection_KF.get_summed_icohp_by_label_list(
@@ -311,6 +370,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [0, -1, -1],
+            "orbitals": None,
         }
         icohpvalue["2"] = {
             "@module": "pymatgen.electronic_structure.cohp",
@@ -324,6 +384,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [-1, 0, -1],
+            "orbitals": None,
         }
         icohpvalue["3"] = {
             "@module": "pymatgen.electronic_structure.cohp",
@@ -337,6 +398,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [0, 0, -1],
+            "orbitals": None,
         }
         icohpvalue["4"] = {
             "@module": "pymatgen.electronic_structure.cohp",
@@ -350,6 +412,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [-1, -1, 0],
+            "orbitals": None,
         }
         icohpvalue["5"] = {
             "@module": "pymatgen.electronic_structure.cohp",
@@ -363,6 +426,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [0, -1, 0],
+            "orbitals": None,
         }
         icohpvalue["6"] = {
             "@module": "pymatgen.electronic_structure.cohp",
@@ -376,6 +440,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "@class": "IcohpValue",
             "atom1": "F1",
             "translation": [-1, 0, 0],
+            "orbitals": None,
         }
 
         dict_KF = self.icohpcollection_KF.get_icohp_dict_by_bondlengths(minbondlength=0.0, maxbondlength=8.0)
@@ -401,6 +466,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "length": 2.83189,
             "@class": "IcohpValue",
             "icohp": {Spin.up: -0.10218, Spin.down: -0.19701},
+            "orbitals": None,
         }
         icohpvalue_spin["2"] = {
             "num": 1,
@@ -414,6 +480,7 @@ def test_get_icohp_dict_by_bondlengths(self):
             "length": 2.45249,
             "@class": "IcohpValue",
             "icohp": {Spin.up: -0.28485, Spin.down: -0.58279},
+            "orbitals": None,
         }
 
         dict_Fe = self.icohpcollection_Fe.get_icohp_dict_by_bondlengths(minbondlength=0.0, maxbondlength=8.0)
@@ -446,6 +513,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "1",
             "icohp": {Spin.up: -0.40075},
+            "orbitals": None,
         }
         icohpvalue["2"] = {
             "translation": [-1, 0, -1],
@@ -459,6 +527,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "2",
             "icohp": {Spin.up: -0.40074},
+            "orbitals": None,
         }
         icohpvalue["3"] = {
             "translation": [0, 0, -1],
@@ -472,6 +541,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "3",
             "icohp": {Spin.up: -0.40079},
+            "orbitals": None,
         }
         icohpvalue["4"] = {
             "translation": [-1, -1, 0],
@@ -485,6 +555,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "4",
             "icohp": {Spin.up: -0.40079},
+            "orbitals": None,
         }
         icohpvalue["5"] = {
             "translation": [0, -1, 0],
@@ -498,6 +569,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "5",
             "icohp": {Spin.up: -0.40074},
+            "orbitals": None,
         }
         icohpvalue["6"] = {
             "translation": [-1, 0, 0],
@@ -511,6 +583,7 @@ def test_get_icohp_dict_of_site(self):
             "num": 1,
             "label": "6",
             "icohp": {Spin.up: -0.40075},
+            "orbitals": None,
         }
 
         dict_KF = self.icohpcollection_KF.get_icohp_dict_of_site(site=0)
@@ -592,6 +665,7 @@ def test_get_icohp_dict_of_site(self):
             "atom1": "Fe8",
             "atom2": "Fe7",
             "label": "1",
+            "orbitals": None,
             "@class": "IcohpValue",
             "num": 2,
         }
@@ -605,6 +679,7 @@ def test_get_icohp_dict_of_site(self):
             "atom1": "Fe8",
             "atom2": "Fe9",
             "label": "2",
+            "orbitals": None,
             "@class": "IcohpValue",
             "num": 1,
         }