DFT+U refactor, step 6

[mohanchen]

DFT+U Slimming: PW occ-path free-function refactor + LCAO occ signature collapse

Overview

This PR continues the DFT+U refactoring by moving the remaining
occupation-matrix member functions of Plus_U_Base/Plus_U out into
namespace free functions, and by collapsing the over-long
cal_occ_mat_k/gamma parameter lists. It further shrinks the class
interface and reduces header coupling.

It consists of four independently compilable commits, each verified by a
full build before committing.

Per-commit breakdown

1. c90489cf2 — slim dftu_lcao.h and turn pot_uterm_HR into free functions

Header slimming (forward declarations)

• dftu_lcao.h: drop 7 transitive
  includes (klist/unitcell/parallel_orbitals/orb_read/hamilt/hcontainer/ density_matrix), replaced with forward declarations for UnitCell,
  Parallel_Orbitals, LCAO_Orbitals, hamilt::HContainer<T> and
  elecstate::DensityMatrix<TK,TR>.
• Add direct includes where the removed transitives were actually used:
  dftu_lcao.cpp, dftu_lcao_pots.cpp, dftu_lcao_energy.cpp,
  force_stress_lcao.h, write_hs_r.h, spar_u.cpp.

Free-function refactor

• Move Plus_U::cal_eff_pot_mat_R_double / cal_eff_pot_mat_R_complex_double
  out of the class into namespace DFTU_LCAO, renamed to pot_uterm_HR_real /
  pot_uterm_HR_complex to match the existing pot_uterm_real/complex family.
• spar_u.cpp now calls DFTU_LCAO::pot_uterm_HR_real/complex(dftu, ...).

Note: this commit merges the originally planned step 1 (free functions) and
step 2 (header slimming). The two are tightly coupled in dftu_lcao.h
(removing the member declarations + forward-declaration change in the same
file, and spar_u.cpp relies on the header's Parallel_Orbitals forward
declaration), so splitting them would have produced a non-compiling
intermediate state. Merging keeps the history clean and buildable.

2. 8f2d19a83 — sink occupation mixing to source_estate and drop forwarding methods

• Add free function elecstate::mix_occ_with_save() in
  occ_matrix.h /
  occ_matrix.cpp:
  occ = beta*occ + (1-beta)*occ_save applied to each atom's correlated
  orbital, nspin-aware (nspin=4 mixes the single Pauli block, nspin=1/2 mixes
  both spin channels).
• LCAO cal_occ_mat_k/gamma: replace the two duplicated hand-written mixing
  loops with a single mix_occ_with_save() call (~60 lines removed).
• Plus_U_Base: delete the pure-forwarding methods copy_occ_mat /
  zero_occ_mat / set_occ_mat and inline their call sites
  (init_base, cal_occ_pw) to call occmat_.copy_to_save / write_save_to_flat / zero / read_from_flat directly.

3. 816f3020d — turn Plus_U_Base PW occ path into DFTU_BASE free functions

Move the remaining PW occupation-matrix members out of Plus_U_Base into
namespace DFTU_BASE (declared in
dftu_base_tools.h, defined in
dftu_base_occ.cpp):

• reduce_occ_mat(cell, nspin, kpar, orbital_corr, occmat): the k-pool reduce
  still goes through Parallel_Reduce::reduce_double_allpool (MPI wrapper
  unchanged).
• compute_eff_pot_and_energy → compute_pot_uterm_and_energy(...), renamed to
  match the pot_uterm_* family; writes energy_u by reference.
• accumulate_occ_one_k<Device>: templated free function, with explicit
  instantiations for DEVICE_CPU (and DEVICE_GPU) moved to the .cpp.
• sync_occ_to_uom removed; cal_occ_pw now calls occmat_.write_to_flat
  directly.

Plus_U_Base::cal_occ_pw remains as the thin orchestration wrapper.

4. 8045f3691 — collapse cal_occ_mat_k/gamma signatures onto Plus_U&

• cal_occ_mat_k: 13 params → 9 (add Plus_U& dftu, drop
  nspin/npol/nlocal/ks_solver/iatlnmipol2iwt/orbital_corr/occ_mat/ occ_mat_save/occ_mat_initialized).
• cal_occ_mat_gamma: 13 params → 7 (add Plus_U& dftu).
• Both now read all occupation-matrix state from dftu.occmat()
  (mat/mat_save/data/data_save/iatlnmipol2iwt/nspin/npol) and the
  Plus_U_Base accessors (get_orbital_corr_vec / is_occ_mat_initialized / mark_occ_mat_initialized), and reuse occmat().copy_to_save()/zero()
  instead of the duplicated hand-written copy/zero loops.
