Finite-size corrections to charged defect supercells with FHI-aims¶

1. Introduction¶

This jupyter notebook workflow has been designed for applying finite-size correction schemes to charged defects when calculating the formation energy for a charged point defect with the supercell method and outputs from the all-electron electronic structure software package FHI-aims. This notebook provides a step-by-step explanation of the processing of the defect data when applying the correction schemes utilised in this workflow as the notebook progresses.

This notebook can be used to perform corrections for one defect at a time. However, if the user has a large set of defect calculations to process, another script has been made available called DefectCorrectionsDataset.py. This uses the NotebookScripter library to run this notebook over a set of defect calculations from the terminal. It is still advisable to first step through the workflow in this notebook to be aware of the steps and assumptions made during the processing of the data. The entire notebook can also be run by selecting 'Kernel --> Restart & Run All'. Outputs generated by this workflow are stored in a directory 'ProcessedDefects' in the same directory as this notebook.

Background reading¶

For background reading on electronic structure calculations for defects in solids, doi: 10.1103/RevModPhys.86.253 provides a thorough review on the topic where it is the 'supercell approach' that is utilised in this workflow. Other important background reading includes doi: 10.1088/0965-0393/17/8/084002, doi: 10.1103/PhysRevLett.102.016402, doi: 10.1016/j.cpc.2018.01.011 and doi: 10.1103/PhysRevB.89.195205. Methods described in all of these works are relevant for this workflow.

Calculating the formation energy of charged defect supercells¶

The formation energy of a charged defect is given by:

$\Delta H_{D,q} = [E_{D,q} + E^{corr}_{q}] - E_{host} + \sum_i n_i (E_i + \mu_i) + q[\epsilon_F + \epsilon_{\nu} + \Delta \nu_{0/b}]$

Where $E_{D,q}$ is the total energy for a defect in charge state q embedded in a supercell from electronic structure calculations, $E^{corr}_{q}$ are finite-size corrections to this defect supercell, $E_{host}$ is the calculated total energy of an equivalent perfect host lattice supercell, $n_i$ is the number of atoms of species i added to ($n_i > 0$) or removed from ($n_i < 0$) the chemical reservoir when the defect is formed. $\mu_i$ is the chemical potential of species i, referenced to the total energy of the pure element in its standard state, $E_i$, $\epsilon_F$ is the position of the Fermi level in the band gap, $q \epsilon_{\nu}$ (energy of bulk VBM) and $\Delta \nu_{0/b}$ (term used to align the electrostatic potential of the VBM for the bulk and defect supercells). Please refer to papers listed in the 'background reading' section for more details.

Finite-size corrections schemes¶

In this workflow, it is possible to obtain the $E^{corr}_{q}$ term for a charged defect in periodic DFT calculations with dielectric screening from the host lattice with either the LZ (doi: 10.1088/0965-0393/17/8/084002) or FNV (doi: 10.1103/PhysRevLett.102.016402) finite-size correction schemes. The FNV scheme is implemented here by using the CoFFEE python code (doi: 10.1016/j.cpc.2018.01.011) and in this workflow functions from pylada-defects (doi: 10.1016/j.commatsci.2016.12.040) are utilised to obtain the $E^{corr}_{q}$ term within the LZ scheme. The option to implement either or both of these schemes to obtain values for $E^{corr}_{q}$ is determined in the user inputs cell by setting 'LZ' or 'FNV' to 'True'. Note that the FNV scheme involves an alignment-like step and this can be performed using either the planar average of the electrostatic potential or averaged atom-site potentials with the sampling region proposed by Kumagai and Oba (doi: 10.1103/PhysRevB.89.195205), see Fig. 2(a) in the publication. The option to implement either or both of these methods is determined in the user inputs cell by setting 'atom_centered_pa' or 'planar_av_pa' to 'True'.

Currently, band filling corrections to $E^{corr}_{q}$ (descriptions of which can be found in in doi: 10.1103/physrevb.72.035211, 10.1103/PhysRevB.78.235104 and 10.1088/0965-0393/17/8/084002) that are most relevant for shallow defects with delocalised defect-induced charge, are not currently available in this workflow (but would make for a nice extension!).

The potential alignment step (towards the end of the workflow) to obtain $\Delta \nu_{0/b}$ can also use either planar averages of the electrostatic potential or averaged atom-site potentials with the sampling region proposed by Kumagai and Oba. The option to implement either or both of these methods is determined in the user inputs cell, also by setting 'atom_centered_pa' or 'planar_av_pa' to 'True'.

Note: For potential alignment, it can be difficult to generate clear plots with planar averages of potentials if the defect supercells have been relaxed. Using averaged atom-site potentials by setting 'atom_centered_pa' to 'True' is usually more robust for this.

2. Assumptions¶

Necessary data¶

1. Components of the static dielectric tensor of the host crystal (calculated electronic + ionic contributions). This must be added to the user inputs cell. It is best to use computed values for dielectric constants for consistent corrections.

2. Relaxed structures for charged defect supercells, neutral defect supercells and equivalent perfect supercell. (Supercells have all been relaxed and then a single one-shot is performed where control.in tags for geometry relaxation have been removed or commented out and the tag for atom site potentials and the planar average has been included output on-site esp and output realspace_esp 100 100 100).

Defect supercells¶

1. The VOLUME of the supercells has NOT been relaxed so lattice vectors are the same for the perfect and defect supercells (this is appropriate for simulating defects in the dilute limit).

2. DFT data being processed is for INTRINSIC defects only (there are no atom types present in the defect supercell that are not also present in the perfect supercell)... but modifying the workflow to allow for extrinsic defects would be a nice extension!

3. A maximum of one point defect is present in each supercell.

4. When vacancies are created in the supercell, the line for the vacancy site has been removed (not just commented out).

5. The supercells have orthogonal axes for the lattice vectors (a small amount of numerical noise in the off-diagonals is okay).

6. Geometry files can be either in fractional or Cartesian coordinates as they are automatically all converted into Cartesian coordinates in the workflow.

The above assumptions are necessary for the methods used to locate the defects in this workflow using functions in DefectSupercellAnalyses.py.

Format of data¶

1. All output filenames from FHI-aims calculations still have their original default names: 'geometry.in', 'On-site_ESP.dat' and 'plane_average_realspace_ESP.out'.

2. The data from all final one-shot calculations of the perfect and defect supercells are stored in one directory. Sub-directories for the different defects within this directory are fine (will be searched through by this notebook) and are to be expected to keep data organised!

3. As no relaxation steps were performed on the final structures, the only geometry files in these directories should be 'geometry.in' files (i.e. no 'geometry.in.next_step' files) and these correspond to the final relaxed structure for each supercell.

Image charge correction for the LZ scheme¶

The calculation of $E^{corr, LZ}_q$ in this workflow makes the approximation of isotropic dielectric screening from the host lattice (from the use of the isotropic average of the static dielectric constant). Extending this workflow to allow for anisotropic dielectric screening with all components of the static dielectric tensor of the host lattice is another good option for future development.

Charge model for the FNV scheme¶

See CoFFEE paper (doi: 10.1016/j.cpc.2018.01.011) and later cells in this notebook for more details on the role of the Gaussian charge model used in this correction scheme. Assumptions made in this workflow:

1. The same sigma value is used for the Gaussian charge model implemented in the CoFFEE package for all defects in the same set for consistency. The maximum sigma is assumed to be the minimum distance of all defects in the set to the supercell boundary, but note that in reality the charge associated with some defects may be more or less disperse than others.

2. This workflow sets the cutoff for the plane waves used in CoFFEE's Poisson solver. Smaller sigma values require more plane waves to achieve convergence and the more plane waves the slower the solver will be. If it turns out that you are unable to achieve convergence with the cutoff value set in the notebook based on sigma, there is also the option to override the value and set this manually. This is the 'manual_cutoff' variable. Tests so far have not shown this to be necessary.

3. This workflow does not currently include band filling corrections. It is assumed that the charge associated with a defect is contained within the supercell.

3. User inputs¶

The cell below is the only cell that must be modified by the user before running the notebook. All subsequent cells do not require user inputs. They just need to be run in order and there is an explanation of the analysis being performed before each cell. The only exception is if manual overrides of certain parameters turn out to be necessary for certain defects, but instructions for how to implement this will be contained in the 'log.info' file if necessary. Otherwise, the user may also want to modify the plots generated by the workflow and instructions on how to do this are contained in the cell with the code to generate each plot. Note that some cells will take longer to run than others.

To run this analysis for a full set of defects data, this information can be inputted as a list in the DefectCorrectionsDataset.py script and the cell below is replaced by an argument list when processing a set of defects with DefectCorrectionsDataset.py (using the NotebookScripter library). Additional instructions are contained in DefectCorrectionsDataset.py.

4. Python libraries used in the workflow¶

5. Defining paths to input and output files and set up log file¶

6. Reading information for the perfect host crystal supercell¶

It has been assumed that the volume of defect supercells has not been relaxed, therefore the lattice vectors of the perfect host supercell are the same as those of the defect supercells. (This was listed in the assumptions list earlier in the notebook)

7. Locating the defect position in the supercell¶

This script locates the coordinates of the defect by comparing the defect supercell to the perfect host supercell and using the functions defined in DefectSupercellAnalyses.py. The defect position is referenced to the species that would be in the perfect host in the case of a vacancy, e.g. the S atom removed from the host to make a S vacancy. The position is referenced to the atom added to the defect supercell in the case of an interstitial or antisite.

These coordinates will be used to define the centre of the Gaussian used in the charge model for the FNV method in the input 'in' file used with the Poisson solver in coffee.py. The z-position of the defect will also be used when generating plots for the potential alignment (pa) correction step towards the end of this workflow.

7.1. Finding smallest distances of the defect from the supercell boundary¶

In the script below the minimum distance of the defect to the supercell boundary is computed. In the plot of the defect supercell, the volume of the region used to represent the defect location is arbitrary and not representative of atomic radii.

8. Lany-Zunger (LZ) finite-size correction scheme¶

For more information on the LZ scheme, please see doi: 10.1088/0965-0393/17/8/084002. This workflow uses modified functions from pylada-defects. For more information on this implementation of this correction scheme, please see doi: 10.1016/j.commatsci.2016.12.040.

As a reminder, the formation energy of a charged defect is given by:

$\Delta H_{D,q} = [E_{D,q} + E^{corr}_{q}] - E_{host} + \sum_i n_i (E_i + \mu_i) + q[\epsilon_F + \epsilon_{\nu} + \Delta \nu_{0/b}]$

For this correction scheme, we will first obtain $E^{corr}_{q}$ as the 'image charge correction' in the LZ scheme (but note that this workflow does not include band filling corrections as mentioned already in the assumptions cell in this notebook). After that, we perform potential alignment to obtain the $\nu_{0/b}$ term.

8.1 Obtain $E^{corr}_{q}$ from the LZ scheme¶

The $E^{corr}_q$ term for the defect formation energy is first obtained in this workflow using the LZ correction scheme (see doi: 10.1088/0965-0393/17/8/084002 and doi: 10.1016/j.commatsci.2016.12.040 for more information). Here it is just the 'image-charge correction' that is being computed first and is denoted as $E^{corr}_q$, while sometimes in the literature $E^{corr}_q$ is used to refer to this correction, band filling corrections and the potential alignment that is performed towards the end of this workflow. To avoid confusion with the term calculated next with the FNV scheme, $E^{corr, LZ}_q$ is how this term will be denoted in this workflow. Within the LZ scheme, this correction term is given by:

$E^{corr, LZ}_q = [ 1 + c_{sh}(1- \frac{1}{\epsilon}) ] \frac{q^2 \alpha_M}{2 \epsilon L}$

