Charged defects are often studied within the periodic density functional theory (DFT), but this introduces strong finite-size artifacts. In this work, we develop an electrostatic image interaction correction (IIC) method based on the direct solution of the Poisson equation for charge models constructed directly from DFT calculations. These IICs are found to be detail-insensitive, depending almost entirely on bulk dielectric properties. As these IICs are not able to fully explain the observed finite-size scaling, we explore potential alignment in detail and introduce a novel decomposition to separate out different contributions. We find that the two main sources of potential alignment are defect image interactions and changes in the number of atoms present in the supercell. This first effect is accurately predicted by the periodic part of our IIC. The second contribution is unrelated to the IIC and justifies the common observation that the magnitude of finite-size dependence can strongly vary between vacancy and interstitial defects. It can be approximately predicted using atomic radius, but is strongly sensitive to the pseudopotential employed. Combined, these developments provide a new justification for known finite-size scaling rules. Our results suggest that for cubic supercells, the Lany-Zunger IIC, combined with simplified potential alignment between neutral systems, can yield accurate corrections in spite of the simplicity of the approach.

## I. INTRODUCTION

Accurate defect formation energies are vital to understanding the properties of materials and electronic devices, as they are required to calculate defect concentrations and charge transition levels.^{1} Defect calculations are often carried out using periodic boundary conditions in classical simulations^{2} or at different levels of many-electron theory, notably density functional theory (DFT). Periodic translations allow defects to interact with their own images in the neighboring cells of the simulation. This leads to strong electrostatic interactions between a defect and its images when the defect has net charge.^{3} The Coulomb interaction is long ranged, and these image interactions may remain significant even when relatively large supercells are used. The use of periodic boundary conditions also creates an additional problem that the bulk electrostatic potential is not properly recovered far from the defect.^{4}

In practice, it is not possible to use sufficiently large supercells in which these effects can be neglected, especially with the high computational costs associated with hybrid functionals, now in common usage.^{5} This problem is typically solved by applying *post hoc* corrections to the final energy. In this way, highly developed standard DFT codes can still be employed at a reasonable computational cost.

*Post hoc* correction methods have been developing over the last 25 years, and there are now many approaches available and compared in the literature (see, e.g., Refs. 3 and 6–13 which is by no means a complete list). Detailed comparison of different methods is complicated by a realization that image interactions contribute to the need of potential alignment, meaning that the two problems must be considered together.^{8–10} It seems likely that some of the suggested techniques are reliable, but many disagreements in the literature mean that the current position is rather unclear, especially for non-experts. This is perhaps one of the reasons why it has recently been suggested that formation energies can only be rigorously calculated for neutral defects.^{14}

By contrast, the actual finite-size scaling of charged defect formation energies has been observed to exhibit some common features across many kinds of materials and defects. For cubic systems, for example, the main errors in the formation energy are observed to scale in terms of the length and volume of the supercell.^{15} Furthermore, vacancy and interstitial defects in the same charge state often appear to demonstrate quite different finite-size scaling.^{10,13} The origins of these observations are not yet fully understood. In this paper, we attempt to separate image interaction corrections and potential alignment corrections and to use this approach to provide a new rational for these observations.

First, we introduce a new image interaction correction (IIC), where we solve the Poisson equation directly for the additional charge present in the DFT supercell. This enables us to explore the degree of finite-size dependence that is introduced by electrostatic interactions with images and the inclusion of a jellium background, without recourse to an approximate description of the image charge distribution. Such an approach predicts that these image interactions introduce a length dependent finite-size error, for the case of cubic supercells.

Next, we consider potential alignment and introduce a new decomposition of the total potential alignment into three parts. We find that the image correction we constructed exactly predicts one of these parts, leaving two other components that are unrelated to image interactions. These alignment components are introduced by the removal (or addition) of atoms from the cell, and ionic relaxation of the cell, respectively. We observe that the first component can be predicted by calculations on isolated atoms and that the potential alignment correction required for isolated atoms scales with the volume of the supercell. Additionally, this potential alignment correction has an opposite sign between vacancies and interstitials. It is observed that inclusion of the second component often introduces additional finite-size dependence.

We apply these methods to vacancy defects in three materials, diamond, MgO, and SnO_{2}, where our examples are chosen both to represent materials with different types of bonding and defect states that are quantitatively different. The relative magnitude of the IICs and potential alignment corrections is substantively different between these examples. However, the two corrections outlined above are shown to be sufficient, without any empirical parameters.

The rest of the paper is organized as follows. In Sec. II, we give a brief overview of finite-size dependence and some of the existing approaches to calculating image charge corrections. Then in Sec. III, we propose and justify a new image charge correction method. In Sec. IV, we introduce a new decomposition of the potential alignment. In Sec. V, we apply these developments to defect calculations in diamond, MgO, and SnO_{2} and compare our results with those obtained using several common correction methods. The discussion and conclusions of our findings are provided in Sec. VI.

## II. CALCULATION OF DEFECT FORMATION ENERGIES

### A. Cell size dependence

One standard definition^{16} of the defect formation energy *E*_{f} is given as

where $EdefectDFT$ and $EbulkDFT$ are the total energies of the supercell containing the defect and an equivalently shaped supercell of the perfect bulk crystal, respectively. Additionally, this raw energy difference is modified by the energy required to exchange *n*_{i} ions with some reservoir of chemical potential *μ*_{i}. The first three terms of Eq. (1) are often fully adequate for calculating the formation energy of neutral defects, although it neglects other sources of errors such as elastic relaxation or defect state hybridization. Additional terms which describe these effects have been suggested in Refs. 10 and 17.

When a defect has net charge *q*, additional terms appear in Eq. (1). First, the energy cost to exchange *q* electrons with an electron chemical potential *μ*_{e} must be applied. We follow the usual convention and use the valence band maximum (VBM) of the bulk crystal as our chemical potential such that *μ*_{e} = *E*_{VBM}. Second, the spurious interaction energy *E*_{IIC} due to the extra interaction of the defect both with its own charged images and the neutralizing jellium background (neutralizing charge is required in all charged periodic calculations to avoid a divergence of the total energy) must be removed. Finally, the average electrostatic potential cannot be directly calculated using periodic DFT and is conventionally set to be zero.^{18–20} This neglect of the average potential (NAP) introduces a numerical error into the formation energy of the charged defect. A potential alignment correction *q*Δ*V*_{NAP} is often applied to account for this problem.

