Incorrect DOS when using the tetrahedron method

[tmueller]

The documentation of the DOSCAR file says the DOS is calculated as follows, which conserves the total number of electrons exactly:



However this does not appear to be the case if the tetrahedron method is used. The DOS for the tetrahedron method appears to be calculated as the continuous derivative (rather than discrete difference) of the integrated DOS. As a result, when adding up the DOS the number of electrons is not conserved, and some sharp peaks in the DOS do not show up at all.

A manifestation of this problem has been posted previously on this forum. The problem is partially addressed in VASP 6, which implements special treatment for the situation in which a tetrahedron is contained entirely between successive energy levels and thus does not show up in the density of states. However this does not fix the broader problem, and adding up the DOS will still not conserve the number of electrons in VASP 6.

[martin.schlipf]

I agree that the tetrahedron method will produce the DOS as you report. If a grid-point has any DOS at all, it will provide the analytic derivative of the integrated DOS. If there are two consecutive points where the integrated DOS changes but the analytic DOS would be zero, one of the points will be arbitrarily assigned to get the DOS.

What I currently do not understand is what you would like to see changed. Would you like an update to the documentation that clarifies this point a bit better or do you want a change to the code do that it is a finite-difference derivative? We agree that this is an artifact of the finite energy grid fixed by increasing NEDOS as needed, don't we?

[tmueller]

Thanks for the reply. I think it would be better to change the code. That would make the behavior consistent with other values of ISMEAR and with what I think most users would expect. Ideally the change would be documented somewhere, so that in the future people are able to appropriately interpret the densities of states generated with earlier versions of VASP.

I'm not sure I'd say that this problem can be fixed by increasing NEDOS, for both theoretical and practical reasons. From a theoretical perspective, increasing NEDOS only reduces the error. Recovering the behavior described in the documentation ("This method conserves the total number of electrons exactly") can only be achieved when NEDOS is infinite.

From a practical perspective, it's hard to know what value of NEDOS would result in a reasonably small error. I've run some tests using VASP 6.3.1 on mp-567498 from the Materials Project, and even when NEDOS is 50,001 the difference between the VASP-calculated integrated-DOS and adding up the DOS is still about 2%. VASP 5.4.4 has a 10% error when NEDOS is 50,001.

Input files for this calculation can be found here, using the CHGCAR from here to initialize the calculation.

When doing these tests, I noticed that for some reason the execution time of VASP 6.3.1 seems to scale poorly with NEDOS. The calculation took about 18 times as long with with NEDOS=50001 as it did with NEDOS=501. I don't have the same problem with VASP 5.4.4.

[martin.schlipf]

Could you comment why you need the DOS to have the property of conserving the number of electrons? I currently don't see where this feature would be required and we are conservative about breaking existing behavior unless there is a very strong reason for it. You always have the integrated DOS as well, so you can compute the numeric derivative if that is necessary.
We updated the documentation on the wiki to state that the number of electrons is only conserved for ISMEAR >= -1.

Regarding the particular large mismatch for VASP 5.4.4. I would suspect that one tetrahedron has all 4 corners degenerate. This can happen for systems with large symmetry and then it would never contribute to the DOS.

[tmueller]

It's possible for a user to work around this problem by writing some post-processing code to extract the correct (electron-conserving) density of states from the integrated density of states, but I suspect many users will not do this. For example, I noticed this problem because I was using DOS data from the Materials Project database. They directly report the VASP outputs without doing any post-processing to correct the densities of states, so their densities of states do not conserve the number of electrons and completely miss some bands. Unfortunately, that makes the data nearly useless for my project.

The Materials Project data was generated using VASP 5, but even with VASP 6 the error can be significant, both qualitatively and quantitatively. Below are a couple of DOS plots for the example I listed above, run at values of NEDOS ranging from 1000 to 5000. The plots labeled "raw" were generated using the raw VASP output, and those labeled "corrected" were extracted from the VASP-generated integrated density of states. The corrected plots are more consistent both qualitatively and quantitatively.

movie_zoomed.gif

movie_zoomed.gif

Of course it is entirely up to you as to whether you change your code, and I understand the hesitancy to break existing behavior. However the VASP 5 behavior has already been broken in VASP 6. The updated DOSCAR documentation is helpful, but it does not mention the changes made in VASP 6, which are significant. Whatever you decide to do about this, I think it's important that the behaviors of the different versions of the code are well-documented.

[martin.schlipf]

We will discuss this internally. Perhaps one could set a flag or additional ISMEAR options to change to the behavior that you want.

Regarding the correctness of the DOS: Both methods have advantages and disadvantages. Taking the numerical derivative guarantees that you can numerically integrate it to get the correct integrated DOS. However, the DOS at a particular energy E depends on the energy mesh. The analytic derivative is of course independent of the mesh, but can "miss" states if your mesh happens not to include them. So you need a mesh spacing smaller than the width of your bands. For the plots that you show, a denser k-point mesh would probably also improve the results.

But ultimately it comes down to what you want to use the DOS for. Perhaps, you could tell us why you need the DOS to be consistent with the integrated DOS? If you want to calculated integrals involving the DOS, either method of obtaining the DOS is valid and which one is better suited depends on the quantity that you want to integrate.

[tmueller]

Thanks. Without going into too much detail, I'm doing some machine learning with densities of states. The fluctuations in the total electron count and the appearance and disappearance of peaks with respect to NEDOS are effectively noise in the data. VASP 6 is somewhat better than VASP 5, but the challenge with VASP 6 is that it's not straightforward to interpret the DOS data. Sometimes it is the analytical derivative, sometimes it is the numerical derivative, and sometimes (I believe) it's a hybrid of both.

[martin.schlipf]

We discussed this today and plan to implement a new INCAR tag that allows you to switch between numeric and analytic derivative for the DOS. Unfortunately, this will likely not make it into the next version 6.4. For your research, we recommend that you take the integrated DOS and take the numeric derivative by hand (or script). But you also need to be aware that the density of the k-point mesh might be something you need to carefully consider. The spikes you see at the moment seem to at least partially caused by a too coarse grid.