CoFFEE: Corrections For Formation Energy and Eigenvalues for charged defect simulations

TL;DR

The paper introduces CoFFEE, a Python-based tool implementing FNV correction scheme to accurately compute formation energies and eigenvalues of charged defects in various material geometries, improving defect simulation accuracy.

Contribution

It presents CoFFEE, a versatile and optimized Python code that applies FNV corrections for charged defect simulations across different material shapes.

Findings

CoFFEE effectively corrects formation energies and eigenvalues.

The code supports diverse geometries including bulk, slab, wires, and nanoribbons.

Parallelization and optimization enhance computational efficiency.

Abstract

Charged point defects in materials are widely studied using Density Functional Theory (DFT) packages with periodic boundary conditions. The formation energy and defect level computed from these simulations need to be corrected to remove the contributions from the spurious long-range interaction between the defect and its periodic images. To this effect, the CoFFEE code implements the Freysoldt-Neugebauer-Van de Walle (FNV) correction scheme. The corrections can be applied to charged defects in a complete range of material shapes and size: bulk, slab (or two-dimensional), wires and nanoribbons. The code is written in Python and features MPI parallelization and optimizations using the Cython package for slow steps.

Equations

$${\mathrm{E}}_0^f[{\vec{\mathrm{R}}}_0]({ϵ}_F)={\mathrm{E}}_0^{\mathrm{tot}}[{\vec{\mathrm{R}}}_0]-{\mathrm{E}}_{\mathrm{pristine}}-n_x{\mu }_x,$$

$${\mathrm{E}}_q^f[{\vec{\mathrm{R}}}_q]({ϵ}_F)=\{{\mathrm{E}}_q^{\mathrm{tot}}[{\vec{\mathrm{R}}}_q]+{\mathrm{E}}_q^{\mathrm{corr}}\}-{\mathrm{E}}_{\mathrm{pristine}}+q\{{ϵ}_{\mathrm{vbm}}^{\mathrm{pristine}}+{ϵ}_F-\Delta V_{0/p}\}-n_x{\mu }_x$$

$$\Delta V_{0/p}=V_0{|}_{\mathrm{far}}-V_p$$

$${\mathrm{E}}_q^{\mathrm{corr}}={\mathrm{E}}_q^{\mathrm{lat}}-q\Delta V_{q-0/m}$$

$${\mathrm{E}}_q^{\mathrm{per},\mathrm{m}}=\frac{1}{2}{\int }_{\Omega }{\rho }^m(\mathbf{r})V_q^{\mathrm{per},\mathrm{m}}(\mathbf{r})d\mathbf{r}$$

$$\Delta V_{q-0/m}=(V_q^{\mathrm{DFT}}-V_0^{\mathrm{DFT}}){|}_{\mathrm{far}}-V_q^{\mathrm{per},\mathrm{m}}{|}_{\mathrm{far}}$$

$$\nabla .[\epsilon (\mathbf{r})\nabla V_q^{\mathrm{per},\mathrm{m}}(\mathbf{r})]=-4\pi {\rho }^m(\mathbf{r})$$

$$\sum _{{\mathbf{G}}^{\prime }}\sum _{i=1}^3{G}_i{G}_i^{\prime }{\epsilon }_{ii}(\mathbf{G}-{\mathbf{G}}^{\prime })V_q^{\mathrm{per},\mathrm{m}}({\mathbf{G}}^{\prime })=4\pi {\rho }^m(\mathbf{G})$$

$$V_q^{\mathrm{per},\mathrm{m}}(\mathbf{G})=\frac{4\pi {\rho }^m(\mathbf{G})}{\epsilon |\mathbf{G}{|}^2}$$

$$\epsilon (z)=\frac{1}{2}({\epsilon }^v-{\epsilon }^m)\mathrm{erf}(\frac{z-t_1}{s})-\frac{1}{2}({\epsilon }^v-{\epsilon }^m)\mathrm{erf}(\frac{z-t_2}{s})+{\epsilon }^v$$

$$\begin{gathered} \sum _{G_z^{\prime }}[{\epsilon }_{\perp }(G_z-G_z^{\prime })]G_z^{\prime }{G}_z{V}_q^{\mathrm{per},\mathrm{m}}(G_x,G_y,G_z^{\prime }) \\ +\sum _{G_z^{\prime }}[{\epsilon }_{\parallel }(G_z-G_z^{\prime })](G_x^2+G_y^2)V_q^{\mathrm{per},\mathrm{m}}(G_x,G_y,G_z^{\prime })=4\pi {\rho }^m(G_x,G_y,G_z) \end{gathered}$$

$$\begin{gathered} \sum _{G_x^{\prime },G_y^{\prime }}[{ϵ}_{\parallel }(G_x-G_x^{\prime },G_y-G_y^{\prime })G_z^2 \\ +{ϵ}_{\perp }(G_x-G_x^{\prime },G_y-G_y^{\prime })(G_x{G}_x^{\prime }+G_y{G}_y^{\prime })]V(G_x^{\prime },G_y^{\prime },G_z) \\ =4\pi \rho (G_x,G_y,G_z) \end{gathered}$$

