Use parametrization to combine tests for different backends, and agai…

…n test actual values

doped/thermodynamics.py

@@ -4660,7 +4660,9 @@ def _get_single_chempot_dict(
     ) -> tuple[dict[str, float], Any]:
         """
         Get the chemical potentials for a single limit (``limit``) from the
-        ``chempots`` (or ``self.defect_thermodynamics.chempots``) dictionary.
+        ``chempots`` (or ``self.defect_thermodynamics.chempots``) dictionary,
+        giving the chemical potentials `with respect to the elemental
+        references` (i.e. from ``"limits_wrt_el_refs"``).
 
         Returns a `single` chemical potential dictionary for the specified limit.
         """
@@ -4705,7 +4707,7 @@ def equilibrium_solve(
                 should be a dictionary of chemical potentials for a single limit
                 (``limit``), in the format: ``{element symbol: chemical potential}``.
                 If ``el_refs`` is provided or set in ``self.defect_thermodynamics.el_refs``
-                then it is the formal chemical potentials (i.e. relative to the elemental
+                then it is the `formal` chemical potentials (i.e. relative to the elemental
                 reference energies) that should be given here, otherwise the absolute
                 (DFT) chemical potentials should be given.
             el_refs (dict[str, float]):
@@ -4808,28 +4810,35 @@ def equilibrium_solve(
                 conc_dict = defect_system.concentration_dict()
             data = []
 
+            df_dict_base = {
+                "Temperature": defect_system.temperature,
+                "Fermi Level": conc_dict["Fermi Energy"],
+                "Electrons (cm^-3)": conc_dict["n0"],
+                "Holes (cm^-3)": conc_dict["p0"],
+            }
             for k, v in conc_dict.items():
                 if k not in ["Fermi Energy", "n0", "p0", "Dopant"]:
                     row = {
                         "Defect": k,
                         "Concentration (cm^-3)": v,
-                        "Temperature": defect_system.temperature,
-                        "Fermi Level": conc_dict["Fermi Energy"],
-                        "Electrons (cm^-3)": conc_dict["n0"],
-                        "Holes (cm^-3)": conc_dict["p0"],
+                        **df_dict_base,
                     }
-                    if "Dopant" in conc_dict:
-                        row["Dopant (cm^-3)"] = conc_dict["Dopant"]
                     data.append(row)
+            if "Dopant" in conc_dict:
+                data.append(
+                    {
+                        "Defect": "Dopant",
+                        "Concentration (cm^-3)": conc_dict["Dopant"],
+                        **df_dict_base,
+                    }
+                )
 
             concentrations = pd.DataFrame(data)
             concentrations = concentrations.set_index("Defect", drop=True)
 
         if append_chempots:
             for key, value in single_chempot_dict.items():
                 concentrations[f"μ_{key}"] = value
-            if effective_dopant_concentration:
-                concentrations["Dopant (cm^-3)"] = effective_dopant_concentration
 
         return concentrations
 
@@ -4910,7 +4919,7 @@ def pseudo_equilibrium_solve(
                 potentials for a single limit (``limit``), in the format:
                 ``{element symbol: chemical potential}``.
                 If ``el_refs`` is provided or set in ``self.defect_thermodynamics.el_refs``
-                then it is the formal chemical potentials (i.e. relative to the elemental
+                then it is the `formal` chemical potentials (i.e. relative to the elemental
                 reference energies) that should be given here, otherwise the absolute
                 (DFT) chemical potentials should be given.
             el_refs (dict[str, float]):
@@ -5062,8 +5071,6 @@ def pseudo_equilibrium_solve(
         if append_chempots:
             for key, value in single_chempot_dict.items():
                 concentrations[f"μ_{key}"] = value
-            if effective_dopant_concentration:
-                concentrations["Dopant (cm^-3)"] = effective_dopant_concentration
 
         return concentrations
 
@@ -6360,7 +6367,7 @@ def _generate_defect_system(
                 should be a dictionary of chemical potentials for a single limit
                 (``limit``), in the format: ``{element symbol: chemical potential}``.
                 If ``el_refs`` is provided or set in ``self.defect_thermodynamics.el_refs``
-                then it is the formal chemical potentials (i.e. relative to the elemental
+                then it is the `formal` chemical potentials (i.e. relative to the elemental
                 reference energies) that should be given here, otherwise the absolute
                 (DFT) chemical potentials should be given.
             el_refs (dict[str, float]):
@@ -6402,6 +6409,9 @@ def _generate_defect_system(
         self._check_required_backend_and_error("py-sc-fermi")
         assert self._DefectSpecies
         assert self._DefectSystem
+        single_chempot_dict, el_refs = self.defect_thermodynamics._get_chempots(
+            single_chempot_dict, el_refs
+        )  # returns self.defect_thermodynamics.chempots/self.defect_thermodynamics.el_refs if None
         dft_chempots = _get_dft_chempots(single_chempot_dict, el_refs)
 
         defect_species = []  # dicts of: {"charge_states": {...}, "nsites": X, "name": label}
@@ -6535,7 +6545,7 @@ def _generate_annealed_defect_system(
                 potentials for a single limit (``limit``), in the format:
                 ``{element symbol: chemical potential}``.
                 If ``el_refs`` is provided or set in
-                ``self.defect_thermodynamics.el_refs`` then it is the formal
+                ``self.defect_thermodynamics.el_refs`` then it is the `formal`
                 chemical potentials (i.e. relative to the elemental reference
                 energies) that should be given here, otherwise the absolute
                 (DFT) chemical potentials should be given.

examples/fermisolver_tutorial.ipynb

@@ -34,8 +34,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2025-01-01T23:45:58.269116Z",
-     "start_time": "2025-01-01T23:45:54.202069Z"
+     "end_time": "2025-01-06T17:34:17.117193Z",
+     "start_time": "2025-01-06T17:34:13.398285Z"
     }
    },
    "source": [
@@ -56,8 +56,8 @@
   {
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2025-01-01T23:45:59.171110Z",
-     "start_time": "2025-01-01T23:45:58.274702Z"
+     "end_time": "2025-01-06T17:34:21.979462Z",
+     "start_time": "2025-01-06T17:34:21.313094Z"
     }
    },
    "cell_type": "code",
@@ -102,8 +102,8 @@
    "cell_type": "code",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2025-01-01T23:46:00.846448Z",
-     "start_time": "2025-01-01T23:45:59.212926Z"
+     "end_time": "2025-01-06T17:34:32.343730Z",
+     "start_time": "2025-01-06T17:34:30.398465Z"
     }
    },
    "source": [
@@ -119,16 +119,16 @@
     "\n",
     "# or they can be specified as a dictionary of chemical potentials referenced to the\n",
     "# elements in the bulk structure\n",
-    "temperature_df_py = py_fs.scan_temperature(\n",
+    "temperature_df_py = py_fs.scan_temperature(  # manually specify Cd-rich chempots\n",
     "    chempots={'Cd': 0.0, 'Te': -1.2513}, temperature_range=temperatures)"
    ],
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "100%|██████████| 23/23 [00:00<00:00, 58.65it/s]\n",
-      "100%|██████████| 23/23 [00:01<00:00, 18.87it/s] \n"
+      "100%|██████████| 23/23 [00:00<00:00, 37.75it/s]\n",
+      "100%|██████████| 23/23 [00:01<00:00, 17.61it/s]\n"
      ]
     }
    ],
@@ -275,8 +275,8 @@
   {
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2025-01-01T23:46:00.900059Z",
-     "start_time": "2025-01-01T23:46:00.896542Z"
+     "end_time": "2025-01-06T17:34:51.896856Z",
+     "start_time": "2025-01-06T17:34:51.890124Z"
     }
    },
    "cell_type": "code",
@@ -288,13 +288,13 @@
     "plt.style.use(f\"{doped.__path__[0]}/utils/doped.mplstyle\")  # use doped style"
    ],
    "outputs": [],
-   "execution_count": 5
+   "execution_count": 4
   },
   {
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2025-01-01T23:46:01.411185Z",
-     "start_time": "2025-01-01T23:46:01.082664Z"
+     "end_time": "2025-01-06T17:34:55.213415Z",
+     "start_time": "2025-01-06T17:34:54.849573Z"
     }
    },
    "cell_type": "code",
