Tutorial updates and minor code cleanup

doped/interface/fermi_solver.py

@@ -23,7 +23,11 @@
 from scipy.interpolate import griddata
 from scipy.spatial import ConvexHull, Delaunay
 
-from doped.thermodynamics import _parse_chempots, get_rich_poor_limit_dict
+from doped.thermodynamics import (
+    _add_effective_dopant_concentration,
+    _parse_chempots,
+    get_rich_poor_limit_dict,
+)
 from doped.utils.parsing import get_neutral_nelect_from_vasprun
 
 if TYPE_CHECKING:
@@ -489,7 +493,8 @@ def _add_effective_dopant_concentration_and_solve_pseudo(
     def scan_chemical_potential_grid(
         self,
         dependent_variable: str,
-        chempots=None,
+        chempots: Optional[dict] = None,
+        el_refs: Optional[dict] = None,
         n_points: int = 10,
         temperature: float = 300,
         annealing_temperature: Optional[float] = None,
@@ -504,7 +509,8 @@ def scan_chemical_potential_grid(
 
         Args:
             dependent_variable (str): the dependent variable to scan
-            chempots (dict): chemical potentials to scan
+            chempots (dict): chemical potentials to scan. Uses ``self.chempots`` if not provided.
+            el_refs (dict): elemental reference energies for the chemical potentials (optional)
             n_points (int): number of points to scan
             temperature (float): temperature to solve at
             annealing_temperature (float): temperature to anneal at
@@ -516,13 +522,10 @@ def scan_chemical_potential_grid(
         Returns:
             pd.DataFrame: DataFrame containing the Fermi energy solutions at the grid points
         """
-        if chempots is None:
-            if self.chempots is None or "limits_wrt_el_refs" not in self.chempots:
-                raise ValueError(
-                    "self.chempots or self.chempots['limits_wrt_el_refs'] is None or missing."
-                )
-            chempots = self.chempots["limits_wrt_el_refs"]
-
+        # Parse chemical potentials, either using input values (after formatting them in the doped format)
+        # or using the class attributes if set:
+        chempots, _el_refs = _parse_chempots(chempots or self.chempots, el_refs or self.el_refs)
+        assert chempots is not None
         grid = ChemicalPotentialGrid.from_chempots(chempots).get_grid(dependent_variable, n_points)
 
         if annealing_temperature is not None and quenching_temperature is not None:
@@ -558,7 +561,8 @@ def min_max_X(
         dependent_chempot,
         target,
         min_or_max,
-        chempots=None,
+        chempots: Optional[dict] = None,
+        el_refs: Optional[dict] = None,
         tolerance=0.01,
         n_points=10,
         temperature=300,
@@ -575,16 +579,17 @@ def min_max_X(
         chemical potential that minimizes or maximizes the target variable
         until the target value no longer changes by more than the tolerance.
         """
-        if chempots is None:
-            chempots = self.chempots["limits_wrt_el_refs"]
-
+        chempots, _el_refs = _parse_chempots(chempots or self.chempots, el_refs or self.el_refs)
+        assert chempots is not None
         starting_grid = ChemicalPotentialGrid.from_chempots(chempots)
         current_vertices = starting_grid.vertices
         chempots_labels = list(current_vertices.columns)
         previous_value = None
 
         while True:
             # Solve and append chemical potentials
+            print(f"Scanning with grid size {starting_grid.vertices.shape}")
+            print(starting_grid.vertices)  # debugging, remove when fixed! TODO
             if annealing_temperature is not None and quenching_temperature is not None:
                 all_data = Parallel(n_jobs=processes)(
                     delayed(self._solve_and_append_chempots_pseudo)(
@@ -762,7 +767,7 @@ def equilibrium_solve(
             per_site=False,
             skip_formatting=False,
         )
-        concentrations = self.defect_thermodynamics._add_effective_dopant_concentration(
+        concentrations = _add_effective_dopant_concentration(
             concentrations, effective_dopant_concentration
         )
         new_columns = {

doped/thermodynamics.py

@@ -1522,7 +1522,7 @@ def _get_total_q(fermi_level):
                 lean=True,
             )
             # add effective dopant concentration if supplied:
-            conc_df = self._add_effective_dopant_concentration(conc_df, effective_dopant_concentration)
+            conc_df = _add_effective_dopant_concentration(conc_df, effective_dopant_concentration)
             qd_tot = (conc_df["Charge"] * conc_df["Concentration (cm^-3)"]).sum()
             qd_tot += get_doping(
                 fermi_dos=self.fermi_dos, fermi_level=fermi_level + self.vbm, temperature=temperature
@@ -1779,7 +1779,7 @@ def get_quenched_fermi_level_and_concentrations(
                 per_charge=False,  # give total concentrations for each defect
                 lean=True,
             )
-            annealing_defect_concentrations = self._add_effective_dopant_concentration(
+            annealing_defect_concentrations = _add_effective_dopant_concentration(
                 annealing_defect_concentrations, effective_dopant_concentration
             )  # add effective dopant concentration if supplied
             total_concentrations = dict(  # {Defect: Total Concentration (cm^-3)}
@@ -1801,7 +1801,7 @@ def _get_constrained_concentrations(
                     skip_formatting=True,
                     lean=lean,
                 )
-                conc_df = self._add_effective_dopant_concentration(conc_df, effective_dopant_concentration)
+                conc_df = _add_effective_dopant_concentration(conc_df, effective_dopant_concentration)
                 conc_df["Total Concentration (cm^-3)"] = conc_df["Defect"].map(total_concentrations)
                 conc_df["Concentration (cm^-3)"] = (  # set total concentration to match annealing conc
                     conc_df["Concentration (cm^-3)"]  # but with same relative concentrations