Merge pull request #2350 from JaGeo/master

Correction of small bugs in Lobster module

pymatgen/electronic_structure/tests/test_cohp.py

@@ -795,7 +795,7 @@ def setUp(self):
 
         # spinpolarized case:
         filepath = os.path.join(test_dir, "environments", "COHPCAR.lobster.mp-190.gz")
-        structure = os.path.join(test_dir, "environments", "POSCAR.mp_190")
+        structure = os.path.join(test_dir, "environments", "POSCAR.mp_190.gz")
         self.cohp_lobster_spin_polarized = CompleteCohp.from_file(
             "lobster", filename=filepath, structure_file=structure
         )
@@ -1212,7 +1212,7 @@ def setUp(self):
         self.cohp_orb = CompleteCohp.from_file("lobster", filename=filepath, structure_file=structure)
 
         filepath = os.path.join(test_dir, "environments", "COHPCAR.lobster.mp-190.gz")
-        structure = os.path.join(test_dir, "environments", "POSCAR.mp_190")
+        structure = os.path.join(test_dir, "environments", "POSCAR.mp_190.gz")
         self.cohp_lobster_spin_polarized = CompleteCohp.from_file(
             "lobster", filename=filepath, structure_file=structure
         )

pymatgen/io/lobster/lobster_basis/BASIS_PBE_54_standard.yaml

@@ -83,7 +83,7 @@ BASIS:
   La: '4f 5d 5p 5s 6s '
   La_s: '4f 5d 5p 6s '
   Li: '2p 2s '
-  Li_sv: '1s 2s '
+  Li_sv: '1s 2s 2p '
   Lu: '4f 5d 5p 5s 6s '
   Lu_3: '5d 5p 6s '
   Mg: '3p 3s '

pymatgen/io/lobster/lobsterenv.py

@@ -294,24 +294,25 @@ def get_light_structure_environment(self, only_cation_environments=False, only_i
 
         return lse
 
-    def get_info_icohps_to_neighbors(self, isites=[], onlycation_isites=True):
+    def get_info_icohps_to_neighbors(self, isites=None, onlycation_isites=True):
         """
         this method will return information of cohps of neighbors
         Args:
-            isites: list of site ids, if isite==[], all isites will be used to add the icohps of the neighbors
+            isites: list of site ids, if isite==None, all isites will be used to add the icohps of the neighbors
             onlycation_isites: will only use cations, if isite==[]
 
 
         Returns:
             sum of icohps of neighbors to certain sites [given by the id in structure], number of bonds to this site,
             labels (from ICOHPLIST) for
             these bonds
-            [the latter is useful for plotting summed COHP plots]
+            [the latter is useful for plotting summed COHP plots],
+            list of the central isite for each label
         """
 
         if self.valences is None and onlycation_isites:
             raise ValueError("No valences are provided")
-        if isites == []:
+        if isites is None:
             if onlycation_isites:
                 isites = [i for i in range(len(self.structure)) if self.valences[i] >= 0.0]
             else:
@@ -322,6 +323,7 @@ def get_info_icohps_to_neighbors(self, isites=[], onlycation_isites=True):
         number_bonds = 0
         labels = []
         atoms = []
+        final_isites = []
         for ival, site in enumerate(self.structure):
             if ival in isites:
                 for keys, icohpsum in zip(self.list_keys[ival], self.list_icohps[ival]):
@@ -335,8 +337,8 @@ def get_info_icohps_to_neighbors(self, isites=[], onlycation_isites=True):
                         ]
                     )
                     number_bonds += 1
-
-        return summed_icohps, list_icohps, number_bonds, labels, atoms
+                    final_isites.append(ival)
+        return summed_icohps, list_icohps, number_bonds, labels, atoms, final_isites
 
     def plot_cohps_of_neighbors(
         self,
@@ -396,7 +398,7 @@ def plot_cohps_of_neighbors(
     def get_info_cohps_to_neighbors(
         self,
         path_to_COHPCAR="COHPCAR.lobster",
-        isites=[],
+        isites=None,
         only_bonds_to=None,
         onlycation_isites=True,
         per_bond=True,
@@ -408,7 +410,7 @@ def get_info_cohps_to_neighbors(
             path_to_COHPCAR: str, path to COHPCAR
             isites: list of int that indicate the number of the site
             only_bonds_to: list of str, e.g. ["O"] to only show cohps of anything to oxygen
-            onlycation_isites: if isites=[], only cation sites will be returned
+            onlycation_isites: if isites=None, only cation sites will be returned
             per_bond: will normalize per bond
             summed_spin_channels: will sum all spin channels
 
@@ -417,7 +419,7 @@ def get_info_cohps_to_neighbors(
 
         """
         # TODO: add options for orbital-resolved cohps
-        summed_icohps, list_icohps, number_bonds, labels, atoms = self.get_info_icohps_to_neighbors(
+        summed_icohps, list_icohps, number_bonds, labels, atoms, final_isites = self.get_info_icohps_to_neighbors(
             isites=isites, onlycation_isites=onlycation_isites
         )
 
@@ -452,21 +454,26 @@ def get_info_cohps_to_neighbors(
             )
 
         else:
-            # TODO: check if this is okay
             # labels of the COHPs that will be summed!
             # iterate through labels and atoms and check which bonds can be included
             new_labels = []
             new_atoms = []
             # print(labels)
             # print(atoms)
-            for label, atompair in zip(labels, atoms):
-                # durchlaufe only_bonds_to=[] und sage ja, falls eines der Labels in atompair ist, dann speichere
-                # new_label
+            for label, atompair, isite in zip(labels, atoms, final_isites):
                 present = False
-                # print(only_bonds_to)
                 for atomtype in only_bonds_to:
-                    if atomtype in (self._split_string(atompair[0])[0], self._split_string(atompair[1])[0]):
-                        present = True
+                    # This is necessary to identify also bonds between the same elements correctly!
+                    if str(self.structure[isite].species.elements[0]) != atomtype:
+                        if atomtype in (self._split_string(atompair[0])[0], self._split_string(atompair[1])[0]):
+                            present = True
+                    else:
+                        if (
+                            atomtype == self._split_string(atompair[0])[0]
+                            and atomtype == self._split_string(atompair[1])[0]
+                        ):
+                            present = True
+
                 if present:
                     new_labels.append(label)
                     new_atoms.append(atompair)
@@ -506,13 +513,13 @@ def _get_plot_label(self, atoms, per_bond):
             plotlabel = plotlabel + " (per bond)"
         return plotlabel
 
-    def get_info_icohps_between_neighbors(self, isites=[], onlycation_isites=True):
+    def get_info_icohps_between_neighbors(self, isites=None, onlycation_isites=True):
 
         """
         will return infos about interactions between neighbors of a certain atom
         Args:
-            isites: list of site ids, if isite==[], all isites will be used
-            onlycation_isites: will only use cations, if isite==[]
+            isites: list of site ids, if isite==None, all isites will be used
+            onlycation_isites: will only use cations, if isite==None
 
         Returns:
 
@@ -523,7 +530,7 @@ def get_info_icohps_between_neighbors(self, isites=[], onlycation_isites=True):
 
         if self.valences is None and onlycation_isites:
             raise ValueError("No valences are provided")
-        if isites == []:
+        if isites is None:
             if onlycation_isites:
                 isites = [i for i in range(len(self.structure)) if self.valences[i] >= 0.0]
             else: