Add a Newton chemical-potential search (mu_search='newton')

[k-yoshimi]

Summary

The chemical-potential adjustment re-evaluates the (expensive) lattice density n(μ) at every step of the existing bracketing + Brent search (calc_mu). Since n(μ) is smooth and monotonic and its derivative is almost free, a Newton iteration usually needs fewer of those density evaluations.

n(μ)    = Σ_k w_k Σ_block [ (1/β) Σ_iω Tr G + 0.5 n_orb ]
dn/dμ   = -(1/β) Σ_k w_k Σ_block Σ_iω Tr[G²]        (dG/dμ = -G², tail term is μ-independent)

• total_density_and_derivative() — returns (n, dn/dμ) in a single k-loop (Matsubara-only; guards iw_or_w).
• calc_mu_newton() — a Numerical-Recipes rtsafe safeguarded Newton: brackets the root (both directions), takes a Newton step only when it stays in the bracket and reduces the residual fast enough, else bisects. In a charge gap (dn/dμ ≈ 0) it falls back to bisection, so it is never less robust than bisection. Returns None (like calc_mu) on failure.
• New [system] mu_search option ('brent' default, or 'newton'), plumbed through dmft_core to SumkDFTWorkerGloc. 'newton' adjusts μ from the Matsubara-summed charge, so it requires no_tail_fit=True (enforced where the search actually runs).

Why the analytic derivative is exact here

The μ-search is the standard DMFT inner step at fixed Σ, so dn/dμ has no ∂Σ/∂μ term and is exact. Both methods solve the identical n(μ; Σ_fixed) − N = 0, so they find the same μ — Newton only changes the path. Σ's μ-dependence is captured by the outer self-consistency.

Measured effect

Square model, n_orb=24, no_tail_fit=True, one Gloc evaluation:

method	density evaluations	wall time	total charge
brent	6	19.8 s	24.000000
newton	1	8.4 s	24.000000

The gain is largest near self-consistency, where μ barely moves between iterations (Newton's early return).

Tests

tests/non-mpi/mu_newton covers the root finder (smooth, charge-gap, wrong-signed initial derivative, already-converged, non-convergence→None), the iw-only guard, parameter defaulting/validation, and the derivative vs a finite difference. The default mu_search='brent' leaves existing behavior unchanged (wannier90 regression tests pass).

🤖 Generated with Claude Code

[k-yoshimi]

Full-run validation: newton vs brent across all iterations

To confirm mu_search='newton' is physically equivalent to 'brent' over a complete self-consistency (not just one Gloc evaluation), I ran an interacting, off-half-filling DMFT and compared every iteration.

• Bethe lattice, norb=1, nelec=0.7, U=4, beta=20, n_iw=1000, prec_mu=1e-4, no_tail_fit=True, sigma_mix=0.5, 6 iterations
• Solver: scipy/sparse (exact diagonalization, n_bath=2) — so Σ is non-trivial and evolves each iteration

iter	μ (brent)	μ (newton)	|Δμ|	max|ΔΣ|
1	−0.47790644	−0.47790658	1.4e-07	7.1e-08
2	0.02806586	0.02806710	1.2e-06	3.2e-07
3	0.13191055	0.13191073	1.8e-07	1.8e-07
4	0.20703274	0.20703290	1.6e-07	1.3e-07
5	0.26180316	0.26180335	1.9e-07	1.8e-07
6	0.30191163	0.30191187	2.4e-07	1.8e-07

Both μ and Σ agree at every iteration to ~1e-7–1e-6, well below prec_mu=1e-4, and the difference does not accumulate over iterations — the two searches follow the same self-consistent trajectory. This is expected: each iteration solves the same fixed-Σ equation n(μ;Σ)=N, so brent and newton converge to the same root (to within the requested tolerance); tightening prec_mu shrinks the difference further.