diff --git a/doped/chemical_potentials.py b/doped/chemical_potentials.py
index 88789b6a..6fb2f778 100644
--- a/doped/chemical_potentials.py
+++ b/doped/chemical_potentials.py
@@ -1509,6 +1509,117 @@ def cplap_input(self, dependent_variable=None, filename="input.dat"):
                     print(int(v), k, end=" ")
                 print(f"{i['formation_energy']}  # num_atoms, element, formation_energy")
 
+    def to_LaTeX_table(self, columns=1, sort_by_energy=False, prune_polymorphs=True):
+        """
+        A very simple function to print out the competing phase formation
+        energies in a LaTeX table format. Needs the mhchem package to work and
+        does *not* use the booktabs package -- change hline to toprule, midrule
+        and bottomrule if you want to use booktabs style.
+
+        Args:
+            columns (int):
+                Number of columns to print out the table in; either 1 or 2
+            sort_by_energy (bool):
+                If True, sorts the table by formation energy (highest to lowest).
+                Default sorting by formula.
+            prune_polymorphs (bool):
+                Whether to only print out the lowest energy polymorphs for each composition.
+                Default is True.
+
+        Returns:
+            str: LaTeX table string
+        """
+        if columns not in [1, 2]:
+            raise ValueError("columns must be either 1 or 2")
+        # done in the pyscfermi report style
+        data = copy.deepcopy(self.data)
+
+        if any("kpoints" not in item for item in data):
+            raise (
+                ValueError(
+                    "kpoints need to be present in data; run CompetingPhasesAnalyzer.from_vaspruns "
+                    "instead of from_csv"
+                )
+            )
+
+        if prune_polymorphs:  # only keep the lowest energy polymorphs
+            formation_energy_df = _calculate_formation_energies(data, self.elemental_energies)
+            indices = formation_energy_df.groupby("formula")["energy_per_atom"].idxmin()
+            pruned_df = formation_energy_df.loc[indices]
+            data = pruned_df.to_dict(orient="records")
+
+        if sort_by_energy:
+            data = sorted(data, key=lambda x: x["formation_energy"], reverse=True)
+        # moves the bulk composition to the top of the list
+        _move_dict_to_start(data, "formula", self.bulk_composition.reduced_formula)
+
+        string = "\\begin{table}[h]\n\\centering\n"
+        string += (
+            "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{"
+            + self.bulk_composition.reduced_formula
+            + "} and all competing phases, with k-meshes used in calculations."
+            + ("}\n" if not prune_polymorphs else " Only the lowest energy polymorphs are included}\n")
+        )
+        string += "\\label{tab:competing_phase_formation_energies}\n"
+        if columns == 1:
+            string += "\\begin{tabular}{ccc}\n"
+            string += "\\hline\n"
+            string += "Formula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\hline \n"
+            for i in data:
+                kpoints = i["kpoints"].split("x")
+                fe = i["formation_energy"]
+                string += (
+                    "\\ce{"
+                    + i["formula"]
+                    + "} & "
+                    + f"{kpoints[0]}$\\times${kpoints[1]}$\\times${kpoints[2]}"
+                    " & " + f"{fe:.3f} \\\\ \n"
+                )
+
+        elif columns == 2:
+            string += "\\begin{tabular}{ccc|ccc}\n"
+            string += "\\hline\n"
+            string += (
+                "Formula & k-mesh & $\\Delta E_f$ (eV) & Formula & k-mesh & $\\Delta E_f$ (eV)\\\\ "
+                "\\hline \n"
+            )
+
+            mid = len(data) // 2
+            first_half = data[:mid]
+            last_half = data[mid:]
+
+            for i, j in zip(first_half, last_half):
+                kpoints = i["kpoints"].split("x")
+                fe = i["formation_energy"]
+                kpoints2 = j["kpoints"].split("x")
+                fe2 = j["formation_energy"]
+                string += (
+                    "\\ce{"
+                    + i["formula"]
+                    + "} & "
+                    + f"{kpoints[0]}$\\times${kpoints[1]}$\\times${kpoints[2]}"
+                    " & "
+                    + f"{fe:.3f} & "
+                    + "\\ce{"
+                    + j["formula"]
+                    + "} & "
+                    + f"{kpoints2[0]}$\\times${kpoints2[1]}$\\times${kpoints2[2]}"
+                    " & " + f"{fe2:.3f} \\\\ \n"
+                )
+
+        string += "\\hline\n"
+        string += "\\end{tabular}\n"
+        string += "\\end{table}"
+
+        return string
+
+
+def _move_dict_to_start(data, key, value):
+    for index, item in enumerate(data):
+        if key in item and item[key] == value:
+            data.insert(0, data.pop(index))
+            return
+
 
 def combine_extrinsic(first, second, extrinsic_species):
     # TODO: Can we just integrate this to `CompetingPhaseAnalyzer`, so you just pass in a list of
