From c8cf3421b712609bcb09bff63bc2364955eb84c3 Mon Sep 17 00:00:00 2001
From: Chen Chengbing <1747193328@qq.com>
Date: Mon, 12 Oct 2026 12:52:26 +0800
Subject: [PATCH 1/5] Fix smeared-occupation Fermi convergence and
 particle-number accuracy

---
 source/source_estate/occupy.cpp               | 194 ++++++++------
 .../test/elecstate_occupy_test.cpp            | 244 +++++++++++++++++-
 2 files changed, 352 insertions(+), 86 deletions(-)

diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp
index 8fd2898b9c2..8d6b1aebc77 100644
--- a/source/source_estate/occupy.cpp
+++ b/source/source_estate/occupy.cpp
@@ -1,5 +1,8 @@
 #include "occupy.h"
 
+#include <limits>
+#include <sstream>
+
 #include "source_base/constants.h"
 #include "source_base/mymath.h"
 #include "source_base/parallel_reduce.h"
@@ -235,7 +238,7 @@ void Occupy::gweights(const int nks,
             continue;
 }
 
-        for (int ib = 0; ib < PARAM.globalv.nbands_l; ib++)
+        for (int ib = 0; ib < nband; ib++)
         {
             //================================
             // Calculate the gaussian weights
@@ -281,104 +284,126 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
                      const int& is,
                      const std::vector<int>& isk)
 {
-    // ModuleBase::TITLE("Occupy","efermig");
-    //==================================================================
-    //  Finds the Fermi energy - Gaussian Broadening (Methfessel-Paxton)
-    //==================================================================
-    // find bounds for the Fermi energy. Very safe choice!
-    const int maxiter = 300;
-    const double eps = 1.0e-10;
-
-    /*
-    for(int ik=0; ik<nks; ik++)
+    // Resolve particle number to floating-point accuracy, independently of the
+    // outer SCF threshold. A fixed absolute tolerance leaves a charge floor.
+    const double electron_magnitude = std::abs(nelec);
+    const double scale = std::max(1.0, electron_magnitude);
+    const double tolerance = 8.0 * std::numeric_limits<double>::epsilon() * scale;
+    if (!(smearing_sigma > 0.0) || !std::isfinite(smearing_sigma) || !std::isfinite(nelec))
+    {
+        ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Invalid smearing width or electron number");
+    }
+
+    // Empty local partitions contribute neutral extrema to the global bounds.
+    double lower = std::numeric_limits<double>::infinity();
+    double upper = -std::numeric_limits<double>::infinity();
+    for (int ik = 0; ik < nks; ++ik)
     {
-        for(int i=0; i<nband; i++)
+        if (is != -1 && is != isk[ik])
         {
-            std::cout << " ib=" << i << " ekb=" << ekb[ik][i] << std::endl;
+            continue;
+        }
+        for (int ib = 0; ib < nband; ++ib)
+        {
+            const double eigenvalue = ekb(ik, ib);
+            lower = std::min(lower, eigenvalue);
+            upper = std::max(upper, eigenvalue);
         }
     }
-    */
-
-    // the first 0 stands for the first k point.
-    double elw = ekb(0, 0);
-    double eup = ekb(0, nband - 1);
-
-    for (int ik = 1; ik < nks; ik++) // do ik = 2, nks
+#ifdef __MPI
+    Parallel_Reduce::reduce_min(lower);
+    Parallel_Reduce::reduce_max(upper);
+#endif
+    if (!std::isfinite(lower) || !std::isfinite(upper))
     {
-        elw = std::min(elw, ekb(ik, 0));
-        eup = std::max(eup, ekb(ik, nband - 1));
+        ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "No finite eigenvalue bounds");
     }
 
-    eup += 2 * smearing_sigma;
-    elw -= 2 * smearing_sigma;
-    // find min and max across pools
-    #ifdef __MPI
-    Parallel_Reduce::reduce_max(eup);
-    Parallel_Reduce::reduce_min(elw);
-    #endif
-    //=================
-    // Bisection method
-    //=================
-    // call sumkg
-    int changetime = 0;
-    while (true)
+    double step = 2.0 * smearing_sigma;
+    lower -= step;
+    upper += step;
+    double low_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, lower, is, isk);
+    double high_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, upper, is, isk);
+    int expansion = 0;
+    while (low_count > nelec + tolerance || high_count < nelec - tolerance)
     {
-        const double sumkup = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, eup, is, isk);
-        const double sumklw = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, elw, is, isk);
-
-        if (changetime > 1000)
-        {
-            std::cout << " SOMETHING WRONG: " << std::endl;
-            std::cout << " is = " << is << std::endl;
-            std::cout << " eup = " << eup << std::endl;
-            std::cout << " elw = " << elw << std::endl;
-            std::cout << " nband = " << nband << std::endl;
-            std::cout << " nelec = " << nelec << std::endl;
-            std::cout << " sumkup = " << sumkup << std::endl;
-            std::cout << " sumklw = " << sumklw << std::endl;
-            std::cout << " sumkup - nelec = " << sumkup - nelec << std::endl;
-            std::cout << " sumklw - nelec = " << sumklw - nelec << std::endl;
-            ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "ERROS in SMEARING");
-            // no need to break; quit directly
-        }
-        else if ((sumkup - nelec) < -eps)
+        if (++expansion > 100 || !std::isfinite(step))
         {
-            eup += 2 * smearing_sigma;
-            ++changetime;
+            ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Cannot bracket the requested electron number");
         }
-        else if ((sumklw - nelec) > eps)
+        step *= 2.0;
+        if (low_count > nelec + tolerance)
         {
-            elw -= 2 * smearing_sigma;
-            ++changetime;
+            lower -= step;
+            low_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, lower, is, isk);
         }
-        else
+        if (high_count < nelec - tolerance)
         {
-            break;
+            upper += step;
+            high_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, upper, is, isk);
         }
     }
 
-    for (int i = 0; i < maxiter; i++)
+    if (!std::isfinite(low_count) || !std::isfinite(high_count)
+        || !std::isfinite(lower) || !std::isfinite(upper))
     {
-        //======================
-        // change ef value
-        //======================
-        ef = (eup + elw) / 2.0;
-        const double sumkmid = occupy_smearing::sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, ef, is, isk);
+        ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Non-finite Fermi bracket");
+    }
 
-        if (std::abs(sumkmid - nelec) < eps)
+    // Keep the best evaluated point, including the saturated endpoints. For
+    // empty/full bands or a resolved gap, the chemical potential is not unique.
+    ef = lower;
+    double best_error = std::abs(low_count - nelec);
+    const double high_error = std::abs(high_count - nelec);
+    if (high_error < best_error)
+    {
+        ef = upper;
+        best_error = high_error;
+    }
+    for (int iteration = 0; iteration < 300; ++iteration)
+    {
+        if (best_error <= tolerance)
         {
             return;
         }
-        else if ((sumkmid - nelec) < -eps)
+        const double middle = 0.5 * lower + 0.5 * upper;
+        if (middle == lower || middle == upper)
         {
-            elw = ef;
+            // No representable interior point remains. Do not silently spin to
+            // maxiter, nor claim that a small energy interval proves charge accuracy.
+            std::ostringstream message;
+            message.precision(17);
+            message << "Fermi bracket reached floating-point resolution; electron residual = "
+                    << best_error << ", requested tolerance = " << tolerance;
+            const std::string diagnostic = message.str();
+            ModuleBase::WARNING_QUIT("occupy_smearing::efermig", diagnostic);
+            return;
+        }
+        const double count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, middle, is, isk);
+        if (!std::isfinite(count))
+        {
+            ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Non-finite electron number");
+        }
+        const double error = std::abs(count - nelec);
+        if (error < best_error)
+        {
+            ef = middle;
+            best_error = error;
+        }
+        if (best_error <= tolerance)
+        {
+            return;
+        }
+        if (count < nelec)
+        {
+            lower = middle;
         }
         else
         {
-            eup = ef;
+            upper = middle;
         }
     }
