Theoretically, to obtain the band edges one should calculate the integrated density of electronic states D(ε). Given the bandstructure {εl(k) | l, k}, where l denotes the band index and k a point inside the first Brillouin zone, 1BZ, one has

(1) D(ε) = Ω-1 Σk Σl δ(ε - εl(k)),

where Ω denotes the volume of the system, Σk sum over all points of the 1BZ and Σl sum over all band indices. Using the equality

(2) Ω-1 Σk f(k) = (2π)-d ∫1BZ ddk f(k)

for microscopic systems, where d denotes the dimensionality of the system (d=3 for three-dimensional systems), one calculates D(ε) by sampling the integral on the right-hand side of Eq. (2), using for instance a sufficiently large set of special points. To smoothen the function D(ε) thus calculated, one may adopt the following expression for the δ function:

(3) δ(x) = (η/π)/(x2 + η2),  η ↓ 0,

where instead of using an infinitesimal positive value for η, in practice one should use a finite positive value. This value should not be too large, since ultimately any non-infinitesimal value for this constant leads to disappearance of the gap. □

I should point out that for the calculation of the band edges one does not need explicitly to calculate the function D(ε). What is needed is calculation of the value of the chemical potential μ, separating the valence-band energies and the conduction-band energies. Given the number of particles per unit cell Nc, the chemical potential is determined from the equality

(4)  Ωc ∫-∞μ dε D(ε) = Nc,

where Ωc is the unit-cell volume. From Eqs (1), (2) and (4) one obtains

(5)  (2π)-d ∫1BZ ddk θ(μ-εl(k)) = Nc / Ωc,

where θ(x) is the unit-step function. In practice one replaces this function by the Fermi-Dirac distribution function fT(μ-εl(k)) for sufficiently small value of the temperature T. This has a similar (but not identical) effect as using a non-infinitesimal value for η in Eq. (3). Solving Eq. (5) for μ after discretising the integral over the 1BZ, one can calculate the band edges as follows.

Given the set S = {εl(k) | l, k} generated for the calculation of μ, one has to separate this set into two sets S1 and S2, with S1 containing all elements of S less than μ, and S2 all elements of S greater than μ. The sought-after band edges coincide with the value of the largest element of S1 and smallest element of S2. □

The Haydock-Heine-Kelly Recursion method provides a great way of obtaining the band edges. Please read Springer-Verlag "Solid State Physics" series of monographs Vol 35, spacially the article by Haydock.  Then go to the Homepage of Haydock in University of Oregon. The "Cambridge Recursion Library"  can be downloaded free of charge.  The relevant subroutines to use are  BNDCRD (crude band edges), BNDEST (better estimate of band edges) and BNDREF (refined band edges).