From 0252acbd0df7c7ff77bac3ef00632fdfe607dd0c Mon Sep 17 00:00:00 2001 From: William Dawson Date: Mon, 16 Dec 2024 14:45:10 +0900 Subject: [PATCH 1/2] 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 From 5d5d264238584ccac5d9f89b89f3da9a97cb313b Mon Sep 17 00:00:00 2001 From: William Dawson Date: Tue, 17 Dec 2024 13:00:58 +0900 Subject: [PATCH 2/2] Include our own ERF --- Source/Fortran/FermiOperatorModule.F90 | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 index 376a0321..ebe343a3 100644 --- a/Source/Fortran/FermiOperatorModule.F90 +++ b/Source/Fortran/FermiOperatorModule.F90 @@ -625,5 +625,24 @@ SUBROUTINE ComputeCStep(X, A, W, pool, threshold, Out) CALL DestructMatrix(XA) END SUBROUTINE ComputeCStep +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> ERF is a Fortran2008 feature, and we just need a crude approximation here + !! so we implement it with the Abramowitz and Stegun approximation + !! (credit ChatGPT). + FUNCTION ERF(x) + REAL(NTREAL) :: x, erf + REAL(NTREAL) :: t, z, tau + REAL(NTREAL), PARAMETER :: a1 = 0.254829592_NTREAL + REAL(NTREAL), PARAMETER :: a2 = -0.284496736_NTREAL + REAL(NTREAL), PARAMETER :: a3 = 1.421413741_NTREAL + REAL(NTREAL), PARAMETER :: a4 = -1.453152027_NTREAL + REAL(NTREAL), PARAMETER :: a5 = 1.061405429_NTREAL + REAL(NTREAL), PARAMETER :: p = 0.3275911_NTREAL + + z = ABS(x) + t = 1.0_NTREAL / (1.0_NTREAL + p * z) + tau = t * (a1 + t * (a2 + t * (a3 + t * (a4 + t * a5)))) + erf = SIGN(1.0_NTREAL, x) * (1.0_NTREAL - tau * EXP(-z * z)) + END FUNCTION ERF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE FermiOperatorModule