Merge pull request #1423 from JaGeo/master

New methods for COHP classes

pymatgen/electronic_structure/cohp.py

@@ -180,6 +180,41 @@ def get_interpolated_value(self, energy, integrated=False):
                 raise ValueError("ICOHP is empty.")
         return inter
 
+    def has_antibnd_states_below_efermi(self, spin=None, limit=0.01):
+        """
+        Returns dict indicating if there are antibonding states below the Fermi level depending on the spin
+            spin: Spin
+            limit: -COHP smaller -limit will be considered.
+
+        """
+        warnings.warn("This method has not been tested on many examples. Check the parameter limit, pls!")
+
+        populations = self.cohp
+        number_energies_below_efermi = len([x for x in self.energies if x <= self.efermi])
+
+        if populations is None:
+            return None
+        elif spin is None:
+            dict_to_return = {}
+            for spin, cohpvalues in populations.items():
+                if (max(cohpvalues[0:number_energies_below_efermi])) > limit:
+                    dict_to_return[spin] = True
+                else:
+                    dict_to_return[spin] = False
+        else:
+            dict_to_return = {}
+            if isinstance(spin, int):
+                spin = Spin(spin)
+            elif isinstance(spin, str):
+                s = {"up": 1, "down": -1}[spin.lower()]
+                spin = Spin(s)
+            if (max(populations[spin][0:number_energies_below_efermi])) > limit:
+                dict_to_return[spin] = True
+            else:
+                dict_to_return[spin] = False
+
+        return dict_to_return
+
     @classmethod
     def from_dict(cls, d):
         """
@@ -385,6 +420,42 @@ def get_summed_cohp_by_label_list(self, label_list, divisor=1):
                     are_coops=first_cohpobject.are_coops,
                     icohp=divided_icohp)
 
+    def get_summed_cohp_by_label_and_orbital_list(self, label_list, orbital_list, divisor=1):
+        """
+        Returns 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
+        Returns:
+            Returns a COHP object including a summed COHP
+        """
+        # check if cohps are spinpolarized or not
+        first_cohpobject = self.get_orbital_resolved_cohp(label_list[0], orbital_list[0])
+        summed_cohp = first_cohpobject.cohp.copy()
+        summed_icohp = first_cohpobject.icohp.copy()
+        for ilabel, label in enumerate(label_list[1:], 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)
+            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)
+            if Spin.down in summed_icohp:
+                summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.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)
+        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)
+
+        return Cohp(efermi=first_cohpobject.efermi, energies=first_cohpobject.energies, cohp=divided_cohp,
+                    are_coops=first_cohpobject.are_coops,
+                    icohp=divided_icohp)
+
     def get_orbital_resolved_cohp(self, label, orbitals):
         """
         Get orbital-resolved COHP.
@@ -944,7 +1015,7 @@ def extremum_icohpvalue(self, summed_spin_channels=True, spin=Spin.up):
 
         if not self._is_spin_polarized:
             if spin == Spin.down:
-                raise Warning("This spin channel does not exist. I am switching to Spin.up")
+                warnings.warn("This spin channel does not exist. I am switching to Spin.up")
             spin = Spin.up
 
         for value in self._icohplist.values():

pymatgen/electronic_structure/tests/test_cohp.py

@@ -63,6 +63,12 @@ def test_str(self):
         self.assertEqual(self.cohp.__str__(), str_cohp)
         self.assertEqual(self.coop.__str__(), str_coop)
 
+    def test_antibnd_states_below_efermi(self):
+        self.assertDictEqual(self.cohp.has_antibnd_states_below_efermi(spin=None), {Spin.up: True, Spin.down: True})
+        self.assertDictEqual(self.cohp.has_antibnd_states_below_efermi(spin=None, limit=0.5),
+                             {Spin.up: False, Spin.down: False})
+        self.assertDictEqual(self.cohp.has_antibnd_states_below_efermi(spin=Spin.up, limit=0.5), {Spin.up: False})
+
 
 class IcohpValueTest(unittest.TestCase):
     def setUp(self):
@@ -576,6 +582,23 @@ def test_get_summed_cohp_by_label_list(self):
         self.assertEqual(self.cohp_orb.get_summed_cohp_by_label_list(["1", "1"]).cohp[Spin.up][300], 0.03392 * 2.0)
         self.assertEqual(self.cohp_orb.get_summed_cohp_by_label_list(["1", "1"], divisor=2).cohp[Spin.up][300], 0.03392)
 
+    def test_get_summed_cohp_by_label_and_orbital_list(self):
+        ref = self.cohp_orb.orb_res_cohp["1"]["4s-4px"]
+        ref2 = self.cohp_orb.orb_res_cohp["1"]["4px-4pz"]
+        cohp_label = self.cohp_orb.get_summed_cohp_by_label_and_orbital_list(["1"], ["4s-4px"])
+        cohp_label2 = self.cohp_orb.get_summed_cohp_by_label_and_orbital_list(["1", "1"], ["4s-4px", "4s-4px"])
+        cohp_label2x = self.cohp_orb.get_summed_cohp_by_label_and_orbital_list(["1", "1"], ["4s-4px", "4s-4px"],
+                                                                               divisor=2)
+        cohp_label3 = self.cohp_orb.get_summed_cohp_by_label_and_orbital_list(["1", "1"], ["4px-4pz", "4s-4px"])
+
+        self.assertArrayEqual(cohp_label.cohp[Spin.up], ref["COHP"][Spin.up])
+        self.assertArrayEqual(cohp_label2.cohp[Spin.up], ref["COHP"][Spin.up] * 2.0)
+        self.assertArrayEqual(cohp_label3.cohp[Spin.up], ref["COHP"][Spin.up] + ref2["COHP"][Spin.up])
+        self.assertArrayEqual(cohp_label.icohp[Spin.up], ref["ICOHP"][Spin.up])
+        self.assertArrayEqual(cohp_label2.icohp[Spin.up], ref["ICOHP"][Spin.up] * 2.0)
+        self.assertArrayEqual(cohp_label2x.icohp[Spin.up], ref["ICOHP"][Spin.up])
+        self.assertArrayEqual(cohp_label3.icohp[Spin.up], ref["ICOHP"][Spin.up] + ref2["ICOHP"][Spin.up])
+
     def test_orbital_resolved_cohp(self):
         # When read from a COHPCAR file, total COHPs are calculated from
         # the orbital-resolved COHPs if the total is missing. This may be