-    return;
+    ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Fermi search exceeded its iteration limit");
 }
 
 /**
@@ -407,21 +432,24 @@ double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
 {
     // ModuleBase::TITLE("Occupy","sumkg");
     double sum2 = 0.0;
-    for (int ik = 0; ik < nks; ik++)
+    double compensation = 0.0;
+    for (int ik = 0; ik < nks; ++ik)
     {
-        if (is != -1 && is != isk[ik]) {
+        if (is != -1 && is != isk[ik])
+        {
             continue;
-}
-
-        double sum1 = 0.0;
-        for (int ib = 0; ib < nband; ib++)
+        }
+        for (int ib = 0; ib < nband; ++ib)
         {
-            //===========================
-            // call wgauss
-            //===========================
-            sum1 += occupy_smearing::wgauss((e - ekb(ik, ib)) / smearing_sigma, ngauss);
+            // Sum the same weighted occupations as gweights, compensating the
+            // local accumulation before the existing pool/band reduction.
+            const double argument = (e - ekb(ik, ib)) / smearing_sigma;
+            const double occupation = wk[ik] * wgauss(argument, ngauss);
+            const double corrected = occupation - compensation;
+            const double next = sum2 + corrected;
+            compensation = (next - sum2) - corrected;
+            sum2 = next;
         }
-        sum2 += wk[ik] * sum1;
     }
 
     // GlobalV::ofs_running << "\n sum2 before reduce = " << sum2 << std::endl;
diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/elecstate_occupy_test.cpp
index fdcf6ac46cf..34d9efc8571 100644
--- a/source/source_estate/test/elecstate_occupy_test.cpp
+++ b/source/source_estate/test/elecstate_occupy_test.cpp
@@ -1,4 +1,7 @@
 #include <string>
+#include <cmath>
+#include <limits>
+#include <sstream>
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 #include "source_io/module_parameter/parameter.h"
@@ -268,7 +271,11 @@ TEST_F(OccupyTest, Efermig)
   std::vector<int> isk = {0, 0};
   double ef = 0.0;
   occupy_smearing::efermig(ekb, 1, 1, 1.0, wk, smearing_sigma, ngauss, ef, is, isk);
-  EXPECT_NEAR(ef, -0.5, 1e-13);
+  // Full occupation fixes particle number, not a unique finite chemical potential.
+  EXPECT_TRUE(std::isfinite(ef));
+  EXPECT_GT(ef, ekb(0, 0));
+  const double count = occupy_smearing::sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, ef, is, isk);
+  EXPECT_NEAR(count, 1.0, 2e-15);
 }
 
 TEST_F(OccupyTest, Gweights)
@@ -283,10 +290,241 @@ TEST_F(OccupyTest, Gweights)
   std::vector<int> isk = {0, 0};
   double ef = 0.0;
   ModuleBase::matrix wg(1, 1);
-  wg(0, 0) = 1.0;
+  wg(0, 0) = -123.0;
   double demet = 0.0;
   occupy.gweights(1, wk, 1, 1.0, smearing_sigma, ngauss, ekb, ef, demet, wg, is, isk);
-  EXPECT_NEAR(ef, -0.5, 1e-13);
+  // Full occupation fixes particle number, not a unique finite chemical potential.
+  EXPECT_TRUE(std::isfinite(ef));
+  EXPECT_GT(ef, ekb(0, 0));
+  const double count = occupy_smearing::sumkg(ekb, 1, 1, wk, smearing_sigma, ngauss, ef, is, isk);
+  EXPECT_NEAR(count, 1.0, 2e-15);
   EXPECT_NEAR(demet, 0.0, 1e-13);
   EXPECT_NEAR(wg(0, 0), 1.0, 1e-13);
 }
+
+TEST_F(OccupyTest, FractionalFermiDiracOccupation)
+{
+    ModuleBase::matrix energies(1, 1);
+    energies(0, 0) = -0.7;
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    const double sigma = 0.13;
+    const double target = 0.3141592653589793;
+    double ef = 0.0;
+    occupy_smearing::efermig(energies, 1, 1, target, weights, sigma, -99, ef, -1, spins);
+    const double expected = energies(0, 0) + sigma * std::log(target / (1.0 - target));
+    EXPECT_NEAR(ef, expected, 2e-15);
+    ModuleBase::matrix occupations(1, 1);
+    occupations(0, 0) = -123.0;
+    double entropy = 0.0;
+    occupy.gweights(1, weights, 1, target, sigma, -99, energies, ef, entropy, occupations, -1, spins);
+    EXPECT_NEAR(occupations(0, 0), target, 2e-15);
+    const double expected_entropy = sigma * (target * std::log(target) + (1.0 - target) * std::log(1.0 - target));
+    EXPECT_NEAR(entropy, expected_entropy, 2e-15);
+}
+
+TEST_F(OccupyTest, WeightedParticleNumberAllSmearings)
+{
+    ModuleBase::matrix energies(3, 4);
+    const std::vector<double> weights{0.17, 0.29, 0.54};
+    const std::vector<int> spins{0, 0, 0};
+    for (int ik = 0; ik < 3; ++ik)
+    {
+        for (int ib = 0; ib < 4; ++ib)
+        {
+            energies(ik, ib) = -0.4 + 0.27 * ib + 0.07 * ik;
+        }
+    }
+    const double target = 1.23456789;
+    for (const int smearing : {0, 1, 2, -1, -99})
+    {
+        double ef = 0.0;
+        double entropy = 0.0;
+        ModuleBase::matrix occupations(3, 4);
+        occupy.gweights(3, weights, 4, target, 0.1, smearing, energies, ef, entropy, occupations, -1, spins);
+        long double count = 0.0L;
+        for (int ik = 0; ik < 3; ++ik)
+        {
+            for (int ib = 0; ib < 4; ++ib)
+            {
+                count += occupations(ik, ib);
+            }
+        }
+        EXPECT_NEAR(count, target, 3e-15) << smearing;
+        EXPECT_TRUE(std::isfinite(entropy));
+    }
+}
+
+TEST_F(OccupyTest, EmptyFullAndGappedSpectra)
+{
+    ModuleBase::matrix energies(1, 2);
+    energies(0, 0) = -1.0;
+    energies(0, 1) = 1.0;
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    for (const double target : {0.0, 1.0, 2.0})
+    {
+        double ef = 0.0;
+        double entropy = 0.0;
+        ModuleBase::matrix occupations(1, 2);
+        occupy.gweights(1, weights, 2, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+        const double count = occupations(0, 0) + occupations(0, 1);
+        EXPECT_NEAR(count, target, 4e-15);
+        EXPECT_NEAR(entropy, 0.0, 1e-14);
+        if (target == 1.0)
+        {
+            EXPECT_GT(ef, energies(0, 0));
+            EXPECT_LT(ef, energies(0, 1));
+            EXPECT_NEAR(occupations(0, 0), 1.0, 2e-15);
+            EXPECT_NEAR(occupations(0, 1), 0.0, 2e-15);
+        }
+    }
+}
+
+TEST_F(OccupyTest, SelectedSpinAndUnsortedBands)
+{
+    ModuleBase::matrix energies(2, 3);
+    energies(0, 0) = 100.0;
+    energies(0, 1) = 200.0;
+    energies(0, 2) = 300.0;
+    energies(1, 0) = 0.4;
+    energies(1, 1) = -0.3;
+    energies(1, 2) = 0.1;
+    const std::vector<double> weights{1.0, 1.0};
+    const std::vector<int> spins{0, 1};
+    ModuleBase::matrix occupations(2, 3);
+    for (int ib = 0; ib < 3; ++ib)
+    {
+        occupations(0, ib) = -123.0;
+    }
+    double ef = 0.0;
+    double entropy = 0.0;
+    occupy.gweights(2, weights, 3, 1.314159, 0.08, 0, energies, ef, entropy, occupations, 1, spins);
+    double count = 0.0;
+    for (int ib = 0; ib < 3; ++ib)
+    {
+        EXPECT_EQ(occupations(0, ib), -123.0);
+        count += occupations(1, ib);
+    }
+    EXPECT_NEAR(count, 1.314159, 3e-15);
+}
+
+TEST_F(OccupyTest, CompensatedWeightedOccupationSum)
+{
+    const int points = 4096;
+    ModuleBase::matrix energies(points, 2);
+    std::vector<double> weights(points, 1.0 / points);
+    const std::vector<int> spins(points, 0);
+    long double expected = 0.0L;
+    for (int ik = 0; ik < points; ++ik)
+    {
+        // Repeated degenerate levels expose drift in ordinary accumulation.
+        energies(ik, 0) = -0.1;
+        energies(ik, 1) = 0.2;
+        for (int ib = 0; ib < 2; ++ib)
+        {
+            const double argument = (0.123 - energies(ik, ib)) / 0.1;
+            expected += weights[ik] * 0.5 * std::erfc(-argument);
+        }
+    }
+    const double count = occupy_smearing::sumkg(energies, 2, points, weights, 0.1, 0, 0.123, -1, spins);
+    EXPECT_NEAR(count, expected, 5e-16);
+}
+
+TEST_F(OccupyTest, FermiSearchAtRepresentableEnergyLimit)
+{
+    ModuleBase::matrix energies(1, 1);
+    energies(0, 0) = 1e8;
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    double ef = 0.0;
+    // The width is smaller than one ULP of the eigenvalue. No double-valued
+    // chemical potential can yield the requested fractional occupation.
+    testing::internal::CaptureStdout();
+    EXPECT_EXIT(occupy_smearing::efermig(energies, 1, 1, 0.3, weights, 1e-8, 0, ef, -1, spins),
+                ::testing::ExitedWithCode(1), "");
+    const std::string diagnostic = testing::internal::GetCapturedStdout();
+    EXPECT_THAT(diagnostic, testing::HasSubstr("Fermi bracket reached floating-point resolution"));
+    EXPECT_THAT(diagnostic, testing::HasSubstr("electron residual"));
+}
+
+TEST_F(OccupyTest, RejectInvalidWidthAndUnbracketedCharge)
+{
+    ModuleBase::matrix energies(1, 1);
+    energies(0, 0) = 0.0;
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    double ef = 0.0;
+    EXPECT_EXIT(occupy_smearing::efermig(energies, 1, 1, 0.5, weights, 0.0, 0, ef, -1, spins),
+                ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(occupy_smearing::efermig(energies, 1, 1, 1.1, weights, 0.1, 0, ef, -1, spins),
+                ::testing::ExitedWithCode(1), "");
+}
+
+TEST_F(OccupyTest, FermiParticleNumberAcrossScales)
+{
+    ModuleBase::matrix energies(3, 64);
+    const std::vector<double> weights{0.17, 0.29, 0.54};
+    const std::vector<int> spins{0, 0, 0};
+    for (int ik = 0; ik < 3; ++ik)
+    {
+        for (int ib = 0; ib < 64; ++ib)
+        {
+            energies(ik, ib) = -1.0 + 0.05 * ib + 0.013 * ik;
+        }
+    }
+    // Exercise efermig directly, independently of gweights' band-count fix.
+    for (const double target : {0.3141592653589793, 33.0, 63.125})
+    {
+        for (const int smearing : {0, 1, 2, -1, -99})
+        {
+            double ef = 0.0;
+            occupy_smearing::efermig(energies, 64, 3, target, weights, 0.1, smearing, ef, -1, spins);
+            long double count = 0.0L;
+            for (int ik = 0; ik < 3; ++ik)
+            {
+                for (int ib = 0; ib < 64; ++ib)
+                {
+                    const double argument = (ef - energies(ik, ib)) / 0.1;
+                    const double occupation = weights[ik] * occupy_smearing::wgauss(argument, smearing);
+                    count += occupation;
+                }
+            }
+            const double scale = std::max(1.0, target);
+            // Allow the independently accumulated result a small rounding margin.
+            const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * scale;
+            EXPECT_NEAR(count, target, tolerance) << "N = " << target << ", smearing = " << smearing;
+            if (target == 33.0 && smearing == 0)
+            {
+                const double residual = std::abs(static_cast<double>(count - target));
+                std::ostringstream value;
+                value.precision(17);
+                value << residual;
+                RecordProperty("N33_gaussian_residual", value.str());
+            }
+        }
+    }
+}
+
+TEST_F(OccupyTest, GweightsHonorsExplicitBandRange)
+{
+    ModuleBase::matrix energies(1, 3);
+    energies(0, 0) = -0.7;
+    energies(0, 1) = 0.1;
+    energies(0, 2) = 0.8;
+    ModuleBase::matrix occupations(1, 3);
+    for (int ib = 0; ib < 3; ++ib)
+    {
+        occupations(0, ib) = -123.0;
+    }
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    const double target = 0.3141592653589793;
+    double ef = 0.0;
+    double entropy = 0.0;
+    // The explicit one-band range must govern both root finding and writing.
+    occupy.gweights(1, weights, 1, target, 0.13, -99, energies, ef, entropy, occupations, -1, spins);
+    EXPECT_NEAR(occupations(0, 0), target, 2e-15);
+    EXPECT_DOUBLE_EQ(occupations(0, 1), -123.0);
+    EXPECT_DOUBLE_EQ(occupations(0, 2), -123.0);
+}

From cded5f4248aa1875a586c9100fee28f41d700ed5 Mon Sep 17 00:00:00 2001
From: Chen Chengbing <1747193328@qq.com>
Date: Mon, 12 Oct 2026 13:19:19 +0800
Subject: [PATCH 2/5] Preserve centered Fermi coordinates when evaluating
 occupations

---
 source/source_estate/occupy.cpp               | 161 +++++++++++++-----
 source/source_estate/occupy.h                 |   5 +-
 .../test/elecstate_occupy_test.cpp            | 115 ++++++++++++-
 3 files changed, 235 insertions(+), 46 deletions(-)

diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp
index 8d6b1aebc77..f00af5664d5 100644
--- a/source/source_estate/occupy.cpp
+++ b/source/source_estate/occupy.cpp
@@ -8,6 +8,38 @@
 #include "source_base/parallel_reduce.h"
 #include "source_io/module_parameter/parameter.h"
 
+namespace
+{
+// A chemical potential with an explicit energy origin. Occupations must use
+// both parts; the absolute double is only the rounded public/output value.
+struct FermiLevel
+{
+    double reference;
+    double offset;
+};
+
+FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
+                           const int nband,
+                           const int nks,
+                           const double& nelec,
+                           const std::vector<double>& wk,
+                           const double& smearing_sigma,
+                           const int ngauss,
+                           const int& is,
+                           const std::vector<int>& isk);
+
+double sum_occupations(const ModuleBase::matrix& ekb,
+                       const int nband,
+                       const int nks,
+                       const std::vector<double>& wk,
+                       const double& smearing_sigma,
+                       const int ngauss,
+                       const double reference,
+                       const double offset,
+                       const int& is,
+                       const std::vector<int>& isk);
+} // namespace
+
 Occupy::Occupy()
 {
 }
@@ -204,7 +236,7 @@ void Occupy::iweights(
  * @param smearing_sigma parameter input by user.
  * @param ngauss which type of smearing.
  * @param ekb band energy.
- * @param ef output: fermi level
+ * @param ef output: rounded absolute Fermi level
  * @param demet output: energy correction for metal
  * @param wg output: weight of each band at each k point
  * @param is spin
@@ -227,8 +259,9 @@ void Occupy::gweights(const int nks,
     //===============================
     //  Calculate the Fermi energy ef
     //===============================
-    //  call efermig
-    occupy_smearing::efermig(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, ef, is, isk);
+    // Keep the centered solution when constructing the occupations.
+    const FermiLevel level = find_fermi_level(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, is, isk);
+    ef = level.reference + level.offset;
     demet = 0.0;
 
     for (int ik = 0; ik < nks; ik++)
@@ -244,7 +277,9 @@ void Occupy::gweights(const int nks,
             // Calculate the gaussian weights
             //================================
             // call wgauss
-            wg(ik, ib) = wk[ik] * occupy_smearing::wgauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);
+            const double relative_energy = ekb(ik, ib) - level.reference;
+            const double argument = (level.offset - relative_energy) / smearing_sigma;
+            wg(ik, ib) = wk[ik] * occupy_smearing::wgauss(argument, ngauss);
 
             //====================================================================
             // The correct form of the band energy is  \int e n(e) de   for e<ef
@@ -252,7 +287,7 @@ void Occupy::gweights(const int nks,
             //====================================================================
             // Mohan fix bug 2010-1-9
             // call w1gauss
-            demet += wk[ik] * smearing_sigma * occupy_smearing::w1gauss((ef - ekb(ik, ib)) / smearing_sigma, ngauss);
+            demet += wk[ik] * smearing_sigma * occupy_smearing::w1gauss(argument, ngauss);
         }
     }
 
@@ -260,7 +295,7 @@ void Occupy::gweights(const int nks,
 } // end subroutine gweights
 
 /**
- * @brief Finds the Fermi energy
+ * @brief Finds the Fermi energy in a centered coordinate
  *
  * @param ekb band energy.
  * @param nband number of bands.
@@ -269,20 +304,21 @@ void Occupy::gweights(const int nks,
  * @param wk wk weight of each k point(symmetry considered).
  * @param smearing_sigma parameter input by user.
  * @param ngauss which type of smearing.
- * @param ef output: fermi level
+ * @return energy reference and converged relative Fermi level
  * @param is spin
  * @param isk array to point out each k belong to which spin
  */
