Fix bug in fermi function: Floating-point exception in debug mode

Author: ares201005

src/prg_densitymatrix_mod.F90

@@ -628,8 +628,12 @@ real(dp) function fermi(e,ef,kbt)
 
     real(dp), intent(in) :: e, ef, kbt
 
-    fermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
-
+    if ((e-ef)/kbt > 100.0_dp) then
+      fermi = 0.0_dp
+    else
+      fermi = 1.0_dp/(1.0_dp+exp((e-ef)/(kbt)))
+    endif
+   
   end function fermi
 
 end module prg_densitymatrix_mod