The two terms *E*_{IIC} and Δ*V*_{NAP} are finite-size corrections, which tend to zero as the size of the DFT simulation cell is increased, and are the focus of this paper. Their performance can be assessed by exploring the convergence of *E*_{f} with increasing supercell size. Ideally, once corrected, there would no longer be a supercell size dependence in *E*_{f}. Hence, the simplest measure of the success of these corrections is the supercell size dependence of the corrected energies.

Additionally, detailed studies on the supercell size dependence of the formation energy of charged defects have been carried out. The study of Castleton *et al.*^{15} has produced the following empirical scaling rule for cubic supercells:

which demonstrates an inverse size dependence on both the length of the considered supercell, *L*, and the volume of the supercell Ω = *L*^{3}. It is particularly useful if the physical origins of terms *a*_{1} and *a*_{3} can be identified, and several plausible suggestions have been offered in the literature.^{10,13} [We note the similarity between Eqs. (2) and (4), which was discovered earlier.^{3} However, we will later argue that these two equations are distinct.]

The expression in Eq. (2) is often used to extrapolate charged formation energies to the dilute limit.^{15} This is then used to test the performance of charge correction methods in several papers.^{6,9–11} We do not report these extrapolations in this paper, as we believe that the uncertainties in these extrapolations are about the same as the differences between the correction methods considered. Extrapolating a three parameter model from four supercells provides one point that can be exuded from the fit, to check the consistency of the extrapolation, but we always observe a poor prediction of this point. As we will later justify (in Table I), this is because the *L*^{−3} scaling of the *a*_{3} term is only approximately obeyed, causing the accuracy of this extrapolation to vary between the different defects considered.

System . | Supercell . | ΔV_{AB} (V)
. | V_{atom} (V)
. | Error (%) . |
---|---|---|---|---|

Diamond (C vacancy) | 1 × 1 × 1 | −1.87 | −2.35 | −25.4 |

2 × 2 × 2 | −0.32 | −0.30 | 8.4 | |

3 × 3 × 3 | −0.10 | −0.09 | 15.7 | |

4 × 4 × 4 | −0.05 | −0.04 | 18.4 | |

MgO (O vacancy) | 1 × 1 × 1 | −0.98 | −1.24 | −26.7 |

2 × 2 × 2 | −0.12 | −0.15 | −32.7 | |

3 × 3 × 3 | −0.03 | 0.05 | −31.5 | |

4 × 4 × 4 | −0.02 | −0.02 | −20.5 | |

5 × 5 × 5 | −0.01 | −0.01 | −15.2 | |

SnO_{2} (Sn vacancy) | 1 × 1 × 1 | −3.73 | −3.44 | 7.9 |

2 × 2 × 2 | −0.66 | −0.59 | 10.0 | |

3 × 3 × 3 | −0.23 | −0.16 | 28.6 | |

4 × 4 × 4 | −0.10 | −0.07 | 32.7 |

System . | Supercell . | ΔV_{AB} (V)
. | V_{atom} (V)
. | Error (%) . |
---|---|---|---|---|

Diamond (C vacancy) | 1 × 1 × 1 | −1.87 | −2.35 | −25.4 |

2 × 2 × 2 | −0.32 | −0.30 | 8.4 | |

3 × 3 × 3 | −0.10 | −0.09 | 15.7 | |

4 × 4 × 4 | −0.05 | −0.04 | 18.4 | |

MgO (O vacancy) | 1 × 1 × 1 | −0.98 | −1.24 | −26.7 |

2 × 2 × 2 | −0.12 | −0.15 | −32.7 | |

3 × 3 × 3 | −0.03 | 0.05 | −31.5 | |

4 × 4 × 4 | −0.02 | −0.02 | −20.5 | |

5 × 5 × 5 | −0.01 | −0.01 | −15.2 | |

SnO_{2} (Sn vacancy) | 1 × 1 × 1 | −3.73 | −3.44 | 7.9 |

2 × 2 × 2 | −0.66 | −0.59 | 10.0 | |

3 × 3 × 3 | −0.23 | −0.16 | 28.6 | |

4 × 4 × 4 | −0.10 | −0.07 | 32.7 |

### B. Existing image charge corrections

As the use of periodic boundary conditions replaces an isolated charged defect with an array of periodic images, calculations using the supercell methodology also include an energy error (if one is interested in the dilute limit) due to the interactions with images, *E*_{IIC}. In order to remove these unwanted interactions, they can be modeled separately using classical electrostatics and then be removed. This is shown pictorially in Fig. 1.

The simplest electrostatic model that can be applied to construct an IIC is approximating the defect as a point charge (PC). In Fig. 1, this represents the charges *q* as ideal points. These corrections were first applied to defect calculations using classical molecular dynamics by Leslie and Gillan.^{2} In the context of DFT, this is the first-order Makov-Payne (MP) correction for cubic supercells,^{3}

where *ε* is the dielectric constant of the material containing the defect, *L* is the separation between defects, and *α*_{M} is the Madelung constant, which depends only on the shape of the supercell containing the defect. This provides the exact interaction energy of an array of point charges immersed in neutralizing jellium. The appropriate dielectric constant is the static dielectric constant *ε*_{0} when ionic relaxation of the defect is performed or the high frequency dielectric constant *ε*_{∞} when the ionic structure is kept fixed. This simple method has been extended to general supercell shapes and to anisotropic dielectrics.^{21} For cubic cells, this produces a *L*^{−1} dependent term that could explain the origin of the *a*_{1} term in Eq. (2).

However, a point charge is a highly idealized model, which is known to overestimate the size of the unwanted image interactions. Charged defects are often not point-like, but instead have charge distributed over some defect dependent regions, such as a vacancy site, dangling bonds, or neighboring ions. Additionally, the finite size of the supercell may affect the degree of dielectric screening.

Makov and Payne^{3} were motivated by these considerations to suggest a multipole expansion to generate higher-order terms that could also be calculated. Such an expansion yields

where *Q* is the quadrupole moment of the defect charge model. This expression contains only odd *L* terms because the even *L* terms are zero for charge distributions with cubic symmetry. Hence, Eq. (4) is correct up to the fifth-order for cubic supercells. However, the MP paper^{3} does not provide a systematic way of calculating *Q*, so this quadrupole term is rarely calculated within the MP method.

The Lany-Zunger (LZ) correction^{6,7} directly attacks this problem by making a practical approximation to *Q*. This method is derived by assuming that the defect draws bound screening charge, as demanded by classical electrostatics^{52,53}. This bound charge is drawn uniformly from the supercell, leading to a uniform charging given by

This charge present across the whole supercell is expected to be the dominant charge component far from the defect (note this is not the whole bound charge), is used to define the quadrupole moment *Q*, and leads to a reduction in the magnitude of the correction. In fact, this is shown to reduce the MP correction by about a third. Charge models constructed from DFT density differences are found to include this bound charge, as shown in Fig. 2(a).

