Fermi-Dirac分布函数对变量E

编程问答 2022-06-13

问题描述

我编写了一个Fortran代码，以计算Fermi-Dirac分布函数对变量E的导数。Fermi-Dirac分布函数本身如下所述。

$$ 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）$$。

但是，在运行Fortran代码以计算此派生函数之后，我获得了'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,df)
WRITE (UNIT=*,FMT=*) 't=0.0d0 ','ef = 0.5 ','e=0.1 ',df

e = -5.0d-1
CALL FE(e,'e=-0.5 ',df

e = 5.0d-1
CALL FE(e,'e=0.5 ',df

e = 1.0d0
CALL FE(e,'e=1.0 ',df

e = 5.0d0
CALL FE(e,'e=5.0 ',df

t = 5.0d0

e = 1.0d-1
CALL FE(e,FMT=*) 't=5.0d0 ',df

STOP

END PROGRAM TEST

SUbroUTINE FE(e,df)
IMPLICIT NONE
DOUBLE PRECISION :: e,kb,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

有人可以给我一些有关如何解决问题的建议吗？预先谢谢你。

解决方法

好吧，我可以说只看代码就有两件事

  1. 代码中没有规定处理T = 0。当您通过T = 0时，您应该期望NaN，即除以0

  2. 看起来您使用的单位是以eV为单位的能量，以K为单位的温度。那么您的Boltzman常数将完全关闭。

k B = 1/11605 eV / K

尝试此代码，可能会更好

    SUBROUTINE FE(e,ef,T,df)
    IMPLICIT NONE

    DOUBLE PRECISION :: e,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
distributiondistributionfortran