Introduce a guard against overflow in the dense FOE

william-dawson · Dec 16, 2024 · 0252acb · 0252acb
1 parent 20f55fe
commit 0252acb
Showing 1 changed file with 6 additions and 2 deletions.
diff --git 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