Tutorial of pydefect

This page explains how to use the pydefect code.

Note 1: Pydefect currently supports only the Vienna Ab-initio Simulation Package (VASP). Therefore, it assumes the use of VASP-specific input and output file names (e.g., POSCAR, POTCAR, OUTCAR) as well as computational methods such as periodic boundary conditions.

Note 2: The units used in Pydefect follow the VASP convention: energy is in electronvolts (eV) and length is in angstroms (Å).

Note 3: At present, only nonmagnetic host materials are supported.

The workflow for point-defect calculations in non-metallic solids is illustrated below. Some tasks can be performed in parallel, while others must be carried out sequentially. Since these procedures are often complex and time-consuming, researchers are prone to making mistakes.

The primary goal of pydefect is to provide an automated system for point-defect calculations in non-metallic materials. This helps researchers save time and minimize human error.

_images/pydefect.png

This cheat sheet is divided into three parts: unitcell, chemical potential diagram (cpd), and defect, each corresponding to a working directory, as illustrated below.

Red text indicates related commands, while blue text denotes the files generated by those commands. Branches represent alternative workflows or options. Filled blue rectangles indicate steps that require VASP calculations. Bolded rectangles represent steps that involve multiple directories related to the considered defects. Green rectangles indicate optional steps.

Here, we assume the following directory structure. The placeholder <project_name> typically represents the name of the target material, optionally including its crystal structure — for example, rutile-TiO:sub:2``.

<project_name>
 │
 ├ pydefect.yaml
 ├ vise.yaml
 │
 ├ unitcell/ ── structure_opt/
 │            ├ band/
 │            ├ dielectric/
 │            └ dos/
 │
 ├ cpd/ ──── <competing_phase 1>
 │       ├── <competing_phase 2>
 │           ....
 │
 └ defect/ ── perfect/
            ├─ Va_X_0/
            ├─ Va_X_1/
            ├─ Va_X_2/
             ...

We recommend that users follow the same directory structure if possible. The details of the workflow are explained step by step, using an example of MgAl₂O₄ calculated with the PBEsol functional.

In :code:pydefect, there are five main command-line tools:

pydefect pydefect_vasp pydefect_util pydefect_vasp_util pydefect_print The commands pydefect_vasp and pydefect_vasp_util are related to VASP, while pydefect and pydefect_util are independent of any specific DFT code. Commands with the suffix _util are not strictly necessary, but they can be useful in certain situations.

Each subcommand supports a variety of arguments. You can always display the help message using the -h flag. For example:

pydefect s -h

to view detailed information.

1. Relaxation of the unit cell

Point-defect calculations are typically performed on a theoretically relaxed structure, using a specified exchange-correlation functional and projector augmented wave (PAW) potentials. This ensures that the structure is free from artificial strain and stress, which can otherwise introduce unwanted dependencies on the supercell size. Therefore, the first step is usually to optimize the lattice constants and the fractional atomic coordinates in the unit cell.

We first prepare the POSCAR file of the pristine bulk unit cell, and create the unitcell/ directory and its subdirectory unitcell/structure_opt/ (using mkdir -p unitcell/structure_opt/), then move into it. (In this tutorial, any name ending with / indicates a directory.) When pydefect needs to generate VASP input files, namely INCAR, POTCAR, and KPOINTS, we use the vise (short for vasp integrated supporting environment), which creates input files for various tasks and exchange-correlation (XC) functionals. vise depends on pymatgen, and therefore, as explained in pymatgen usage page 1 or page 2, you need to set the PMG_DEFAULT_FUNCTIONAL and PMG_VASP_PSP_DIR in the .pmgrc.yaml file located in your home directory. For example:

PMG_DEFAULT_FUNCTIONAL: PBE_54
PMG_MAPI_KEY: xxxxxxxxxxxxxxxx
PMG_VASP_PSP_DIR: /home/kumagai/potcars/

In pydefect, PMG_MAPI_KEY is required to query POSCAR files and the total energies of competing materials.

Input files for optimizing a unit cell using the PBEsol functional can be generated with the following command:

vise vasp_set -x pbesol

Here, vasp_set, or its abbreviation vs, is a subcommand of the main vise function. Similarly, all subcommands in pydefect and vise have their own abbreviations.

Since the PBE functional is the default in vise, we use the -x argument to switch the XC functional to PBEsol.

Note that structure optimization generally needs to be repeated with a 1.3× larger cutoff energy until both the forces and stresses are well converged at the first ionic step. See the vasp manual for further details.