diff --git a/examples/dope_chemical_potentials.ipynb b/examples/dope_chemical_potentials.ipynb
index 47bad80b..72af34fd 100644
--- a/examples/dope_chemical_potentials.ipynb
+++ b/examples/dope_chemical_potentials.ipynb
@@ -377,7 +377,9 @@
     "\n",
     "You may want to change the default `ENCUT` (350 eV) or k-point densities that the convergence tests span (5 - 60 kpoints/Å<sup>3</sup> for semiconductors & insulators and 40 - 120 kpoints/Å<sup>3</sup> for metals in steps of 5 kpoints/Å<sup>3</sup>). Note that `ISMEAR = -5` is used for metals by default (better kpoint convergence for energies but should not be used during metal geometry relaxations) and k-point convergence testing is not required for molecules (Γ-point sampling is sufficient).\n",
     "\n",
-    "Note that `doped` generates \"molecule in a box\" structures for the gaseous elemental phases H<sub>2</sub>, O<sub>2</sub>, N<sub>2</sub>, F<sub>2</sub> and Cl<sub>2</sub>. The molecule is placed in a 30 Å x 30 Å x 30 Å box, and relaxed with Γ-point-only k-point sampling.\n",
+    "Note that `doped` generates \"molecule in a box\" structures for the gaseous elemental phases\n",
+    "H<sub>2</sub>, O<sub>2</sub>, N<sub>2</sub>, F<sub>2</sub> and Cl<sub>2</sub>. The molecule is placed in\n",
+    " a slightly-symmetry-broken (to avoid metastable electronic solutions) 30 Å cuboid box, and relaxed with Γ-point-only k-point sampling.\n",
     "\n",
     "The kpoints convergence calculations are set up with:"
    ]
