fortran – Derivative of Fermi-Dirac 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.
kB = 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