• The intermediate cal_occ_mat template forwarders shrink to a single direct
  call.

Note: the MPI_Allreduce(MPI_COMM_WORLD) inside the accumulation loop is
intentionally left unchanged — it has no exact Parallel_Reduce counterpart
and was kept per the "only replace what maps cleanly" rule.

Verification

Each commit was verified before committing with:

This PR reorganizes the DFT+U force and stress implementation in the LCAO module, with two main goals: splitting the monolithic force/stress file into focused, well-documented units, and renaming all LCAO-specific files to explicitly indicate the NAO (numerical atomic orbital) basis set.

First, the original dftu_fs.cpp, which contained the entire force and stress calculation in a single file, is split into three new files. dftu_nao_fs_r.h/cpp provides the unified entry point cal_fs_nao_r(), handling OpenMP parallelization over atom pairs and the final MPI reduction. dftu_nao_for_r.h/cpp contains the single atom-pair force kernel cal_for_IJR_nao_r(). dftu_nao_str_r.h/cpp contains the single atom-pair stress kernel cal_str_IJR_nao_r(). Each file header includes the full mathematical formulas for the Hubbard correction energy, forces, and stress tensor, and explicitly documents that the implementation operates in real space using the real-space density matrix DMR. Because DMR already contains the Brillouin-zone integration, this implementation is unified for both gamma-only and multi-k-point calculations, eliminating the need for separate gamma and k-space branches.

[lanshuyue]

There seems to be an existing out-of-bounds write for nspin == 1 in write_save_to_flat().

For nspin == 1, init_base() allocates only one occupation-matrix block in uom_save: pot_index is doubled only when nspin == 2. However, write_save_to_flat() handles nspin == 1 and nspin == 2 together and unconditionally writes both spin channels:

const int half_size = uom_save.size() / 2;

uom_save[index[iat] + mm] = ...[0]...;
uom_save[half_size + index[iat] + mm] = ...[1]...;

For example, for a single correlated d orbital, uom_save.size() == 25 when nspin == 1. The first matrix is written to [0, 24], while the second write goes to [12, 36], so part of the first block is overwritten and [25, 36] is out of bounds. Since operator[] is used, this results in undefined behavior rather than a guaranteed crash.

This does not appear to be newly introduced by this PR: the legacy Plus_U_Base::copy_occ_mat() contains the same logic. However, now that the flat serialization has been extracted into OccupationMatrix, I think it would be good to fix the existing inconsistency here.

In particular, write_to_flat() and read_from_flat() already treat the layout differently: they access the second flat block only for nspin == 2. Should write_save_to_flat() follow the same logic, i.e. always write spin 0 for nspin == 1/2, but write spin 1 only when nspin == 2?

[lanshuyue]

This is probably not a blocker, but I think it would be worth fixing while touching this code.

Before this refactor, init_base() always updated the Yukawa state explicitly:

use_yukawa_ = yukawa_potential;

Therefore, after every call to init_base(), the internal state was guaranteed to match the argument, including a possible true -> false re-initialization.

After this refactor, the state is inferred from the lifetime of the object:

bool use_yukawa() const { return yukawa_ != nullptr; }

which is a cleaner design in itself. However, init_base() currently only updates yukawa_ when yukawa_potential is true:

if (yukawa_potential)
{
    yukawa_.reset(new YukawaScreening());
    yukawa_->init(...);
}

When yukawa_potential == false, an existing yukawa_ is not cleared. Therefore the following sequence behaves differently from the old implementation:

init_base(..., true, ...);
init_base(..., false, ...);

Previously, use_yukawa() would become false; after this refactor it would remain true because the old YukawaScreening object is still alive.

I checked the current call path as well. For the PW solver, init_base() is indeed called repeatedly on the same Plus_U object through before_scf() -> setup_pot(). In the normal workflow, however, yukawa_potential comes from the same input configuration, so it appears to stay constant across these repeated calls. Therefore I don't currently see a normal production path that actually performs a true -> false transition, and this may not cause an observable problem today.

Still, this is a small behavioral regression relative to the previous implementation, and init_base() is not strictly a one-shot initialization function. Since the fix is essentially free, I think it would be safer to preserve the old semantics by clearing the pointer in the disabled case, e.g.

if (yukawa_potential)
{
    yukawa_ = std::make_unique<YukawaScreening>();
    yukawa_->init(...);
}
else
{
    yukawa_.reset();
}

or simply reset it before the conditional.

So this is more of a low-cost robustness fix than a confirmed runtime bug in the current workflow, but I think it is worth fixing opportunistically in this refactor.