Generally, we may express this reduction in the IIC by defining an internal polarization ratio *f*. The MP correction is too large and needs to be reduced by a scaling factor. Then, the complete interaction correction is given as

The LZ method provides a practical approximation to *f* in the form

where the shape factor, *c*_{sh}, depends only on the shape of the defect containing the supercell. These factors have been calculated for common supercells.^{7} As the LZ method is derived from Eq. (4), it is also partly underpinned by an assumption of a cubic supercell. For cubic cells with a high dielectric constant (*ε* → *∞*), this yields *f*^{LZ} = 0.631. By contrast, the original MP correction is recovered when *ε* → 1. As the LZ method reduces to a scaling of the original MP correction, it also leads to a *L*^{−1} dependent correction that could produce the *a*_{1} term in Eq. (2).

The interpretation of the LZ IIC is simple—the finite size of the supercell prevents ideal screening behavior being observed. Instead, the expected dielectric screening is enhanced by the electrostatic bound charge present in the supercell. This additional screening depends on the magnitude of the bound charge, which in turn depends on the dielectric constant. The LZ method requires the dielectric of the system to be constant and isotropic, presently limiting its application to bulk systems.

We note that the LZ method was also introduced with a potential alignment technique based on atomic site averages, performed over a large volume of the cell. Hence, when LZ results are presented in the literature, they may also include this alignment procedure. More complex methods for performing potential alignment fall beyond the scope of this paper, and we limit ourselves to planar averages. Additionally, we will later argue that the LZ IIC does well describe the image electrostatics, which encourages simpler alignment procedures.

Another distinct conceptual approach is the defect correction method of Freysoldt, Neugebauer, and Van de Walle (FNV).^{8} The FNV method is different to the other approaches as both the IIC and the potential alignment problem are tackled simultaneously. In this method, the charge induced by the defect is modeled as a Gaussian distribution like that shown in Fig. 2(b) (or other simple localized distribution), which normally make only a very minor change to the interaction strength *E*_{IIC}, and produces an energy almost identical to $EIICMP$. Then, the electrostatic potentials of the model system and the DFT calculation are matched far from the defect, introducing extra polarization effects into the correction.

In this method, the correction is given as

where $EIICFNV$ is the interaction energy of the defect charge model used and Δ*Ṽ*_{q/0} is the shift in potential between the charged and neutral defect, beyond that introduced by the calculation of $EIICFNV$. (Potential alignment is discussed in more detail in Sec. IV.) Equation (8) contains two terms, where $EIICFNV$ can be identified with the *a*_{1} term and $\u2212q\Delta V\u0303q/0$ can be associated with the *a*_{3} term in Eq. (2).

Out of the considered methods, the FNV approach offers the most flexibility because both the charge model used and the region over which potential alignment is calculated can be varied, effectively leading to a whole family of corrections. These options can allow some fine tuning of the magnitude of the correction to the specific defect being studied. More recent work on the method has looked at fitting these charge models to defect state orbitals directly.^{22}

The FNV method has been extended to surfaces and interfaces.^{23} However, it can be complicated in practice by the effects of atomic relaxation, as the movement of ions introduces long range changes in the potential not related to alignment effects.^{9,11} The Kumagai-Oba (KO) method^{9} has been suggested to circumvent these limitations.

One important question is whether we should expect IICs to have different magnitudes for positively and negatively charged defects. The corrections based on analytic theory (MP and LZ) are identical for positively and negatively charged defects. The FNV method is observed to break this symmetry when it is applied, a flexibility introduced by the Δ*Ṽ*_{q/0} term.

## III. PROPOSED IMAGE CHARGE CORRECTION

Above, we have seen electrostatic methods for improving the MP point charge description of defect-defect interactions. In this work, we calculate interaction corrections from DFT electronic densities directly, rather than approximating the defect charge distribution.

To do this, we follow very similar conceptual lines to the density countercharge (DCC) method of Dabo *et al.*,^{24} where the charge due to the defect is modeled by an explicit charge distribution and a full electrostatic model system is solved (in the DCC method, corrections are then applied self-consistently within the DFT calculation, in contrast to our approach where corrections are applied after the calculation). In the DCC method, the spurious interaction with images is defined as

where the interaction energy *E*_{IIC} is defined as being caused by the model charge $qr$ interacting electrostatically with its own images. This interaction potential $vIIC$ is broken down into two contributions, which can be computed more easily than evaluating $vIIC$ directly. These two components are

These two terms are the potentials produced by the charge distribution $qr$, when it is interacting with itself in the absence of other charges, $viso$, and when it is surrounded by periodic images, $vPBC$. We note that the periodic potential $vPBC$ also includes the contribution of a jellium background.

Previously, the DCC method has been developed for molecules in vacuum, where the dielectric *ε* = 1 through all space. In order to produce a *post hoc* correction for bulk materials, we need to include the effect of the dielectric of the host material. The appropriate version of the Poisson equation to use is

where the model charge distribution *q* is embedded in a host crystal with dielectric tensor *ε*, producing the potential $v$.

In the original DCC paper, the whole charge distribution is used in the corrective electrostatic calculation. But to apply similar ideas to a periodic material, we must distinguish between electronic charge which is present in the bulk system (which in turn contributes to the dielectric properties) from that induced by the defect itself. In order to do this, we need to define the model charge distribution $qr$. When ions are not allowed to relax, this can be simply defined as the electronic density difference between the neutral defect and the charged defect in state *q*. This yields a charge distribution with total magnitude *q*,

where $nq=qdefect$ is the electronic density of the supercell containing the charged defect and $nq=0defect$ is the electronic density when the defect is in a charge neutral state. If relaxation around the defect is not allowed, the atomic structure is identical in both calculations.

When the defect geometry is allowed to relax, there is more flexibility in the construction of the defect charge model, *q*. In this paper, we fully relax the atomic geometry of the charged defect in the calculation of $nq=qdefect$. Then, we use this same ionic geometry in the calculation of $nq=0defect$. This matches the ionic structure that is present in the image cells. Effectively, we replace charged defects in the image cells with neutral defects of the same ionic structure. If the geometry is also relaxed in the calculation of $nq=0defect$, large fluctuations are introduced into *q*, reflecting the movement of ions.

It must be stressed that these charge models do not reflect just the electrons added (or removed) from the system, but also contain information about the polarization of the supercell as a whole.

