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): \begin{itemize} \item 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),
\item semiconductors with a rather small dielectric constant $\varepsilon$, and \item high-charge states like 3+, 4+, 3- or 4- since the spurious interactions
scale (approximately) proportional to the {\em square} of the total cell
charge, e.g., for a 4+ state the error is about 16 times larger than for
a 1+ state!
\end{itemize} 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: \begin{itemize} \item Unfortunately a {\em 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).
\item 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 {\em 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)!
\item 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
{\em without} any corrections and later you have to add the corrections
``by hand using the output printed in OUTCAR (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 {\em 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
{\em first-principles} calculations usually agree very well (often within
the order of 1-3\%).
\item \index{INCAR!I!IDIPOL|textit}\index{INCAR!D!DIPOL|textit}
The dipole-dipole plus quadrupole-monopole corrections printed in OUTCAR
are meaningless in their original form! We have to calculate a correction
for the {\em 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 {\tt IDIPOL}=4 and {\tt DIPOL}=same position as in defect
cell (search for the line containing Tr[quadrupol] ... in file OUTCAR).
The corresponding quadrupole moment has to be subtracted from the
quadrupole moment printed for the charged defect cell. The difference
corresponds to the {\em 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 ..." of
file OUTCAR 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 file OUTCAR). Add now the
corrected value for the quadrupole-monopole interaction to the calculated
monopole-monopole interaction energy (and finally scale the sum with
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!) {\em 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 {\em defect
cell} you can now calculate the quadrupole-monopole contribution alone and
hence, the dipole-dipole contribution is then known too. The dipole-dipole
contribution will be kept and the {\em defect-induced} quadrupole-monopole
contribution has to be added to this (just multiply the proportionality
constant with the {\em 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.
\item \index{INCAR!L!LDIPOL|textit}
Any potential correction ({\tt LDIPOL}=.TRUE.) is currently {\em impossible}!
Hence you can only use {\tt LDIPOL}=.FALSE.! The reasons are: first the
downscaling with $\varepsilon$ is missing and second the correction is
not calculated from the {\em 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 routine dipol.F!). However, one has to expect that
the potential corrections do not change the results dramatically ... .
\end{itemize}
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 {\em acceptor-donor pairs}. The most critical cases are again given for semiconductors with rather small lattice constants, rather small dielectric constants of 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, {\tt 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: since in cubic cells any dipole moment can only be {\em defect-induced} no additional corrections are necessary (in contrast to the monopole-quadrupole energies of charged cells). However, the other bad news is: for such defect complexes it may sometimes be hard to find the correct ``center of mass (input {\tt DIPOL}=... in INCAR!) for the {\em defect induced} charge perturbation (it's usually more easy for single point defects since usually {\tt DIPOL}=position of the point defect is the correct choice). This introduces some uncertainties and one might try different values for {\tt DIPOL} (the one giving the minimum correction should be the correct one). But also note: {\tt DIPOL} is internally aligned to the position of the closest FFT-grid point in real space. Hence, the position {\tt DIPOL} is only determined within distances corresponding to the FFT-grid spacing (controlled by {\tt NGXF}, {\tt NGYF}, and {\tt NGZF}). As an additional note 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 {\tt 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 {\em 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 {\bf 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 {\bf k}-point sampling (though one should note that this may then hold for the electrostatic corrections too!).