-void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
-                     const int nband,
-                     const int nks,
-                     const double& nelec,
-                     const std::vector<double>& wk,
-                     const double& smearing_sigma,
-                     const int ngauss,
-                     double& ef,
-                     const int& is,
-                     const std::vector<int>& isk)
+namespace
+{
+FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
+                           const int nband,
+                           const int nks,
+                           const double& nelec,
+                           const std::vector<double>& wk,
+                           const double& smearing_sigma,
+                           const int ngauss,
+                           const int& is,
+                           const std::vector<int>& isk)
 {
     // Resolve particle number to floating-point accuracy, independently of the
     // outer SCF threshold. A fixed absolute tolerance leaves a charge floor.
@@ -319,11 +355,17 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
         ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "No finite eigenvalue bounds");
     }
 
+    // Keep the energy origin separate from the small chemical-potential offset.
+    // Rounding reference + offset before evaluating occupations can destroy
+    // particle-number accuracy even for a harmless constant shift of all levels.
+    const double reference = 0.5 * lower + 0.5 * upper;
+    lower -= reference;
+    upper -= reference;
     double step = 2.0 * smearing_sigma;
     lower -= step;
     upper += step;
-    double low_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, lower, is, isk);
-    double high_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, upper, is, isk);
+    double low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
+    double high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
     int expansion = 0;
     while (low_count > nelec + tolerance || high_count < nelec - tolerance)
     {
@@ -335,12 +377,12 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
         if (low_count > nelec + tolerance)
         {
             lower -= step;
-            low_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, lower, is, isk);
+            low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
         }
         if (high_count < nelec - tolerance)
         {
             upper += step;
-            high_count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, upper, is, isk);
+            high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
         }
     }
 
@@ -352,19 +394,19 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
 
     // Keep the best evaluated point, including the saturated endpoints. For
     // empty/full bands or a resolved gap, the chemical potential is not unique.
-    ef = lower;
+    FermiLevel result{reference, lower};
     double best_error = std::abs(low_count - nelec);
     const double high_error = std::abs(high_count - nelec);
     if (high_error < best_error)
     {
-        ef = upper;
+        result.offset = upper;
         best_error = high_error;
     }
     for (int iteration = 0; iteration < 300; ++iteration)
     {
         if (best_error <= tolerance)
         {
-            return;
+            return result;
         }
         const double middle = 0.5 * lower + 0.5 * upper;
         if (middle == lower || middle == upper)
@@ -377,9 +419,9 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
                     << best_error << ", requested tolerance = " << tolerance;
             const std::string diagnostic = message.str();
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", diagnostic);
-            return;
+            return result;
         }
-        const double count = sumkg(ekb, nband, nks, wk, smearing_sigma, ngauss, middle, is, isk);
+        const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, middle, is, isk);
         if (!std::isfinite(count))
         {
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Non-finite electron number");
@@ -387,12 +429,12 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
         const double error = std::abs(count - nelec);
         if (error < best_error)
         {
-            ef = middle;
+            result.offset = middle;
             best_error = error;
         }
         if (best_error <= tolerance)
         {
-            return;
+            return result;
         }
         if (count < nelec)
         {
@@ -404,6 +446,25 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
         }
     }
     ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Fermi search exceeded its iteration limit");
+    return result;
+}
+
+} // namespace
+
+void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
+                             const int nband,
+                             const int nks,
+                             const double& nelec,
+                             const std::vector<double>& wk,
+                             const double& smearing_sigma,
+                             const int ngauss,
+                             double& ef,
+                             const int& is,
+                             const std::vector<int>& isk)
+{
+    const FermiLevel level = find_fermi_level(ekb, nband, nks, nelec, wk, smearing_sigma, ngauss, is, isk);
+    // Public output is rounded to an absolute double; gweights retains both parts.
+    ef = level.reference + level.offset;
 }
 
 /**
@@ -415,20 +476,24 @@ void occupy_smearing::efermig(const ModuleBase::matrix& ekb,
  * @param wk weight of each k point(symmetry considered).
  * @param smearing_sigma parameter input by user
  * @param ngauss which type of smearing.
- * @param e a givern energy
+ * @param reference common energy origin
+ * @param offset Fermi level relative to the common origin
  * @param is spin
  * @param isk array to point out each k belong to which spin
  * @return (double) the number of states
  */