@@ -357,7 +357,7 @@
      "output_type": "display_data"
     }
    ],
-   "execution_count": 6
+   "execution_count": 5
   },
   {
    "cell_type": "markdown",

tests/test_fermisolver.py

@@ -462,12 +462,12 @@ def test_equilibrium_solve(self, backend):
             "Electrons (cm^-3)",
             "Holes (cm^-3)",
             "Temperature",
-            "Dopant (cm^-3)",
         ]:
             assert i in concentrations.columns, f"Missing column: {i}"
 
         # Check that concentrations are reasonable numbers
         assert np.all(concentrations["Concentration (cm^-3)"] >= 0)
+        assert np.isclose(concentrations["Temperature"].iloc[0], 300)
         # Check appended chemical potentials
         for element in single_chempot_dict:
             assert f"μ_{element}" in concentrations.columns
@@ -480,9 +480,25 @@ def test_equilibrium_solve(self, backend):
         doped_e_h = get_e_h_concs(self.CdTe_fermi_dos, expected_fermi_level + self.example_thermo.vbm, 300)
         assert np.isclose(concentrations["Electrons (cm^-3)"].iloc[0], doped_e_h[0], rtol=1e-3)
         assert np.isclose(concentrations["Holes (cm^-3)"].iloc[0], doped_e_h[1], rtol=1e-3)