@@ -525,7 +527,11 @@
    "id": "7e57cb66",
    "metadata": {},
    "source": [
-    "Remember that the final `ENCUT` used for the energy calculations should be the same as for your host material & defects, and that you may still need to account for Pulay stress by increasing `ENCUT` for the geometry relaxations (usual rule of thumb being 1.3*converged `ENCUT`) or re-relaxing each structure until the volume change is minimal (roughly <0.3%)."
+    "Remember that the final `ENCUT` used for the energy calculations should be the same as for your host\n",
+    "material & defects, and that you may still need to account for Pulay stress by increasing `ENCUT` for\n",
+    "the geometry relaxations (a typical rule of thumb being 1.3*converged `ENCUT`) or re-relaxing each\n",
+    "structure until the volume change is minimal (roughly <0.3%). This is not the case for the\n",
+    "molecule-in-a-box competing phases however, due to the large simulation box size and fixed volume."
    ]
   },
   {
@@ -730,16 +736,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 1,
    "id": "e273f8fa",
    "metadata": {
-    "ExecuteTime": {
-     "end_time": "2023-08-21T11:54:14.575272Z",
-     "start_time": "2023-08-21T11:54:14.548811Z"
-    },
     "collapsed": false,
     "jupyter": {
      "outputs_hidden": false
+    },
+    "ExecuteTime": {
+     "end_time": "2023-08-23T21:18:44.960159Z",
+     "start_time": "2023-08-23T21:18:39.250023Z"
     }
    },
    "outputs": [],
@@ -750,12 +756,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 2,
    "id": "dd65cb61-bd10-41c8-bf6d-15398841f226",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2023-08-21T11:54:16.335532Z",
-     "start_time": "2023-08-21T11:54:15.916656Z"
+     "end_time": "2023-08-23T21:18:45.128309Z",
+     "start_time": "2023-08-23T21:18:44.961171Z"
     }
    },
    "outputs": [
@@ -789,16 +795,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 3,
    "id": "a3d1b363",
    "metadata": {
-    "ExecuteTime": {
-     "end_time": "2023-08-21T11:51:09.266661Z",
-     "start_time": "2023-08-21T11:51:09.052205Z"
-    },
     "collapsed": false,
     "jupyter": {
      "outputs_hidden": false
+    },
+    "ExecuteTime": {
+     "end_time": "2023-08-23T21:18:45.377891Z",
+     "start_time": "2023-08-23T21:18:45.129641Z"
     }
    },
    "outputs": [
@@ -844,6 +850,73 @@
     "cpa.from_vaspruns(all_paths)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "### Print Formation Energies as a `LaTeX` Table\n",
+    "You may want to print out the parsed calculation results for your competing phases as a LaTeX table, to\n",
+    "be included in supplementary information of journal articles or in theses, to aid open-science and reproducibility.\n",
+    "The formation energies are automatically saved to a `csv` file, which can be converted to LaTeX (see\n",
+    "https://www.tablesgenerator.com/latex_tables) or\n",
+    "Word formats, but you can also print out a nicely-formatted LaTeX table using the `to_LaTeX_table()`\n",
+    "method as shown belo –\n",
+    "either as a one- or two-column table; ordered alphabetially or in descending order of energies."
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\\begin{table}[h]\n",
+      "\\centering\n",
+      "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs are included}\n",
+      "\\label{tab:competing_phase_formation_energies}\n",
+      "\\begin{tabular}{ccc}\n",
+      "\\hline\n",
+      "Formula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\hline \n",
+      "\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 \\\\ \n",
+      "\\ce{O2} & 2$\\times$2$\\times$2 & 0.000 \\\\ \n",
+      "\\ce{Zr} & 9$\\times$9$\\times$5 & 0.000 \\\\ \n",
+      "\\ce{Zr2O} & 5$\\times$5$\\times$2 & -5.729 \\\\ \n",
+      "\\ce{Zr3O} & 5$\\times$5$\\times$5 & -5.987 \\\\ \n",
+      "\\hline\n",
+      "\\end{tabular}\n",
+      "\\end{table}\n"
+     ]
+    }
+   ],
+   "source": [
+    "table = cpa.to_LaTeX_table()\n",
+    "print(table)"
+   ],
+   "metadata": {
+    "collapsed": false,
+    "ExecuteTime": {
+     "end_time": "2023-08-23T21:18:45.385574Z",
+     "start_time": "2023-08-23T21:18:45.382991Z"
+    }
+   }
+  },
+  {
+   "cell_type": "markdown",
+   "source": [
+    "```{note}\n",
+    "The table does not use `booktabs` by default but you can change the horizontal lines (`hline`) to\n",
+    "toprule, midrule and bottomrule to get the desired style effect. The output also relies on using the\n",
+    "`\\ce{}` command from the `mhchem` LaTeX package to get the correct formatting for the chemical formulae.\n",
+    "```"
+   ],
+   "metadata": {
+    "collapsed": false
+   }
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
diff --git a/tests/test_chemical_potentials.py b/tests/test_chemical_potentials.py
index c8f91f87..bcbd0dc7 100644
--- a/tests/test_chemical_potentials.py
+++ b/tests/test_chemical_potentials.py
@@ -217,7 +217,46 @@ def test_cplap_input(self):
         assert (
             "1 Zr 2 O -10.951109995000005\n" not in contents
         )  # shouldn't be included as is a higher energy polymorph of the bulk phase
-        assert not any("2 Zr 1 O" in i for i in contents)  # non-bordering phase
+        assert all("2 Zr 1 O" not in i for i in contents)  # non-bordering phase
+
+    def test_latex_table(self):
+        cpa = chemical_potentials.CompetingPhasesAnalyzer(self.stable_system)
+        path = self.path / "ZrO2"
+        cpa.from_vaspruns(path=path, folder="relax", csv_fname=self.csv_path)
+
+        with self.assertRaises(ValueError):
+            cpa.to_LaTeX_table(columns="M")
+
+        string = cpa.to_LaTeX_table(columns=1)
+        assert (
+            string[28:209]
+            == "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all "
+            "competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs "
+            "are included"
+        )
+        assert len(string) == 586
+        assert string.split("hline")[1] == "\nFormula & k-mesh & $\\Delta E_f$ (eV) \\\\ \\"
+        assert string.split("hline")[2][2:45] == "\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 \\"
+
+        string = cpa.to_LaTeX_table(columns=2)
+        assert (
+            string[28:209]
+            == "\\caption{Formation energies ($\\Delta E_f$) per formula unit of \\ce{ZrO2} and all "
+            "competing phases, with k-meshes used in calculations. Only the lowest energy polymorphs "
+            "are included"
+        )
+        assert (
+            string.split("hline")[1]
+            == "\nFormula & k-mesh & $\\Delta E_f$ (eV) & Formula & k-mesh & $\\Delta E_f$ (eV)\\\\ \\"
+        )
+
+        assert string.split("hline")[2][2:45] == "\\ce{ZrO2} & 3$\\times$3$\\times$3 & -10.975 &"
+        assert len(string) == 579
+
+        cpa_csv = chemical_potentials.CompetingPhasesAnalyzer(self.stable_system)
+        cpa_csv.from_csv(self.csv_path)
+        with self.assertRaises(ValueError):
+            cpa_csv.to_LaTeX_table(columns=2)
 
 
 class BoxedMoleculesTestCase(unittest.TestCase):