Persistent Delocalized States in Vacuum Region for Bi Surface Calculation

[akuklin]

Hello VASP community,

I am encountering an unusual issue while calculating the ELFCAR of a Bi surface using the PBE functional. The results reveal a significant presence of delocalized states throughout the vacuum region. Delocalized electronic states appear distinctly in the vacuum region, not extending directly from the surface, but rather appearing independently within the vacuum. This behaviour suggests these states may be artefacts rather than physical features.
I also calculated LOCPOT which changes quite a lot in the vacuum region (see attachments)

Here are the specific details of my setup and observations:
PREC= high
GGA = PE
IVDW = 11
Electronic Relaxation 1
NELM = 80;
EDIFF = 1E-04 stopping-criterion for ELM
Ionic Relaxation c-
NSW = 0 number of steps for IOM
NBLOCK = 1;
IBRION = -1 ionic relax: 0-MD 1-quasi-New 2-CG
ISIF = 2 stress and relaxation
ISYM = 2
LCORR = T Harris-correction to forces
LDIAG = T
LCHARG = .T.
EDIFFG = -1E-02
ISPIN = 1
LWAVE =.T.
LVTOT = T
LELF = T
ISMEAR = -5;
IALGO = 38 algorithm
SIGMA = .05 broadening in eV -4-tet -1-fermi 0-gaus
LREAL = A real-space projection
NCORE = 4
LSCALU = .F.
LASPH = .TRUE.
LORBIT = 11
NEDOS = 5000
ENCUT = 500
A vacuum thickness was over 20 Å, ensuring minimal interaction between periodic images.

I believe that states could be a numerical artefact, perhaps related to the functional but not sure 100%. I attempted calculations using SCAN, which completely eliminated these vacuum states, suggesting they may indeed be artefacts of PBE. However, when using HSE06, the delocalized states persist, though they are less pronounced.
Would it be worthwhile to test a GW calculation for this system, or are there adjustments within DFT that could better address this issue and make me 100% sure that this is an artefact?
Any insights or suggestions would be highly appreciated, especially from those who have experience with surface calculations and dealing with spurious vacuum states in similar setups.

Thank you!

[henrique_miranda]

Could you share the INCAR, KPOINTS, OUTCAR and OSZICAR of the calculations with the different functionals: PBE, HSE06 and SCAN?
Calculations with a lot of vacuum tend to be hard to converge, before trying to answer your question we should make sure that the calculations converged properly.

[akuklin]

Please check the attachments

[henrique_miranda]

OK, thank you for sharing the OUTCAR and OSZICAR.
Your calculations seem to be converging well. At least, I don't see any obvious problem.
EDIFF is maybe a bit large, but I doubt lowering it will change the final results significantly.

Do you also see this delocalization of the states in the CHGCAR by means of a large charge density in the vaccum region?
Could you share the ELFCAR.jpg for the different functionals too?

[akuklin]

Yes, the states can also be seen in the CHGCAR when setting the appropriate isosurface level
The ELFCARs are attached