remake dopedxpy-sc-fermi

.github/workflows/pip_install_test.yml

@@ -19,7 +19,7 @@ jobs:
 
       matrix:
         os: [ ubuntu-latest, macos-14 ]
-        python-version: [ '3.9', '3.10', '3.11', '3.12' ]
+        python-version: [ '3.9', '3.10', '3.11', '3.12' ]  # 3.13 released in October 2024
         exclude:
           - os: macos-14
             python-version: '3.9'  # Exclude Python 3.9 on macOS, not supported for macOS-14 tests

.github/workflows/test.yml

@@ -14,7 +14,7 @@ jobs:
 
       matrix:
         os: [ ubuntu-latest, macos-14 ]
-        python-version: [ '3.9', '3.10', '3.11', '3.12' ]
+        python-version: [ '3.9', '3.10', '3.11', '3.12' ]  # 3.13 released in October 2024
         exclude:
           - os: macos-14
             python-version: '3.9'  # Exclude Python 3.9 on macOS, not supported for macOS-14 tests

docs/Dev_ToDo.md

@@ -91,7 +91,7 @@
 - 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
 - Test chempot grid plotting tool.
-- `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)
+- `dist_tol` should also group defects for the concentration etc functions, both in `thermodynamics.py` and `fermi_solver.py`, 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)
 - Plotting lines colour updates.
 - Add example to chemical potentials / thermodynamics analysis tutorials of varying chemical potentials as a function of temperature/pressure (i.e. gas phases), using the `Spinney` functions detailed here (https://spinney.readthedocs.io/en/latest/tutorial/chemipots.html#including-temperature-and-pressure-effects-through-the-gas-phase-chemical-potentials) or possibly `DefAP` functions otherwise.
 - Add note to chempots tutorial that if your bulk phase is lower energy than its version on the MP (e.g. distorted perovskite), then you should use this for your bulk competing phase calculation.

docs/Tutorials.rst

@@ -15,5 +15,6 @@ tutorials:
    thermodynamics_tutorial
    chemical_potentials_tutorial
    advanced_analysis_tutorial
+   fermisolver_tutorial
    plotting_customisation_tutorial
    GGA_workflow_tutorial

docs/fermisolver_tutorial.ipynb

@@ -0,0 +1 @@
+../examples/fermisolver_tutorial.ipynb

doped/chemical_potentials.py

@@ -3367,3 +3367,142 @@ def combine_extrinsic(first, second, extrinsic_species):
         "limits": new_limits,
         "limits_wrt_el_refs": new_limits_wrt_el_refs,
     }
+
+
+class ChemicalPotentialGrid:
+    """
+    A class to represent a grid in chemical potential space and to perform
+    operations such as generating a grid within the convex hull of given
+    vertices (chemical potential limits).
+
+    This class provides methods for handling and manipulating chemical
+    potential data, including the creation of a grid that spans a specified
+    chemical potential space.
+    """
+
+    def __init__(self, chempots: dict[str, Any]):
+        r"""
+        Initializes the ``ChemicalPotentialGrid`` with chemical potential data.
+
+        This constructor takes a dictionary of chemical potentials and sets up
+        the initial vertices of the grid.
+
+        Args:
+            chempots (dict):
+                Dictionary of chemical potentials for the grid. This can have
+                the form of ``{"limits": [{'limit': [chempot_dict]}], ...}``
+                (the format generated by ``doped``\'s chemical potential parsing
+                functions), or alternatively can be a dictionary of the form
+                ``{'limit': [chempot_dict, ...]}`` (i.e. matching the format of
+                ``chempots["limits_wrt_el_refs"]`` from the ``doped`` ``chempots``
+                dict) where the keys are the limit names (e.g. "Cd-CdTe", "Cd-rich"
+                etc) and the values are dictionaries of a single chemical potential
+                limit in the format: ``{element symbol: chemical potential}``.
+
+                If ``chempots`` in the ``doped`` format is supplied, then the
+                chemical potentials `with respect to the elemental reference
+                energies` will be used (i.e. ``chempots["limits_wrt_el_refs"]``)!
+        """
+        unformatted_chempots_dict = chempots.get("limits_wrt_el_refs", chempots)
+        test_elt = Element("H")
+        formatted_chempots_dict = {
+            limit: {
+                f"μ_{k}" if test_elt.is_valid_symbol(k) else k: v
+                for (k, v) in unformatted_chempots_subdict.items()
+            }
+            for limit, unformatted_chempots_subdict in unformatted_chempots_dict.items()
+        }
+
+        self.vertices = pd.DataFrame.from_dict(formatted_chempots_dict, orient="index")
+
+    def get_grid(self, n_points: int = 100) -> pd.DataFrame:
+        """
+        Generates a grid within the convex hull of the vertices and
+        interpolates the dependent variable values.
+
+        This method creates a grid of points that spans the chemical potential
+        space defined by the vertices. It ensures that the generated points lie
+        within the convex hull of the provided vertices and interpolates the
+        chemical potential values at these points.
+
+        Args:
+            n_points (int):
+                The number of points to generate along each axis of the grid.
+                Note that this may not always be the final number of points in the grid,
+                as points lying outside the convex hull are excluded.
+                Default is 100.
+
+        Returns:
+            pd.DataFrame:
+                A ``DataFrame`` containing the points within the convex hull,
+                along with their corresponding interpolated chemical potential values.
+                Each row represents a point in the grid with associated chemical
+                potential values.
+        """
+        return self.grid_from_dataframe(self.vertices, n_points)
+
+    @staticmethod
+    def grid_from_dataframe(mu_dataframe: pd.DataFrame, n_points: int = 100) -> pd.DataFrame:
+        r"""
+        Generates a grid within the convex hull of the vertices.
+
+        This method creates a grid of points within the convex hull
+        defined by the input ``DataFrame``\.
+        It interpolates the values of chemical potentials over this
+        grid, ensuring that all generated points lie within the convex
+        hull of the given vertices.
+
+        Args:
+            mu_dataframe (pd.DataFrame):
+                A ``DataFrame`` containing the chemical potential data,
+                with the last column representing the dependent variable
+                and the preceding columns representing the independent
+                variables.
+            n_points (int):
+                The number of points to generate along each axis of the
+                grid. Note that this may not always be the final number
+                of points in the grid, as points lying outside the convex
+                hull are excluded. Defaults to 100.
+
+        Returns:
+            pd.DataFrame:
+                A ``DataFrame`` containing the points within the convex
+                hull along with their corresponding interpolated values of
+                the dependent variable. Each row represents a point in the
+                grid.
+        """
+        dependent_variable = mu_dataframe.columns[-1]
+        dependent_var = mu_dataframe[dependent_variable].to_numpy()
+        independent_vars = mu_dataframe.drop(columns=dependent_variable)
+
+        n_dims = independent_vars.shape[1]  # Get the number of independent variables (dimensions)
+
+        # Get the convex hull of the vertices
+        hull = ConvexHull(independent_vars.values)
+
+        # Create a dense grid that covers the entire range of the vertices
+        grid_ranges = [
+            np.linspace(independent_vars.iloc[:, i].min(), independent_vars.iloc[:, i].max(), n_points)
+            for i in range(n_dims)
+        ]
+        grid = np.meshgrid(*grid_ranges, indexing="ij")  # Create N-dimensional grid
+        grid_points = np.vstack([g.ravel() for g in grid]).T  # Flatten the grid to points
+
+        # Delaunay triangulation to get points inside the convex hull
+        delaunay = Delaunay(hull.points[hull.vertices])
+        inside_hull = delaunay.find_simplex(grid_points) >= 0
+        points_inside = grid_points[inside_hull]
+
+        # Interpolate the values to get the dependent chemical potential
+        values_inside = griddata(independent_vars.values, dependent_var, points_inside, method="linear")
+
+        # Combine points with their corresponding interpolated values
+        grid_with_values = np.hstack((points_inside, values_inside.reshape(-1, 1)))
+
+        # Add vertices to the grid
+        grid_with_values = np.vstack((grid_with_values, mu_dataframe.to_numpy()))
+
+        return pd.DataFrame(
+            grid_with_values,
+            columns=[*list(independent_vars.columns), dependent_variable],
+        )

doped/core.py

@@ -1361,7 +1361,9 @@ def _no_chempots_warning(property="Formation energies (and concentrations)"):
     )
 
 
-def _get_dft_chempots(chempots, el_refs, limit):
+def _get_dft_chempots(
+    chempots: Optional[dict], el_refs: Optional[dict] = None, limit: Optional[str] = None
+):
     """
     Parse the DFT chempots from the input chempots and limit.
     """