Where $L = \Omega^{\frac{-1}{3}}$ is the linear supercell dimension, which has volume $\Omega$, $\epsilon$ is the isotropic average of the calculated static dielectric constant inputed in the 'user inputs' cell, $\alpha_M$ is the Madelung constant and $c_{sh}$ is the shape factor. The latter two terms are related to the supercell geometry. The approximation of isotropic dielectric screening from the host lattice has been made here, this was listed in the assumptions cell earlier in this notebook.

Note: The cell below may take a while to run.

8.2. Potential alignment to obtain $\Delta \nu_{0/b}$ with the LZ scheme¶

As a reminder, the formation energy of a charged defect is given by:

$\Delta H_{D,q} = [E_{D,q} + E^{corr}_{q}] - E_{host} + \sum_i n_i (E_i + \mu_i) + q[\epsilon_F + \epsilon_{\nu} + \Delta \nu_{0/b}]$

In the cell above, $E^{corr}_{q}$ was calculated within the LZ finite-size correction schemes. This part of the workflow obtains the $\Delta \nu_{0/b}$ potential alignment term for the expression above for $\Delta H_{D,q}$ by comparing the electrostatic potentials from electronic structure calculations of a perfect host supercell to that of a charge neutral defect, in a region far from the defect.

Here, this term can be obtained by either using planar averages of the electrostatic potential or averaged atom-site potentials with the sampling region proposed by Kumagai and Oba, depending on settings specified in the user inputs cell.

Note: The atom corresponding to the defect is omitted in the following potential alignment plots.

8.2.1. Using atom centres and KO sampling region¶

The script below is used to generate a potential alignment plot between the perfect host supercell and the charge neutral defect supercell using potentials from the electronic structure calculations and the sampling region proposed by Kumagai and Oba for performing the potential alignment step with averaged atom centre potentials. For more information see doi: 10.1103/PhysRevB.89.195205) and Fig. 2(a) in the paper.

8.2.2. Using planar averages of calculated potentials¶

Here the alignment term for the LZ scheme is obtained by aligning planar averages of the potentials from electronic structure calculations of the perfect host supercell and charge neutral defect supercell in a region far from the defect.

Note: It is usually better to use atom potentials to align the potential if the defect supercell has been relaxed.

9. The Freysoldt-Neugebauer-Van de Walle (FNV) finite-size correction scheme¶

For further reading on this finite size correction scheme see doi: 10.1103/PhysRevLett.102.016402. The FNV scheme is implemented here by using the CoFFEE python code. For further information on this specific implementation see doi: 10.1016/j.cpc.2018.01.011.

As a reminder, the formation energy of a charged defect is given by:

$\Delta H_{D,q} = [E_{D,q} + E^{corr}_{q}] - E_{host} + \sum_i n_i (E_i + \mu_i) + q[\epsilon_F + \epsilon_{\nu} + \Delta \nu_{0/b}]$

The finite-size electrostatic correction term from the FNV correction scheme is given by: $E^{corr}_{q} = E_{q}^{lat} - q\Delta V_{q-0/m}$

$E_{q}^{lat}$ is obtained from a Gaussian charge model and the next term is determined later in this workflow through a 'potential alignment like' procedure which compares the model potential to the difference in potentials for the defect in a charged and neutral states from the electronic structure calculations. This will be discussed more in later cells when this procedure is applied.

Charge model for the FNV scheme¶

This workflow uses the CoFFEE software package when implementing the FNV correction scheme (doi: 10.1103/PhysRevLett.102.016402) for finite-size corrections to charged defect supercells. This correction scheme involves solving the Poisson equation for a charge model system (to infer the electrostatic interactions between the defect and its periodic images) and aligning potentials to derive the correction to the formation energy. The model system used to simulate defect charge interactions with periodic images is a Gaussian within a periodicly repeated supercell where the calculated static dielectric tensor of the host lattice is used to represent the dielectric screening of the defect charge by the host lattice. The schematic below shows a charge model used to simulate the interaction between the defect charge in the central supercell and some of its periodic images. For more details on the charge model used when performing corrections to charged defect supercells with the scheme used in this workflow, please refer to the CoFFEE paper (doi: 10.1016/j.cpc.2018.01.011).

