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 . 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 , only the contributions scaling may become dangerous for small cells),
- semiconductors with a rather small dielectric constant , and
- 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!

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 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
! Yet, there is no possibility to enter any dielectric
constant, all corrections are calculated and printed for
.
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 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*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 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*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/ ). 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 hence, 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/ ). 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 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 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 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 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, `LDIPOL`=.FALSE. and you have to rescale
the correction printed in OUTCAR by the bulk dielectric constant
(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 *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 `DIPOL`=... in INCAR!) for the
*defect induced* charge perturbation (it's usually more easy for single
point defects since usually `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).
But also note: `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`).
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 `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
(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 (hopefully a linear function
- if not try to plot it against and look whether it matches a linear
function) and taking the value for
(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
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!).

N.B. Requests for support are to be addressed to: vasp.materialphysik@univie.ac.at