dav_subspace aborts on a small correction direction after normalization (FeRh GPU reproducer) #8127

Describe the bug

A noncollinear constrained-moment FeRh PW calculation exits with Diago_DavSubspace::cal_grad: psi_norm <= 0 for band 0 late in SCF. The diagnostic suggests npwx < nbands, but the recorded values are npwx = 3294 and NBANDS = 53.

The failure was observed with official develop commit f8118a0de8c8818c7c35f32ed443c18ad70a486d on one V100 GPU. Source inspection of develop at cd2c7ce5b2c8b24893d3426922a20fb7d4b7fd3c confirms that the relevant check/order is still present; that newer commit has not been runtime-tested for this report.

The last printed SCF row is DS46, with DRHO = 5.1186e-11 and EDIFF = -1.24903194e-9 eV, followed by:

Diago_DavSubspace::cal_grad: psi_norm <= 0 for band 0
This may be due to npwx < nbands: the number of plane waves is less than
the number of bands, leading to a rank-deficient problem.
Please increase ecutwfc or reduce nbands.

Exit status is 1. The requested total-energy criterion (scf_ene_thr = 1e-10 eV) was not satisfied; I am not suggesting accepting the whole SCF as converged.

Expected behavior

Check a correction direction before dividing by its norm. Distinguish an already converged eigenpair from an unusable direction or genuine breakdown, and retain the outer SCF stopping criteria. When termination is necessary, report the measured norm, threshold, physical band/k-point and residual information instead of suggesting an insufficient basis without checking it.

To Reproduce

Reproducer archive prepared: ferh-dav-subspace-reproducer.zip (389 KB). The initial attachment URL failed the public download check, so archive attachment is pending. Exact text inputs are included below in the meantime.

The archive contains original INPUT/STRU/KPT, both pseudopotentials, both orbital files, a run wrapper, effective INPUT, stdout/running/warning logs, commit/binary provenance and SHA256 checksums. All seven calculation input files match the recorded failing run byte-for-byte. Archive SHA256:

41bdb9c40aa9a82a7d698f6e0bd0d8376c3ef162a88d9e9baa7f0c9d3b8a2492

1. Build a CUDA-enabled ABACUS from the tested commit above, using your normal supported toolchain.
2. Extract the archive and allocate one GPU.
3. From the extracted directory:

sha256sum -c SHA256SUMS
export ABACUS_BIN=/absolute/path/to/abacus
bash case/run.sh

The wrapper runs OMP_NUM_THREADS=1 mpirun -np 1 and records stdout and exit status. This is the observed failing case, not yet a minimized reproducer; no new calculation was launched while packaging it.

Key settings: two-atom FeRh (3.06 angstrom cell), noncollinear/SOC with all moment components constrained, PW/PBE (gga_grad=2), 100 Ry, Gamma-centered 4x4x4 mesh, Gaussian 0.01 Ry, dav_subspace, pw_diag_ndim=2, pw_diag_nmax=20, pw_seed=1, scf_thr=1e-10, scf_ene_thr=1e-10 eV and sc_thr=2e-10. Full parameters are included.

Environment

• Linux HPC environment; NVIDIA Tesla V100-SXM2-16GB, one GPU, one MPI rank, one OpenMP thread.
• Runtime modules: NVHPC 26.3 / CUDA 12 toolchain, OpenMPI 5.0.10, FFTW 3.3.10, Libxc 7.0.0, ELPA 2026.02.001 and site BLAS.
• Executable SHA256: 004371ca6b070e985ac639beeaa3d627a0269ce6baa0a2e7e2069216c5d49324.
• GPU runtime failure is recorded. CPU behavior below is source inspection, not a reproduced CPU failure.

Additional Context

What the code checks

Here psi_norm (denoted by $\nu_n$ below) is the norm of a preconditioned residual/search direction, conceptually

$$ r_n = H\psi_n-\epsilon_n S\psi_n,\qquad p_n=P^{-1}r_n,\qquad \nu_n=\lVert p_n\rVert. $$

It is not the norm of the physical occupied wavefunction. An already converged eigenpair can have a vanishing correction direction, so a small norm alone does not establish a rank-deficient PW basis.

Pinned source evidence:

• cal_grad, CPU/GPU normalization and guards: normalize_op runs before psi_norm <= 1.0e-12 is tested. The message says <= 0, although a small positive value also triggers it.
• CUDA kernel: divides the vector by norm before returning that norm to the caller.
• CPU kernel: similarly normalizes before the caller's guard; an earlier assertion on the squared norm is build-dependent.
• diag_once ordering: all bands initially remain unconverged, and the first cal_grad precedes the iterative eigenvalue-change convergence check. This permits attempting to expand an eigenpair already converged by the initial reduced solve.

The actual norm is not logged, so an already converged eigenpair is a hypothesis, not an established root cause of this FeRh run. The evidence excludes the particular npwx < nbands suggestion, not every possible numerical breakdown.

A possible fix should separate norm measurement/finite checks from normalization, assess the unpreconditioned eigen-residual with appropriate tolerances, and handle converged or dependent directions through convergence/deflation/restart as appropriate. Simply lowering the constant or removing the guard would not resolve the ordering problem. Useful CPU/GPU tests include an exactly converged initial Ritz pair, a tiny positive correction, zero/nonfinite directions and this real case.

Related: #4919 reported pre-loop Davidson breakdown for an older solver and is closed; #7727 discusses broader GPU/paged-memory Diago_DavSubspace issues. This report supplies a regular GPU FeRh case and the current norm-check ordering evidence.

INPUT_PARAMETERS

device          gpu
gga_grad        2

suffix          ABACUS
calculation     scf
ecutwfc         100
scf_thr_type    2
scf_thr 1e-10
scf_ene_thr 1e-10
scf_nmax 200

smearing_method gaussian
smearing_sigma  0.01

mixing_type     broyden
mixing_beta 0.15
mixing_ndim     30

ks_solver       dav_subspace
basis_type      pw
pw_diag_ndim    2
pw_diag_nmax    20
pw_seed 1

kpar 1

symmetry        -1
noncolin        1
nspin           4
lspinorb        1

onsite_radius   3.0
sc_mag_switch 1
sc_thr 2e-10
nsc 50
nsc_min         2
alpha_trial     0.01
sccut           3.0

sc_scf_thr      10
sc_os_ndim      5

scf_os_stop 0
scf_os_thr      -0.01
scf_os_ndim     20

cal_stress 0
cal_force 0

out_chg 0

init_chg atomic

mixing_beta_mag 0.8

ATOMIC_SPECIES
Fe 55.845 Fe.upf
Rh 102.9055 Rh.upf
NUMERICAL_ORBITAL
Fe_gga_8au_200.0Ry_4s2p2d1f.orb
Rh_gga_8au_300.0Ry_4s2p2d1f.orb
LATTICE_CONSTANT
1.889726125458
LATTICE_VECTORS
3.06 0 0
0 3.06 0
0 0 3.06
ATOMIC_POSITIONS
Direct
Fe
0.0
1
0 0 0 1 1 1 mag 0 0 3 sc 1 1 1
Rh
0.0
1
0.517 0.491 0.506 1 1 1 mag 0.5 0 0 sc 1 1 1

K_POINTS
0
Gamma
4 4 4 0 0 0

b9b7671a2c667a6f87a7228d30fcc205e339b69f5568fe3f3dba567fdce59a19  case/Fe.upf
8ff6a789b1c7636a1d8ce857edb2477aa094aa36c25cdad11fce646b852c484e  case/Fe_gga_8au_200.0Ry_4s2p2d1f.orb
11803a8b7994745ce0b4b1b2f2be31d17a7d26fbdbd8097e0269ffb9f1aac532  case/INPUT
1ad281c0f88e2f4269b222db91fc59e3365188569a96e8629b57c8f65b0613e7  case/KPT
871688eff20d3f3555f1f7372f302d831582db9d87860979a068daf5ccfc111e  case/Rh.upf
c6912dc929a6613cd927e1b205b3bb90de6e48686eb16aef5cf3ec9fd16828c3  case/Rh_gga_8au_300.0Ry_4s2p2d1f.orb
4b890534c1943cf2bcb61a8b0702a86473f75fe90cb35132a096a4f8c3390b64  case/STRU