$${ϵ}_q^{\mathrm{corr}}=\frac{-2}{q}{\mathrm{E}}_q^{\mathrm{corr}}$$

$$p(\Omega )=f_1+f_2/{\Omega }^{1/3}+f_3/\Omega$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: CoFFEE

Program obtainable from: http://www.physics.iisc.ernet.in/~mjain/pages/software.html

Journal Reference:

Catalogue identifier:

Licensing provisions: Open source BSD License

Programming language: Python

Computer: Any computer with Python installed. The code has been tested with Python2.7 and Python3.6.

Operating system: Unix/Linux/Windows

RAM: 5-1000 MB (dependent on system size)

Keywords: Density Functional Theory, Defect formation energy, Charged defects, Eigenvalue corrections, Electronic Structure, GW

Classification: 7.1, 7.3

External routines/libraries: numpy, scipy, mpi4py, matplotlib

Nature of problem: Most electronic structure codes based on Density Functional Theory use periodic boundary conditions. This leads to spurious electrostatic interactions during simulation of charged defects, which affects the computed defect formation energy and the defect eigenvalue. Solution method: We implement the Freysoldt-Neugebauer-Van de Walle (FNV) correction scheme to correct the defect formation energy and eigenvalues. Our implementation can be applied to charged defects in 3D bulk materials as well as materials having 2D and 1D geometries.

Running time: 1-600 minutes (depends on the number of processors and system size)

1 Introduction