This iterative process is not supported by pydefect, but it can easily be implemented using a simple shell script.

In pydefect, users can control various default parameters via the vise.yaml file used by vise. See the vise documentation for more information. The set of configurable parameters can be found in defaults.py in pydefect.

2. Calculation of band, DOS, and dielectric tensor

We then calculate the band structure (BS), density of states (DOS), and dielectric constants. In defect calculations, the BS is used to determine the valence band maximum (VBM) and conduction band minimum (CBM), while the dielectric constant, defined as the sum of the electronic (ion-clamped) and ionic dielectric tensors, is required for correcting the defect formation energies.

First, create the band/, dos/, and dielectric/ directories inside unitcell/, copy the POSCAR file from unitcell/structure_opt/, and run the following command in each directory:

vise vs -x pbesol -t band -pd ../structure_opt
vise vs -x pbesol -t dos -pd ../structure_opt -uis LVTOT True LAECHG True KPAR 1
vise vs -x pbesol -t dielectric_dfpt -pd ../structure_opt

The additional user_incar_settings (abbreviated as uis) in the dos directory are used to generate volumetric data for the electrostatic potential and all-electron charge density.

vise also provides plotting utilities for both the band structure (BS) and density of states (DOS). See the vise documentation for details.

3. Gathering unit cell information for point-defect calculations

We next collect the bulk information, namely the band edges and electronic and ionic dielectric tensors using the unitcell (= u) sub-command.

pydefect_vasp u -vb band/vasprun-finish.xml -ob band/OUTCAR-finish -odc dielectric/OUTCAR-finish -odi dielectric/OUTCAR-finish -n MgAl2O4

Here, the electronic and ionic dielectric constants can be set with different OUTCAR files. Then, unitcell.yaml is generated, which will be used for analyzing defect calculations later.

system: MgAl2O4
vbm: 4.0183
cbm: 9.2376
ele_dielectric_const:
- - 3.075988
  - 0.0
  - -0.0
- - 0.0
  - 3.075988
  - 0.0
- - -0.0
  - -0.0
  - 3.075988
ion_dielectric_const:
- - 5.042937
  - -0.0
  - -0.0
- - -0.0
  - 5.042937
  - 0.0
- - -0.0
  - 0.0
  - 5.042937

Of course, the users can also create it by hand.

4. Calculation of competing phases

When a defect is introduced, atoms are exchanged with hypothetical atomic reservoirs within a thermodynamic framework. To calculate the free energy of defect formation, often approximated by the defect formation energy, we need to determine the chemical potentials of the atoms involved in the defect. Typically, these chemical potentials are chosen under the condition where the host material coexists with its competing phases, as determined from a chemical potential diagram (CPD).

To set up this process, we create directories under cpd/. We retrieve POSCAR files of stable or slightly unstable competing phases from the Materials Project database (MPD). This requires an MP API key, as described earlier.

For example, to obtain competing phases for MgAl2O4 with energies above the convex hull less than 0.5 meV/atom, use:

pydefect_vasp mp -e Mg Al O --e_above_hull 0.0005

This command generates directories like:

Al2O3_mp-1143/  Al_mp-134/  Mg149Al_mp-1185596/  Mg17Al12_mp-2151/
MgAl2O4_mp-3536/  MgAl2_mp-1094116/  MgO_mp-1265/  Mg_mp-1056702/  mol_O2/

To reduce computational time in this tutorial, we remove Mg149Al_mp-1185596/.

Each directory contains a POSCAR file and a prior_info.yaml file. The prior_info.yaml file includes information retrieved from the MPD that is useful for setting VASP calculation conditions via vise.

For example, Mg_mp-1056702/prior_info.yaml contains:

band_gap: 0.0
data_source: mp-1056702
total_magnetization: 0.0007357

indicating that Mg is a non-magnetic metal. vise uses this file to determine the k-point density in KPOINTS and whether spin polarization is needed via the ISPIN tag in INCAR. If any value is clearly incorrect, the user can modify it manually.

The molecules O2, H2, N2, NH3, and NO2 are not retrieved from MPD but are instead generated by pydefect, since MPD often treats them as solids, which is unsuitable for defect calculations.

We then generate INCAR, POTCAR, and KPOINTS files for these systems. A common cutoff energy, ENCUT, is needed

5. Construction of a supercell and defect species

We have finished the calculations for the unit cell and the CPD, and are now ready to proceed with point-defect calculations. First, create the defect/ directory.