This charge model includes a jellium background only implicitly and hence is of magnitude *q*. When the interaction potential $vIIC$ is calculated, a jellium background is effectively introduced in exactly the same way as for DFT calculations. This is because we follow the same convention as the DFT calculation that generated the charge model^{19,20} and set the average value of the periodic part of the interaction correction model $vPBC$ to zero, in the same way as the average electrostatic potential *V* within the DFT calculation is set to zero.

Once $\epsilon r$ and $qr$ are defined, the interaction energy can be calculated through Eq. (11). In this paper, we use the open source libraries dl_mg^{25,26} and pspfft^{27} to solve these electrostatic problems directly, without the introduction of further approximations. dl_mg is an iterative multigird solver.^{28} pspfft is an FFT based solver for isolated boundary conditions.

## IV. POTENTIAL ALIGNMENT

In Secs. II B and III, we have discussed interaction corrections *E*_{IIC}. However, it is readily observed that the use of simple interaction corrections (such as MP) is not always successful by itself. Significant improvements in cell size convergence can be seen if a potential alignment correction, *q*Δ*V*_{NAP}, is also applied.^{4,11} How much of this is due to failures in interaction corrections *E*_{IIC} is not clear. This is a question of how much of *V*_{NAP} is introduced by the same image interactions treated in *E*_{IIC}. It has been suggested that no potential alignment should be required as long as *E*_{IIC} can be determined with sufficient accuracy.^{9}

The need for potential alignment is introduced into periodic DFT calculations by NAP, the convention of setting the average electrostatic potential to zero,^{18–20} $V=0$. NAP can introduce non-physical shifts into the electrostatic potential, which effectively shift the position of the bulk VBM [*μ*_{e} in Eq. (1)] and hence lead to errors in the formation energy of charge defects. The average potential of the bulk supercell and the supercell containing the charged defect should be different, so NAP results in the electrostatic potential of the bulk crystal not being recovered far from the defect. An example of this behavior is shown in Fig. 3.

Generally speaking, when using PBC DFT, total energies of neutral systems can be directly compared between different DFT calculations, but the eigenvalues can always be shifted by an unknown amount, with this unknown amount being introduced by NAP. When the system is charged, both the eigenvalues and the total energy can be affected by NAP.

If the shift in the electrostatic potentials introduced by NAP can be identified, its non-physical contribution to the formation energy can be removed. In order to determine the shift, a potential difference between a supercell containing the charged defect and a supercell of the same size containing the bulk crystal can be calculated

where the electrostatic potential of the supercell containing the charged defect, *V*_{q}, is compared against the potential of a bulk supercell of the same shape, *V*_{b}. As the average of both the electrostatic potentials *V*_{q} and *V*_{b} is zero, the average of Δ*V*_{q/b} over the entire supercell will also be zero. To determine a scalar potential shift, Δ*V*_{q/b}, between the two configurations, it must be averaged over some region of the supercell that is chosen to reproduce bulk properties,^{22}

If the “far region” where Δ*V*_{q/b} is averaged over does reproduce bulk properties in the defective cell, then the difference between the potentials must be the unphysical potential shift introduced into the calculation by NAP. Hence,

Equation (14) inevitably introduces a lot of choices in the way the far region that is averaged over is defined, potentially leading to different values of the scalar alignment shift Δ*V*_{q/b}. In this paper, we limit ourselves to one of the simplest averaging techniques, calculating the xy planar average of $\Delta Vq/br$ and then taking the z-value furthest from the defect. In this case, the far region is the xy plane furthest from the defect. Note that Eq. (15) is approximate, as the scalar potential shift introduced by NAP will only be correctly identified if the “far” region where the average is performed over does in fact reproduce bulk properties.

Possibly, some of the total potential alignment shifts present in Fig. 3 are caused by the same electrostatic interactions with image defects that were discussed in Secs. II B and III. Hence, we should be careful that the potential changes included in the calculation of *E*_{IIC} are not double counted in the potential alignment correction *q*Δ*V*_{NAP}. Avoiding this double counting is the main motivation of the FNV method.^{8} We will demonstrate that the IIC we introduced above can be used to predict part of the total potential alignment and that significant contributions to the potential alignment remain that are caused by effects entirely separate to IICs.

### A. Decomposition of potential alignment

In the standard approach, Δ*V*_{NAP} is calculated by comparing the electrostatic potentials of a supercell of the bulk material and a supercell containing the charged defect, via Eqs. (13)–(15). However, when the problem is approached in this way, it remains unclear how NAP is introducing a potential shift. In order to further explore the causes of potential alignment, we split the process of forming a charged defect into steps and compare the potentials of four DFT calculations. Such a procedure is illustrated in Fig. 4. As the start and end DFT calculations are identical to the standard approach, splitting up the alignment in this way produces a total potential alignment $\Delta Vr$ which is identical. This leads to the following expression for the total alignment:

where the total change in potential of the bulk and defective systems Δ*V*_{q/b} is composed of the three alignment terms between the four configurations introduced via Fig. 4. Evaluating the potential in this way does require two additional single-point DFT calculations that would not otherwise be carried out (configurations B and C), whereas traditional approaches use only configurations that have already been calculated to evaluate the formation energy. We note that the charge model we defined in Eq. (12) is the charge density difference between configurations C and D. We show an example of this decomposition of the total potential difference in Fig. 5, for a $VC\u22122$ defect in diamond.

### B. The Δ*V*_{AB} component

Interestingly, in our example in Fig. 5, we observe that a significant potential shift is introduced far from the defect by the removal of a carbon atom from the supercell alone, producing the Δ*V*_{AB} component. This particular alignment component is shown in Fig. 6, for the diamond supercell.

The Δ*V*_{AB} component is extremely local to the defect site, so bulk properties are quickly recovered far from the defect. This is shown by how quickly Δ*V*_{AB} takes a constant value far from the defect. This ensures that the average in Eq. (14) can be carried out easily and without ambiguity—any sensible averaging technique will return the same value. As a result, for this component of potential alignment, we expect the approximate Eq. (15) to be very accurate.

Δ*V*_{AB} can be readily approximated by simply placing an atom of the type being exchanged in a periodic box of the same size, as justified in Fig. 7. This approximation is accurate to ∼30% for the examples we considered, as shown in Table I. Each atom present in the periodic DFT supercell makes a contribution to the average potential (which is conventionally neglected and set to zero), and hence changing the number of atoms in the supercell modifies the potential shift introduced by NAP.

The electrostatic potentials depend on the exact distribution of electric charge within the atom. As a consequence, any component of a DFT calculation that modifies the electron density can also modify the potential shift associated with each atom. The atomic species is clearly the most important, but the chemical environment in which the atom is present in also makes an important contribution to the distribution of its electric charge.

