Add strict_anions option to MaterialsProject2020Compatibility (#3803

)

* Add `strict_anions` option to `MaterialsProject2020Compatibility`

* Add `require_exact` option and additional documentation

---------

Co-authored-by: Matthew Horton <matthew.horton@microsoft.com>

pymatgen/entries/compatibility.py

@@ -50,6 +50,22 @@
 MP2020_COMPAT_CONFIG = loadfn(f"{MODULE_DIR}/MP2020Compatibility.yaml")
 MP_COMPAT_CONFIG = loadfn(f"{MODULE_DIR}/MPCompatibility.yaml")
 
+# This was compiled by cross-referencing structures in Materials Project from exp_compounds.json.gz
+# used in the fitting of the MP2020 correction scheme, and applying the BVAnalyzer algorithm to
+# determine oxidation state. O and S are not included since these are treated separately.
+MP2020_ANION_OXIDATION_STATE_RANGES = {
+    "Br": (-1, -1),
+    "Cl": (-1, -1),
+    "F": (-1, -1),
+    "H": (-1, -1),
+    "I": (-1, -1),
+    "N": (-3, -2),
+    "Sb": (-3, -2),
+    "Se": (-2, -1),
+    "Si": (-4, -1),
+    "Te": (-2, -1),
+}
+
 assert (  # ping @janosh @rkingsbury on GitHub if this fails
     MP2020_COMPAT_CONFIG["Corrections"]["GGAUMixingCorrections"]["O"]
     == MP2020_COMPAT_CONFIG["Corrections"]["GGAUMixingCorrections"]["F"]
@@ -830,6 +846,10 @@ class MaterialsProject2020Compatibility(Compatibility):
     Materials Project input set parameters (see pymatgen.io.vasp.sets.MPRelaxSet). Using
     this compatibility scheme on calculations with different parameters is not valid.
 
+    The option `strict_anions` was added due to a bug. See PR #3803 (May 2024) for
+    related discussion. This behavior may change in subsequent versions as a more comprehensive
+    fix for this issue may be found.
+
     Note: While the correction scheme is largely composition-based, the energy corrections
     applied to ComputedEntry and ComputedStructureEntry can differ for O and S-containing
     structures if entry.data['oxidation_states'] is not populated or explicitly set. This
@@ -843,6 +863,7 @@ def __init__(
         self,
         compat_type: str = "Advanced",
         correct_peroxide: bool = True,
+        strict_anions: Literal["require_exact", "require_bound", "no_check"] = "require_bound",
         check_potcar: bool = True,
         check_potcar_hash: bool = False,
         config_file: str | None = None,
@@ -869,6 +890,15 @@ def __init__(
             correct_peroxide: Specify whether peroxide/superoxide/ozonide
                 corrections are to be applied or not. If false, all oxygen-containing
                 compounds are assigned the 'oxide' correction. Default: True
+            strict_anions: only apply the anion corrections to anions. The option
+                "require_exact" will only apply anion corrections in cases where the
+                anion oxidation state is between the oxidation states used
+                in the experimental fitting data. The option "require_bound" will
+                define an anion as any species with an oxidation state value of <= -1.
+                This prevents the anion correction from being applied to unrealistic
+                hypothetical structures containing large proportions of very electronegative
+                elements, thus artificially over-stabilizing the compound. Set to "no_check"
+                to restore the original behavior described in the associated publication. Default: True
             check_potcar (bool): Check that the POTCARs used in the calculation are consistent
                 with the Materials Project parameters. False bypasses this check altogether. Default: True
                 Can also be disabled globally by running `pmg config --add PMG_POTCAR_CHECKS false`.
@@ -893,6 +923,7 @@ def __init__(
 
         self.compat_type = compat_type
         self.correct_peroxide = correct_peroxide
+        self.strict_anions = strict_anions
         self.check_potcar = check_potcar
         self.check_potcar_hash = check_potcar_hash
 
@@ -955,7 +986,7 @@ def get_adjustments(self, entry: AnyComputedEntry) -> list[EnergyAdjustment]:
         comp = entry.composition
         rform = comp.reduced_formula
         # sorted list of elements, ordered by electronegativity
-        elements = sorted((el for el in comp.elements if comp[el] > 0), key=lambda el: el.X)
+        sorted_elements = sorted((el for el in comp.elements if comp[el] > 0), key=lambda el: el.X)
 
         # Skip single elements
         if len(comp) == 1:
@@ -1048,18 +1079,35 @@ def get_adjustments(self, entry: AnyComputedEntry) -> list[EnergyAdjustment]:
                 "only the most electronegative atom."
             )
 
-        for anion in "Br I Se Si Sb Te H N F Cl".split():
+        for anion in ("Br", "I", "Se", "Si", "Sb", "Te", "H", "N", "F", "Cl"):
             if Element(anion) in comp and anion in self.comp_correction:
                 apply_correction = False
+                oxidation_state = entry.data["oxidation_states"].get(anion, 0)
                 # if the oxidation_states key is not populated, only apply the correction if the anion
                 # is the most electronegative element
-                if entry.data["oxidation_states"].get(anion, 0) < 0:
+                if oxidation_state < 0:
                     apply_correction = True
+                    if self.strict_anions == "require_bound" and oxidation_state > -1:
+                        # This is not an anion. Noting that the rare case of a fractional
+                        # oxidation state in range [-1, 0] might be considered an anionic.
+                        # This could include suboxides or metal-rich pnictides, chalcogenides etc.
+                        # However! these cases are not included in the experimental fitting data
+                        # used for the correction scheme, and so there is no information for
+                        # whether the corrections are appropriate in this instance, and likely
+                        # may.
+                        apply_correction = False
                 else:
-                    most_electroneg = elements[-1].symbol
+                    most_electroneg = sorted_elements[-1].symbol
                     if anion == most_electroneg:
                         apply_correction = True
 
+                if self.strict_anions == "require_exact":
+                    apply_correction = False
+                    if (oxi_range := MP2020_ANION_OXIDATION_STATE_RANGES.get(anion)) and (
+                        oxi_range[0] <= oxidation_state <= oxi_range[1]
+                    ):
+                        apply_correction = True
+
                 if apply_correction:
                     adjustments.append(
                         CompositionEnergyAdjustment(
@@ -1074,7 +1122,7 @@ def get_adjustments(self, entry: AnyComputedEntry) -> list[EnergyAdjustment]:
         # GGA / GGA+U mixing scheme corrections
         calc_u = entry.parameters.get("hubbards")
         calc_u = defaultdict(int) if calc_u is None else calc_u
-        most_electroneg = elements[-1].symbol
+        most_electroneg = sorted_elements[-1].symbol
         u_corrections = self.u_corrections.get(most_electroneg, defaultdict(float))
         u_settings = self.u_settings.get(most_electroneg, defaultdict(float))
         u_errors = self.u_errors.get(most_electroneg, defaultdict(float))

tests/entries/test_compatibility.py

@@ -551,6 +551,20 @@ def setUp(self):
         self.compat = MaterialsProject2020Compatibility(check_potcar_hash=False)
         self.gga_compat = MaterialsProject2020Compatibility("GGA", check_potcar_hash=False)
 
+        self.entry_many_anions = ComputedEntry(
+            "CuI9",
+            -1,
+            0.0,
+            parameters={
+                "is_hubbard": False,
+                "run_type": "GGA",
+                "potcar_spec": [
+                    {"titel": "PAW_PBE Cu_pv 06Sep2000", "hash": "2d718b6be91068094207c9e861e11a89"},
+                    {"titel": "PAW_PBE I 08Apr2002", "hash": "f4ff16a495dd361ff5824ee61b418bb0"},
+                ],
+            },
+        )
+
     def test_process_entry(self):
         # Correct parameters
         assert self.compat.process_entry(self.entry1) is not None
@@ -1048,6 +1062,19 @@ def test_process_entry_with_oxidation_state(self):
         assert processed_entry.correction == approx(-6.084)
         assert processed_entry.energy == approx(-58.97 + -6.084)
 
+    def test_many_anions(self):
+        compat = MaterialsProject2020Compatibility(strict_anions="require_bound")
+        processed_entry = compat.process_entry(self.entry_many_anions)
+        assert processed_entry.energy == -1
+
+        compat = MaterialsProject2020Compatibility(strict_anions="no_check")
+        processed_entry = compat.process_entry(self.entry_many_anions)
+        assert processed_entry.energy == approx(-4.411)
+
+        compat = MaterialsProject2020Compatibility(strict_anions="require_exact")
+        processed_entry = compat.process_entry(self.entry_many_anions)
+        assert processed_entry.energy == -1
+
 
 class TestMITCompatibility(TestCase):
     def setUp(self):