-double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
-                     const int nband,
-                     const int nks,
-                     const std::vector<double>& wk,
-                     const double& smearing_sigma,
-                     const int ngauss,
-                     const double& e,
-                     const int& is,
-                     const std::vector<int>& isk)
+namespace
+{
+double sum_occupations(const ModuleBase::matrix& ekb,
+                       const int nband,
+                       const int nks,
+                       const std::vector<double>& wk,
+                       const double& smearing_sigma,
+                       const int ngauss,
+                       const double reference,
+                       const double offset,
+                       const int& is,
+                       const std::vector<int>& isk)
 {
     // ModuleBase::TITLE("Occupy","sumkg");
     double sum2 = 0.0;
@@ -443,8 +508,9 @@ double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
         {
             // Sum the same weighted occupations as gweights, compensating the
             // local accumulation before the existing pool/band reduction.
-            const double argument = (e - ekb(ik, ib)) / smearing_sigma;
-            const double occupation = wk[ik] * wgauss(argument, ngauss);
+            const double relative_energy = ekb(ik, ib) - reference;
+            const double argument = (offset - relative_energy) / smearing_sigma;
+            const double occupation = wk[ik] * occupy_smearing::wgauss(argument, ngauss);
             const double corrected = occupation - compensation;
             const double next = sum2 + corrected;
             compensation = (next - sum2) - corrected;
@@ -463,6 +529,21 @@ double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
 
     return sum2;
 }
+} // namespace
+
+double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
+                             const int nband,
+                             const int nks,
+                             const std::vector<double>& wk,
+                             const double& smearing_sigma,
+                             const int ngauss,
+                             const double& e,
+                             const int& is,
+                             const std::vector<int>& isk)
+{
+    const double reference = 0.0;
+    return sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, e, is, isk);
+}
 
 double occupy_smearing::wgauss(const double& x, const int n)
 {
diff --git a/source/source_estate/occupy.h b/source/source_estate/occupy.h
index abb30b8d26c..087686f1493 100644
--- a/source/source_estate/occupy.h
+++ b/source/source_estate/occupy.h
@@ -36,7 +36,10 @@ double sumkg(const ModuleBase::matrix& ekb,
              const int& is,
              const std::vector<int>& isk);
 
-/// @brief find the Fermi level by bisection on sumkg()
+/// @brief Find the Fermi level using a compensated, energy-centered electron count.
+/// @note ef is the rounded absolute chemical potential. Re-evaluating sumkg(ef)
+///       can lose precision for a large energy offset; gweights retains the
+///       centered chemical potential internally when forming occupations.
 void efermig(const ModuleBase::matrix& ekb,
              const int nbnd,
              const int nks,
diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/elecstate_occupy_test.cpp
index 34d9efc8571..f6e70589c9a 100644
--- a/source/source_estate/test/elecstate_occupy_test.cpp
+++ b/source/source_estate/test/elecstate_occupy_test.cpp
@@ -433,15 +433,16 @@ TEST_F(OccupyTest, CompensatedWeightedOccupationSum)
 
 TEST_F(OccupyTest, FermiSearchAtRepresentableEnergyLimit)
 {
-    ModuleBase::matrix energies(1, 1);
-    energies(0, 0) = 1e8;
+    ModuleBase::matrix energies(1, 2);
+    energies(0, 0) = 0.0;
+    energies(0, 1) = 1e8;
     const std::vector<double> weights{1.0};
     const std::vector<int> spins{0};
     double ef = 0.0;
-    // The width is smaller than one ULP of the eigenvalue. No double-valued
-    // chemical potential can yield the requested fractional occupation.
+    // Unlike a common energy shift, this enormous spectral span still leaves
+    // the root unresolved in the centered coordinate. It must not pass silently.
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(occupy_smearing::efermig(energies, 1, 1, 0.3, weights, 1e-8, 0, ef, -1, spins),
+    EXPECT_EXIT(occupy_smearing::efermig(energies, 2, 1, 0.3, weights, 1e-8, 0, ef, -1, spins),
                 ::testing::ExitedWithCode(1), "");
     const std::string diagnostic = testing::internal::GetCapturedStdout();
     EXPECT_THAT(diagnostic, testing::HasSubstr("Fermi bracket reached floating-point resolution"));
@@ -528,3 +529,107 @@ TEST_F(OccupyTest, GweightsHonorsExplicitBandRange)
     EXPECT_DOUBLE_EQ(occupations(0, 1), -123.0);
     EXPECT_DOUBLE_EQ(occupations(0, 2), -123.0);
 }
+
+TEST_F(OccupyTest, DegenerateSpectrumEnergyShift)
+{
+    // Reproduce ElecStateTest.CalculateWeightsGWeights without changing its fixture.
+    const int points = 5;
+    const int bands = 6;
+    ModuleBase::matrix energies(points, bands);
+    ModuleBase::matrix occupations(points, bands);
+    const std::vector<double> weights(points, 2.0);
+    const std::vector<int> spins(points, 0);
+    const double target = 10.0;
+    const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * target;
+    for (const double shift : {0.0, 100.0, 200.0, 1e8})
+    {
+        for (int ik = 0; ik < points; ++ik)
+        {
+            for (int ib = 0; ib < bands; ++ib)
+            {
+                energies(ik, ib) = shift;
+            }
+        }
+        double ef = 0.0;
+        double entropy = 0.0;
+        occupy.gweights(points, weights, bands, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+        long double count = 0.0L;
+        for (int ik = 0; ik < points; ++ik)
+        {
+            for (int ib = 0; ib < bands; ++ib)
+            {
+                EXPECT_NEAR(occupations(ik, ib), 1.0 / 3.0, 2e-15);
+                count += occupations(ik, ib);
+            }
+        }
+        EXPECT_NEAR(count, target, tolerance) << shift;
+        if (shift == 100.0)
+        {
+            const double residual = std::abs(static_cast<double>(count - target));
+            std::ostringstream value;
+            value.precision(17);
+            value << residual;
+            RecordProperty("shift100_particle_residual", value.str());
+        }
+    }
+}
+
+TEST_F(OccupyTest, OccupationAndEntropyEnergyShift)
+{
+    ModuleBase::matrix energies(1, 3);
+    ModuleBase::matrix occupations(1, 3);
+    ModuleBase::matrix reference(1, 3);
+    const std::vector<double> weights{2.0};
+    const std::vector<int> spins{0};
+    const double target = 2.3;
+    for (const int smearing : {0, 1, 2, -1, -99})
+    {
+        double reference_entropy = 0.0;
+        double reference_fermi = 0.0;
+        for (const double shift : {0.0, 100.0, 200.0, 1e8})
+        {
+            for (int ib = 0; ib < 3; ++ib)
+            {
+                // Binary fractions ensure the shifted input spectrum is exact.
+                energies(0, ib) = shift - 0.5 + 0.5 * ib;
+            }
+            double ef = 0.0;
+            double entropy = 0.0;
+            occupy.gweights(1, weights, 3, target, 0.125, smearing, energies, ef, entropy, occupations, -1, spins);
+            if (shift == 0.0)
+            {
+                reference_entropy = entropy;
+                reference_fermi = ef;
+                for (int ib = 0; ib < 3; ++ib)
+                {
+                    reference(0, ib) = occupations(0, ib);
+                }
+            }
+            double count = 0.0;
+            for (int ib = 0; ib < 3; ++ib)
+            {
+                EXPECT_NEAR(occupations(0, ib), reference(0, ib), 2e-15);
+                count += occupations(0, ib);
+            }
+            EXPECT_NEAR(count, target, 8e-15);
+            EXPECT_NEAR(entropy, reference_entropy, 2e-15);
+            const double expected_fermi = reference_fermi + shift;
+            EXPECT_DOUBLE_EQ(ef, expected_fermi);
+        }
+    }
+}
+
+TEST_F(OccupyTest, NarrowSmearingAtLargeEnergyOffset)
+{
+    ModuleBase::matrix energies(1, 1);
+    energies(0, 0) = 1e8;
+    ModuleBase::matrix occupations(1, 1);
+    const std::vector<double> weights{1.0};
+    const std::vector<int> spins{0};
+    double ef = 0.0;
+    double entropy = 0.0;
+    occupy.gweights(1, weights, 1, 0.3, 1e-8, 0, energies, ef, entropy, occupations, -1, spins);
+    EXPECT_NEAR(occupations(0, 0), 0.3, 2e-15);
+    EXPECT_TRUE(std::isfinite(ef));
+    EXPECT_TRUE(std::isfinite(entropy));
+}

From 7eb5b12d33090ee64bd747d99a6a78afa341c16f Mon Sep 17 00:00:00 2001
From: Chen Chengbing <1747193328@qq.com>
Date: Mon, 12 Oct 2026 13:25:49 +0800
Subject: [PATCH 3/5] Keep the Fermi energy origin near the shrinking bracket

---
 source/source_estate/occupy.cpp               | 15 +++++++-
 .../test/elecstate_occupy_test.cpp            | 38 ++++++++++++++++---
 2 files changed, 47 insertions(+), 6 deletions(-)

diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp
index f00af5664d5..fcbd726e25b 100644
--- a/source/source_estate/occupy.cpp
+++ b/source/source_estate/occupy.cpp
@@ -358,7 +358,7 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
     // Keep the energy origin separate from the small chemical-potential offset.
     // Rounding reference + offset before evaluating occupations can destroy
     // particle-number accuracy even for a harmless constant shift of all levels.
-    const double reference = 0.5 * lower + 0.5 * upper;
+    double reference = 0.5 * lower + 0.5 * upper;
     lower -= reference;
     upper -= reference;
     double step = 2.0 * smearing_sigma;
@@ -408,6 +408,18 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         {
             return result;
         }
+        // Follow the shrinking bracket: a fixed spectral midpoint can be far
+        // from the Fermi surface when there are high unoccupied bands. Keep
+        // the small interval in offsets even when the absolute origin rounds.
+        const double trial_offset = 0.5 * lower + 0.5 * upper;
+        const double next_reference = reference + trial_offset;
+        const double shift = next_reference - reference;
+        if (shift != 0.0 && std::isfinite(next_reference))
+        {
+            lower -= shift;
+            upper -= shift;
+            reference = next_reference;
+        }
         const double middle = 0.5 * lower + 0.5 * upper;
         if (middle == lower || middle == upper)
         {
@@ -429,6 +441,7 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         const double error = std::abs(count - nelec);
         if (error < best_error)
         {
+            result.reference = reference;
             result.offset = middle;
             best_error = error;
         }
diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/elecstate_occupy_test.cpp
index f6e70589c9a..984fb221ee6 100644
--- a/source/source_estate/test/elecstate_occupy_test.cpp
+++ b/source/source_estate/test/elecstate_occupy_test.cpp
@@ -433,16 +433,16 @@ TEST_F(OccupyTest, CompensatedWeightedOccupationSum)
 
 TEST_F(OccupyTest, FermiSearchAtRepresentableEnergyLimit)
 {
-    ModuleBase::matrix energies(1, 2);
+    ModuleBase::matrix energies(1, 1);
     energies(0, 0) = 0.0;
-    energies(0, 1) = 1e8;
     const std::vector<double> weights{1.0};
     const std::vector<int> spins{0};
     double ef = 0.0;
-    // Unlike a common energy shift, this enormous spectral span still leaves
-    // the root unresolved in the centered coordinate. It must not pass silently.
+    // A subnormal width leaves too few representable offsets to resolve this
+    // occupation, even after moving the energy origin close to the root.
+    const double sigma = 1024.0 * std::numeric_limits<double>::denorm_min();
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(occupy_smearing::efermig(energies, 2, 1, 0.3, weights, 1e-8, 0, ef, -1, spins),
+    EXPECT_EXIT(occupy_smearing::efermig(energies, 1, 1, 0.3, weights, sigma, 0, ef, -1, spins),
                 ::testing::ExitedWithCode(1), "");
     const std::string diagnostic = testing::internal::GetCapturedStdout();
     EXPECT_THAT(diagnostic, testing::HasSubstr("Fermi bracket reached floating-point resolution"));
@@ -633,3 +633,31 @@ TEST_F(OccupyTest, NarrowSmearingAtLargeEnergyOffset)
     EXPECT_TRUE(std::isfinite(ef));
     EXPECT_TRUE(std::isfinite(entropy));
 }
+
+TEST_F(OccupyTest, FermiSurfaceAwayFromSpectralMidpoint)
+{
+    ModuleBase::matrix energies(1, 7);
+    ModuleBase::matrix occupations(1, 7);
+    const std::vector<double> weights{2.0};
+    const std::vector<int> spins{0};
+    const double target = 10.0;
+    const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * target;
+    for (const double shift : {0.0, 100.0, 1e8})
+    {
+        for (int ib = 0; ib < 6; ++ib)
+        {
+            energies(0, ib) = shift;
+        }
+        energies(0, 6) = shift + 100.0;
+        double ef = 0.0;
+        double entropy = 0.0;
+        occupy.gweights(1, weights, 7, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+        long double count = 0.0L;
+        for (int ib = 0; ib < 7; ++ib)
+        {
+            count += occupations(0, ib);
+        }
+        EXPECT_NEAR(count, target, tolerance) << shift;
+        EXPECT_DOUBLE_EQ(occupations(0, 6), 0.0);
+    }
+}

From 70866e2f6cd498e1c83c27344b9c5d006ba79f45 Mon Sep 17 00:00:00 2001
From: Chen Chengbing <1747193328@qq.com>
Date: Tue, 13 Oct 2026 16:54:06 +0800
Subject: [PATCH 4/5] Bound Fermi stagnation by an explicit roundoff budget

---
 source/source_estate/occupy.cpp               | 168 ++++++++++++++++--
 source/source_estate/occupy.h                 |   5 +
 .../test/elecstate_occupy_test.cpp            |  73 ++++++++
 3 files changed, 229 insertions(+), 17 deletions(-)

diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp
index fcbd726e25b..c783a4bd441 100644
--- a/source/source_estate/occupy.cpp
+++ b/source/source_estate/occupy.cpp
@@ -37,7 +37,85 @@ double sum_occupations(const ModuleBase::matrix& ekb,
                        const double reference,
                        const double offset,
                        const int& is,
-                       const std::vector<int>& isk);
+                       const std::vector<int>& isk,
+                       double* roundoff);
+
+// First-order scalar error estimate, used only when the Fermi search stalls.
+// The 8u/4u allowances for erfc/exp are tested libm assumptions, not a portable
+// correctly-rounded guarantee. Keep MP component magnitudes before cancellation.
+double smearing_roundoff(const double x, const int n, const double value, const double x_error)
+{
+    const double u = 0.5 * std::numeric_limits<double>::epsilon();
+    const double tiny = std::numeric_limits<double>::denorm_min();
+    if (n == -99)
+    {
+        const double magnitude = std::abs(x);
+        if (magnitude > 200.0)
+        {
+            const double tail_argument = -magnitude;
+            return std::exp(tail_argument) + tiny;
+        }
+        const double slope = value * (1.0 - value);
+        return u * (2.0 * std::abs(value) + 4.0 * slope) + slope * x_error + tiny;
+    }
+    const double negative_x = -x;
+    const double gaussian = 0.5 * std::erfc(negative_x);
+    const double x2 = x * x;
+    const double negative_x2 = -x2;
+    const double density = std::exp(negative_x2) / std::sqrt(ModuleBase::PI);
+    if (n == 0)
+    {
+        return 8.0 * u * gaussian + density * x_error + tiny;
+    }
+    if (n == -1)
+    {
+        const double shift = 1.0 / ModuleBase::SQRT2;
+        const double xp = x - shift;
+        const double xp2 = xp * xp;
+        const double negative_xp = -xp;
+        const double base = 0.5 * std::erfc(negative_xp);
+        const double argument = std::min(200.0, xp2);
+        const double negative_argument = -argument;
+        const double correction = std::exp(negative_argument) / std::sqrt(ModuleBase::TWO_PI);
+        if (xp2 > 200.0)
+        {
+            // The existing exponential clamp leaves a nonzero tail.
+            return 8.0 * u * base + correction + tiny;
+        }
+        const double derivative = (ModuleBase::SQRT2 - 2.0 * xp) * correction;
+        const double slope = std::abs(derivative);
+        const double shifted_error = x_error + u * (std::abs(xp) + 2.0 * shift);
+        // Square, exp, normalization and addition; include the shift's rounding.
+        return u * (8.0 * base + (xp2 + 10.0) * correction + std::abs(value)) + slope * shifted_error + tiny;
+    }
+    if (n != 1 && n != 2)
+    {
+        return std::numeric_limits<double>::infinity();
+    }
+    const double magnitude = std::abs(x);
+    double polynomial = 0.5 * magnitude;
+    double derivative = 1.5 - x2;
+    double slope = std::abs(derivative) * density;
+    double operations = 8.0;
+    if (n == 2)
+    {
+        polynomial = magnitude * (0.875 + 0.25 * x2);
+        derivative = 1.875 - 2.5 * x2 + 0.5 * x2 * x2;
+        slope = std::abs(derivative) * density;
+        operations = 12.0;
+    }
+    const double argument = std::min(200.0, x2);
+    const double negative_argument = -argument;
+    const double correction_scale = polynomial * std::exp(negative_argument) / std::sqrt(ModuleBase::PI);
+    if (x2 > 200.0)
+    {
+        return 8.0 * u * gaussian + correction_scale + tiny;
+    }
+    // 8/12 count first-order rounding allowances in the MP1/MP2 recurrence,
+    // including the exponential and normalization, not an electron tolerance.
+    return u * (8.0 * gaussian + (operations + x2) * correction_scale + n * std::abs(value))
+           + slope * x_error + tiny;
+}
 } // namespace
 
 Occupy::Occupy()
@@ -364,8 +442,8 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
     double step = 2.0 * smearing_sigma;
     lower -= step;
     upper += step;
-    double low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
-    double high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
+    double low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk, nullptr);
+    double high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk, nullptr);
     int expansion = 0;
     while (low_count > nelec + tolerance || high_count < nelec - tolerance)
     {
@@ -377,12 +455,12 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         if (low_count > nelec + tolerance)
         {
             lower -= step;
-            low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
+            low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk, nullptr);
         }
         if (high_count < nelec - tolerance)
         {
             upper += step;
-            high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
+            high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk, nullptr);
         }
     }
 
@@ -402,7 +480,10 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         result.offset = upper;
         best_error = high_error;
     }
