diff --git a/source/source_estate/occupy.cpp b/source/source_estate/occupy.cpp index 8fd2898b9c..c573052fe1 100644 --- a/source/source_estate/occupy.cpp +++ b/source/source_estate/occupy.cpp @@ -1,10 +1,45 @@ #include "occupy.h" +#include +#include + #include "source_base/constants.h" #include "source_base/mymath.h" #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& wk, + const double& smearing_sigma, + const int ngauss, + const int& is, + const std::vector& isk); + +double sum_occupations(const ModuleBase::matrix& ekb, + const int nband, + const int nks, + const std::vector& wk, + const double& smearing_sigma, + const int ngauss, + const double reference, + const double offset, + const int& is, + const std::vector& isk); +} // namespace + Occupy::Occupy() { } @@ -201,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 @@ -224,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++) @@ -235,13 +271,15 @@ 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 //================================ // 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& wk, - const double& smearing_sigma, - const int ngauss, - double& ef, - const int& is, - const std::vector& isk) +namespace { - // 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; +FermiLevel find_fermi_level(const ModuleBase::matrix& ekb, + const int nband, + const int nks, + const double& nelec, + const std::vector& wk, + const double& smearing_sigma, + const int ngauss, + const int& is, + const std::vector& isk) +{ + // Resolve particle number to floating-point accuracy, independently of the + // 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 = 16.0 * std::numeric_limits::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"); + } - /* - for(int ik=0; ik::infinity(); + double upper = -std::numeric_limits::infinity(); + for (int ik = 0; ik < nks; ++ik) { - for(int i=0; i 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) + if (++expansion > 100 || !std::isfinite(step)) { - 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 + ModuleBase::WARNING_QUIT("occupy_smearing::efermig", "Cannot bracket the requested electron number"); } - else if ((sumkup - nelec) < -eps) + step *= 2.0; + if (low_count > nelec + tolerance) { - eup += 2 * smearing_sigma; - ++changetime; + lower -= step; + low_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, lower, is, isk); } - else if ((sumklw - nelec) > eps) + if (high_count < nelec - tolerance) { - elw -= 2 * smearing_sigma; - ++changetime; - } - else - { - break; + upper += step; + high_count = sum_occupations(ekb, nband, nks, wk, smearing_sigma, ngauss, reference, 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. + 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) + { + result.offset = upper; + best_error = high_error; + } + // Inherited safety cap; exhausting it must not report convergence. + const int maxiter = 300; + for (int iteration = 0; iteration < maxiter; ++iteration) + { + if (best_error <= tolerance) + { + 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) + { + // 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 result; + } + 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"); + } + const double error = std::abs(count - nelec); + if (error < best_error) { - return; + result.reference = reference; + result.offset = middle; + best_error = error; } - else if ((sumkmid - nelec) < -eps) + if (best_error <= tolerance) { - elw = ef; + return result; + } + if (count < nelec) + { + lower = middle; } else { - eup = ef; + upper = middle; } } - return; + 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& wk, + const double& smearing_sigma, + const int ngauss, + double& ef, + const int& is, + const std::vector& 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; } /** @@ -390,38 +492,46 @@ 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& wk, - const double& smearing_sigma, - const int ngauss, - const double& e, - const int& is, - const std::vector& isk) +namespace +{ +double sum_occupations(const ModuleBase::matrix& ekb, + const int nband, + const int nks, + const std::vector& wk, + const double& smearing_sigma, + const int ngauss, + const double reference, + const double offset, + const int& is, + const std::vector& isk) { // 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 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; + sum2 = next; } - sum2 += wk[ik] * sum1; } // GlobalV::ofs_running << "\n sum2 before reduce = " << sum2 << std::endl; @@ -435,6 +545,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& wk, + const double& smearing_sigma, + const int ngauss, + const double& e, + const int& is, + const std::vector& 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) { @@ -481,14 +606,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="< 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 +289,346 @@ TEST_F(OccupyTest, Gweights) std::vector 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 weights{1.0}; + const std::vector 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, EmptyFullAndGappedSpectra) +{ + ModuleBase::matrix energies(1, 2); + energies(0, 0) = -1.0; + energies(0, 1) = 1.0; + const std::vector weights{1.0}; + const std::vector 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 weights{1.0, 1.0}; + const std::vector 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 weights(points, 1.0 / points); + const std::vector 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) = 0.0; + const std::vector weights{1.0}; + const std::vector spins{0}; + double ef = 0.0; + // 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::denorm_min(); + testing::internal::CaptureStdout(); + 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")); + EXPECT_THAT(diagnostic, testing::HasSubstr("electron residual")); +} + +TEST_F(OccupyTest, RejectInvalidWidthAndUnbracketedCharge) +{ + ModuleBase::matrix energies(1, 1); + energies(0, 0) = 0.0; + const std::vector weights{1.0}; + const std::vector 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 weights{0.17, 0.29, 0.54}; + const std::vector 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; + } + } + // 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; + 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) + { + count += occupations(ik, ib); + } + } + const double scale = std::max(1.0, target); + // Add one epsilon for the independent sum of the returned doubles. + const double tolerance = 17.0 * std::numeric_limits::epsilon() * scale; + EXPECT_NEAR(count, target, tolerance) << "N = " << target << ", smearing = " << smearing; + } + } +} + +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 weights{1.0}; + const std::vector 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); +} + +TEST_F(OccupyTest, OccupationAndEntropyEnergyShift) +{ + ModuleBase::matrix energies(1, 3); + ModuleBase::matrix occupations(1, 3); + ModuleBase::matrix reference(1, 3); + const std::vector weights{2.0}; + const std::vector 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, 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 weights{1.0}; + const std::vector 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); + const double tolerance = 16.0 * std::numeric_limits::epsilon(); + EXPECT_NEAR(occupations(0, 0), 0.3, tolerance); + 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 weights{2.0}; + const std::vector spins{0}; + const double target = 10.0; + const double tolerance = 16.0 * std::numeric_limits::epsilon() * target; + for (int ib = 0; ib < 6; ++ib) + { + 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, LowFillingDegenerateSpectrum) +{ + // 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 weights{2.0}; + const std::vector spins{0}; + for (const Case& test : cases) + { + ModuleBase::matrix energies(1, test.bands); + ModuleBase::matrix occupations(1, test.bands); + double ef = 0.0; + double entropy = 0.0; + 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 < test.bands; ++ib) + { + count += occupations(0, ib); + } + const double scale = std::max(1.0, test.target); + // Include one epsilon for the independent accumulation. + const double tolerance = 17.0 * std::numeric_limits::epsilon() * scale; + EXPECT_NEAR(count, test.target, tolerance) << test.smearing; + 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, -8.0, 5.61214858649146354e-30}, + {-1, -8.0, 5.1183964952735562742e-34}, + {1, -2.0, -0.007994625186522394009569}, + {2, -0.5, 0.06124734974721450735834}, + {-1, 2.0, 1.041238947228640924942}, + {-99, -0.5, 0.3775406687981454353611}}; + for (const Reference& reference : references) + { + const double value = occupy_smearing::wgauss(reference.x, reference.order); + const double tolerance = 32.0 * std::numeric_limits::epsilon() * std::abs(reference.occupation); + EXPECT_NEAR(value, reference.occupation, tolerance) << reference.order << ", x = " << reference.x; + } +}