Fermi level in VASP and all the other ab-initio codes is calculated using the Fermi-Dirac distribution. This is because electrons in solids and molecules obey the Fermion statistics, which prohibits double occupation of any state with particles of the same quantum character. Accordingly, in any configuration, the electronic states are filled in ascending energetical order by at most two electrons (one for each spin orientation), until all the available electrons are accommodated. The energy of the highest occupied state then becomes the Fermi energy of the whole system.

In codes such as VASP, to properly treat the step-like behaviour of the Fermi-Dirac distribution, different methods have been implemented. The simplest one is the so-called Gaussian smearing, which defines an effective broadening for the otherwise desecrate electronic states, allowing smooth sampling of the entire Brillouin zone using a finite number of k-points. There are more advanced methods such as Methfessel-Paxton, which can provide a more accurate description for the partial occupancies and the resulting Fermi energy.

Thanks, but this does not really help me to understand how the energies work out in detail.

Here's are some example numbers for Pd:

When I put an isolated Pd atom in a large box (edge length 15A), I get a final energy value of -1.479 eV (Using the Pd_pv pseudopotential).

When I look at the "band structure" of this, the Fermi energy is given as

-4 eV in the DOSCAR. The DOS of course jumps at this energy because this is not really a band, but an isolated energy level.

When I create a Pd fcc crystal, the final energy I get is

-5.189eV/atom.

The Fermi energy in the DOSCAR is now given as

8.916 eV, but relative to what?

It is not the bottom of the band because that sits at about 2.2eV (integrated number of electrons is 6/atom (4p electrons) below this value and increases from there.)

So what sets the zero points of these energies? Why is the Fermi energy is 12eV lower in the isolated case?

Sorry for misunderstanding your question. In the case of single Pd that you mentioned, the Fermi level is basically the atomic energy of the topmost orbital which is 5s orbital which is -4.3 for the Pd_pv pseudopotential. You can find this information in the POTCAR as well. This is for example what you get for Pd_pv:

Atomic configuration

11 entries

n l j E occ.

1 0 0.50 -24164.6632 2.0000

2 0 0.50 -3530.9218 2.0000

2 1 1.50 -3165.6822 6.0000

3 0 0.50 -639.4578 2.0000

3 1 1.50 -519.4532 6.0000

3 2 2.50 -328.9758 10.0000

4 0 0.50 -88.7014 2.0000

4 1 1.50 -54.4541 6.0000

4 2 2.50 -6.6878 9.0000

5 0 0.50 -4.3543 1.0000

4 3 2.50 -1.3606 0.0000

The situation in a crystal is drastically different from a single atom. In this case the bonding-anti bonding interactions between atoms and their underlying Coulombic interactions can greatly change the electronic structure. It is thus not a surprise to see such a shift in the Fermi level.

Thanks, but I am still confused whether the energies in the DOSCAR can actually be compared directly or whether they only make sense after some energy shift.

For example, in the calculation with the Pd_pv-potential, even the energy of the 4p-electrons shifts by about 10 eV between the crystal structure and the isolated atom.

(In the crystal, the p-level starts to become occupied at -40.8eV, in the isolated atom at -50.9eV)

If I were to simply add up the energy contributions from the DOSCAR (multiply each energy value with the DOS at this point and the bin width and sum it all up), the energy of the electrons would be much lower for the isolated atom than for the crystal which makes not much sense to me.

Again, this difference comes form the coulomb interaction which is obviously different for an atom in a crystal as compared with that in vacuum. To clarify this, I just did two test calculations; one for isolated Pd (Gamma-only) and one for FCC-Pd (10x10x10 k-mesh) using Pd_pv-potential. Comparing their OUTCAR files, I get this:

For isolated Pd

E-fermi : -3.7225 XC(G=0): -0.8575 alpha+bet : -0.2438

For FCC-Pd

E-fermi : 8.8276 XC(G=0): -13.0534 alpha+bet :-15.7412

The term XC(G=0) is the G=0 component of Coulomb interaction excluded from the energy bands. Now if you add this value to E-fermi you get these:

For isolated Pd

E-fermi+XC(G=0)= -4.58 eV

For FCC-Pd

E-fermi+XC(G=0)= -4.22 eV

Which are pretty much the same. Also, the difference between the XC(G=0) is (I hope you agree) very close to the 10 eV difference, puzzling you.

Great, that was what I was looking for.

I did not realise that VASP reports the numbers this way.

Thanks a lot.

My pleasure. I'm glad I could be of help.

May I add another question?

Could you also briefly tell me what the alpha+beta - Term is?

In older versions of VASP; after the E-Fermi-line, it is stated:

"add alpha+bet to get absolute eigenvalues".

Or is there a detailed documentation of this output somewhere?

alpha+bet term is in-principle the reference energy point of vacuum (Zero energy), but in practice it only works for isolated atoms. The alpha part corresponds to the divergent (G=0) component of the electrostatic ion-ion interaction per electron. You can get this value by dividing the PSCENC term by NELECT (total number of electrons). The value of PSENCE can be found in the following line in the OUTCAR file:

alpha Z PSCENC = 2.03904736

As for the bet term, it corresponds to the compensating G=0 component of Pseudopotential. In the case of an isolated atom, the corresponding alhpha+bet can give a reliable reference energy for the vacuum. For a slab, this value becomes less and less reliable, depending on whether the slab is charged or not and how the vacuum is treated. A better way to find the vacuum reference energy in this case is to use LVTOT=True option. Finally, for a bulk this value has no real physical meaning.

Great, thanks again, that was super helpful