We will further develop these ideas in Sec. IV F, in order to predict how Δ*V*_{AB} scales with the size of the supercell.

### C. The Δ*V*_{BC} component

The Δ*V*_{BC} component is the most challenging to interpret. This potential change is introduced when the final ionic configuration of the charged defect is introduced. Hence, this component is zero by definition when the formation energy of an unrelaxed defect is considered. An example of all three components along with the total change in potential is shown in Fig. 5. For the example of diamond, this last alignment component makes the smallest contribution to the total potential shift far from the defect. In this example, the relaxation of ions makes almost no change to the electrostatic potential far from the defect. In this case, it appears that NAP does not introduce a meaningful potential shift between configurations B and C.

In some other materials, such as MgO, $\Delta VBCr$ is large, far from the defect. However, it never tends to a constant value, making it difficult to select a far region for Eq. (14). In turn, this will damage the approximation made in Eq. (15).

When the contribution of Δ*V*_{BC} to Δ*V*_{NAP} is entirely neglected, our corrections appear to perform as reliably for both examples where Δ*V*_{BC} is large and small. This suggests that this component does not in fact introduce a significant NAP shift. We can justify this observation using the following argument—as the number of electrons and ions in the cell is kept fixed, as well as their total density, there is not a large change in the neglected average potential. We note that the Δ*V*_{BC} is large in materials where the original formulation of the FNV method is believed to perform poorly and is the alignment component most strongly altered by the use of the KO far averaging method.^{9} It would be interesting to explore in detail how the KO sampling method affects the different alignment components we define.

### D. The Δ*V*_{CD} component

Finally, we consider the potential difference Δ*V*_{CD} between the charged defect and a neutral reference structure (as defined in Fig. 4). We observe the following relationship between our corrective potential, as introduced in Eq. (10), and the Δ*V*_{CD} potential alignment component:

The electrostatic potential difference present in the DFT calculation is experienced by all the electrons in the calculation and therefore is not screened, whereas the interaction model $vPBC$ is experienced only by the additional defect charge and hence is screened by the dielectric constant *ε* (which is *ε*_{0} when the ionic structure is relaxed and *ε*_{∞} when it is held fixed). Although this factor of *ε* is justified by this electrostatic argument, it is not included in the derivation of the FNV method.^{8} We note that the factor of two in this expression is due to double counting of the electron Hartree potential in castep and may not be required in other codes depending on the electrostatic convention used (this factor is not required in cp2k). Finally, only the periodic part of the interaction model, $vPBC$, contributes. The periodic DFT calculation does not contain the isolated component of the image interaction, $viso$—as a result, Δ*V*_{CD} does not contain enough information to infer *E*_{IIC}.

An example of this behavior is shown in Fig. 8. Across the defects we considered in this paper in isotropic dielectrics, we find Eq. (17) is satisfied to within 1% (for a material with an anisotropic dielectric, *ε* should be replaced with a tensor, but we did not explore how this could be done in this paper). This correspondence demonstrates that the potential alignment component Δ*V*_{CD} is already fully incorporated in our interaction correction and need not be considered separately. Note that this is the only part of potential alignment that is influenced by the jellium background present in charged calculations.

### E. Discussion of all components

We show an example of the finite-size scaling resulting when the various components of potential alignment are used to perform alignment corrections, in Fig. 9. We switch to the example of an oxygen vacancy in MgO, as Δ*V*_{BC} is significant for this system. By contrast, the Δ*V*_{AB} component is relatively small in this case. These results demonstrate that alignment with Δ*V*_{AB} makes an improvement to the size scaling, but that the Δ*V*_{BC} and Δ*V*_{CD} components appear to do more harm than good. As we have argued, inclusion of Δ*V*_{CD} is expected to be harmful as this component is already fully included in the IIC. The negative contribution of the Δ*V*_{BC} contribution requires some further discussion.

The accurate determination of Δ*V*_{BC} is challenging for some materials, where we cannot evaluate Eq. (14) reliably. This equation requires the selection of a region of the supercell that reproduces bulk properties, but if relaxation of the defect ionic structure leads to long-range effects, we cannot identify such a region. For the example of diamond, we recover a bulk-like region (i.e., Δ*V*_{BC} tends to a constant far from the defect) and find that Δ*V*_{BC} ∼ 0. In the case of MgO, we cannot identify such a region, and the planar average we apply is certainly not sophisticated enough to accurately determine Δ*V*_{BC}. However, the results in Fig. 9 very strongly suggest that Δ*V*_{BC} ∼ 0. If this is the case, it is best not to include the Δ*V*_{BC} component in the potential alignment correction.

Taken together, all these results suggest a tempting premise. The Δ*V*_{AB} component is conveniently both the most significant source of the potential alignment error (once an image correction has been consistently applied) and the most well-behaved component. Almost any far averaging procedure used in Eq. (14) will unambiguously determine this component because the modification to the electrostatic potential is local to the defect site and does not lead to long-range changes in the material.

### F. Finite-size scaling of the Δ*V*_{AB} component

Returning to the Δ*V*_{AB} component, simple model calculations can demonstrate that the exchange of atoms in periodic DFT leads to a potential shift. Consider an isolated atom in a large cubic box. If the atom is neutral, we expect that the same energy will be calculated under isolated and periodic boundary conditions—as well separated neutral atoms do not interact. In fact, this is the behavior we observe. Using cp2k,^{29–32} the same DFT calculation can be performed with both boundary conditions.^{33,34} For a calculation on an Sn atom in a 10 Å cubic box, the total energy is only modified by 0.5 meV. However, as shown in the first panel of Fig. 10, the electrostatic potential contains a significant shift of 0.33 V. This shift is observed to be sensitive to the details of the DFT calculation, depending on the atom type, pseudopotential, and xc functional. This sensitivity to the pseudopotential used is also observed in Δ*V*_{AB}, as shown in Table II for the example of diamond. This shift is introduced by NAP.

^{a}

This PP is distributed with castep.

^{b}

Castep PP definition: $1|1.2|17|20|23|20N:21L(qc=8)$.

^{c}

castep PP definition: $2|1.4|8|9|10|20:21(qc=5)$.

This same argument that the change in boundary conditions between the two calculations introduces an alignment shift can also be made mathematically. First, consider the constraints that the boundary conditions place on the electrostatic potential of a neutral atom,

where $VatomPBCr$ is the average electrostatic potential of the isolated atom under periodic boundary conditions. For the same system under isolated boundary conditions, $Vatomisor$, we do not require an average potential of zero and instead have the average value *C*.

However, in both cases, we are describing the same physical situation and should still calculate the same energy. This requires that the only difference between the two potentials is a uniform shift,

