Fix smeared-occupation Fermi convergence and particle-number accuracy

[aboys-cb]

Reminder

• Read AGENTS.md and docs/developers_guide/agent_governance.md.
• Included reproducing tests, verification commands, behavior changes and core-module impact.
• Explained INPUT/documentation scope and governance findings.

Linked Issue

No separate issue; numerical reproductions are included in the tests. Base: 457d9dd4a4d1715d4d0c772a77d33f08d7e36ca9.

What's changed?

The smeared-occupation Fermi search could return after exhausting its iterations without reaching the requested electron number. Its fixed 1e-10 charge tolerance also limited higher-precision calculations, and gweights used different band counts for searching and writing occupations.

This PR uses the explicit band count consistently, selected-spin eigenvalue bounds, compensated occupation sums, and centered Fermi coordinates retained when evaluating occupations and entropy. Gaussian/MP and cold-smearing bases use erfc to avoid lower-tail cancellation. The public absolute ef remains rounded to one double.

Following review, the convergence rule is now simply:

abs(computed_count - nelec) <= 16 * epsilon * max(1, abs(nelec))

All runtime error-estimation code, its extra MPI reductions and the separate stagnation fallback have been removed. Unresolved stagnation, invalid inputs, unbracketable targets and exhaustion of the named maxiter = 300 remain failures. The factor 16 is a validated allowance, not a universal error bound for arbitrary spectra or MPI layouts.

Tests now retain one case per distinct failure mechanism: the duplicate weighted-count and degenerate-shift tests were removed, repeated shifts were trimmed, and scalar references were reduced to six cases covering tails, signed MP values and overshoot. LowFillingDegenerateSpectrum describes the numerical case and covers both the 74-level Gaussian reproduction and a 26-level MP1 cancellation case. The unchanged ElecStateTest.CalculateWeightsGWeights still covers the original 100-Ry CI regression.

Unit Tests and/or Case Tests for my changes

Fresh validation with OMP_NUM_THREADS=1:

Check	Result
Occupation suite, GCC 13.3 / GoogleTest 1.14	38/38 passed
Entire MODULE_ESTATE group	19/19 passed
Serial 16-epsilon sweep, macOS/Clang and Linux/GCC	1,442/1,442 passed on each platform
Identical sweep with only the tolerance changed to 8 epsilon	macOS: 1,441/1,442; Linux: 1,442/1,442
Actual MPI implementation, ranks/pools = 1/1, 4/2, 56/2, 56/56	120/120 passed, including replicated pools and empty partitions

The serial sweep covers five smearings, 1/3 weighted k points, 1/26/74/256 bands, 5%/30%/90% filling, widths 0.01/1e-8, energy shifts 0/100/1e8, and degenerate/dispersed spectra, plus the two low-filling regressions. On macOS, the 8-epsilon MP1 case stalls at a computed residual of 2.387e-15 against 1.776e-15; 16 epsilon resolves it without a fallback.

Independent 80-digit summation of the returned double occupations gives a maximum residual of 16.114 epsilon times max(1, |N|) across both platforms. Independent 80-digit evaluation of the analytic smearing functions in the saved centered coordinates gives 19.672 epsilon times max(1, |N|). These distinguish the computed stopping criterion from function-evaluation/accumulation error. The independent-sum tests allow one extra epsilon; the external analytic-reference check uses a preselected 20-epsilon limit. Neither margin changes the solver threshold.

The narrow-width/large-offset test's old 2e-15 assertion was updated to the new 16*epsilon contract: its observed residual is 3.053e-15, below 3.553e-15. The subnormal-width failure test remains, with the removed budget text assertion deleted. Other analytic occupation/entropy checks retain their existing limits.

Commands and validation scope

The existing Sai build configuration was reused: Release, GNU 13.3, OpenMPI 5.0.10, CMake 3.31.6, CUDA sm70, LCAO/ELPA/Libxc and testing enabled. Changed sources were rebuilt and SHA256-checked against the submitted files.