Next, we generate files related to the supercell and defect species using the supercell (abbreviated as s) and defect_set (abbreviated as ds) subcommands.

pydefect recommends using a nearly isotropic (and, if possible, cubic-like) supercell with a moderate number of atoms. The following command creates a SPOSCAR file:

pydefect s -p ../unitcell/structure_opt/CONTCAR-finish

If the input structure is not a standardized primitive cell, a NotPrimitiveError will be raised.

Currently, pydefect constructs the supercell by expanding the conventional unit cell. In principle, it is possible to construct supercells with lattice angles that differ from those of the conventional cell— for example, one could make a supercell with mutually orthogonal a-, b-, and c-axes for hexagonal systems.

However, this is generally discouraged, because such modifications break the original symmetry, which can reduce the accuracy of point-defect calculations and make symmetry analysis at the defect site more difficult. Therefore, pydefect expands the lattice vectors along their original directions by default.

One exception is the tetragonal system, where rotated supercells (e.g., rotated by 45 degrees in the basal plane) can preserve the original symmetry and are therefore acceptable.

In pydefect, users can also specify the supercell transformation matrix. For example:

pydefect s -p ../unitcell/structure_opt/CONTCAR-finish --matrix 2 1 1

This matrix is applied to the conventional unit cell.

Note again that using an anisotropic supercell may change the symmetry of the structure, which can result in incorrect symmetry analysis in subsequent steps.

To view the conventional unit cell, use the following command:

vise si -p ../unitcell/structure_opt/CONTCAR-finish -c

Refer to the help message for more options and details.

Since JSON files are generally less human-readable than YAML files, pydefect provides the pydefect_print command to extract and display human-readable information from JSON files. For example:

pydefect_print supercell_info.json

This prints a summary of the supercell_info.json file as shown below:

Space group: Fd-3m
Transformation matrix: [-1, 1, 1]  [1, -1, 1]  [1, 1, -1]
Cell multiplicity: 4

   Irreducible element: Mg1
        Wyckoff letter: b
         Site symmetry: -43m
         Cutoff radius: 2.528
          Coordination: {'O': [1.94, 1.94, 1.94, 1.94]}
      Equivalent atoms: 0..7
Fractional coordinates: 0.2500000  0.2500000  0.2500000
     Electronegativity: 1.31
       Oxidation state: 2

   Irreducible element: Al1
        Wyckoff letter: c
         Site symmetry: .-3m
         Cutoff radius: 2.5
          Coordination: {'O': [1.92, 1.92, 1.92, 1.92, 1.92, 1.92]}
      Equivalent atoms: 8..23
Fractional coordinates: 0.6250000  0.3750000  0.3750000
     Electronegativity: 1.61
       Oxidation state: 3

   Irreducible element: O1
        Wyckoff letter: e
         Site symmetry: .3m
         Cutoff radius: 2.5
          Coordination: {'Mg': [1.94], 'Al': [1.92, 1.92, 1.92]}
      Equivalent atoms: 24..55
Fractional coordinates: 0.8614957  0.8614957  0.8614957
     Electronegativity: 3.44
       Oxidation state: -2

Using the defect_set (abbreviated as ds) subcommand, we can generate the defect_in.yaml file. An example of defect_in.yaml for MgAl2O4 is shown below:

Al_Mg1: [-1, 0, 1]
Al_O1: [-1, 0, 1, 2, 3, 4, 5]
Mg_Al1: [-1, 0, 1]
O_Al1: [-5, -4, -3, -2, -1, 0, 1]
Va_Al1: [-3, -2, -1, 0, 1]
Va_Mg1: [-2, -1, 0]
Va_O1: [0, 1, 2]

This file lists the combinations of defect species and their charge states. It can be edited manually if needed, or modified using the --keywords option.

To add dopants, for example Ca, use:

pydefect ds -d Ca