Point defects, ubiquitous in materials, influence their electrical and optical properties. First principles electronic structure calculations have proven to be vital in understanding and predicting the role of defects [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. Defect engineering to enhance or suppress certain characteristics of materials often rely on such simulations for inputs [19, 20, 21, 22, 23, 24, 25].

Intrinsic point defects in materials are generally found in low concentrations, about one in a million atoms in 3D solids to one in a thousand atoms in 2D materials. Simulations thus attempt to study isolated defects in materials. However, a large number of first principles calculation codes based on Density Functional Theory (DFT) employ periodic boundary conditions. As a result, in order to simulate and understand isolated defects, super cells are constructed. The size of the super cell is chosen to minimize the overlap of the defect wavefunction with its periodic image in neighboring cells [1, 26, 27, 28, 29, 30, 31]. The computed formation energy of non-shallow neutral defects in such calculations is found to converge quickly with the super cell size.

Charged point defects are simulated by artificially introducing a compensating uniform background charge to avoid the divergence in the electrostatic energy. The formation energy of charged defects show slow convergence with super cell size due to the Coulomb interaction between the defect charge and its periodic images. The super cell sizes necessary to completely overcome these spurious defect-defect interactions are computationally intractable. The defect energy levels in the gap are also affected similarly, and shows slow convergence with super cell size [27, 32, 33]. Furthermore, since the uniform background charge is artificially introduced, an absolute reference for the electrostratic potential is undefined. The potential in the defect supercell thus needs to be aligned with the bulk potential, in accordance with the chosen reference of the Fermi level with respect to the bulk VBM.

Several a posteriori correction schemes have been developed to tackle this issue [26, 28, 34, 35, 29, 36, 37]. Most of the correction schemes involve solving the Poisson equation for a model system, and aligning of potentials, to derive the correction to the formation energy. In particular, the scheme proposed by Freysoldt, Neugebauer and Van de Walle (FNV) has gained a lot of popularity owing to its consistency in deriving accurate corrections for charged defects in numerous materials [26]. This scheme has been extended to low-dimensional systems as well, and shown to perform well [38, 39, 40, 41, 42, 43, 44]. However, a generalized correction scheme implementation that works with bulk as well as low-dimensional systems is absent in the various DFT packages, or as an independent package [45, 46]. The recent independent packages are restricted to compute the corrections for charged defects in bulk systems alone [45, 46].

We present a complete electrostatic corrections package, CoFFEE: Corrections For Formation Energy and Eigenvalues for charged defect simulations. The package is applicable to charged defects in materials ranging bulk solids, interfaces, surfaces/slabs, two-dimensional (2D) materials, nanowires and nanoribbons. These materials can be classified according to the number of periodic directions, into 3D (bulk) , 2D (slabs, 2D materials) and 1D (nanowires, nanoribbons) systems. We implement a generalized Poisson solver based on the FNV correction scheme with a gaussian model charge distribution. Tools to compute the potential alignment terms in the FNV correction scheme are also provided with the package. The code is written entirely in Python [47]. We use Message Passing Interface (MPI) to parallelize the code and Cython [48] to accelerate slow steps. Our implementation can be used alongside any DFT package to obtain an a posteriori correction for the formation energy and the defect level position in the gap, for the charged defect being simulated.

2 Theoretical framework

The formation energy of a neutral defect in a material is given by:

[TABLE]

where $\epsilon_{F}$ is the Fermi level, $\mathrm{E}^{\mathrm{tot}}_{0}[\vec{\mathrm{R}}_{0}]$ is the total energy of a system containing a neutral defect with atom positions at $\vec{\mathrm{R}}_{0}$. $\mathrm{E}_{\mathrm{pristine}}$ represents the total energy of a pristine super cell of the same size. $n_{x}$ refers to the number of atoms of type $x$ added (positive) or removed (negative) from the pristine system, and $\mu_{x}$ is the atom’s chemical potential [27]. The neutral defect formation energy needs no electrostatic correction term due to the absence of any long range electrostatic defect-defect interactions.

The formation energy of a defect in charge state $q$ is given by [27]:

[TABLE]

The first term on the right hand side is the total energy of a system containing a defect in charge state q with the requisite finite-size electrostatic correction $\mathrm{E}^{\mathrm{corr}}_{q}$ (described below). The formation energy is now a function of the Fermi level in the system, $\epsilon_{F}$, with respect to the pristine valence band maximum (VBM), $\epsilon_{\mathrm{vbm}}^{\mathrm{pristine}}$.

[TABLE]

is a potential alignment term found by comparing the electrostatic potentials from a pristine calculation and far from the defect in a neutral defect calculation.

2.1 FNV correction scheme

The electrostatic correction term, $\mathrm{E}^{\mathrm{corr}}_{q}$, is incorporated to correct the spurious interaction between the defect charge and its periodic images. In the FNV scheme, this term is given by [27, 26]:

[TABLE]

where $\mathrm{E}^{\mathrm{lat}}_{q}=\mathrm{E}^{\mathrm{iso,m}}_{q}-\mathrm{E}^{\mathrm{per,m}}_{q}$ is obtained from a model calculation. It involves solving the Poisson equation using a model charge distribution, $\rho^{m}(\mathbf{r})$, and model dielectric profile under periodic boundary conditions to obtain the potential, $V^{\mathrm{per,m}}_{q}(\mathbf{r})$ [27, 26]. $\mathrm{E}^{\mathrm{per,m}}_{q}$ is then given by:

[TABLE]

where the integral is over the super cell volume, $\Omega$. $\mathrm{E}^{\mathrm{per,m}}_{q}$ is evaluated for larger super cells and extrapolated to obtain $\mathrm{E}^{\mathrm{iso,m}}_{q}$. Performing larger super cell calculations within this model is computationally inexpensive compared to a DFT calculation on such systems. $\mathrm{E}^{\mathrm{lat}}_{q}$ accounts for the long range interactions. The dependence of $\mathrm{E}^{\mathrm{corr}}_{q}$ on the model charge distribution is eliminated by the second term in 4. $\Delta V_{q-0/m}$ is a potential alignment term found by comparing the model potential to the DFT difference potential:

[TABLE]

2.2 Model calculation

The model super cell calculation for $\mathrm{E}^{\mathrm{per,m}}_{q}$ involves solving the Poisson equation for the periodic model potential, $V^{\mathrm{per,m}}_{q}(\mathbf{r})$ [39, 38, 42]:

[TABLE]

where $\mathbf{\varepsilon}(\mathbf{r})$ is the dielectric tensor profile of the material, $\rho^{m}(\mathbf{r})$ is the model charge distribution. The dielectric tensor can be obtained using Density Functional Perturbation Theory (DFPT) [49, 50]. Eqn 7 can be effectively solved in the reciprocal space [39, 38, 42]:

[TABLE]

where $\varepsilon_{ii}$ are the diagonal terms of the dielectric tensor. $V^{\mathrm{per,m}}_{q}(\mathbf{G^{\prime}=0})$ is set to zero, which is equivalent to introducing a uniform, neutralizing background charge. The number of G vectors used in the calculation is determined by an energy cut off.

As one is interested in the long range corrections, $\varepsilon$ in the material can be assumed to have no spatial profile; however, in general, it is a 3$\times$3 tensor. In a 3D bulk system with an isotropic dielectric, $V^{\mathrm{per,m}}_{q}$ can be obtained from Eqn 8 as:

[TABLE]

For slab or two-dimensional systems, on the other hand, the dielectric tensor, while not having any profile inside the material, would still have a spatial profile in the aperiodic direction, say $z$ [39, 38].

$\begin{bmatrix}\varepsilon^{1}_{\parallel}(z)&0.0&0.0\\ 0.0&\varepsilon^{2}_{\parallel}(z)&0.0\\ 0.0&0.0&\varepsilon_{\perp}(z)\\ \end{bmatrix}$

where $\varepsilon^{1}_{\parallel}(z)$ and $\varepsilon^{2}_{\parallel}(z)$ is the dielectric profile for the in-plane dielectric constants $\vec{a}_{1}$ and $\vec{a}_{2}$ and $\varepsilon_{\perp}(z)$ is the dielectric profile for the out-of-plane dielectric constant of the material.

For slab systems or a 2D material with more than one atom thickness, like transition metal dichalcogenides, phosphorene, etc., the dielectric profile used is of the form, [39, 38]:

[TABLE]

where $t_{1}$ and $t_{2}$ are the edges of the slab in the $z$-direction in the simulation cell. A parameter, $s$, is used with the error function (erf) to smoothen the profile at the slab edges. $\varepsilon^{m}$ is the dielectric constant of the material, found using DFPT [39]. $\varepsilon^{v}$ is the dielectric constant of the space surrounding the material, set to 1 for vacuum. A sample slab profile is shown in Fig 1 (a). For a 2D material with one atom thickness, like Boron Nitride (BN), Silicon Carbide (SiC), graphene, etc., the dielectric profile used is a gaussian to mimic the planar-averaged charge density of the material in the $z$-direction. A sample profile is shown in Fig 1 (b).

For slab or 2D systems, Eqn 8 can then be written as:

[TABLE]

The average potential, $V^{\mathrm{per,m}}_{q}(\mathbf{G}=0)$, is set to zero, to introduce a neutralizing bakground charge.

For a quasi-1D system like nanowires or nanoribbons, with perodicity along the $z$ direction, the dielectric tensor profile is of the form [42] :

$\begin{bmatrix}\varepsilon_{\perp}(x,y)&0.0&0.0\\ 0.0&\varepsilon_{\perp}(x,y)&0.0\\ 0.0&0.0&\varepsilon_{\parallel}(x,y)\\ \end{bmatrix}$

where $\varepsilon_{\perp}(x,y)$ is the profile for the dielectric constant perpendicular to the wire and $\varepsilon_{\parallel}(x,y)$ is the profile for the dielectric constant along the wire. A sample profile for a silicon nanowire along the [111] direction is shown in Fig. 2. The wire, oriented along the $z$ direction, has a hexagonal cross-section [42]. The profile is constructed such that the points inside the wire have the dielectric constant of the material, $\varepsilon^{m}$, and the points outside have $\varepsilon^{v}$. A sample profile for a BN nanoribbon is shown in Fig. 3. The profile is constructed to mimic the electron charge density of the material by combining a slab like profile (Eqn 10) along the width of the nanoribbon and a gaussian along the out-of-plane direction. Eqn 8, for 1D systems, can then be written as [42]:

[TABLE]

We set $V(\mathbf{G}=0)=0$ here, to introduce a neutralizing background charge.

2.3 Defect eigenvalue correction

The defect eigenvalues, like the formation energy, show slow convergence with the super cell size [32, 33]. The correction to the eigenvalue is given by [33]:

[TABLE]

where $\mathrm{E}_{q}^{\mathrm{corr}}$ is the FNV correction term, Eqn 4. The correction is negative for positively charged defects and positive for negatively charged defects. Defect level corrections are crucial in studying defects using the DFT+GW formalism [32, 33, 2, 3, 1] and also in interpreting absorption experiments.

3 Test Systems

We show the performance of our code in correcting the formation energies and defect levels in the following systems:

Bulk: Carbon vacancy in diamond in charge state -2, $V_{C}^{-2}$. 2. 2.

2D: Sulfur vacancy in monolayer MoS2 in charge state -1, $V_{S}^{-1}$. 3. 3.

1D: Boron vacancy in BN nanoribbon in charge state -1, $V_{B}^{-1}$.

All the DFT calculations are performed using the plane-wave, pseudopotential package, Quantum Espresso [51]. We perform simulations on a range of super cell sizes and compare the extrapolated formation energy and defect level (to the isolated limit) with the corrected formation energy and defect level.

3.1 Vacancy in diamond, $V_{C}^{-2}$

We perform DFT calculations on different cubic super cell sizes with number of C atoms ranging from 64 to 1024. We use norm-conserving pseudopotentials and the LDA exchange correlation functional [52]. The wavefunctions are expanded in plane-waves upto an energy cut off of 60 Ry. For the 64 atom, 2$\times$2$\times$2 super cell, a k-point sampling of 5$\times$5$\times$5 was used. An equivalent sampling is used for the other super cell sizes. We did not include any atomic relaxations.

The model calculation is performed using the CoFFEE code, solving Eqn 9. The defect is modelled with a gaussian of width 2.6 bohr, a dielectric constant of 5.76 and a plane wave energy cut off of 16 Ry. Fig. 4 (a) shows the defect wavefunction charge density, $|\psi_{d}|^{2}$, and the model charge density. 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. Note also that the plane wave energy cut off required to converge the model calculation is inversely proportional to the width of the Gaussian. Fig. 4 (c) shows the extrapolation of the model energy, $\mathrm{E}^{\mathrm{per,m}}_{-2}$ (Eqn 5), with super cell size, $\Omega$. The energy is fit with a polynomial of the form:

[TABLE]

$f_{1}$ here corresponds to the isolated model energy, $\mathrm{E}^{\mathrm{iso,m}}_{-2}$. The fitting parameters here are: $f_{1}=2.04$, $f_{2}=-14.19$ and $f_{3}=120.85$. The lattice correction is then given by: $\mathrm{E}^{\mathrm{iso,m}}_{-2}-\mathrm{E}^{\mathrm{per,m}}_{-2}$ (Eqn 4). Fig. 4 (b) compares the planar averaged model potential and the DFT difference potential which contribute to the potential alignment term, $\Delta V_{q-0/m}$ (Eqn 6) for a 3$\times$3$\times$3 super cell. The defect in this calculation is at the center of the cell. The model potential matches well with the DFT difference potential far from the defect, leading to a very small correction. On the other hand, the potential alignment term, $\Delta V_{0/p}$ in Eqn 2, which compares the potential far from the neutral defect to the pristine, has a substantial contribution.

The uncorrected formation energy of $V_{C}^{-2}$ for various supercell sizes is fit using Eqn 14 to extrapolate to the isolated limit as shown in Fig. 4 (d). The chemical potential of carbon is taken from bulk diamond. The FNV corrected formation energy shows good agreement with the extrapolated value for super cell sizes larger than 2$\times$2$\times$2. The correction scheme performs well when the defect wavefunction is localized well within the super cell; this is not the case with the 2$\times$2$\times$2 super cell. The defect level accomodating the -2 charge is shown in Fig. 5 with respect to the pristine VBM eigenvalue in the charged super cell. It shows a similar scaling with super cell size and is corrected using Eqn 13.

3.2 Vacancy in monolayer MoS2, $V_{S}^{-1}$

We study sulfur vacancies in monolayer MoS2 to demonstrate the electrostatic corrections in 2D or slab systems. The DFT calculations are performed using PAW pseudoptentials [53], PBE scheme for the exchange correlation functional [54]. A wavefunction cut off of 50 Ry, and charge density cut off of 500 Ry is used. We perform calculations on four super cell sizes: $\alpha\times\alpha\times\alpha$, for $\alpha$ = 4, 5, 6 and 8. The $\alpha$ refers to the scaling of the cell dimension. It has been reported that scaling just the in-plane supercell, keeping the amount of vacuum fixed, leads to a divergence in the model energy [38]. We hence uniformly scale the vacuum with the in-plane super cell. For $\alpha=6$, the vacuum thickness is $\sim 16\mathrm{\AA}$ and thickness of MoS2 is $3.2\mathrm{\AA}$. The k-point sampling for these is chosen to effectively sample the unit cell Brillouin zone with a grid better than 12$\times$12$\times$1. The atoms are relaxed to their equilibrium positions such that the force on each atom is $\leq 10^{-2}$ eV/Å.

The model calculation is performed solving Eqn 11. The dielectric profile is chosen as in Eqn 10, with $\varepsilon_{m}=15$ in the direction parallel to the material and $\varepsilon_{m}=2$ in the perpendicular direction. The smoothening parameter, s, is chosen to be 0.38 bohr. The dielectric constants are determined using DFPT as prescribed in reference [39]. The plane wave energy cut off is set to 30.0 Ryd. The defect is modelled with a Gaussian of width 1.9 bohr as shown in Fig. 6 (a). It is important here to choose a Gaussian width that mimics the defect charge density. Choosing a width too large that spills beyond the gray shaded region in Fig. 6 (a) would affect the results. The defect wavefunction charge density can be fit to a Gaussian using the script g_fit.py supplied with the CoFFEE code.

Fig 6 (b) compares the planar averaged model potential with the planar averaged DFT difference potential, $V_{-1}^{\mathrm{DFT}}-V_{0}^{\mathrm{DFT}}$, for the $\alpha=6$ supercell. The DFT difference potential used here is from unrelaxed calculations and shows good agreement with the model potential far from the defect. Fig 6 (c) shows the variation of the model electrostatic energy with two different scalings of the super cell size: $\alpha\times\alpha\times\alpha$ and $\alpha\times\alpha\times 1.6\alpha$. For $\alpha\times\alpha\times 1.6\alpha$ super cell scaling, the vacuum thickness is 1.6 times larger than that of $\alpha\times\alpha\times\alpha$ scaling. The interaction between the periodic charges is unscreened in the out-of-plane direction, and screened by the material in the in-plane direction. The model energy is thus strongly dependent on the vacuum thickness, for small super cell sizes. As expected, at the infinite limit ($\alpha=\infty$), the two scalings extrapolate to the same value. The model scaling plots clearly do not follow a simple fit as in the case of bulk solids due to the spatially varying dielectric profile. The energy is fit with a polynomial of degree five in $\alpha^{-1}$: $\mathrm{E}^{\mathrm{per,m}}_{-1}(\alpha)=f_{0}+f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3}+f_{4}\alpha^{-4}+f_{5}\alpha^{-5}$ [39]. The isolated model energy is given by $f_{0}=\mathrm{E}^{\mathrm{iso,m}}_{-1}$, and the lattice correction is given by: $\mathrm{E}^{\mathrm{iso,m}}_{-1}-\mathrm{E}^{\mathrm{per,m}}_{-1}$. The fitting parameters, for $\alpha\times\alpha\times\alpha$ scaling, are found to be: $f_{0}=0.66$, $f_{1}=-4.04$, $f_{2}=54.95$, $f_{3}=-388.69$, $f_{4}=1267.38$ and $f_{5}=-1579.27$.