A uniform shift in the electrostatic potential of a charge neutral system does not change its energy. Following this reasoning, we would like to understand how this potential shift introduced by PBC scales with the size of our simulation supercell. In the case of isolated boundary conditions, the vacuum level remains well defined, so this shift must be a non-physical change in the PBC potential. We can calculate this shift as

where the average is performed by taking an integral over the simulation cell of volume Ω. Note, however, that we can infer a lot about $Vatomisor$ using our knowledge of isolated atoms—it must be some kind of local well that binds electrons. As a local potential, it is only non-zero near the atom being considered. As long as our isolated box is large enough for the potential to properly decay to zero at the edges, then the actual size of the box does not affect the shape of this potential at all. This is equivalent to stating that

where the integral given in Eq. (21) evaluates to a constant *α*_{atom}, which is a property of the atom independent of the size of the simulation cell. Different choices of atoms, pseudopotentials and xc functionals will all have their own consistent value of *α*_{atom} across any large simulation supercell considered. By combining Eqs. (20)–(22), we discover the following size dependence for the non-physical potential shift:

Hence, this potential shift scales as a constant over the volume Ω. When we recall Eq. (1), this leads to an alignment correction of *q*Δ*V* and in fact would then explain the *a*_{3} term in Eq. (2). Clearly this example is simpler than the actual situation of removing an atom from a lattice and leads to an ideal volume-dependent contribution to the total energy. When the defect is present within a material, this volume scaling is only approximately maintained.

As this potential shift is associated with isolated atoms, it is worth exploring if any atomic properties can be used to predict it. We find that atomic radius is the best predictor, as illustrated in Fig. 11. We present other descriptors we considered in Sec. S2 of the supplementary material. (We use atomic radii calculated by Ghosh and Biswas,^{35} using Slater orbitals.) This suggests that this kind of potential alignment is likely to be more significant in calculations involving heavier elements. An important consequence of this argument is that we expect this component of potential alignment to swap signs when considering interstitials rather than vacancies. This behavior can be observed in the results of Shim *et al.*^{13} and Taylor and Bruneval,^{10} although these authors ascribe this behavior to different mechanisms (change in screening behavior and defect state hybridization, respectively).

## V. APPLICATIONS TO DEFECTS IN DIFFERENT MATERIALS

### A. Computational details

In order to test the introduced IIC, we studied three example vacancy defects using DFT. All these simulations were performed within the plane-wave *ab inito* density functional theory code castep.^{36} We used a modified version of release 16.1.1 that correctly calculated the non-coulombic pseudopotential energy of systems with net charge.^{20} As a result of this error, the total energy of charged systems in castep was calculated incorrectly in versions before 16.1.2. All ions were represented with norm-conserving pseudopotentials. We used the following Castep on-the-fly pseudopotential definitions: For diamond, C 2|1.4|10|12|13|20N:21N(qc = 7). For MgO, O 1|1.2|23|26|31|20N:21L(qc = 9), Mg 1|1.8|3|4|4|30N:31L:32N. For SnO_{2}, O 2|1.1|14|16|19|20N:21N(tm,nonlcc), Sn 3|2.2|6|8|8|50N:51N:42N(tm,nonlcc). Geometry optimizations were performed with the BFGS algorithm.^{37}

The dielectric constants used in this paper were calculated using castep’s implementation of the density functional perturbation theory (DFPT).^{38}

For diamond, we used high quality cutoffs, in light of difficulties encountered by previous studies of carbon vacancies in this material. In particular, previous work has demonstrated the importance of **k**-point sampling in calculating accurate formation energies in this material.^{39} Brillouin zone sampling was achieved using the Monkhorst-Pack (MP) scheme.^{40} The **k**-point grids used were chosen to maintain a constant sampling of Δ*k* = 0.0234 Å^{−1} across the sizes of the supercell used. We also used a planewave cutoff of 2000 eV. This was equivalent to a 24 × 24 × 24 sampling in the conventional unit cell of 8 atoms. Additionally, the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional was used,^{41} in order to compare directly with previous studies.^{13,39}

Using the above cutoffs, the structural properties of diamond were well reproduced. The lattice constant was calculated as 3.56 Å, 0.3% smaller than the experimental value of 3.57 Å. As is typical, the bandstructure was not as well replicated, with a bandgap of 4.26 eV. This is a 22% underestimation of the experimental value of 5.49 eV.

For MgO, a more coarse parameterization was used, as this ionic material is computationally cheaper to accurately simulate. The PBEsol exchange correlation functional was used,^{42} as it led to improved agreement with MgO’s experimental properties. The planewave cutoff was 843.55 eV, and a 7 × 7 × 7 **k**-point grid was used to sample the unit cell.

The experimental lattice constant of MgO is 4.21 Å, compared to the value of 4.24 Å calculated here.^{43} This represents a 0.6% overestimate. This leads to a bandgap of 4.5 eV, a 43% underestimation compared to the experiential gap of 7.7 eV.^{44}

For SnO_{2}, a planewave basis set with a cutoff energy of 2100 eV was used, with a 4 × 4 × 8 **k**-point sampling for the unit cell. We used the PBE functional^{41} to directly compare with a previous study. We calculated lattice constants of *a* = 4.78 Å, *b* = 3.20 Å compared to the experimental values of *a* = 4.74 Å, *b* = 3.19 Å.^{45} These were 0.98% and 0.39% overestimates, respectively.

The bandgap was calculated as 2.04 eV, which was a 43% underestimate of the experimental value of 3.59 eV.^{46}

We note that the semi-local functionals used in this study do not produce accurate formation energies,^{47} but do allow the use of larger supercells, for which size dependence can be seen more clearly. The use of hybrid functionals enables an improved description of defect level positions relative to the bulk bandstructure, but does not strongly affect the electronic charge densities which form the focus of this paper.

### B. Supercell size convergence in diamond

In diamond, the dielectric screening is isotropic and principally due to electronic contributions. In fact, ionic relaxation makes a negligible contribution to screening, and the static and high frequency dielectric constants calculated using DFPT are identical (*ε*_{0} = *ε*_{∞} = 5.76). We used a carbon chemical potential calculated from bulk diamond, as in the previous studies. Previously the formation energy of the $VC0$ defect was reported as 7.31 eV,^{13} but in good agreement, we calculated it as 7.18 eV.

