Laplace correction was used for solutions of Poisson equation in density functional theory (DFT) slab simulations of GaN (0001) surface. In accordance to the scenario of Meyer and Marx [Phys. Rev. B 67, 35403 (2003)], the electric field, emerging within the slab, is determined by the condition that Fermi energy is equal to the energies of valence band maximum and conduction band minimum at both sides of the slab. In contrast to their predictions that the field is identical for different slab thicknesses, it was found from DFT calculations that the electric potential difference between two sides of the slab is kept constant. Therefore, the energy of the slab does not diverge for large slab thickness. It was also found that, in most cases, the field in the slab is determined by pinning by surface states. Thus the change in opposite side slab termination, both by different types of atoms and by their location, can be used to change electric field in the slab, creating a tool that can be used to simulate the change in surface properties due to doping of the bulk. It was shown that, depending on the electric field, the energy of surface states changes in the way different from the band states. This change could be large, comparable to the bandgap. Despite the fact that the Fermi energy can be pinned to surface states, it may change its energy with respect to valence band and conduction band as a function of the doping of bulk semiconductor. These results are in agreement with the experimental data showing variation in Fermi energy at the GaN (0001) surface in function of the doping in the bulk. It was also shown using the electron density and potential profiles that ten Ga–N atomic layers are necessary to achieve small quantum overlap of the surface states of both sides of the slab, which is required for high precision simulations of the surface.