-        # doped_defect_concs = self.example_thermo.get_equilibrium_concentrations(
-        #     fermi_level=expected_fermi_level, limit="Te-rich", temperature=300
-        # )
+        doped_defect_concs = self.example_thermo.get_equilibrium_concentrations(
+            fermi_level=expected_fermi_level,
+            limit="Te-rich",
+            temperature=300,
+            per_charge=False,
+            skip_formatting=True,
+        )  # currently FermiSolver only supports per charge=False
+
+        fermisolver_concentrations = concentrations["Concentration (cm^-3)"]
+        fermisolver_concentrations = fermisolver_concentrations.drop("Dopant")
+        # drop Dopant row, not included with ``DefectThermodynamics.get_equilibrium_concentrations()``
+        if backend == "py-sc-fermi":
+            fermisolver_concentrations["Te_i_Td_Te2.83"] = fermisolver_concentrations["Te_i_Td_Te2.83_a"]
+            fermisolver_concentrations = fermisolver_concentrations.drop(
+                ["Te_i_Td_Te2.83_a", "Te_i_Td_Te2.83_b"],
+            )
+        pd.testing.assert_series_equal(
+            doped_defect_concs["Concentration (cm^-3)"], fermisolver_concentrations, rtol=1e-3
+        )  # also checks the index and ordering
 
     def test_equilibrium_solve_mocked_py_sc_fermi_backend(self):
         """
@@ -521,7 +537,6 @@ def test_equilibrium_solve_mocked_py_sc_fermi_backend(self):
             "Electrons (cm^-3)",
             "Holes (cm^-3)",
             "Temperature",
-            "Dopant (cm^-3)",
         ]:
             assert i in concentrations.columns, f"Missing column: {i}"
 
@@ -535,14 +550,16 @@ def test_equilibrium_solve_mocked_py_sc_fermi_backend(self):
 
     # Tests for pseudo_equilibrium_solve
 
-    def test_pseudo_equilibrium_solve_doped_backend(self):
+    @parameterize_backend()
+    def test_pseudo_equilibrium_solve_doped_backend(self, backend):
         """