First, we looked at the unrelaxed carbon vacancy, as charge corrections have been studied in detail for this defect.^{13,22} The V_{C} defect supports both positive and negative charge states. Shim *et al.*^{39} considered several correction methods and reported that the MP method appeared to perform reliably for negatively charged carbon vacancies, but seemed to fail for positively charged vacancies. They suggested that this behavior was due to different screening responses for the two types of defect. Later work by Freysoldt *et al.*^{22} attempted to capture this physics using the FNV method, with exponentially decaying Gaussians as the model charge. This motivated the present study—changes in the screening response of the material should be naturally included in our method. We present the internal polarization ratios *f* we calculated using various image corrections in Sec. S1 of the supplementary material.

Using the procedure defined in Eq. (12), charge models were constructed for the investigated defects. The positive and negative charge states of the defect distribute the charge in very similar ways, as shown by the comparison in Fig. 12. This leads to similar interaction corrections for both states and suggests that the polarization response of diamond does not strongly depend on the sign of the charge.

This result does not explain the behavior observed by the previous Shim *et al.* study^{13} and suggests that we need to find another reason to understand the difference between the correction of the positive and negative states they reported.

The supercell size convergence of two of these charge states is shown in Fig. 13. The 1 × 1 × 1 cells are excluded, as these cells provide a poor description of the defect. We avoid these small cells in all the convergence plots shown in this paper.

The first general point to make is that applying a Δ*V*_{AB} alignment correction does make a significant contribution to the total corrections calculated. However, this alignment increases the formation energy of negative states and lowers the formation energy of positive states. This is in contrast to the IIC, which always increases the defect formation energy. This is a general result for cubic cells, where the jellium background decreases the formation energy more than the repulsive interaction with images which increases it. It can be seen that the unaligned corrections (dashed lines) show a stronger size dependence. Introducing alignment makes the solid lines flatter and makes the calculated formation energy less dependent on the supercell size. But also, only including alignment without the IIC appears to make the formation energy of positive states worse. This changes the way the IIC methods compare for smaller cells (≤216 ions), when the Δ*V*_{AB} potential alignment is neglected. Note that the finite-size scaling of results with only the MP correction appears to be improved for negative defects over positive defects.

This shows how potential alignment can directly cloud comparison of IICs. The potential alignment component Δ*V*_{AB} is unrelated to the image interactions present in *E*_{IIC}. If the positive alignment correction is neglected for the negative vacancies, then an IIC that is too large appears more accurate. In the opposite way, weaker IICs appear to be more effective when the positive states do not include the negative *V*_{AB} alignment correction.

When Δ*V*_{AB} potential alignment is included in the PC corrections [yielding MP(A) and LZ(A), respectively], we can see more clearly that the MP IIC includes stronger interaction than the Poisson IIC, just as the LZ method contains a slightly weaker interaction.

We suggest that the above discussion helps understand the results reported by Shim *et al.*^{13,39} The MP method appears to yield an accurate correction for the negative charge states because the MP method’s overcorrection is compensated by the missing potential alignment term. For the positive states, this favorable error cancellation does not occur.

Finally, the FNV correction falls between the corrections of MP(A) and Poisson(A). The exact ratio is dependent on the charge state in question. There is freedom in the construction of the FNV correction, which allows other models to return different ratios. For these diamond calculations, we have used the default parameters. Diamond vacancy specific values were suggested in a previous publication,^{22} but we find that these modifications do not affect the overall trends observed.

When we allow full atomic relaxation, the results are quite similar to those shown previously. As shown in Fig. 14, the FNV method is the most strongly affected—as discussed in the literature, and it is sensitive to the changes in potential alignment caused by atomic relaxation.^{9,11} Atomic relaxation increases the volume in the cell that is perturbed by the defect and hence makes it harder to define a region far from the defect. As we use a Δ*V*_{AB} alignment correction based on the unrelaxed cells, this complication can be completely avoided in the other results we present. Additionally, the corrections calculated using the Poisson method change only very slightly. As the static and high frequency dielectric constants are identical, the MP and LZ corrections do not change on relaxation.

Some examples of the charge models used in our correction are shown in Fig. 15, which shows the minor role of relaxation in this material.

Our final defect formation energies are summarized in Table III.

### C. Supercell size convergence in MgO

In MgO, electronic and ionic contributions to the dielectric response are both important and hence allowing relaxation triples the dielectric constant (*ε*_{∞} = 3.19, *ε*_{0} = 9.41). For the exchange of oxygen atoms, we used a chemical potential set to an isolated O atom in its triplet ground state.

We looked at the oxygen vacancy defect in MgO, which supports a + and +2 charge state. The supercell size convergence of the relaxed defect is shown in Fig. 17. These results show very similar trends to those we presented for diamond.

Previous work found that the formation energy of the $VO0$ defect with respect to the chemical potential of an isolated oxygen atom is 10.08 eV.^{48} In good agreement, we calculated a formation energy of 10.27 eV, using this chemical potential (the majority of this energy difference is due to our use of the PBEsol functional—using PW91, we calculate a formation energy of 10.05 eV).

The calculated defect charge models are shown in Fig. 16. These plots show clear fringes of positive and negative polarization in our charge model, as expected from an ionic crystal.

Corrected formation energies for the relaxed MgO cells are shown in Fig. 17. The FNV corrections also become larger than the MP corrections when relaxation is allowed. The difficulty in applying the FNV method to materials that display large atomic relaxations has already been reported, as the method can become very sensitive to the region in which potential alignment is performed over.^{9}

Our final defect formation energies for MgO are summarized in Table IV.

### D. Supercell size convergence in SnO_{2}

In SnO_{2}, the ionic contribution to screening is stronger than the electronic contribution. Additionally, the unit cell is a rectangular cuboid and the dielectric constants are anisotropic. The calculated dielectric responses were

where *ε*_{∞} is almost isotropic and *ε*_{0} is anisotropic. We used a chemical potential for Sn calculated from an isolated atom in its triplet ground state.

The defect formation energy of the relaxed $VSn0$ defect has been reported as 16.0 eV.^{49} This is higher than the 15.56 eV we calculated. Additionally, we find that a net magnetic moment of 4*μ*_{B} on the neutral defect is energetically favorable, in agreement with Ref. 50 and disagreement with the value of 3.63*μ*_{B} previously reported.^{49} A more detailed study of the V_{Sn} defect will follow in a future publication.

The new method we have outlined can still be applied to materials with an anisotropic dielectric. Additionally, the MP method has already been extended to materials with anisotropic dielectrics.^{21} However, the LZ method has not been extended. Instead, we apply *f*_{LZ} (calculated using the knowledge of the non-cubic cell, but not the anisotropic dielectric) to the anisotropic MP correction. We have not applied the FNV method, although Kumagai and Oba’s^{9} work generalizes the FNV method to anisotropic media.

