diff --git a/doped/interface/fermi_solver.py b/doped/interface/fermi_solver.py index beaff8ce..9716791f 100644 --- a/doped/interface/fermi_solver.py +++ b/doped/interface/fermi_solver.py @@ -24,11 +24,8 @@ from scipy.spatial import ConvexHull, Delaunay from tqdm import tqdm -from doped.thermodynamics import ( - _add_effective_dopant_concentration, - _parse_chempots, - get_rich_poor_limit_dict, -) +from doped.core import _get_dft_chempots +from doped.thermodynamics import _add_effective_dopant_concentration, _parse_chempots, _parse_limit if TYPE_CHECKING: from pymatgen.util.typing import PathLike @@ -55,12 +52,10 @@ def _get_label_and_charge(name: str) -> tuple[str, int]: label = name[:last_underscore] if last_underscore != -1 else name charge_str = name[last_underscore + 1 :] if last_underscore != -1 else None - # Initialize charge with a default value - charge = 0 + charge = 0 # Initialize charge with a default value if charge_str is not None: with contextlib.suppress(ValueError): charge = int(charge_str) - # If charge_str cannot be converted to int, charge remains its default value return label, charge @@ -148,6 +143,12 @@ def __init__( necessary attributes, which includes loading the bulk density of states (DOS) data from the provided VASP output. + If using the ``py-sc-fermi`` backend (currently required for the + ``free_defects`` and ``fix_charge_states`` options in the scanning + functions), please cite the code paper: + Squires et al., (2023). Journal of Open Source Software, 8(82), 4962 + https://doi.org/10.21105/joss.04962 + Args: defect_thermodynamics (DefectThermodynamics): A ``DefectThermodynamics`` object, providing access to defect @@ -176,7 +177,7 @@ def __init__( If ``None`` (default), will use ``DefectThermodynamics.chempots`` (= 0 for all chemical potentials by default, if not set). - This can have the form of ``{"limits": [{'limit': [chempot_dict]}]}`` + This can have the form of ``{"limits": [{'limit': [chempot_dict]}], ...}`` (the format generated by ``doped``\'s chemical potential parsing functions (see tutorials)) and specific limits (chemical potential limits) can then be chosen using ``limit``. @@ -188,6 +189,13 @@ def __init__( in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + If provided here, then ``defect_thermodynamics.chempots`` is set to this + input. + + Chemical potentials can also be supplied later in each analysis function, or + set using ``FermiSolver.defect_thermodynamics.chempots = ...`` or + ``DefectThermodynamics.chempots = ...`` (with the same input options). el_refs (Optional[dict]): Dictionary of elemental reference energies for the chemical potentials in the format: @@ -196,7 +204,15 @@ def __init__( ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in ``DefectThermodynamics`` in the format generated by ``doped`` (see tutorials), or if ``DefectThermodynamics.el_refs`` has been - set. Default is None. + set. + + If provided here, then ``defect_thermodynamics.el_refs`` is set to this + input. + + Elemental reference energies can also be supplied later in each analysis + function, or set using ``FermiSolver.defect_thermodynamics.el_refs = ...`` + or ``DefectThermodynamics.el_refs = ...`` (with the same input options). + Default is ``None``. backend (Optional[str]): The code backend to use for the thermodynamic calculations, which can be either ``"doped"`` or ``"py-sc-fermi"``. ``"py-sc-fermi"`` @@ -209,12 +225,33 @@ def __init__( Whether to skip the warning about the DOS VBM differing from ``DefectThermodynamics.vbm`` by >0.05 eV. Should only be used when the reason for this difference is known/acceptable. Default is ``False``. + + Key attributes: + defect_thermodynamics (DefectThermodynamics): + The ``DefectThermodynamics`` object used for the thermodynamic + calculations. + backend (str): + The code backend used for the thermodynamic calculations (``"doped"`` + or ``"py-sc-fermi"``). + volume (float): + Volume of the unit cell in the bulk DOS calculation (stored in + ``self.defect_thermodynamics.bulk_dos``). + skip_check (bool): + Whether to skip the warning about the DOS VBM differing from the defect + entries VBM by >0.05 eV. Should only be used when the reason for this + difference is known/acceptable. + multiplicity_scaling (int): + Scaling factor to account for the difference in volume between the + defect supercells and the bulk DOS calculation cell, when using the + ``py-sc-fermi`` backend. + py_sc_fermi_dos (DOS): + A ``py-sc-fermi`` ``DOS`` object, generated from the input + ``FermiDos`` object, for use with the ``py-sc-fermi`` backend. """ # Note: In theory multiprocessing could be introduced to make thermodynamic calculations # over large grids faster, but with the current code there seems to be issues with thread # locking / synchronisation (which shouldn't be necessary...). Worth keeping in mind if # needed in future. - # TODO: List class attributes in docstring self.defect_thermodynamics = defect_thermodynamics self.skip_check = skip_check if bulk_dos is not None: @@ -250,11 +287,11 @@ def __init__( # Parse chemical potentials, either using input values (after formatting them in the doped format) # or using the class attributes if set: - self.chempots, self.el_refs = _parse_chempots( + self.defect_thermodynamics._chempots, self.defect_thermodynamics._el_refs = _parse_chempots( chempots or self.defect_thermodynamics.chempots, el_refs or self.defect_thermodynamics.el_refs ) - if self.chempots is None: + if self.defect_thermodynamics.chempots is None: raise ValueError( "You must supply a chemical potentials dictionary or have them present in the " "DefectThermodynamics object." @@ -313,9 +350,8 @@ def _check_required_backend_and_error(self, required_backend: str): if required_backend.lower() == "doped": if not isinstance(self.defect_thermodynamics.bulk_dos, FermiDos): raise_error = True - else: - if self._DOS is None: - raise_error = True + elif self._DOS is None: + raise_error = True if raise_error: raise RuntimeError( @@ -325,19 +361,33 @@ def _check_required_backend_and_error(self, required_backend: str): def _get_fermi_level_and_carriers( self, - chempots: dict[str, float], + single_chempot_dict: dict[str, float], + el_refs: Optional[dict[str, float]] = None, temperature: float = 300, effective_dopant_concentration: Optional[float] = None, ) -> tuple[float, float, float]: """ - Calculate the Fermi level and carrier concentrations under a given - chemical potential regime and temperature. + Calculate the equilibrium Fermi level and carrier concentrations under + a given chemical potential regime and temperature. Args: - chempots (dict[str, float]): - A dictionary of chemical potentials for the elements. - The keys are element symbols, and the values are their - corresponding chemical potentials. + single_chempot_dict (dict[str, float]): + Dictionary of chemical potentials to use for calculating the equilibrium + Fermi level position. Here, this 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 + reference energies) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + el_refs (dict[str, float]): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}``. Unnecessary if + ``self.defect_thermodynamics.el_refs`` is set (i.e. if ``chempots`` was + provided to ``self.defect_thermodynamics`` in the format generated by + ``doped``). + (Default: None) temperature (float): The temperature at which to solve for the Fermi level and carrier concentrations, in Kelvin. @@ -360,24 +410,50 @@ def _get_fermi_level_and_carriers( """ self._check_required_backend_and_error("doped") fermi_level, electrons, holes = self.defect_thermodynamics.get_equilibrium_fermi_level( # type: ignore - chempots=chempots, - limit=None, + chempots=single_chempot_dict, + el_refs=el_refs, temperature=temperature, return_concs=True, effective_dopant_concentration=effective_dopant_concentration, ) # use already-set bulk dos return fermi_level, electrons, holes - def _get_limits(self, limit): - limit_dict = get_rich_poor_limit_dict(self.chempots) - limit = next(v for k, v in limit_dict.items() if k == limit) - return self.chempots["limits_wrt_el_refs"][limit] + def _get_single_chempot_dict( + self, limit: Optional[str] = None, chempots: Optional[dict] = None, el_refs: Optional[dict] = None + ) -> tuple[dict[str, float], Any]: + """ + Get the chemical potentials for a single limit (``limit``) from the + ``chempots`` (or ``self.defect_thermodynamics.chempots``) dictionary. + + Returns a `single` chemical potential dictionary for the specified limit. + """ + chempots, el_refs = self.defect_thermodynamics._get_chempots( + chempots, el_refs + ) # returns self.defect_thermodynamics.chempots if chempots is None + if chempots is None: + raise ValueError( + "No chemical potentials supplied or present in self.defect_thermodynamics.chempots!" + ) + limit = _parse_limit(chempots, limit) + + if ( + limit not in chempots["limits_wrt_el_refs"] + and "User Chemical Potentials" not in chempots["limits_wrt_el_refs"] + ): + raise ValueError( + f"Limit '{limit}' not found in the chemical potentials dictionary! You must specify " + f"an appropriate limit or provide an appropriate chempots dictionary." + ) + + return chempots["limits_wrt_el_refs"][limit or "User Chemical Potentials"], el_refs def equilibrium_solve( self, - chempots: dict[str, float], + single_chempot_dict: dict[str, float], + el_refs: Optional[dict[str, float]] = None, temperature: float = 300, effective_dopant_concentration: Optional[float] = None, + append_chempots: bool = True, ) -> pd.DataFrame: """ Calculate the Fermi level and defect/carrier concentrations under @@ -387,12 +463,27 @@ def equilibrium_solve( Typically not intended for direct usage, as the same functionality is provided by ``DefectThermodynamics.get_equilibrium_concentrations/fermi_level()``. + Rather this is used internally to provide a unified interface for both + backends, within the ``scan_...`` functions. Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the elements. - The keys are element symbols, and the values are their corresponding - chemical potentials. + single_chempot_dict (dict[str, float]): + Dictionary of chemical potentials to use for calculating the equilibrium + Fermi level position and defect/carrier concentrations. Here, this + 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 + reference energies) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + el_refs (dict[str, float]): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}``. Unnecessary if + ``self.defect_thermodynamics.el_refs`` is set (i.e. if ``chempots`` was + provided to ``self.defect_thermodynamics`` in the format generated by + ``doped``). + (Default: None) temperature (float): The temperature at which to solve for defect concentrations, in Kelvin. Defaults to 300 K. @@ -405,12 +496,16 @@ def equilibrium_solve( value corresponds to acceptor doping. Defaults to ``None``, corresponding to no additional extrinsic dopant. + append_chempots (bool): + Whether to append the chemical potentials (and effective dopant + concentration, if provided) to the results ``DataFrame``. + Default is ``True``. Returns: # TODO: Check this output format matches for both backends! pd.DataFrame: A DataFrame containing the defect and carrier concentrations, as well - as the self-consistent Fermi energy and additional properties such as + as the self-consistent Fermi level and additional properties such as electron and hole concentrations. The columns include: - "Fermi Level": The self-consistent Fermi level in eV. - "Electrons (cm^-3)": The electron concentration. @@ -420,15 +515,22 @@ def equilibrium_solve( - "Defect": The defect type. - "Concentration (cm^-3)": The concentration of the defect in cm^-3. Additional columns may include concentrations for specific defects and carriers. + + Returns: + pd.DataFrame: A DataFrame containing the calculated defect and carrier concentrations, + along with the self-consistent Fermi level. The DataFrame also includes the provided + chemical potentials as additional columns. """ if self.backend == "doped": fermi_level, electrons, holes = self._get_fermi_level_and_carriers( - chempots=chempots, + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, temperature=temperature, effective_dopant_concentration=effective_dopant_concentration, ) concentrations = self.defect_thermodynamics.get_equilibrium_concentrations( - chempots=chempots, + chempots=single_chempot_dict, + el_refs=el_refs, fermi_level=fermi_level, temperature=temperature, per_charge=False, @@ -449,43 +551,53 @@ def equilibrium_solve( excluded_columns = ["Defect", "Charge", "Charge State Population"] for column in concentrations.columns.difference(excluded_columns): concentrations[column] = concentrations[column].astype(float) - return concentrations - # else py-sc-fermi backend: - defect_system = self._generate_defect_system( - chempots=chempots, - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - ) + else: # py-sc-fermi backend: + defect_system = self._generate_defect_system( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) - with np.errstate(all="ignore"): - conc_dict = defect_system.concentration_dict() - data = [] - - for k, v in conc_dict.items(): - if k not in ["Fermi Energy", "n0", "p0", "Dopant"]: - row = { - "Temperature": defect_system.temperature, - "Fermi Level": conc_dict["Fermi Energy"], - "Holes (cm^-3)": conc_dict["p0"], - "Electrons (cm^-3)": conc_dict["n0"], - } - if "Dopant" in conc_dict: - row["Dopant (cm^-3)"] = conc_dict["Dopant"] - row.update({"Defect": k, "Concentration (cm^-3)": v}) - data.append(row) - - results_df = pd.DataFrame(data) - return results_df.set_index("Defect", drop=True) + with np.errstate(all="ignore"): + conc_dict = defect_system.concentration_dict() + data = [] + + for k, v in conc_dict.items(): + if k not in ["Fermi Energy", "n0", "p0", "Dopant"]: + row = { + "Temperature": defect_system.temperature, + "Fermi Level": conc_dict["Fermi Energy"], + "Holes (cm^-3)": conc_dict["p0"], + "Electrons (cm^-3)": conc_dict["n0"], + } + if "Dopant" in conc_dict: + row["Dopant (cm^-3)"] = conc_dict["Dopant"] + row.update({"Defect": k, "Concentration (cm^-3)": v}) + data.append(row) + + 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 def pseudo_equilibrium_solve( self, - chempots: dict[str, float], annealing_temperature: float, + single_chempot_dict: dict[str, float], + el_refs: Optional[dict[str, float]] = None, quenched_temperature: float = 300, effective_dopant_concentration: Optional[float] = None, fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, + append_chempots: bool = True, ) -> pd.DataFrame: """ Calculate the self-consistent Fermi level and corresponding @@ -496,6 +608,8 @@ def pseudo_equilibrium_solve( Typically not intended for direct usage, as the same functionality is provided by ``DefectThermodynamics.get_quenched_fermi_level_and_concentrations()``. + Rather this is used internally to provide a unified interface for both + backends, within the ``scan_...`` functions. 'Pseudo-equilibrium' here refers to the use of frozen defect and dilute limit approximations under the constraint of charge neutrality (i.e. @@ -542,16 +656,30 @@ def pseudo_equilibrium_solve( ``defect_entry.degeneracy_factors`` attributes. Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the elements. - The keys are element symbols, and the values are their corresponding - chemical potentials. annealing_temperature (float): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to the highest temperature during annealing/synthesis of the material (at which we assume equilibrium defect concentrations) within the frozen defect approach. + single_chempot_dict (dict[str, float]): + Dictionary of chemical potentials to use for calculating the + pseudo-equilibrium Fermi level position and defect/carrier + concentrations. Here, this 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 + reference energies) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + el_refs (dict[str, float]): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}``. Unnecessary if + ``self.defect_thermodynamics.el_refs`` is set (i.e. if ``chempots`` was + provided to ``self.defect_thermodynamics`` in the format generated by + ``doped``). + (Default: None) quenched_temperature (float): Temperature in Kelvin at which to calculate the self-consistent (constrained equilibrium) Fermi level and carrier concentrations, @@ -577,6 +705,10 @@ def pseudo_equilibrium_solve( fixing, useful for highly mobile defects that are not expected to be "frozen-in" upon quenching. Defaults to ``None``. # TODO: Allow matching of substring + append_chempots (bool): + Whether to append the chemical potentials (and effective dopant + concentration, if provided) to the results ``DataFrame``. + Default is ``True``. Returns: # TODO: Check output format for both backends! pd.DataFrame: @@ -584,7 +716,7 @@ def pseudo_equilibrium_solve( under pseudo-equilibrium conditions, along with the self-consistent Fermi level. The columns include: - - "Fermi Level": The self-consistent Fermi energy. + - "Fermi Level": The self-consistent Fermi level. - "Electrons (cm^-3)": The electron concentration. - "Holes (cm^-3)": The hole concentration. - "Annealing Temperature": The annealing temperature. @@ -595,6 +727,11 @@ def pseudo_equilibrium_solve( Additional columns may include concentrations for specific defects and other relevant data. + + Returns: + pd.DataFrame: A DataFrame containing the calculated defect and carrier concentrations + using the pseudo-equilibrium approach, along with the provided chemical potentials + as additional columns. """ py_sc_fermi_required = fix_charge_states or free_defects is not None if py_sc_fermi_required and self._DOS is None: @@ -610,8 +747,8 @@ def pseudo_equilibrium_solve( holes, concentrations, ) = self.defect_thermodynamics.get_quenched_fermi_level_and_concentrations( # type: ignore - chempots=chempots, - limit=None, + chempots=single_chempot_dict, + el_refs=el_refs, annealing_temperature=annealing_temperature, quenched_temperature=quenched_temperature, effective_dopant_concentration=effective_dopant_concentration, @@ -641,42 +778,51 @@ def pseudo_equilibrium_solve( for column in concentrations.columns.difference(excluded_columns): concentrations[column] = concentrations[column].astype(float) - renamed_concentrations = trimmed_concentrations_sub_duplicates.rename( + concentrations = trimmed_concentrations_sub_duplicates.rename( columns={"Total Concentration (cm^-3)": "Concentration (cm^-3)"}, ) - return renamed_concentrations.set_index("Defect", drop=True) + concentrations = concentrations.set_index("Defect", drop=True) - # else py-sc-fermi: - defect_system = self._generate_annealed_defect_system( - chempots=chempots, - annealing_temperature=annealing_temperature, - quenched_temperature=quenched_temperature, - effective_dopant_concentration=effective_dopant_concentration, - free_defects=free_defects, - fix_charge_states=fix_charge_states, - ) + else: # py-sc-fermi + defect_system = self._generate_annealed_defect_system( + annealing_temperature=annealing_temperature, + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, + quenched_temperature=quenched_temperature, + effective_dopant_concentration=effective_dopant_concentration, + free_defects=free_defects, + fix_charge_states=fix_charge_states, + ) - with np.errstate(all="ignore"): - conc_dict = defect_system.concentration_dict() - - data = [] - for k, v in conc_dict.items(): - if k not in ["Fermi Energy", "n0", "p0"]: - row = { - "Annealing Temperature": annealing_temperature, - "Quenched Temperature": quenched_temperature, - "Fermi Level": conc_dict["Fermi Energy"], - "Holes (cm^-3)": conc_dict["p0"], - "Electrons (cm^-3)": conc_dict["n0"], - } - row.update({"Defect": k, "Concentration (cm^-3)": v}) - if "Dopant" in conc_dict: - row["Dopant (cm^-3)"] = conc_dict["Dopant"] - row.update({"Defect": k, "Concentration (cm^-3)": v}) - data.append(row) - - results_df = pd.DataFrame(data) - return results_df.set_index("Defect", drop=True) + with np.errstate(all="ignore"): + conc_dict = defect_system.concentration_dict() + + data = [] + for k, v in conc_dict.items(): + if k not in ["Fermi Energy", "n0", "p0"]: + row = { + "Annealing Temperature": annealing_temperature, + "Quenched Temperature": quenched_temperature, + "Fermi Level": conc_dict["Fermi Energy"], + "Holes (cm^-3)": conc_dict["p0"], + "Electrons (cm^-3)": conc_dict["n0"], + } + row.update({"Defect": k, "Concentration (cm^-3)": v}) + if "Dopant" in conc_dict: + row["Dopant (cm^-3)"] = conc_dict["Dopant"] + row.update({"Defect": k, "Concentration (cm^-3)": v}) + data.append(row) + + 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 def scan_temperature( self, @@ -685,13 +831,14 @@ def scan_temperature( temperature_range: Union[float, list[float]] = 300, chempots: Optional[dict[str, float]] = None, limit: Optional[str] = None, + el_refs: Optional[dict[str, float]] = None, effective_dopant_concentration: Optional[float] = None, fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: - """ + r""" Scan over a range of temperatures and solve for the defect - concentrations, carrier concentrations, and Fermi energy at each + concentrations, carrier concentrations, and Fermi level at each temperature / annealing & quenched temperature pair. If ``annealing_temperature_range`` (and ``quenched_temperature_range``; @@ -722,21 +869,56 @@ def scan_temperature( Temperature range to solve over, under thermodynamic equilibrium (if ``annealing_temperature_range`` is not specified). Defaults to just 300 K. - chempots (Optional[dict[str, float]]): Dictionary of chemical potentials to use for - calculating the defect formation energies (and thus concentrations and Fermi level). - This can be a dictionary of chemical potentials for a single limit (limit), - in the format: {element symbol: chemical potential}. If manually specifying chemical - potentials this way, you can set the `el_refs` option with the DFT reference energies - of the elemental phases, in which case it is the formal chemical potentials - (i.e., relative to the elemental references) that should be given here, - otherwise the absolute (DFT) chemical potentials should be given. - limit (Optional[str]): The chemical potential limit for which to determine the equilibrium - Fermi level. Can be either: - - `None`, if `chempots` corresponds to a single chemical potential limit - otherwise - will use the first chemical potential limit in the `chempots` dict. - - `"X-rich"/"X-poor"` where X is an element in the system, in which case the most - X-rich/poor limit will be used (e.g., "Li-rich"). - - A key in the `chempots["limits"]` dictionary. + chempots (dict): + Dictionary of chemical potentials to use for calculating the defect + formation energies (and thus concentrations). If ``None`` (default), + will use ``self.defect_thermodynamics.chempots`` (= 0 for all + chemical potentials by default). + This can have the form of ``{"limits": [{'limit': [chempot_dict]}], ...}`` + (the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials)) and specific limits (chemical potential + limits) can then be chosen using ``limit``. + + Alternatively this can be a dictionary of chemical potentials for a + single limit (``limit``), in the format: + ``{element symbol: chemical potential}``. + If manually specifying chemical potentials this way, you can set the + ``el_refs`` option with the DFT reference energies of the elemental phases, + in which case it is the formal chemical potentials (i.e. relative to the + elemental references) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + + Note that you can also set ``FermiSolver.defect_thermodynamics.chempots = ...`` + or ``DefectThermodynamics.chempots = ...`` (with the same input options) + to set the default chemical potentials for all calculations. + (Default: None) + limit (str): + The chemical potential limit for which to calculate the Fermi level + positions and defect/carrier concentrations. Can be either: + + - ``None``, if ``chempots`` corresponds to a single chemical potential + limit - otherwise will use the first chemical potential limit in the + ``chempots`` dict. + - ``"X-rich"/"X-poor"`` where X is an element in the system, in which + case the most X-rich/poor limit will be used (e.g. "Li-rich"). + - A key in the ``(self.defect_thermodynamics.)chempots["limits"]`` + dictionary. + + The latter two options can only be used if ``chempots`` is in the + ``doped`` format (see chemical potentials tutorial). + (Default: None) + el_refs (dict): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}`` (to determine the formal chemical + potentials, when ``chempots`` has been manually specified as + ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is + provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``FermiSolver.defect_thermodynamics.el_refs = ...`` + or ``DefectThermodynamics.el_refs = ...`` (with the same input options) + to set the default elemental reference energies for all calculations. + (Default: None) effective_dopant_concentration (Optional[float]): The fixed concentration (in cm^-3) of an arbitrary dopant or impurity in the material. This value is included in the charge @@ -767,16 +949,14 @@ def scan_temperature( if isinstance(quenched_temperature_range, float): quenched_temperature_range = [quenched_temperature_range] - if chempots is None and limit is not None: - chempots = self._get_limits(limit) - elif chempots is None: - raise ValueError("You must specify a limit or chempots dictionary.") + single_chempot_dict, el_refs = self._get_single_chempot_dict(limit, chempots, el_refs) if annealing_temperature_range is not None: - all_data_df = pd.concat( + return pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots, + self.pseudo_equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, quenched_temperature=quench_temp, annealing_temperature=anneal_temp, effective_dopant_concentration=effective_dopant_concentration, @@ -789,33 +969,33 @@ def scan_temperature( ] ) - else: - all_data_df = pd.concat( - [ - self._solve_and_append_chempots( - chempots=chempots, - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - ) - for temperature in tqdm(temperature_range) - ] - ) - - return all_data_df + return pd.concat( + [ + self.equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) + for temperature in tqdm(temperature_range) + ] + ) - def scan_chempots( + def scan_dopant_concentration( self, - chempots: list[dict[str, float]], + effective_dopant_concentration_range: Union[float, list[float]], + chempots: Optional[dict[str, float]] = None, + limit: Optional[str] = None, + el_refs: Optional[dict[str, float]] = None, annealing_temperature: Optional[float] = None, quenched_temperature: float = 300, temperature: float = 300, - effective_dopant_concentration: Optional[float] = None, fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: - """ - Scan over a range of chemical potentials and solve for the defect - concentrations and Fermi energy at each set of chemical potentials. + r""" + Calculate the defect concentrations under a range of effective dopant + concentrations. If ``annealing_temperature`` (and ``quenched_temperature``; 300 K by default) are specified, then the frozen defect approximation is @@ -828,10 +1008,60 @@ def scan_chempots( assuming thermodynamic equilibrium at that temperature. Args: - chempots (list[dict[str, float]]): - A list of dictionaries where each dictionary represents a set - of chemical potentials to scan over. The keys are element symbols, - and the values are their corresponding chemical potentials. + effective_dopant_concentration_range (Union[float, list[float]]): + The range of effective dopant concentrations to solver over. + This can be a single value or a list of values representing + different concentrations. + chempots (dict): + Dictionary of chemical potentials to use for calculating the defect + formation energies (and thus concentrations). If ``None`` (default), + will use ``self.defect_thermodynamics.chempots`` (= 0 for all + chemical potentials by default). + This can have the form of ``{"limits": [{'limit': [chempot_dict]}], ...}`` + (the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials)) and specific limits (chemical potential + limits) can then be chosen using ``limit``. + + Alternatively this can be a dictionary of chemical potentials for a + single limit (``limit``), in the format: + ``{element symbol: chemical potential}``. + If manually specifying chemical potentials this way, you can set the + ``el_refs`` option with the DFT reference energies of the elemental phases, + in which case it is the formal chemical potentials (i.e. relative to the + elemental references) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + + Note that you can also set ``FermiSolver.defect_thermodynamics.chempots = ...`` + or ``DefectThermodynamics.chempots = ...`` (with the same input options) + to set the default chemical potentials for all calculations. + (Default: None) + limit (str): + The chemical potential limit for which to calculate the Fermi level + positions and defect/carrier concentrations. Can be either: + + - ``None``, if ``chempots`` corresponds to a single chemical potential + limit - otherwise will use the first chemical potential limit in the + ``chempots`` dict. + - ``"X-rich"/"X-poor"`` where X is an element in the system, in which + case the most X-rich/poor limit will be used (e.g. "Li-rich"). + - A key in the ``(self.defect_thermodynamics.)chempots["limits"]`` + dictionary. + + The latter two options can only be used if ``chempots`` is in the + ``doped`` format (see chemical potentials tutorial). + (Default: None) + el_refs (dict): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}`` (to determine the formal chemical + potentials, when ``chempots`` has been manually specified as + ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is + provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``FermiSolver.defect_thermodynamics.el_refs = ...`` + or ``DefectThermodynamics.el_refs = ...`` (with the same input options) + to set the default elemental reference energies for all calculations. + (Default: None) annealing_temperature (Optional[float]): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to @@ -847,18 +1077,9 @@ def scan_chempots( Default is 300 K. temperature (float): The temperature at which to solve for defect concentrations - and Fermi energy, under thermodynamic equilibrium (if + and Fermi level, under thermodynamic equilibrium (if ``annealing_temperature`` is not specified). Defaults to 300 K. - effective_dopant_concentration (Optional[float]): - The fixed concentration (in cm^-3) of an arbitrary dopant or - impurity in the material. This value is included in the charge - neutrality condition to analyze the Fermi level and doping - response under hypothetical doping conditions. - A positive value corresponds to donor doping, while a negative - value corresponds to acceptor doping. - Defaults to ``None``, corresponding to no additional extrinsic - dopant. fix_charge_states (bool): Whether to fix the concentrations of individual defect charge states (``True``) or allow charge states to vary while keeping total defect @@ -870,51 +1091,62 @@ def scan_chempots( to be "frozen-in" upon quenching. Defaults to ``None``. Returns: - pd.DataFrame: A DataFrame containing defect and carrier concentrations for each set - of chemical potentials. Each row corresponds to a different set of chemical potentials. + pd.DataFrame: + A ``DataFrame`` containing the defect and carrier concentrations for each + effective dopant concentration. Each row represents the concentrations + for a different dopant concentration. """ - # TODO: This code needs to be either fixed or cut + if isinstance(effective_dopant_concentration_range, float): + effective_dopant_concentration_range = [effective_dopant_concentration_range] + + single_chempot_dict, el_refs = self._get_single_chempot_dict(limit, chempots, el_refs) + if annealing_temperature is not None: - all_data_df = pd.concat( + return pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots, + self.pseudo_equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, quenched_temperature=quenched_temperature, annealing_temperature=annealing_temperature, effective_dopant_concentration=effective_dopant_concentration, fix_charge_states=fix_charge_states, free_defects=free_defects, ) + for effective_dopant_concentration in tqdm(effective_dopant_concentration_range) ] ) - else: - all_data_df = pd.concat( - [ - self._solve_and_append_chempots( - chempots=chempots, - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - ) - ] - ) - - return all_data_df + return pd.concat( + [ + self.equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) + for effective_dopant_concentration in tqdm(effective_dopant_concentration_range) + ] + ) - def scan_dopant_concentration( + def interpolate_chempots( self, - effective_dopant_concentration_range: Union[float, list[float]], - chempots: Optional[dict[str, float]] = None, - limit: Optional[str] = None, + n_points: int, + chempots: Optional[Union[list[dict], dict]] = None, + limits: Optional[list[str]] = None, + el_refs: Optional[dict[str, float]] = None, annealing_temperature: Optional[float] = None, quenched_temperature: float = 300, temperature: float = 300, + effective_dopant_concentration: Optional[float] = None, fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: """ - Calculate the defect concentrations under a range of effective dopant - concentrations. + Interpolate between two sets of chemical potentials and solve for the + defect concentrations and Fermi level at each interpolated point. + Chemical potentials can be interpolated between two sets of chemical + potentials or between two limits. If ``annealing_temperature`` (and ``quenched_temperature``; 300 K by default) are specified, then the frozen defect approximation is @@ -927,18 +1159,53 @@ def scan_dopant_concentration( assuming thermodynamic equilibrium at that temperature. Args: - effective_dopant_concentration_range (Union[float, list[float]]): - The range of effective dopant concentrations to solver over. - This can be a single value or a list of values representing - different concentrations. - chempots (Optional[dict[str, float]]): - A dictionary of chemical potentials for the elements, where - the keys are element symbols and the values are the - corresponding chemical potentials. - limit (Optional[str]): - The chemical potential limit to use for calculations. This can - specify a particular limit such as "X-rich" or "X-poor", or can - be a key from the chemical potentials dictionary. + n_points (int): + The number of points to generate between chemical potential + end points. + chempots (Optional[list[dict]]): + The chemical potentials to interpolate between. This can be + either a list containing two dictionaries, each representing + a set of chemical potentials for a single limit (in the format: + ``{element symbol: chemical potential}``) to interpolate between, + or can be a single chemical potentials dictionary in the + ``doped`` format (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``) + -- in which case ``limits`` must be specified to pick the + end-points to interpolate between. + + If ``None`` (default), will use ``self.defect_thermodynamics.chempots``. + Note that you can also set + ``FermiSolver.defect_thermodynamics.chempots = ...`` + or ``DefectThermodynamics.chempots = ...`` (with the same input + options) to set the default chemical potentials for all calculations. + + If manually specifying chemical potentials with a list of two + dictionaries, you can also set the ``el_refs`` option with the DFT + reference energies of the elemental phases if desired, in which case + it is the formal chemical potentials (i.e. relative to the elemental + references) that should be given here, otherwise the absolute (DFT) + chemical potentials should be given. + limits (Optional[list[str]]): + The chemical potential limits to interpolate between, as a list + containing two strings. Each string should be in the format + ``"X-rich"/"X-poor"``, where X is an element in the system, or + be a key in the ``(self.defect_thermodynamics.)chempots["limits"]`` + dictionary. + + If not provided, ``chempots`` must be specified as a list of two + single chemical potential dictionaries for single limits. + el_refs (dict): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}`` (to determine the formal chemical + potentials, when ``chempots`` has been manually specified as + ``[{element symbol: chemical potential}, ...]``). Unnecessary if + ``chempots`` is provided/present in format generated by ``doped`` + (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``). + + Note that you can also set ``FermiSolver.defect_thermodynamics.el_refs = ...`` + or ``DefectThermodynamics.el_refs = ...`` (with the same input options) + to set the default elemental reference energies for all calculations. + (Default: None) annealing_temperature (Optional[float]): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to @@ -954,9 +1221,18 @@ def scan_dopant_concentration( Default is 300 K. temperature (float): The temperature at which to solve for defect concentrations - and Fermi energy, under thermodynamic equilibrium (if + and Fermi level, under thermodynamic equilibrium (if ``annealing_temperature`` is not specified). Defaults to 300 K. + effective_dopant_concentration (Optional[float]): + The fixed concentration (in cm^-3) of an arbitrary dopant or + impurity in the material. This value is included in the charge + neutrality condition to analyze the Fermi level and doping + response under hypothetical doping conditions. + A positive value corresponds to donor doping, while a negative + value corresponds to acceptor doping. + Defaults to ``None``, corresponding to no additional extrinsic + dopant. fix_charge_states (bool): Whether to fix the concentrations of individual defect charge states (``True``) or allow charge states to vary while keeping total defect @@ -968,54 +1244,117 @@ def scan_dopant_concentration( to be "frozen-in" upon quenching. Defaults to ``None``. Returns: - pd.DataFrame: A DataFrame containing the defect and carrier concentrations for each - effective dopant concentration. Each row represents the concentrations for a different - dopant concentration. + pd.DataFrame: + A ``DataFrame`` containing the defect and carrier concentrations for + each interpolated set of chemical potentials. Each row represents the + concentrations for a different interpolated point. """ - if isinstance(effective_dopant_concentration_range, float): - effective_dopant_concentration_range = [effective_dopant_concentration_range] + if isinstance(chempots, list): # should be two single chempot dictionaries + if len(chempots) != 2: + raise ValueError( + f"If `chempots` is a list, it must contain two dictionaries representing the starting " + f"and ending chemical potentials. The provided list has {len(chempots)} entries!" + ) + single_chempot_dict_1, single_chempot_dict_2 = chempots - if chempots is None and limit is not None: - chempots = self._get_limits(limit) - elif chempots is None: - raise ValueError("You must specify a limit or chempots dictionary.") + else: # should be a dictionary in the ``doped`` format or ``None``: + chempots, el_refs = self.defect_thermodynamics._get_chempots( + chempots, el_refs + ) # returns self.defect_thermodynamics.chempots if chempots is None + if chempots is None: + raise ValueError( + "No chemical potentials supplied or present in self.defect_thermodynamics.chempots!" + ) + + if limits is None or len(limits) != 2: + raise ValueError( + f"If `chempots` is not provided as a list, then `limits` must be a list containing " + f"two strings representing the chemical potential limits to interpolate between. " + f"The provided `limits` is: {limits}." + ) + + assert isinstance(chempots, dict) # typing + single_chempot_dict_1, el_refs = self._get_single_chempot_dict(limits[0], chempots, el_refs) + single_chempot_dict_2, el_refs = self._get_single_chempot_dict(limits[1], chempots, el_refs) + + interpolated_chempots = self._get_interpolated_chempots( + single_chempot_dict_1, single_chempot_dict_2, n_points + ) if annealing_temperature is not None: - all_data_df = pd.concat( + return pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots, + self.pseudo_equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, quenched_temperature=quenched_temperature, annealing_temperature=annealing_temperature, effective_dopant_concentration=effective_dopant_concentration, fix_charge_states=fix_charge_states, free_defects=free_defects, ) - for effective_dopant_concentration in tqdm(effective_dopant_concentration_range) + for single_chempot_dict in tqdm(interpolated_chempots) ] ) - else: - all_data_df = pd.concat( - [ - self._solve_and_append_chempots( - chempots=chempots, - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - fix_charge_states=fix_charge_states, - free_defects=free_defects, - ) - for effective_dopant_concentration in tqdm(effective_dopant_concentration_range) - ] - ) - - return all_data_df + return pd.concat( + [ + self.equilibrium_solve( + single_chempot_dict=single_chempot_dict, + el_refs=el_refs, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) + for single_chempot_dict in tqdm(interpolated_chempots) + ] + ) - def interpolate_chempots( + def _get_interpolated_chempots( self, + chempot_start: dict, + chempot_end: dict, n_points: int, - chempots: Optional[list[dict]] = None, + ) -> list: + """ + Generate a set of interpolated chemical potentials between two points. + + Here, these should be dictionaries of chemical potentials for `single` + limits, in the format: ``{element symbol: chemical potential}``. + If ``el_refs`` is provided (to the parent function) or set in + ``self.defect_thermodynamics.el_refs``, then it is the formal chemical + potentials (i.e. relative to the elemental reference energies) that + should be used here, otherwise the absolute (DFT) chemical potentials + should be used. + + Args: + chempot_start (dict): + A dictionary representing the starting chemical potentials. + chempot_end (dict): + A dictionary representing the ending chemical potentials. + n_points (int): + The number of interpolated points to generate, `including` + the start and end points. + + Returns: + list: + A list of dictionaries, where each dictionary contains a `single` + set of interpolated chemical potentials. The length of the list + corresponds to `n_points`, and each dictionary corresponds to an + interpolated state between the starting and ending chemical potentials. + """ + return [ + { + key: chempot_start[key] + (chempot_end[key] - chempot_start[key]) * i / (n_points - 1) + for key in chempot_start + } + for i in range(n_points) + ] + + def scan_chempots( + self, + chempots: Union[list[dict[str, float]], dict[str, dict]], limits: Optional[list[str]] = None, + el_refs: Optional[dict[str, float]] = None, annealing_temperature: Optional[float] = None, quenched_temperature: float = 300, temperature: float = 300, @@ -1024,10 +1363,8 @@ def interpolate_chempots( free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: """ - Interpolate between two sets of chemical potentials and solve for the - defect concentrations and Fermi energy at each interpolated point. - Chemical potentials can be interpolated between two sets of chemical - potentials or between two limits. + Scan over a range of chemical potentials and solve for the defect + concentrations and Fermi level at each set of chemical potentials. If ``annealing_temperature`` (and ``quenched_temperature``; 300 K by default) are specified, then the frozen defect approximation is @@ -1040,17 +1377,51 @@ def interpolate_chempots( assuming thermodynamic equilibrium at that temperature. Args: - n_points (int): - The number of points to generate between chemical potential - end points. - chempots (Optional[list[dict]]): - A list containing two dictionaries, each representing a set - of chemical potentials to interpolate between. If not provided, - ``limits`` must be specified. + chempots (Optional[Union[list[dict], dict]]): + The chemical potentials to scan over. This can be either a list + containing dictionaries of a set of chemical potentials for a + `single` limit (in the format: ``{element symbol: chemical potential}``), + or can be a single chemical potentials dictionary in the + ``doped`` format (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``) + -- in which case ``limits`` can be specified to pick the chemical + potential limits to scan over (otherwise scans over all limits in the + ``chempots`` dictionary). + + If ``None`` (default), will use ``self.defect_thermodynamics.chempots``. + Note that you can also set + ``FermiSolver.defect_thermodynamics.chempots = ...`` + or ``DefectThermodynamics.chempots = ...`` (with the same input + options) to set the default chemical potentials for all calculations. + + If manually specifying chemical potentials with a list of dictionaries, + you can also set the ``el_refs`` option with the DFT reference energies + of the elemental phases if desired, in which case it is the formal + chemical potentials (i.e. relative to the elemental references) that + should be given here, otherwise the absolute (DFT) chemical potentials + should be given. limits (Optional[list[str]]): - A list containing two strings, each representing a chemical - potential limit (e.g., "X-rich" or "X-poor") to interpolate - between. If not provided, ``chempots`` must be specified. + The chemical potential limits to scan over, as a list of strings, if + ``chempots`` was provided / is present in the ``doped`` format. Each + string should be in the format ``"X-rich"/"X-poor"``, where X is an + element in the system, or be a key in the + ``(self.defect_thermodynamics.)chempots["limits"]`` dictionary. + + If ``None`` (default) and ``chempots`` is in the ``doped`` format + (rather than a list of single chemical potential limits), will scan + over all limits in the ``chempots`` dictionary. + el_refs (dict): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}`` (to determine the formal chemical + potentials, when ``chempots`` has been manually specified as + ``[{element symbol: chemical potential}, ...]``). Unnecessary if + ``chempots`` is provided/present in format generated by ``doped`` + (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``). + + Note that you can also set ``FermiSolver.defect_thermodynamics.el_refs = ...`` + or ``DefectThermodynamics.el_refs = ...`` (with the same input options) + to set the default elemental reference energies for all calculations. + (Default: None) annealing_temperature (Optional[float]): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to @@ -1066,7 +1437,7 @@ def interpolate_chempots( Default is 300 K. temperature (float): The temperature at which to solve for defect concentrations - and Fermi energy, under thermodynamic equilibrium (if + and Fermi level, under thermodynamic equilibrium (if ``annealing_temperature`` is not specified). Defaults to 300 K. effective_dopant_concentration (Optional[float]): @@ -1089,153 +1460,61 @@ def interpolate_chempots( to be "frozen-in" upon quenching. Defaults to ``None``. Returns: - pd.DataFrame: A DataFrame containing the defect and carrier concentrations for each - interpolated set of chemical potentials. Each row represents the concentrations for - a different interpolated point. + pd.DataFrame: + A ``DataFrame`` containing defect and carrier concentrations for each + set of chemical potentials. Each row corresponds to a different set + of chemical potentials. """ - if chempots is None and limits is not None: - chempots_1 = self._get_limits(limits[0]) - chempots_2 = self._get_limits(limits[1]) - elif chempots is not None: - chempots_1 = chempots[0] - chempots_2 = chempots[1] + if isinstance(chempots, dict): # should be a dictionary in the ``doped`` format or ``None``: + chempots, el_refs = self.defect_thermodynamics._get_chempots( + chempots, el_refs + ) # returns self.defect_thermodynamics.chempots if chempots is None + if chempots is None: + raise ValueError( + "No chemical potentials supplied or present in self.defect_thermodynamics.chempots!" + ) + assert isinstance(chempots, dict) # typing + chempots, el_refs = self._parse_and_check_grid_like_chempots(chempots) - interpolated_chem_pots = self._get_interpolated_chempots(chempots_1, chempots_2, n_points) + if limits is None: + limits = list(chempots["limits_wrt_el_refs"].keys()) + elif not isinstance(limits, list): + raise ValueError( + "`limits` must be either a list of limits (as strings) or `None` for `scan_chempots`!" + ) + + chempots = [self._get_single_chempot_dict(limit, chempots, el_refs)[0] for limit in limits] if annealing_temperature is not None: - all_data_df = pd.concat( + return pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots, + self.pseudo_equilibrium_solve( + single_chempot_dict=single_chempot_dict, quenched_temperature=quenched_temperature, annealing_temperature=annealing_temperature, effective_dopant_concentration=effective_dopant_concentration, fix_charge_states=fix_charge_states, free_defects=free_defects, ) - for chempots in tqdm(interpolated_chem_pots) + for single_chempot_dict in chempots ] ) - else: - all_data_df = pd.concat( - [ - self._solve_and_append_chempots( - chempots=chem_pots, - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - ) - for chem_pots in tqdm(interpolated_chem_pots) - ] - ) - - return all_data_df - - def _get_interpolated_chempots( - self, - chem_pot_start: dict, - chem_pot_end: dict, - n_points: int, - ) -> list: - """ - Generate a set of interpolated chemical potentials between two points. - - Args: - chem_pot_start (dict): A dictionary representing the starting chemical potentials. - The keys are element symbols and the values are their corresponding chemical potentials. - chem_pot_end (dict): A dictionary representing the ending chemical potentials. - The keys are element symbols and the values are their corresponding chemical potentials. - n_points (int): The number of interpolated points to generate, - including the start and end points. - - Returns: - list: A list of dictionaries, where each dictionary contains a set of interpolated chemical - potentials. The length of the list corresponds to `n_points`, and each dictionary corresponds - to an interpolated state between the starting and ending chemical potentials. - """ - return [ - { - key: chem_pot_start[key] + (chem_pot_end[key] - chem_pot_start[key]) * i / (n_points - 1) - for key in chem_pot_start - } - for i in range(n_points) - ] - - def _solve_and_append_chempots( - self, chempots: dict[str, float], temperature: float = 300, **kwargs - ) -> pd.DataFrame: - """ - Solve for the defect concentrations at a given temperature and set of - chemical potentials. - - Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the elements. - The keys are element symbols, and the values are their - corresponding chemical potentials. - temperature (float): - The temperature at which to solve for defect concentrations, - in Kelvin. Defaults to 300 K. - kwargs: - Additional keyword arguments that may be passed to the solver, - such as options for specific defect calculations or solver configurations. - - Returns: - pd.DataFrame: A DataFrame containing the calculated defect and carrier concentrations, - along with the self-consistent Fermi energy. The DataFrame also includes the provided - chemical potentials as additional columns. - """ - results_df = self.equilibrium_solve(chempots, temperature, **kwargs) # type: ignore - for key, value in chempots.items(): - results_df[f"μ_{key}"] = value - if kwargs.get("effective_dopant_concentration"): - results_df["Dopant (cm^-3)"] = kwargs["effective_dopant_concentration"] - return results_df - - def _solve_and_append_chempots_pseudo( - self, - chempots: dict[str, float], - annealing_temperature: float, - quenched_temperature: float = 300, - **kwargs, - ) -> pd.DataFrame: - """ - Solve for the defect concentrations using a pseudo-equilibrium - approach, given a range of chemical potentials and temperatures. - - Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the - elements. The keys are element symbols, and the values are - their corresponding chemical potentials. - annealing_temperature (float): - The temperature at which the system is annealed, in Kelvin. - quenched_temperature (float): - The temperature at which the system is quenched, in Kelvin. - Defaults to 300 K. - kwargs: - Additional keyword arguments that may be passed to the solver, - such as options for specific defect calculations or solver - configurations. - - Returns: - pd.DataFrame: A DataFrame containing the calculated defect and carrier concentrations - using the pseudo-equilibrium approach, along with the provided chemical potentials - as additional columns. - """ - results_df = self.pseudo_equilibrium_solve( - chempots, quenched_temperature, annealing_temperature, **kwargs - ) # type: ignore - for key, value in chempots.items(): - results_df[f"μ_{key}"] = value - if kwargs.get("effective_dopant_concentration"): - results_df["Dopant (cm^-3)"] = kwargs["effective_dopant_concentration"] - return results_df + return pd.concat( + [ + self.equilibrium_solve( + single_chempot_dict=single_chempot_dict, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) + for single_chempot_dict in chempots + ] + ) def scan_chemical_potential_grid( self, chempots: Optional[dict] = None, - n_points: Optional[int] = 10, + n_points: int = 10, annealing_temperature: Optional[float] = None, quenched_temperature: float = 300, temperature: float = 300, @@ -1243,7 +1522,7 @@ def scan_chemical_potential_grid( fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: - """ + r""" Given a ``doped``-formatted chemical potential dictionary, generate a ``ChemicalPotentialGrid`` object and calculate the Fermi level positions and defect/carrier concentrations at the grid points. @@ -1260,13 +1539,24 @@ def scan_chemical_potential_grid( Args: chempots (Optional[dict]): - A dictionary of chemical potentials to scan. If not provided, - the default chemical potentials from `self.chempots` will be - used, if available. - n_points (Optional[int]): + Dictionary of chemical potentials to scan over, in the ``doped`` + format (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``) + -- the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials). + + If ``None`` (default), will use ``self.defect_thermodynamics.chempots``. + + Note that you can also set + ``FermiSolver.defect_thermodynamics.chempots = ...`` or + ``DefectThermodynamics.chempots = ...`` to set the default chemical + potentials for all calculations, and you can set + ``FermiSolver.defect_thermodynamics.el_refs = ...`` or + ``DefectThermodynamics.el_refs = ...`` if you want to update the + elemental reference energies for any reason. + n_points (int): The number of points to generate along each axis of the grid. The actual number of grid points may be less, as points outside - the convex hull are excluded. Defaults to 10. + the convex hull are excluded. Default is 10. annealing_temperature (Optional[float]): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to @@ -1282,7 +1572,7 @@ def scan_chemical_potential_grid( Default is 300 K. temperature (float): The temperature at which to solve for defect concentrations - and Fermi energy, under thermodynamic equilibrium (if + and Fermi level, under thermodynamic equilibrium (if ``annealing_temperature`` is not specified). Defaults to 300 K. effective_dopant_concentration (Optional[float]): @@ -1305,53 +1595,88 @@ def scan_chemical_potential_grid( to be "frozen-in" upon quenching. Defaults to ``None``. Returns: - pd.DataFrame: A DataFrame containing the Fermi energy solutions at the grid + pd.DataFrame: A DataFrame containing the Fermi level solutions at the grid points, based on the provided chemical potentials and conditions. """ - 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 - - grid = ChemicalPotentialGrid.from_chempots(chempots).get_grid(n_points) + chempots, el_refs = self._parse_and_check_grid_like_chempots(chempots) + grid = ChemicalPotentialGrid(chempots).get_grid(n_points) if annealing_temperature is not None: - all_data_df = pd.concat( + return pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots[1].to_dict(), + self.pseudo_equilibrium_solve( + single_chempot_dict=chempot_series.to_dict(), + el_refs=el_refs, annealing_temperature=annealing_temperature, quenched_temperature=quenched_temperature, effective_dopant_concentration=effective_dopant_concentration, fix_charge_states=fix_charge_states, free_defects=free_defects, ) - for chempots in tqdm(grid.iterrows()) + for _idx, chempot_series in tqdm(grid.iterrows()) ] ) - else: - all_data_df = pd.concat( - [ - self._solve_and_append_chempots( - chempots=chempots[1].to_dict(), - temperature=temperature, - effective_dopant_concentration=effective_dopant_concentration, - ) - for chempots in tqdm(grid.iterrows()) - ] + return pd.concat( + [ + self.equilibrium_solve( + single_chempot_dict=chempot_series.to_dict(), + el_refs=el_refs, + temperature=temperature, + effective_dopant_concentration=effective_dopant_concentration, + ) + for _idx, chempot_series in tqdm(grid.iterrows()) + ] + ) + + def _parse_and_check_grid_like_chempots(self, chempots: Optional[dict] = None) -> tuple[dict, dict]: + r""" + Parse a dictionary of chemical potentials for the chemical potential + scanning functions, checking that it is in the correct format. + + Args: + chempots (Optional[dict]): + Dictionary of chemical potentials to scan over, in the ``doped`` + format (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``) + -- the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials). + + If ``None`` (default), will use ``self.defect_thermodynamics.chempots``. + + Note that you can also set + ``FermiSolver.defect_thermodynamics.chempots = ...`` or + ``DefectThermodynamics.chempots = ...`` to set the default chemical + potentials for all calculations, and you can set + ``FermiSolver.defect_thermodynamics.el_refs = ...`` or + ``DefectThermodynamics.el_refs = ...`` if you want to update the + elemental reference energies for any reason. + + Returns: + dict: + The chemical potentials in the correct format, along with the + elemental reference energies. + """ + # returns self.defect_thermodynamics.chempots if chempots is None: + chempots, el_refs = self.defect_thermodynamics._get_chempots(chempots) + + if chempots is None: + raise ValueError( + "No chemical potentials supplied or present in `self.defect_thermodynamics.chempots`!" + ) + if len(chempots["limits"]) == 1: + raise ValueError( + "Only one chemical potential limit is present in " + "`chempots`/`self.defect_thermodynamics.chempots`, which makes no sense for a chemical " + "potential grid scan (with `scan_chemical_potential_grid`/`min_max_X`/`scan_chempots`)!" ) - return all_data_df + return chempots, el_refs def min_max_X( self, target: str, min_or_max: str = "max", chempots: Optional[dict] = None, - el_refs: Optional[dict] = None, annealing_temperature: Optional[float] = None, quenched_temperature: float = 300, temperature: float = 300, @@ -1361,7 +1686,7 @@ def min_max_X( fix_charge_states: bool = False, free_defects: Optional[list[str]] = None, ) -> pd.DataFrame: - """ + r""" Search for the chemical potentials that minimize or maximize a target variable, such as electron concentration, within a specified tolerance. @@ -1383,13 +1708,25 @@ def min_max_X( Args: target (str): The target variable to minimize or maximize, e.g., "Electrons (cm^-3)". + # TODO: Allow this to match a substring of the column name. min_or_max (str): Specify whether to "minimize" ("min") or "maximize" ("max"; default) the target variable. - chempots (Optional[dict]): A dictionary of initial chemical potentials to use. - If not provided, default potentials from `self.chempots` are used. - el_refs (Optional[dict]): A dictionary of elemental reference energies used - for calculating chemical potentials relative to these references. + chempots (Optional[dict]): + Dictionary of chemical potentials to scan over, in the ``doped`` + format (i.e. ``{"limits": [{'limit': [chempot_dict]}], ...}``) + -- the format generated by ``doped``\'s chemical potential parsing + functions (see tutorials). + + If ``None`` (default), will use ``self.defect_thermodynamics.chempots``. + + Note that you can also set + ``FermiSolver.defect_thermodynamics.chempots = ...`` or + ``DefectThermodynamics.chempots = ...`` to set the default chemical + potentials for all calculations, and you can set + ``FermiSolver.defect_thermodynamics.el_refs = ...`` or + ``DefectThermodynamics.el_refs = ...`` if you want to update the + elemental reference energies for any reason. annealing_temperature (Optional[float]): Temperature in Kelvin at which to calculate the high temperature (fixed) total defect concentrations, which should correspond to @@ -1405,7 +1742,7 @@ def min_max_X( Default is 300 K. temperature (float): The temperature at which to solve for defect concentrations - and Fermi energy, under thermodynamic equilibrium (if + and Fermi level, under thermodynamic equilibrium (if ``annealing_temperature`` is not specified). Defaults to 300 K. tolerance (float): @@ -1445,9 +1782,8 @@ def min_max_X( If neither ``chempots`` nor ``self.chempots`` is provided, or if ``min_or_max`` is not ``"minimize"/"min"`` or ``"maximize"/"max"``. """ - 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) + chempots, el_refs = self._parse_and_check_grid_like_chempots(chempots) + starting_grid = ChemicalPotentialGrid(chempots) current_vertices = starting_grid.vertices chempots_labels = list(current_vertices.columns) previous_value = None @@ -1456,38 +1792,40 @@ def min_max_X( if annealing_temperature is not None: results_df = pd.concat( [ - self._solve_and_append_chempots_pseudo( - chempots=chempots[1].to_dict(), + self.pseudo_equilibrium_solve( + single_chempot_dict=chempot_series.to_dict(), + el_refs=el_refs, annealing_temperature=annealing_temperature, quenched_temperature=quenched_temperature, effective_dopant_concentration=effective_dopant_concentration, fix_charge_states=fix_charge_states, free_defects=free_defects, ) - for chempots in tqdm(starting_grid.get_grid(n_points).iterrows()) + for _idx, chempot_series in tqdm(starting_grid.get_grid(n_points).iterrows()) ] ) else: results_df = pd.concat( [ - self._solve_and_append_chempots( - chempots=chempots[1].to_dict(), + self.equilibrium_solve( + single_chempot_dict=chempot_series.to_dict(), + el_refs=el_refs, temperature=temperature, effective_dopant_concentration=effective_dopant_concentration, ) - for chempots in tqdm(starting_grid.get_grid(n_points).iterrows()) + for _idx, chempot_series in tqdm(starting_grid.get_grid(n_points).iterrows()) ] ) # Find chemical potentials value where target is lowest or highest if target in results_df.columns: if "min" in min_or_max: - target_chem_pot = results_df[results_df[target] == results_df[target].min()][ + target_chempot = results_df[results_df[target] == results_df[target].min()][ chempots_labels ] target_dataframe = results_df[results_df[target] == results_df[target].min()] elif "max" in min_or_max: - target_chem_pot = results_df[results_df[target] == results_df[target].max()][ + target_chempot = results_df[results_df[target] == results_df[target].max()][ chempots_labels ] target_dataframe = results_df[results_df[target] == results_df[target].max()] @@ -1501,20 +1839,20 @@ def min_max_X( # Find the row where "Concentration (cm^-3)" is at its minimum or maximum if "min" in min_or_max: min_value = filtered_df["Concentration (cm^-3)"].min() - target_chem_pot = results_df.loc[results_df["Concentration (cm^-3)"] == min_value][ + target_chempot = results_df.loc[results_df["Concentration (cm^-3)"] == min_value][ chempots_labels ] target_dataframe = results_df[ - results_df[chempots_labels].eq(target_chem_pot.iloc[0]).all(axis=1) + results_df[chempots_labels].eq(target_chempot.iloc[0]).all(axis=1) ] elif "max" in min_or_max: max_value = filtered_df["Concentration (cm^-3)"].max() - target_chem_pot = results_df.loc[results_df["Concentration (cm^-3)"] == max_value][ + target_chempot = results_df.loc[results_df["Concentration (cm^-3)"] == max_value][ chempots_labels ] target_dataframe = results_df[ - results_df[chempots_labels].eq(target_chem_pot.iloc[0]).all(axis=1) + results_df[chempots_labels].eq(target_chempot.iloc[0]).all(axis=1) ] current_value = ( filtered_df["Concentration (cm^-3)"].min() @@ -1529,85 +1867,61 @@ def min_max_X( ): break previous_value = current_value - target_chem_pot = target_chem_pot.drop_duplicates(ignore_index=True) + target_chempot = target_chempot.drop_duplicates(ignore_index=True) - new_vertices = [ # get midpoint between current vertices and target_chem_pot - (current_vertices + row[1]) / 2 for row in target_chem_pot.iterrows() + new_vertices = [ # get midpoint between current vertices and target_chempot + (current_vertices + row[1]) / 2 for row in target_chempot.iterrows() ] - # Generate a new grid around the target_chem_pot that + # Generate a new grid around the target_chempot that # does not go outside the bounds of the starting grid new_vertices_df = pd.DataFrame(new_vertices[0], columns=chempots_labels) starting_grid = ChemicalPotentialGrid(new_vertices_df.to_dict("index")) return target_dataframe - def _handle_chempots_for_py_sc_fermi(self, chempots: dict[str, float]) -> dict[str, float]: - """ - Adjust the provided chemical potentials for use with the py-sc-fermi - backend. - - This method ensures that the chemical potentials are correctly referenced using - the elemental reference energies stored in the class. - - Args: - chempots (dict[str, float]): A dictionary containing the chemical potentials - for the elements. The keys are element symbols, and the values are their - corresponding chemical potentials. - - Returns: - dict[str, float]: A dictionary containing the adjusted chemical potentials, - where each potential is shifted by the corresponding elemental reference energy. - - Raises: - ValueError: If `self.chempots` or `self.chempots['elemental_refs']` is None, indicating - that the necessary elemental reference energies are not available. - """ - if self.chempots is not None and self.chempots.get("elemental_refs") is not None: - elemental_refs = self.chempots["elemental_refs"] - else: - # Handle the case where self.chempots or self.chempots["elemental_refs"] is None - raise ValueError("self.chempots or self.chempots['elemental_refs'] is None") - return {k: v + elemental_refs.get(k, 0) for k, v in chempots.items()} - def _generate_dopant_for_py_sc_fermi(self, effective_dopant_concentration: float) -> "DefectSpecies": """ - Generate a dopant defect charge state object. + Generate a dopant defect charge state object, for use with the ``py-sc- + fermi`` functions. - This method creates a defect charge state object representing an arbitrary dopant or - impurity in the material, used to include in the charge neutrality condition and - analyze the Fermi level/doping response under hypothetical doping conditions. + This method creates a defect charge state object representing an + arbitrary dopant or impurity in the material, used to include in the + charge neutrality condition and analyze the Fermi level/doping + response under hypothetical doping conditions. Args: - effective_dopant_concentration (float): The fixed concentration of the dopant or impurity - in the material, specified in cm^-3. A positive value indicates donor doping (positive - defect charge state), while a negative value indicates acceptor doping (negative defect - charge state). + effective_dopant_concentration (float): + The fixed concentration of the dopant or impurity in the + material, specified in cm^-3. A positive value indicates + donor doping (positive defect charge state), while a negative + value indicates acceptor doping (negative defect charge state). Returns: - DefectSpecies: An instance of the `DefectSpecies` class, representing the generated dopant - with the specified charge state and concentration. + DefectSpecies: + An instance of the ``DefectSpecies`` class, representing the + generated dopant with the specified charge state and concentration. Raises: - ValueError: If `effective_dopant_concentration` is zero or if there is an issue with generating - the dopant. + ValueError: + If ``effective_dopant_concentration`` is zero or if there is an + issue with generating the dopant. """ self._check_required_backend_and_error("py-sc-fermi") assert self._DefectChargeState assert self._DefectSpecies - if effective_dopant_concentration > 0: - charge = 1 - effective_dopant_concentration = abs(effective_dopant_concentration) / 1e24 * self.volume - elif effective_dopant_concentration < 0: - charge = -1 - effective_dopant_concentration = abs(effective_dopant_concentration) / 1e24 * self.volume dopant = self._DefectChargeState( - charge=charge, fixed_concentration=effective_dopant_concentration, degeneracy=1 + charge=np.sign(effective_dopant_concentration), + fixed_concentration=abs(effective_dopant_concentration) / 1e24 * self.volume, + degeneracy=1, + ) + return self._DefectSpecies( + nsites=1, charge_states={np.sign(effective_dopant_concentration): dopant}, name="Dopant" ) - return self._DefectSpecies(nsites=1, charge_states={charge: dopant}, name="Dopant") def _generate_defect_system( self, - chempots: dict[str, float], + single_chempot_dict: dict[str, float], + el_refs: Optional[dict[str, float]] = None, temperature: float = 300, effective_dopant_concentration: Optional[float] = None, ) -> "DefectSystem": @@ -1621,10 +1935,23 @@ def _generate_defect_system( concentration. Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the elements. - The keys are element symbols, and the values are their corresponding - chemical potentials. + single_chempot_dict (dict[str, float]): + Dictionary of chemical potentials to use for calculating the equilibrium + Fermi level position and defect/carrier concentrations. Here, this + 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 + reference energies) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + el_refs (dict[str, float]): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}``. Unnecessary if + ``self.defect_thermodynamics.el_refs`` is set (i.e. if ``chempots`` was + provided to ``self.defect_thermodynamics`` in the format generated by + ``doped``). + (Default: None) temperature (float): The temperature at which to perform the calculations, in Kelvin. Defaults to 300 K. @@ -1650,14 +1977,14 @@ def _generate_defect_system( defect_species: dict[str, Any] = { label: {"charge_states": {}, "nsites": None, "name": label} for label in labels } - chempots = self._handle_chempots_for_py_sc_fermi(chempots) + dft_chempots = _get_dft_chempots(single_chempot_dict, el_refs) for entry in entries: label, charge = _get_label_and_charge(entry.name) defect_species[label]["nsites"] = entry.defect.multiplicity / self.multiplicity_scaling formation_energy = self.defect_thermodynamics.get_formation_energy( - entry, chempots=chempots, fermi_level=0 + entry, chempots=dft_chempots, fermi_level=0 ) total_degeneracy = np.prod(list(entry.degeneracy_factors.values())) defect_species[label]["charge_states"][charge] = { @@ -1681,8 +2008,9 @@ def _generate_defect_system( def _generate_annealed_defect_system( self, - chempots: dict[str, float], annealing_temperature: float, + single_chempot_dict: dict[str, float], + el_refs: Optional[dict[str, float]] = None, quenched_temperature: float = 300, fix_charge_states: bool = False, effective_dopant_concentration: Optional[float] = None, @@ -1701,13 +2029,26 @@ def _generate_annealed_defect_system( charge states to vary while keeping total defect concentrations fixed. Args: - chempots (dict[str, float]): - A dictionary containing the chemical potentials for the elements. - The keys are element symbols, and the values are their corresponding - chemical potentials. annealing_temperature (float): The higher temperature (in Kelvin) at which the system is annealed to set initial defect concentrations. + single_chempot_dict (dict[str, float]): + Dictionary of chemical potentials to use for calculating the equilibrium + Fermi level position and defect/carrier concentrations. Here, this + 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 + reference energies) that should be given here, otherwise the absolute + (DFT) chemical potentials should be given. + el_refs (dict[str, float]): + Dictionary of elemental reference energies for the chemical potentials + in the format: + ``{element symbol: reference energy}``. Unnecessary if + ``self.defect_thermodynamics.el_refs`` is set (i.e. if ``chempots`` was + provided to ``self.defect_thermodynamics`` in the format generated by + ``doped``). + (Default: None) quenched_temperature (float): The lower temperature (in Kelvin) to which the system is quenched. Defaults to 300 K. @@ -1734,13 +2075,14 @@ def _generate_annealed_defect_system( free_defects = [] defect_system = self._generate_defect_system( - chempots=chempots, # chempots handled in _generate_defect_system() + single_chempot_dict=single_chempot_dict, # chempots handled in _generate_defect_system() + el_refs=el_refs, temperature=annealing_temperature, effective_dopant_concentration=effective_dopant_concentration, ) initial_conc_dict = defect_system.concentration_dict() # concentrations at initial temperature - # Exclude the free_defects, carrier concentrations and Fermi energy from fixing + # Exclude the free_defects, carrier concentrations and Fermi level from fixing all_free_defects = ["Fermi Energy", "n0", "p0"] all_free_defects.extend(free_defects) @@ -1775,9 +2117,9 @@ def _generate_annealed_defect_system( class ChemicalPotentialGrid: """ - A class to represent a grid of chemical potentials and to perform + 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. + 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 @@ -1786,21 +2128,35 @@ class ChemicalPotentialGrid: properties. """ - def __init__(self, chempots: dict[str, Any]) -> None: - """ - Initializes the `ChemicalPotentialGrid` with chemical potential data. + 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[str, Any]): A dictionary containing chemical potential - information. The keys are element symbols or other identifiers, and - the values are the corresponding chemical potentials. + 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"]``)! """ - self.vertices = pd.DataFrame.from_dict(chempots, orient="index") + if "limits_wrt_el_refs" in chempots: + self.vertices = pd.DataFrame.from_dict(chempots["limits_wrt_el_refs"], orient="index") + else: + self.vertices = pd.DataFrame.from_dict(chempots, orient="index") - def get_grid(self, n_points: Optional[int] = None) -> pd.DataFrame: + 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. @@ -1811,65 +2167,51 @@ def get_grid(self, n_points: Optional[int] = None) -> pd.DataFrame: chemical potential values at these points. Args: - n_points (int): The number of points to generate along each axis of the grid. + 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. + 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. + 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. """ - if n_points is None: - n_points = 100 return self.grid_from_dataframe(self.vertices, n_points) - @classmethod - def from_chempots(cls, chempots: dict[str, Any]) -> "ChemicalPotentialGrid": - """ - Initializes the ChemicalPotentialGrid with chemical potential data. - - This class method creates an instance of the `ChemicalPotentialGrid` class - using a dictionary of chemical potentials. It extracts the relevant data - from the provided dictionary and initializes the grid with this data. - - Args: - chempots (dict[str, Any]): A dictionary containing chemical potential - information. The keys are identifiers for the chemical potential limits, - and the values provide the corresponding chemical potential values relative - to elemental references. - - Returns: - ChemicalPotentialGrid: An instance of the `ChemicalPotentialGrid` class - initialized with the provided chemical potential data. - """ - return cls(chempots["limits_wrt_el_refs"]) - @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's independent variables. It interpolates the values of the - dependent variable over this grid, ensuring that all generated points lie - within the convex hull of the given vertices. + This method creates a grid of points within the convex hull + defined by the input ``DataFrame``\'s independent variables. + It interpolates the values of the dependent variable 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. + 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, with the last column containing the interpolated - dependent variable values. + 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, with the last column containing the interpolated + dependent variable values. """ # Exclude the dependent variable from the vertices dependent_variable = mu_dataframe.columns[-1] diff --git a/doped/thermodynamics.py b/doped/thermodynamics.py index 10991b16..a89d0349 100644 --- a/doped/thermodynamics.py +++ b/doped/thermodynamics.py @@ -458,9 +458,9 @@ def __init__( (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. - If None (default), sets all chemical potentials to zero. Chemical - potentials can also be supplied later in each analysis function. - (Default: None) + If ``None`` (default), sets all chemical potentials to zero. Chemical + potentials can also be supplied later in each analysis function, or set + using ``DefectThermodynamics.chempots = ...`` (with the same input options). el_refs (dict): Dictionary of elemental reference energies for the chemical potentials in the format: @@ -468,7 +468,10 @@ def __init__( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided in format generated by ``doped`` (see tutorials). - (Default: None) + + If ``None`` (default), sets all elemental reference energies to zero. Reference + energies can also be supplied later in each analysis function, or set + using ``DefectThermodynamics.el_refs = ...`` (with the same input options). vbm (float): VBM eigenvalue to use as Fermi level reference point for analysis. If None (default), will use ``"vbm"`` from the ``calculation_metadata`` @@ -1372,12 +1375,18 @@ def get_equilibrium_concentrations( limits) can then be chosen using ``limit``. Alternatively this can be a dictionary of chemical potentials for a - single limit (``limit``), in the format: ``{element symbol: chemical potential}``. + single limit (``limit``), in the format: + ``{element symbol: chemical potential}``. If manually specifying chemical potentials this way, you can set the ``el_refs`` option with the DFT reference energies of the elemental phases, in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. + (Default: None) limit (str): The chemical potential limit for which to obtain the equilibrium concentrations. Can be either: @@ -1399,6 +1408,10 @@ def get_equilibrium_concentrations( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) fermi_level (float): Value corresponding to the electron chemical potential, referenced @@ -1621,6 +1634,11 @@ def get_equilibrium_fermi_level( in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. + (Default: None) limit (str): The chemical potential limit for which to determine the equilibrium Fermi level. Can be either: @@ -1642,6 +1660,10 @@ def get_equilibrium_fermi_level( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) temperature (float): Temperature in Kelvin at which to calculate the equilibrium Fermi level. @@ -1814,6 +1836,11 @@ def get_quenched_fermi_level_and_concentrations( in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. + (Default: None) limit (str): The chemical potential limit for which to determine the Fermi level and concentrations. Can be either: @@ -1835,6 +1862,10 @@ def get_quenched_fermi_level_and_concentrations( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) annealing_temperature (float): Temperature in Kelvin at which to calculate the high temperature @@ -2109,6 +2140,9 @@ def get_formation_energy( (DFT) chemical potentials should be given. If None (default), sets all chemical potentials to zero. + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. (Default: None) limit (str): The chemical potential limit for which to @@ -2131,6 +2165,10 @@ def get_formation_energy( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) fermi_level (float): Value corresponding to the electron chemical potential, referenced @@ -2243,6 +2281,11 @@ def get_dopability_limits( in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. + (Default: None) limit (str): The chemical potential limit for which to calculate formation energies (and thus dopability limits). Can be either: @@ -2264,6 +2307,10 @@ def get_dopability_limits( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) Returns: @@ -2421,6 +2468,11 @@ def get_doping_windows( in which case it is the formal chemical potentials (i.e. relative to the elemental references) that should be given here, otherwise the absolute (DFT) chemical potentials should be given. + + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. + (Default: None) limit (str): The chemical potential limit for which to calculate formation energies (and thus doping windows). Can be either: @@ -2442,6 +2494,10 @@ def get_doping_windows( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) Returns: @@ -2595,6 +2651,9 @@ def plot( otherwise the absolute (DFT) chemical potentials should be given. If None (default), sets all chemical potentials to zero. + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. (Default: None) limit (str): The chemical potential limit for which to @@ -2615,6 +2674,10 @@ def plot( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) chempot_table (bool): Whether to print the chemical potential table above the plot. @@ -2943,6 +3006,9 @@ def get_formation_energies( (DFT) chemical potentials should be given. If None (default), sets all chemical potentials to zero. + Note that you can also set ``DefectThermodynamics.chempots = ...`` + (with the same input options) to set the default chemical potentials + for all calculations. (Default: None) limit (str): The chemical potential limit for which to @@ -2963,6 +3029,10 @@ def get_formation_energies( potentials, when ``chempots`` has been manually specified as ``{element symbol: chemical potential}``). Unnecessary if ``chempots`` is provided/present in format generated by ``doped`` (see tutorials). + + Note that you can also set ``DefectThermodynamics.el_refs = ...`` + (with the same input options) to set the default elemental reference + energies for all calculations. (Default: None) fermi_level (float): Value corresponding to the electron chemical potential, referenced