-    for (int iteration = 0; iteration < 300; ++iteration)
+    // Inherited safety cap, not an accuracy criterion. Stagnation is checked
+    // separately; exhausting this cap must never report convergence.
+    const int maxiter = 300;
+    for (int iteration = 0; iteration < maxiter; ++iteration)
     {
         if (best_error <= tolerance)
         {
@@ -423,17 +504,30 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         const double middle = 0.5 * lower + 0.5 * upper;
         if (middle == lower || middle == upper)
         {
-            // No representable interior point remains. Do not silently spin to
-            // maxiter, nor claim that a small energy interval proves charge accuracy.
+            // Audit only the best candidate, in its original centered coordinates.
+            // A large error estimate must not buy a loose electron tolerance.
+            double roundoff = 0.0;
+            const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss,
+                                                 result.reference, result.offset, is, isk, &roundoff);
+            const double count_error = count - nelec;
+            const double residual = std::abs(count_error);
             std::ostringstream message;
             message.precision(17);
+            if (std::isfinite(roundoff) && roundoff <= tolerance && residual <= tolerance + roundoff)
+            {
+                message << "Fermi search reached rounding accuracy; electron residual = " << residual
+                        << ", target = " << tolerance << ", estimated roundoff = " << roundoff;
+                const std::string diagnostic = message.str();
+                ModuleBase::WARNING("occupy_smearing::efermig", diagnostic);
+                return result;
+            }
             message << "Fermi bracket reached floating-point resolution; electron residual = "
-                    << best_error << ", requested tolerance = " << tolerance;
+                    << residual << ", requested tolerance = " << tolerance << ", estimated roundoff = " << roundoff;
             const std::string diagnostic = message.str();
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", diagnostic);
             return result;
         }
-        const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, middle, is, isk);
+        const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, middle, is, isk, nullptr);
         if (!std::isfinite(count))
         {
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Non-finite electron number");
@@ -506,11 +600,15 @@ double sum_occupations(const ModuleBase::matrix& ekb,
                        const double reference,
                        const double offset,
                        const int& is,
-                       const std::vector<int>& isk)
+                       const std::vector<int>& isk,
+                       double* roundoff)
 {
     // ModuleBase::TITLE("Occupy","sumkg");
     double sum2 = 0.0;
     double compensation = 0.0;
+    const double u = 0.5 * std::numeric_limits<double>::epsilon();
+    double absolute_sum = 0.0;
+    double evaluation_error = 0.0;
     for (int ik = 0; ik < nks; ++ik)
     {
         if (is != -1 && is != isk[ik])
@@ -523,24 +621,58 @@ double sum_occupations(const ModuleBase::matrix& ekb,
             // local accumulation before the existing pool/band reduction.
             const double relative_energy = ekb(ik, ib) - reference;
             const double argument = (offset - relative_energy) / smearing_sigma;
-            const double occupation = wk[ik] * occupy_smearing::wgauss(argument, ngauss);
+            const double value = occupy_smearing::wgauss(argument, ngauss);
+            const double occupation = wk[ik] * value;
             const double corrected = occupation - compensation;
             const double next = sum2 + corrected;
             compensation = (next - sum2) - corrected;
             sum2 = next;
+            if (roundoff != nullptr)
+            {
+                const double magnitude = std::abs(occupation);
+                absolute_sum += magnitude;
+                const double difference = offset - relative_energy;
+                const double coordinate_scale = std::abs(relative_energy) / smearing_sigma
+                                                + std::abs(difference) / smearing_sigma + std::abs(argument);
+                const double x_error = u * coordinate_scale;
+                const double scalar_error = smearing_roundoff(argument, ngauss, value, x_error);
+                evaluation_error += std::abs(wk[ik]) * scalar_error + u * magnitude;
+            }
         }
     }
 
     // GlobalV::ofs_running << "\n sum2 before reduce = " << sum2 << std::endl;
 
+    double totals[3] = {sum2, absolute_sum, evaluation_error};
+    double reduction_error = 0.0;
 #ifdef __MPI
     const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, sum2);
+    const int count = roundoff == nullptr ? 1 : 3;
+    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, totals, count);
+    if (roundoff != nullptr && npool > 1)
+    {
+        int participants = 1;
+        Parallel_Reduce::reduce_all(participants);
+        const double additions = (participants - 1) * u;
+        // Pool scaling plus an order-independent reduction bound; used only
+        // for the capped stagnation check, never for early termination.
+        reduction_error = (u + additions / (1.0 - additions)) * totals[1];
+    }
 #endif
+    if (roundoff != nullptr)
+    {
+        // First-order Kahan bound, function/argument evaluation, MPI and the
+        // final count subtraction. Subnormal widths invalidate this model.
+        *roundoff = totals[2] + 2.0 * u * totals[1] + reduction_error + u * std::abs(totals[0]);
+        if (smearing_sigma < std::numeric_limits<double>::min())
+        {
+            *roundoff = std::numeric_limits<double>::infinity();
+        }
+    }
 
     // GlobalV::ofs_running << "\n sum2 after reduce = " << sum2 << std::endl;
 
