Make method_init_density='pp' work for psp8 and add a generic atomic-density fallback

[MZKC]

Summary

method_init_density='pp' builds the initial electron density from the
pseudopotential's atomic valence density (rho_pp_tbl). That density is smooth
and positive everywhere, unlike the default 'wf' guess that is built from the
gaussian initial orbitals — whose near-zero regions make meta-GGA functionals
(e.g. TBmBJ) diverge on the first iterations. So 'pp' is the natural,
non-pathological initial density (cf. issue #1253).

Two gaps are fixed so that 'pp' works for more pseudopotentials:

1. psp8 stores the pseudo valence density in the file, but the reader read
   it only to print a diagnostic and then discarded it. It is now stored
   into rho_pp_tbl as 4*pi*r^2*rho (= r^2 * rho_tmp, since psp8 stores
   4*pi*rho, the same convention as its NLCC data).
2. Generic fallback: if a pseudopotential still provides no tabulated atomic
   density, 'pp' previously aborted with "radial density is not available".
   It now builds a generic atomic density from the valence charge
   (set_generic_atomic_density: an exponential model normalised to pp%zps)
   instead of aborting.

'wf'/'read_dns_cube' and UPF (which already reads <PP_RHOATOM>) are
unaffected.

Verification (Wisteria-O, A64FX)

1. psp8 'pp' works and is correctly normalised. bulk Si (8 atoms,
Si.psp8 ONCVPSP, zion=4, so 32 valence electrons), nproc_k=4:

• 'pp' no longer aborts; calc_density_pp reports Int(rho) = 32.0000000 —
  i.e. the initial density integrates to exactly the number of valence
  electrons, confirming the normalisation.
• the GS converges (iter 60: total energy -951.4 eV, gap +0.03 eV). 'pp' is in
  fact a better starting guess here than 'wf' (which still had a negative gap
  at iter 60).

2. 'pp' cures the meta-GGA initial-density pathology (issue #1253). bulk Si
(8 atoms, Si_rps.dat), xc='TBmBJ':

method_init_density	result
'wf' (gaussian)	diverges — iter 4 total energy +1580 eV, gap overflows
'pp'	converges — total energy -862.8 eV, gap 2.64 eV

🤖 Generated with Claude Code

[jenkins-diana]

Can one of the admins verify this patch?

[MZKC]

psp8 valence-density support and the generic fallback look correct: the psp8 reader runs before the empty-table check, so the fallback only fires when the file genuinely has no atomic density, and the exponential model integrates to zps. One robustness concern on the storage bounds, inline.

[MZKC]

if (i <= pp%nrmax) pp%rho_pp_tbl(i,ik) = pp%rad(i,ik)**2 * rho_tmp

The i <= pp%nrmax guard avoids an out-of-bounds write, but it does so by silently truncating the valence density when the psp8 grid has more points than pp%nrmax: the loop runs i=1..pp%mr(ik)+1, yet sum_rho_pp on the next line still accumulates all of them, so the stored table and the reported norm would disagree. The tail-fill just below (do i=pp%mr(ik)+2,pp%nrmax) and the other readers assume pp%mr(ik)+1 <= pp%nrmax. Consider asserting that explicitly (error out on an oversized grid) rather than partial storage, so the case is caught instead of silently dropped.

[MZKC]

Fixed in b29617d: the reader now errors out when pp%mr(ik)+1 > pp%nrmax instead of silently truncating (which would have left the stored density inconsistent with sum_rho_pp). Builds cleanly on Fugaku (A64FX).

---

✅ Verified on Fugaku (A64FX), 2026-06-24: a psp8 pseudopotential (Be) with method_init_density='pp' reads the valence density and runs the SCF, confirming the psp8 initial-density path is intact after the change.