-        Test pseudo_equilibrium_solve method for doped backend.
+        Test ``pseudo_equilibrium_solve`` method for both backends.
         """
-        single_chempot_dict, el_refs = self.solver_doped._get_single_chempot_dict(limit="Te-rich")
+        solver = self.solver_doped if backend == "doped" else self.solver_py_sc_fermi
+        single_chempot_dict, el_refs = solver._get_single_chempot_dict(limit="Te-rich")
 
         # Call the method
-        concentrations = self.solver_doped.pseudo_equilibrium_solve(
+        concentrations = solver.pseudo_equilibrium_solve(
             annealing_temperature=800,
             single_chempot_dict=single_chempot_dict,
             el_refs=el_refs,
@@ -551,19 +568,45 @@ def test_pseudo_equilibrium_solve_doped_backend(self):
             append_chempots=True,
         )
 
-        # Assertions
-        assert "Fermi Level" in concentrations.columns
-        assert "Electrons (cm^-3)" in concentrations.columns
-        assert "Holes (cm^-3)" in concentrations.columns
-        assert "Annealing Temperature" in concentrations.columns
-        assert "Quenched Temperature" in concentrations.columns
+        for i in [
+            "Fermi Level",
+            "Electrons (cm^-3)",
+            "Holes (cm^-3)",
+            "Annealing Temperature",
+            "Quenched Temperature",
+        ]:
+            assert i in concentrations.columns, f"Missing column: {i}"
+
         # Check that concentrations are reasonable numbers
         assert np.all(concentrations["Concentration (cm^-3)"] >= 0)
+        assert np.isclose(concentrations["Quenched Temperature"].iloc[0], 300)
+        assert np.isclose(concentrations["Annealing Temperature"].iloc[0], 800)
         # Check appended chemical potentials
         for element in single_chempot_dict:
             assert f"μ_{element}" in concentrations.columns
             assert concentrations[f"μ_{element}"].iloc[0] == single_chempot_dict[element]
 
+        fermi_level, e_conc, h_conc, conc_df = self.example_thermo.get_fermi_level_and_concentrations(
+            annealing_temperature=800,
+            effective_dopant_concentration=1e16,
+            limit="Te-rich",
+            per_charge=False,  # currently FermiSolver only supports per charge=False
+            skip_formatting=True,
+        )
+        assert np.isclose(concentrations["Fermi Level"].iloc[0], fermi_level)
+        assert np.isclose(concentrations["Electrons (cm^-3)"].iloc[0], e_conc, rtol=1e-3)
+        assert np.isclose(concentrations["Holes (cm^-3)"].iloc[0], h_conc, rtol=1e-3)
+
+        fermisolver_concentrations = concentrations["Concentration (cm^-3)"]
+        if backend == "py-sc-fermi":
+            fermisolver_concentrations["Te_i_Td_Te2.83"] = fermisolver_concentrations["Te_i_Td_Te2.83_a"]
+            fermisolver_concentrations = fermisolver_concentrations.drop(
+                ["Te_i_Td_Te2.83_a", "Te_i_Td_Te2.83_b"],
+            )
+        pd.testing.assert_series_equal(
+            conc_df["Concentration (cm^-3)"], fermisolver_concentrations, rtol=1e-3
+        )  # also checks the index and ordering
+
     def test_pseudo_equilibrium_solve_py_sc_fermi_backend(self):
         """
         Test pseudo_equilibrium_solve method for py-sc-fermi backend with
@@ -600,14 +643,22 @@ def test_pseudo_equilibrium_solve_py_sc_fermi_backend(self):
         )
 
         # Assertions
-        assert "Fermi Level" in concentrations.columns
-        assert "Electrons (cm^-3)" in concentrations.columns
-        assert "Holes (cm^-3)" in concentrations.columns
-        assert "Annealing Temperature" in concentrations.columns
-        assert "Quenched Temperature" in concentrations.columns
+        for i in [
+            "Fermi Level",
+            "Electrons (cm^-3)",
+            "Holes (cm^-3)",
+            "Annealing Temperature",
+            "Quenched Temperature",
+        ]:
+            assert i in concentrations.columns, f"Missing column: {i}"
+
         # Check defects are included
         assert "defect1" in concentrations.index
         assert "defect2" in concentrations.index
+
+        assert np.isclose(concentrations["Quenched Temperature"].iloc[0], 300)
+        assert np.isclose(concentrations["Annealing Temperature"].iloc[0], 800)
+
         # Check appended chemical potentials
         for element in single_chempot_dict:
             assert f"μ_{element}" in concentrations.columns