Fig 7 (a) shows the scaling of the uncorrected formation energy of neutral S vacancy, $V_{S}^{0}$, and -1 charged S vacancy, $V_{S}^{-1}$, as a function of super cell size. The chemical potential of sulfur is taken from the cyclo-S8 allotrope of sulfur. The neutral vacancy formation energy is well converged and shows no scaling with system size. The scaling of the charged vacancy formation energy is fit with a polynomial of the form [39]: $\mathrm{E}_{-1}^{f}(\alpha)=t_{0}+(f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3}+f_{4}\alpha^{-4}+f_{5}\alpha^{-5})+t_{3}\alpha^{-3}$. The additional terms are necessary to capture the scaling of the potential alignment term with system size. $t_{0}$ and $t_{3}$ are found to be 4.71 and -8.69 respectively. The FNV corrected formation energy shows excellent agreement with the extrapolated value for all the super cell sizes considered. Fig 7 (b) shows the dependence of the defect level accomodating the -1 charge with system size. The uncorrected values are fit with a polynomial of the form: $\epsilon_{-1}(\alpha)=t_{0}+2\times(f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3}+f_{4}\alpha^{-4}+f_{5}\alpha^{-5})+t_{3}\alpha^{-3}$. The corrected eigenvalues are in good agreement with the extrapolated value.