Some tips related to supercell_info.json and defect_in.yaml are provided below:

  1. Antisite and substitutional defects are determined based on the difference in electronegativity between the replacing and replaced atoms. The default threshold for the maximum allowed difference is defined in defaults.py, but it can be changed via pydefect.yaml. Since pydefect uses DefaultsBase from vise, the configuration rules follow those in vise.yaml (with different keywords).

  2. Oxidation states determine the possible defect charge states. For example, a vacancy of Sn2+ may have charge states of 0, -1, or -2, while Sn4+ may have 0 through -4. For antisite and substitutional defects, pydefect considers all possible vacancy and interstitial combinations. For example, Sn2+ on an S2- site may yield charge states of 0, +1, +2, +3, and +4.

    Oxidation states are automatically inferred using the oxi_state_guesses method of the Composition class in pymatgen. However, users may also set them manually. For example:

    pydefect ds --oxi_states Mg 4
    

    Note that the recommended charge states may sometimes miss important cases. For example, Zn vacancies in ZnO are known to exhibit a +1 charge state due to the formation of polarons at neighboring O sites. See Frodason et al., Phys. Rev. B (2017). Such cases must be added manually by the user.

  3. By default, atomic positions near the defect are perturbed to lower the symmetry to P1. While this can help stabilize defect configurations, it also increases the number of irreducible k-points. To suppress this, set displace_distance to 0 in pydefect.yaml.

  4. To restrict calculations to specific defect types, such as oxygen vacancies only, use the -k option with a Python regular expression. For example:

    pydefect ds -k "Va_O[0-9]?_[0-9]+"
    

    This creates the following directories:

    perfect/  Va_O1_0/  Va_O1_1/  Va_O1_2/
    

    For information on regular expressions, see the Python re module documentation.

6. Decision of interstitial sites

In addition to vacancies and substitutional defects, interstitials may also be important. While most people identify interstitial sites by visually inspecting the host crystal structure, there are some methods that attempt to systematically predict them.

However, identifying the most likely interstitial sites is generally difficult, as it depends strongly on the specific combination of the host and the inserted elements.

For example, the largest vacant space is often a favorable interstitial site for positively charged, closed-shell cations (e.g., Mg2+, Al3+), which tend not to form strong bonds with neighboring atoms.

In contrast, a proton (H+) tends to prefer positions near O2- or N3- to form strong O–H or N–H bonds. Conversely, a hydride ion (H-) is expected to occupy very different sites.

Therefore, careful consideration is required when selecting interstitial positions.

pydefect includes a utility to assist in identifying interstitial sites based on volumetric data, such as the all-electron charge density of the unit cell. This feature uses the ChargeDensityAnalyzer class from pymatgen.

To use it, you need to generate volumetric data such as AECCAR and LOCPOT based on the standardized primitive cell. This step has already been performed during the DOS calculation in this tutorial.

Note that it is generally inappropriate to use band structure calculation outputs for this purpose, since the cell may not match the standardized primitive cell.

After completing the VASP calculation, run the following command in the directory containing AECCAR{0,2}:

pydefect_vasp le -v AECCAR{0,2} -i all_electron_charge

This will print local minima in the charge density as follows:

       a      b      c      value  ave_value
0  0.125  0.125  0.125   3.387930   0.029844
1  0.625  0.125  0.125   3.383668   0.029845
2  0.125  0.625  0.125   3.383668   0.029845
3  0.125  0.125  0.625   3.383445   0.029845
4  0.500  0.500  0.500  16.501119   0.155178
5  0.750  0.750  0.750  16.501119   0.155178

More details are stored in volumetric_data_local_extrema.json, which can be viewed using:

pydefect_print volumetric_data_local_extrema.json

Example output:

info: all_electron_charge
min_or_max: min
extrema_points:
#  site_sym  coordination                                                     frac_coords               quantity
1  -3m       {'Mg': [1.75, 1.75], 'O': [2.14, 2.14, 2.14, 2.14, 2.14, 2.14]}  ( 0.125,  0.125,  0.125)  0.03
2  -43m      {'Al': [1.75, 1.75, 1.75, 1.75], 'O': [1.57, 1.57, 1.57, 1.57]}  ( 0.500,  0.500,  0.500)  0.16

Keep in mind that the identified local minima may not always be suitable for all interstitials. Users should understand the limitations of this method.

To add interstitial sites #1 and #2, use the add_interstitials_from_local_extrema (abbreviated ai) subcommand:

pydefect_util ai --local_extrema ../unitcell/dos/volumetric_data_local_extrema.json -i 1 2

The selected interstitial sites are added to the supercell_info.json file:

...
-- interstitials
#1
                  Info: all_electron_charge #1
Fractional coordinates: 0.1250000  0.1250000  0.1250000
         Site symmetry: -3m
          Coordination: {'Mg': [1.75, 1.75], 'O': [2.14, 2.14, 2.14, 2.14, 2.14, 2.14]}

#2
                  Info: all_electron_charge #2
Fractional coordinates: 0.5000000  0.5000000  0.5000000
         Site symmetry: -43m
          Coordination: {'Al': [1.75, 1.75, 1.75, 1.75], 'O': [1.57, 1.57, 1.57, 1.57]}

