Add PhaseDiagram method get_reference_energy

pymatgen/analysis/phase_diagram.py

@@ -543,15 +543,27 @@ def _get_stable_entries_in_space(self, space) -> list[Entry]:
         """
         return [e for e, s in zip(self._stable_entries, self._stable_spaces) if space.issuperset(s)]
 
-    def get_reference_energy_per_atom(self, comp: Composition) -> float:
+    def get_reference_energy(self, comp: Composition) -> float:
+        """Sum of elemental reference energies over all elements in a composition.
+
+        Args:
+            comp (Composition): Input composition.
+
+        Returns:
+            float: Reference energy
         """
+        return sum(comp[el] * self.el_refs[el].energy_per_atom for el in comp.elements)
+
+    def get_reference_energy_per_atom(self, comp: Composition) -> float:
+        """Sum of elemental reference energies over all elements in a composition.
+
         Args:
             comp (Composition): Input composition.
 
         Returns:
-            Reference energy of the terminal species at a given composition.
+            float: Reference energy per atom
         """
-        return sum(comp[el] * self.el_refs[el].energy_per_atom for el in comp.elements) / comp.num_atoms
+        return self.get_reference_energy(comp) / comp.num_atoms
 
     def get_form_energy(self, entry: PDEntry) -> float:
         """
@@ -565,7 +577,7 @@ def get_form_energy(self, entry: PDEntry) -> float:
             float: Formation energy from the elemental references.
         """
         comp = entry.composition
-        return entry.energy - sum(comp[el] * self.el_refs[el].energy_per_atom for el in comp.elements)
+        return entry.energy - self.get_reference_energy(comp)
 
     def get_form_energy_per_atom(self, entry: PDEntry) -> float:
         """

tests/analysis/test_phase_diagram.py

@@ -210,39 +210,53 @@ def test_ordering(self):
 
     def test_stable_entries(self):
         stable_formulas = [ent.composition.reduced_formula for ent in self.pd.stable_entries]
-        expected_stable = [
-            "Fe2O3",
-            "Li5FeO4",
-            "LiFeO2",
-            "Fe3O4",
-            "Li",
-            "Fe",
-            "Li2O",
-            "O2",
-            "FeO",
-        ]
+        expected_stable = "Fe2O3 Li5FeO4 LiFeO2 Fe3O4 Li Fe Li2O O2 FeO".split()
         for formula in expected_stable:
-            assert formula in stable_formulas, formula + " not in stable entries!"
+            assert formula in stable_formulas, f"{formula} not in stable entries!"
 
-    def test_get_formation_energy(self):
-        stable_formation_energies = {
-            ent.composition.reduced_formula: self.pd.get_form_energy(ent) for ent in self.pd.stable_entries
-        }
+    def test_get_form_energy(self):
         expected_formation_energies = {
-            "Li5FeO4": -164.8117344866667,
-            "Li2O2": -14.119232793333332,
-            "Fe2O3": -16.574164339999996,
-            "FeO": -5.7141519966666685,
+            "Li5FeO4": -164.8117344,
+            "Li2O2": -14.119232793,
+            "Fe2O3": -16.57416433,
+            "FeO": -5.71415199,
             "Li": 0.0,
-            "LiFeO2": -7.732752316666666,
-            "Li2O": -6.229303868333332,
+            "LiFeO2": -7.73275231,
+            "Li2O": -6.229303868,
             "Fe": 0.0,
-            "Fe3O4": -22.565714456666683,
-            "Li2FeO3": -45.67166036000002,
+            "Fe3O4": -22.56571445,
+            "Li2FeO3": -45.67166036,
             "O2": 0.0,
         }
-        for formula, energy in expected_formation_energies.items():
-            assert energy == approx(stable_formation_energies[formula])
+
+        for entry in self.pd.stable_entries:
+            formula = entry.composition.reduced_formula
+            expected = expected_formation_energies[formula]
+            n_atoms = entry.composition.num_atoms
+            # test get_form_energy
+            assert self.pd.get_form_energy(entry) == approx(expected), formula
+            # test get_form_energy_per_atom
+            assert self.pd.get_form_energy_per_atom(entry) == approx(expected / n_atoms), formula
+
+    def test_get_reference_energy(self):
+        expected_ref_energies = {
+            "Li5FeO4": -265.5546721,
+            "Li2O2": -24.685172046,
+            "Fe2O3": -51.93425725,
+            "FeO": -21.70885048,
+            "Li": -1.91301487,
+            "LiFeO2": -17.02571825,
+            "Li2O": -8.084307881,
+            "Fe": -6.5961471,
+            "Fe3O4": -73.6431077,
+            "Li2FeO3": -92.78804506,
+            "O2": -25.54966885,
+        }
+        for entry in self.pd.stable_entries:
+            formula = entry.composition.reduced_formula
+            actual = self.pd.get_reference_energy(entry.composition)
+            expected = expected_ref_energies[formula]
+            assert actual == approx(expected), formula
 
     def test_all_entries_hulldata(self):
         assert len(self.pd.all_entries_hulldata) == 490