Tutorial¶
A complete worked example: building a DefectSystem from defect formation energies and a density of states, solving for the self-consistent Fermi level and the equilibrium defect, electron, and hole concentrations, and then working through temperature dependence, concentration and site constraints, metastable charge states, and temperature-dependent band-edge corrections.
1. Building a defect system¶
The DefectChargeState object¶
The smallest “unit” of a DefectSystem object in py-sc-fermi is a DefectChargeState object: this is where the formation energies of the defect (which you have likely calculated from a density functional theory study) are defined. Let’s take the example of NaCl, a rock salt structured material, and assume the only defects that can form are vacancies, \(V_\mathrm{Na}\) and \(V_\mathrm{Cl}\). We’ll assume each vacancy has two possible charge states, \(V^{0}_\mathrm{Na}\) and
\(V^{-1}_\mathrm{Na}\) and \(V^{0}_\mathrm{Cl}\) and \(V^{+1}_\mathrm{Cl}\). The cell below shows how we might define these defect objects.
[1]:
from py_sc_fermi.defect_charge_state import DefectChargeState
v_Na_0 = DefectChargeState(charge= 0, # the charge of the defect
degeneracy= 1, # see the usage notes page in the documentation for a discussion of degeneracy values
energy = 1.0 # defect formation energy in eV, we will just use a dummy value here
)
v_Na_minus_1 = DefectChargeState(charge= -1, degeneracy= 1, energy = 2.0)
# repeat for the Cl vacancy
v_Cl_0 = DefectChargeState(charge= 0, degeneracy= 1, energy = 1.0)
v_Cl_1 = DefectChargeState(charge= 1, degeneracy= 1, energy = -1.0)
The DefectSpecies object¶
The DefectSpecies object collects DefectChargeStates together under a single “kind” of defect, this is useful for assigning consistent site degeneracies (and specifying different constraints, which we will come to later). These are defined like so:
[2]:
from py_sc_fermi.defect_species import DefectSpecies
# defining the Na vacancy
v_Na = DefectSpecies(name = "v_Na", # the name can be any string, it only needs to be unique in the set of all species
nsites = 1, # the number of sites on which the defect can form in the unit cell (see usage notes)
charge_states = [v_Na_0, v_Na_minus_1]) # list of DefectChargeState objects
# defining the Cl vacancy
v_Cl = DefectSpecies(name = "v_Cl", nsites = 1, charge_states = [v_Cl_0, v_Cl_1])
The DOS object¶
The DOS object contains the density of states information calculated for the unit cell (see usage notes in the documentation for a strict definition of the unit cell) we take data calculated from the Materials Project for NaCl (mp-22862) which we have dumped into a csv file.
[3]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format='retina'
from py_sc_fermi.dos import DOS
energy = np.loadtxt('dos.csv', delimiter=',')[0] # the energy range of the DOS information
dos = np.loadtxt('dos.csv', delimiter=',')[1] # the (total) DOS
dos = DOS(dos = dos, # the total density of states for the unit cell
edos = energy, # the energy range of the calculated density of states
bandgap = 5, # the bandgap in eV (taken from the Materials Project, mp-22862)
nelect = 20, # the number of electrons in the unit cell calculation
spin_polarised=False) # This dos calculation is not spin polarised
Putting it all together: the DefectSystem object¶
The DefectSystem object collects our DefectSpecies together, and allows for us to search for the self-consistent Fermi energy. It is defined as below:
[4]:
from py_sc_fermi.defect_system import DefectSystem
defect_system = DefectSystem(defect_species = [v_Na, v_Cl], # list of py_sc_fermi.defect_species.DefectSpecies objects
dos = dos, # py_sc_fermi.dos.DOS object
volume = 46.10, # the volume of the unit cell in angstrom^3 (taken from the Materials Project, mp-22862)
temperature = 500) # the temperature in K
Printing defect_system.result gives a summary of the solved system: the self-consistent Fermi energy, the electron and hole concentrations, and the equilibrium defect concentrations broken down by charge state.
[5]:
print(defect_system.result)
SC Fermi energy: 1.500000 eV
Carriers:
p0 (holes): 1.138e+06 cm⁻³
n0 (electrons): 3.921e-18 cm⁻³
Defect concentrations:
v_Na 1.979e+17 cm⁻³
q+0 1.806e+12 cm⁻³ ( 0.00%)
q-1 1.979e+17 cm⁻³ (100.00%)
v_Cl 1.979e+17 cm⁻³
q+0 1.806e+12 cm⁻³ ( 0.00%)
q+1 1.979e+17 cm⁻³ (100.00%)
Snapshots at other temperatures: DefectSystemFactory¶
DefectSystem is an immutable, fixed-temperature snapshot: defect_system.temperature is fixed at construction, and the defect_species it stores are independent copies of whatever was passed in — construction-time options such as fixed_concentrations or formation_energy_corrections act on those copies and never affect v_Na, v_Cl, or any other DefectSystem built from them.
To get a DefectSystem for the same defect species at a different temperature, build a DefectSystemFactory once from the same defect_species, dos and volume, then call .at(T) — this returns a fresh, independent DefectSystem for that temperature. We’ll use this throughout the rest of the tutorial whenever we need a snapshot at a different temperature, instead of mutating defect_system.temperature or copying it by hand.
[6]:
from py_sc_fermi.defect_system import DefectSystemFactory
# v_Na and v_Cl exceed the 1% site-occupancy threshold at the high end of the
# temperature sweeps below (about 14% at 3000 K), so the DiluteLimitWarning is
# disabled for this toy system. In a real calculation it should be heeded: high
# occupancy means the dilute, non-interacting-defect assumption is breaking down.
defect_system_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume,
occupancy_warning_threshold=None,
)
2. Misc. DefectSystem analysis¶
There are various analyses we can perform once our defect system is defined, some non-exhaustive examples are included in the cells below:
Plot transition level diagrams¶
[7]:
# plot transition levels
# get a dictionary of the lines for a transition level (TL) diagram where the key is the name of the DefectSpecies
# the values is a list of lists, where the first list is the x coordinates of the TL profile and the second list
# is the y coordinates of the TL profile
transition_levels = defect_system.get_transition_levels()
# get the self consistent Fermi energy to plot on the TL diagram
sc_fermi = defect_system.get_sc_fermi()[0]
# plot each defect species TL profile on the full diagram, and label as
# the defect species name
for k, v in transition_levels.items():
plt.plot(v[0], v[1], "-o", label=k)
# plot the self consistent Fermi energy
plt.vlines(sc_fermi, 0, 2, linestyles="dashed", label = "SC Fermi level")
plt.ylim(0, 2)
plt.xlim(0, 5)
plt.legend()
plt.ylabel("Formation energy / eV")
plt.xlabel("Fermi Energy / eV")
plt.show()
Plot defect concentrations as a function of temperature¶
[8]:
# plot the defect concentration as a function of temperature
temperatures = np.linspace(400, 3000, 100)
v_Cl_conc = []
v_Na_conc = []
for t in temperatures:
concentrations = defect_system_factory.at(t).result.concentrations
v_Cl_conc.append(concentrations["v_Cl"])
v_Na_conc.append(concentrations["v_Na"])
plt.plot(temperatures, v_Cl_conc, label=r"$[V_\mathrm{Cl}]$", ls = "--")
plt.plot(temperatures, v_Na_conc, label=r"$[V_\mathrm{Na}]$", ls = "--")
plt.yscale("log")
plt.legend()
plt.ylabel("Concentration / cm^-3")
plt.xlabel("temperature / K")
plt.show()
Plot defect concentrations as a function of Fermi energy¶
DefectSystem.defect_concentrations_at_fermi_level(e_fermi) gives the site-excluded (Langmuir) concentration of each defect at a specified Fermi energy, in cm^-3 (the same statistics as the solve, capped at each species’ available sites). Sweeping the Fermi energy traces the compensation behaviour behind a material’s doping limits: as an external dopant pushes the Fermi level toward a band edge, the native defects that rise to compensate are read straight off these curves. The equilibrium
concentrations, at the self-consistent Fermi level, are on defect_system.result.
Where a compensating defect approaches one per available site its curve flattens at that site-density ceiling. Read the curves for where compensation drives a defect towards saturation, not for the plateau height itself: at such high occupancy the dilute, non-interacting assumption underlying the model breaks down — the regime that DiluteLimitWarning flags for a solved system.
[9]:
# plot defect concentrations as a function of Fermi energy
fermi_energies = np.linspace(0, 3, 100)
# plot the concentrations of the defects as a function of Fermi energy
v_Cl_conc = []
v_Na_conc = []
for e_fermi in fermi_energies:
conc = defect_system.defect_concentrations_at_fermi_level(e_fermi) # cm^-3, site-excluded
v_Cl_conc.append(conc["v_Cl"])
v_Na_conc.append(conc["v_Na"])
plt.plot(fermi_energies, v_Cl_conc, label=r"$[V_\mathrm{Cl}]$")
plt.plot(fermi_energies, v_Na_conc, label=r"$[V_\mathrm{Na}]$")
plt.yscale("log")
# plot the concentration of electrons as a function of Fermi energy
h_conc = []
e_conc = []
for e_fermi in fermi_energies:
h_conc.append(dos.carrier_concentrations(e_fermi, 500)[0] * 1e24 / defect_system.volume)
e_conc.append(dos.carrier_concentrations(e_fermi, 500)[1] * 1e24 / defect_system.volume)
plt.plot(fermi_energies, h_conc, label="$[h^+]$")
plt.plot(fermi_energies, e_conc, label="$[e^-]$")
plt.legend()
plt.ylabel("Concentration / cm^-3")
plt.xlabel("Fermi Energy / eV")
plt.show()
3. Applying concentration constraints¶
Fixing the concentrations of DefectSpecies¶
There are a range of constraints we can apply to the concentrations of the defects in a DefectSystem, the first example is to fix the concentration of the DefectSpecies. This is purely a thought experiment, but may represent a system where barriers for ionic motion are very large, but electronic transport is facile. Here we will look at what happens to the hole concentration if we run our NaCl system at high temperatures, then rerun at low temperatures.
[10]:
# calculate defect and electronic carrier concentrations at high T.
live_defect_system = defect_system_factory.at(1500)
high_t = live_defect_system.result
# repeat at low T
low_t_defect_system = defect_system_factory.at(300)
low_t = low_t_defect_system.result
# quench: fix the species totals solved at high T, and re-solve the
# electron holes at low T
high_t_concs = high_t.concentrations_per_cell
fixed_species_defect_system = defect_system_factory.at(
300,
fixed_concentrations={
"v_Na": high_t_concs["v_Na"],
"v_Cl": high_t_concs["v_Cl"],
},
)
fixed_species = fixed_species_defect_system.result
# plot the results
plt.bar([0, 1, 2], [high_t.p0, low_t.p0, fixed_species.p0])
plt.xticks([0, 1, 2], ["high T", "low T", "fixed species"])
plt.ylabel("Concentration / cm^-3")
plt.yscale("log")
We find the hole concentration in the fixed system (high temperature defect concentrations with electrons and holes recalculated at lower temperatures) is inflated relative to the low temperature system.
Fixing the concentrations of DefectChargeStates¶
What about if we have a system where we imagine neither ionic defects or their charge states can re-equilibrate on the timescale of interest, here we take high temperature charged defects and fix their concentrations, rerunning at lower temperatures, and find this time, that the hole concentration in the system with fixed concentration charge states ~ to the system calculated at high temperatures.
[11]:
# calculate defect and electronic carrier concentrations at high T.
live_defect_system = defect_system_factory.at(1500)
high_t = live_defect_system.result
# repeat at low T
low_t_defect_system = defect_system_factory.at(300)
low_t = low_t_defect_system.result
# quench at the charge-state level: fix the per-state concentrations
# solved at high T, keyed by (species_name, charge_state_name)
high_t_cs_concs = high_t.charge_state_concentrations_per_cell
fixed_species_defect_system = defect_system_factory.at(
300,
fixed_concentrations={
("v_Na", "q-1"): high_t_cs_concs["v_Na"]["q-1"],
("v_Cl", "q+1"): high_t_cs_concs["v_Cl"]["q+1"],
},
)
fixed_species = fixed_species_defect_system.result
plt.bar([0, 1, 2], [high_t.p0, low_t.p0, fixed_species.p0])
plt.xticks([0, 1, 2], ["high T", "low T", "fixed charge states"])
plt.ylabel("Concentration / cm^-3")
plt.yscale("log")
Introducing an additional fixed concentration defect¶
Finally, we can introduce an additional defect to model how the defect concentrations will respond to e.g. the introduction of a charged dopant.
[12]:
dopant_concentrations = np.logspace(17,20)
v_Na_conc = []
v_Cl_conc = []
p_0 = []
n_0 = []
for dopant_concentration in dopant_concentrations:
dopant_charge_state = DefectChargeState(1, fixed_concentration = dopant_concentration / 1e24 * defect_system.volume)
dopant_species = DefectSpecies("dopant", 1, charge_states = [dopant_charge_state])
live_defect_system = DefectSystem(defect_species = [dopant_species,
v_Na, v_Cl], volume = defect_system.volume, temperature = 500, dos = defect_system.dos)
result = live_defect_system.result
v_Na_conc.append(result.concentrations["v_Na"])
v_Cl_conc.append(result.concentrations["v_Cl"])
p_0.append(result.p0)
n_0.append(result.n0)
plt.plot(dopant_concentrations, v_Na_conc, label=r"$V_\mathrm{Na}$")
plt.plot(dopant_concentrations, v_Cl_conc, label=r"$V_\mathrm{Cl}$")
plt.plot(dopant_concentrations, p_0, label="$p_0$")
plt.plot(dopant_concentrations, n_0, label="$n_0$")
plt.ylabel("Concentration / cm^-3")
plt.xlabel("Dopant concentration / cm^-3")
plt.yscale("log")
plt.xscale("log")
plt.legend()
plt.show()
Site competition with site_pools¶
By default, every DefectSpecies is assigned its own site-exclusion group of nsites sites per unit cell, and its variable charge states follow Langmuir-type statistics within that group: \(c_i = n_\mathrm{free} \cdot w_i / (1 + \sum_j w_j)\).
This breaks down when multiple defect species occupy the same physical sites — for example, a vacancy and a substitutional impurity cannot simultaneously occupy the same lattice site, so they should compete for one shared budget rather than each getting their own.
The site_pools parameter groups DefectSpecies that compete for the same pool of sites. Each pool is a SitePool(n_sites=N_pool, species=[...]), where n_sites is the number of sites per unit cell shared by the listed species.
Here we add a Mg dopant substituting on Na sites (Mg_Na), which competes with v_Na for the available Na sites. Because Mg is 2+ and occupies a 1+ site, Mg_Na carries a formal charge of +1.
To see this competition directly, we fix Mg_Na at half of the Na-site pool and compare the resulting v_Na concentration with and without the shared site_pools constraint.
Fixing Mg_Na at half the Na-site pool puts it at 50% site occupancy, well outside py-sc-fermi’s dilute, non-interacting-defect assumption, so solving these systems emits a DiluteLimitWarning. That is the expected signal here — the fixed 50% occupancy is exactly what drives the site competition being demonstrated — so this comparison is illustrative rather than a quantitative prediction.
[13]:
from py_sc_fermi.pools import SitePool
# Mg substituting on the Na site (charge +1, since Mg is 2+ on a 1+ site)
Mg_Na_1 = DefectChargeState(charge=1, degeneracy=1, energy=-0.5)
Mg_Na_0 = DefectChargeState(charge=0, degeneracy=1, energy=0.3)
Mg_Na = DefectSpecies(name="Mg_Na", nsites=1, charge_states=[Mg_Na_0, Mg_Na_1])
# A copy of Mg_Na with its concentration fixed at half the Na-site pool
Mg_Na_fixed = DefectSpecies(
name="Mg_Na", nsites=1, charge_states=[Mg_Na_0, Mg_Na_1], fixed_concentration=0.5
)
# v_Na and Mg_Na_fixed each get their own 1-site budget
separate_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl, Mg_Na_fixed],
dos=dos,
volume=46.10,
)
# v_Na and Mg_Na_fixed share the 1 Na site per unit cell
pool_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl, Mg_Na_fixed],
dos=dos,
volume=46.10,
site_pools={"Na_sites": SitePool(n_sites=1.0, species=[v_Na, Mg_Na_fixed])},
)
# Mg_Na_fixed sits at 50% occupancy, so each solve below emits a DiluteLimitWarning
# (see the markdown above) — that high occupancy is the point of this demonstration.
separate_concs = separate_factory.at(1500).result.concentrations
pool_concs = pool_factory.at(1500).result.concentrations
print(f"[v_Na], own site budget: {separate_concs['v_Na']:.3e} cm⁻³")
print(f"[v_Na], shared with Mg_Na: {pool_concs['v_Na']:.3e} cm⁻³")
[v_Na], own site budget: 4.765e+20 cm⁻³
[v_Na], shared with Mg_Na: 3.343e+20 cm⁻³
/var/folders/np/pvnqvv3958gfjdtbgbbpjqk00000gn/T/ipykernel_40677/4273524033.py:30: DiluteLimitWarning: 'Mg_Na' (50%), 'v_Na' (2.2%) and 'v_Cl' (1.99%) reach high site occupancy at the self-consistent Fermi level. py-sc-fermi assumes dilute, non-interacting defects, so results at this occupancy may be non-physical.
separate_concs = separate_factory.at(1500).result.concentrations
/var/folders/np/pvnqvv3958gfjdtbgbbpjqk00000gn/T/ipykernel_40677/4273524033.py:31: DiluteLimitWarning: 'Mg_Na' (50%), 'v_Na' (1.54%) and 'v_Cl' (1.42%) reach high site occupancy at the self-consistent Fermi level. py-sc-fermi assumes dilute, non-interacting defects, so results at this occupancy may be non-physical.
pool_concs = pool_factory.at(1500).result.concentrations
[14]:
import warnings
from py_sc_fermi.warnings import DiluteLimitWarning
# These factories hold Mg_Na fixed at 50% occupancy, so every .at(t) repeats the
# DiluteLimitWarning shown in the cell above. Silence those per-temperature repeats
# here so the 60-point sweep stays readable.
temperatures = np.linspace(300, 1500, 60)
separate_v_Na = []
pool_v_Na = []
with warnings.catch_warnings():
warnings.simplefilter("ignore", DiluteLimitWarning)
for t in temperatures:
separate_v_Na.append(separate_factory.at(t).result.concentrations["v_Na"])
pool_v_Na.append(pool_factory.at(t).result.concentrations["v_Na"])
plt.plot(temperatures, separate_v_Na, label=r"$v_\mathrm{Na}$, own site budget")
plt.plot(temperatures, pool_v_Na, label=r"$v_\mathrm{Na}$, shared with $\mathrm{Mg}_\mathrm{Na}$", ls="--")
plt.yscale("log")
plt.xlabel("Temperature / K")
plt.ylabel(r"$[v_\mathrm{Na}]$ / cm$^{-3}$")
plt.legend()
plt.title(r"Site competition between $v_\mathrm{Na}$ and $\mathrm{Mg}_\mathrm{Na}$")
plt.show()
Fixing elemental content with element_pools¶
site_pools controls where a defect species can sit. element_pools controls how much of a given element is present in total — distributing that element across multiple defect species according to their statistical weights.
For example, a Mg dopant might sit substitutionally on Na sites (Mg_Na) or form an interstitial (Mg_i). If the total Mg content is fixed (e.g. by the synthesis conditions), we can model this with element_pools:
element_pools = {
# target = total Mg content; members map each species to its Mg atoms per defect
"Mg": ElementPool(target=N_Mg, members={Mg_Na: 1.0, Mg_i: 1.0})
}
The defect statistics determine the relative distribution of Mg between sites; the pool constraint fixes the absolute total, by solving for the element’s chemical potential. The constraint composes with site_pools: the chemical potential enters each state’s weight inside its site-exclusion group, so site competition and the element constraint are satisfied together in a single solve.
[15]:
from py_sc_fermi.pools import ElementPool
# Mg interstitial — competes with Mg_Na for the total Mg budget.
# At low Fermi energies (p-type), Mg_i_1 (E=-0.7 eV) is the stable state and
# dominates the Mg budget. As temperature rises and E_F shifts upward, the +1
# charge penalty grows and Mg_Na_0 (charge-neutral, E=0.3 eV) takes over —
# producing a clear crossover in the fraction plot below.
Mg_i_2 = DefectChargeState(charge=2, degeneracy=1, energy=0.5)
Mg_i_1 = DefectChargeState(charge=1, degeneracy=1, energy=-0.7)
Mg_i = DefectSpecies(name="Mg_i", nsites=1, charge_states=[Mg_i_1, Mg_i_2])
# Fix total Mg content at 1e-3 per unit cell (~2e19 cm-3)
N_Mg = 1e-3
# Both factories below cross the 1% site-occupancy threshold at the temperatures used
# here (unconstrained Mg_Na grows with T; the vacancies exceed it near 1500 K), so the
# DiluteLimitWarning is disabled for these toy systems.
# No element constraint: Mg_Na and Mg_i take their unconstrained equilibrium concentrations
free_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl, Mg_Na, Mg_i],
dos=dos,
volume=46.10,
occupancy_warning_threshold=None,
)
# Element pool: total Mg content fixed at N_Mg per unit cell
elem_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl, Mg_Na, Mg_i],
dos=dos,
volume=46.10,
element_pools={"Mg": ElementPool(target=N_Mg, members={Mg_Na: 1.0, Mg_i: 1.0})},
occupancy_warning_threshold=None,
)
free_system = free_factory.at(1000)
elem_system = elem_factory.at(1000)
print(elem_system)
print()
scale = 1e24 / 46.10
free_cd = free_system.result.concentrations
elem_cd = elem_system.result.concentrations
print(f"{'':20s} {'free':>12s} {'element pool':>12s}")
print(f"{'Mg_Na':20s} {free_cd['Mg_Na']:12.3e} {elem_cd['Mg_Na']:12.3e} cm⁻³")
print(f"{'Mg_i':20s} {free_cd['Mg_i']:12.3e} {elem_cd['Mg_i']:12.3e} cm⁻³")
print(f"{'total Mg':20s} {free_cd['Mg_Na']+free_cd['Mg_i']:12.3e} {elem_cd['Mg_Na']+elem_cd['Mg_i']:12.3e} cm⁻³")
print(f"{'target':20s} {'—':>12s} {N_Mg * scale:12.3e} cm⁻³")
DefectSystem
bandgap: 5 eV nelect: 20
volume: 46.1 ų temperature: 1000 K
4 defect species:
v_Na [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = -1 E = 2.0000 eV (deg. 1)
v_Cl [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = +1 E = -1.0000 eV (deg. 1)
Mg_Na [nsites: 1]
q = +0 E = 0.3000 eV (deg. 1)
q = +1 E = -0.5000 eV (deg. 1)
Mg_i [nsites: 1]
q = +1 E = -0.7000 eV (deg. 1)
q = +2 E = 0.5000 eV (deg. 1)
element pools:
Mg: 0.001 per cell → [Mg_Na ×1, Mg_i ×1]
free element pool
Mg_Na 6.476e+20 2.163e+19 cm⁻³
Mg_i 1.982e+18 6.534e+16 cm⁻³
total Mg 6.496e+20 2.169e+19 cm⁻³
target — 2.169e+19 cm⁻³
[16]:
# Show how Mg is distributed between substitutional and interstitial sites as a function of temperature
temperatures = np.linspace(300, 1500, 60)
Mg_Na_conc = []
Mg_i_conc = []
Mg_total = []
for t in temperatures:
cd = elem_factory.at(t).result.concentrations
Mg_Na_conc.append(cd["Mg_Na"])
Mg_i_conc.append(cd["Mg_i"])
Mg_total.append(cd["Mg_Na"] + cd["Mg_i"])
plt.plot(temperatures, Mg_Na_conc, label=r"$\mathrm{Mg}_\mathrm{Na}$")
plt.plot(temperatures, Mg_i_conc, label=r"$\mathrm{Mg}_i$", ls="--")
plt.plot(temperatures, Mg_total, label="total Mg", ls=":", color="k")
plt.xlabel("Temperature / K")
plt.ylabel("Concentration / cm⁻³")
plt.yscale("log")
plt.legend()
plt.title("Distribution of fixed Mg content between sites")
plt.show()
4. Metastable defect charge states¶
A DefectSpecies can hold more than one DefectChargeState with the same formal charge, to represent metastable configurations of the same defect — for example, an interstitial dopant that can sit on two different sites with the same charge but different formation energies (and possibly different degeneracies).
DefectSpecies.get_formation_energies, tl_profile and get_transition_level_and_energy group DefectChargeStates by their formal charge. By default (temperature=0) each charge is represented by its lowest-energy form. Passing temperature>0 instead combines same-charge forms via a Boltzmann-weighted effective formation energy, \(-k_BT\ln\sum_i g_i e^{-E_i/k_BT}\).
DefectSpecies.effective_formation_energy applies the same Boltzmann weighting across all of a species’ charge states and metastable forms at once, giving the formation energy at any Fermi energy — not just the kinks returned by tl_profile.
As an example, here’s a hypothetical interstitial dopant Al_i, with two distinct +1 configurations (“tetrahedral” and “octahedral”) in addition to its 0 and +3 states. The two +1 configurations share a formal charge, so they are given explicit names: DefectSpecies requires charge-state names to be unique, and unnamed states default to their charge string ("q+1"), which would collide.
[17]:
Al_i_0 = DefectChargeState(charge=0, degeneracy=1, energy=2.0)
Al_i_1_tet = DefectChargeState(charge=1, degeneracy=1, energy=0.8, name="Al_i_1_tet") # tetrahedral interstitial site
Al_i_1_oct = DefectChargeState(charge=1, degeneracy=2, energy=1.0, name="Al_i_1_oct") # octahedral interstitial site (metastable)
Al_i_3 = DefectChargeState(charge=3, degeneracy=1, energy=-1.0)
Al_i = DefectSpecies(name="Al_i", nsites=1,
charge_states=[Al_i_0, Al_i_1_tet, Al_i_1_oct, Al_i_3])
print("formation energies at E_F = 0, T = 0 K: ", Al_i.get_formation_energies(0))
print("formation energies at E_F = 0, T = 1500 K:", Al_i.get_formation_energies(0, temperature=1500))
formation energies at E_F = 0, T = 0 K: {0: 2.0, 1: 0.8, 3: -1.0}
formation energies at E_F = 0, T = 1500 K: {0: 2.0, 1: 0.7541601883699446, 3: -1.0}
At T = 0 K, the +1 charge takes the lower of its two formation energies (0.8 eV, the tetrahedral form). At T = 1500 K, the more degenerate octahedral form is also thermally populated, pulling the effective +1 formation energy slightly below 0.8 eV.
The same effect is visible across the whole Fermi-energy range: effective_formation_energy(e_fermi, temperature=0) reproduces the tl_profile lower envelope at every e_fermi (not just its kinks), while effective_formation_energy(e_fermi, temperature=1500) is the smoothed, finite-temperature curve.
[18]:
e_fermi_grid = np.linspace(0, 5, 200)
eff_0 = [Al_i.effective_formation_energy(e, temperature=0) for e in e_fermi_grid]
eff_1500 = [Al_i.effective_formation_energy(e, temperature=1500) for e in e_fermi_grid]
tl = Al_i.tl_profile(0, 5)
plt.plot(e_fermi_grid, eff_0, label="effective formation energy, T = 0 K")
plt.plot(e_fermi_grid, eff_1500, label="effective formation energy, T = 1500 K", ls="--")
plt.plot(tl[:, 0], tl[:, 1], "o", color="k", label="tl_profile kinks, T = 0 K")
plt.ylabel("Formation energy / eV")
plt.xlabel("Fermi energy / eV")
plt.legend()
plt.show()
To inspect per-charge-state concentrations at the DefectSystem level, use result.charge_state_concentrations. This is a {species_name: {name: conc}} dict with one entry per DefectChargeState, in cm\(^{-3}\). The inner key is the charge state’s name: an explicit name where one was set, otherwise the charge-derived default such as "q+3". result.charge_state_concentrations_per_cell gives the same breakdown per unit cell.
The self-consistent Fermi energy and the electron and hole concentrations are on the same result object too, as result.fermi_energy, result.p0 and result.n0.
[19]:
Al_i_system = DefectSystem(
defect_species=[Al_i],
dos=dos,
volume=46.10,
temperature=1500,
)
result = Al_i_system.result
cs_concs = result.charge_state_concentrations_per_cell
print("Al_i charge-state concentrations at 1500 K (per unit cell):")
for key, conc in cs_concs["Al_i"].items():
print(f" {key}: {conc:.3e}")
# the same breakdown in cm⁻³, with the Fermi energy alongside on the result object
print(f"\nFermi energy: {result.fermi_energy:.3f} eV")
print("Al_i charge-state concentrations (cm⁻³):")
for key, conc in result.charge_state_concentrations["Al_i"].items():
print(f" {key}: {conc:.3e}")
Al_i charge-state concentrations at 1500 K (per unit cell):
q+0: 1.907e-07
Al_i_1_tet: 5.091e-13
Al_i_1_oct: 2.167e-13
q+3: 3.498e-26
Fermi energy: 2.859 eV
Al_i charge-state concentrations (cm⁻³):
q+0: 4.136e+15
Al_i_1_tet: 1.104e+10
Al_i_1_oct: 4.700e+09
q+3: 7.588e-04
5. Temperature-dependent band-edge corrections¶
Quantities like the band gap and defect formation energies can themselves depend on temperature, e.g. through thermal expansion or electron-phonon coupling. DefectSystem accepts pre-evaluated corrections for a single temperature:
vbm_shift/cbm_shift: shifts (in eV) of the valence-band maximum and conduction-band minimum. They scissor the DOS so the band gap changes bycbm_shift - vbm_shift(the VBM is pinned at E=0), which moves the carrier concentrations and the self-consistent Fermi level, and sets the displayed band gap. A shifted system stores its own scissored copy; theDOSyou pass is never modified.rigid_shift(defaultTrue): gates the formation-energy channel only; the DOS scissor always applies. If the band structure and defect levels move together as a rigid body, formation energies are unchanged. IfFalse, the defect levels are fixed in absolute energy while the bands move, so each variable-concentration charge state’s formation energy shifts bycharge * vbm_shift.formation_energy_corrections: adictof one-off corrections (e.g. from a hybrid-functional or GW calculation) for specific charge states, keyed per charge state — as(species_name, charge_state_name)pairs or theDefectChargeStateobjects themselves — so that metastable states sharing a charge can be corrected independently.
As we saw earlier, DefectSystem is an immutable snapshot: corrections are applied once, to copies of defect_species, at construction — the original objects passed in are never modified.
[20]:
shifted = DefectSystem(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume, temperature=500,
vbm_shift=0.05, cbm_shift=-0.1,
)
print(shifted)
DefectSystem
bandgap: 4.85 eV nelect: 20
volume: 46.1 ų temperature: 500 K
2 defect species:
v_Na [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = -1 E = 2.0000 eV (deg. 1)
v_Cl [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = +1 E = -1.0000 eV (deg. 1)
With the default rigid_shift=True, the band gap has shifted from 5 eV to 4.85 eV: this scissors the DOS, so the carrier concentrations and self-consistent Fermi level change, while the formation energies of v_Na’s and v_Cl’s charge states are unchanged. With rigid_shift=False, the same vbm_shift additionally shifts each variable charge state’s formation energy by charge * vbm_shift (and the original v_Na/v_Cl charge states are still untouched):
[21]:
nonrigid = DefectSystem(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume, temperature=500,
vbm_shift=0.05, cbm_shift=-0.1, rigid_shift=False,
)
print(nonrigid)
print("\noriginal v_Na charge-state energies:", [cs.energy for cs in v_Na.charge_states])
DefectSystem
bandgap: 4.85 eV nelect: 20
volume: 46.1 ų temperature: 500 K
2 defect species:
v_Na [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = -1 E = 1.9500 eV (deg. 1)
v_Cl [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = +1 E = -0.9500 eV (deg. 1)
original v_Na charge-state energies: [1.0, 2.0]
formation_energy_corrections applies a one-off correction to a specific DefectChargeState — for example, a +0.3 eV correction to v_Cl’s +1 state from a higher-level calculation:
[22]:
corrected = DefectSystem(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume, temperature=500,
formation_energy_corrections={v_Cl_1: 0.3},
)
print(corrected)
print("\noriginal v_Cl_1.energy:", v_Cl_1.energy)
DefectSystem
bandgap: 5 eV nelect: 20
volume: 46.1 ų temperature: 500 K
2 defect species:
v_Na [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = -1 E = 2.0000 eV (deg. 1)
v_Cl [nsites: 1]
q = +0 E = 1.0000 eV (deg. 1)
q = +1 E = -0.7000 eV (deg. 1)
original v_Cl_1.energy: -1.0
Temperature-dependent corrections with DefectSystemFactory¶
We’ve already used defect_system_factory.at(T) to get plain temperature snapshots. DefectSystemFactory can also take temperature-dependent shift functions: factory.at(T) evaluates vbm_shift_fn(T)/cbm_shift_fn(T)/formation_energy_correction_fns[cs](T) at T and returns an independent DefectSystem with those corrections applied.
As a toy example, suppose the band gap shrinks linearly with temperature (vbm_shift_fn(T) = 1e-4 * T, cbm_shift_fn(T) = -1e-4 * T), and (with rigid_shift=False) the defect levels are fixed in absolute energy while the bands move. We compare the self-consistent Fermi energy with and without these corrections across a temperature range:
[23]:
def vbm_shift_fn(temperature):
return 1e-4 * temperature
def cbm_shift_fn(temperature):
return -1e-4 * temperature
# v_Na and v_Cl exceed the 1% site-occupancy threshold at the high end of the sweep,
# so the DiluteLimitWarning is disabled for both toy factories below.
band_edge_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume,
vbm_shift_fn=vbm_shift_fn, cbm_shift_fn=cbm_shift_fn, rigid_shift=False,
occupancy_warning_threshold=None,
)
baseline_factory = DefectSystemFactory(
defect_species=[v_Na, v_Cl], dos=dos, volume=defect_system.volume,
occupancy_warning_threshold=None,
)
temperatures = np.linspace(300, 1500, 60)
corrected_ef = []
baseline_ef = []
for t in temperatures:
corrected_ef.append(band_edge_factory.at(t).get_sc_fermi()[0])
baseline_ef.append(baseline_factory.at(t).get_sc_fermi()[0])
plt.plot(temperatures, baseline_ef, label="no corrections")
plt.plot(temperatures, corrected_ef, label="with band-edge corrections", ls="--")
plt.xlabel("Temperature / K")
plt.ylabel("SC Fermi Energy / eV")
plt.legend()
plt.show()
[ ]: