Dipole corrections for defects in solids

Similar to the case of charged atoms and molecules in a large cubic box also charged defects in semiconductors impose the problem of potentially slow convergence of the results with respect to the supercell size due to spurious electrostatic interaction between defects in neighboring supercells. Generally, the errors are less dramatic than for charged atoms or molecules since the charged defect is embedded in a dielectric medium (bulk) and all spurious interactions between neighboring cells are scaled down by the bulk dielectric constant $\varepsilon$. Hence, the total error might remain small (order of 0.1 eV) and one has not to worry too much about spurious electrostatic interactions between neighboring cells. However, there exist three critical cases where one should definitely start to worry (and to apply dipole corrections):

• Semiconductors containing first-row elements since they possess rather small lattice constants and hence the distance between two neighboring defects is smaller than in most other semiconductor materials (though one should note that the smaller lattice constant alone must not yet increase the errors dramatically since the leading scaling is $1/L$, only the contributions scaling $1/L^3$ may become dangerous for small cells).
• Semiconductors with a rather small dielectric constant $\varepsilon$.
• High-charge states like 3+, 4+, 3- or 4- since the spurious interactions scale (approximately) proportional to the square of the total cell charge, e.g., for a 4+ state the error is about 16 times larger than for a 1+ state!

The worst case one can ever think of is that all three conditions mentioned above are fulfilled simultaneously. In this case the corrections can amount to the order of several eV (instead of the otherwise typical order of few 0.1 eV)!

In principle it is possible to apply the same procedure as in the case of charged atoms and molecules in vacuum. However, with the current implementation one has to care about following things and following restrictions apply:

