Add formal_chempots option to ChemicalPotentialDiagram to plot the formal chemical potentials rather than the DFT energies (materialsproject#2916)

* Add `formal_chempots` option to `ChemicalPotentialDiagram` to plot the formal chemical potentials rather than the DFT energies

* Add tests for `formal_chempots` option in `ChemicalPotentialDiagram` to all applicable tests in `test_chempot_diagram.py`

* remove useless dict.copy() in _renormalize_entry(), fix some doc strings and typos

* replace assertArrayAlmostEqual() with pytest.approx() in test_chempot_diagram.py

* use dot access for plot attributes

---------

Co-authored-by: Janosh Riebesell <[email protected]>

Author: 2 people

pymatgen/analysis/chempot_diagram.py

@@ -65,24 +65,42 @@ class ChemicalPotentialDiagram(MSONable):
     def __init__(
         self,
         entries: list[PDEntry],
-        limits: dict[Element, float] | None = None,
+        limits: dict[Element, tuple[float, float]] | None = None,
         default_min_limit: float = -50.0,
-    ):
+        formal_chempots: bool = True,
+    ) -> None:
         """
         Args:
-            entries: List of PDEntry-like objects containing a composition and
+            entries (list[PDEntry]): PDEntry-like objects containing a composition and
                 energy. Must contain elemental references and be suitable for typical
                 phase diagram construction. Entries must be within a chemical system
                 of with 2+ elements.
-            limits: Bounds of elemental chemical potentials (min, max), which are
-                used to construct the border hyperplanes used in the
-                HalfSpaceIntersection algorithm; these constrain the space over which the
-                domains are calculated and also determine the size of the plotted
-                diagram. Any elemental limits not specified are covered in the
-                default_min_limit argument. e.g., {Element("Li"): [-12.0, 0.0], ...}
+            limits (dict[Element, float] | None): Bounds of elemental chemical potentials (min, max),
+                which are used to construct the border hyperplanes used in the HalfSpaceIntersection
+                algorithm; these constrain the space over which the domains are calculated and also
+                determine the size of the plotted diagram. Any elemental limits not specified are
+                covered in the default_min_limit argument. e.g., {Element("Li"): [-12.0, 0.0], ...}
             default_min_limit (float): Default minimum chemical potential limit (i.e.,
                 lower bound) for unspecified elements within the "limits" argument.
+            formal_chempots (bool): Whether to plot the formal ('reference') chemical potentials
+                (i.e. μ_X - μ_X^0) or the absolute DFT reference energies (i.e. μ_X(DFT)).
+                Default is True (i.e. plot formal chemical potentials).
         """
+        entries = sorted(entries, key=lambda e: e.composition.reduced_composition)
+        _min_entries, _el_refs = self._get_min_entries_and_el_refs(entries)
+
+        if formal_chempots:
+            # renormalize entry energies to be relative to the elemental references
+            renormalized_entries = []
+            for entry in entries:
+                comp_dict = entry.composition.as_dict()
+                renormalization_energy = sum(
+                    [comp_dict[el] * _el_refs[Element(el)].energy_per_atom for el in comp_dict]
+                )
+                renormalized_entries.append(_renormalize_entry(entry, renormalization_energy / sum(comp_dict.values())))
+
+            entries = renormalized_entries
+
         self.entries = sorted(entries, key=lambda e: e.composition.reduced_composition)
         self.limits = limits
         self.default_min_limit = default_min_limit
@@ -622,7 +640,7 @@ def __repr__(self):
 
 def simple_pca(data: np.ndarray, k: int = 2) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     """
-    A barebones implementation of principal component analysis (PCA) used in the
+    A bare-bones implementation of principal component analysis (PCA) used in the
     ChemicalPotentialDiagram class for plotting.
 
     Args:
@@ -645,15 +663,15 @@ def simple_pca(data: np.ndarray, k: int = 2) -> tuple[np.ndarray, np.ndarray, np
 
 def get_centroid_2d(vertices: np.ndarray) -> np.ndarray:
     """
-    A barebones implementation of the formula for calculating the centroid of a 2D
+    A bare-bones implementation of the formula for calculating the centroid of a 2D
     polygon. Useful for calculating the location of an annotation on a chemical
     potential domain within a 3D chemical potential diagram.
 
