Derivative of Fermi-Dirac Distribution Function to Variable E

Question

I wrote a Fortran code to compute the derivative of Fermi-Dirac distribution function to variable E. The Fermi-Dirac distribution function itself is written below.

$$f_{(E)} =\frac{1}{e^{\frac{E-E_f}{k_BT}}+1}$$

where, $$E$$ and $$T$$ are variables; $$E_f$$ and $$k_B$$ are constant

The derivative of this function to variable is written below.

$$\frac{\partial f_{(E)}}{\partial E} =\frac{1}{k_BT}\frac{e^{\frac{E-E_f}{k_BT}}}{(e^{\frac{E-E_f}{k_BT}}+1)^2}$$

I suppose this derivative function should act as a delta function $$\delta (E-E_F)$$, when $$T$$ is zero.

However, I obtain 'NAN' sign after I ran my Fortran code to compute this derivative function. Here are the results after I ran my code.

t=0.0d0 ef = 0.5 e=-0.5                      NaN
t=0.0d0 ef = 0.5 e=0.5                      NaN
t=0.0d0 ef = 0.5 e=1.0                      NaN
t=0.0d0 ef = 0.5 e=5.0                      NaN
t=5.0d0 ef = 0.5 e=0.1   0.000000000000000E+000
t=5.0d0 ef = 0.5 e=-0.5   0.000000000000000E+000
t=5.0d0 ef = 0.5 e=0.5  -3.621485258019960E+021
t=5.0d0 ef = 0.5 e=1.0                      NaN
t=5.0d0 ef = 0.5 e=5.0                      NaN

Here are my code.

PROGRAM TEST
IMPLICIT NONE
DOUBLE PRECISION :: e, ef, t,df

t = 0.0d0
ef = 5.0d-1

e = 1.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=0.0d0 ', 'ef = 0.5 ', 'e=0.1 ', df

e = -5.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=0.0d0 ', 'ef = 0.5 ', 'e=-0.5 ', df

e = 5.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=0.0d0 ', 'ef = 0.5 ', 'e=0.5 ', df

e = 1.0d0
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=0.0d0 ', 'ef = 0.5 ', 'e=1.0 ', df

e = 5.0d0
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=0.0d0 ', 'ef = 0.5 ', 'e=5.0 ', df

t = 5.0d0

e = 1.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=5.0d0 ', 'ef = 0.5 ', 'e=0.1 ', df

e = -5.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=5.0d0 ', 'ef = 0.5 ', 'e=-0.5 ', df

e = 5.0d-1
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=5.0d0 ', 'ef = 0.5 ', 'e=0.5 ', df

e = 1.0d0
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=5.0d0 ', 'ef = 0.5 ', 'e=1.0 ', df

e = 5.0d0
CALL FE(e,ef,t,df)
WRITE (UNIT=*, FMT=*) 't=5.0d0 ', 'ef = 0.5 ', 'e=5.0 ', df

STOP

END PROGRAM TEST

SUBROUTINE FE(e,ef,t,df)
IMPLICIT NONE
DOUBLE PRECISION :: e, ef, kb, t,df

kb = 1.380649d-23

df = 0.0d0

df = 1.0d0 / kb / t * EXP(1.0d0 / kb / t * (e - ef))
df = df / (EXP(1.0d0 / kb / t * (e - ef)) + 1.0d0) ** 2
df = df * -1.0d0

RETURN
END SUBROUTINE FE

Would anyone please give me some suggestions on how to sort out the problem? Thank you in advance.

Answer 1

Well, I could say there are two thing wrong with just eye viewing the code

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

2. 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