• Unfortunately a full correction is only possible for cubic cells, the only contribution which can always be corrected for any arbitrary cell shape, is the monopole-monopole interaction. However, for intermediate cell sizes the quadrupole-monopole interaction is not always negligible (it can reach the order of minus 30-40 % of the monopole-monopole term). Therefore, whenever possible the use of cubic cells is recommended. Otherwise one should try to use as large as possible cells (the dipole-dipole and monopole-quadrupole interactions scale like $1/L^3$ and therefore, for larger cells a monopole-monopole correction alone becomes more and more reliable).
• The corrections are only reasonable if the defect-induced perturbation of the charge density is strictly localized around the defect, i.e., if only the occupation of localized defect states is changed. Whenever the problem occurs that (partially) wrong bands (e.g. delocalized conduction band or valence band states instead of defect states) are occupied the calculated corrections become meaningless (the correction formulas are not valid for overlapping charges)! Therefore one should first calculate the difference between the charge densities of the charged defect cell and the ideal unperturbed bulk cell and check the localization of this difference charge (in between the defects the difference must vanish within the numerical error bars for the charge densities)!
• Don't forget to scale down all results by the bulk dielectric constant $\varepsilon$! Yet, there is no possibility to enter any dielectric constant, all corrections are calculated and printed for $\varepsilon=1$. Therefore, the corrected total energies printed after the final electronic iteration are meaningless! Hence, you should first calculate the energies without any corrections and later you have to add the corrections by hand using the output printed in the OUTCAR file. You must search for a line DIPCOR: dipole corrections for dipole and following lines. There you find the dipole moment, the quadrupole moment and the energy corrections. One should note that strictly one has to take the dielectric constant calculated by first-principles methods. Since VASP does not yet allow a simple calculation of dielectric constants, however, you have to use the experimental value (or values taken from other calculations). This empirism introduces slight uncertainties in your energy corrections. However, one can expect that the uncertainty should rarely exceed 5-10% since dielectric constants taken from experiment and those obtained from first-principles calculations usually agree very well (often within the order of 1-3%).
• The dipole-dipole plus quadrupole-monopole corrections printed in OUTCAR are meaningless in their original form! We have to calculate a correction for the defect-induced multipoles, but since we have also included the surrounding bulk a quadrupole moment associated with the corresponding charge (extending over the whole cell!) is also included in the printed quadrupole moment (and in the corresponding energy corrections). Since in systems with cubic symmetry dipoles are forbidden by symmetry a dipole moment can only be defect induced (and only if the cubic symmetry is broken by atomic relaxations). In order to obtain the correct (usually quadrupole-monopole interaction only) energy correction, one has to proceed as follows:
  • One has to calculate the quadrupole moment for an ideal bulk cell (neutral!) by setting IDIPOL=4 and DIPOL= same as in the defect cell (search for the line containing Tr[quadrupol] ... in the OUTCAR file).
  • The corresponding quadrupole moment has to be subtracted from the quadrupole moment printed for the charged defect cell. The difference corresponds to the defect-induced part of the quadrupole moment.
  • If no dipole-dipole interaction is present you can now simply scale down the energy printed on the line dipol+quadrupol energy correction ... in the OUTCAR file by the ratio defect-induced quadrupole/total cell quadrupole since this interaction is proportional to the quadrupole moment. After this scaling you should end up with reasonable numbers (usually smaller than the monopole-monopole correction printed on the line containing energy correction for charged system ... in the OUTCAR file).
  • Add now the corrected value for the quadrupole-monopole interaction to the calculated monopole-monopole interaction energy (and finally scale the sum wit $1/\varepsilon$).
  • The whole procedure is even more complicated if a dipole moment occurs also, since then only the quadrupole-monopole term has to be corrected but the dipole-dipole term is already correct!
  • But you can easily help yourself: Take simply a cell of the same dimension and calculate a free ion (does not matter which one) of the same charge state (if this causes trouble try the opposite state, e.g. 4+ instead of 4-, but don't forget then to take the opposite sign for the printed monopole quadrupole energy since this energy is proportional to the cell charge). The calculation will provide a quadrupole moment and a certain quadrupole-monopole interaction energy. Since this energy is proportional to the quadrupole moment (times total cell charge) you can estimate the proportionality constant with which one has to multiply the quadrupole moment in order to obtain the corresponding monople-quadrupole interaction for the given cell size by dividing the energy by the quadrupole moment. Multiplying this constant by the quadrupole moment of the defect cell you can now calculate the quadrupole-monopole contribution alone and the dipole-dipole contribution is then known too. The dipole-dipole contribution will be kept and the defect-induced quadrupole-monopole contribution has to be added to this (just multiply the proportionality constant with the defect-induced} quadrupole moment).
  • Then you finally end up with the correct values for all interactions (which have to be summed again and rescaled with $1/\varepsilon$). It's currently a clumsy procedure but it works satisfactorily.
• Any potential correction (LDIPOL=.TRUE.) is currently impossible! Hence you can only use LDIPOL=.FALSE.! The reasons are:
  • First the downscaling with $\varepsilon$ is missing and second the correction is not calculated from the defect-induced multipoles but from the total monopoles of the defect cell containing at least a meaningless quadrupole contribution (one had to subtract the quadrupole moment of the ideal cell before calculating any correction potential, but this is not yet implemented in the routine dipol.F). However, one has to expect that the potential corrections do not change the results dramatically ...

Besides charged defects there's another critical type of defects which may cause serious trouble (and for which one should also apply dipole corrections): neutral defects or defect complexes of low symmetry. For such defects a dipole moment may occur leading to considerable dipole-dipole interactions. Though they fall off like $1/L^3$ they might not be negligible (even for somewhat larger cells) if the induced dipole moment is rather large. The worst case that can happen is a defect complex with two (or more) rather distant defects (separated by distances of the order of nearest-neighbor bond lengths or larger) with a strong charge transfer between the defects forming the complex (e.g., one defect might possess the charge state 2+ and the other one the charge state 2-). This can easily happen for defect complexes representing acceptor-donor pairs. The most critical cases are again given for semiconductors with rather small lattice constants, rather small dielectric constants and for any defect complex causing strong charge transfers. Again the same restrictions and comments hold as stated above for charged cells:

• You may currently only use cubic cells, LDIPOL=.FALSE. and you have to rescale the correction printed in OUTCAR by the bulk dielectric constant $\varepsilon$ (i.e., the printed energies are again meaningless and have to be corrected by hand). There is only one point which might help, that is since in cubic cells any dipole moment can only be defect-induced no additional corrections are necessary (in contrast to the monopole-quadrupole energies of charged cells). *The other bad news is that for such defect complexes it may sometimes be hard to find the correct center of mass (DIPOL in INCAR file) for the defect induced charge perturbation (it's usually more easy for single point defects since DIPOL=position of the point defect is the correct choice). This introduces some uncertainties and one might try different values for DIPOL (the one giving the minimum correction should be the correct one).
• DIPOL is internally aligned to the position of the closest FFT-grid point in real space. Hence, the position DIPOL is only determined within distances corresponding to the FFT-grid spacing (controlled by NGXF, NGYF, and NGZF). This might also play a certain role if for charged single point defects the position of the defect is not chosen to be (0,0,0)! In this case DIPOL might correspond to a position lying slightly off the position of the defect what may also introduces inaccuracies in the calculation of the electrostatic interactions (i.e., apparent dipole moments may occur which should be zero if the correct position DIPOL would have been chosen). In this case you should whenever possible try to adjust your FFT-grid in such a way that the position of the defect matches exactly some FFT-grid point in real space or otherwise never use any other (point) defect position than (0,0,0)...

A final note has to be made: besides the electrostatic interactions there exist also spurious elastic interactions between neighboring cells which (according to a simple elastic dipole lattice model) should scale like $1/L^3$ (leading order). Therefore, the corrected values may still show a certain variation with respect to the supercell size. One can check the relaxation energies (elastic energies) separately by calculating (and correcting) also unrelaxed cells (defect plus remaining atoms in their ideal bulk positions). If the k-point sampling is sufficient to obtain well-converged results (with respect to the BZ integration) one might even try to extrapolate the elastic interaction energies empirically by plotting the relaxation energies versus $1/L^3$ (hopefully a linear function; if not try to plot it against $1/L^5$ and look whether it matches a linear function) and taking the value for $1/L \rightarrow 0$ (i.e. the axis offset). However, usually the remaining errors due to spurious elastic interactions can be expected to be small (rarely larger than about 0.1 eV) and the extrapolation towards $L \rightarrow \infty$ may also be rather unreliable if the results are not perfectly converged with respect to the k-point sampling (though one should note that this may then hold for the electrostatic corrections too).