CdSe surface calculations

[glab]

Dear all,

I am a new VASP user. I am trying to perform CdSe surface gap calculations by the DFT method with a hybrid exchange-correlation functional. The hybrid functional has been chosen to fit the experimental gap value for the bulk. I have obtained very strange results for the surface calculations. I would expect the band gap values of the surface reaching the bulk band gap with increasing the number of layers. (Thus, if I have, let's say, 2.5 eV for 1 layer , then the gap value should be less for 2 layers, and so on until it reaches the bulk band gap.) However, my values for the gap oscillate in the wide range of energies (0.7-5.0 eV for 1-5 layers) with a number of layers if I keep ions unrelaxed. On the other hand, if I relax ions, cell shape and cell volume, the gap values go down how it is supposed to be but the gap value for 2 layers is already lower than the bulk gap value. I also tried GW calculations which give similar results. Here is my input:

INCAR

System = zinc-blende CdSe
ISMEAR = 0
SIGMA = 0.01
GGA = PE
#NSW = 100
#POTIM = 0.2
#IBRION = 2
#ISIF = 3
## Hybrid function
#LHFCALC = .TRUE. ; GGA = B3 ; AEXX = 0.2 ; AGGAX = 0.72
#AGGAC = 0.81 ; ALDAC = 0.19
#ALGO = D ; TIME = 0.4

First, I perform a standard PBE calculation and then HF+DFT.

The POSCAR file was created with the VESTA software. The distance between layers is about 35 A for all calculations.
For KPOINTS file I used 4x4x1 mesh for preliminary results.

What can be the problem? Do I do something wrong?

Thank you!

[zzhlax]

It is quite complicated for surface calculations I think. Maybe there are some surface oscillations overlapping with the bulk states, leading to a smaller band bap for the surface than for the bulk.
In addition, it is also important whether the surface is nonpolar or polar.

chears,