3.3 Vacancy in BN nanoribbon, $V_{B}^{-1}$

To show the application of the code for 1D systems, we study B vacancy in a BN nanoribbon of thickness 13.39 Å[13]. We perform calculations on three super cell sizes: $\alpha$ (1$\times$2.31$\times$1), for $\alpha=6,8,10$. For $\alpha=6$, the simulation cell dimensions are 15$\mathrm{\AA}$, 35$\mathrm{\AA}$ and 15$\mathrm{\AA}$, in the $x$, $y$ and $z$ directions repectively. The $z$ direction is the periodic direction. A vacuum of 15$\mathrm{\AA}$ and 22$\mathrm{\AA}$ has been introduced in the out-of-plane $x$ direction and lateral $y$ direction. The thickness of the ribbon is fixed in these cell sizes, with the number of atoms in the periodic direction and the vacuum padding scaling with $\alpha$. The ribbon is passivated on either side with H atoms. Ribbon in $\alpha=6$ super cell with a B vacancy is shown in Fig 8 (a). The DFT calculations are performed using norm conserving pseudopotentials and the PBE scheme is used for the exchange correlation functional [54]. A plane wave energy cut off of 70 Ry is used for the wavefunctions. The Brillouin zone is sampled with a 1$\times$1$\times$2 grid for $\alpha=6,8$ and with the $\Gamma$ point for $\alpha=10$. We relaxed the atoms in the super cell containing the defect until the force on each atom is $\leq 10^{-2}$ eV/Å.

