{'fraction_non_alphanumeric': 0.0230859360697155, 'fraction_numerical': 0.010490278642305298, 'mean_word_length': 4.530322972562519, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite size scaling as a cure for supercell approximation errors in calculations of neutral native defects in InP

The relaxed and unrelaxed formation energies of neutral antisites and interstitial defects in InP are calculated using ab initio density functional theory and simple cubic supercells of up to 512 atoms. The finite size errors in the formation energies of all the neutral defects arising from the supercell approximation are examined and corrected for using finite size scaling methods, which are shown to be a very promising approach to the problem. Elastic errors scale linearly, whilst the errors arising from charge multipole interactions between the defect and its images in the periodic boundary conditions have a linear plus a higher order term, for which a cubic provides the best fit. These latter errors are shown to be significant even for neutral defects. Instances are also presented where even the 512 atom supercell is not sufficiently converged. Instead, physically relevant results can be obtained only by finite size scaling the results of calculations in several supercells, up to and including the 512 atom cell and in extreme cases possibly even including the 1000 atom supercell.


I. INTRODUCTION
Over the past decade or so first principles density functional theory (DFT) 1 has become a powerful tool for studying the properties of defects in semiconductors. It is possible to calculate formation energies, binding energies and migration barriers, to predict local structure and, up to a certain point, defect levels and electrical activity. Problems and limitations remain, however, and one of the greatest is the very limited size of the systems for which calculations are feasible: 10s or 100s of atoms, even when we wish to describe physical problems involving 1,000s or 10,000s. This leaves the results heavily influenced by errors arising from the boundary conditions. These errors must therefore be carefully studied so that their effects can be understood and accounted for when results are interpreted. There are two types of boundary conditions commonly used: open and periodic. Open boundary conditions are usually encountered in cluster calculations. The surface atoms are "terminated" with hydrogen to use up spare electrons, but are otherwise surrounded by empty space. Periodic boundary conditions (PBCs), meanwhile, are found in supercell calculations, in which a block of atoms is surrounded not by empty space but by an infinite array of copies of itself. Both approaches have strengths and weaknesses. We present here a detailed study of the problems arising from the use of the supercell approximation and we propose a method to overcome them.
Finite size errors in supercell calculations come from two main sources. Elastic errors often arise because the supercell is simply not large enough to contain all of the local relaxation around the defect, leaving calculated formation energies too high. In addition, the defect interacts with an infinite array of spurious images of itself seen in the PBCs, via both elastic and electrostatic interactions. The direct elastic interactions can easily be truncated by freezing all atoms lying on the surface of the cell at their ideal lattice positions. The electrostatic interactions, on the other hand, cannot be truncated or removed. They result in errors in the calculated formation, binding and migration energies, errors which can be on the same order as the energies themselves. For practical supercell sizes they need not even be negligible for neutral defects, since dipolar and quadrupolar interactions can remain significant. These latter can even result in errors in the calculated structures, since they favour certain symmetries and local relaxation modes over others. Hence indirect elastic errors cannot be avoided either. Finally, a third source of finite size errors is also present: the defect state wavefunctions can overlap with their images in the PBCs leading to a spurious dispersion of the defect levels which in turn can affect the formation energies, especially if the defect level is only partially filled. The errors related to this dispersion (or tunnelling) are expected to have only a fairly small and rather short ranged (exponential) effect.
Recently, various correction schemes have been suggested 2,3 to compensate for at least the leading terms in the errors arising from the electrostatic interactions. They are usually based upon fits to quasi-classical models and/or multipole expansions of the electrostatic interactions. They have met with varying levels of success but, so far at least, are generally considered insufficiently reliable. There are more direct approaches, however. Probert and Payne 4 recently presented a detailed ab initio study of the neutral vacancy in Si, considering all aspects of convergence, from basis set and k-point sampling to size and symmetry of supercells. They demonstrated that the use of "large" supercells (200+ atoms) can be essential for obtaining the correct physical results. Meanwhile, we recently presented 5 a study of the neutral vacancies in InP. We demonstrated the advantages of not only using large supercells but also finite size scaling the results obtained. We evaluated both the relaxed and unrelaxed formation energies of the phosphorus (V P ) and indium (V In ) vacancies in simple cubic supercells of 8, 64, 216 and 512 atoms. We then showed that the variation in the formation energy with supercell size, L, follows rather closely the form where E C d (L) is the formation energy in supercell "C" and a 1 and a 3 are fitting parameters. E ∞ d is the finite size scaled formation energy corresponding to an infinitely large supercell. Equation 1 has, in fact, the same form (L − 1 plus L − 3 ) as the corrections proposed by Makov and Payne 2 . We will return to the issue of whether or not this is the correct form in section IV. For the vacancies we showed that the error bars on the values obtained for E ∞ d are usually rather small: O(0.1 eV) or less, depending also upon the level of k-point convergence in individual cells. Care must be taken, however, when scaling relaxed formation energies of strongly Jahn-Teller 6 active defects, such as V P . In such cases rather wider error bars are obtained if the symmetry of the relaxed structures varies with supercell size. Numerically, we also found that (for example) E ∞ d for the unrelaxed V In is actually ∼0.2 eV above the value obtained in the 512 atom supercell, suggesting that there are cases for which even the 512 atom cell is not large enough to be converged, so that scaling becomes essential.
In the present paper we extend the study to that of the other neutral native defects: the antisites and interstitials. It should be noted, however, that our primary purpose is not the study of the defects themselves but of the finite size errors which arise when calculating their formation energies. InP has the zinc-blend structure, with two antisites, P In and In P , both with tetragonal (T d ) symmetry when unrelaxed. For interstitials there are three high symmetry sites: two tetragonal, with four In or four P nearest neighbours (X i(In) or X i(P) ) and a quasihexagonal site (X i(hex) ) with six nearest neighbours (three P and three In) and C 3v symmetry. Previous work for the isolated neutral cases has mostly been limited to the 64 atom supercell. Here, formation energies were around 5-6 eV for the tetrahedral interstitials 7 and around 3 eV for the antisites 8 . In P displayed some Jahn-Teller behaviour whilst P In did not. We will examine how these results change with larger supercells and with scaling. In the next section we describe the method to be used. In section III we will describe the basic scaling results for the various neutral defects and will examine in sections IV and V the form of the scaling and when it does and does not work. In section VI we consider briefly other charge states of the defects. In Section VII we discuss the origin of the surprising linear scaling we find in certain cases. In section VIII we estimate the size of the other non-finite size dependent errors for the purpose of comparison, before concluding in IX.

II. METHOD AND K-POINT CONVERGENCE.
We use planewave ab initio DFT 1 within the local density approximation (LDA) together with ultrasoft pseudopotentials 9 (US-PP) using the VASP code 10 . We recently presented 11  memory limitations in even the largest parallel computer facilities available to us mean that we must treat the 4d electrons of In as core, even though they are comparatively shallow. (A calculation for bulk InP in its FCC primitive cell places the In 4d states about 14.5 eV below the valence band edge at the Γ point when they are treated as valence -fairly deep but still close enough to potentially make some contributions to formation energies.) It would be preferable to treat them as valence, but then properly k-point converged calculations in the 512 atom cell -which our analysis requires -would be impossible. The size of the resulting error will be examined in section VIII.
The optimized LDA lattice constant using these pseudopotentials is 5.827Å and the band gap is 0.667 eV, compared to 5.869Å and 1.344 eV in experiment. As stated above, we will only use simple cubic supercells of 8, 64, 216 and 512 atoms, since it is important to keep to a single supercell symmetry since the scaling is different for different symmetries. No restrictions are placed upon the symmetry of relaxations, but we do not allow atoms located on the surface of the cell to relax. The exception is interstitials at the quasi-hexagonal sites in the 8 atom cell. Here, three of the nearest neighbours lie on the surface. Initial tests showed that the relaxation was not stable if all of these were allowed to relax at once so relaxations were done in stages: firstly one group of neighbours was relaxed whilst the others were kept fixed, then vice verse.
The key quantity is the formation energy where E C T (defect) and E C T (bulk) are the total energy of the supercell with and without the defect, calculated using the same values of planewave cutoff, k-point grid, etc, to make use of the cancellation of errors. The defect is formed by adding/removing n i atoms of chemical potential µ i . We use µ P = 3.485 eV and µ In = 6.243 eV, 11 corresponding to stoiciometric conditions. A planewave cutoff energy of 200 eV and a Monkhorst-Pack 4x4x4 k-point grid 13 was found sufficient 11 for converged non-relaxed calculations in the 64 atom supercell with errors around O(0.01 eV). Hence grids of 8x8x8 in the 8 atom cell and 2x2x2 in the larger cells should suffice. However, in this study we examine specifically the errors arising from the supercell approximation itself. K-point convergence is different in different supercells and we do not wish to include any significant errors due to this. Hence we need to ensure even higher levels of convergence during this studymuch higher than is normally required or practical. (Examining the scaled properties of the native InP defects themselves, over all relevant charge states, is left to future work, however, as it will be done with rather different levels of k-point convergence, pseudopotentials etc, once we have in the present work established the behaviour of the finite size errors.) We use k-point grids of up to and including 12x12x12, 8x8x8, 4x4x4 and 4x4x4 in the 8, 64 216 and 512 atom cells respectively. (Exceptions are: P In -only 2x2x2 needed in the 512 atom cell, and In i(In) and In i(P) -6x6x6 used in the 216 atom cell to check convergence.) To improve convergence further we use weighted averages over E C d values calculated using a series of grids.
helps dramatically, as it effectively increases the k-point density: the points in a 4x4x4 grid are not contained in a 6x6x6 grid, for example. There are two obvious choices for w N : the best is w N = N 3 , the number of k-points in the full Brillouin zone, but setting w N equal to the number in the irreducible wedge is not bad. (Solid and dashed lines respectively.) All subsequent results will be weighted averages using w N = N 3 . (Incidentally, for the unrelaxed neutral vacancies we find 5 errors of 0.36 eV and 0.06 eV respectively for V P and V In in the 512 atom cell when comparing Γ point only calculations to converged values, so the Γ point is never sufficient.) In reference 11 we also showed that the relaxation energy (where E C:Rx The relaxation energies used will be weighted averages using 6x6x6 and 8x8x8 Monkhorst-Pack k-point grids in the 8 . This amounts to assuming that the distortion in the bandstructure due to the presence of the defect is either localized (thus important only very near Γ) or symmetric. It introduces a small error whose significance disappears in the large supercell limit.

III. SCALING RESULTS.
In Figs 2, 3 and 4 we show the formation energies for the antisites, phosphorus interstitials and indium interstitials respectively, both relaxed (minimum energy configurations) and unrelaxed (atoms at ideal bulk lattice sites) plotted against inverse supercell size. The solid lines are fits of the four points to equation 1. The y-axis intersect of each of these fits is the scaled formation energy E ∞ d corresponding to the formation energy of a single isolated defect in an otherwise perfect, infinite lattice. The inclusion of formation energies from the 8 atom supercell could be questioned, since in itself it is so far from being converged. However, the results shown in the figures clearly support our expectation that the correct form for the scaling has (at least) three parameters. We therefore need results from at least four supercells of the same symmetry. It would have been preferable to use the 1000 atom simple cubic cell, but sufficiently converged ab initio calculations for InP defects in this cell are not possible with current facilities. In addition, the results themselves justify the use of the 8 atom cell: the scaling mostly works very well, producing small error bars and with E 008 d lying on or near the curves. This, incidentally, also tells us that the k-point convergence in the individual cells was sufficient. The cases where the scal-ing does not work so well turn out to be due to other problems, see section V.
To get an idea of the accuracy of the fitting and of the derived E ∞ d values (and the individual E C d ) four more fits (dotted lines) are added in each case. For each of these, one of the four data points has been omitted. The spread in the resulting y-axis intersects gives an error bar for the scaled formation energy E ∞ d . (It should be emphasized that the dotted lines -three parameters to three points -are not intended to be meaningful in themselves but are merely indications of the scale of the error in the real solid line fits.) The scaled formation energies are listed in table I, together with those of the vacancies taken from Ref [5].
The scaling curves also have a more general and rather practical meaning: they are essentially predictions of the formation energies in all simple cubic supercells from 8 atoms to infinity. For example, table II shows the predicted formation energies in the 1000 and 8000 atom supercells, as they would arise from k-point and basis set converged LDA calculations. (The 64,000 atom cell would be just as easy, but these predictions are more likely to be tested and -hopefully -confirmed within our lifetimes!) As a more immediate test, table II also shows the formation energies in the 512 atom cell predicted by scaling the results from only the 8, 64 and 216 atom cells: following along that dotted line which doesn't pass exactly through the 512 atom value for each case in the figures. The error in this prediction (as compared to the actual calculated values) is also shown. The errors are pleasantly small, especially considering only 3 cells have been used, including the 8 atom one. The errors in the other predictions are expected to be much smaller still.
All the relaxed structures turn out to be symmetric, with just breathing mode relaxations. The exception is In P , which shows some moderate Jahn-Teller behaviour (see section V A. ) All the interstitials apart from P i(P) relax outwards, as does In P , whilst P In relaxes inwards. This is all as expected, since P is smaller than In. Fig 5 shows the scaling of the percentage volume change upon relaxation. (The volume shown is that of the polyhedron defined by the nearest neighbours.) The fits are described in section IV and the scaled results are in table I. The error bars are derived in the same way as for the formation energies, although the dotted lines are omitted here for clarity.
The most stable neutral native defect in stoiciometric InP turns out to be the phosphorus antisite P In closely followed by the vacancy V P and then the indium antisite In P . Then come the interstitials, of which indium is the more stable. The least stable neutral defect is V In . The most stable site for the phosphorus interstitial is the quasi-hexagonal site P i(HEX) . For indium, the two tetrahedral sites are degenerate to within the error bars but numerically the phosphorus surrounded site In i(P) is 0.05 eV lower -which is probably correct, since the error bars are asymmetric. (See Fig 4.) This is a result which is only apparent in very large supercells and the degeneracy only appears at the 512 atom supercell. In smaller cells In i(In) seems more stable. What is more, scaling shows that simply taking the 512 atom result would give a relaxed formation energy about 0.4 eV too high (0.8 eV too high if we stopped at 64 atoms) which can be large enough to make real differences in the predicted physics of InP. As with the formation energy 5 of V In this emphasizes the value of finite size scaling the results of supercell calculations, since the largest cells for which we can actually do calculations can still be too small to be fully converged.
There are complications with the stable interstitial sites as a function of supercell size. Whilst both interstitials are stable at all three locations in the two smallest cells they are not so in the larger cells. For indium the quasi-hexagonal site lies about ∼ 1 1 2 eV above the tetrahedral sites and is not stable in the larger cells. An indium atom placed here migrates to a tetrahedral site upon relaxation, indicating that the quasihexagonal site is probably not even metastable for indium in the real material. For the phosphorus interstitial, on the other hand, it is the tetrahedral sites which are unstable. In this case, however, we can still obtain (hypothetical) relaxed formation energies at the tetrahedral sites by only allowing T d breathing mode relaxations. The tetrahedral sites are 1.2 eV higher than the quasi-hexagonal site. The reason for this difference in stabilities may be (partially) stearic: the unrelaxed nearest neighbour distances are shorter at the hexagonal site, where the smaller phosphorus interstitials are stable and longer at the tetrahedral site where the larger indium interstitial sits.

A. Formation and relaxation energies.
So far it has been assumed that the correct functional form for scaling is that of equation 1. This need not be the case, however. Equation 1 is based upon approximations and predictions for the form of the errors arising from electrostatic charge multipole interactions between the defect and it's im-  ages in the PBCs and these could be incorrect. In addition, we include here relaxations within finite sized supercells, so we have both elastic errors and, potentially, cross terms between the elastic and electrostatic errors, so other possible scalings should be considered. There is clearly a linear term present, plus at least one higher order term, so we consider scaling of the form: with n = 2, 3 and 4. The fit quality is assessed using the "χ 2 " test. This is not easy to do reliably with only four points, so we average over different defects. Simply summing or averaging  the χ 2 values would bias the conclusion towards the worst data sets. Instead, for each data set we find χ 2 for each value of n, select the n giving the best fit (χ 2 best ) and then calculate a quality factor Q n =χ 2 n /χ 2 best for each n. We then compare the averaged Q n .
We first examine the scaling of the elastic errors. To do this we have performed a series of relaxations in the 216 atom cell for each of six defects. In these relaxations, the number of shells of atoms permitted to relax around the defect is varied from 1→4 for interstitials or 1→5 for antisites and vacancies. (4 and 5 are the maximum numbers of complete shells which fit inside the supercell.) Since the cell size is kept constant the electrostatic interactions will be (almost) constant so any variation in E d is due to elastic effects. In the left panel of almost purely linear (solid lines in the figure) even -within the bounds of error -for the Jahn-Teller active defects like In P and V P . Adding a higher order term, such as L −3 (dashed lines) clearly only improves the fit very slightly: χ 2 is reduced by about 30% on average. Having established that the scaling of the relaxation energy, and hence of the elastic errors in E ∞ d , is linear we turn to the scaling of the formation energies with supercell size. Unfortunately there is too much scatter in several of the curves so that whilst n = 3 is best in some cases, in others it comes an often close second to n = 2 or n = 4. However, since the elastic errors are now known to be linear in relaxation radius and hence in supercell size, we can assume that it will simply add to the linear term in equation 6, so that we can calculate Q n using all 19 data sets (Figs 2 to 4 and Ref. [5]). The result obtained is that the scaling does indeed fit best with n = 3, with both Q 2 and Q 4 being 4 times larger than Q 3 . In fact, if we calculate the {Q n } using only the unrelaxed formation energies we still find n = 3 fits best and if we use the relaxed formation energies minus that of P In we find the same result. Unfortunately if we use all of the relaxed energies including P In we get a different result: n = 2 provides a better fit. This suggests the (faint) possibility of a cubic scaling for unrelaxed formation energies but quadratic for relaxed energies, which shifts the predicted E ∞ d by 0.01→0.2 eV depending on the defect. A quadratic scaling would seem odd, since it comes from neither the elastic nor the electrostatic errors but it could possibly arise from the cross terms between them. An alternative and perhaps more likely explanation is that the uncertainty regarding n is a result of the spurious dispersion of the defect levels mentioned above. This adds an exponential term to the formation energies in the smallest supercells, blurring the picture slightly. It seems most likely that the scaling should be cubic even for the relaxed formation energies: even now it is only one defect which seems to particularly disagree. This should be confirmed once reliable calculations involving the 1000 atom supercell become feasible.  It is very clear from table I that some formation energies scale better than others. When considering the scaling for the neutral vacancies 5 we noted that the scaling becomes more difficult to do reliably for strongly Jahn-Teller active defects such as V P . For a Jahn-Teller active defect there is a partially filled degenerate state at the Fermi level, which the Jahn-Teller theorem 6 says will be lifted by symmetry reducing relaxations (if no other effect achieves this first). This leads to poor scaling since the symmetry of the most stable relaxed structure can vary with supercell size, so that data points from some cells scale differently to those from others. In order to get good error bars for scaled formation energies each possible reduced symmetry structure must be scaled separately. Amongst the current defects there is a further example of a Jahn-Teller active defect, In P . The distortions here are much weaker, so the error bar is ±0.06 eV even when symmetry differences are ignored, which is still reasonable. Never-the-less, we have done a search for the various stable and metastable structures in the four supercells. Their various formation energies are shown scaled in Fig 7, together with a) the scaling of the lowest lying formation energy irrespective of symmetry (labelled Min) and b) the formation energy when only T d symmetry breathing mode relaxations are permitted (labelled Td).
In the 8 atom cell the lowest lying structure has the full T d symmetry of the unrelaxed antisite. Indeed, this is the only stable structure we find in this cell. In the larger cells relaxation breaks the T d symmetry to give (primarily) C 3v , D 2h or C 2v point groups at the defect site. C 3v symmetry is reached by the antisite atom moving either towards the midpoint of three of its nearest neighbours (away from the fourth) in a DX like structure (DX in the figure). or towards one neighbour and away from the other three -an anti-DX (ADX) structure. C 2v arises when the antisite moves towards one pair of neighbours and away from the other pair, i.e. towards a bond centre site (BCR). The degeneracy can also be lifted when the antisite stays still with the four antisite-neighbour distances equal, but the angles between the bonds changes. D 2h structures occur if the neighbours rotate to form two opposing pairs of either shorter or longer neighbour-neighbour distances (DDM -double dimer or DAD -double antidimer, respectively). In the 64 and 512 atom cells the lowest lying structure is a DAD structure with two neighbour-neighbour distances respectively 5% and 8% shorter than the others. The most stable structure in the 216 cell is a 7% DAD-like structures but with an additional 4% DX-like distortion, although a 4% BCR structure and a 7% pure DAD structure (with no DX component) both come a close second. (The distortion quoted for the BCR, DX and ADX structures is the % variation in antisite-neighbour distances.) In the 512 atom cell the potential energy surface for small (up to ∼2%) DX distortions from the DAD structure is also very flat. Overall, these results suggest that a lone In P in an infinite supercell would have a DAD structure with a formation energy lying about 0.4 eV below the T d structure found when only breathing mode relaxations are allowed, and 0.1 eV below the formation energy found by scaling the minimum formation energy irrespective of Jahn-Teller structure.
The changes in relative stability of the different structures are due to one or a combination of two things: a) stabilizing/destabilizing dipolar, quadrupolar or higher interactions, which can in certain cases lift the symmetry without distortion, (in the 8 atom cell for example,) or favour certain Jahn- Teller structures over others. These effects become weaker as the cells grow. b) The lack of shells of atoms in the smaller cells to absorb the elastic strain, which favours more symmetric structures. In the right hand panel of Fig 6 the variation in the degree of Jahn-Teller distortion for the ADX and DAD structures was plotted versus the number of shells permitted to relax within the 216 atom cell. For the ADX structure there is virtually no variation at all and the same was previously 5 found for V P . For the DAD structure, on the other hand, a rather strong variation is found. This suggests that elastic effects are involved for some local distortion symmetries but not for others, thus further complicating the possible variations in lowest symmetry structure with supercell size. At least for the ADX structure the charge multipolar interactions act essentially in competition with the normal Jahn-Teller mechanism, making it uncertain if the correct structure has been found at all for smaller supercells. It is sometimes pointed out (Ref [14] and elsewhere) that one way around this problem would be to use k-point integration at the Γ point only, since this restores the degeneracy of the degenerate levels (prior to relaxation). However, since for InP (and doubtless many other materials) the Γ point does not give sufficiently converged formation energies even in the 512 atom supercell this will simply result in unconverged results: errors arising from the use of just the Γ point can be tenths of eV, often the same size or larger than the splittings between various stable and metastable Jahn-Teller distorted structures. Here we study the relative stability of all the possible Jahn-Teller distorted structures with converged k-point grids and then scale them to predict which structure will be most stable in the infinite supercell limit, where these spurious degeneracy lifting interactions become zero.

B. Defect level dispersion and defect states outside the bandgap.
The scaling is also rather bad for the interstitials, even though no Jahn-Teller behaviour is anticipated or detected what-so-ever. One might have expected that the reason for this was the defect level dispersion which was also identified in section IV as a possible source of the slight uncertainty in the correct form for the scaling equation. The dispersion can lead to errors in the formation energies of individual defects and since it should decrease exponentially with cell size it could affect the error bars on the scaled infinite supercell formation energies. In a fully occupied defect level the dispersion should have no direct effect as long as the defect level lies within the gap at all k-points, since the average energy of the band should equal the energy of the level in the absence of dispersion. On the other hand, the dispersion can lead to hybridization of the defect level with conduction band states of the same symmetry, thus artificially lowering the mean value of the defect level and hence of the defect formation energies in the smaller cells. (The same holds for empty defect levels hybridizing with valence band states, allowing the latter to be lowered in energy.) For partially filled defect levels the effect upon the formation energy is more direct, since only the lower parts of the defect level (band) will be filled, again leading to too low a value for the formation energy. However, the connection between the amount of defect level dispersion and the size of the errors in the formation energy in a particular supercell is not simple. Indeed, checks on the bandwidth of the defect levels find no correlation at all between them and the scaling error bars: for example the unrelaxed structures for both P In and In i(In) have a defect level in the lower part of the gap with a rather large dispersion (∼0.6 eV in the 64 atom supercell) and yet the scaling error bar for P In is 0.02 eV whilst that for In i(In) is 0.35 eV.
The main difference is actually that In i(In) also has a partially filled defect level inside the conduction band, resulting in an electron occupying a delocalized state at the conduction band edge. This occurs for all of the defects which have poor scaling error bars. This is seen clearly in Fig 8 where we present the level diagrams for all the neutral native defects, shown at the Γ point in the 64 and 512 atom cells. While the antisites have no delocalized states, the interstitials always have at least one electron occupying a state at the conduction band edge. For phosphorus at the tetragonal sites not one but three electrons lie above the conduction band edge, which may be why they are not stable at all in the larger supercells. The level diagrams for the vacancies seem at first glance to contradict the rule: V P has one delocalized electron in the smaller cells and V In has three delocalized holes in all four cells, despite the fact that the scaling is rather good for these defects. This is itself a finite size effect, however. For V P a triply degenerate defect level lies just above the conduction band edge. This level Jahn-Teller splits, with one state moving downwards. In the smaller cells it does not make it to the band gap so the metastable distorted structures found must arise from hybridization between the localized Jahn-Teller split levels and levels at the conduction band edge. For the largest cells, however, the lower Jahn-Teller split level drops into the gap, so that the charge state becomes stable. In the case of V In we again see a triply degenerate localized defect level, this time inside the valence band leaving three holes at the valence band edge, and even in the 512 atom cell this is all we see. However, as the supercell size grows the defect level moves rapidly towards the valence band edge. Clearly, for an isolated defect in a real material we would expect this level to lie inside the band gap, even though the 512 atom cell is not sufficiently large to show it. This again underlines the need for using comparison and scaling rather than just large supercells. (Similar effects are visible for P i(In) and P i(P) , although they remain unstable.)

VI. POSSIBLE CHARGE STATES ACROSS THE BAND GAP.
Although we only consider the neutral defects in this paper, the level diagrams in Fig 8 also give a rough indication of which charge states are most likely in different parts of the band gap. For the indium interstitials the +1 charge state seems the most stable in the upper part of the band gap, but two transition levels will lie in the mid-gap, giving a +3 charge state for p-type material. For P i(HEX) we expect a +1 charge state across most of the gap. With increasing supercell size a second defect level approaches the band edge from below so there may be one or two transition levels near the bottom of the gap. This is a results which is only apparent when one compares or scales the results from different supercell sizes: even in the 216 or 512 atom cells it is not apparent. P In has a single, filled level in the upper half of the gap, suggesting 0 for n-type material and +1 or +2 otherwise. In P has a threefold degenerate level containing 4 electrons in the middle of the gap. A large number of transition levels across the entire gap are thus expected, ranging from -1 or -2 in n-type material to potentially even +4 for strongly p-type material. V P should have two transition levels in the upper half of the band gap, +1 and -1 being the most stable charge states in strongly p and n-type material respectively. For V In we expect the -3 charge state to be stable from the mid-gap upwards. However, the movement of the threefold degenerate defect level up into the band gap (section V B) leads us to expect six transition levels all lying near (above or below) the valence band edge. The most stable charge state at the valence band edge itself could be anything from 0 to +3.

VII. THE ORIGIN OF THE LINEAR SCALING TERM FOR NEUTRAL DEFECTS.
We expect linear scaling terms to arise from two sources: from incomplete relaxation (elastic errors) and from electrostatic errors. For electrostatic errors linear scaling comes 2 from the monopole term in the multipole expansion of the electrostatic interactions: roughly speaking, it is the Madelung energy of the localized charge on the defect interacting with the compensating jellium background in the image cells. This is shown schematically in Fig 9. Part a) shows the typical situation for a positively charged defect: An infinite array of very tightly localized positive core charges (one in each periodic image of the supercell) are only partially compensated by the localized electrons in the (less tightly) localized defect levels. The remaining charge is compensated by jellium. Makov and Payne 2 extracted the part of the system shown in Fig 9 b) consisting of an infinite array of positive delta functions interacting with jellium and showed that is gives rise to a formation energy error which is linear in the supercell size, with a strength proportional to the square of the charge on the defect. This error should therefore be absent for unrelaxed neutral defects. The fact that we see a clear linear contribution to the unrelaxed formation energy of, for example, the neutral In P and tetrahedral In interstitials is thus somewhat unexpected.
The explanation lies in the localization/delocalization of the states which become occupied and unoccupied when the cell contains a neutral defect. In the cases of the neutral interstitials and vacancies one or more defect level(s) lie outside the bandgap for all or some supercell sizes. As a result there are either electrons in delocalized band states at the LDA conduction band edge or holes at the valence band edge. This is shown for the example of In i(P) in Fig 10. To quantify the localization of a particular Kohn-Sham eigenstate we start by projecting it onto Wigner-Seitz cells 15 around each atom. We then select some radius r l around the defect and find the average projection per atom ρ > over the part of the cell where r > r l , and similarly the average ρ < for the r < r l region. The localization in then given by the ratio For an isolated defect in an otherwise perfect infinite lattice, α ρ → ∞ as r l → ∞ if the state is localized, since ρ > → 0. On the other hand, α ρ → 1 if the state is delocalized as the two averages ρ > and ρ < then tend to the same value.
Here we have finite L sized supercells, so we choose r l such that ρ > is calculated over the outermost complete shell of atoms around the defect, plus the atoms in the cell corners. . To a first approximation we can replace the charge density of the delocalized electron by its average value, thus recovering the situation of Fig 9 a). We thus predict a linear term in the formation energy scaling. However, the distribution of the delocalized charge is in reality far from uniform on the scale of the atomic spacing, so it is not clear what the prefactor and effective charge should be. Hence, without performing the detailed mathematical derivation required (which lies beyond the scope of the current paper) we cannot predict what gradient the linear term should have. The case of In P is different. There are no electrons in delocalized conduction band states or holes in the valence band. However, the localization of the defect level at the Γ point is rather weak, certainly compared to P In and In i(P) . The defect level of In P is three-fold degenerate and partially filled. Away from the Γ point this degeneracy is split by the interaction of the defect with its images, leading to different dispersions of the defect levels in different directions in k-space, with one level having lower energy (away from Γ). At most k-points this state is thus completely filled, whilst the other two are partially empty. Hence, averaging the Wigner-Seitz projection over the whole Brillouin zone (using 8x8x8, 4x4x4 or 2x2x2 The other two levels are more occupied at the origin in k-space than elsewhere and hence in real space are only partially localized. Clearly the Jahn-Teller structural relaxations are required in order to properly localize these states. In the unrelaxed structure we effectively have two localized electrons and two partially delocalized ones in the middle of the band gap, leading to the spatial charge distribution shown schematically in Fig 9d). Electrostatically, this again behaves to first order like a jellium background charge surrounding a positively charged, localized defect, resulting in the linear scaling observed. As with the interstitials, predicting the specific gradient expected must be left to future work. Here, we simply note that this mechanism should also affect the unrelaxed formation energies of other defects with partially filled degenerate levels in the band gap. It seems reasonable to presume that it is involved with the neutral V In and perhaps V P , although they also have linear terms arising from partially filled band states.

VIII. OTHER SOURCES OF ERROR.
Finite size errors are not the only types of errors present in these calculations: errors also arise from the truncation of the basis set and the k-point integration as well as from the use of both pseudopotentials and LDA. Furthermore, memory size limitations forced us to use pseudopotentials in which the indium 4d shell is treated as core rather than valence. Although the central aim of this paper is to study the treatment of finite size errors it is still informative to estimate the size of these other sorts of errors for comparison. In principle, errors arising from the pseudopotentials and the LDA should be independent of supercell size, though some short range dependence may still be anticipated since changes may affect the amount of defect band dispersion. This contribution should disappear exponentially with supercell size.
In fig 11 we show the change in the unrelaxed formation energies when (left panel) the In 4d electrons are treated as core rather than valence, (centre panel) the US-PP are replaced by the Projector Augmented Wave method 16 (PAW) or (right panel) the LDA exchange correlation functional is replaced by the Generalized Gradient Approximation (GGA) of Perdue and Wang 17 . These formation energy differences are shown as a function of supercell size for the seven stable point defects. (V In , In P , In i(In) , In i(P) , V P , P In and P i(HEX) .) For the GGA and for LDA with the In 4d as valence a planewave cutoff of 200 eV gives the same level of accuracy as it does with the In 4d as core. For PAW on the other hand a planewave cutoff of 300 eV is required for the same accuracy level. The relaxed lattice constants and bandgaps also change. For LDA with In 4d as valence we find 5.833Å and 0.581 eV respectively, for LDA with PAW we find 5.830Å and 0.597 eV whilst for GGA (with US-PP and In 4d as core) we find 5.956Å and 0.473 eV.
The errors coming from the pseudopotentials (left and centre panels of Fig 11) are, as anticipated largely independent of the supercell size. The slight dependence is well described by a two parameter fit to the general exponential (Solid lines in the figure.) Changing the In 4d electrons from core to valence has, not surprisingly, only a fairly small effect (well below 0.1 eV) upon the phosphorus related defects, but a much more significant effect upon the indium related ones, particularly In i(In) and In P where the effect is on the same order of magnitude as the finite size errors. PAW produces more accurate results than any pseudopotential method since it reconstructs the exact valence wavefunction with all nodes in the core region. The replacement of US-PP with PAW produces only small (O(0.1 eV) or less) changes in the formation energies with virtually no size dependence: the largest difference between the correction in the 8 atom cell and that in the 64 atom cell is only 0.01 eV, so the US-PP to PAW changes have not been calculated in the 216 atom cell. This demonstrates that the widely used US-PP 9 are perfectly reasonable for this type of calculation. It should also be noted that in most cases the small correction when PAW is introduced partially cancel that made when moving the In 4d electrons from core to valence, at least for the examples studied here. The errors arising from the exchange-correlation functional have two main forms. Firstly, the band gap is (usually) strongly reduced compared to experiment. (For InP GGA does worse than LDA once lattice parameter optimization has been included.) This reduction leads to ambiguities in the definition of the formation energy for charged defects, in turn leading to large uncertainties in predictions. For some semiconductors the bandgap can even be reduced to zero making the material appear metallic and dramatically altering the properties of many defects. However, neither of these ef-fects occurs here, since we consider neutral defects and obtain a non-zero band gap. Never-the-less, a second exchangecorrelation related error is present: LDA overbinds all bonds, moving some defect formation energies up and others down. Hence to assess the errors involved we here compare with GGA, which is known to have the opposite effect: underbinding. An exact solution to the DFT formation energies would lie somewhere in between -and probably closer to the LDA results since LDA gives a better lattice constant for bulk InP. We note also that we have not allowed spin polarization in our calculations, using LDA instead of the Local Spin Density Approximation (LSDA). However, for most of our defects the use of LSDA would simply cause further finite size errors since it would introduce spurious magnetic interactions between the defect and it's PBC images. The exceptions occur for the Jahn-Teller active defects (In P , V P and V In ) where the degeneracies can in certain cases be lifted by Hund's rule couplings. However, this is only important for materials 18 with more tightly localized bonds or dangling bonds than we have here, so we may safely omit it.
The results from replacing LDA with GGA are shown in the right panel of Fig 11 and are somewhat surprising: there is a clear finite size error associated with the choice of exchangecorrelation functional for certain defects, namely In i(P) , In i(In) , V P and P i(HEX) . The antisites on the other hand show no significant supercell size dependence. (V In shows only a small size dependence but is still fitted better by a polynomial than by an exponential.) The defects for which the LDA/GGA error difference depends on supercell size are those for which the defect levels lie outside the band gap. This indicates that the spurious delocalized states are treated differently by different exchange correlation functionals. This is confirmed by checking the degree of localization of the electrons/holes at the band edges: when GGA is used the values of α ρ are 20 to 400% larger than with LDA. These changes are linked to the change in band gap and also to the fact that the defect bands move closer to the band edges (at least near the Γ point) for these particular defects. Hence we may expect large, supercell size dependent errors of O(0.5 eV) from the choice of exchange correlation functional when partially filled defect bands lie outside the bandgap but we expect smaller O(0.1 eV) size independent errors when the defect levels lie within the band gap. Table IV compares the size of the errors arising from different sources. The basis set errors are estimated from the difference in formation energy in the 64 atom supercell when the planewave cutoff is raised from 200 eV to 400 eV. The kpoint errors are different in each supercell: the errors shown in the table are the largest occurring for any of the supercells used, estimated as the difference between the mean value E C d and the value E C d (N max ) obtained using the largest k-point grid actually calculated for the defect and supercell in question. All other sources of error -such as defect band dispersion etc -are contained within the finite size errors shown. Most of the errors listed are on the 0.1 eV scale or below, which in practical calculations is usually acceptable. The finite size errors and some of those arising from treating the In 4d electrons as core are larger, lying on the ∼0.5 eV level.  We note, however, that for charged defects we anticipate even larger finite size errors -up to 1-2 eV in many cases 12 -whilst we do not anticipate the errors from the In pseudopotentials being significantly different from those here. The presence of defect level dispersion effects is confirmed by the existence of the exponentially shrinking supercell size dependence in the errors related to pseudopotentials and to the LDA for the antisites. These results also confirm that they are short ranged. Their energy scale is rather hard to estimate directly, but the size of the exponential components in Fig 11 suggests that the errors involved are only on the 0.01-0.1 eV scale even for the 8 atom cell. An alternative estimation comes from the error bars on the scalings themselves (see table I) which are 0.01-0.35 eV, only a small part of which comes from the defect level dispersion.
The fact that the errors coming from the pseudopotentials have only a small and exponentially decaying cell size dependence means that it is perfectly reasonable to make the approximations we needed to make in order to be able to do calculations in sufficiently large supercells to correctly assess the finite size errors. The amplitude and sign of these nonsize dependent errors can be calculated separately in a smaller cell -the 64 atom cell for example -and simply added or subtracted from the final finite size scaled results in order to produce much more accurate and reliable defect formation energies than those normally published. This has been done for the LDA formation energies of the present examples in the last column of table IV, but is equally valid for the errors in the relaxed formation energies and those in the GGA (at least if no delocalized hole/electron states have appeared at the valence/conduction band edges.)