Determine sigma for the Gaussian charge models for the defect set¶

The purpose here is to determine the value of sigma for the Gaussian charge model that is small enough to ensure that the charge is contained within the supercell for all supercells in the defect set, including the ones where the defect is closest to the boundary of the supercell. For consistency, the same value of sigma is used for all defects in the set of calculations. As stated in the CoFFEE paper (doi: 10.1016/j.cpc.2018.01.011): 'For bulk systems, it is not necessary that the width of the Gaussian model charge match the defect wavefunction charge density. It suffices if the width is appropriately small to keep the model charge inside the cell.' As already mentioned, band-filling corrections are not currently applied in this workflow. All of these assumptions were listed earlier in the notebook in the assumptions cell.

This step needs to only be run once for your set of defects, so if it turns out to be slow for a large set of defects just run it once and for subsequent runs replace the cell below with a cell that just contains:

sigma = computed_value_from_first_run

However, tests so far have not shown this to be necessary.

User option: manually override sigma value¶

Determine cutoff for the Poisson solver based on sigma of Gaussian charge model¶

The script below sets the cutoff for the plane waves used in CoFFEE's Poisson solver. Smaller sigma values require more plane waves to achieve convergence and the more plane waves the slower the solver will be.

If it turns out that you are unable to achieve convergence with the cutoff value set based on sigma, there is also the option to override the value and set this manually.

User option: manually override cutoff value¶

9.1. Determine $E_{q}^{lat}$ from the Gaussian charge model¶

The CoFFEE software package has been developed to be compatible with 1D, 2D and 3D systems. The long range model potential for bulk (3D) systems can be analytically derived, but in the case of spatially varying dielectric profiles (in 1D and 2D systems) it is practical to solve the Poisson equation numerically and extrapolate the electrostatic energy to the isolated limit $E_q^{iso,m}$. Although this workflow is for 3D defect supercells, it follows the general formalism in the CoFFEE workflow that is applicable to 1D, 2D and 3D systems. Therefore, the electrostatic energy is determined for a supercell that is the same size as that used in the electronic structure calculations with FHI-aims (this is denoted by $E_q^{per,m}$) and also for 2x2x2 and 3x3x3 supercells ($E_q^{per,2m}$ and $E_q^{per,3m}$ respectively) of that supercell to allow $E_q^{iso,m}$ to be determined by extrapolating the fit to the energies for the three supercell sizes. $E_q^{lat}$ is then obtained from the equation below.

$E_q^{lat} = E_q^{iso,m} - E_q^{per,m}$

See the CoFFEE paper for more details.

Run coffee.py Possion solver for charge models¶

This is used to obtain $E_q^{per,m}$ for the charge model that is the same size as the supercell used in the DFT defect calculation and also to calculate $E_q^{per,2m}$ and $E_q^{per,3m}$ for 2x2x2 and 3x3x3 supercells of the DFT supercell.

The values computed for all three supercells will be used to obtain $E_q^{iso,m}$ for the charge model from a fit to the three data points and extrapolating to the limit of an infinite supercell.

Note: The cell below may take a while to run.

Obtain $E_q^{iso}$ from fit to $E_q^{per,m}$ for defect supercell, 2x2x2 and 3x3x3 of this supercell¶

Below is a plotting script adapted from plot_fit.py from the CoFFEE_1.1 'Examples' directory. It uses a fit to $E_q^{per,m}$, $E_q^{per,2m}$, $E_q^{per,3m}$ to extrapolate to the infinite supercell limit to get obtain $E_q^{iso,m}$.

$E_q^{iso,m}$ from this fit is used after to compute $E_q^{lat}$.

Obtain $E_q^{lat}$ term from calculated values¶

$E_q^{lat} = E_q^{iso,m} - E_q^{per,m}$

9.2. Align potentials of charge model and defect supercell to obtain $q\Delta V_{q-0/m}$¶

As a reminder, the finite-size electrostatic correction term from the FNV correction scheme is given by: $E_q^{corr} = E_q^{lat} - q\Delta V_{q-0/m}$

