Charge Transition Levels¶
In section Thermodynamic limits for the chemical potentials we have discussed how Spinney can be used for obtaining the values of the chemical potentials in given conditions.
The defect formation energy, equation (2), depends also on the chemical potential of the electron, usually expressed as \(\epsilon_{VMB} + E_F\).
For this reason, the defect formation energy is often plotted as a function of
the Fermi level \(E_F\). Spinney offers the class
Diagram for this purpose.
Such diagrams are also useful for visualizing the defect thermodynamic charge transition levels, which define the value of \(E_F\) for which two defect charge states have the same formation energy.
We can take again the defect chemistry of intrinsic GaN as an illustrative example. In particular, the electronic energy of the intrinsic defects of GaN in several charge states was calculated using a supercell containing 96 atoms and the PBE exchange-correlation functional using a 4x2x2 reciprocal-point mesh. The formation energies will be calculated in the Ga-rich limit.
Contents
Using the class Diagram¶
Step 1: Calculate the Defect Formation Energies¶
The easiest way to calculate defect formation energies, including electrostatic finite-size effect corrections,
is through the PointDefect class. Since many defects in several charge states
have to be considered, a straightforward way for processing the DFT data is collecting the results in a proper directory hierarchy.
For example, using output from VASP, if we decide to use the correction scheme of Kumagai and Oba, which only requires the OUTCAR files, the
directory tree might look like this (see also Manage a defective system with the class DefectiveSystem):
data_path
├── data_defects
│ ├── Ga_int
│ │ ├── 0
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ ├── 1
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ ├── 2
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ └── 3
│ │ ├── OUTCAR
│ │ └── position.txt
│ ├── Ga_N
│ │ ├── 0
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ ├── 1
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ ├── -1
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ ├── 2
│ │ │ ├── OUTCAR
│ │ │ └── position.txt
│ │ └── 3
│ │ ├── OUTCAR
│ │ └── position.txt
...
├── Ga
│ └── OUTCAR
├── N2
│ └── OUTCAR
└── pristine
└── OUTCAR
where the calculations result are collected for each charge state of any specific defect, as well for the pristine system and parent compounds.
Now the defect formation energy can be calculated for each defect by walking the directory tree. This code snippet also writes a file containing the formation energies for each defect in each charge state.
import os
import ase.io
from spinney.structures.pointdefect import PointDefect
from spinney.io.vasp import extract_potential_at_core_vasp
from spinney.tools.formulas import count_elements
path_defects = os.path.join('data', 'data_defects')
path_pristine = os.path.join('data', 'pristine', 'OUTCAR')
path_ga = os.path.join('data', 'Ga', 'OUTCAR')
# prepare preliminary data
ase_pristine = ase.io.read(path_pristine, format='vasp-out')
pot_pristine = extract_potential_at_core_vasp(path_pristine)
ase_ga = ase.io.read(path_ga, format='vasp-out')
vbm = 5.009256 # valence band maximum pristine system
e_rx = 5.888338 + 4.544304
e_rz = 6.074446 + 5.501630
# dielectric tensor
e_r = [[e_rx, 0, 0], [0, e_rx, 0], [0, 0, e_rz]]
# store defect positions in dictionary for access convenience
defect_positions = {'Ga_int' : (0.25, 0.55, 0.55),
'N_int' : (0.5, 0.41667, 0.50156),
'Vac_Ga' : (0.5, 0.6666666666, 0.50156),
'Vac_N' : (0.5, 0.3333333333, 0.46045),
'Ga_N' : (0.5, 0.3333333333, 0.46045),
'N_Ga' : (0.5, 0.6666666666, 0.50156)
}
# get the chemical potential of Ga
chem_pot_ga = ase_ga.get_total_energy()/ase_ga.get_number_of_atoms()
# get the chemical potential of N in the Ga-rich conditions
elements = count_elements(ase_pristine.get_chemical_formula())
chem_pot_n = ase_pristine.get_total_energy() - elements['Ga']*chem_pot_ga
chem_pot_n /= elements['N']
# write the formation energy for Fermi level = 0 to a file
energy_file = open('formation_energies_GaN_Ga_rich.txt', 'w')
energy_file.write('#{:<8} {:8} {:>30}\n'.format('system',
'charge',
'formation energy (eV)'))
for root, dirs, files in os.walk(path_defects):
path = root.split(os.sep)
if 'OUTCAR' in files:
def_path = os.path.join(root, 'OUTCAR')
ase_outcar = ase.io.read(def_path)
pot_defective = extract_potential_at_core_vasp(def_path)
defect_name = path[-2]
charge_state = int(path[-1])
print('Processing defect {} in charge state {}'.format(defect_name,
charge_state))
print('Data in: {}'.format(path))
# prepare PointDefect object
pdf = PointDefect(ase_outcar)
pdf.set_pristine_system(ase_pristine)
pdf.set_chemical_potential_values({'N':chem_pot_n, 'Ga':chem_pot_ga})
pdf.set_vbm(vbm)
pdf.set_defect_charge(charge_state)
pdf.set_defect_position(defect_positions[defect_name])
pdf.set_dielectric_tensor(e_r)
pdf.set_finite_size_correction_scheme('ko')
pdf.add_correction_scheme_data(potential_pristine=pot_pristine,
potential_defective=pot_defective)
corrected_energy = pdf.get_defect_formation_energy(True)
energy_file.write('{:<8} {:>8} {:30.10f}\n'.format(defect_name,
charge_state,
corrected_energy))
print('Done.\n------------------------------------------------------')
energy_file.close()
A completely similar approach can be used for WIEN2k or any other first-principle code.
Step 2: Get Thermodynamic Charge Transition Levels and plot the Diagram¶
The file formation_energies_GaN_Ga_rich.txt, containing the formation energies calculated for each point defect, has the general format of a text file, with optional headers starting with a #:
# description or other comments
defect_name charge_state formation_energy_at_vbm
formation_energy_at_vbm is the calculated value of the formation energy for an electron chemical potential equal to the valence band maximum, i.e. for \(E_F = 0\).
Thermodynamic charge transition levels can be easily calculated by initializing a
Diagram object.
The class constructor takes two mandatory parameters: a dictionary containing information about point defects formation energies (at \(E_F = 0\)) and their charge state and a 2-ple indicating the system valence band maximum and conduction band minimum. The latter argument will give the range within which possible transition levels are considered. Absolute values are intended for \(\mu_e\) as independent variables; setting the valence band maximum to 0, will give \(E_F\) as independent variable.
The former argument, might be obtained from a file analogous to formation_energies_GaN_Ga_rich.txt using the
function extract_formation_energies_from_file().
The following code snipper will calculate the charge transition levels and print them in the text file transition_levels.txt.
from spinney.defects.diagrams import extract_formation_energies_from_file, Diagram
data_file = 'formation_energies_GaN_Ga_rich.txt'
defect_dictionary = extract_formation_energies_from_file(data_file)
dgm = Diagram(defect_dictionary, (0, 1.713)) # E_F from 0 to PBE band gap value
dgm.write_transition_levels('transition_levels.txt')
The content of transition_levels.txt will be:
#Defect type q/q'
Ga_N 2/3 0.441121
1/2 0.731836
0/1 1.250109
Ga_int 2/3 0.892796
N_Ga 1/2 0.192495
0/1 0.806926
...
Transition level values can also be accessed in a Python session using the attribute
transition_levels:
dgm.transition_levels
returns a Pandas Series object with the transition levels.
Extending the band gap¶
It is well known that standard local and semilocal functionals severely underestimate the band gap of semiconductor materials.
For this purpose the valence band maxima are often aligned to those obtained from more accurate functionals, like hybrids.
Transition levels can be calculated in this case by specifying an extended_gap parameter in the constructor of
Diagram.
For example, the paper if reference [LVdW17] finds that the valence band maximum of GaN calculated with HSE with the 31% of Hartree-Fock exchange lies 0.85 eV below the valence band maximum found by PBE and predicts a band gap of 3.51 eV.
We can calculate the transition levels in this extended gap by using the following snippet:
# Fermi level will range in the HSE range
dgm = Diagram(defect_dictionary, (0, 1.713), (-0.85, - 0.85 + 3.51))
dgm.write_transition_levels('transition_levels_extended.txt')
The file transition_levels_extended.txt will contain also transition levels located outside the PBE gap:
#Defect type q/q'
Ga_N 2/3 0.441121
1/2 0.731836
0/1 1.250109
-1/0 1.825432
Ga_int 2/3 0.892796
1/2 2.008318
0/1 2.374206
...
A plot is very useful to visualize these results. The following snippet shows how to make the plot with Spinney.
from spinney.defects.diagrams import extract_formation_energies_from_file, Diagram
data_file = 'formation_energies_GaN_Ga_rich.txt'
defect_dictionary = extract_formation_energies_from_file(data_file)
dgm = Diagram(defect_dictionary, (0, 1.713), (-0.85, - 0.85 + 3.51))
# use some prettier labels in the plot
dgm.labels = {'Ga_int' : r'$Ga_i$',
'N_int' : r'$N_i$',
'Vac_Ga' : r'$Vac_{Ga}$',
'Vac_N' : r'$Vac_N$',
'Ga_N' : r'$Ga_N$',
'N_Ga' : r'$N_{Ga}$'}
# personalize colors to use in the plot
colors = {'Ga_int' : 'red',
'N_int' : 'blue',
'Vac_Ga' : 'orange',
'Vac_N' : 'gray',
'Ga_N' : 'magenta',
'N_Ga' : 'green'}
# save the plot
dgm.plot(save_flag=True, save_title='diagram',
title='Intrinsic GaN, Ga-rich conditions', legend=True,
colors_dict=colors, x_label=r'$E_F$ (eV)')
This code will save the file diagram.pdf with the plot:
Using the class DefectiveSystem¶
The directories structure used in the previous section is compatible with the
one required by an instance of DefectiveSystem.
Such object has a Diagram instance accessible through
the attribute diagram. In order to create this Diagram object,
one needs to specify the parameters gap_range and eventually extended_gap_range
needed by Diagram.
These values can be added using the attributes gap_range and extended_gap_range of
a DefectiveSystem.
For example, suppose that a DefectiveSystem
instance as been initialized as defective_system (see Manage a defective system with the class DefectiveSystem). Then for our example,
one can type:
defective_system.gap_range = (0, 1.713)
defective_system.extended_gap_range = (-0.85, - 0.85 + 3.51)
dgm = defective_system.diagram
dgm is an initialized Diagram object, which we can use as
shown in the previous sections.