-    return sum2;
+    return totals[0];
 }
 } // namespace
 
@@ -555,7 +687,7 @@ double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
                              const std::vector<int>& isk)
 {
     const double reference = 0.0;
-    return sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, e, is, isk);
+    return sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, e, is, isk, nullptr);
 }
 
 double occupy_smearing::wgauss(const double& x, const int n)
@@ -603,14 +735,16 @@ double occupy_smearing::wgauss(const double& x, const int n)
     {
         const double xp = x - 1.00 / ModuleBase::SQRT2;
         const double arg = std::min(maxarg, xp * xp);
-        wga = 0.50 * erf(xp) + 1.00 / sqrt(ModuleBase::TWO_PI) * std::exp(-arg) + 0.50;
+        const double negative_xp = -xp;
+        wga = 0.50 * std::erfc(negative_xp) + 1.00 / sqrt(ModuleBase::TWO_PI) * std::exp(-arg);
         return wga;
     }
 
     //====================
     // Methfessel-Paxton       //pengfei 2014-10-13
     //====================
-    wga = 0.5 * (1 - erf(-x));
+    const double negative_x = -x;
+    wga = 0.5 * std::erfc(negative_x);
     // wga = gauss_freq(x * ModuleBase::SQRT2);
     //	std::cout<<"\n x="<<x<<" wga="<<wga;
     if (n == 0)
diff --git a/source/source_estate/occupy.h b/source/source_estate/occupy.h
index 087686f1493..d63abe53632 100644
--- a/source/source_estate/occupy.h
+++ b/source/source_estate/occupy.h
@@ -40,6 +40,11 @@ double sumkg(const ModuleBase::matrix& ekb,
 /// @note ef is the rounded absolute chemical potential. Re-evaluating sumkg(ef)
 ///       can lose precision for a large energy offset; gweights retains the
 ///       centered chemical potential internally when forming occupations.
+///       The usual electron-count target is 8*epsilon*max(1,abs(nelec)). At
+///       floating-point stagnation only, a diagnosed rounding-limited result
+///       may exceed it by an estimated error no larger than the target itself.
+///       The estimate is a first-order libm/arithmetic model, not a certified
+///       physical-accuracy bound; a larger estimate still causes failure.
 void efermig(const ModuleBase::matrix& ekb,
              const int nbnd,
              const int nks,
diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/elecstate_occupy_test.cpp
index 984fb221ee6..dea764a43b4 100644
--- a/source/source_estate/test/elecstate_occupy_test.cpp
+++ b/source/source_estate/test/elecstate_occupy_test.cpp
@@ -447,6 +447,7 @@ TEST_F(OccupyTest, FermiSearchAtRepresentableEnergyLimit)
     const std::string diagnostic = testing::internal::GetCapturedStdout();
     EXPECT_THAT(diagnostic, testing::HasSubstr("Fermi bracket reached floating-point resolution"));
     EXPECT_THAT(diagnostic, testing::HasSubstr("electron residual"));
+    EXPECT_THAT(diagnostic, testing::HasSubstr("estimated roundoff = inf"));
 }
 
 TEST_F(OccupyTest, RejectInvalidWidthAndUnbracketedCharge)
@@ -661,3 +662,75 @@ TEST_F(OccupyTest, FermiSurfaceAwayFromSpectralMidpoint)
         EXPECT_DOUBLE_EQ(occupations(0, 6), 0.0);
     }
 }
+
+TEST_F(OccupyTest, ReviewerDegenerateGaussianSpectrum)
+{
+    // PR #8130: the former 1-erf evaluation stalls at 3.99680288865e-15.
+    const int bands = 74;
+    ModuleBase::matrix energies(1, bands);
+    ModuleBase::matrix occupations(1, bands);
+    const std::vector<double> weights{2.0};
+    const std::vector<int> spins{0};
+    const double target = 2.0;
+    const double tolerance = 8.0 * std::numeric_limits<double>::epsilon() * target;
+    for (const double shift : {0.0, 100.0, 1e8})
+    {
+        for (int ib = 0; ib < bands; ++ib)
+        {
+            energies(0, ib) = shift;
+        }
+        double ef = 0.0;
+        double entropy = 0.0;
+        occupy.gweights(1, weights, bands, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+        long double count = 0.0L;
+        for (int ib = 0; ib < bands; ++ib)
+        {
+            count += occupations(0, ib);
+        }
+        EXPECT_NEAR(count, target, tolerance) << shift;
+        EXPECT_TRUE(std::isfinite(entropy));
+    }
+}
+
+TEST_F(OccupyTest, SmearingHighPrecisionReference)
+{
+    // Independent 100-decimal-digit erfc/exp evaluations of the closed forms.
+    // Include signed MP tails and cold-smearing overshoot; do not clamp to [0,1].
+    struct Reference
+    {
+        int order;
+        double x;
+        double occupation;
+    };
+    const Reference references[] = {
+        {0, -6.0, 1.0759868356249456558e-17},
+        {0, -8.0, 5.61214858649146354e-30},
+        {-1, -6.0, 1.2797773847197835025e-20},
+        {-1, -8.0, 5.1183964952735562742e-34},
+        {0, -2.0, 0.002338867490523632918965},
+        {0, -0.5, 0.2397500610934767311586},
+        {0, 0.5, 0.7602499389065232688414},
+        {0, 2.0, 0.997661132509476367081},
+        {1, -2.0, -0.007994625186522394009569},
+        {1, -0.5, 0.1299022387265461318969},
+        {1, 0.5, 0.8700977612734538681031},
+        {1, 2.0, 1.00799462518652239401},
+        {2, -2.0, 0.004922240659785139651099},
+        {2, -0.5, 0.06124734974721450735834},
+        {2, 0.5, 0.9387526502527854926417},
+        {2, 2.0, 0.9950777593402148603489},
+        {-1, -2.0, 0.0003264301114214816118753},
+        {-1, -0.5, 0.1368184140991744579818},
+        {-1, 0.5, 0.766994049270187616701},
+        {-1, 2.0, 1.041238947228640924942},
+        {-99, -2.0, 0.1192029220221175559403},
+        {-99, -0.5, 0.3775406687981454353611},
+        {-99, 0.5, 0.6224593312018545646389},
+        {-99, 2.0, 0.8807970779778824440597}};
+    for (const Reference& reference : references)
+    {
+        const double value = occupy_smearing::wgauss(reference.x, reference.order);
+        const double tolerance = 32.0 * std::numeric_limits<double>::epsilon() * std::abs(reference.occupation);
+        EXPECT_NEAR(value, reference.occupation, tolerance) << reference.order << ", x = " << reference.x;
+    }
+}

From d31dec557698499781ffee982bff0d37d7012520 Mon Sep 17 00:00:00 2001
From: Chen Chengbing <1747193328@qq.com>
Date: Tue, 13 Oct 2026 22:21:34 +0800
Subject: [PATCH 5/5] Simplify Fermi convergence tolerance and focused
 regressions

---
 source/source_estate/occupy.cpp               | 165 ++-------------
 source/source_estate/occupy.h                 |   7 +-
 .../test/elecstate_occupy_test.cpp            | 188 ++++--------------
 3 files changed, 63 insertions(+), 297 deletions(-)

diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp
index c783a4bd441..c573052fe1a 100644
--- a/source/source_estate/occupy.cpp
+++ b/source/source_estate/occupy.cpp
@@ -37,85 +37,7 @@ double sum_occupations(const ModuleBase::matrix& ekb,
                        const double reference,
                        const double offset,
                        const int& is,
-                       const std::vector<int>& isk,
-                       double* roundoff);
-
-// First-order scalar error estimate, used only when the Fermi search stalls.
-// The 8u/4u allowances for erfc/exp are tested libm assumptions, not a portable
-// correctly-rounded guarantee. Keep MP component magnitudes before cancellation.
-double smearing_roundoff(const double x, const int n, const double value, const double x_error)
-{
-    const double u = 0.5 * std::numeric_limits<double>::epsilon();
-    const double tiny = std::numeric_limits<double>::denorm_min();
-    if (n == -99)
-    {
-        const double magnitude = std::abs(x);
-        if (magnitude > 200.0)
-        {
-            const double tail_argument = -magnitude;
-            return std::exp(tail_argument) + tiny;
-        }
-        const double slope = value * (1.0 - value);
-        return u * (2.0 * std::abs(value) + 4.0 * slope) + slope * x_error + tiny;
-    }
-    const double negative_x = -x;
-    const double gaussian = 0.5 * std::erfc(negative_x);
-    const double x2 = x * x;
-    const double negative_x2 = -x2;
-    const double density = std::exp(negative_x2) / std::sqrt(ModuleBase::PI);
-    if (n == 0)
-    {
-        return 8.0 * u * gaussian + density * x_error + tiny;
-    }
-    if (n == -1)
-    {
-        const double shift = 1.0 / ModuleBase::SQRT2;
-        const double xp = x - shift;
-        const double xp2 = xp * xp;
-        const double negative_xp = -xp;
-        const double base = 0.5 * std::erfc(negative_xp);
-        const double argument = std::min(200.0, xp2);
-        const double negative_argument = -argument;
-        const double correction = std::exp(negative_argument) / std::sqrt(ModuleBase::TWO_PI);
-        if (xp2 > 200.0)
-        {
-            // The existing exponential clamp leaves a nonzero tail.
-            return 8.0 * u * base + correction + tiny;
-        }
-        const double derivative = (ModuleBase::SQRT2 - 2.0 * xp) * correction;
-        const double slope = std::abs(derivative);
-        const double shifted_error = x_error + u * (std::abs(xp) + 2.0 * shift);
-        // Square, exp, normalization and addition; include the shift's rounding.
-        return u * (8.0 * base + (xp2 + 10.0) * correction + std::abs(value)) + slope * shifted_error + tiny;
-    }
-    if (n != 1 && n != 2)
-    {
-        return std::numeric_limits<double>::infinity();
-    }
-    const double magnitude = std::abs(x);
-    double polynomial = 0.5 * magnitude;
-    double derivative = 1.5 - x2;
-    double slope = std::abs(derivative) * density;
-    double operations = 8.0;
-    if (n == 2)
-    {
-        polynomial = magnitude * (0.875 + 0.25 * x2);
-        derivative = 1.875 - 2.5 * x2 + 0.5 * x2 * x2;
-        slope = std::abs(derivative) * density;
-        operations = 12.0;
-    }
-    const double argument = std::min(200.0, x2);
-    const double negative_argument = -argument;
-    const double correction_scale = polynomial * std::exp(negative_argument) / std::sqrt(ModuleBase::PI);
-    if (x2 > 200.0)
-    {
-        return 8.0 * u * gaussian + correction_scale + tiny;
-    }
-    // 8/12 count first-order rounding allowances in the MP1/MP2 recurrence,
-    // including the exponential and normalization, not an electron tolerance.
-    return u * (8.0 * gaussian + (operations + x2) * correction_scale + n * std::abs(value))
-           + slope * x_error + tiny;
-}
+                       const std::vector<int>& isk);
 } // namespace
 
 Occupy::Occupy()