The corrections calculated in Fig. 18 are similar to those for the other materials. The anisotropic dielectric further increases the screening factor *f*, calculated using our method, by approximately 0.008. Hence, the anisotropic dielectric only has a minor effect on defect formation energies in this material. Potentially, this could be because the anisotropy in the cell vectors and dielectric is balanced against each other. A stronger effect could be seen in materials where this is not the case.

## VI. DISCUSSION AND CONCLUSIONS

In this paper, we have developed an IIC which directly calculates the electrostatic image interaction introduced by the additional charge associated with a charged defect. When applying this method to cubic bulk crystals, we found strong agreement with the Lany-Zunger^{6,7} IIC. We also demonstrated the extension of the developed method to SnO_{2}, as an example of a material with an anisotropic dielectric constant, beyond the scope of the LZ method.

However, even with this new correction, we find that an appropriately constructed IIC still requires a potential alignment procedure. We introduced a decomposition of potential alignment into three components that can, in principle, contribute to the failure of the bulk potential to be recovered far from the defect. We find that the first component Δ*V*_{AB} is due to a change in the atom number. This can introduce significant energy errors when the atoms added (or removed) have a large atomic radius. Second, the Δ*V*_{BC} component is introduced by ionic relaxation. In cases where we can still find a good reference region that reproduces bulk properties inside the defect cell, this potential shift is near zero. For examples where the relaxation around the charged defect leads to atomic displacements across the entire supercell, it remains a challenge to evaluate this potential alignment term. However, we present indirect evidence that it remains close to zero. Finally, we find that the Δ*V*_{CD} component also introduces a large alignment shift, but that this is already included in the image charge interaction *E*_{IIC}.

Although the arguments look rather different, the conclusions we have reached readily support what is normally termed the traditional potential alignment technique,^{11} where alignment is calculated between the bulk material and a neutral defect. The atomic relaxation present in the traditional approach can lead to differences to the Δ*V*_{AB} component, but these modifications are very minor for many examples (see Sec. S3 of the supplementary material). This correspondence could be expected to fail when ionic relaxations around the neutral defect are long range.

We have presented a plausible argument that Δ*V*_{BC} alignment shifts are small, but we have not been able to prove this. If a defect can be identified where this alignment term can both be unambiguously determined and is significant, this would provide an example where the methods endorsed by this paper would not be reliable. In such an example, the FNV method could be anticipated to be more effective—as the inclusion of a Δ*V*_{BC} could become beneficial rather than harmful. The application of the potential alignment decomposition introduced in this paper to a wide range of defects and materials would help to further understand the role of the Δ*V*_{BC} component.

The above results lead us to the conclusion that the Δ*V*_{AB} component we define, caused by the removal (or addition) of atoms from the supercell alone, is the dominant remaining finite-size error once the *E*_{IIC} correction has been applied. As a consequence, we expect that potential alignment corrections are not requited for calculations on self-trapped polarons. An advantage of this conclusion is that determination of the potential alignment correction is no longer hindered by relaxation of the defect geometry.

It is already known that IICs are positive for cubic cells, as the presence of a jellium background over stabilizes defects.^{12} However, in this work, we show that the sign of the potential alignment correction introduced by the exchange of atoms (the dominant *V*_{AB} component) can also be predicted, as in Table V. Generally, the total finite-size corrections required are smaller for positive vacancies and negative interstitials, as the IICs and potential alignment corrections compete and partly cancel—in the opposite case, negative vacancies and positive interstitials require larger corrections as the two sources of finite size error add together and reinforce the size dependence. This behavior has already been noted in the literature,^{10,13} but our results provide a strong justification for this observed behavior. (We note that the FNV correction does not always obey these rules for the sign of the potential alignment, as the Δ*V*_{AB} component is not always the dominant component of the FNV alignment scheme.)

. | Positive charge state . | Negative charge state . |
---|---|---|

Vacancy defect | − | + |

Interstitial defect | + | − |

. | Positive charge state . | Negative charge state . |
---|---|---|

Vacancy defect | − | + |

Interstitial defect | + | − |

All of the above conclusions combine to justify the origin of the empirical scaling rule suggested in Ref. 15, which predicts finite size dependence that scales inversely with supercell length (*L*^{−1}) and supercell volume (*L*^{−3}). The *L*^{−1} scaling is found to be introduced for cubic supercells by the defect electrostatics of a charge interacting with its images in a jellium background (the IIC). This behavior has already been predicted in Refs. 6 and 7, but we maintain the same result without using an approximate multipole expansion. We find that the second volume dependent *L*^{−3} term is caused by a potential alignment effect (the Δ*V*_{AB} component), due to the removal (or addition) of ions effectively modifying the average potential in a periodic DFT calculation. The Δ*V*_{AB} is found to be approximated by calculations on neutral atoms, which are simple enough for an analytic *L*^{−3} dependence to be derived. This term is predicted to be larger when atoms with larger atomic radii are exchanged to form a defect.

Taken together, these developments enable robust corrections across various materials and defect charge states. For cubic systems with isotropic dielectrics, the LZ IIC combined with a potential alignment between the bulk material (the A configuration) and with atoms added (or removed, the B configuration) can provide reliable corrections. Our conclusions on potential alignment are fully transferable to anisotropic dielectrics and non-cubic cells, where the *V*_{AB} component will always produce an inverse volume scaling finite size error. For the case of IICs, we present a method that can be applied to anisotropic dielectrics and non-cubic cells, where an inverse-length scaling finite-size dependence can no longer be expected.

Finally we note that the IIC introduced in this paper can be applied to more complex dielectric environments, such as systems containing multiple interfaces between materials, and such systems will be considered in a future paper. It would also be possible to use a similar solver to apply self-consistent corrections to DFT calculations, as in the DCC method.^{24} The relationship discovered in Eq. (17) is an important ingredient of such a scheme.

## SUPPLEMENTARY MATERIAL

See supplementary material for the following additional data: Internal polarization ratios *f* calculated using Eq. (6), the various correction methods and materials, and additional fits to the atomic descriptors we considered for the alignment constants of neutral atoms. We introduce the agreement between traditional alignment constants and the Δ*V*_{AB} alignment.

## ACKNOWLEDGMENTS

We acknowledge funding from EPSRC via a doctoral training partnership (No. EP/M506448). Via our membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC (No. EP/L000202), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). The authors acknowledge the use of the UCL Grace High Performance Computing Facility (Grace@UCL) and associated support services in the completion of this work. We thank Leeor Kronik, Philip Hasnip, Stephan Lany, Gilberto Teobaldi, and Keith McKenna for helpful discussions and Lucian Anton for assistance with the dl_mg library. The matplotlib library was used to produce figures.^{51}