From 0252acbd0df7c7ff77bac3ef00632fdfe607dd0c Mon Sep 17 00:00:00 2001
From: William Dawson <william.dawson@riken.jp>
Date: Mon, 16 Dec 2024 14:45:10 +0900
Subject: [PATCH] Introduce a guard against overflow in the dense FOE

---
 Source/Fortran/FermiOperatorModule.F90 | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

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