This step is used to calculate $q\Delta V_{q-0/m}$ to obtain the full expression for $E_q^{corr}$ within the FNV finite-size correction scheme (not accounting for any band filling effects for shallow defects). A 'potential alignment like' procedure is used to obtain $q\Delta V_{q-0/m}$ where the potential of the charge model is compared to the difference in potentials for the defect in a charged and neutral states from the electronic structure calculations.

As a reminder, the formation energy of a charged defect is given by:

$\Delta H_{D,q} = [E_{D,q} + E^{corr}_{q}] - E_{host} + \sum_i n_i (E_i + \mu_i) + q[\epsilon_F + \epsilon_{\nu} + \Delta \nu_{0/b}]$

When implementing the FNV scheme, it is possible to obtain $E^{corr}_{q}$ and $\Delta \nu_{0/b}$ at the same time by also aligning the potentials of the perfect host supercell to the charge neutral defect in the region far from the defect. All terms will be generated in the plots below.

Note: The atom corresponding to the defect is omitted in the following potential alignment plots.

9.2.1. Using atom centres and KO sampling region¶

The script below is used to generate the sampling region proposed by Kumagai and Oba for performing the potential alignment step with averaged atom centre potentials. For more information see doi: 10.1103/PhysRevB.89.195205) and Fig. 2(a) in the paper.

First calculate atom centre potentials from the CoFFEE charge model:¶

Generate plot to align atom centre potentials of charge model and defect supercells:¶

9.2.1. Using planar average of electrostatic potential¶

This is the default method for performing the potential aligment step in the CoFFEE implementation of the FNV scheme.

Note: It is usually better to use atom potentials to align the potential if the defect supercell has been relaxed.

Obtain planar average of the potential of the charge model:¶

This is needed to align potentials between the charge model and the defect supercells when using planar averages for the alignment method. This is used with planar averaged potentials along the z-direction (referred to as 'a3' in CoFFEE) from FHI-aims calculations for the perfect host, charged defect and neutral defect supercells.

First we write an 'in_V' file for the supercell to use with the utility script from CoFFEE_1.1/PotentialAlignment/Utilities/plavg.py

The script below uses V_r.npy generated from the coffee.py Poisson solver for the charge model with the same supercell size as the DFT calculations with the in_V file written in cell above to obtain planar average of potential for the charge model along the z-direction.

This produces to following file: plavg_a3.plot, which is in units of Bohr.

Generate plot to align planar average of potentials of charge model and defect supercells:¶

The final plotting script below uses plavg_a3.plot for the planar average of the charge model from CoFFEE with 'plane_average_realspace_ESP.out' files from FHI-aims calculations for the perfect host, neutral defect and charged defect supercells.

TZ: Potential Alignment from planar avearge of potentials

CoFFEE format: divided into two plots: 
            q((V_Charge-V_0) - V_model) (FNV_planar_pa_C2.png)     q(V_0-V_host) (FNV_planar_pa_C1.png)

In the Original Freysoldt Correction PRL paper: 
        Just   q((V_charge-V_host) - V_model)  (FNV_planar_pa_Frey.png)

TZ: Several Enhancement. We can either contact Mit Naik or change that by ourself. 1. The plane average electrostatic potential of the model from CoFFEE seems not correct, the center are not exact in the defect position, a little shift.
2. The grids of plane average electrostatic potential of the model from the CoFFEE are not the same as that in the FHI-aims. It's easy to let the user input the grids they want to produce the electrostatic potential, so that it can produce the Model potential at the exact same grids as that in the FHI-aims. Based on these Model potential, we can plot the Difference easily (V(Defect) -V(Host)) - V(Model). Otherwise, it's not easy to plot this Diff. directly in the plot.

10. Write file containing computed values for corrections¶

11. Cleanup unnecessary outputs¶

The cell below removes unnecessary output files. See 'log.info' for an overview of the outputs from this notebook and 'corrections_summary.dat' for the final correction terms that have been computed. Have a nice day!