Update over-dependent chempot handling and TODOs/notes

docs/Dev_ToDo.md

@@ -2,10 +2,8 @@
 ## Chemical potential
 - Check through chemical potential TO-DOs. Need to recheck validity of approximations used for extrinsic competing phases (and code for this). Proper `vasp_std` setup (with `NKRED` folders like for defect calcs) and `vasp_ncl` generation.
 - Efficient generation of competing phases for which there are many polymorphs? See SK notes from CdTe competing phases.
-- Update chemical potential tools to work with new Materials Project API. Currently, supplying an API key for the new Materials Project API returns entries which do not have `e_above_hull` as a property, and so crashes. Ideally would be good to be compatible with both the legacy and new API, which should be fairly straightforward (try importing MPRester from mp_api client except ImportError import from pmg then will need to make a whole separate query/search because `band_gap` and `total_magnetisation` no longer accessible from `get_entries`). See https://docs.materialsproject.org/downloading-data/using-the-api
 - Publication ready chemical potential diagram plotting tool as in Adam Jackson's `plot-cplap-ternary` (3D) and Sungyhun's `cplapy` (4D) (see `doped_chempot_plotting_example.ipynb`; code there, just needs to be implemented in module functions). `ChemicalPotentialGrid` in `py-sc-fermi` interface could be quite useful for this? (Worth moving that part of code out of `interface` subpackage?) Also can see https://github.com/materialsproject/pymatgen-analysis-defects/blob/main/pymatgen/analysis/defects/plotting/phases.py for 2D chempot plotting.
   - Also see `Cs2SnTiI6` notebooks for template code for this.
-- Note in tutorial that LaTeX table generator website can also be used with the `to_csv()` function to generate LaTeX tables for the competing phases.
 
 ## Post-processing / analysis / plotting
 - Better automatic defect formation energy plot colour handling (auto-change colormap based on number of defects, set similar colours for similar defects (types and inequivalent sites)) – and more customisable?
@@ -88,6 +86,7 @@
 
 ## SK To-Do for next update:
 - `doped` repo/docs cleanup `TODO`s above
+- Note in chempots tutorial that LaTeX table generator website can also be used with the `to_csv()` function to generate LaTeX tables for the competing phases.
 - Quick-start tutorial suggested by Alex G
 - Add chempot grid plotting tool, shown in `JOSS_plots` using Alex's chemical potential grid, and test (and remove TODO from JOSS plots notebook).
 - `dist_tol` should also group defects for the concentration etc functions, currently doesn't (e.g. `CdTe_thermo.get_equilibrium_concentrations(limit="Te-rich", per_charge=False, fermi_level=0.5)` and `CdTe_thermo.dist_tol=10; CdTe_thermo.get_equilibrium_concentrations(limit="Te-rich", per_charge=False, fermi_level=0.5)`, same output)

doped/chemical_potentials.py

@@ -744,9 +744,8 @@ def __init__(
         self.bulk_composition = Composition(composition)
         self.chemsys = list(self.bulk_composition.as_dict().keys())
 
-        # TODO: Update installation pages and any mentions of 'legacy' in code
+        # TODO: Update installation pages, docs and tutorials
         # TODO: Add tests with new API keys
-        # TODO: Update docs and tutorials
 
         # get all entries in the chemical system:
         self.MP_full_pd_entries, self.property_key_dict, self.property_data_fields = (
@@ -1351,16 +1350,16 @@ def __init__(
                     MP_extrinsic_bordering_phases: list[str] = []
 
                     for limit in MP_extrinsic_gga_chempots:
-                        # note that the the number of phases in equilibria at each vertex (limit) is
-                        # equal to the number of elements in the chemical system (here being the
-                        # host composition plus the extrinsic species)
+                        # note that the number of phases in equilibria at each vertex (limit) is equal
+                        # to the number of elements in the chemical system (here being the host
+                        # composition plus the extrinsic species)
                         extrinsic_bordering_phases = {
                             phase for phase in limit.split("-") if sub_el in phase
                         }
                         # only add to MP_extrinsic_bordering_phases when only 1 extrinsic bordering phase
                         # (i.e. ``full_sub_approach=False`` behaviour):
-                        if len(
-                            extrinsic_bordering_phases
+                        if len(  # this should always give the same number of facets as the bulk PD
+                            extrinsic_bordering_phases  # TODO: Explicitly test this for all cases in tests
                         ) == 1 and not extrinsic_bordering_phases.issubset(MP_extrinsic_bordering_phases):
                             MP_extrinsic_bordering_phases.extend(extrinsic_bordering_phases)
 
@@ -1371,19 +1370,25 @@ def __init__(
                         or (entry.is_element and sub_el in entry.name)
                     ]
 
-                    # check that extrinsic competing phases list is not empty (can happen with
-                    # 'over-dependent' limits); if so then set full_sub_approach = True and re-run
-                    # the extrinsic phase addition process
+                    # check that extrinsic competing phases list is not empty (according to PyCDT
+                    # chemical potential handling this can happen (despite purposely neglecting these
+                    # "over-dependent" facets above), but no known cases... (apart from when `extrinsic`
+                    # actually contains an intrinsic element, which we handle above anyway)
                     if not single_bordering_sub_el_entries:
-                        warnings.warn(
+                        # warnings.warn(
+                        #     f"Determined chemical potentials to be over-dependent on the extrinsic "
+                        #     f"species {sub_el}, meaning we need to revert to `full_sub_approach = True` "
+                        #     f"for this species."
+                        # )  # Revert to this handling if we ever find a case of this actually happening
+                        # self.entries += sub_el_entries
+                        raise RuntimeError(
                             f"Determined chemical potentials to be over-dependent on the extrinsic "
-                            f"species {sub_el}, meaning we need to revert to `full_sub_approach = True` "
-                            f"for this species."
+                            f"species {sub_el} despite `full_sub_approach=False`, which shouldn't happen. "
+                            f"Please report this to the developers on the GitHub issues page: "
+                            f"https://github.com/SMTG-Bham/doped/issues"
                         )
-                        self.entries += sub_el_entries
 
-                    else:
-                        self.entries += single_bordering_sub_el_entries
+                    self.entries += single_bordering_sub_el_entries
 
         self.MP_full_pd_entries.sort(  # sort by energy above hull, num_species, then alphabetically
             key=lambda x: _pd_entries_sorting_func(x, self.legacy_MP)