The defect formation energy in the supercell approach

Within the supercell approach, the formation energy of a point defect d in charge state q can be expressed as the difference in grand potential that the formation of the defect in an otherwise pristine crystal entails [ZN91]:

(2)\[\Delta E_f(d; q) = G(d; q) - G(\mathrm{bulk}) - \sum_i n_i \mu_i + q \left(\epsilon_{VBM} + E_F\right) + E_{corr}\]

where \(\Delta E_f(d; q)\) represents the defect formation energy, \(G(d; q)\) and \(G(\mathrm{bulk})\) are the free energies of the defective and pristine supercells, respectively, \(n_i\) are the amount of atoms added or removed from the host material in order to create the defect, \(\mu_i\) are the chemical potentials of these atoms, \(\epsilon_{VBM}\) is the valence band maximum eigenvalue, \(E_F\) is the Fermi level, which can range within the material band gap, and \(E_{corr}\) is an energy term that corrects for finite-size effects originating from the use of the supercell approach. As a common approximation, instead of the free energies we will use the electronic energies, the main quantity calculated by first-principles codes.

Note

In the following tutorials, to illustrate how \(\Delta E_f(d; q)\) can be calculated with Spinney, we will use an example the formation energy of a B vacancy in the charge state -3 in cubic BN.

Read the quickstart guide for a succint tutorial explaining how to calculate the defect formation energy with Spinney using output files from VASP and WIEN2k.

Otherwise, read the In-depth Tutorial.

In-depth Tutorial

Initialize a PointDefect object

The easiest way to calculate \(\Delta E_f(d; q)\) in Spinney is by using the PointDefect class. To initialize a PointDefect object, we need a Atoms object of the ASE library representing the defective system. Such object must have attached a Calculator implementing a get_total_energy() method that returns the energy of the defective system.

Using VASP

ASE supports several first-principles codes; for example, if VASP is used, and the result of the calculation of the defective system is stored in the OUTCAR_def file, then an Atoms object suitable for initializing a PointDefect instance can be simply obtained using the ase.io.read() function. The PointDefect object can in this case be initialized with few lines of Python code:

from spinney.structures.pointdefect import PointDefect
import ase.io

# use ASE to read the OUTCAR file of the defective system
outcar = ase.io.read('OUTCAR_def', format='vasp-out')
# initialize a Spinney PointDefect object
pd = PointDefect(outcar)

Using WIEN2k or other codes not fully supported by ASE

While ASE can read the outputs of several first-principles codes, not all of them are supported. For example, up to the version 3.17, the last ASE version compatible with Spinney, there is no function that can read the .scf files produced by WIEN2k. In this case, one can write an ad hoc function serving this purpose.

Warning

When writing the helper functions needed to create an ASE Atoms object to be used for initialize a Spinney PointDefect instance from the output files of a first-principles suite not yet supported by ASE, you should convert lengths in Angstrom and energies in eV. This is the defalut in ASE and what is assumed by the PointDefect class.

In the case of WIEN2k, ASE can read the .struct files of WIEN2k and return the proper Atoms object. In this case create the ASE Atoms object required by PointDefect is easy.

Suppose you have calculated the defective supercell using WIEN2k and have the files defective.scf and defective.struct. In this case, the general way for initializing a PointDefect object is:

  1. read the total energy calculated by WIEN2k in the defective.scf file, for example using a grep command on Linux:

    grep :ENE defective.scf

  2. use ASE to read the defective.struct file and initialize an Atoms object:

    struct=ase.io.read('defective.struct', format='struct')

  3. use the DummyAseCalculator class to initialize a dummy calculator and use the method set_total_energy() to insert the calculated energy found in point 1 (the energy has to be converted into units of eV, which is the energy unit used in ASE),

    from spinney.structures.pointdefect import PointDefect, DummyAseCalculator
from spinney.constants import conversion_table
# convert the energy read from the .scf file to eV
energy = -5048.05765754*conversion_table['Ry']['eV']
calc = DummyAseCalculator(struct)
calc.set_total_energy(energy)

  4. attach the dummy calculator to the Atoms object created in step 2,

    struct.set_calculator(calc)

  5. Use this Atoms object to initialize a new PointDefect object.

    pd = PointDefect(struct)