-    **NOTE**: vertices must be ordered circumfrentially!
+    **NOTE**: vertices must be ordered circumferentially!
 
     Args:
         vertices: array of 2-d coordinates corresponding to a polygon, ordered
-            circumfrentially
+            circumferentially
 
     Returns:
         Array giving 2-d centroid coordinates
@@ -690,7 +708,7 @@ def get_2d_orthonormal_vector(line_pts: np.ndarray) -> np.ndarray:
             coordinates of a line
 
     Returns:
-
+        np.ndarray: A length-2 vector that is orthonormal to the line.
     """
     x = line_pts[:, 0]
     y = line_pts[:, 1]
@@ -703,3 +721,15 @@ def get_2d_orthonormal_vector(line_pts: np.ndarray) -> np.ndarray:
     vec = np.array([np.sin(theta), np.cos(theta)])
 
     return vec
+
+
+def _renormalize_entry(entry: PDEntry, renormalization_energy_per_atom: float) -> PDEntry:
+    """
+    Regenerate the input entry with an energy per atom decreased by renormalization_energy_per_atom
+    """
+    renormalized_entry_dict = entry.as_dict()
+    renormalized_entry_dict["energy"] = entry.energy - renormalization_energy_per_atom * sum(
+        entry.composition.values()
+    )  # entry.energy includes MP corrections as desired
+    renormalized_entry = PDEntry.from_dict(renormalized_entry_dict)
+    return renormalized_entry

pymatgen/analysis/tests/test_chempot_diagram.py

@@ -5,6 +5,7 @@
 from pathlib import Path
 
 import numpy as np
+import pytest
 from plotly.graph_objects import Figure
 
 from pymatgen.analysis.chempot_diagram import (
@@ -23,15 +24,18 @@
 class ChemicalPotentialDiagramTest(PymatgenTest):
     def setUp(self):
         self.entries = EntrySet.from_csv(str(module_dir / "pdentries_test.csv"))
-        self.cpd_ternary = ChemicalPotentialDiagram(entries=self.entries, default_min_limit=-25)
+        self.cpd_ternary = ChemicalPotentialDiagram(entries=self.entries, default_min_limit=-25, formal_chempots=False)
+        self.cpd_ternary_formal = ChemicalPotentialDiagram(
+            entries=self.entries, default_min_limit=-25, formal_chempots=True
+        )
         elements = [Element("Fe"), Element("O")]
         binary_entries = list(
             filter(
                 lambda e: set(e.composition.elements).issubset(elements),
                 self.entries,
             )
         )
-        self.cpd_binary = ChemicalPotentialDiagram(entries=binary_entries, default_min_limit=-25)
+        self.cpd_binary = ChemicalPotentialDiagram(entries=binary_entries, default_min_limit=-25, formal_chempots=False)
         warnings.simplefilter("ignore")
 
     def tearDown(self):
@@ -40,6 +44,7 @@ def tearDown(self):
     def test_dim(self):
         assert self.cpd_binary.dim == 2
         assert self.cpd_ternary.dim == 3
+        assert self.cpd_ternary_formal.dim == 3
 
     def test_el_refs(self):
         el_refs = {elem: entry.energy for elem, entry in self.cpd_ternary.el_refs.items()}
@@ -48,17 +53,26 @@ def test_el_refs(self):
         energies = [-1.91301487, -6.5961471, -25.54966885]
         correct_el_refs = dict(zip(elems, energies))
 
-        self.assertDictsAlmostEqual(el_refs, correct_el_refs)
+        assert el_refs == pytest.approx(correct_el_refs)
+
+    def test_el_refs_formal(self):
+        el_refs = {elem: entry.energy for elem, entry in self.cpd_ternary_formal.el_refs.items()}
+        elems = [Element("Li"), Element("Fe"), Element("O")]
+        energies = [0, 0, 0]
+        correct_el_refs = dict(zip(elems, energies))
+        assert el_refs == pytest.approx(correct_el_refs)
 
     def test_border_hyperplanes(self):
         desired = np.array(
             [[-1, 0, 0, -25], [1, 0, 0, 0], [0, -1, 0, -25], [0, 1, 0, 0], [0, 0, -1, -25], [0, 0, 1, 0]]
         )
-        self.assertArrayAlmostEqual(self.cpd_ternary.border_hyperplanes, desired)
+        assert self.cpd_ternary.border_hyperplanes == pytest.approx(desired)
+        assert self.cpd_ternary_formal.border_hyperplanes == pytest.approx(desired)
 
     def test_lims(self):
         desired_lims = np.array([[-25, 0], [-25, 0], [-25, 0]])
-        self.assertArrayAlmostEqual(self.cpd_ternary.lims, desired_lims)
+        assert self.cpd_ternary.lims == pytest.approx(desired_lims)
+        assert self.cpd_ternary_formal.lims == pytest.approx(desired_lims)
 
     def test_pca(self):
         points_3d = np.array(
@@ -80,7 +94,7 @@ def test_pca(self):
 
         points_2d, _, _ = simple_pca(points_3d, k=2)
 
-        self.assertArrayAlmostEqual(points_2d, points_2d_desired)
+        assert points_2d == pytest.approx(points_2d_desired)
 
     def test_centroid(self):
         vertices = np.array(
@@ -97,7 +111,7 @@ def test_centroid(self):
         centroid = get_centroid_2d(vertices)
         centroid_desired = np.array([-0.00069433, -0.00886174])
 
-        self.assertArrayAlmostEqual(centroid, centroid_desired)
+        assert centroid == pytest.approx(centroid_desired, abs=1e-6)
 
     def test_get_2d_orthonormal_vector(self):
         pts_1 = np.array([[1, 1], [2, 2]])
@@ -109,18 +123,25 @@ def test_get_2d_orthonormal_vector(self):
         vec_1_desired = np.array([0.70710678, 0.70710678])
         vec_2_desired = np.array([0.98386991, 0.17888544])
 
-        self.assertArrayAlmostEqual(vec_1, vec_1_desired)
-        self.assertArrayAlmostEqual(vec_2, vec_2_desired)
+        assert vec_1 == pytest.approx(vec_1_desired)
+        assert vec_2 == pytest.approx(vec_2_desired)
 
     def test_get_plot(self):
         fig_2d = self.cpd_binary.get_plot()
         fig_3d = self.cpd_ternary.get_plot()
+        fig_3d_formal = self.cpd_ternary_formal.get_plot()
 
         assert isinstance(fig_2d, Figure)
-        assert fig_2d["data"][0]["type"] == "scatter"
+        assert fig_2d.data[0].type == "scatter"
 
         assert isinstance(fig_3d, Figure)
-        assert fig_3d["data"][0]["type"] == "scatter3d"
+        assert fig_3d.data[0].type == "scatter3d"
+
+        assert isinstance(fig_3d_formal, Figure)
+        assert fig_3d_formal.data[0].type == "scatter3d"
+        assert fig_3d_formal.data[0].mode == "lines"
+        assert fig_3d_formal.layout.plot_bgcolor == "rgba(0,0,0,0)"
+        assert fig_3d_formal.layout.scene.annotations[0].text == "FeO"
 
     def test_domains(self):
         correct_domains = {
@@ -229,7 +250,114 @@ def test_domains(self):
             d = self.cpd_ternary.domains[formula]
             d = d.round(6)  # to get rid of numerical errors from qhull
             actual_domain_sorted = d[np.lexsort((d[:, 2], d[:, 1], d[:, 0]))]
-            self.assertArrayAlmostEqual(actual_domain_sorted, domain)
+            assert actual_domain_sorted == pytest.approx(domain)
+
+        formal_domains = {
+            "FeO": np.array(
+                [
+                    [-2.50000000e01, 3.55271368e-15, -2.85707600e00],
+                    [-2.01860032e00, 3.55271368e-15, -2.85707600e00],
+                    [-2.50000000e01, -1.45446765e-01, -2.71162923e00],
+                    [-2.16404709e00, -1.45446765e-01, -2.71162923e00],
+                ]
+            ),
+            "Fe2O3": np.array(
+                [
+                    [-25.0, -4.14354109, 0.0],
+                    [-3.637187, -4.14354108, 0.0],
+                    [-3.49325969, -3.85568646, -0.19190308],
+                    [-25.0, -0.70024301, -2.29553205],
+                    [-2.44144521, -0.70024301, -2.29553205],
+                ]
+            ),
+            "Fe3O4": np.array(
+                [
+                    [-25.0, -0.70024301, -2.29553205],
+                    [-25.0, -0.14544676, -2.71162923],
+                    [-2.44144521, -0.70024301, -2.29553205],
+                    [-2.16404709, -0.14544676, -2.71162923],
+                ]
+            ),
+            "LiFeO2": np.array(
+                [
+                    [-3.49325969e00, -3.85568646e00, -1.91903083e-01],
+                    [-2.01860032e00, 3.55271368e-15, -2.85707600e00],
+                    [-2.44144521e00, -7.00243005e-01, -2.29553205e00],
+                    [-2.16404709e00, -1.45446765e-01, -2.71162923e00],
+                    [-1.71198739e00, 3.55271368e-15, -3.01038246e00],
+                    [-2.74919447e00, -3.11162124e00, -9.35968300e-01],
+                ]
+            ),
+            "Li2O": np.array(
+                [
+                    [0.00000000e00, -2.50000000e01, -6.22930387e00],
+                    [-2.69949567e00, -2.50000000e01, -8.30312528e-01],
+                    [3.55271368e-15, 3.55271368e-15, -6.22930387e00],
+                    [-1.43858289e00, 3.55271368e-15, -3.35213809e00],
+                    [-2.69949567e00, -3.78273835e00, -8.30312528e-01],
+                ]
+            ),
+            "Li2O2": np.array(
+                [
+                    [-3.52980820e00, -2.50000000e01, 0.00000000e00],
+                    [-2.69949567e00, -2.50000000e01, -8.30312528e-01],
+                    [-3.52980820e00, -4.35829869e00, 3.55271368e-15],
+                    [-2.69949567e00, -3.78273835e00, -8.30312528e-01],
+                    [-2.82687176e00, -3.65536226e00, -7.02936437e-01],
+                ]
+            ),
+            "Li2FeO3": np.array(
+                [
+                    [-3.52980820e00, -4.35829869e00, 3.55271368e-15],
+                    [-3.63718700e00, -4.14354108e00, 0.00000000e00],
+                    [-3.49325969e00, -3.85568646e00, -1.91903083e-01],
+                    [-2.74919447e00, -3.11162124e00, -9.35968300e-01],
+                    [-2.82687176e00, -3.65536226e00, -7.02936437e-01],
+                ]
+            ),
+            "Li5FeO4": np.array(
+                [
+                    [-1.43858289e00, 3.55271368e-15, -3.35213809e00],
+                    [-1.71198739e00, 3.55271368e-15, -3.01038246e00],
+                    [-2.74919447e00, -3.11162124e00, -9.35968300e-01],
+                    [-2.69949567e00, -3.78273835e00, -8.30312528e-01],
+                    [-2.82687176e00, -3.65536226e00, -7.02936437e-01],
+                ]
+            ),
+            "O2": np.array(
+                [
+                    [-2.50000000e01, -2.50000000e01, 3.55271368e-15],
+                    [-3.52980820e00, -2.50000000e01, 0.00000000e00],
+                    [-2.50000000e01, -4.14354109e00, 0.00000000e00],
+                    [-3.52980820e00, -4.35829869e00, 3.55271368e-15],
+                    [-3.63718700e00, -4.14354108e00, 0.00000000e00],
+                ]
+            ),
+            "Fe": np.array(
+                [
+                    [0.00000000e00, 0.00000000e00, -2.50000000e01],
+                    [-2.50000000e01, 0.00000000e00, -2.50000000e01],
+                    [3.55271368e-15, 3.55271368e-15, -6.22930387e00],
+                    [-2.50000000e01, 3.55271368e-15, -2.85707600e00],
+                    [-2.01860032e00, 3.55271368e-15, -2.85707600e00],
+                    [-1.43858289e00, 3.55271368e-15, -3.35213809e00],
+                    [-1.71198739e00, 3.55271368e-15, -3.01038246e00],
+                ]
+            ),
+            "Li": np.array(
+                [
+                    [3.55271368e-15, -2.50000000e01, -2.50000000e01],
+                    [0.00000000e00, -2.50000000e01, -6.22930387e00],
+                    [0.00000000e00, 0.00000000e00, -2.50000000e01],
+                    [3.55271368e-15, 3.55271368e-15, -6.22930387e00],
+                ]
+            ),
+        }
+
+        for formula, domain in formal_domains.items():
+            d = self.cpd_ternary_formal.domains[formula]
+            d = d.round(6)  # to get rid of numerical errors from qhull
+            assert d == pytest.approx(domain, abs=1e-5)
 
 
 if __name__ == "__main__":