This question is better suited for the Computational Science site. Still here are some points for you to think about.
First, the range of integration is the intersection of
(-oo, -eV-gap) U (-eV+gap, +oo) and
(-oo, -gap) U (gap, +oo). There are two possible cases:
eV < 2*gap then the allowed energy values are in
(-oo, -eV-gap) U (gap, +oo);
eV > 2*gap then the allowed energy values are in
(-oo, -eV-gap) U (-eV+gap, -gap) U (gap, +oo).
Second, you are working in a very low temperature region. With
t equal to 0.02 K, the denominator in the Boltzmann factor is 1.7 ÂµeV, while the energy gap is 400 ÂµeV. In this case the value of the exponent is huge for positive energies and it soon goes off the limits of the double precision floating point numbers, used by Python. As this is the minimum possible positive energy, things would not get any better at higher energies. With negative energies the value would always be very close to zero. Note that at this temperature, the Fermi-Dirac distribution has a very sharp edge and resembles a reflected theta function. At
E = gap you would have
exp(E/kT) of approximately 6.24E+100. You would run out of range when
E/kT > 709.78 or
E > 3.06*gap.
Yet it makes no sense to go to such energies since at that temperature the difference between the two Fermi functions very quickly becomes zero outside the
[-eV, 0] interval which falls entirely inside the gap for the given temperature when
V < (2*gap)/e (0.8 mV). That's why one would expect that the current would be very close to zero when the bias voltage is less than 0.8 mV. When it is more than 0.8 mV, then the main value of the integral would come from the integrand in
(-eV+gap, -gap), although some non-zero value would come from the region near the singularity at
E = gap and some from the region near the singularity at
E = -eV-gap. You should not avoid the singularities in the DoS, otherwise you would not get the expected discontinuities (vertical lines) in the I(V) curve (image taken from Wikipedia):
Rather, you have to derive equivalent approximate expressions in the vicinity of each singularity and integrate them instead.
As you can see, there are many special cases for the value of the integrand and you have to take them all into account when computing numerically. If you don't want to do that, you should probably turn to some other mathematical package like Maple or Mathematica. These have much more sophisticated numerical integration routines and might be able to directly handle your formula.
Note that this is not an attempt to answer your question but rather a very long comment that would not fit in any comment field.