@@ -399,10 +321,11 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
                            const std::vector<int>& isk)
 {
     // Resolve particle number to floating-point accuracy, independently of the
-    // outer SCF threshold. A fixed absolute tolerance leaves a charge floor.
+    // outer SCF threshold. The modest 16-epsilon allowance is validated against
+    // independent occupation sums; it is not an all-spectra rounding bound.
     const double electron_magnitude = std::abs(nelec);
     const double scale = std::max(1.0, electron_magnitude);
-    const double tolerance = 8.0 * std::numeric_limits<double>::epsilon() * scale;
+    const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * scale;
     if (!(smearing_sigma > 0.0) || !std::isfinite(smearing_sigma) || !std::isfinite(nelec))
     {
         ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Invalid smearing width or electron number");
@@ -442,8 +365,8 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
     double step = 2.0 * smearing_sigma;
     lower -= step;
     upper += step;
-    double low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk, nullptr);
-    double high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk, nullptr);
+    double low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
+    double high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
     int expansion = 0;
     while (low_count > nelec + tolerance || high_count < nelec - tolerance)
     {
@@ -455,12 +378,12 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         if (low_count > nelec + tolerance)
         {
             lower -= step;
-            low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk, nullptr);
+            low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk);
         }
         if (high_count < nelec - tolerance)
         {
             upper += step;
-            high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk, nullptr);
+            high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, upper, is, isk);
         }
     }
 
@@ -480,8 +403,7 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         result.offset = upper;
         best_error = high_error;
     }
-    // Inherited safety cap, not an accuracy criterion. Stagnation is checked
-    // separately; exhausting this cap must never report convergence.
+    // Inherited safety cap; exhausting it must not report convergence.
     const int maxiter = 300;
     for (int iteration = 0; iteration < maxiter; ++iteration)
     {
@@ -504,30 +426,17 @@ FermiLevel find_fermi_level(const ModuleBase::matrix& ekb,
         const double middle = 0.5 * lower + 0.5 * upper;
         if (middle == lower || middle == upper)
         {
-            // Audit only the best candidate, in its original centered coordinates.
-            // A large error estimate must not buy a loose electron tolerance.
-            double roundoff = 0.0;
-            const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss,
-                                                 result.reference, result.offset, is, isk, &roundoff);
-            const double count_error = count - nelec;
-            const double residual = std::abs(count_error);
+            // No representable interior point remains. Do not silently spin to
+            // maxiter, nor claim that a small energy interval proves charge accuracy.
             std::ostringstream message;
             message.precision(17);
-            if (std::isfinite(roundoff) && roundoff <= tolerance && residual <= tolerance + roundoff)
-            {
-                message << "Fermi search reached rounding accuracy; electron residual = " << residual
-                        << ", target = " << tolerance << ", estimated roundoff = " << roundoff;
-                const std::string diagnostic = message.str();
-                ModuleBase::WARNING("occupy_smearing::efermig", diagnostic);
-                return result;
-            }
             message << "Fermi bracket reached floating-point resolution; electron residual = "
-                    << residual << ", requested tolerance = " << tolerance << ", estimated roundoff = " << roundoff;
+                    << best_error << ", requested tolerance = " << tolerance;
             const std::string diagnostic = message.str();
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", diagnostic);
             return result;
         }
-        const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, middle, is, isk, nullptr);
+        const double count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, middle, is, isk);
         if (!std::isfinite(count))
         {
             ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Non-finite electron number");
@@ -600,15 +509,11 @@ double sum_occupations(const ModuleBase::matrix& ekb,
                        const double reference,
                        const double offset,
                        const int& is,
-                       const std::vector<int>& isk,
-                       double* roundoff)
+                       const std::vector<int>& isk)
 {
     // ModuleBase::TITLE("Occupy","sumkg");
     double sum2 = 0.0;
     double compensation = 0.0;
-    const double u = 0.5 * std::numeric_limits<double>::epsilon();
-    double absolute_sum = 0.0;
-    double evaluation_error = 0.0;
     for (int ik = 0; ik < nks; ++ik)
     {
         if (is != -1 && is != isk[ik])
@@ -621,58 +526,24 @@ double sum_occupations(const ModuleBase::matrix& ekb,
             // local accumulation before the existing pool/band reduction.
             const double relative_energy = ekb(ik, ib) - reference;
             const double argument = (offset - relative_energy) / smearing_sigma;
-            const double value = occupy_smearing::wgauss(argument, ngauss);
-            const double occupation = wk[ik] * value;
+            const double occupation = wk[ik] * occupy_smearing::wgauss(argument, ngauss);
             const double corrected = occupation - compensation;
             const double next = sum2 + corrected;
             compensation = (next - sum2) - corrected;
             sum2 = next;
-            if (roundoff != nullptr)
-            {
-                const double magnitude = std::abs(occupation);
-                absolute_sum += magnitude;
-                const double difference = offset - relative_energy;
-                const double coordinate_scale = std::abs(relative_energy) / smearing_sigma
-                                                + std::abs(difference) / smearing_sigma + std::abs(argument);
-                const double x_error = u * coordinate_scale;
-                const double scalar_error = smearing_roundoff(argument, ngauss, value, x_error);
-                evaluation_error += std::abs(wk[ik]) * scalar_error + u * magnitude;
-            }
         }
     }
 
     // GlobalV::ofs_running << "\n sum2 before reduce = " << sum2 << std::endl;
 
-    double totals[3] = {sum2, absolute_sum, evaluation_error};
-    double reduction_error = 0.0;
 #ifdef __MPI
     const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-    const int count = roundoff == nullptr ? 1 : 3;
-    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, totals, count);
-    if (roundoff != nullptr && npool > 1)
-    {
-        int participants = 1;
-        Parallel_Reduce::reduce_all(participants);
-        const double additions = (participants - 1) * u;
-        // Pool scaling plus an order-independent reduction bound; used only
-        // for the capped stagnation check, never for early termination.
-        reduction_error = (u + additions / (1.0 - additions)) * totals[1];
-    }
+    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, sum2);
 #endif
-    if (roundoff != nullptr)
-    {
-        // First-order Kahan bound, function/argument evaluation, MPI and the
-        // final count subtraction. Subnormal widths invalidate this model.
-        *roundoff = totals[2] + 2.0 * u * totals[1] + reduction_error + u * std::abs(totals[0]);
-        if (smearing_sigma < std::numeric_limits<double>::min())
-        {
-            *roundoff = std::numeric_limits<double>::infinity();
-        }
-    }
 
     // GlobalV::ofs_running << "\n sum2 after reduce = " << sum2 << std::endl;
 
-    return totals[0];
+    return sum2;
 }
 } // namespace
 
@@ -687,7 +558,7 @@ double occupy_smearing::sumkg(const ModuleBase::matrix& ekb,
                              const std::vector<int>& isk)
 {
     const double reference = 0.0;
-    return sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, e, is, isk, nullptr);
+    return sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, e, is, isk);
 }
 
 double occupy_smearing::wgauss(const double& x, const int n)
diff --git a/source/source_estate/occupy.h b/source/source_estate/occupy.h
index d63abe53632..4e89f6b1102 100644
--- a/source/source_estate/occupy.h
+++ b/source/source_estate/occupy.h
@@ -40,11 +40,8 @@ double sumkg(const ModuleBase::matrix& ekb,
 /// @note ef is the rounded absolute chemical potential. Re-evaluating sumkg(ef)
 ///       can lose precision for a large energy offset; gweights retains the
 ///       centered chemical potential internally when forming occupations.
-///       The usual electron-count target is 8*epsilon*max(1,abs(nelec)). At
-///       floating-point stagnation only, a diagnosed rounding-limited result
-///       may exceed it by an estimated error no larger than the target itself.
-///       The estimate is a first-order libm/arithmetic model, not a certified
-///       physical-accuracy bound; a larger estimate still causes failure.
+///       The electron-count target is 16*epsilon*max(1,abs(nelec)). Unresolved
+///       floating-point stagnation remains an error; there is no relaxed fallback.
 void efermig(const ModuleBase::matrix& ekb,
              const int nbnd,
              const int nks,
diff --git a/source/source_estate/test/elecstate_occupy_test.cpp b/source/source_estate/test/elecstate_occupy_test.cpp
index dea764a43b4..9689d794a35 100644
--- a/source/source_estate/test/elecstate_occupy_test.cpp
+++ b/source/source_estate/test/elecstate_occupy_test.cpp
@@ -1,7 +1,6 @@
 #include <string>
 #include <cmath>
 #include <limits>
-#include <sstream>
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 #include "source_io/module_parameter/parameter.h"
@@ -323,38 +322,6 @@ TEST_F(OccupyTest, FractionalFermiDiracOccupation)
     EXPECT_NEAR(entropy, expected_entropy, 2e-15);
 }
 
-TEST_F(OccupyTest, WeightedParticleNumberAllSmearings)
-{
-    ModuleBase::matrix energies(3, 4);
-    const std::vector<double> weights{0.17, 0.29, 0.54};
-    const std::vector<int> spins{0, 0, 0};
-    for (int ik = 0; ik < 3; ++ik)
-    {
-        for (int ib = 0; ib < 4; ++ib)
-        {
-            energies(ik, ib) = -0.4 + 0.27 * ib + 0.07 * ik;
-        }
-    }
-    const double target = 1.23456789;
-    for (const int smearing : {0, 1, 2, -1, -99})
-    {
-        double ef = 0.0;
-        double entropy = 0.0;
-        ModuleBase::matrix occupations(3, 4);
-        occupy.gweights(3, weights, 4, target, 0.1, smearing, energies, ef, entropy, occupations, -1, spins);
-        long double count = 0.0L;
-        for (int ik = 0; ik < 3; ++ik)
-        {
-            for (int ib = 0; ib < 4; ++ib)
-            {
-                count += occupations(ik, ib);
-            }
-        }
-        EXPECT_NEAR(count, target, 3e-15) << smearing;
-        EXPECT_TRUE(std::isfinite(entropy));
-    }
-}
-
 TEST_F(OccupyTest, EmptyFullAndGappedSpectra)
 {
     ModuleBase::matrix energies(1, 2);
@@ -447,7 +414,6 @@ TEST_F(OccupyTest, FermiSearchAtRepresentableEnergyLimit)
     const std::string diagnostic = testing::internal::GetCapturedStdout();
     EXPECT_THAT(diagnostic, testing::HasSubstr("Fermi bracket reached floating-point resolution"));
     EXPECT_THAT(diagnostic, testing::HasSubstr("electron residual"));
-    EXPECT_THAT(diagnostic, testing::HasSubstr("estimated roundoff = inf"));
 }
 
 TEST_F(OccupyTest, RejectInvalidWidthAndUnbracketedCharge)
@@ -475,35 +441,27 @@ TEST_F(OccupyTest, FermiParticleNumberAcrossScales)
             energies(ik, ib) = -1.0 + 0.05 * ib + 0.013 * ik;
         }
     }