This procedure can be generally used for any first-principles code for which ASE can create an Atoms object with the structural information but without information about the system energy.

For WIEN2k, Spinney offers some helping functions. The proper Atoms object that can be used to initialize a PointDefect object can be obtained using the function prepare_ase_atoms_wien2k() which takes as arguments the .struct and .scf files and returns the Atoms object with a dummy calculator.

from spinney.io.wien2k import prepare_ase_atoms_wien2k
ase_scf = prepare_ase_atoms_wien2k('defective.struct', 'defective.scf')
pd = PointDefect(ase_scf)

Similar helping functions can be easily written for any other ab initio code.

Calculating the defect formation energy

Once a new PointDefect instance has been created, we can feed it the data necessary to calculate the defect formation energy. Then needed data come from the very definition of \(\Delta E_f(d; q)\) of equation (2).

Feed the required data to the PointDefect instance

  1. Energy of the host material.

    \(G(\mathrm{bulk})\) is read from an ASE Atoms object representing the pristine system. Such object must thus have attached a calculator. The necessary steps are analogous to the one described in the previous section for the defective system. Once the instance has been created, we can feed it to the PointDefect object:

    # for example, if VASP was used
ase_pristine = ase.io.read('OUTCAR_pristine', format='vasp-out')
pd.set_pristine_system(ase_pristine)

  2. Chemical potential of the elements involved in the creation of the defect.

    Suppose we are interested in the defect formation energy in the B-rich limit. In order to find \(\mu_\mathrm{B}\), we need to calculate the reference state for the B atom. We can take it as the \(\alpha-\mathrm{B}\) phase in the rhombohedral crystal family. \(\mu_\mathrm{B}\) is then equal to the system energy per atom. \(\mu_\mathrm{N}\) is then simply obtained as \(\mu_\mathrm{N} = E_\mathrm{BN} - \mu_\mathrm{B}\). Suppose \(\mu_\mathrm{B}\) and \(\mu_\mathrm{N}\) have been stored in the variables chem_pot_B and chem_pot_N, respectively. Then we can inform the PointDefect object that we want to use these values of the chemical potentials by typing:

    pd.set_chemical_potential_values({'N':chem_pot_N, 'B':chem_pot_B})

  3. Information about the pristine system.

    From equation (2), we see that we also need \(\epsilon_{VBM}\) and \(E_F\). If the values are stored in the variables e_vbm and fermi_level, respectively, then we can feed them to the PointDefect object using:

    pd.set_vbm(e_vbm)
pd.set_fermi_level_value_from_vbm(fermi_level)

  4. Information about the defective system.

    This information is needed in order to compute \(E_{corr}\). This term mainly corrects for electrostatic finize-size effects due to the presence of charged defects. As other finize-size errors can generelly be made negligibly small by increasing the supercell size [FGH+14], the Spinney package implements two popular correction schemes for correcting electrostatic finite-size effects, as explained in detail in Correction schemes for electrostatic finite-size effects. Each of correction scheme requires some specific data in order to compute \(E_{corr}\); additionally, some data are required by every scheme.

    General data:

    • The defect charge state \(q\), q=-3 in our example:

      pd.set_defect_charge(q)

    • The defect position with respect to the supercell basis vectors, def_position=(0.5, 0.5, 0.5), if the vacancy was placed in the centre of the supercell:

      pd.set_defect_position(def_position)

    • The dielectric constant/tensor of the host material dielectric_tensor:

      pd.set_dielectric_tensor(dielectric_tensor)

    Correction-scheme-specific data:

    • First we need to inform the PointDefect instance about which correction scheme we intend to use. Spinney implements two state-of-the-art methods: the scheme proposed by Freysoldt, Neugebauer and Van de Walle, [FNVdW09] and the improved version proposed by Kumagai and Oba [KO14]. More information about such schemes can be found in Correction schemes for electrostatic finite-size effects.

      To specify the correction scheme of Freysoldt, Neugebauer and Van de Walle:

      pd.set_finite_size_correction_scheme('fnv')

      To specify the correction scheme of Kumagai and Oba:

      pd.set_finite_size_correction_scheme('ko')

    • Next, we need to add the required scheme-specific data. In this example we will consider the schem of Kumagai and Oba. We refer to Correction schemes for electrostatic finite-size effects for details about the scheme of Freysoldt, Neugebauer and Van de Walle.

      The method of Kumagai and Oba requires the electrostatic potential calculated at the ionic sites for the pristine and defective systems. Such data have to be stored in two arrays, potential_pristine, potential_defective. The way these data are obtained is specific of the employed first-principle package. We provide functions to extract such information from VASP and WIEN2k outputs.

      • VASP: the required output files are the OUTCAR files for the pristine and defective supercells, OUTCAR_pristine and OUTCAR_def, respectively.

        from spinney.io.vasp import extract_potential_at_core_vasp
