Repository: materialsproject/pymatgen

Issue:
Mistake in regex string of ion_or_solid_comp_object leads to example string formula "Na[+]" failing
### Python version

Python 3.12.3

### Pymatgen version

 2024.11.13

### Operating system version

Windows 10

### Current behavior

The following code leads to an error (ValueError: + is an invalid formula!)
```
from pymatgen.analysis.pourbaix_diagram import ion_or_solid_comp_object
fails = ion_or_solid_comp_object('Na[+]')
```
Even though the definition of ion_or_solid_comp_object states that the argument string should be formatted as 'Na[+]'
```
def ion_or_solid_comp_object(formula):
    """Get an Ion or Composition object given a formula.

    Args:
        formula: String formula. Eg. of ion: NaOH(aq), Na[+];
            Eg. of solid: Fe2O3(s), Fe(s), Na2O

    Returns:
        Composition/Ion object
    """
    if re.match(r"\[([^\[\]]+)\]|\(aq\)", formula):
        comp_obj = Ion.from_formula(formula)
    elif re.search(r"\(s\)", formula):
        comp_obj = Composition(formula[:-3])
    else:
        comp_obj = Composition(formula)
    return comp_obj
```
In it's current form, the string needs to be '[+]Na' in order for it to work:
```
from pymatgen.analysis.pourbaix_diagram import ion_or_solid_comp_object
works = ion_or_solid_comp_object('[+]Na')
```

### Expected Behavior

If 'Na[+]' is the desired string format, the definition of ion_or_solid_comp_object should be modified

L398 of pourbaix_diagram.py is currently:
```
if re.match(r"\[([^\[\]]+)\]|\(aq\)", formula):
```
https://github.com/materialsproject/pymatgen/blob/da607e86f9ce8aec942067c6c1a4fda6e04915dd/src/pymatgen/analysis/pourbaix_diagram.py#L398C18-L398C19

I suggest replacement with 
```
if re.match(r".*\[([^\[\]]+)\]|\(aq\)", formula):
```

I am not super familiar with git workflows and pull requests so I have sent this as an issue. Please advise if you think this is something that should be fixed and if you wish me to do a pull request :)

### Minimal example

```Python
from pymatgen.analysis.pourbaix_diagram import ion_or_solid_comp_object

works = ion_or_solid_comp_object('[+]Na')
fails = ion_or_solid_comp_object('Na[+]')
```


### Relevant files to reproduce this bug

_No response_

Hints from issue comments:
@DanielYang59: Hi @IsaacGentle thanks for reporting and I could confirm this issue. For a quick guide on making contributions to a project, you could checkout [the guide here](https://github.com/materialsproject/pymatgen/blob/master/docs/contributing.md).

I could confirm your fix indeed work, but things looks more complex that this as some other formulas still don't work (for example `Cl[-]`).

Also the unit test for this helper method seems missing, so I might take over from here and fix this :) Thanks again for reporting and investigating.

---