For the model calculation, we use a dielectric profile as shown in Fig 3. A slab like profile, similar to the one used for MoS2, is used along the width of the nanoribbon and a gaussian profile in the out-of-plane, $x$ direction. The width of the slab is taken to be that of the width of the nanoribbon, with a smoothening parameter of 0.37 bohr. The gaussian width is taken to be 0.86 bohr, chosen to mimic the electron charge density in that direction. $\epsilon^{m}$ used is 2.9 along the $x$ direction and 12 along the other two directions. These are computed from DFPT calculations on a monolayer BN sheet. We use a gaussian model charge with width 1.89 bohr and integrated charge -1. The planar averaged model charge and the defect wavefunction charge density are shown in Fig. 8 (b). The calculation is performed with a plane wave energy cut off of 16.0 Ry. Fig 8 (e) shows the scaling of the computed $\mathrm{E}^{\mathrm{per,m}}_{-1}$ as a function of 1/$\alpha$. The points are fit with a third degree polynomial in $\alpha^{-1}$: $\mathrm{E}^{\mathrm{per,m}}_{-1}(\alpha)=f_{0}+f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3}$. $f_{0}$ then corresponds to the isolated model energy, $\mathrm{E}^{\mathrm{iso,m}}_{-1}$ , and the lattice correction is given by: $\mathrm{E}^{\mathrm{iso,m}}_{-1}-\mathrm{E}^{\mathrm{per,m}}_{-1}$. The fitting parameters are found to be: $f_{0}=1.90$, $f_{1}=-1.85$, $f_{2}=-4.35$ and $f_{3}=0.58$. The potential alignment terms, $\Delta V_{q-0/m}$ and $\Delta V_{0/p}$, are calculated by comparing the potentials in the out-of-plane, $x$ direction, far from the defect. Both these terms are found to be small ($\textless$ 20 meV). Fig 8 (c) and (d) show the DFT difference potential in Eqn 6 and the model potential, respectively, for $\alpha=6$ super cell size. The defect is at (x, y) = (7.5, 9) in the figure. The model potential matches well with the DFT difference potential far from the defect.

Fig 9 (a) shows the scaling of the uncorrected formation energy of neutral vacancy, $V_{B}^{0}$, and -1 charged vacancy, $V_{B}^{-1}$, as a function of $1/\alpha$. The formation energy is computed for nitrogen-rich conditions. The chemical potential for N is taken from N2 molecule. Formation energy of the neutral vacancy is well converged and shows no scaling with the system size. The uncorrected formation energy of the charged vacancy is fit with a polynomial of the form: $\mathrm{E}_{-1}^{f}(\alpha)=t_{0}+(f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3})+t_{3}\alpha^{-3}$. $t_{0}$ and $t_{3}$ are found to be 9.39, -3.92. The FNV corrected formation energy shows excellent agreement with the extrapolated value for all the super cell sizes considered. Fig 9 (b) shows the dependence of the defect level accomodating the -1 charge with system size. The uncorrected values are fit with a polynomial of the form: $\epsilon_{-1}(\alpha)=t_{0}+2\times(f_{1}\alpha^{-1}+f_{2}\alpha^{-2}+f_{3}\alpha^{-3})+t_{3}\alpha^{-3}$. The corrected eigenvalues are in good agreement with the extrapolated value.