potential_pristine = extract_potential_at_core_vasp('OUTCAR_pristine')
potential_defective = extract_potential_at_core_vasp('OUTCAR_def')

      • WIEN2k: the required output files are the .struct and .vcoul files for the pristine and defective systems: pristine.struct/.vcoul and defective.struct/.vcoul, respectively.

        from spinney.io.wien2k import extract_potential_at_core_wien2k
potential_pristine = extract_potential_at_core_wien2k('pristine.struct', 'pristine.vcoul')
potential_defective = extract_potential_at_core_wien2k('defective.struct', 'defective.vcoul')
# convert to eV
potential_pristine *= conversion_table['Ry']['eV']
potential_defective *= conversion_table['Ry']['eV']

      • Other first-principles codes: for other codes, one needs to write an output-specific helper function in order to extract the ionic-site electrostatic potential and initialize the arrays: potential_pristine and potential_defective.

      Once the arrays containing the electrostatic potential at the ionic sites have been created, we can feed the to the PointDefect instance:

      pd.add_correction_scheme_data(potential_pristine=potential_pristine,
                              potential_defective=potential_defective)

      the add_correction_scheme_data() takes correction-scheme-specific keyword arguments.

Obtaining the defect formation energy

Once the PointDefect object has been initializated and the needed data has been fed to it, we can obtain the defect formation energy.

# energy without adding corrections for electrostatic finite-size effects
uncorrected_energy = pd.get_defect_formation_energy()
# corrected energy
corrected_energy = pd.get_defect_formation_energy(True)

Quickstart Guide

Calculate the formation energy of charged point defects using Spinney is easy and fast. In this getting-started guide we will explain how to obtain this quantity for a Boron vacancy in charge state -3 in the B-rich limit. For correcting for electrostatic finite-size effects we apply the method of Kumagai and Oba [KO14].

We give explicit examples using two popular first-principles code:

Each ab initio code produces different outputs; however the information necessary to calculate the defect formation energy is generally present somewhere in these output files. To use Spinney with a different first-principle code, you will need to create some ad hoc helper functions. Read the whole guide for more information.

Using VASP

The snippet of code below shows how the formation energy of the defect can be calculated from VASP output files. It is assumed that the script working directory contains the following files:

  • OUTCAR files of the pristine and defective supercells: OUTCAR_pristine and OUTCAR_defective, respectively.

  • OUTCAR file for the reference state of bulk B: OUTCAR_B

Moreover, the following variables need to be created:

  • vbm: the valence-band-maximum eigenvalue of the host material.

  • e_r: the dielectric tensor/constant of the host material.

  • defect_position: the defect position in the supercell in fractional coordinates.

  • q: the defect charge state (-3 in this example).

from spinney.structures.pointdefect import PointDefect
from spinney.io.vasp import extract_potential_at_core_vasp
from spinney.tools.formulas import count_elements

import ase.io

### initialize a point defect object
ase_defective = ase.io.read('OUTCAR_defective', format='vasp-out')
pd = PointDefect(ase_defective)

### Feed it the data
ase_pristine = ase.io.read('OUTCAR_pristine', format='vasp-out')
pd.set_pristine_system(ase_pristine)
# get the chemical potential of Boron
ase_boron = ase.io.read('OUTCAR_B', format='vasp-out')
chem_pot_B = ase_boron.get_total_energy()/ase_boron.get_number_of_atoms()
# get the chemical potential of Nitrogen in the B-rich conditions
elements = count_elements(ase_pristine.get_chemical_formula())
chem_pot_N = ase_pristine.get_total_energy() - elements['B']*chem_pot_B
chem_pot_N /= elements['N']
pd.set_chemical_potential_values({'N':chem_pot_N, 'B':chem_pot_B})
# set valence band maximum, the Fermi level is set by default to zero
pd.set_vbm(vbm)
# if one wants a Fermi level not equal to 0, uncomment:
# pd.set_fermi_level_value_from_vbm(E_F)