IX. CONCLUSIONS.
We have studied the finite size errors which occur when the supercell approximation is used in the calculation of the formation energies and structures of point defects in semiconductors, using the neutral native defects of InP as an example. We have calculated the relaxed and unrelaxed formation energies using planewave ab initio DFT in simple cubic supercells containing 8, 64, 216 and 512 atoms -the largest currently computable. To examine and correct for these errors we have used finite size scaling with inverse supercell size, which we consider to be the most reliable and accurate way to treat the finite size supercell approximation errors. Unlike other methods this does not rely upon any modelling or assumptions about the errors, other than that they are primarily long ranged (polynomial rather than exponential) and decrease with increasing supercell size. This method requires the results of calculations in at least 3 supercells, and at least 4 if we are to have an idea of the accuracy of the resulting scaled energies. Hence for some difficult cases in which the 8 atom cell is simply too unreliable it may occasionally be necessary to use supercells with up to 1000 atoms.
Three sources of finite size error have been examined: in the case of relaxed formation energies there are elastic errors due to the finite volume available for relaxation. We showed that, as they should, these scale linearly with inverse supercell size (L −1 ) with very little hint of any higher order error term arising, even when Jahn-Teller distortions are taken into account. The second type of error is the dispersion of the defect levels, which has only a relatively small effect upon the formation energies, at least when the defect levels are completely filled. These effects appear to shrink exponentially with increasing supercell size, as anticipated, but appear to slightly increase the uncertainties in the final scaled formation energies. The third source of error is much more significant and arises for both relaxed and unrelaxed formation energies. It is due to charged multipole interactions between a defect and its images in the PBCs. We have shown that these errors are present and not always negligible, even for neutral defects. Linear errors can arise if fully or partially filled defect levels lie outside the band gap in the neutral charge state, leading to delocalized holes at the valence band edge or electrons at the conduction band edge. Linear errors can also appear in unrelaxed formation energies due to the way in which the defect/image interactions lift degeneracies for partially filled degenerate defect states within the band gap. Both of these ways for linear scaling errors to arise can apply for other non-neutral defects. They indicate that the calculation is in some way unphysical, for example that the charge state involved is not actually stable. However, in practise it is often necessary to calculate formation energies of such unstable charge states in order to check which transitions levels do actually lie inside the gap. It is thus important to be aware that large finite size errors, such as those reported here, can occur even in calculations for neutral defects, as this can lead to transition levels calculated in individual supercells (without scaling) to appear to lie inside the gap when they should lie outside and vice versa. Indeed, electrostatic errors are still present even for physically rea-sonable cases, such as V P and (probably) V In , due to defect states which enter the gap as the cell size grows. They are even present and may remain significant for defects such as the neutral P In which has all of its defect levels within the gap in all supercells. This is because neutral defects in crystalline solids still have higher order charge multipole moments, especially if the system is made up of more than one type of atom with different electronegativities. We have shown that for unrelaxed formation energies these errors scale as the inverse cube of the supercell size (L −3 ). For relaxed formation energies this is almost certainly the case also, although we do find possible indications that the leading non-linear error term may sometimes scale as the inverse square (L −2 ). Further work with more defects is needed to confirm or definitively rule out this possibility. It could clearly be answered by removing the 8 atom cell from the scaling and replacing it with the 1000 atom cell, but that must wait for improved computing facilities. In the meantime, an answer may be obtained from scaling studies of formation energies for charged defects 12 since both the electrostatic errors themselves and the cross terms between them and the elastic errors will then be stronger.
To summarize: the use of large (up to 500 or occasionally even 1000 atom) supercells with finite size scaling has been shown to be a very promising route around the errors which arise in the use of the supercell approximation to calculate for-mation energies of defects in III-V (and other) semiconductors. Errors scale with a linear plus a higher order term, most probably cubic. We have also found several instances where scaling recovers physically relevant results that are even hidden in calculations on the 512 atom cell: formation energies which are wrong by ∼ 1 2 eV and defect levels which appear inside the valence or conduction bands in supercell calculations when they should actually lie inside the band gap.