What are some approaches to modeling charged point defects in materials?

Question

From what I understand, one major issue in modeling charged point defects in materials is the existence of artificial interactions between charged defects and their periodic images. Even when using large supercells, as is typical when modeling uncharged point defects, the problem persists.

In the past, I briefly used the PyCDT Python package, which uses sxdefectalign, and VASP to calculate the formation energies of charged point defects. Besides using plane-wave DFT and the correction scheme in sxdefectalign, are there other approaches to modeling charged defects that might trade speed for accuracy or vice-versa? What other challenges are associated with modeling charged defects?

Answer 1

This is an important question that isn't answered yet, and I've been doing some work with charged defects recently so I'll attempt an answer--though I readily admit I am not really an expert (i.e. I don't work on developing these correction approaches). I'm not familiar with many approaches apart from DFT, so I'll limit my answer to that area.

In DFT at least, a charged cell is compensated by a uniform jellium background of equal and opposite charge, to keep the energy finite. The problem is that the long-range potential of the charged defect in this medium decays very slowly, as you mentioned. It depends on the defect and the dielectric properties of the material under study, but can be significant (on the order of eV). So something needs to be done to correct for this spurious interaction. There are also several potential-alignment terms that come into play when calculating energy differences between different systems. A good place to start reading is references 1 and 2 below. I'll keep this answer less theoretical since the references do a much better job of explaining. I'll summarize the concepts quickly and mention some practical issues.

One approach you may have heard of is the Makov-Payne correction, which was derived based on a Madelung-type sum in a cubic cell. I haven't used this correction but I understand it can be quite inaccurate in realistic systems, typically overcorrecting 3. In the situations I've encountered it, it's also been limited to systems with cubic symmetry, though I'm not sure if this is always the case.

Freysoldt's scheme is a common approach. The basic idea is that you use a simple model to describe the defect's charge, such that you can calculate it's isolated energy interacting with the jellium background, as well as the energy of a periodic system that is interacting with it's images, using Poisson's equation. If you align the potential of your model with the DFT calculations, you can use the difference between the isolated and periodic energies of your model as a correction term. There is also related correction by Kumegai and Oba 4.

The third correction I've seen is from Lany and Zunger, which goes beyond the Makov-Payne correction by using DFT-calculated difference in charge between a charged and neutral defect to calculate a higher-order term in the correction. I have the least amount of experience with this method.

When it comes down to implementing these schemes I have the most experience with the Freysoldt method. I have used sxdefectalign and CoFFEE. One thing I learned early on was that all the nice plots in the papers and code examples are typically unrelaxed calculations where the atomic positions are the same in all the calculations. This makes differences in potentials smooth and easy to analyze. When you want an accurate formation energy you need to allow the system to relax; the differences in the atomic positions of the different calculations (pristine vs. defect) cause significant variation in the potential. You typically need to do some smoothing or averaging, which can be challenging to get right.

I've also encountered challenges when using these codes with cells that have non-orthogonal lattice vectors (e.g. monoclinic). Planar averaging and model calculations can be especially difficult. It can be hard to determine if you're doing something wrong or if there is a bug in the code somewhere. You need to do a lot of experimenting.

There are other codes that automate the process more. You mentioned PyCDT; PyDEF and Pylada are also interesting, but these three only really support VASP at the moment. I know PyCDT has a wrapper for sxdefectalign, but I thought they had a separate module for performing the calculation without sxdefectalign. Pymatgen can perform various corrections taken from the PyCDT code, but I haven't had much time to spend with it yet.

I'm always interested to hear if there are other approaches and codes out there.

References

1. C. Freysoldt et al., Rev. Mod. Phys. 86, 253 (2014).
2. H.-P. Komsa, T. T. Rantala, and A. Pasquarello, Phys. Rev. B 86, 045112 (2012).
3. H.-P. Komsa, T. Rantala, and A. Pasquarello, Physica B: Condensed Matter 407, 3063 (2012).
4. Y. Kumagai and F. Oba, Phys. Rev. B 89, 195205 (2014).