费米-狄拉克分布函数对变量 E 的导数

Question

我写了一个 Fort运行代码来计算费米-狄拉克分布函数对变量 E 的导数。费米-狄拉克分布函数本身写在下面。

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

其中，$$E$$和$$T$$是变量； $$E_f$$ 和 $$k_B$$ 是常量

这个函数对变量的导数写在下面。

$$\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}$$

我想当 $$T$$ 为零时，这个导数函数应该充当增量函数 $$\delta (E-E_F)$$。

然而，我在运行我的 Fort运行代码后获得 'NAN' 符号来计算这个导数函数。这是我运行我的代码后的结果。

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

这是我的代码。

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

谁能给我一些解决问题的建议？提前谢谢你。

Answer 1

嗯，我可以说只用眼睛看代码有两点不对

代码中没有规定处理T=0。当你通过 T=0 你应该期待 NaN，你除以 0
看起来你使用的单位是能量（eV），温度（K）。那么你的玻尔兹曼常数就完全不对了。

k_B = 1/11605 eV/K

试试这个代码，可能效果更好

    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

费米-狄拉克分布函数对变量 E 的导数

Derivative of Fermi-Dirac Distribution Function to Variable E

fortran

distribution