-    // Exercise efermig directly, independently of gweights' band-count fix.
+    // Check the returned occupations across electron-number scales and smearings.
     for (const double target : {0.3141592653589793, 33.0, 63.125})
     {
         for (const int smearing : {0, 1, 2, -1, -99})
         {
             double ef = 0.0;
-            occupy_smearing::efermig(energies, 64, 3, target, weights, 0.1, smearing, ef, -1, spins);
+            double entropy = 0.0;
+            ModuleBase::matrix occupations(3, 64);
+            occupy.gweights(3, weights, 64, target, 0.1, smearing, energies, ef, entropy, occupations, -1, spins);
             long double count = 0.0L;
             for (int ik = 0; ik < 3; ++ik)
             {
                 for (int ib = 0; ib < 64; ++ib)
                 {
-                    const double argument = (ef - energies(ik, ib)) / 0.1;
-                    const double occupation = weights[ik] * occupy_smearing::wgauss(argument, smearing);
-                    count += occupation;
+                    count += occupations(ik, ib);
                 }
             }
             const double scale = std::max(1.0, target);
-            // Allow the independently accumulated result a small rounding margin.
-            const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * scale;
+            // Add one epsilon for the independent sum of the returned doubles.
+            const double tolerance = 17.0 * std::numeric_limits<double>::epsilon() * scale;
             EXPECT_NEAR(count, target, tolerance) << "N = " << target << ", smearing = " << smearing;
-            if (target == 33.0 && smearing == 0)
-            {
-                const double residual = std::abs(static_cast<double>(count - target));
-                std::ostringstream value;
-                value.precision(17);
-                value << residual;
-                RecordProperty("N33_gaussian_residual", value.str());
-            }
         }
     }
 }
@@ -531,50 +489,6 @@ TEST_F(OccupyTest, GweightsHonorsExplicitBandRange)
     EXPECT_DOUBLE_EQ(occupations(0, 2), -123.0);
 }
 
-TEST_F(OccupyTest, DegenerateSpectrumEnergyShift)
-{
-    // Reproduce ElecStateTest.CalculateWeightsGWeights without changing its fixture.
-    const int points = 5;
-    const int bands = 6;
-    ModuleBase::matrix energies(points, bands);
-    ModuleBase::matrix occupations(points, bands);
-    const std::vector<double> weights(points, 2.0);
-    const std::vector<int> spins(points, 0);
-    const double target = 10.0;
-    const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * target;
-    for (const double shift : {0.0, 100.0, 200.0, 1e8})
-    {
-        for (int ik = 0; ik < points; ++ik)
-        {
-            for (int ib = 0; ib < bands; ++ib)
-            {
-                energies(ik, ib) = shift;
-            }
-        }
-        double ef = 0.0;
-        double entropy = 0.0;
-        occupy.gweights(points, weights, bands, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
-        long double count = 0.0L;
-        for (int ik = 0; ik < points; ++ik)
-        {
-            for (int ib = 0; ib < bands; ++ib)
-            {
-                EXPECT_NEAR(occupations(ik, ib), 1.0 / 3.0, 2e-15);
-                count += occupations(ik, ib);
-            }
-        }
-        EXPECT_NEAR(count, target, tolerance) << shift;
-        if (shift == 100.0)
-        {
-            const double residual = std::abs(static_cast<double>(count - target));
-            std::ostringstream value;
-            value.precision(17);
-            value << residual;
-            RecordProperty("shift100_particle_residual", value.str());
-        }
-    }
-}
-
 TEST_F(OccupyTest, OccupationAndEntropyEnergyShift)
 {
     ModuleBase::matrix energies(1, 3);
@@ -587,7 +501,7 @@ TEST_F(OccupyTest, OccupationAndEntropyEnergyShift)
     {
         double reference_entropy = 0.0;
         double reference_fermi = 0.0;
-        for (const double shift : {0.0, 100.0, 200.0, 1e8})
+        for (const double shift : {0.0, 1e8})
         {
             for (int ib = 0; ib < 3; ++ib)
             {
@@ -630,7 +544,8 @@ TEST_F(OccupyTest, NarrowSmearingAtLargeEnergyOffset)
     double ef = 0.0;
     double entropy = 0.0;
     occupy.gweights(1, weights, 1, 0.3, 1e-8, 0, energies, ef, entropy, occupations, -1, spins);
-    EXPECT_NEAR(occupations(0, 0), 0.3, 2e-15);
+    const double tolerance = 16.0 * std::numeric_limits<double>::epsilon();
+    EXPECT_NEAR(occupations(0, 0), 0.3, tolerance);
     EXPECT_TRUE(std::isfinite(ef));
     EXPECT_TRUE(std::isfinite(entropy));
 }
@@ -643,51 +558,52 @@ TEST_F(OccupyTest, FermiSurfaceAwayFromSpectralMidpoint)
     const std::vector<int> spins{0};
     const double target = 10.0;
     const double tolerance = 16.0 * std::numeric_limits<double>::epsilon() * target;
-    for (const double shift : {0.0, 100.0, 1e8})
+    for (int ib = 0; ib < 6; ++ib)
     {
-        for (int ib = 0; ib < 6; ++ib)
-        {
-            energies(0, ib) = shift;
-        }
-        energies(0, 6) = shift + 100.0;
-        double ef = 0.0;
-        double entropy = 0.0;
-        occupy.gweights(1, weights, 7, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
-        long double count = 0.0L;
-        for (int ib = 0; ib < 7; ++ib)
-        {
-            count += occupations(0, ib);
-        }
-        EXPECT_NEAR(count, target, tolerance) << shift;
-        EXPECT_DOUBLE_EQ(occupations(0, 6), 0.0);
+        energies(0, ib) = 0.0;
     }
+    energies(0, 6) = 100.0;
+    double ef = 0.0;
+    double entropy = 0.0;
+    occupy.gweights(1, weights, 7, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+    long double count = 0.0L;
+    for (int ib = 0; ib < 7; ++ib)
+    {
+        count += occupations(0, ib);
+    }
+    EXPECT_NEAR(count, target, tolerance);
+    EXPECT_DOUBLE_EQ(occupations(0, 6), 0.0);
 }
 
-TEST_F(OccupyTest, ReviewerDegenerateGaussianSpectrum)
+TEST_F(OccupyTest, LowFillingDegenerateSpectrum)
 {
-    // PR #8130: the former 1-erf evaluation stalls at 3.99680288865e-15.
-    const int bands = 74;
-    ModuleBase::matrix energies(1, bands);
-    ModuleBase::matrix occupations(1, bands);
+    // Gaussian lower-tail cancellation and MP1 cancellation near zero filling.
+    struct Case
+    {
+        int bands;
+        int smearing;
+        double target;
+    };
+    const Case cases[] = {{74, 0, 2.0}, {26, 1, 0.3}};
     const std::vector<double> weights{2.0};
     const std::vector<int> spins{0};
-    const double target = 2.0;
-    const double tolerance = 8.0 * std::numeric_limits<double>::epsilon() * target;
-    for (const double shift : {0.0, 100.0, 1e8})
+    for (const Case& test : cases)
     {
-        for (int ib = 0; ib < bands; ++ib)
-        {
-            energies(0, ib) = shift;
-        }
+        ModuleBase::matrix energies(1, test.bands);
+        ModuleBase::matrix occupations(1, test.bands);
         double ef = 0.0;
         double entropy = 0.0;
-        occupy.gweights(1, weights, bands, target, 0.01, 0, energies, ef, entropy, occupations, -1, spins);
+        occupy.gweights(1, weights, test.bands, test.target, 0.01, test.smearing,
+                        energies, ef, entropy, occupations, -1, spins);
         long double count = 0.0L;
-        for (int ib = 0; ib < bands; ++ib)
+        for (int ib = 0; ib < test.bands; ++ib)
         {
             count += occupations(0, ib);
         }
-        EXPECT_NEAR(count, target, tolerance) << shift;
+        const double scale = std::max(1.0, test.target);
+        // Include one epsilon for the independent accumulation.
+        const double tolerance = 17.0 * std::numeric_limits<double>::epsilon() * scale;
+        EXPECT_NEAR(count, test.target, tolerance) << test.smearing;
         EXPECT_TRUE(std::isfinite(entropy));
     }
 }
@@ -703,30 +619,12 @@ TEST_F(OccupyTest, SmearingHighPrecisionReference)
         double occupation;
     };
     const Reference references[] = {
-        {0, -6.0, 1.0759868356249456558e-17},
         {0, -8.0, 5.61214858649146354e-30},
-        {-1, -6.0, 1.2797773847197835025e-20},
         {-1, -8.0, 5.1183964952735562742e-34},
-        {0, -2.0, 0.002338867490523632918965},
-        {0, -0.5, 0.2397500610934767311586},
-        {0, 0.5, 0.7602499389065232688414},
-        {0, 2.0, 0.997661132509476367081},
         {1, -2.0, -0.007994625186522394009569},
-        {1, -0.5, 0.1299022387265461318969},
-        {1, 0.5, 0.8700977612734538681031},
-        {1, 2.0, 1.00799462518652239401},
-        {2, -2.0, 0.004922240659785139651099},
         {2, -0.5, 0.06124734974721450735834},
-        {2, 0.5, 0.9387526502527854926417},
-        {2, 2.0, 0.9950777593402148603489},
-        {-1, -2.0, 0.0003264301114214816118753},
-        {-1, -0.5, 0.1368184140991744579818},
-        {-1, 0.5, 0.766994049270187616701},
         {-1, 2.0, 1.041238947228640924942},
-        {-99, -2.0, 0.1192029220221175559403},
-        {-99, -0.5, 0.3775406687981454353611},
-        {-99, 0.5, 0.6224593312018545646389},
-        {-99, 2.0, 0.8807970779778824440597}};
+        {-99, -0.5, 0.3775406687981454353611}};
     for (const Reference& reference : references)
     {
         const double value = occupy_smearing::wgauss(reference.x, reference.order);