Ok. Your integrand has indeed a singularity at the origin. See below. Are you sure that your model / equations are correct?
m = 0.02
J = 50
muB = 0.67
Dst = 7.1
U = 313
g = 2.06
mu = muB*g*Hext - Dst
eps[kx_, ky_, kz_] := J*(3 - Cos[kx] - Cos[ky] - Cos[kz])
Ek[Delta_, kx_, ky_, kz_] :=
eps[kx, ky, kz]^(1/2)*(eps[kx, ky, kz] + Delta)^(1/2)
eps[0, 0, 0]
Ek[del, 0, 0, 0]
So the expression containing Ek will tend to infinitiy whern Ek is evaulated at the origin
fac = .2;
Plot3D[((eps[kx, ky, 0] + 1)/Ek[2, kx, ky, 0] - 1) Boole[ kx^2 + ky^2 + 0 > .0004], {kx, -Pi fac, Pi fac}, {ky, -Pi fac,
Pi fac}, PlotRange -> All, PlotPoints -> 20]
Plot[((eps[kx, 0, 0] + 1)/Ek[2, kx, 0, 0] - 1), {kx, -Pi fac, Pi fac},
PlotRange -> All, PlotStyle -> {Thick, Red}, AxesOrigin -> {0, 0}]