diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index fafeb1c3..376a0321 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -122,8 +122,12 @@ SUBROUTINE ComputeDenseFOE(H, ISQ, trace, K, inv_temp_in, & chemical_potential = left + (right - left) / 2 DO II = 1, num_eigs sval = eigs(II) - chemical_potential - ! occ(II) = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) - occ(II) = 1.0_NTREAL / (1.0_NTREAL + EXP(inv_temp * sval)) + ! Guard against overflow + IF (inv_temp * sval .GT. 30) THEN + occ(II) = 0.5_NTREAL * (1.0_NTREAL - ERF(inv_temp * sval)) + ELSE + occ(II) = 1.0_NTREAL / (1.0_NTREAL + EXP(inv_temp * sval)) + END IF END DO sv = SUM(occ) IF (ABS(trace - sv) .LT. 1E-8_NTREAL) THEN