fortran – Derivative of FermiDirac Distribution Function to Variable E
fortran – Derivative of FermiDirac Distribution Function to Variable E
Well, I could say there are two thing wrong with just eye viewing the code

There is no provision to handle T=0 in the code. When you pass T=0 you should expect NaN, you divide by 0

Looks like units you use are energy in eV, temperature in K. Then your Boltzman constant is completely off.
k_{B} = 1/11605 eV/K
Try this code, might work better
SUBROUTINE FE(e, ef, T, df)
IMPLICIT NONE
DOUBLE PRECISION :: e, ef, T, df
DOUBLE PRECISION :: kB, q
kB = 1.0d0/11605.0d0 ! eV/K
df = 0.0d0
if (T .ne. 0.0d0) then ! 0 at T=0 ?
q = (e  ef)/(kB*T)
q = EXP( q ) + 1.0d0
df = (q  1.0d0) / q / q
df = df / (kB*T)
else ! T = 0
if (e .eq. ef) then
df = 1.0d0
endif
endif
RETURN
END SUBROUTINE FE