export OMP_NUM_THREADS=1
ctest --test-dir "$BUILD" --show-only=json-v1 -R MODULE_ESTATE
# Build all 19 target names returned above:
cmake --build "$BUILD" -j8 --target <the 19 MODULE_ESTATE targets>
ctest --test-dir "$BUILD" --output-on-failure -R MODULE_ESTATE
"$BUILD/source/source_estate/test/MODULE_ESTATE_elecstate_occupy" --gtest_output="xml:$OUT/fixed.xml"

# Standalone serial probes extract the current solver unchanged,
# apart from the explicit 8-versus-16 control; compile with Clang and GCC:
clang++ -std=c++11 -O2 -I"$SRC/source" "$OUT/probe16.cpp" -o "$OUT/probe16"
g++ -std=c++11 -O2 -I"$SRC/source" "$OUT/probe16.cpp" -o "$OUT/probe16"
"$OUT/probe16" < "$OUT/probe-input.txt"
# Repeat with probe8.cpp; verify_counts.py uses mpmath at 80 digits.

mpirun --oversubscribe --mca pml ob1 --mca btl self,sm,tcp -np 1 "$OUT/mpi_probe_lean" 1
mpirun --oversubscribe --mca pml ob1 --mca btl self,sm,tcp -np 4 "$OUT/mpi_probe_lean" 2
mpirun --oversubscribe --mca pml ob1 --mca btl self,sm,tcp -np 56 "$OUT/mpi_probe_lean" 2
mpirun --oversubscribe --mca pml ob1 --mca btl self,sm,tcp -np 56 "$OUT/mpi_probe_lean" 56

git diff --cached --check
python3 tools/03_code_analysis/agent_governance_check.py --staged --event-path "$OUT/pr-event.json"
python3 tools/03_code_analysis/code_quality_score.py source/source_estate/occupy.cpp source/source_estate/occupy.h source/source_estate/test/elecstate_occupy_test.cpp

MPI probes include the current occupation implementation, link the real pool/reduction wrappers, and independently evaluate closed-form occupations in long double. All layouts pass; their maximum reference residual is 19.663 epsilon times N. These are single-node correctness probes with oversubscription, not performance or multi-node results. The first 56-rank probe, linked to the full CUDA-enabled test dependencies, exceeded the allocation's 15-GiB limit; the successful probe links only the required occupation/MPI sources and uses CPU transports. No product change was needed for that harness failure.

Limits: no new full SCF/force/D campaign or guarantee for arbitrarily ill-conditioned spectra, other libm implementations or larger MPI deployments. CUDA-enabled build/test success is not CPU/GPU physical-accuracy validation. External sweeps/probes are validation artifacts, not additional production code or repetitive committed tests.

Governance Notes

• INPUT/docs: No documentation update required for parameter metadata; no INPUT keywords, defaults or public signatures change. The numerical criterion and rounded absolute-output limitation are documented in the header.
• Core impact: ElecState Fermi search and occupation evaluation. Entropy/force/XC formulas and MPI wrappers are unchanged. No new globals, public interfaces or direct MPI calls; the global dependency budget decreases by one versus the PR base.
• Exceptions requested: none. The remaining non-increasing global-dependency move uses the existing pool-size dependency for the original scalar reduction.
• File-wide quality advisory: source 22 (previous revision 7), header 70, test 86. The source remains below the advisory 60-point line; this is disclosed rather than reported as a passing quality check.

[mohanchen]

there are still some issues I need to consider

[AsTonyshment]

Thanks for the improvement. However, I still have a couple of concerns noted inline.

[AsTonyshment]

Thanks for the updates. However, I think nowadays AI agents tend to overengineer solutions (e.g., long, repetitive tests), so we should weigh the added complexity against its actual benefit.

To clarify my main points:

1. I still think modestly relaxing the main tolerance, for example to twice its current value, would be simpler and more robust than adding a separate fallback that seems unlikely to accept a result (and hard to maintain!!).
2. I would also suggest trimming the tests to remove redundant cases.

[AsTonyshment]

Thanks! Now it looks much cleaner!