Add keyword check_occu: bool = True to CifParser.get_structures() (#2836)

jonathanjdenney · pre-commit-ci[bot] · janosh · web-flow · commit 0dca5b1c035c · 2023-08-21T18:34:15.000-07:00
* Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Add files via upload * merging edits * Update test_fstar.py * Addressed comments from Matt * Update test_fstar.py * Add skip_checks to CifParser * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add tests to skip_checks add tests to skip_checks * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update test_skip_checks.py * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update test_skip_checks.py * make recommended changes * pre-commit auto-fixes * remove Unnecessary `else` after `return` statement * pre-commit auto-fixes * Update cif.py * pre-commit auto-fixes * removed folder test skip checks this folder must have gotten restored at some point. It has been deleted again. * Made recommended changes Shortened some if statements to contain less duplicates. Moved the skip_checks test to the right place. * pre-commit auto-fixes * Update test_cif.py * Update test_cif.py * deleted this version of fstar * pre-commit auto-fixes * pre-commit auto-fixes * Update cif.py * pre-commit auto-fixes * Create Skip_checks_test.cif * made requested edits * pre-commit auto-fixes * Update cif.py * add doc str Warning: Aphysical site occupancies are incompatible with many pymatgen features. * fix test_skip_checks and improve var names * Update cif.py * pre-commit auto-fixes * fix errors * pre-commit auto-fixes * Update test_cif.py * Update cif.py * Update cif.py * Update cif.py * Update cif.py * not sure why comp_dict construction was so complicated? * ensure test_skip_checks raises ValueError without occupancy_tolerance=1.5 * Update cif.py * pre-commit auto-fixes * remove code duplication in repeated if skip_occu_checks blocks * fix typos * fix mypy * reduce git diff * rename skip_occu_checks to check_occu --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Janosh Riebesell <janosh.riebesell@gmail.com>
diff --git a/pymatgen/apps/battery/plotter.py b/pymatgen/apps/battery/plotter.py
@@ -108,7 +108,7 @@ def get_plot(self, width=8, height=8, term_zero=True, ax: plt.Axes = None):
             ax.plot(x, y, "-", linewidth=2, label=label)
 
         ax.legend()
-        ax.set_xlabel(self._choose_best_x_lable(formula=formula, wion_symbol=wion_symbol))
+        ax.set_xlabel(self._choose_best_x_label(formula=formula, wion_symbol=wion_symbol))
         ax.set_ylabel("Voltage (V)")
         plt.tight_layout()
         return ax
@@ -157,7 +157,7 @@ def get_plotly_figure(
                 width=width,
                 height=height,
                 font=font_dict,
-                xaxis={"title": self._choose_best_x_lable(formula=formula, wion_symbol=wion_symbol)},
+                xaxis={"title": self._choose_best_x_label(formula=formula, wion_symbol=wion_symbol)},
                 yaxis={"title": "Voltage (V)"},
                 **kwargs,
             ),
@@ -166,7 +166,7 @@ def get_plotly_figure(
         fig.update_layout(template="plotly_white", title_x=0.5)
         return fig
 
-    def _choose_best_x_lable(self, formula, wion_symbol):
+    def _choose_best_x_label(self, formula, wion_symbol):
         if self.xaxis in {"capacity", "capacity_grav"}:
             return "Capacity (mAh/g)"
         if self.xaxis == "capacity_vol":
diff --git a/pymatgen/io/cif.py b/pymatgen/io/cif.py
@@ -13,7 +13,7 @@
 from io import StringIO
 from itertools import groupby
 from pathlib import Path
-from typing import TYPE_CHECKING, Literal
+from typing import TYPE_CHECKING, Any, Literal
 
 import numpy as np
 from monty.io import zopen
@@ -24,6 +24,7 @@
 from pymatgen.core.lattice import Lattice
 from pymatgen.core.operations import MagSymmOp, SymmOp
 from pymatgen.core.periodic_table import DummySpecies, Element, Species, get_el_sp
+from pymatgen.core.sites import PeriodicSite
 from pymatgen.core.structure import Structure
 from pymatgen.electronic_structure.core import Magmom
 from pymatgen.symmetry.analyzer import SpacegroupAnalyzer, SpacegroupOperations
@@ -629,9 +630,9 @@ def get_lattice(
 
         except KeyError:
             # Missing Key search for cell setting
-            for lattice_lable in ["_symmetry_cell_setting", "_space_group_crystal_system"]:
-                if data.data.get(lattice_lable):
-                    lattice_type = data.data.get(lattice_lable).lower()
+            for lattice_label in ["_symmetry_cell_setting", "_space_group_crystal_system"]:
+                if data.data.get(lattice_label):
+                    lattice_type = data.data.get(lattice_label).lower()
                     try:
                         required_args = getargspec(getattr(Lattice, lattice_type)).args
 
@@ -910,8 +911,10 @@ def _parse_symbol(self, sym):
 
         return parsed_sym
 
-    def _get_structure(self, data, primitive, symmetrized) -> Structure | None:
-        """Generate structure from part of the CIF."""
+    def _get_structure(
+        self, data: dict[str, Any], primitive: bool, symmetrized: bool, check_occu: bool = False
+    ) -> Structure | None:
+        """Generate structure from part of the cif."""
 
         def get_num_implicit_hydrogens(sym):
             num_h = {"Wat": 2, "wat": 2, "O-H": 1}
@@ -960,7 +963,7 @@ def get_matching_coord(coord):
             if oxi_states is not None:
                 o_s = oxi_states.get(symbol, 0)
                 # use _atom_site_type_symbol if possible for oxidation state
-                if "_atom_site_type_symbol" in data.data:
+                if "_atom_site_type_symbol" in data.data:  # type: ignore[attr-defined]
                     oxi_symbol = data["_atom_site_type_symbol"][idx]
                     o_s = oxi_states.get(oxi_symbol, o_s)
                 try:
@@ -979,19 +982,19 @@ def get_matching_coord(coord):
                 occu = str2float(data["_atom_site_occupancy"][idx])
             except (KeyError, ValueError):
                 occu = 1
-
-            if occu > 0:
+            # If check_occu is True or the occupancy is greater than 0, create comp_d
+            if check_occu or occu > 0:
                 coord = (x, y, z)
                 match = get_matching_coord(coord)
-                comp_d = {el: occu}
+                comp_dict = {el: max(occu, 1e-8)}
+
                 if num_h > 0:
-                    comp_d["H"] = num_h  # type: ignore
+                    comp_dict["H"] = num_h  # type: ignore
                     self.warnings.append(
-                        "Structure has implicit hydrogens defined, "
-                        "parsed structure unlikely to be suitable for use "
-                        "in calculations unless hydrogens added."
+                        "Structure has implicit hydrogens defined, parsed structure unlikely to be "
+                        "suitable for use in calculations unless hydrogens added."
                     )
-                comp = Composition(comp_d)
+                comp = Composition(comp_dict)
 
                 if not match:
                     coord_to_species[coord] = comp
@@ -1002,7 +1005,6 @@ def get_matching_coord(coord):
                     # disordered magnetic not currently supported
                     coord_to_magmoms[match] = None
                     labels[match] = label
-
         sum_occu = [
             sum(c.values()) for c in coord_to_species.values() if set(c.elements) != {Element("O"), Element("H")}
         ]
@@ -1060,7 +1062,7 @@ def get_matching_coord(coord):
 
                 # The following might be a more natural representation of equivalent indices,
                 # but is not in the format expect by SymmetrizedStructure:
-                #   equivalent_indices.append(list(range(len(allcoords), len(coords)+len(allcoords))))
+                #   equivalent_indices.append(list(range(len(all_coords), len(coords)+len(all_coords))))
                 # The above gives a list like:
                 #   [[0, 1, 2, 3], [4, 5, 6, 7, 8, 9, 10, 11]] where the
                 # integers are site indices, whereas the version used below will give a version like:
@@ -1076,6 +1078,7 @@ def get_matching_coord(coord):
                 all_labels.extend(new_labels)
 
             # rescale occupancies if necessary
+            all_species_noedit = all_species[:]  # save copy before scaling in case of check_occu=True, used below
             for idx, species in enumerate(all_species):
                 total_occu = sum(species.values())
                 if 1 < total_occu <= self._occupancy_tolerance:
@@ -1105,12 +1108,21 @@ def get_matching_coord(coord):
                 # TODO: extract Wyckoff labels (or other CIF attributes) and include as site_properties
                 wyckoffs = ["Not Parsed"] * len(struct)
 
-                # Names of space groups are likewise not parsed (again, not all CIFs will contain this information)
+                # space groups names are likewise not parsed (again, not all CIFs will contain this information)
                 # What is stored are the lists of symmetry operations used to generate the structure
                 # TODO: ensure space group labels are stored if present
                 sg = SpacegroupOperations("Not Parsed", -1, self.symmetry_operations)
+                struct = SymmetrizedStructure(struct, sg, equivalent_indices, wyckoffs)
 
-                return SymmetrizedStructure(struct, sg, equivalent_indices, wyckoffs)
+            if check_occu:
+                struct = Structure(lattice, all_species, all_coords, site_properties=site_properties, labels=all_labels)
+                for idx in range(len(struct)):
+                    struct[idx] = PeriodicSite(
+                        all_species_noedit[idx], all_coords[idx], lattice, properties=site_properties, skip_checks=True
+                    )
+
+            if symmetrized or check_occu:
+                return struct
 
             struct = struct.get_sorted_structure()
 
@@ -1124,7 +1136,11 @@ def get_matching_coord(coord):
         return None
 
     def get_structures(
-        self, primitive: bool = True, symmetrized: bool = False, on_error: Literal["ignore", "warn", "raise"] = "warn"
+        self,
+        primitive: bool = True,
+        symmetrized: bool = False,
+        check_occu: bool = True,
+        on_error: Literal["ignore", "warn", "raise"] = "warn",
     ) -> list[Structure]:
         """Return list of structures in CIF file.
 
@@ -1140,12 +1156,18 @@ def get_structures(
                 currently Wyckoff labels and space group labels or numbers are
                 not included in the generated SymmetrizedStructure, these will be
                 notated as "Not Parsed" or -1 respectively.
+            check_occu (bool): If False, site occupancy will not be checked, allowing unphysical
+                occupancy != 1. Useful for experimental results in which occupancy was allowed
+                to refine to unphysical values. Warning: unphysical site occupancies are incompatible
+                with many pymatgen features. Defaults to True.
             on_error ('ignore' | 'warn' | 'raise'): What to do in case of KeyError or ValueError
                 while parsing CIF file. Defaults to 'warn'.
 
         Returns:
             list[Structure]: All structures in CIF file.
         """
+        if check_occu:  # added in https://github.com/materialsproject/pymatgen/pull/2836
+            warnings.warn("Structures with unphysical site occupancies are not compatible with many pymatgen features.")
         if primitive and symmetrized:
             raise ValueError(
                 "Using both 'primitive' and 'symmetrized' arguments is not currently supported "
@@ -1155,7 +1177,7 @@ def get_structures(
         structures = []
         for idx, dct in enumerate(self._cif.data.values()):
             try:
-                struct = self._get_structure(dct, primitive, symmetrized)
+                struct = self._get_structure(dct, primitive, symmetrized, check_occu=check_occu)
                 if struct:
                     structures.append(struct)
             except (KeyError, ValueError) as exc:
diff --git a/tests/io/test_cif.py b/tests/io/test_cif.py
@@ -914,33 +914,46 @@ def test_empty_deque(self):
         with pytest.raises(ValueError, match="Invalid CIF file with no structures"):
             parser.get_structures()
 
+    def test_skip_checks(self):
+        with open(f"{TEST_FILES_DIR}/site_type_symbol_test.cif") as c:
+            cif_str = c.read()
+        cif_str = cif_str.replace("Te    Te 1.0000", "Te    Te 1.5000", 1)
+
+        with pytest.raises(ValueError, match="Invalid CIF file with no structures"):
+            # should fail without setting custom occupancy tolerance
+            CifParser.from_string(cif_str).get_structures()
+
+        parser = CifParser.from_string(cif_str, occupancy_tolerance=1.5)
+        structs = parser.get_structures(primitive=False, symmetrized=True, check_occu=False)[0]
+        assert structs[0].species.as_dict()["Te"] == 1.5
+
 
 class TestMagCif(PymatgenTest):
     def setUp(self):
         self.mcif = CifParser(f"{TEST_FILES_DIR}/magnetic.example.NiO.mcif")
         self.mcif_ncl = CifParser(f"{TEST_FILES_DIR}/magnetic.ncl.example.GdB4.mcif")
-        self.mcif_incom = CifParser(f"{TEST_FILES_DIR}/magnetic.incommensurate.example.Cr.mcif")
-        self.mcif_disord = CifParser(f"{TEST_FILES_DIR}/magnetic.disordered.example.CuMnO2.mcif")
+        self.mcif_incommensurate = CifParser(f"{TEST_FILES_DIR}/magnetic.incommensurate.example.Cr.mcif")
+        self.mcif_disordered = CifParser(f"{TEST_FILES_DIR}/magnetic.disordered.example.CuMnO2.mcif")
         self.mcif_ncl2 = CifParser(f"{TEST_FILES_DIR}/Mn3Ge_IR2.mcif")
 
     def test_mcif_detection(self):
         assert self.mcif.feature_flags["magcif"]
         assert self.mcif_ncl.feature_flags["magcif"]
-        assert self.mcif_incom.feature_flags["magcif"]
-        assert self.mcif_disord.feature_flags["magcif"]
+        assert self.mcif_incommensurate.feature_flags["magcif"]
+        assert self.mcif_disordered.feature_flags["magcif"]
         assert not self.mcif.feature_flags["magcif_incommensurate"]
         assert not self.mcif_ncl.feature_flags["magcif_incommensurate"]
-        assert self.mcif_incom.feature_flags["magcif_incommensurate"]
-        assert not self.mcif_disord.feature_flags["magcif_incommensurate"]
+        assert self.mcif_incommensurate.feature_flags["magcif_incommensurate"]
+        assert not self.mcif_disordered.feature_flags["magcif_incommensurate"]
 
     def test_get_structures(self):
         # incommensurate structures not currently supported
         with pytest.raises(NotImplementedError, match="Incommensurate structures not currently supported"):
-            self.mcif_incom.get_structures()
+            self.mcif_incommensurate.get_structures()
 
         # disordered magnetic structures not currently supported
         with pytest.raises(NotImplementedError, match="Disordered magnetic structures not currently supported"):
-            self.mcif_disord.get_structures()
+            self.mcif_disordered.get_structures()
 
         # taken from self.mcif_ncl, removing explicit magnetic symmops
         # so that MagneticSymmetryGroup() has to be invoked