A small tip, to share a code block, you could use the [permalink](https://stackoverflow.com/a/31282559/24021108) instead of copying the code block.

Another small tip, to show differences in code, you could [add `diff`](https://stackoverflow.com/a/40883538/24021108), which would look like:
```diff
- if re.match(r"\[([^\[\]]+)\]|\(aq\)", formula):
+ if re.match(r".*\[([^\[\]]+)\]|\(aq\)", formula):
```

---

Update: looks like [`ion_or_solid_comp_object` is not used across the code base at all](https://github.com/search?q=repo%3Amaterialsproject%2Fpymatgen%20ion_or_solid_comp_object&type=code), not sure about its use case, and it looks like a private helper method to me
@IsaacGentle: Hi @DanielYang59, thanks a lot for your comprehensive reply and helpful tips! Interesting that [ion_or_solid_comp_object](https://github.com/search?q=repo%3Amaterialsproject%2Fpymatgen%20ion_or_solid_comp_object&type=code) is not used across the code base. To add extra information about how I discovered this minor bug, I was use the pymatgen scripts to generate Pourbaix diagrams, but I was manually adding some of the ionic species using [ion_or_solid_comp_object](https://github.com/search?q=repo%3Amaterialsproject%2Fpymatgen%20ion_or_solid_comp_object&type=code).@DanielYang59: No problem at all. Thanks for sharing the context (though it still looks more like a private function to me, but I wouldn't rename it to private for now as I'm not super familiar with `pourbaix_diagram` module)

Happy coding in the new year!


Solution patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index d57d1fa020..5dbb8143a6 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -8,7 +8,7 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.8.1
+    rev: v0.9.0
     hooks:
       - id: ruff
         args: [--fix, --unsafe-fixes]
@@ -22,7 +22,7 @@ repos:
       - id: trailing-whitespace
 
   - repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v1.13.0
+    rev: v1.14.1
     hooks:
       - id: mypy
 
@@ -65,6 +65,6 @@ repos:
         args: [--drop-empty-cells, --keep-output]
 
   - repo: https://github.com/RobertCraigie/pyright-python
-    rev: v1.1.389
+    rev: v1.1.391
     hooks:
       - id: pyright
diff --git a/docs/pymatgen.md b/docs/pymatgen.md
index a72f8e125b..90f1e934b2 100644
--- a/docs/pymatgen.md
+++ b/docs/pymatgen.md
@@ -2454,7 +2454,6 @@ nav_order: 6
       * [`ElementBase.is_noble_gas`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_noble_gas)
       * [`ElementBase.is_post_transition_metal`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_post_transition_metal)
       * [`ElementBase.is_quadrupolar`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_quadrupolar)
-      * [`ElementBase.is_rare_earth_metal`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_rare_earth_metal)
       * [`ElementBase.is_transition_metal`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_transition_metal)
       * [`ElementBase.is_valid_symbol()`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.is_valid_symbol)
       * [`ElementBase.iupac_ordering`](pymatgen.core.md#pymatgen.core.periodic_table.ElementBase.iupac_ordering)
diff --git a/pyproject.toml b/pyproject.toml
index cc7ff40208..3f0e090ec1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -242,6 +242,7 @@ docstring-code-format = true
 "src/pymatgen/analysis/*" = ["D"]
 "src/pymatgen/io/*" = ["D"]
 "dev_scripts/*" = ["D"]
+"src/pymatgen/util/string.py" = ["A005"]
 
 [tool.pytest.ini_options]
 addopts = "--durations=30 --quiet -r xXs --color=yes --import-mode=importlib"
diff --git a/src/pymatgen/analysis/phase_diagram.py b/src/pymatgen/analysis/phase_diagram.py
index 8e07887b11..2bf5271f3d 100644
--- a/src/pymatgen/analysis/phase_diagram.py
+++ b/src/pymatgen/analysis/phase_diagram.py
@@ -1805,7 +1805,7 @@ class PatchedPhaseDiagram(PhaseDiagram):
         Raises:
             ValueError: If no suitable PhaseDiagram is found for the entry.
         """
-        entry_space = frozenset(entry.elements) if isinstance(entry, Composition) else frozenset(entry.elements)
+        entry_space = frozenset(entry.elements)
 
         try:
             return self.pds[entry_space]
diff --git a/src/pymatgen/analysis/pourbaix_diagram.py b/src/pymatgen/analysis/pourbaix_diagram.py
index e5e9c003e7..15c9a407f2 100644
--- a/src/pymatgen/analysis/pourbaix_diagram.py
+++ b/src/pymatgen/analysis/pourbaix_diagram.py
@@ -31,11 +31,16 @@ from pymatgen.util.plotting import pretty_plot
 from pymatgen.util.string import Stringify
 
 if TYPE_CHECKING:
-    from typing import Any
+    from collections.abc import Sequence
+    from typing import Any, ClassVar, Literal
 
     import matplotlib.pyplot as plt
+    from numpy.typing import NDArray
     from typing_extensions import Self
 
+    from pymatgen.core import DummySpecies, Species
+    from pymatgen.entries.computed_entries import ComputedStructureEntry
+
 __author__ = "Sai Jayaraman"
 __copyright__ = "Copyright 2012, The Materials Project"
 __version__ = "0.4"
@@ -62,12 +67,15 @@ due.cite(
 
 logger = logging.getLogger(__name__)
 
-PREFAC = 0.0591
+PREFAC: float = 0.0591  # ln(10) * RT/nF in Nernst Equation
 
 
 # TODO: Revise to more closely reflect PDEntry, invoke from energy/composition
+
 # TODO: PourbaixEntries depend implicitly on having entry energies be
 # formation energies, should be a better way to get from raw energies
+
+
 # TODO: uncorrected_energy is a bit of a misnomer, but not sure what to rename
 class PourbaixEntry(MSONable, Stringify):
     """
@@ -82,7 +90,12 @@ class PourbaixEntry(MSONable, Stringify):
     work. This may be changed to be more flexible in the future.
     """
 
-    def __init__(self, entry, entry_id=None, concentration=1e-6):
+    def __init__(
+        self,
+        entry: ComputedEntry | ComputedStructureEntry,
+        entry_id: str | None = None,
+        concentration: float = 1e-6,
+    ) -> None:
         """
         Args:
             entry (ComputedEntry | ComputedStructureEntry | PDEntry | IonEntry): An entry object
@@ -100,29 +113,42 @@ class PourbaixEntry(MSONable, Stringify):
             self.charge = 0
         self.uncorrected_energy = entry.energy
         if entry_id is not None:
-            self.entry_id = entry_id
+            self.entry_id: str | None = entry_id
         elif getattr(entry, "entry_id", None):
             self.entry_id = entry.entry_id
         else:
             self.entry_id = None
 
+    def __repr__(self) -> str:
+        energy, npH, nPhi, nH2O, entry_id = (
+            self.energy,
+            self.npH,
+            self.nPhi,
+            self.nH2O,
+            self.entry_id,
+        )
+        return (
+            f"{type(self).__name__}({self.entry.composition} with {energy=:.4f}, {npH=}, "
+            f"{nPhi=}, {nH2O=}, {entry_id=})"
+        )
+
     @property
-    def npH(self):
+    def npH(self) -> float:
         """The number of H."""
         return self.entry.composition.get("H", 0) - 2 * self.entry.composition.get("O", 0)
 
     @property
-    def nH2O(self):
+    def nH2O(self) -> float:
         """The number of H2O."""
         return self.entry.composition.get("O", 0)
 
     @property
-    def nPhi(self):
+    def nPhi(self) -> float:
         """The number of electrons."""
         return self.npH - self.charge
 
     @property
-    def name(self):
+    def name(self) -> str:
         """The entry's name."""
         if self.phase_type == "Solid":
             return f"{self.entry.reduced_formula}(s)"
@@ -130,22 +156,22 @@ class PourbaixEntry(MSONable, Stringify):
         return self.entry.name
 
     @property
-    def energy(self):
+    def energy(self) -> float:
         """Total energy of the Pourbaix entry (at pH, V = 0 vs. SHE)."""
         # Note: this implicitly depends on formation energies as input
         return self.uncorrected_energy + self.conc_term - (MU_H2O * self.nH2O)
 
     @property
-    def energy_per_atom(self):
+    def energy_per_atom(self) -> float:
         """Energy per atom of the Pourbaix entry."""
         return self.energy / self.composition.num_atoms
 
     @property
-    def elements(self):
+    def elements(self) -> list[Element | Species | DummySpecies]:
         """Elements in the entry."""
         return self.entry.elements
 
-    def energy_at_conditions(self, pH, V):
+    def energy_at_conditions(self, pH: float, V: float) -> float:
         """Get free energy for a given pH and V.
 
         Args:
@@ -157,7 +183,7 @@ class PourbaixEntry(MSONable, Stringify):
         """
         return self.energy + self.npH * PREFAC * pH + self.nPhi * V
 
-    def get_element_fraction(self, element):
+    def get_element_fraction(self, element: Element | str) -> float:
         """Get the elemental fraction of a given non-OH element.
 
         Args:
@@ -170,13 +196,13 @@ class PourbaixEntry(MSONable, Stringify):
         return self.composition.get(element) * self.normalization_factor
 
     @property
-    def normalized_energy(self):
+    def normalized_energy(self) -> float:
         """Energy normalized by number of non H or O atoms, e.g.
         for Zn2O6, energy / 2 or for AgTe3(OH)3, energy / 4.
         """
         return self.energy * self.normalization_factor
 
-    def normalized_energy_at_conditions(self, pH, V):
+    def normalized_energy_at_conditions(self, pH: float, V: float) -> float:
         """Energy at an electrochemical condition, compatible with
         numpy arrays for pH/V input.
 
@@ -190,7 +216,7 @@ class PourbaixEntry(MSONable, Stringify):
         return self.energy_at_conditions(pH, V) * self.normalization_factor
 
     @property
-    def conc_term(self):
+    def conc_term(self) -> float:
         """The concentration contribution to the free energy. Should only be present
         when there are ions in the entry.
         """
@@ -202,14 +228,17 @@ class PourbaixEntry(MSONable, Stringify):
         Note that the pH, voltage, H2O factors are always calculated when
         constructing a PourbaixEntry object.
         """
-        dct = {"@module": type(self).__module__, "@class": type(self).__name__}
+        dct = {
+            "@module": type(self).__module__,
+            "@class": type(self).__name__,
+            "entry": self.entry.as_dict(),
+            "concentration": self.concentration,
+            "entry_id": self.entry_id,
+        }
         if isinstance(self.entry, IonEntry):
             dct["entry_type"] = "Ion"
         else:
             dct["entry_type"] = "Solid"
-        dct["entry"] = self.entry.as_dict()
-        dct["concentration"] = self.concentration
-        dct["entry_id"] = self.entry_id
         return dct
 
     @classmethod
@@ -224,17 +253,17 @@ class PourbaixEntry(MSONable, Stringify):
         return cls(entry, entry_id, concentration)
 
     @property
-    def normalization_factor(self):
+    def normalization_factor(self) -> float:
         """Sum of number of atoms minus the number of H and O in composition."""
         return 1.0 / (self.num_atoms - self.composition.get("H", 0) - self.composition.get("O", 0))
 
     @property
-    def composition(self):
+    def composition(self) -> Composition:
         """Composition."""
         return self.entry.composition
 
     @property
-    def num_atoms(self):
+    def num_atoms(self) -> float:
         """Number of atoms in current formula. Useful for normalization."""
         return self.composition.num_atoms
 
@@ -245,40 +274,27 @@ class PourbaixEntry(MSONable, Stringify):
 
         return self.entry.name
 
-    def __repr__(self):
-        energy, npH, nPhi, nH2O, entry_id = (
-            self.energy,
-            self.npH,
-            self.nPhi,
-            self.nH2O,
-            self.entry_id,
-        )
-        return (
-            f"{type(self).__name__}({self.entry.composition} with {energy=:.4f}, {npH=}, "
-            f"{nPhi=}, {nH2O=}, {entry_id=})"
-        )
-
 
 class MultiEntry(PourbaixEntry):
     """PourbaixEntry-like object for constructing multi-elemental Pourbaix diagrams."""
 
-    def __init__(self, entry_list, weights=None):
+    def __init__(self, entry_list: Sequence[PourbaixEntry], weights: list[float] | None = None) -> None:
         """Initialize a MultiEntry.
 
         Args:
-            entry_list ([PourbaixEntry]): List of component PourbaixEntries
-            weights ([float]): Weights associated with each entry. Default is None
+            entry_list (Sequence[PourbaixEntry]): Component PourbaixEntries.
+            weights (list[float]): Weights associated with each entry. Default is None
         """
         self.weights = weights or [1.0] * len(entry_list)
-        self.entry_list = entry_list
+        self.entry_list = list(entry_list)
 
-    def __getattr__(self, attr):
+    def __getattr__(self, attr: str) -> Any:
         """
         Because most of the attributes here are just weighted averages of the entry_list,
         we save some space by having a set of conditionals to define the attributes.
         """
         # Attributes that are weighted averages of entry attributes
-        if attr in [
+        if attr in {
             "energy",
             "npH",
             "nH2O",
@@ -287,7 +303,7 @@ class MultiEntry(PourbaixEntry):
             "composition",
             "uncorrected_energy",
             "elements",
-        ]:
+        }:
             # TODO: Composition could be changed for compat with sum
             start = Composition() if attr == "composition" else 0
             weighted_values = (
@@ -296,18 +312,13 @@ class MultiEntry(PourbaixEntry):
             return sum(weighted_values, start)
 
         # Attributes that are just lists of entry attributes
-        if attr in ["entry_id", "phase_type"]:
+        if attr in {"entry_id", "phase_type"}:
             return [getattr(entry, attr) for entry in self.entry_list]
 
         # normalization_factor, num_atoms should work from superclass
         return self.__getattribute__(attr)
 
-    @property
-    def name(self):
-        """MultiEntry name, i.e. the name of each entry joined by ' + '."""
-        return " + ".join(entry.name for entry in self.entry_list)
-
-    def __repr__(self):
+    def __repr__(self) -> str:
         energy, npH, nPhi, nH2O, entry_id = (
             self.energy,
             self.npH,
@@ -318,7 +329,12 @@ class MultiEntry(PourbaixEntry):
         cls_name, species = type(self).__name__, self.name
         return f"Pourbaix{cls_name}({energy=:.4f}, {npH=}, {nPhi=}, {nH2O=}, {entry_id=}, {species=})"
 
-    def as_dict(self):
+    @property
+    def name(self) -> str:
+        """MultiEntry name, i.e. the name of each entry joined by ' + '."""
+        return " + ".join(entry.name for entry in self.entry_list)
+
+    def as_dict(self) -> dict[str, Any]:
         """Get MSONable dict."""
         return {
             "@module": type(self).__module__,
@@ -367,6 +383,9 @@ class IonEntry(PDEntry):
         name = name or self.ion.reduced_formula
         super().__init__(composition=ion.composition, energy=energy, name=name, attribute=attribute)
 
+    def __repr__(self) -> str:
+        return f"IonEntry : {self.composition} with energy = {self.energy:.4f}"
+
     @classmethod
     def from_dict(cls, dct: dict) -> Self:
         """Get an IonEntry object from a dict."""
@@ -377,34 +396,30 @@ class IonEntry(PDEntry):
             dct.get("attribute"),
         )
 
-    def as_dict(self):
+    def as_dict(self) -> dict[Literal["ion", "energy", "name"], Any]:
         """Create a dict of composition, energy, and ion name."""
         return {"ion": self.ion.as_dict(), "energy": self.energy, "name": self.name}
 
-    def __repr__(self):
-        return f"IonEntry : {self.composition} with energy = {self.energy:.4f}"
-
 
-def ion_or_solid_comp_object(formula):
-    """Get an Ion or Composition object given a formula.
+def ion_or_solid_comp_object(formula: str) -> Composition | Ion:
+    """Get an Ion or Composition from a formula.
 
     Args:
-        formula: String formula. Eg. of ion: NaOH(aq), Na[+];
-            Eg. of solid: Fe2O3(s), Fe(s), Na2O
+        formula (str): Formula. E.g. of ion: NaOH(aq), Na[+], Na+;
+            E.g. of solid: Fe2O3(s), Fe(s), Na2O.
 
     Returns:
-        Composition/Ion object
+        Composition/Ion object.
     """
-    if re.match(r"\[([^\[\]]+)\]|\(aq\)", formula):
-        comp_obj = Ion.from_formula(formula)
-    elif re.search(r"\(s\)", formula):
-        comp_obj = Composition(formula[:-3])
-    else:
-        comp_obj = Composition(formula)
-    return comp_obj
+    # Formula for ion
+    if formula.endswith("(aq)") or re.search(r"\[.*\]$", formula) or "-" in formula or "+" in formula:
+        return Ion.from_formula(formula)
 
+    # Formula for solid
+    if formula.endswith("(s)"):
+        return Composition(formula[:-3])
 
-ELEMENTS_HO = {Element("H"), Element("O")}
+    return Composition(formula)
 
 
 # TODO: the solids filter breaks some of the functionality of the
@@ -418,6 +433,8 @@ ELEMENTS_HO = {Element("H"), Element("O")}
 class PourbaixDiagram(MSONable):
     """Create a Pourbaix diagram from entries."""
 
+    elements_ho: ClassVar[set[Element]] = {Element("H"), Element("O")}
+
     def __init__(
         self,
         entries: list[PourbaixEntry] | list[MultiEntry],
@@ -425,7 +442,7 @@ class PourbaixDiagram(MSONable):
         conc_dict: dict[str, float] | None = None,
         filter_solids: bool = True,
         nproc: int | None = None,
-    ):
+    ) -> None:
         """
         Args:
             entries ([PourbaixEntry] or [MultiEntry]): Entries list
@@ -451,7 +468,7 @@ class PourbaixDiagram(MSONable):
 
         # Get non-OH elements
         self.pbx_elts = list(
-            set(itertools.chain.from_iterable([entry.composition.elements for entry in entries])) - ELEMENTS_HO
+            set(itertools.chain.from_iterable([entry.composition.elements for entry in entries])) - self.elements_ho
         )
         self.dim = len(self.pbx_elts) - 1
 
@@ -463,7 +480,7 @@ class PourbaixDiagram(MSONable):
             self._unprocessed_entries = single_entries
             self._filtered_entries = single_entries
             self._conc_dict = None
-            self._elt_comp = {k: v for k, v in entries[0].composition.items() if k not in ELEMENTS_HO}
+            self._elt_comp = {k: v for k, v in entries[0].composition.items() if k not in self.elements_ho}
             self._multi_element = True
 
         # Process single entry inputs
@@ -483,7 +500,7 @@ class PourbaixDiagram(MSONable):
 
             # If a conc_dict is specified, override individual entry concentrations
             for entry in ion_entries:
-                ion_elts = list(set(entry.elements) - ELEMENTS_HO)
+                ion_elts = list(set(entry.elements) - self.elements_ho)
                 # TODO: the logic here for ion concentration setting is in two
                 # places, in PourbaixEntry and here, should be consolidated
                 if len(ion_elts) == 1:
@@ -512,15 +529,15 @@ class PourbaixDiagram(MSONable):
 
         self._stable_domains, self._stable_domain_vertices = self.get_pourbaix_domains(self._processed_entries)
 
-    def _convert_entries_to_points(self, pourbaix_entries):
+    def _convert_entries_to_points(self, pourbaix_entries: list[PourbaixEntry]) -> NDArray:
         """
         Args:
-            pourbaix_entries ([PourbaixEntry]): list of Pourbaix entries
+            pourbaix_entries (list[PourbaixEntry]): Pourbaix entries
                 to process into vectors in nph-nphi-composition space.
 
         Returns:
-            list of vectors, [[nph, nphi, e0, x1, x2, ..., xn-1]]
-            corresponding to each entry in nph-nphi-composition space
+            NDAarray: vectors as [[nph, nphi, e0, x1, x2, ..., xn-1]]
+                corresponding to each entry in nph-nphi-composition space
         """
         vecs = [
             [entry.npH, entry.nPhi, entry.energy] + [entry.composition.get(elt) for elt in self.pbx_elts[:-1]]
@@ -531,15 +548,18 @@ class PourbaixDiagram(MSONable):
         vecs *= norms
         return vecs
 
-    def _get_hull_in_nph_nphi_space(self, entries) -> tuple[list[PourbaixEntry], list[Simplex]]:
+    def _get_hull_in_nph_nphi_space(
+        self,
+        entries: list[PourbaixEntry],
+    ) -> tuple[list[PourbaixEntry], list[Simplex]]:
         """Generate convex hull of Pourbaix diagram entries in composition,
         npH, and nphi space. This enables filtering of multi-entries
         such that only compositionally stable combinations of entries
         are included.
 
         Args:
-            entries ([PourbaixEntry]): list of PourbaixEntries to construct
-                the convex hull
+            entries (list[PourbaixEntry]): PourbaixEntries to construct
+                the convex hull.
 
         Returns:
             tuple[list[PourbaixEntry], list[Simplex]]: PourbaixEntry list and stable
@@ -588,11 +608,15 @@ class PourbaixDiagram(MSONable):
 
         return min_entries, valid_facets
 
-    def _preprocess_pourbaix_entries(self, entries, nproc=None):
+    def _preprocess_pourbaix_entries(
+        self,
+        entries: list[PourbaixEntry],
+        nproc: int | None = None,
+    ) -> list[MultiEntry]:
         """Generate multi-entries for Pourbaix diagram.
 
         Args:
-            entries ([PourbaixEntry]): list of PourbaixEntries to preprocess
+            entries (list[PourbaixEntry]): PourbaixEntries to preprocess
                 into MultiEntries
             nproc (int): number of processes to be used in parallel
                 treatment of entry combos
@@ -605,23 +629,23 @@ class PourbaixDiagram(MSONable):
 
         min_entries, valid_facets = self._get_hull_in_nph_nphi_space(entries)
 
-        combos = []
+        combos: list[list[frozenset]] = []
         for facet in valid_facets:
             for idx in range(1, self.dim + 2):
-                these_combos = []
+                these_combos: list[frozenset] = []
                 for combo in itertools.combinations(facet, idx):
                     these_entries = [min_entries[i] for i in combo]
                     these_combos.append(frozenset(these_entries))
                 combos.append(these_combos)
 
-        all_combos = set(itertools.chain.from_iterable(combos))
+        all_combos: set | list = set(itertools.chain.from_iterable(combos))
 
-        list_combos = []
-        for idx in all_combos:
-            list_combos.append(list(idx))
+        list_combos: list = []
+        for combo in all_combos:
+            list_combos.append(list(combo))
         all_combos = list_combos
 
-        multi_entries = []
+        multi_entries: list = []
 
         # Parallel processing of multi-entry generation
         if nproc is not None:
@@ -637,7 +661,11 @@ class PourbaixDiagram(MSONable):
 
         return multi_entries
 
-    def _generate_multielement_entries(self, entries, nproc=None):
+    def _generate_multielement_entries(
+        self,
+        entries: list[PourbaixEntry],
+        nproc: int | None = None,
+    ) -> list[MultiEntry]:
         """
         Create entries for multi-element Pourbaix construction.
 
@@ -655,13 +683,13 @@ class PourbaixDiagram(MSONable):
         total_comp = Composition(self._elt_comp)
 
         # generate all combinations of compounds that have all elements
-        entry_combos = [itertools.combinations(entries, idx + 1) for idx in range(n_elems)]
-        entry_combos = itertools.chain.from_iterable(entry_combos)
+        entry_combos: list = [itertools.combinations(entries, idx + 1) for idx in range(n_elems)]
+        entry_combos = list(itertools.chain.from_iterable(entry_combos))
 
-        entry_combos = filter(lambda x: total_comp < MultiEntry(x).composition, entry_combos)
+        entry_combos = list(filter(lambda x: total_comp < MultiEntry(x).composition, entry_combos))
 
         # Generate and filter entries
-        processed_entries = []
+        processed_entries: list = []
         total = sum(comb(len(entries), idx + 1) for idx in range(n_elems))
         if total > 1e6:
             warnings.warn(
@@ -676,15 +704,17 @@ class PourbaixDiagram(MSONable):
             processed_entries = list(filter(bool, processed_entries))
         # Serial processing of multi-entry generation
         else:
-            for entry_combo in entry_combos:
-                processed_entry = self.process_multientry(entry_combo, total_comp)
+            for combo in entry_combos:
+                processed_entry = self.process_multientry(combo, total_comp)
                 if processed_entry is not None:
                     processed_entries.append(processed_entry)
 
         return processed_entries
 
     @staticmethod
-    def process_multientry(entry_list, prod_comp, coeff_threshold=1e-4):
+    def process_multientry(
+        entry_list: Sequence, prod_comp: Composition, coeff_threshold: float = 1e-4
+    ) -> MultiEntry | None:
         """Static method for finding a multientry based on
         a list of entries and a product composition.
         Essentially checks to see if a valid aqueous
@@ -693,7 +723,7 @@ class PourbaixDiagram(MSONable):
         with weights according to the coefficients if so.
 
         Args:
-            entry_list ([Entry]): list of entries from which to
+            entry_list (Sequence[Entry]): Entries from which to
                 create a MultiEntry
             prod_comp (Composition): composition constraint for setting
                 weights of MultiEntry
@@ -722,7 +752,10 @@ class PourbaixDiagram(MSONable):
             return None
 
     @staticmethod
-    def get_pourbaix_domains(pourbaix_entries, limits=None):
+    def get_pourbaix_domains(
+        pourbaix_entries: list[PourbaixEntry],
+        limits: list[list[float]] | None = None,
+    ) -> tuple[dict, dict]:
         """Get a set of Pourbaix stable domains (i.e. polygons) in
         pH-V space from a list of pourbaix_entries.
 
@@ -736,25 +769,27 @@ class PourbaixDiagram(MSONable):
         points.
 
         Args:
-            pourbaix_entries ([PourbaixEntry]): Pourbaix entries
+            pourbaix_entries (list[PourbaixEntry]): Pourbaix entries
                 with which to construct stable Pourbaix domains
-            limits ([[float]]): limits in which to do the pourbaix
+            limits (list[list[float]]): limits in which to do the pourbaix
                 analysis
 
         Returns:
-            Returns a dict of the form {entry: [boundary_points]}.
-            The list of boundary points are the sides of the N-1
-            dim polytope bounding the allowable ph-V range of each entry.
+            tuple[dict[PourbaixEntry, list], dict[PourbaixEntry, NDArray]:
+                The first dict is of form: {entry: [boundary_points]}.
+                The list of boundary points are the sides of the N-1
+                dim polytope bounding the allowable ph-V range of each entry.
         """
         if limits is None:
             limits = [[-2, 16], [-4, 4]]
 
         # Get hyperplanes
-        hyperplanes = [
-            np.array([-PREFAC * entry.npH, -entry.nPhi, 0, -entry.energy]) * entry.normalization_factor
-            for entry in pourbaix_entries
-        ]
-        hyperplanes = np.array(hyperplanes)
+        hyperplanes = np.array(
+            [
+                np.array([-PREFAC * entry.npH, -entry.nPhi, 0, -entry.energy]) * entry.normalization_factor
+                for entry in pourbaix_entries
+            ]
+        )
         hyperplanes[:, 2] = 1
 
         max_contribs = np.max(np.abs(hyperplanes), axis=0)
@@ -773,7 +808,7 @@ class PourbaixDiagram(MSONable):
         hs_int = HalfspaceIntersection(hs_hyperplanes, np.array(interior_point))
 
         # organize the boundary points by entry
-        pourbaix_domains = {entry: [] for entry in pourbaix_entries}
+        pourbaix_domains: dict[PourbaixEntry, list] = {entry: [] for entry in pourbaix_entries}
         for intersection, facet in zip(hs_int.intersections, hs_int.dual_facets, strict=True):
             for v in facet:
                 if v < len(pourbaix_entries):
@@ -782,18 +817,21 @@ class PourbaixDiagram(MSONable):
 
         # Remove entries with no Pourbaix region
         pourbaix_domains = {k: v for k, v in pourbaix_domains.items() if v}
-        pourbaix_domain_vertices = {}
+        pourbaix_domain_vertices: dict[PourbaixEntry, NDArray[float]] = {}
 
         for entry, points in pourbaix_domains.items():
             points = np.array(points)[:, :2]
             # Initial sort to ensure consistency
             points = points[np.lexsort(np.transpose(points))]
-            center = np.mean(points, axis=0)
-            points_centered = points - center
+            center: NDArray[float] = np.mean(points, axis=0)
+            points_centered: NDArray[float] = points - center
 
             # Sort points by cross product of centered points,
             # isn't strictly necessary but useful for plotting tools
-            points_centered = sorted(points_centered, key=cmp_to_key(lambda x, y: x[0] * y[1] - x[1] * y[0]))
+            points_centered = sorted(
+                points_centered,
+                key=cmp_to_key(lambda x, y: x[0] * y[1] - x[1] * y[0]),  # type: ignore[index]
+            )
             points = points_centered + center
 
             # Create simplices corresponding to Pourbaix boundary
@@ -803,7 +841,7 @@ class PourbaixDiagram(MSONable):
 
         return pourbaix_domains, pourbaix_domain_vertices
 
-    def find_stable_entry(self, pH, V):
+    def find_stable_entry(self, pH: float, V: float) -> PourbaixEntry:
         """Find stable entry at a pH,V condition.
 
         Args:
@@ -816,15 +854,20 @@ class PourbaixDiagram(MSONable):
         energies_at_conditions = [entry.normalized_energy_at_conditions(pH, V) for entry in self.stable_entries]
         return self.stable_entries[np.argmin(energies_at_conditions)]
 
-    def get_decomposition_energy(self, entry, pH, V):
+    def get_decomposition_energy(
+        self,
+        entry: PourbaixEntry,
+        pH: float,
+        V: float,
+    ) -> NDArray:
         """Find decomposition to most stable entries in eV/atom,
         supports vectorized inputs for pH and V.
 
         Args:
             entry (PourbaixEntry): PourbaixEntry corresponding to
                 compound to find the decomposition for
-            pH (float, list[float]): pH at which to find the decomposition
-            V (float, list[float]): voltage at which to find the decomposition
+            pH (float): pH at which to find the decomposition
+            V (float): voltage at which to find the decomposition
 
         Returns:
             Decomposition energy for the entry, i.e. the energy above
@@ -833,7 +876,7 @@ class PourbaixDiagram(MSONable):
         # Check composition consistency between entry and Pourbaix diagram:
         pbx_comp = Composition(self._elt_comp).fractional_composition
         entry_pbx_comp = Composition(
-            {elt: coeff for elt, coeff in entry.composition.items() if elt not in ELEMENTS_HO}
+            {elt: coeff for elt, coeff in entry.composition.items() if elt not in self.elements_ho}
         ).fractional_composition
         if entry_pbx_comp != pbx_comp:
             raise ValueError("Composition of stability entry does not match Pourbaix Diagram")
@@ -846,7 +889,7 @@ class PourbaixDiagram(MSONable):
         decomposition_energy /= entry.composition.num_atoms
         return decomposition_energy
 
-    def get_hull_energy(self, pH: float | list[float], V: float | list[float]) -> np.ndarray:
+    def get_hull_energy(self, pH: float | list[float], V: float | list[float]) -> NDArray:
         """Get the minimum energy of the Pourbaix "basin" that is formed
         from the stable Pourbaix planes. Vectorized.
 
@@ -860,7 +903,7 @@ class PourbaixDiagram(MSONable):
         all_gs = np.array([entry.normalized_energy_at_conditions(pH, V) for entry in self.stable_entries])
         return np.min(all_gs, axis=0)
 
-    def get_stable_entry(self, pH, V):
+    def get_stable_entry(self, pH: float, V: float) -> PourbaixEntry | MultiEntry:
         """Get the stable entry at a given pH, V condition.
 
         Args:
@@ -875,26 +918,26 @@ class PourbaixDiagram(MSONable):
         return self.stable_entries[np.argmin(all_gs)]
 
     @property
-    def stable_entries(self):
+    def stable_entries(self) -> list:
         """The stable entries in the Pourbaix diagram."""
         return list(self._stable_domains)
 
     @property
-    def unstable_entries(self):
+    def unstable_entries(self) -> list:
         """All unstable entries in the Pourbaix diagram."""
         return [entry for entry in self.all_entries if entry not in self.stable_entries]
 
     @property
-    def all_entries(self):
+    def all_entries(self) -> list:
         """All entries used to generate the Pourbaix diagram."""
         return self._processed_entries
 
     @property
-    def unprocessed_entries(self):
+    def unprocessed_entries(self) -> list:
         """Unprocessed entries."""
         return self._unprocessed_entries
 
-    def as_dict(self):
+    def as_dict(self) -> dict[str, Any]:
         """Get MSONable dict."""
         return {
             "@module": type(self).__module__,
@@ -926,14 +969,33 @@ class PourbaixDiagram(MSONable):
 class PourbaixPlotter:
     """A plotter class for phase diagrams."""
 
-    def __init__(self, pourbaix_diagram):
+    def __init__(self, pourbaix_diagram: PourbaixDiagram) -> None:
         """
         Args:
             pourbaix_diagram (PourbaixDiagram): A PourbaixDiagram object.
         """
         self._pbx = pourbaix_diagram
 
-    def show(self, *args, **kwargs):
+    @staticmethod
+    def _generate_entry_label(entry: PourbaixEntry | MultiEntry) -> str:
+        """
+        Generates a label for the Pourbaix plotter.
+
+        Args:
+            entry (PourbaixEntry or MultiEntry): entry to get a label for
+        """
+        if isinstance(entry, MultiEntry):
+            return " + ".join(entry.name for entry in entry.entry_list)
+
+        # TODO - a more elegant solution could be added later to Stringify
+        # for example, the pattern re.sub(r"([-+][\d\.]*)", r"$^{\1}$", )
+        # will convert B(OH)4- to B(OH)$_4^-$.
+        # for this to work, the ion's charge always must be written AFTER
+        # the sign (e.g., Fe+2 not Fe2+)
+        string = entry.to_latex_string()
+        return re.sub(r"()\[([^)]*)\]", r"\1$^{\2}$", string)
+
+    def show(self, *args, **kwargs) -> None:
         """Show the Pourbaix plot.
 
         Args:
@@ -1003,7 +1065,7 @@ class PourbaixPlotter:
 
             if label_domains:
                 ax.annotate(
-                    generate_entry_label(entry),
+                    self._generate_entry_label(entry),
                     center,
                     ha="center",
                     va="center",
@@ -1032,12 +1094,12 @@ class PourbaixPlotter:
 
         Args:
             entry (Any): The entry to plot stability for.
-            pH_range (tuple[float, float], optional): pH range for the plot. Defaults to (-2, 16).
-            pH_resolution (int, optional): pH resolution. Defaults to 100.
-            V_range (tuple[float, float], optional): Voltage range for the plot. Defaults to (-3, 3).
-            V_resolution (int, optional): Voltage resolution. Defaults to 100.
-            e_hull_max (float, optional): Maximum energy above the hull. Defaults to 1.
-            cmap (str, optional): Colormap for the plot. Defaults to "RdYlBu_r".
+            pH_range (tuple[float, float]): pH range for the plot. Defaults to (-2, 16).
+            pH_resolution (int): pH resolution. Defaults to 100.
+            V_range (tuple[float, float]): Voltage range for the plot. Defaults to (-3, 3).
+            V_resolution (int): Voltage resolution. Defaults to 100.
+            e_hull_max (float): Maximum energy above the hull. Defaults to 1.
+            cmap (str): Colormap for the plot. Defaults to "RdYlBu_r".
             ax (Axes, optional): Existing matplotlib Axes object for plotting. Defaults to None.
             **kwargs (Any): Additional keyword arguments passed to `get_pourbaix_plot`.
 
@@ -1056,7 +1118,7 @@ class PourbaixPlotter:
         # Plot stability map
         cax = ax.pcolor(pH, V, stability, cmap=cmap, vmin=0, vmax=e_hull_max)
         cbar = ax.figure.colorbar(cax)
-        cbar.set_label(f"Stability of {generate_entry_label(entry)} (eV/atom)")
+        cbar.set_label(f"Stability of {self._generate_entry_label(entry)} (eV/atom)")
 
         # Set ticklabels
         # ticklabels = [t.get_text() for t in cbar.ax.get_yticklabels()]
@@ -1065,7 +1127,7 @@ class PourbaixPlotter:
 
         return ax
 
-    def domain_vertices(self, entry):
+    def domain_vertices(self, entry) -> list:
         """Get the vertices of the Pourbaix domain.
 
         Args:
@@ -1075,22 +1137,3 @@ class PourbaixPlotter:
             list of vertices
         """
         return self._pbx._stable_domain_vertices[entry]
-
-
-def generate_entry_label(entry):
-    """
-    Generates a label for the Pourbaix plotter.
-
-    Args:
-        entry (PourbaixEntry or MultiEntry): entry to get a label for
-    """
-    if isinstance(entry, MultiEntry):
-        return " + ".join(entry.name for entry in entry.entry_list)
-
-    # TODO - a more elegant solution could be added later to Stringify
-    # for example, the pattern re.sub(r"([-+][\d\.]*)", r"$^{\1}$", )
-    # will convert B(OH)4- to B(OH)$_4^-$.
-    # for this to work, the ion's charge always must be written AFTER
-    # the sign (e.g., Fe+2 not Fe2+)
-    string = entry.to_latex_string()
-    return re.sub(r"()\[([^)]*)\]", r"\1$^{\2}$", string)
diff --git a/src/pymatgen/core/periodic_table.py b/src/pymatgen/core/periodic_table.py
index d3ae2e44d4..7bddfbe6df 100644
--- a/src/pymatgen/core/periodic_table.py
+++ b/src/pymatgen/core/periodic_table.py
@@ -14,7 +14,6 @@ from pathlib import Path
 from typing import TYPE_CHECKING, overload
 
 import numpy as np
-from monty.dev import deprecated
 from monty.json import MSONable
 
 from pymatgen.core.units import SUPPORTED_UNIT_NAMES, FloatWithUnit, Ha_to_eV, Length, Mass, Unit
@@ -754,19 +753,6 @@ class ElementBase(Enum):
         """
         return self.is_lanthanoid or self.is_actinoid or self.symbol in {"Sc", "Y"}
 
-    @property
-    @deprecated(
-        is_rare_earth,
-        message="is_rare_earth is corrected to include Y and Sc.",
-        deadline=(2025, 1, 1),
-    )
-    def is_rare_earth_metal(self) -> bool:
-        """True if element is a rare earth metal, Lanthanides (La) series and Actinides (Ac) series.
-
-        This property is Deprecated, and scheduled for removal after 2025-01-01.
-        """
-        return self.is_lanthanoid or self.is_actinoid
-
     @property
     def is_metal(self) -> bool:
         """True if is a metal."""
diff --git a/src/pymatgen/io/ase.py b/src/pymatgen/io/ase.py
index f3cf8516d6..308128c798 100644
--- a/src/pymatgen/io/ase.py
+++ b/src/pymatgen/io/ase.py
@@ -6,7 +6,6 @@ Atoms object and pymatgen Structure objects.
 from __future__ import annotations
 
 import warnings
-from collections.abc import Iterable
 from importlib.metadata import PackageNotFoundError
 from typing import TYPE_CHECKING
 
@@ -18,7 +17,7 @@ from pymatgen.core.structure import Lattice, Molecule, Structure
 try:
     from ase.atoms import Atoms
     from ase.calculators.singlepoint import SinglePointDFTCalculator
-    from ase.constraints import FixAtoms
+    from ase.constraints import FixAtoms, FixCartesian
     from ase.io.jsonio import decode, encode
     from ase.spacegroup import Spacegroup
 
@@ -26,7 +25,7 @@ try:
 
 except ImportError:
     NO_ASE_ERR = PackageNotFoundError("AseAtomsAdaptor requires the ASE package. Use `pip install ase`")
-    encode = decode = FixAtoms = SinglePointDFTCalculator = Spacegroup = None
+    encode = decode = FixAtoms = FixCartesian = SinglePointDFTCalculator = Spacegroup = None
 
     class Atoms:  # type: ignore[no-redef]
         def __init__(self, *args, **kwargs):
@@ -157,31 +156,38 @@ class AseAtomsAdaptor:
         # Get the oxidation states from the structure
         oxi_states: list[float | None] = [getattr(site.specie, "oxi_state", None) for site in structure]
 
-        # Read in selective dynamics if present. Note that the ASE FixAtoms class fixes (x,y,z), so
-        # here we make sure that [False, False, False] or [True, True, True] is set for the site selective
-        # dynamics property. As a consequence, if only a subset of dimensions are fixed, this won't get passed to ASE.
+        # Read in selective dynamics if present.
+        # Note that FixCartesian class uses an opposite notion of
+        # "fix" and "not fix" flags: in ASE True means fixed and False
+        # means not fixed.
+        fix_atoms: dict | None = None
         if "selective_dynamics" in structure.site_properties:
-            fix_atoms = []
+            fix_atoms = {
+                str([xc, yc, zc]): ([xc, yc, zc], [])
+                for xc in [True, False]
+                for yc in [True, False]
+                for zc in [True, False]
+            }
+            # [False, False, False] is free to move - no constraint in ASE.
+            del fix_atoms[str([False, False, False])]
             for site in structure:
                 selective_dynamics: ArrayLike = site.properties.get("selective_dynamics")  # type: ignore[assignment]
-                if (
-                    isinstance(selective_dynamics, Iterable)
-                    and True in selective_dynamics
-                    and False in selective_dynamics
-                ):
-                    raise ValueError(
-                        "ASE FixAtoms constraint does not support selective dynamics in only some dimensions. "
-                        f"Remove the {selective_dynamics=} and try again if you do not need them."
-                    )
-                is_fixed = bool(~np.all(site.properties["selective_dynamics"]))
-                fix_atoms.append(is_fixed)
-
+                for cmask_str in fix_atoms:
+                    cmask_site = (~np.array(selective_dynamics)).tolist()
+                    fix_atoms[cmask_str][1].append(cmask_str == str(cmask_site))
         else:
             fix_atoms = None
 
-        # Set the selective dynamics with the FixAtoms class.
+        # Set the selective dynamics with the FixCartesian class.
         if fix_atoms is not None:
-            atoms.set_constraint(FixAtoms(mask=fix_atoms))
+            atoms.set_constraint(
+                [
+                    FixAtoms(indices) if cmask == [True, True, True] else FixCartesian(indices, mask=cmask)
+                    for cmask, indices in fix_atoms.values()
+                    # Do not add empty constraints
+                    if any(indices)
+                ]
+            )
 
         # Add any remaining site properties to the ASE Atoms object
         for prop in structure.site_properties:
@@ -254,21 +260,40 @@ class AseAtomsAdaptor:
         oxi_states = atoms.get_array("oxi_states") if atoms.has("oxi_states") else None
 
         # If the ASE Atoms object has constraints, make sure that they are of the
-        # kind FixAtoms, which are the only ones that can be supported in Pymatgen.
+        # FixAtoms or FixCartesian kind, which are the only ones that
+        # can be supported in Pymatgen.
         # By default, FixAtoms fixes all three (x, y, z) dimensions.
         if atoms.constraints:
             unsupported_constraint_type = False
-            constraint_indices = []
+            constraint_indices: dict = {
+                str([xc, yc, zc]): ([xc, yc, zc], [])
+                for xc in [True, False]
+                for yc in [True, False]
+                for zc in [True, False]
+            }
             for constraint in atoms.constraints:
                 if isinstance(constraint, FixAtoms):
-                    constraint_indices.extend(constraint.get_indices().tolist())
+                    constraint_indices[str([False] * 3)][1].extend(constraint.get_indices().tolist())
+                elif isinstance(constraint, FixCartesian):
+                    cmask = (~np.array(constraint.mask)).tolist()
+                    constraint_indices[str(cmask)][1].extend(constraint.get_indices().tolist())
                 else:
                     unsupported_constraint_type = True
             if unsupported_constraint_type:
                 warnings.warn(
-                    "Only FixAtoms is supported by Pymatgen. Other constraints will not be set.", stacklevel=2
+                    "Only FixAtoms and FixCartesian is supported by Pymatgen. Other constraints will not be set.",
+                    stacklevel=2,
                 )
-            sel_dyn = [[False] * 3 if atom.index in constraint_indices else [True] * 3 for atom in atoms]
+            sel_dyn = []
+            for atom in atoms:
+                constrained = False
+                for mask, indices in constraint_indices.values():
+                    if atom.index in indices:
+                        sel_dyn.append(mask)
+                        constrained = True
+                        break  # Assume no duplicates
+                if not constrained:
+                    sel_dyn.append([False] * 3)
         else:
             sel_dyn = None
 
diff --git a/src/pymatgen/io/exciting/inputs.py b/src/pymatgen/io/exciting/inputs.py
index 7ef99482b0..41bc7ef09a 100644
--- a/src/pymatgen/io/exciting/inputs.py
+++ b/src/pymatgen/io/exciting/inputs.py
@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import itertools
+import warnings
 import xml.etree.ElementTree as ET
 from typing import TYPE_CHECKING
 
@@ -16,10 +17,13 @@ from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
 from pymatgen.symmetry.bandstructure import HighSymmKpath
 
 if TYPE_CHECKING:
-    from pathlib import Path
+    from typing import ClassVar, Literal
 
+    from numpy.typing import ArrayLike, NDArray
     from typing_extensions import Self
 
+    from pymatgen.util.typing import PathLike
+
 __author__ = "Christian Vorwerk"
 __copyright__ = "Copyright 2016"
 __version__ = "1.0"
@@ -37,10 +41,18 @@ class ExcitingInput(MSONable):
     Attributes:
         structure (Structure): Associated Structure.
         title (str): Optional title string.
-        lockxyz (numpy.ndarray): Lockxyz attribute for each site if available. A Nx3 array of booleans.
+        lockxyz (Nx3 NDArray of booleans): Lockxyz attribute for each site if available.
     """
 
-    def __init__(self, structure: Structure, title=None, lockxyz=None):
+    # Conversion factor between Bohr radius and Angstrom
+    bohr2ang: ClassVar[float] = const.value("Bohr radius") / const.value("Angstrom star")
+
+    def __init__(
+        self,
+        structure: Structure,
+        title: str | None = None,
+        lockxyz: ArrayLike | None = None,
+    ) -> None:
         """
         Args:
             structure (Structure): Structure object.
@@ -52,22 +64,19 @@ class ExcitingInput(MSONable):
         if structure.is_ordered:
             site_properties = {}
             if lockxyz:
-                site_properties["selective_dynamics"] = lockxyz
+                site_properties["selective_dynamics"] = np.asarray(lockxyz)
             self.structure = structure.copy(site_properties=site_properties)
             self.title = structure.formula if title is None else title
         else:
             raise ValueError("Structure with partial occupancies cannot be converted into exciting input!")
 
-    # define conversion factor between Bohr radius and Angstrom
-    bohr2ang = const.value("Bohr radius") / const.value("Angstrom star")
-
     @property
-    def lockxyz(self):
+    def lockxyz(self) -> NDArray:
         """Selective dynamics site properties."""
         return self.structure.site_properties.get("selective_dynamics")
 
     @lockxyz.setter
-    def lockxyz(self, lockxyz):
+    def lockxyz(self, lockxyz: ArrayLike) -> NDArray:
         self.structure.add_site_property("selective_dynamics", lockxyz)
 
     @classmethod
@@ -164,7 +173,7 @@ class ExcitingInput(MSONable):
         return cls(structure_in, title_in, lockxyz)
 
     @classmethod
-    def from_file(cls, filename: str | Path) -> Self:
+    def from_file(cls, filename: PathLike) -> Self:
         """
         Args:
             filename: Filename.
@@ -178,32 +187,27 @@ class ExcitingInput(MSONable):
 
     def write_etree(
         self,
-        celltype,
-        cartesian=False,
-        bandstr=False,
+        celltype: Literal["unchanged", "conventional", "primitive"],
+        cartesian: bool = False,
+        bandstr: bool = False,
         symprec: float = 0.4,
-        angle_tolerance=5,
+        angle_tolerance: float = 5,
         **kwargs,
-    ):
-        """Write the exciting input parameters to an XML object.
+    ) -> ET.Element:
+        """Convert the exciting input parameters to an XML Element object.
 
         Args:
             celltype (str): Choice of unit cell. Can be either the unit cell
                 from self.structure ("unchanged"), the conventional cell
                 ("conventional"), or the primitive unit cell ("primitive").
-
             cartesian (bool): Whether the atomic positions are provided in
                 Cartesian or unit-cell coordinates. Default is False.
-
             bandstr (bool): Whether the bandstructure path along the
                 HighSymmKpath is included in the input file. Only supported if the
                 celltype is set to "primitive". Default is False.
-
             symprec (float): Tolerance for the symmetry finding. Default is 0.4.
-
             angle_tolerance (float): Angle tolerance for the symmetry finding.
-            Default is 5.
-
+                Default is 5.
             **kwargs: Additional parameters for the input file.
 
         Returns:
@@ -297,77 +301,68 @@ class ExcitingInput(MSONable):
 
     def write_string(
         self,
-        celltype,
-        cartesian=False,
-        bandstr=False,
+        celltype: Literal["unchanged", "conventional", "primitive"],
+        cartesian: bool = False,
+        bandstr: bool = False,
         symprec: float = 0.4,
-        angle_tolerance=5,
+        angle_tolerance: float = 5,
         **kwargs,
-    ):
-        """Write exciting input.xml as a string.
+    ) -> str:
+        """Convert exciting input to a string.
 
         Args:
             celltype (str): Choice of unit cell. Can be either the unit cell
-            from self.structure ("unchanged"), the conventional cell
-            ("conventional"), or the primitive unit cell ("primitive").
-
+                from self.structure ("unchanged"), the conventional cell
+                ("conventional"), or the primitive unit cell ("primitive").
             cartesian (bool): Whether the atomic positions are provided in
-            Cartesian or unit-cell coordinates. Default is False.
-
+                Cartesian or unit-cell coordinates. Default is False.
             bandstr (bool): Whether the bandstructure path along the
-            HighSymmKpath is included in the input file. Only supported if the
-            celltype is set to "primitive". Default is False.
-
+                HighSymmKpath is included in the input file. Only supported if the
+                celltype is set to "primitive". Default is False.
             symprec (float): Tolerance for the symmetry finding. Default is 0.4.
-
             angle_tolerance (float): Angle tolerance for the symmetry finding.
-            Default is 5.
-
+                Default is 5.
             **kwargs: Additional parameters for the input file.
 
         Returns:
-            String
+            str
         """
         try:
             root = self.write_etree(celltype, cartesian, bandstr, symprec, angle_tolerance, **kwargs)
             self._indent(root)
             # output should be a string not a bytes object
             string = ET.tostring(root).decode("UTF-8")
+
         except Exception:
             raise ValueError("Incorrect celltype!")
+
         return string
 
     def write_file(
         self,
-        celltype,
-        filename,
-        cartesian=False,
-        bandstr=False,
+        celltype: Literal["unchanged", "conventional", "primitive"],
+        filename: PathLike,
+        cartesian: bool = False,
+        bandstr: bool = False,
         symprec: float = 0.4,
-        angle_tolerance=5,
+        angle_tolerance: float = 5,
         **kwargs,
-    ):
+    ) -> None:
         """Write exciting input file.
 
         Args:
             celltype (str): Choice of unit cell. Can be either the unit cell
-            from self.structure ("unchanged"), the conventional cell
-            ("conventional"), or the primitive unit cell ("primitive").
-
-            filename (str): Filename for exciting input.
-
+                from self.structure ("unchanged"), the conventional cell
+                ("conventional"), or the primitive unit cell ("primitive").
+            filename (PathLike): Filename for exciting input.
             cartesian (bool): Whether the atomic positions are provided in
-            Cartesian or unit-cell coordinates. Default is False.
-
+                Cartesian or unit-cell coordinates. Default is False.
             bandstr (bool): Whether the bandstructure path along the
-            HighSymmKpath is included in the input file. Only supported if the
-            celltype is set to "primitive". Default is False.
-
+                HighSymmKpath is included in the input file. Only supported if the
+                celltype is set to "primitive". Default is False.
             symprec (float): Tolerance for the symmetry finding. Default is 0.4.
-
             angle_tolerance (float): Angle tolerance for the symmetry finding.
-            Default is 5.
-
+                Default is 5.
             **kwargs: Additional parameters for the input file.
         """
         try:
@@ -378,15 +373,15 @@ class ExcitingInput(MSONable):
         except Exception:
             raise ValueError("Incorrect celltype!")
 
-    # Missing PrettyPrint option in the current version of xml.etree.cElementTree
     @staticmethod
-    def _indent(elem, level=0):
+    def _indent(elem: ET.Element, level: int = 0) -> None:
         """
-        Helper method to indent elements.
+        Helper method to indent elements, as missing PrettyPrint option
+        in the current version of xml.etree.cElementTree.
 
         Args:
-            elem:
-            level:
+            elem (ET.Element): The Element to process.
+            level (int): The indentation level.
         """
         i = "\n" + level * "  "
         if len(elem):
@@ -401,19 +396,23 @@ class ExcitingInput(MSONable):
         elif level and (not elem.tail or not elem.tail.strip()):
             elem.tail = i
 
-    def _dicttoxml(self, paramdict_, element):
+    def _dicttoxml(self, paramdict_: dict, element: ET.Element) -> None:
         for key, value in paramdict_.items():
             if isinstance(value, str) and key == "text()":
                 element.text = value
+
             elif isinstance(value, str):
                 element.attrib[key] = value
+
             elif isinstance(value, list):
                 for item in value:
                     self._dicttoxml(item, ET.SubElement(element, key))
+
             elif isinstance(value, dict):
                 if element.findall(key) == []:
                     self._dicttoxml(value, ET.SubElement(element, key))
                 else:
                     self._dicttoxml(value, element.findall(key)[0])
+
             else:
-                print("cannot deal with", key, "=", value)
+                warnings.warn(f"cannot deal with {key} = {value}", stacklevel=2)
diff --git a/src/pymatgen/io/vasp/MP24RelaxSet.yaml b/src/pymatgen/io/vasp/MP24RelaxSet.yaml
new file mode 100644
index 0000000000..72de8bef92
--- /dev/null
+++ b/src/pymatgen/io/vasp/MP24RelaxSet.yaml
@@ -0,0 +1,31 @@
+# Default VASP settings for calculations in the Materials Project
+# using the Regularized-Restored Strongly Constrained and Appropriately
+# Normed functional (r2SCAN).
+PARENT: PBE64Base
+INCAR:
+  ALGO: Normal
+  EDIFF: 1.e-05
+  EDIFFG: -0.02
+  ENAUG: 1360
+  ENCUT: 680
+  GGA_COMPAT: False
+  IBRION: 2
+  ISIF: 3
+  ISMEAR: 0 # included to have some reasonable default
+  ISPIN: 2
+  KSPACING: 0.22 # included to have some reasonable default
+  LAECHG: True
+  LASPH: True
+  LCHARG: True
+  LELF: False # LELF = True restricts calculation to KPAR = 1
+  LMAXMIX: 6 # per benchmark
+  LMIXTAU: True
+  LORBIT: 11
+  LREAL: False # per benchmark
+  LVTOT: True
+  LWAVE: False
+  METAGGA: R2SCAN
+  NELM: 200
+  NSW: 99
+  PREC: Accurate
+  SIGMA: 0.05 # included to have some reasonable default
diff --git a/src/pymatgen/io/vasp/PBE64Base.yaml b/src/pymatgen/io/vasp/PBE64Base.yaml
index af0236253e..358944a2b2 100644
--- a/src/pymatgen/io/vasp/PBE64Base.yaml
+++ b/src/pymatgen/io/vasp/PBE64Base.yaml
@@ -11,7 +11,7 @@ POTCAR:
   At: At
   Au: Au
   B: B
-  Ba: Ba_sv
+  Ba: Ba_sv_GW
   Be: Be_sv
   Bi: Bi
   Br: Br
@@ -26,8 +26,8 @@ POTCAR:
   Cr: Cr_pv
   Cs: Cs_sv
   Cu: Cu_pv
-  Dy: Dy_3
-  Er: Er_3
+  Dy: Dy_h
+  Er: Er_h
   Eu: Eu
   F: F
   Fe: Fe_pv
@@ -39,7 +39,7 @@ POTCAR:
   He: He
   Hf: Hf_pv
   Hg: Hg
-  Ho: Ho_3
+  Ho: Ho_h
   I: I
   In: In_d
   Ir: Ir
@@ -54,7 +54,7 @@ POTCAR:
   N: N
   Na: Na_pv
   Nb: Nb_pv
-  Nd: Nd_3
+  Nd: Nd_h
   Ne: Ne
   Ni: Ni_pv
   Np: Np
@@ -64,9 +64,9 @@ POTCAR:
   Pa: Pa
   Pb: Pb_d
   Pd: Pd
-  Pm: Pm_3
+  Pm: Pm_h
   Po: Po_d
-  Pr: Pr_3
+  Pr: Pr_h
   Pt: Pt
   Pu: Pu
   Ra: Ra_sv
@@ -80,22 +80,22 @@ POTCAR:
   Sc: Sc_sv
   Se: Se
   Si: Si
-  Sm: Sm_3
+  Sm: Sm_h
   Sn: Sn_d
   Sr: Sr_sv
   Ta: Ta_pv
-  Tb: Tb_3
+  Tb: Tb_h
   Tc: Tc_pv
   Te: Te
   Th: Th
   Ti: Ti_pv
   Tl: Tl_d
-  Tm: Tm_3
+  Tm: Tm_h
   U: U
   V: V_pv
   W: W_sv
-  Xe: Xe
+  Xe: Xe_GW
   Y: Y_sv
-  Yb: Yb_3
+  Yb: Yb_h
   Zn: Zn
   Zr: Zr_sv
diff --git a/src/pymatgen/io/vasp/sets.py b/src/pymatgen/io/vasp/sets.py
index 11e5e55581..d23a861ca3 100644
--- a/src/pymatgen/io/vasp/sets.py
+++ b/src/pymatgen/io/vasp/sets.py
@@ -84,7 +84,10 @@ MODULE_DIR = os.path.dirname(__file__)
 
 
 def _load_yaml_config(fname):
-    config = loadfn(f"{MODULE_DIR}/{fname}.yaml")
+    fname = f"{MODULE_DIR}/{fname}"
+    if not fname.endswith(".yaml"):
+        fname += ".yaml"
+    config = loadfn(fname)
     if "PARENT" in config:
         parent_config = _load_yaml_config(config["PARENT"])
         for k, v in parent_config.items():
@@ -623,7 +626,11 @@ class VaspInputSet(InputGenerator, abc.ABC):
             elif key == "KSPACING" and self.auto_kspacing:
                 # Default to metal if no prev calc available
                 bandgap = 0 if self.bandgap is None else self.bandgap
-                incar[key] = auto_kspacing(bandgap, self.bandgap_tol)
+                if new_kspacing := getattr(self, "kspacing_update", None):
+                    # allow custom KSPACING update
+                    incar[key] = new_kspacing
+                else:
+                    incar[key] = auto_kspacing(bandgap, self.bandgap_tol)
 
             else:
                 incar[key] = setting
@@ -1294,7 +1301,7 @@ class MPRelaxSet(VaspInputSet):
 
 @due.dcite(
     Doi("10.1021/acs.jpclett.0c02405"),
-    description="AccurAccurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation",
+    description="Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation",
 )
 @due.dcite(
     Doi("10.1103/PhysRevLett.115.036402"),
@@ -1352,7 +1359,7 @@ class MPScanRelaxSet(VaspInputSet):
 
         James W. Furness, Aaron D. Kaplan, Jinliang Ning, John P. Perdew, and Jianwei Sun.
         Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation.
-        The Journal of Physical Chemistry Letters 0, 11 DOI: 10.1021/acs.jpclett.0c02405
+        The Journal of Physical Chemistry Letters 11, 8208-8215 (2022) DOI: 10.1021/acs.jpclett.0c02405
     """
 
     bandgap: float | None = None
@@ -1375,6 +1382,149 @@ class MPScanRelaxSet(VaspInputSet):
                 self._config_dict["INCAR"].pop(k, None)
 
 
+@dataclass
+class MP24RelaxSet(VaspInputSet):
+    """
+    Materials Project relax set after a 2023-2024 benchmarking effort.
+
+    By default, this uses r2SCAN as the xc functional.
+    For discussion around this benchmarking effort, see:
+    https://github.com/materialsproject/foundation/pull/26
+
+    Citation info:
+    - r2SCAN:
+        James W. Furness, Aaron D. Kaplan, Jinliang Ning, John P. Perdew, and Jianwei Sun.
+        Accurate and Numerically Efficient r2SCAN Meta-Generalized Gradient Approximation.
+        J. Phys. Chem. Lett. 11, 8208-8215 (2022) DOI: 10.1021/acs.jpclett.0c02405
+    - r2SCAN-rVV10:
+        Jinliang Ning, Manish Kothakonda, James W. Furness, Aaron D. Kaplan, Sebastian Ehlert,
+            Jan Gerit Brandenburg, John P. Perdew, and Jianwei Sun.
+        Workhorse minimally empirical dispersion-corrected density functional with tests for
+            weakly bound systems: r2SCAN+rVV10.
+        Phys. Rev. B 106, 075422 (2022) DOI: 10.1103/PhysRevB.106.075422
+    - r2SCAN-D4:
+        Sebastian Ehlert, Uwe Huniar, Jinliang Ning, James W. Furness, Jianwei Sun,
+            Aaron D. Kaplan, John P. Perdew, and Jan Gerit Brandenburg.
+        r2SCAN-D4: Dispersion corrected meta-generalized gradient approximation for general chemical applications.
+        J. Chem. Phys. 154, 061101 (2021) DOI: 10.1063/5.0041008
+    - PBE:
+        John P. Perdew, Kieron Burke, and Matthias Ernzerhof,
+        Generalized Gradient Approximation Made Simple,
+        Phys. Rev. Lett. 77, 3865 (1996) DOI: 10.1103/PhysRevLett.77.3865
+    - PBEsol:
+        John P. Perdew, Adrienn Ruzsinszky, Gábor I. Csonka, Oleg A. Vydrov,
+            Gustavo E. Scuseria, Lucian A. Constantin, Xiaolan Zhou, and Kieron Burke.
+        Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces.
+        Phys. Rev. Lett. 100, 136406 (2009) DOI: 10.1103/PhysRevLett.100.136406
+    - PBE-D4 and PBEsol-D4:
+        Eike Caldeweyher,Jan-Michael Mewes, Sebastian Ehlert, and Stefan Grimme.
+        Extension and evaluation of the D4 London-dispersion model for periodic systems.
+        Phys. Chem. Chem. Phys. 22, 8499-8512 (2020) DOI: 10.1039/D0CP00502A
+    """
+
+    xc_functional: Literal["r2SCAN", "PBE", "PBEsol"] = "r2SCAN"
+    dispersion: Literal["rVV10", "D4"] | None = None
+    CONFIG = _load_yaml_config("MP24RelaxSet")
+    auto_ismear: bool = False
+    auto_kspacing: bool = True
+    inherit_incar: bool = False
+
+    def __post_init__(self) -> None:
+        super().__post_init__()
+
+        to_func = {
+            "r2SCAN": "R2SCAN",
+            "PBE": "PE",
+            "PBEsol": "PS",
+        }
+
+        config_updates: dict[str, Any] = {}
+        if self.xc_functional == "r2SCAN":
+            config_updates |= {"METAGGA": to_func[self.xc_functional]}
+            self._config_dict["INCAR"].pop("GGA", None)
+        elif self.xc_functional in ["PBE", "PBEsol"]:
+            config_updates |= {"GGA": to_func[self.xc_functional]}
+            self._config_dict["INCAR"].pop("METAGGA", None)
+        else:
+            raise ValueError(f"Unknown XC functional {self.xc_functional}!")
+
+        if self.dispersion == "rVV10":
+            if self.xc_functional == "r2SCAN":
+                config_updates |= {"BPARAM": 11.95, "CPARAM": 0.0093, "LUSE_VDW": True}
+            else:
+                raise ValueError(
+                    "Use of rVV10 with functionals other than r2 / SCAN is not currently supported in VASP."
+                )
+
+        elif self.dispersion == "D4":
+            d4_pars = {
+                "r2SCAN": {
+                    "S6": 1.0,
+                    "S8": 0.60187490,
+                    "A1": 0.51559235,
+                    "A2": 5.77342911,
+                },
+                "PBE": {
+                    "S6": 1.0,
+                    "S8": 0.95948085,
+                    "A1": 0.38574991,
+                    "A2": 4.80688534,
+                },
+                "PBEsol": {
+                    "S6": 1.0,
+                    "S8": 1.71885698,
+                    "A1": 0.47901421,
+                    "A2": 5.96771589,
+                },
+            }
+            config_updates |= {"IVDW": 13, **{f"VDW_{k}": v for k, v in d4_pars[self.xc_functional].items()}}
+
+        if len(config_updates) > 0:
+            self._config_dict["INCAR"].update(config_updates)
+
+    @staticmethod
+    def _sigmoid_interp(
+        bg, min_dk: float = 0.22, max_dk: float = 0.5, shape: float = 0.43, center: float = 4.15, fac: float = 12.0
+    ):
+        delta = shape * (bg - center)
+        sigmoid = delta / (1.0 + delta**fac) ** (1.0 / fac)
+        return 0.5 * (min_dk + max_dk + (max_dk - min_dk) * sigmoid)
+
+    def _multi_sigmoid_interp(
+        self,
+        bandgap: float,
+        dks: tuple[float, ...] = (0.22, 0.44, 0.5),
+        shape: tuple[float, ...] = (1.1, 2.5),
+        center: tuple[float, ...] = (2.35, 6.1),
+        fac: tuple[float, ...] = (8, 8),
+        bg_cut: tuple[float, ...] = (4.5,),
+    ):
+        if bandgap < self.bandgap_tol:
+            return dks[0]
+
+        min_bds = [self.bandgap_tol, *bg_cut]
+        max_bds = [*bg_cut, np.inf]
+
+        for icut, min_bd in enumerate(min_bds):
+            if min_bd <= bandgap < max_bds[icut]:
+                return self._sigmoid_interp(
+                    bandgap,
+                    min_dk=dks[icut],
+                    max_dk=dks[icut + 1],
+                    shape=shape[icut],
+                    center=center[icut],
+                    fac=fac[icut],
+                )
+
+        return None
+
+    @property
+    def kspacing_update(self):
+        if self.bandgap is None:
+            return 0.22
+        return self._multi_sigmoid_interp(self.bandgap)
+
+
 @dataclass
 class MPMetalRelaxSet(VaspInputSet):
     """
@@ -1571,7 +1721,8 @@ class MatPESStaticSet(VaspInputSet):
             )
 
         if self.xc_functional.upper() == "R2SCAN":
-            self._config_dict["INCAR"].update({"METAGGA": "R2SCAN", "ALGO": "ALL", "GGA": None})
+            self._config_dict["INCAR"].update({"METAGGA": "R2SCAN", "ALGO": "ALL"})
+            self._config_dict["INCAR"].pop("GGA", None)
         if self.xc_functional.upper().endswith("+U"):
             self._config_dict["INCAR"]["LDAU"] = True
 
@@ -1601,6 +1752,7 @@ class MPScanStaticSet(MPScanRelaxSet):
     def incar_updates(self) -> dict[str, Any]:
         """Updates to the INCAR config for this calculation type."""
         updates: dict[str, Any] = {
+            "ALGO": "Fast",
             "LREAL": False,
             "NSW": 0,
             "LORBIT": 11,
@@ -1626,6 +1778,49 @@ class MPScanStaticSet(MPScanRelaxSet):
         return updates
 
 
+@dataclass
+class MP24StaticSet(MP24RelaxSet):
+    """Create input files for a static calculation using MP24 parameters.
+
+    For class information, refer to `MP24RelaxSet`.
+    If you use this set, please consider citing the appropriate papers in `MP24RelaxSet`.
+
+    Args:
+        structure (Structure): Structure from previous run.
+        bandgap (float): Bandgap of the structure in eV. The bandgap is used to
+            compute the appropriate k-point density and determine the smearing settings.
+        lepsilon (bool): Whether to add static dielectric calculation
+        lcalcpol (bool): Whether to turn on evaluation of the Berry phase approximations
+            for electronic polarization.
+        **kwargs: Keywords supported by MP24RelaxSet.
+    """
+
+    lepsilon: bool = False
+    lcalcpol: bool = False
+    inherit_incar: bool = False
+    auto_kspacing: bool = True
+
+    @property
+    def incar_updates(self) -> dict[str, Any]:
+        """Updates to the INCAR config for this calculation type."""
+        updates: dict[str, Any] = {
+            "NSW": 0,
+            "LORBIT": 11,
+            "ISMEAR": -5,
+        }
+
+        if self.lepsilon:
+            # LPEAD=T: numerical evaluation of overlap integral prevents
+            # LRF_COMMUTATOR errors and can lead to better expt. agreement
+            # but produces slightly different results
+            updates |= {"IBRION": 8, "LEPSILON": True, "LPEAD": True, "NSW": 1, "NPAR": None}
+
+        if self.lcalcpol:
+            updates["LCALCPOL"] = True
+
+        return updates
+
+
 @dataclass
 class MPHSEBSSet(VaspInputSet):
     """Implementation of a VaspInputSet for HSE band structure computations.


Test patch:
diff --git a/tests/analysis/test_pourbaix_diagram.py b/tests/analysis/test_pourbaix_diagram.py
index 2206649925..30150f32aa 100644
--- a/tests/analysis/test_pourbaix_diagram.py
+++ b/tests/analysis/test_pourbaix_diagram.py
@@ -8,7 +8,14 @@ import numpy as np
 from monty.serialization import dumpfn, loadfn
 from pytest import approx
 
-from pymatgen.analysis.pourbaix_diagram import IonEntry, MultiEntry, PourbaixDiagram, PourbaixEntry, PourbaixPlotter
+from pymatgen.analysis.pourbaix_diagram import (
+    IonEntry,
+    MultiEntry,
+    PourbaixDiagram,
+    PourbaixEntry,
+    PourbaixPlotter,
+    ion_or_solid_comp_object,
+)
 from pymatgen.core.composition import Composition
 from pymatgen.core.ion import Ion
 from pymatgen.entries.computed_entries import ComputedEntry
@@ -315,3 +322,34 @@ class TestPourbaixPlotter(TestCase):
         binary_plotter = PourbaixPlotter(pd_binary)
         ax = binary_plotter.plot_entry_stability(self.test_data["Ag-Te"][53])
         assert isinstance(ax, plt.Axes)
+
+
+class TestIonOrSolidCompObject:
+    def test_ion(self):
+        # Test cations
+        assert ion_or_solid_comp_object("Li+").charge == 1
+        assert ion_or_solid_comp_object("Li[+]").charge == 1
+        assert ion_or_solid_comp_object("Ca[2+]").charge == 2
+        assert ion_or_solid_comp_object("Ca[+2]").charge == 2
+        assert ion_or_solid_comp_object("Ca++").charge == 2
+        assert ion_or_solid_comp_object("Ca[++]").charge == 2
+        assert ion_or_solid_comp_object("Ca2+").charge == 1
+        assert ion_or_solid_comp_object("C2O4-2").charge == -2
+
+        # Test anions
+        assert ion_or_solid_comp_object("Cl-").charge == -1
+        assert ion_or_solid_comp_object("Cl[-]").charge == -1
+        assert ion_or_solid_comp_object("SO4[-2]").charge == -2
+        assert ion_or_solid_comp_object("SO4-2").charge == -2
+        assert ion_or_solid_comp_object("SO42-").charge == -1
+        assert ion_or_solid_comp_object("SO4--").charge == -2
+        assert ion_or_solid_comp_object("SO4[--]").charge == -2
+        assert ion_or_solid_comp_object("N3-").charge == -1
+
+    def test_solid(self):
+        # Test end with "(s)"
+        assert ion_or_solid_comp_object("Fe2O3(s)") == Composition("Fe2O3")
+        assert ion_or_solid_comp_object("Fe(s)") == Composition("Fe1")
+
+        # Test end without "(s)"
+        assert type(ion_or_solid_comp_object("Na2O")) is Composition
diff --git a/tests/io/exciting/test_inputs.py b/tests/io/exciting/test_inputs.py
index ccbac25394..0f119448c0 100644
--- a/tests/io/exciting/test_inputs.py
+++ b/tests/io/exciting/test_inputs.py
@@ -1,5 +1,7 @@
 from __future__ import annotations
 
+import math
+import re
 from xml.etree import ElementTree as ET
 
 from numpy.testing import assert_allclose
@@ -165,4 +167,12 @@ class TestExcitingInput(PymatgenTest):
         root = tree.getroot()
         ref_str = ET.tostring(root, encoding="unicode")
 
-        assert ref_str.strip() == test_str.strip()
+        ref_list = ref_str.strip().split()
+        test_list = test_str.strip().split()
+
+        # "scale" is float, direct compare might give surprising results
+        ref_scale = float(re.search(r'scale="([-+]?\d*\.\d+|\d+)"', ref_list.pop(7))[1])
+        test_scale = float(re.search(r'scale="([-+]?\d*\.\d+|\d+)"', test_list.pop(7))[1])
+
+        assert ref_list == test_list
+        assert math.isclose(ref_scale, test_scale)
diff --git a/tests/io/test_ase.py b/tests/io/test_ase.py
index 4459b4ea5f..4328d0f695 100644
--- a/tests/io/test_ase.py
+++ b/tests/io/test_ase.py
@@ -91,6 +91,18 @@ def test_get_atoms_from_structure_dyn():
     STRUCTURE.add_site_property("selective_dynamics", [[False] * 3] * len(STRUCTURE))
     atoms = AseAtomsAdaptor.get_atoms(STRUCTURE)
     assert atoms.constraints[0].get_indices().tolist() == [atom.index for atom in atoms]
+    STRUCTURE.add_site_property("selective_dynamics", [[True] * 3] * len(STRUCTURE))
+    atoms = AseAtomsAdaptor.get_atoms(STRUCTURE)
+    assert len(atoms.constraints) == 0
+    rng = np.random.default_rng(seed=1234)
+    sel_dyn = [[rng.random() < 0.5, rng.random() < 0.5, rng.random() < 0.5] for _ in STRUCTURE]
+    STRUCTURE.add_site_property("selective_dynamics", sel_dyn)
+    atoms = AseAtomsAdaptor.get_atoms(STRUCTURE)
+    for c in atoms.constraints:
+        # print(c)
+        assert isinstance(c, ase.constraints.FixAtoms | ase.constraints.FixCartesian)
+        ase_mask = c.mask if isinstance(c, ase.constraints.FixCartesian) else [True, True, True]
+        assert len(c.index) == len([mask for mask in sel_dyn if np.array_equal(mask, ~np.array(ase_mask))])
 
 
 def test_get_atoms_from_molecule():
@@ -178,8 +190,10 @@ def test_get_structure_mag():
     "select_dyn",
     [
         [True, True, True],
+        [True, False, True],
         [False, False, False],
         np.array([True, True, True]),
+        np.array([True, False, True]),
         np.array([False, False, False]),
     ],
 )
diff --git a/tests/io/vasp/test_sets.py b/tests/io/vasp/test_sets.py
index 4509e1b95f..cef111d7d7 100644
--- a/tests/io/vasp/test_sets.py
+++ b/tests/io/vasp/test_sets.py
@@ -29,6 +29,8 @@ from pymatgen.io.vasp.sets import (
     MITMDSet,
     MITNEBSet,
     MITRelaxSet,
+    MP24RelaxSet,
+    MP24StaticSet,
     MPAbsorptionSet,
     MPHSEBSSet,
     MPHSERelaxSet,
@@ -116,13 +118,14 @@ class TestSetChangeCheck(PymatgenTest):
             "MPHSERelaxSet.yaml": "1779cb6a6af43ad54a12aec22882b9b8aa3469b764e29ac4ab486960d067b811",
             "VASPIncarBase.yaml": "8c1ce90d6697e45b650e1881e2b3d82a733dba17fb1bd73747a38261ec65a4c4",
             "MPSCANRelaxSet.yaml": "ad652ea740d06f9edd979494f31e25074b82b9fffdaaf7eff2ae5541fb0e6288",
-            "PBE64Base.yaml": "3434c918c17706feae397d0852f2224e771db94d7e4c988039e8658e66d87494",
+            "PBE64Base.yaml": "40e7e42159f59543b17f512666916001045f7644f422ccc45b8466d6a1cf0c48",
             "MPRelaxSet.yaml": "c9b0a519588fb3709509a9f9964632692584905e2961a0fe2e5f657561913083",
             "MITRelaxSet.yaml": "0b4bec619fa860dac648584853c3b3d5407e4148a85d0e95024fbd1dc315669d",
             "vdW_parameters.yaml": "7d2599a855533865335a313c043b6f89e03fc2633c88b6bc721723d94cc862bd",
             "MatPESStaticSet.yaml": "4ec60ad4bbbb9a756f1b3fea8ca4eab8fc767d8f6a67332e7af3908c910fd7c5",
             "MPAbsorptionSet.yaml": "e49cd0ab87864f1c244e9b5ceb4703243116ec1fbb8958a374ddff07f7a5625c",
             "PBE54Base.yaml": "cdffe123eca8b19354554b60a7f8de9b8776caac9e1da2bd2a0516b7bfac8634",
+            "MP24RelaxSet.yaml": "35a5d4456f01d644cf41218725c5e0896c59e1658045ecd1544579cbb1ed7b85",
         }
 
         for input_set, hash_str in hashes.items():
@@ -2285,3 +2288,98 @@ def test_dict_set_alias():
     ):
         DictSet()
         assert isinstance(DictSet(), VaspInputSet)
+
+
+class TestMP24Sets(PymatgenTest):
+    def setUp(self):
+        self.relax_set = MP24RelaxSet
+        self.static_set = MP24StaticSet
+
+        filepath = f"{VASP_IN_DIR}/POSCAR"
+        self.structure = Structure.from_file(filepath)
+
+    @staticmethod
+    def matches_ref(test_val, ref_val) -> bool:
+        if isinstance(ref_val, float):
+            return test_val == approx(ref_val)
+        return test_val == ref_val
+
+    def test_kspacing(self):
+        bandgaps = [0.0, 0.1, 0.5, 1.0, 2.0, 3, 5, 10, 1e4]
+        expected_kspacing = [
+            0.22,
+            0.22000976174867864,
+            0.2200466614246148,
+            0.22056799311325073,
+            0.2876525546497567,
+            0.40800309817134106,
+            0.44000114629141485,
+            0.4999999999540808,
+            0.5,
+        ]
+        for i, bandgap in enumerate(bandgaps):
+            assert self.relax_set()._multi_sigmoid_interp(bandgap) == approx(expected_kspacing[i])
+
+    def test_default(self):
+        vis = self.relax_set(structure=self.structure)
+        expected_incar_relax = {
+            "ALGO": "Normal",
+            "EDIFF": 1.0e-05,
+            "EDIFFG": -0.02,
+            "ENAUG": 1360,
+            "ENCUT": 680,
+            "GGA_COMPAT": False,
+            "KSPACING": 0.22,
+            "ISMEAR": 0,
+            "SIGMA": 0.05,
+            "METAGGA": "R2scan",
+            "LMAXMIX": 6,
+            "LREAL": False,
+        }
+
+        assert all(self.matches_ref(vis.incar[k], v) for k, v in expected_incar_relax.items())
+
+        assert self.relax_set(self.structure, xc_functional="r2SCAN")._config_dict == vis._config_dict
+        assert vis.inherit_incar is False
+        assert vis.dispersion is None
+        assert vis.potcar_functional == "PBE_64"
+        assert vis.potcar_symbols == ["Fe_pv", "P", "O"]
+        assert vis.kpoints is None
+
+        vis = self.static_set(structure=self.structure, dispersion="rVV10")
+
+        expected_incar_static = {
+            "ISMEAR": -5,
+            "NSW": 0,
+            "LORBIT": 11,
+            "METAGGA": "R2scan",
+            "LUSE_VDW": True,
+            "BPARAM": 11.95,
+            "CPARAM": 0.0093,
+        }
+
+        assert all(self.matches_ref(vis.incar[k], v) for k, v in expected_incar_static.items())
+        assert (
+            self.static_set(self.structure, xc_functional="r2SCAN", dispersion="rVV10")._config_dict == vis._config_dict
+        )
+
+    def test_non_default_xc_func(self):
+        for xc_functional, vasp_name in {"PBE": "Pe", "PBEsol": "Ps"}.items():
+            vis = self.relax_set(structure=self.structure, xc_functional=xc_functional)
+            assert vis.incar.get("METAGGA") is None
+            assert vis.incar["GGA"] == vasp_name
+
+            vis = self.static_set(structure=self.structure, xc_functional=xc_functional, dispersion="D4")
+            assert vis.incar.get("METAGGA") is None
+            print(self.relax_set(structure=self.structure, xc_functional=xc_functional).incar)
+            assert vis.incar["GGA"] == vasp_name
+            assert all(
+                isinstance(vis.incar[k], float | int)
+                for k in (
+                    "IVDW",
+                    "VDW_A1",
+                    "VDW_A2",
+                    "VDW_S6",
+                    "VDW_S8",
+                )
+            )