4 Workflow

The general steps involved in computing corrections for the formation energy of a charged defect are the following (Fig. 10):

Compute the total energy of the pristine super cell of the same size. Save the DFT potential in cube/xsf format. 2. 2.

Compute the total energy of the super cell (say n$\times$n$\times$n) containing the neutral defect. Save the DFT potential in cube/xsf format. 3. 3.

Compute the total energy of the super cell containing the charged defect. Save the DFT potential in cube/xsf format. 4. 4.

Compute $\mathrm{E}^{\mathrm{lat}}_{q}$ term: Compute the model energy for various super cell sizes and extrapolate to obtain $\mathrm{E}^{\mathrm{iso,m}}_{q}$. $\mathrm{E}^{\mathrm{lat}}_{q}$ is then given by: $\mathrm{E}^{\mathrm{iso,m}}_{q}-\mathrm{E}^{\mathrm{per,m}}_{q}$ (n$\times$n$\times$n). 5. 5.

Compute the potential alignment term $\Delta V_{0/p}$, Eqn. 3. The utility script dV_0p.py can be used to compute this. 6. 6.

Compute the potential alignment term $\Delta V_{q-0/m}$, Eqn. 6. The utility dV_mD.py can be used to compute this.

The first three steps are performed by the user with the DFT electronic structure code of his/her choice. The xsf/cube file formats are commonly used to visualize data. Most electronic structure codes provide utilities to convert the DFT potential after a self consistent calculation into these formats. These formats act as an interface between the DFT calculation and the CoFFEE code. The CoFFEE code aids in computing steps 4, 5 and 6. The model calculations are performed by solving the Poisson equation, as detailed above. Details on preparing the input file and running the CoFFEE code to perform these calculations are provided to some extent in the next section and exhaustively in the user guide of the code. The user guide is available for download with the distribution.

5 CoFFEE code framework

5.1 Layout

On unzipping the tar file after download, the CoFFEE folder contains the following directories: PoissonSolver, PotentialAlignment and Examples, and a script coffee.py. The script coffee.py is the main executable. It can be called from the user’s working directory and it performs the model calculation by sequentially calling the relevant functions as depicted in Fig 11. The script coffee.py reads input from a file. The input file contains parameters pertaining to the dimensions of the super cell, the dielectric profile and the model Gaussian charge. The next subsection describes the input file parameters. On running the script, the model total energy is printed out. The model potential, dielectric profile and model charge distribution are written as numpy save files on providing the appropriate flags in the input. The script is run for different supercell sizes and the model energy extrapolated as shown in Fig. 4 (c), 6 (c) and 8 (e) to obtain the isolated model energy.

The PoissonSolver folder holds classes.py which defines three classes: cell, diel_profile and gaussian. Each class has a function to read the parameters pertaining to it from the input file. The cell class has information regarding the cell parameters and the plane wave cut off to be used in the calculation. The diel_profile class has attributes regarding the dielectric profile to be used and functions to construct and Fourier transform the requisite profile. The gaussian class has attributes regarding the construction and Fourier transform of the model gaussian charge. The PoissonSolver folder also consists of three routines, Solver.py, construct_eps.py and PS_main.py. construct_eps.py is used to select and construct the appropriate dielectric profile, based on the user input. Solver.py is used to select the poisson solver to be run. PS_main.py contains the solvers for 1D (PS_1D), 2D (PS_2D) and bulk (PS_3D) systems.

5.2 Input file

The input file is divided into three sections designated with ”&CELL_PARAMETERS”, ”&DIELECTRIC_PARAMETERS” and ”&GAUSSIAN_PARAMETERS”. Each section is ended with a ”/” and contains, in the intervening space, the parameters relevant to the respective section. The following is an example of an input file for $V_{C}^{-2}$ in diamond.

&CELL_PARAMETERS

Normalized lattice vectors: a1, a2 and a3

Lattice_Vectors(normalized): 1.000000000 0.000000000 0.000000000 0.000000000 1.000000000 0.000000000 0.000000000 0.000000000 1.000000000

Cell dimensions. Provide "angstrom" in place

of "bohr" if you wish to specify

these values in angstrom units.

These are multiplied to a1, a2 and a3 respectively.

Cell_dimensions bohr 26.594331775231996 26.594331775231996 26.594331775231996

G-vectors will be used upto this kinetic energy cut off.

Provide "Rydberg" in place of "Hartree"

if you wish to specify the cut off in Rydberg

atomic units.

Ecut=20.0 Hartree /

Set "Bulk" here for 3D, bulk systems.

&DIELECTRIC_PARAMETERS Bulk

Sets the value of the dielectric constant along a1, a2, a3.

Epsilon1_a1 = 5.76 Epsilon1_a2 = 5.76 Epsilon1_a3 = 5.76 /