# data for the correction scheme
pd.pd.set_defect_charge(q)
pd.set_defect_position(defect_position)
pd.set_dielectric_tensor(e_r)
# scheme of Kumagai and Oba
pd.set_finite_size_correction_scheme('ko')
# get the required data and feed them to the instance
pot_pristine = extract_potential_at_core_vasp('OUTCAR_pristine')
pot_defective = extract_potential_at_core_vasp('OUTCAR_defective')
pd.add_correction_scheme_data(potential_pristine=pot_pristine,
                              potential_defective=pot_defective)

### calculate the defect formation energy
uncorrected_energy = pd.get_defect_formation_energy()
corrected_energy = pd.get_defect_formation_energy(True)

Using WIEN2k

The snippet of code below shows how the formation energy of the defect can be calculated from WIEN2k output files. It is assumed that the script working directory contains the following files:

  • .struct, .scf and .vcoul files of the pristine and defective supercells: pristine.struct/.scf/.vcoul and defective.struct/.scf/.vcoul, respectively.

  • .struct and .scf files for the reference state of bulk B: boron.struct/.scf

Moreover, the following variables need to be created:

  • vbm: the valence-band-maximum eigenvalue of the host material.

  • e_r: the dielectric tensor/constant of the host material.

  • defect_position: the defect position in the supercell in fractional coordinates.

  • q: the defect charge state (-3 in this example).

from spinney.structures.pointdefect import PointDefect
from spinney.io.wien2k import prepare_ase_atoms_wien2k
from spinney.io.wien2k import extract_potential_at_core_wien2k
from spinney.tools.formulas import count_elements
from spinney.constants import conversion_table

### initialize a point defect object
ase_defective = prepare_ase_atoms_wien2k('defective.struct', 'defective.scf')
pd = PointDefect(ase_defective)

### Feed it the data
ase_pristine = prepare_ase_atoms_wien2k('pristine.struct', 'pristine.scf')
pd.set_pristine_system(ase_pristine)
# get the chemical potential of Boron
ase_boron = prepare_ase_atoms_wien2k('boron.struct', 'boron.scf')
chem_pot_B = ase_boron.get_total_energy()/ase_boron.get_number_of_atoms()
# get the chemical potential of Nitrogen in the B-rich conditions
elements = count_elements(ase_pristine.get_chemical_formula())
chem_pot_N = ase_pristine.get_total_energy() - elements['B']*chem_pot_B
chem_pot_N /= elements['N']
pd.set_chemical_potential_values({'N':chem_pot_N, 'B':chem_pot_B})
# set valence band maximum, the Fermi level is set by default to zero
pd.set_vbm(vbm)
# if one wants a Fermi level not equal to 0, uncomment:
# pd.set_fermi_level_value_from_vbm(E_F)

# data for the correction scheme
pd.pd.set_defect_charge(q)
pd.set_defect_position(defect_position)
pd.set_dielectric_tensor(e_r)
# scheme of Kumagai and Oba
pd.set_finite_size_correction_scheme('ko')
# get the required data and feed them to the instance
pot_pristine = extract_potential_at_core_wien2k('pristine.struct',
                                                'pristine.vcoul')
pot_defective = extract_potential_at_core_wien2k('defective.struct',
                                                 'defective.vcoul')
# convert to eV
pot_pristine *= conversion_table['Ry']['eV']
pot_defective *= conversion_table['Ry']['eV']
pd.add_correction_scheme_data(potential_pristine=pot_pristine,
                              potential_defective=pot_defective)

### calculate the defect formation energy
uncorrected_energy = pd.get_defect_formation_energy()
corrected_energy = pd.get_defect_formation_energy(True)

Manage a defective system with the class DefectiveSystem

One is usually interested in the study of several point defects in a defective system. In this case, one could insert the code snippets shown above into a loop and create a PointDefect instance for each point defect.