To remove (pop) an interstitial site (e.g., site #1), use:

pydefect pi -i 1 -s supercell_info.json

Then, the first site is removed, and only the second remains:

...
-- interstitials
#1
                  Info: all_electron_charge #2
Fractional coordinates: 0.5000000  0.5000000  0.5000000
         Site symmetry: -43m
          Coordination: {'Al': [1.75, 1.75, 1.75, 1.75], 'O': [1.57, 1.57, 1.57, 1.57]}

To incorporate the added interstitial sites into defect_in.yaml, run the defect_set subcommand again.

7. Creation of defect calculation directories

Next, we create directories for the point-defect calculations using the defect_entries (abbreviated as de) subcommand:

pydefect_vasp de

This command creates directories for each defect, including perfect/.

If you run the same command again, you will see output like the following:

INFO: perfect dir exists, so skipped...
INFO: Al_i1_1 dir exists, so skipped...
INFO: O_i1_0 dir exists, so skipped...
INFO: Mg_i1_1 dir exists, so skipped...
INFO: Mg_i1_2 dir exists, so skipped...
INFO: Al_i1_3 dir exists, so skipped...
INFO: Mg_i1_0 dir exists, so skipped...
INFO: Al_i1_-1 dir exists, so skipped...
INFO: Al_i1_2 dir exists, so skipped...
INFO: Va_O1_1 dir exists, so skipped...
INFO: Va_Al1_-3 dir exists, so skipped...
INFO: Va_O1_0 dir exists, so skipped...
INFO: Va_Al1_0 dir exists, so skipped...
INFO: Va_O1_2 dir exists, so skipped...
INFO: O_i1_-2 dir exists, so skipped...
INFO: O_i1_-1 dir exists, so skipped...
INFO: Va_Mg1_0 dir exists, so skipped...
INFO: Va_Al1_-1 dir exists, so skipped...
INFO: Va_Al1_-2 dir exists, so skipped...
INFO: Va_Mg1_-2 dir exists, so skipped...
INFO: Va_Mg1_-1 dir exists, so skipped...
INFO: Al_i1_0 dir exists, so skipped...
INFO: Va_Al1_1 dir exists, so skipped...

No new directories are created in this case. This fail-safe mechanism prevents accidental overwriting or deletion of existing directories and calculated results.

If you want to recreate the directories, you must remove the existing ones manually.

Each defect directory contains a defect_entry.json file, which holds information about the point defect, prior to running first-principles calculations.

To inspect this file, use the pydefect_print command:

pydefect_print defect_entry.json

Example output:

 -- defect entry info
name: Va_O1_0
site symmetry: .3m
defect center: ( 0.861,  0.861,  0.861)
perturbed sites:
elem dist   initial_coords             perturbed_coords         displacement
  Al 1.92 ( 0.625,  0.875,  0.875) -> ( 0.637,  0.885,  0.864)    0.15
  Al 1.92 ( 0.875,  0.625,  0.875) -> ( 0.879,  0.622,  0.884)    0.09
  Al 1.92 ( 0.875,  0.875,  0.625) -> ( 0.875,  0.875,  0.624)    0.01
  Mg 1.94 ( 1.000,  1.000,  1.000) -> ( 0.004,  0.001,  0.006)    0.06

8. Generation of defect_entry.json

Sometimes, users may need to treat complex defects. For example, O2 molecules can act as anions in MgO2, where vacancies of O2 molecules may be sufficiently stable. Other important cases include methylammonium lead halides (MAPI), where methylammonium ions act as singly charged cations (CH3NH3+), and DX centers, where anion vacancies and interstitial cations coexist.

In such cases, users must prepare the input files and run VASP calculations manually. However, pydefect requires the defect_entry.json file for post-processing, which is not straightforward to create manually.

To assist with this, pydefect provides a subcommand that generates defect_entry.json by analyzing the structural difference between the defect structure and the perfect supercell. The defect charge is determined based on the INCAR, POSCAR, and POTCAR files.

pydefect_vasp_util de -d . -p ../perfect/POSCAR -n complex_defect

This subcommand is also useful when applying pydefect to analyze previously completed defect calculations.

9. Parsing supercell calculation results

Next, run the VASP calculations for each point defect. To generate the VASP input files, use the following command:

for i in */; do cd $i; vise vs -t defect; cd ../; done

Be sure to include the -t defect option to generate input files suitable for defect calculations.

When performing VASP calculations on large supercells with only the Gamma point sampled, we recommend using the Gamma-only version of VASP for efficiency.

After completing the calculations, you can generate the calc_results.json file, which contains the key results from the first-principles calculations relevant to defect properties.

Use the calc_results (abbreviated as cr) subcommand to generate calc_results.json in all defect directories:

pydefect_vasp cr -d *_*/ perfect

If any of the calculations are still in progress, their directories will be automatically skipped during parsing.

10. Corrections of defect formation energies in finite-size supercells

When using supercells under periodic boundary conditions, the total energies of charged defects are not correctly estimated due to interactions between the defect, its periodic images, and the compensating background charge.

Therefore, it is necessary to apply finite-size corrections to the total energies to approximate the dilute limit more accurately.

These corrections can be performed using the extended_fnv_correction (abbreviated as efnv) subcommand:

pydefect efnv -d *_*/ -pcr perfect/calc_results.json -u ../unitcell/unitcell.yaml

This command requires the static dielectric constants and the atomic site potentials of the perfect supercell. Thus, you must provide the paths to both unitcell.yaml and calc_results.json in the perfect/ directory.

Note that the correction process may take some time to complete.

The correction in pydefect is based on the extended Freysoldt–Neugebauer–Van de Walle (eFNV) method.

If you use this correction scheme, please cite the following references:

This command generates a file named correction.pdf, which includes the defect-induced and point-charge potentials as well as their difference at atomic sites. For example:

_images/correction.png

The height of the horizontal line indicates the averaged potential difference between the point-charge potential and the defect-induced potential— specifically, the potential in the defective supercell minus that in the perfect one.

The length of the line indicates the spatial averaging range. For more details, refer to the above publication by Kumagai and Oba (2014).

It is strongly recommended to visually inspect all correction.pdf files for the calculated defects to avoid mistakes in applying corrections.

11. Check defect structures

We analyze the local structures around defects using the defect_structure_info (abbreviated as dsi) subcommand:

pydefect dsi -d *_*/

This generates defect_structure_info.json files for each defect directory. You can view their contents using the pydefect_print command. For example:

 -- defect structure info
Defect type: vacancy
Site symmetry: -3m -> -3m (same)
Has same configuration from initial structure: True
Drift distance: 0.022
Defect center: ( 0.625,  0.375,  0.372)
Removed atoms:
8  Al  ( 0.625,  0.375,  0.375)

Neighbor max distance 2.472
Displacements
Elem  Dist  Displace  Angle  Index  Initial site                  Final site                Neighbor
O     1.9   0.26      160    46     ( 0.639,  0.361,  0.139)  ->  ( 0.647,  0.353,  0.108)  T
O     1.92  0.26      160    42     ( 0.639,  0.139,  0.361)  ->  ( 0.647,  0.110,  0.351)  T
O     1.92  0.26      160    26     ( 0.861,  0.361,  0.361)  ->  ( 0.890,  0.353,  0.351)  T
O     1.92  0.24      170    47     ( 0.389,  0.389,  0.389)  ->  ( 0.360,  0.396,  0.394)  T
O     1.93  0.24      170    31     ( 0.611,  0.611,  0.389)  ->  ( 0.604,  0.639,  0.394)  T
O     1.94  0.23      160    35     ( 0.611,  0.389,  0.611)  ->  ( 0.604,  0.396,  0.637)  T
Al    2.85  0.03      30     14     ( 0.625,  0.125,  0.125)  ->  ( 0.627,  0.128,  0.127)
Al    2.85  0.03      30     21     ( 0.875,  0.375,  0.125)  ->  ( 0.872,  0.373,  0.127)
Al    2.87  0.05      40     16     ( 0.875,  0.125,  0.375)  ->  ( 0.872,  0.128,  0.372)

We can also generate VESTA files to visualize the defect structures using the defect_vesta_file (abbreviated as dvf) subcommand in pydefect_util:

pydefect_util dvf -d *_*/

This creates defect.vesta files in each defect directory.

12. Check defect eigenvalues and band-edge states in supercell calculations

Note: This section is optional.

In general, point defects can be classified into three categories:

(1) Defects with deep localized states inside the band gap. These defects are often detrimental to device performance because they can trap carriers. They may also act as color centers, such as vacancies in NaCl. Thus, it is important to identify the location and character of such defect states.

(2) Defects with hydrogenic carrier states, or perturbed host states (PHS). Here, the carriers are weakly bound to the charged defect center near the band edges. Examples include B-on-Si (p-type) and P-on-Si (n-type) substitutional dopants in silicon. These defects are typically not harmful to devices but may introduce carriers or eliminate counter carriers. The wavefunctions of PHS can extend over millions of atoms, so first-principles calculations of their thermodynamic transition levels require very large supercells, which are computationally infeasible. Thus, their charge transition levels are often not calculated directly but assumed to lie near the band edges based on qualitative assessments.

(3) Defects without any states inside the band gap or near the band edges. These defects have little impact on electronic properties as long as their concentrations remain low.

For examples of each case, see our publications:

To distinguish among these types, one must examine the defect energy levels and wavefunctions to determine whether PHS or localized states are present.

pydefect provides tools to analyze eigenvalues and band-edge states as follows:

First, generate perfect_band_edge_state.json using the perfect_band_edge_state subcommand:

pydefect_vasp pbes -d perfect

Then, generate band_edge_orbital_infos.json files in the defect directories using the band_edge_orbital_infos subcommand (abbreviated as beoi):

pydefect_vasp beoi -d *_* -pbes perfect/perfect_band_edge_state.json

This also creates eigenvalues.pdf plots.

_images/eigenvalues.png

These plots show the single-particle levels and their occupations for spin-up and spin-down channels. The x-axis corresponds to the fractional coordinates of k-points, and the y-axis shows the absolute energy. Each circle represents a state; red, green, and blue denote occupied, partially occupied (0.2–0.8), and unoccupied levels, respectively.

The two horizontal dashed lines indicate the VBM and CBM of the perfect supercell. Band indices are labeled numerically.

Next, generate band_edge_states.json using the band_edge_states (abbreviated as bes) subcommand:

pydefect bes -d *_*/ -pbes perfect/perfect_band_edge_state.json

To examine this file, use:

pydefect_print band_edge_states.json

Example output:

 -- band-edge states info
Spin-up
     Index  Energy  P-ratio  Occupation  OrbDiff  Orbitals                          K-point coords
VBM  128    3.727   0.33     1.00        0.02     O-p: 0.76                         ( 0.250,  0.250,  0.250)
CBM  129    9.980   0.06     0.00        0.03     Al-s: 0.11, O-s: 0.21, O-p: 0.11  ( 0.250,  0.250,  0.250)
vbm has acceptor phs: False (0.000 vs. 0.2)
cbm has donor phs: False (0.000 vs. 0.2)
---
Localized Orbital(s)
Index  Energy  P-ratio  Occupation  Orbitals

Spin-down
     Index  Energy  P-ratio  Occupation  OrbDiff  Orbitals              K-point coords
VBM  126    4.083   0.67     0.88        0.06     O-p: 0.72             ( 0.250,  0.250,  0.250)
CBM  129    10.010  0.06     0.00        0.01     O-s: 0.21, O-p: 0.11  ( 0.250,  0.250,  0.250)
vbm has acceptor phs: False (0.120 vs. 0.2)
cbm has donor phs: False (0.000 vs. 0.2)
---
Localized Orbital(s)
Index  Energy  P-ratio  Occupation  Orbitals
127    4.277   0.62     0.06        O-p: 2.98
128    4.278   0.62     0.06        O-p: 2.98

The analysis is given for both spin channels.

P-ratio is the participation ratio, which measures how much of the orbital density is localized on neighboring atoms (as listed in defect_structure_info.json) relative to the total over the entire cell.

As mentioned, the formation energies of defects with occupied donor phs or unoccupied acceptor phs should typically be excluded from formation energy diagrams.

In pydefect, these classifications are based on the eigenvalues and the similarity of the wavefunctions to those of the VBM and CBM. For details, please refer to our upcoming publication.

Note: Automatic classification of PHS and localized states is not always reliable. Users are encouraged to visually inspect the band-edge states and examine band-decomposed charge densities when necessary.

13. Plot defect formation energies

Here, we demonstrate how to plot defect formation energies.

Plotting defect formation energies requires multiple types of information: the band edges, chemical potentials of the competing phases, and the total energies of both perfect and defective supercells.

First, use the defect_energy_infos (= dei) sub-command:

pydefect dei -d *_*/ -pcr perfect/calc_results.json -u ../unitcell/unitcell.yaml -s ../cpd/standard_energies.yaml

This command generates the defect_energy_info.yaml files within each defect directory. An example is shown below:

name: Va_O1
charge: 0
formation_energy: 6.803585744999999
atom_io:
  O: -1
energy_corrections:
  pc term: 0.0
  alignment term: -0.0
is_shallow: False

Caveats: (1) The formation_energy is calculated assuming the chemical potentials of the elements are at their standard states, and the Fermi level is at 0 eV. (2) Although two energy corrections are included (pc term and alignment term), you may add custom correction terms if needed. (3) The is_shallow value will be empty if the previous section (on PHS identification) is skipped. Here, “shallow” refers to the presence of a perturbed host state (PHS). You may modify the classification (e.g., shallow/deep) manually if desired.

Next, create the defect_energy_summary.json file using the defect_energy_summary (= des) sub-command:

pydefect des -d *_*/ -u ../unitcell/unitcell.yaml -pbes perfect/perfect_band_edge_state.json -t ../cpd/target_vertices.yaml

This sub-command collects information from each defect_energy_info.yaml file and compiles it into defect_energy_summary.json.

title: MgAl₂O₄
rel_chem_pots:
 -A Al: 0.00 Mg: -0.69 O: -5.25
 -B Al: 0.00 Mg: -0.32 O: -5.34
 -C Al: -7.87 Mg: -5.94 O: 0.00
 -D Al: -8.01 Mg: -5.66 O: 0.00
vbm: 0.00, cbm: 5.22, supercell vbm: -0.03, supercell_cbm: 6.14

name    atom_io         charge    energy    correction  is_shallow
------  ------------  --------  --------  ------------  ------------
Al_i1   Al: 1               -1    14.958        -0.049  False
                             0     9.896         0.000  True
                             1     3.376         0.527  False
                             2    -0.872         1.814  True
                             3    -6.994         3.370  False
Mg_i1   Mg: 1                0     8.208         0.000  True
                             1     2.015         0.555  True
                             2    -4.117         1.534  False
O_i1    O: 1                -2     8.122         0.721  False
                            -1     7.055         0.149  False
                             0     4.566         0.000  False
Va_Al1  Al: -1              -3    13.527         3.402  False
                            -2    12.968         1.789  True
                            -1    12.672         0.660  True
                             0    12.595         0.000  True
                             1    12.710        -0.232  True
Va_Mg1  Mg: -1              -2    10.093         1.524  False
                            -1     9.580         0.578  True
                             0     9.311         0.000  False
Va_MgO  O: -1 Mg: -1         0     9.887         0.000  False
Va_O1   O: -1                0     6.804         0.000  False
                             1     3.808         0.166  False
                             2     0.883         0.845  False

You can also create a calc_summary.json file using the calc_summary (= cs) sub-command:

pydefect cs -d *_*/ -pcr perfect/calc_results.json

This outputs:

|:---------:|:----------:|:-----------:|:-----------------:|:------------:|:------------------:|:--------------:|
|   name    | Ele. conv. | Ionic conv. | Is energy strange | Same config. |    Defect type     | Symm. Relation |
| Al_i1_-1  |     .      |      .      |         .         |    False     |      unknown       |    subgroup    |
|  Al_i1_0  |     .      |      .      |         .         |    False     |      unknown       |    subgroup    |
|  Al_i1_1  |     .      |      .      |         .         |    False     |      unknown       |    subgroup    |
|  Al_i1_2  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Al_i1_3  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Mg_i1_0  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Mg_i1_1  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Mg_i1_2  |     .      |      .      |         .         |      .       |         .          |       .        |
|  O_i1_-1  |     .      |      .      |         .         |      .       |         .          |       .        |
|  O_i1_-2  |     .      |      .      |         .         |    False     | interstitial_split |    subgroup    |
|  O_i1_0   |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Al1_-1 |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Al1_-2 |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Al1_-3 |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Al1_0  |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Al1_1  |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Mg1_-1 |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Mg1_-2 |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_Mg1_0  |     .      |      .      |         .         |      .       |         .          |       .        |
| Va_MgO_0  |     .      |      .      |         .         |    False     |      unknown       |      same      |
|  Va_O1_0  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Va_O1_1  |     .      |      .      |         .         |      .       |         .          |       .        |
|  Va_O1_2  |     .      |      .      |         .         |      .       |         .          |       .        |

Note: this feature is still in beta.

Finally, plot the defect formation energies as a function of the Fermi level using the plot_defect_formation_energy (= pe) sub-command:

pydefect pe -d defect_energy_summary.json -l A

This generates a plot like the following:

_images/energy_A.pdf

To change the chemical potential condition—i.e., select a different vertex from the chemical potential diagram—use the -l option.