&GAUSSIAN_PARAMETERS:

The charge state of the defect being simulated.

Total_charge = -2

The width of the model Gaussian charge being used.

(default: bohr units)

Sigma = 2.614

These set the center of the Gaussian in crystal

units.

Centre_a1 = 0.5 Centre_a2 = 0.5 Centre_a3 = 0.5 /

The default unit of length is bohr, unless ”angstrom” is explicitly specified after the number. Ecut stands for the plane wave energy cut off used in solving the Poisson equation. The string after &DIELECTRIC_PARAMETERS determines the dielectric profile and which Poisson solver is to be called. The dielectric parameters are provided in the input file along the crystal axes. The center of the model charge is specified in crystal units in the &GAUSSIAN_PARAMETERS, along with the charge it carries and the gaussian width. For the nanowire profile, the user can provide a list of points that form the boundary of the cross-section of the wire in a file, based on which the profile is created. The user guide lists all the possible input parameters and the properties they control. Also, the Examples folder has sample profiles for different materials.

5.3 Potential alignment

The PotentialAlignment directory holds two scripts: dV_0p.py and dV_mD.py to compute the two potential alignment terms, $\Delta V_{0/p}$ and $\Delta V_{q-0/m}$ in Eqn. 2 and 4. These scripts read the DFT potential from ”cube” or ”xsf” formats. dV_0p.py plots the planar averaged $V_{0}-V_{p}$, in the desired direction. The difference far from the defect can then be read from the plot. This script reads input from a file. The format of the input file:

&dV_0p file_type = cube # No quotes. Takes cube/xsf file_neutral = n.cube # No quotes. Path to the neutral DFT potential file file_pristine = p.cube # No quotes. Path to the pristine DFT potential file plt_dir = a1 # No quotes. Takes a1/a2/a3. If a1 # is specified, the data is averaged along a2 # and a3 directions and the planar averaged data is plotted # along a1 in a file pa_dv0p_a1.plot factor = Ryd # factor to be multiplied to the cube/xsf data. If the data is in # rydberg and the plot is needed in eV, specify # factor = Ryd. If the data is in Hartree units, # specify factor = Hartree /

dV_mD.py plots the planar averaged $(V_{q}^{\mathrm{DFT}}-V_{0}^{\mathrm{DFT}})$ and $V^{\mathrm{per,m}}_{q}|_{\mathrm{far}}$ along the desired direction. These can be plotted as shown in Fig. 4 b, 6 b. This script reads input from a file. The format of the input file:

&dV_mD file_type = cube # No quotes. Format of the DFT potential files: cube/xsf file_model = m.npy # No quotes. Path to the model potential file (.npy) file_charged = q.cube # No quotes. Path to the charged DFT potential file file_neutral = n.cube # No quotes. Path to the neutral DFT potential file plt_dir = a1 # No quotes. Takes a1/a2/a3. If a1 # is specified, the data is averaged along a2 # and a3 directions and the planar averaged data is plotted # along a1 in files DFTdiff_a1.plot, model_a1.plot factor = Ryd # factor to be multiplied to the cube/xsf data. If the data is in # rydberg and the plot is needed in eV, specify # factor = Ryd. If the data is in Hartree units, # specify factor = Hartree /

5.4 Solver parallelization and optimization

PS_1D, PS_2D and PS_3D are the poisson solver functions in PS_main.py for 1D, 2D and bulk systems. We do not use any parallelization for bulk systems, where the potential is obtained using Eqn 9, since the computation time is fairly small. We have, however, optimized this function with the help of Cython [48]. For 2D systems, linear equations of the form in Eqn 11 are solved. The number of reciprocal lattice vectors depends on the plane wave energy cut off set for the calculation. For $N_{G_{1}}$ and $N_{G_{2}}$ reciprocal lattice vectors along the in-plane reciprocal lattice directions, a set of $N_{G_{1}}\times N_{G_{2}}$ linear equations are solved. We use the iterative solver bicgstab() which is a part of the SciPy linear algebra package to solve these linear equations. These calculations are computationally intensive for large system sizes. We parallelize the total number of linear equations so that each process solves $N_{G_{2}}\times N_{G_{1}}/N_{p}$ equations, where $N_{p}$ is the number of processes. The parallelization is done via MPI, using the package mpi4py [55]. For 1D systems, $N_{G_{z}}$ linear equations are solved as shown in Eqn 14 . These are again parallelized in the package in a similar manner. Furthermore, the efficiency of the bicgstab() routine is primarily determined by the time taken to compute matrix vector products. We optimized these products using Cython [48] to reduce the computation time.

6 Conclusion

We present a complete package, CoFFEE, for electrostatic corrections in charged defect simulations. We demonstrate the application of this code on three test systems, namely: bulk diamond, 2D MoS2 and 1D BN nanoribbon. The corrected formation energy and defect eigenvalues for these systems are found to be in good agreement with the isolated limit extrapolation of the corresponding uncorrected quantities. The code, written completely in Python [47], is parallelized using MPI and the slowest steps accelarated using Cython [48].

7 Acknowledgments

We thank the Supercomputer Education and Research Centre (SERC) at IISc for providing the computational facilities.