feat: r-space parallelization for nonorthogonal lattices (nproc_rgrid>1)

[MZKC]

Summary

Enables real-space domain decomposition (nproc_rgrid > 1) for nonorthogonal lattices, which was previously disabled by a guard in initialization.f90. The changes are additive and limited to the nonorthogonal r-space-parallel path; orthogonal stencils and serial (nproc_rgrid = 1) nonorthogonal runs take the existing code path unchanged.

Background

The nonorthogonal kinetic operator (zstencil_nonorthogonal) forms the off-diagonal-metric cross derivatives by differentiating intermediate per-axis first-derivative fields (df/dx, df/dy) held in a work array. The serial path obtains the neighbors of that work array through the periodic index wrap; under r-space decomposition the work array's halo is not exchanged, which is why the case was guarded off. Access is in the face directions only, so the existing face halo is sufficient.

Changes

• stencil.f90 — zstencil_nonorth_rspace_pass1/2: a two-pass split of zstencil_nonorthogonal (algebraically equivalent) that exposes df/dx, df/dy for a halo exchange between the passes.
• hamiltonian.f90 — hpsi_nonorth_rspace: pass1 → update_overlap_complex8(df/dx, df/dy) → pass2, selected when .not. stencil%if_orthogonal .and. info%if_divide_rspace. Reuses the wavefunction srg.
• preconditioning.f90 + conjugate_gradient.f90 — zstencil_nonorthogonal_preconditioning_diag: a diagonal-metric CG preconditioner for the r-space-parallel case (the cross term would require the same halo). This changes only the preconditioner, not the converged state; for strongly nonorthogonal cells a larger ncg may be needed for convergence.
• initialization.f90 — round-robin atom partition over the r-space ranks. The Ewald ion-ion sums loop over info%ia_mg and reduce over icomm_r; with all atoms replicated on each rank the real-space Ewald was counted nproc_rgrid times. The partition assigns each atom to a single rank so the reduction yields the full sum. The guards are removed.

Validation

Primitive Si (nonorthogonal FCC, nr=16, Γ), nproc_rgrid=(1,1,1) vs (1,1,2), ncg=16:

quantity	difference (r-space-parallel vs serial)
total energy	1.8e-8 eV
stress (P_total)	2e-7 GPa
Ewald terms (E_ewa_R, P_ewa, P_ewa_R)	none
eigenvalues (8 states), forces	none

Built with the Fujitsu compiler (A64FX).

Notes

The branch is based on feature/stress-tensor; the changes are independent of the stress / meta-GGA work and can be rebased onto develop-2.0.0. The new routines are CPU/OpenMP only (no GPU/OpenACC path).

🤖 Generated with Claude Code

[jenkins-diana]

Can one of the admins verify this patch?

[MZKC]

Re-validation on develop-2.0.0 (after rebase off feature/stress-tensor)

This branch was rebased from feature/stress-tensor onto the latest develop-2.0.0. The three commits are unchanged in content (added/removed lines are byte-identical; only diff context shifted), and the head no longer carries the stress-tensor history. The serial-vs-r-space-parallel consistency check was re-run on the new base.

Build: compiles cleanly on develop-2.0.0 with the Fujitsu A64FX compiler (fujitsu-a64fx-ea), no errors.

Case: primitive-FCC (nonorthogonal) Si, Γ-only, num_rgrid = 16,16,16, xc = 'PZ', ncg = 16, nstate = 8. Two runs differing only in nproc_rgrid (1,1,1 vs 1,1,2) — i.e. the serial nonorthogonal kinetic operator vs the new r-space-parallel two-pass operator.

quantity	serial 1,1,1	r-space-parallel 1,1,2	difference
total energy (eV)	−199.2029957187	−199.2029957195	8.1e-10
eigenvalues, all 8 states (Ha)			≤ 6e-10
1-particle energies (eV, 4 dp)	identical	identical	0
force on each atom (au)	~2.62e-5	~2.62e-5	~4e-10
SCF iterations to threshold = 1e-9	32	35	—

Both occupied and empty states agree to ~1e-10 Ha; total energy and forces reproduce the serial result to sub-1e-9. This is consistent with the original Validation table and confirms the operator behaves identically on develop-2.0.0.

Scope note: the original table also compared stress (P_total), which depends on yn_out_stress / the stress-tensor work that is not present on develop-2.0.0 (it lives in feature/stress-tensor); stress is therefore out of scope for this base. The same applies to r2scan / yn_tau_nlcc, so this re-validation uses the LDA (PZ) path. The nonorthogonal r-space-parallel kinetic operator is XC-independent, so the LDA check exercises exactly the code this PR adds.

[MZKC]

Checked the two-pass kinetic operator against the serial zstencil_nonorthogonal and it is equivalent: pass1 builds the diagonal-metric Laplacian + lap0 + Bk*nabla and stores df/dx, df/dy; pass2 adds the off-diagonal cross terms as d/dz(df/dy), d/dz(df/dx), d/dy(df/dx), so only the FACE halos of wrk1/wrk2 are needed (no edge/corner ghosts), and the F(4/5/6) pairing, signs and 0.5 factor all match. The diagonal-metric preconditioner only changes CG convergence, not the converged state, and the round-robin atom split is a clean per-atom partition for the icomm_r-reduced Ewald sums. One performance/memory concern, inline.

[MZKC]

allocate(wrk1(mg%is_array(1):mg%ie_array(1),mg%is_array(2):mg%ie_array(2),mg%is_array(3):mg%ie_array(3) &
             ,nspin,io_s:io_e,ik_s:ik_e,im_s:im_e))

This allocates two full complex orbital-block work arrays (~2x the wavefunction-block size) on every hpsi call, and hpsi is in the CG / time-propagation hot path. The repeated allocate/deallocate adds allocator churn and a sizable transient memory spike per call. Consider persistent module/scratch workspace (allocated once, sized to the orbital block), or processing one (ispin,io,ik) tile at a time so only a single 3D df/dx, df/dy pair is live — the latter also shrinks the halo-exchange working set. Functionally correct as-is; this is about hot-path cost on large runs.

[MZKC]

Filed as #1269 - the per-call allocation of the work arrays is a performance follow-up, not blocking this PR; the two-pass kinetic math is correct (verified equivalent to the serial zstencil_nonorthogonal).