This however might result in too much visual noise. Spinney offers the class DefectiveSystem to manage several point defects in the same host material. To initialize a new instance, one needs to specify the directory containing the results of the first-principles calculations (data_path) and the code used to produce them. The class expects this kind of directory tree structure:

"data_path"
 ├── data_defects
 │   ├── "defect_name"
 │   │   ├── "charge_state"
 │   │   │   └── "files"
 │   │   ├── "charge_state"
 │   │   │   └── "files"
 │   │   ├── "charge_state"
 │   │   │   └── "files"
 │   │   └── "charge_state"
 │   │       └── "files"
 │   ├── "defect_name"
 ...
 └── pristine
     └── "files"

  • the names data_defect and pristine are directories and are mandatory, as the class will look in these places in order to find the required information about the defective and pristine systems.

  • files are the output of the first-principles calculations required to calculate the defect formation energy. The number and type of required files depends on the employed first-principles code and on whether one wants to apply finite-size-effects corrections.

    • For VASP: one needs at least the OUTCAR file. If the correction scheme fnv has to be used, then also the LOCPOT file must be present.

    • For WIEN2k: one needs at least the case.struct and case.scf files. If the correction scheme ko has to be used, then also the case.vcoul file must be present. The fnv method has not yet been implemented in WIEN2k.

    • In all cases, if any correction scheme has to be used, one needs to add a file named position.txt with the fractional coordinates of the defective site on a single line, each entry separated by a white space. For example:

    0.5 0.5 0.5

  • charge_state are directories and must be integers representing the defect charge state. The value specify by charge_state will be considered to be the charge state of the defect.

  • defect_name are directories, they can have any name that can describe the defective system.

  • Any other directories can appears below data_path, these will be ignored by the class.

We can take as an example the study of the defect chemistry of intrinsic GaN using the code VASP. In this case the data_path directory tree would look like this:

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

We can now initialize a new instance of DefectiveSystem:

defective_system = DefectiveSystem('data', 'vasp')

Similarly to the initialization of a PointDefect object, in order to calculate the defect formation energies we need some data relative to the pristine system and the chemical potential values. These can be assigned as object attributes:

# calculate chemical potential values in the Ga-rich limit
path_pristine = os.path.join('data', 'pristine', 'OUTCAR')
path_ga = os.path.join('data', 'Ga', 'OUTCAR')
ase_ga = ase.io.read(path_ga, format='vasp-out')
ase_pristine = ase.io.read(path_pristine, format='vasp-out')
chem_pot_ga = ase_ga.get_total_energy()/ase_ga.get_number_of_atoms()
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']

# valence band maximum and dielectric tensor
vbm = 5.009256
e_rx = 5.888338 + 4.544304
e_rz = 6.074446 + 5.501630
e_r = [[e_rx, 0, 0], [0, e_rx, 0], [0, 0, e_rz]]

# add the data
defective_system.vbm = vbm
defective_system.dielectric_tensor = e_r
defective_system.chemical_potentials = {'Ga':chem_pot_ga, 'N':chem_pot_n}

# use the correction scheme of Kumagai and Oba
defective_system.correction_scheme = 'ko'
# calculate defect formation energies for each point defect
# and print output to the screen
defective_system.calculate_energies(verbose=True)

The calculated defect formation energy are stored in a Pandas dataframe and can be accessed from the attribute data:

df = defective_system.data
printf(df)

will output:

               Form Ene (eV)  Form Ene Corr. (eV)
Defect Charge
N_Ga    2           7.022139             7.272291
        3           7.435422             8.143598
       -1           9.977779            10.237273
        0           8.271712             8.271712
        1           7.461967             7.464786
Ga_int  2           3.190741             3.956642
        3           1.487580             3.063845
        0           8.339166             8.339166
        1           5.710694             5.964960
...

they can be written to a text file, in a format readable by other modules of Spinney by calling write_formation_energies():

defective_system.write_formation_energies('formation_energies_GaN_Ga_rich.txt')

The PointDefect objects for each point defect processed by the DefectiveSystem instance are accessible in a list through the attribute point_defects:

for pdf in defective_system.point_defects:
    print(pdf.my_name)

will print:

N_Ga 2
N_Ga 3
N_Ga -1
N_Ga 0
N_Ga 1
Ga_int 2
Ga_int 3
Ga_int 0
Ga_int 1
...