Skip to content

Feature: LR-TDDFT analytical forces (gradients) by Z-vector method and excited-state geometry relaxation - #8069

Open
maki49 wants to merge 74 commits into
deepmodeling:developfrom
maki49:LR-Grad-relax
Open

maki49 wants to merge 74 commits into
deepmodeling:developfrom
maki49:LR-Grad-relax

Conversation

@maki49

@maki49 maki49 commented Oct 2, 2026 •

Copy link
Copy Markdown
Collaborator

Features

  1. Analytic LR-TDDFT excited-state forces (Z-vector method, gamma-only)
  • Implements the Lagrangian/Z-vector equations to solve for the relaxed and energy-weighted transition density matrices, giving analytic excited-state forces matching well with finite differences (see results.md).
  • Supports RPA, LDA, GGA, Hartree-Fock, hybrid functionals (HSE/range-separated, including a sigma-clamping fix for the wpbeh kernel's short-range divergence) and the long-range-corrected (LRC) kernel.
  • Supports both spin-restricted (singlet/triplet) and spin-unrestricted (open-shell) excited states.
  • For a degenerate excited multiplet, computes and prints the full (non-diagonal) gradient matrix over the degenerate subspace, with its eigen-decomposition, instead of a single force vector.
  1. Excited-state geometry relaxation
  • Drives the existing Relax_Driver on an excited-state PES: the ground state is re-converged and the Z-vector gradient is recomputed at every ionic step.
  • Tracks the target excited state across ionic steps by maximum wavefunction overlap, so relaxation follows the same diabatic state through near-degeneracies/state crossings.
  • For degenerate states, adds an isotropic (symmetry-averaged) relax mode and a Jahn-Teller mode that follows the distortion direction extracted from the degenerate gradient matrix (lr_relax_degen_mode).

Refactors

  • Generalizes OperatorLRHxc/OperatorLREXX (previously built only for the Casida A-matrix) so the same Hxc/EXX operators are reused to build the Z-vector/multiplier equations.
  • Keeps the ground-state KS solver as a persistent member of ESolver_LR (rather than a one-shot temporary) so it can be re-run every ionic step during relaxation.

Copilot AI balanced review requested due to automatic review settings October 2, 2026 10:06

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

@maki49
maki49 marked this pull request as draft October 2, 2026 10:11
Comment thread source/source_lcao/module_lr/cal_edm.cpp
Comment thread source/source_lcao/module_lr/lr_force.cpp
@mohanchen mohanchen added Features Needed The features are indeed needed, and developers should have sophisticated knowledge EXX and lr-TDDFT Related to EXX or lr-TDDFT labels Oct 2, 2026
@maki49
maki49 force-pushed the LR-Grad-relax branch 6 times, most recently from e5a8e55 to cd8346b Compare October 5, 2026 07:56
@maki49
maki49 marked this pull request as ready for review October 6, 2026 13:40
@mohanchen

Copy link
Copy Markdown
Collaborator

Thanks a lot for your hard work on this PR — the overall direction looks great! Below are a few suggestions, ordered from the easiest to the more involved ones. Most should be quick to address, and I'm happy to discuss any of them.

  1. File naming: please use lowercase for source file names. For consistency with the rest of the codebase, new source files should use all-lowercase names.
  2. INPUT parameter names must not exceed 15 characters. The name lr_relax_degen_mode violates this rule — please shorten it. This limit exists so that INPUT files stay nicely aligned and readable.
  3. Keep new files under 500 lines. A single file shouldn't exceed 500 lines. I know some existing code is over this limit too, but that's exactly why we need everyone's help to gradually break things down — at the very least, new code should not add to the problem. Splitting it up will also make this PR much easier to review and maintain.
  4. No mutable members. Please avoid mutable, such as the members dmr and dmr_save declared mutable so that the const function cal_dmr can modify them. A const function that mutates state is a sign the design needs adjusting — you're welcome to refactor the density matrix class itself instead, for example by making cal_dmr non-const, or by restructuring where these buffers live.
  5. Don't obtain data indirectly through the density matrix. I've already removed Parallel_Orbitals and get_kvec_d() from the density matrix. If you need these quantities, please pass them in as explicit parameters rather than borrowing them through the density matrix — that creates unnecessary data dependencies. Convenience is nice, but dependency hygiene matters more.
  6. ESolver functions should stay high-level. There's a strict convention for ESolver functions: they should only contain the top-level flow — the skeleton of the algorithm. All detailed work should go into free functions (not class member functions); if that's truly impractical, a class member function is acceptable as a fallback. For example, something at the level of print_force is a detail, and its definition should not live inside the ESolver. Think of the ESolver as an outline that directs the flow — it shouldn't manage any details itself.
  7. Merge with the existing LR ESolver rather than adding a new one. LR already has an ESolver — why introduce another one? If the new functionality can't be merged into the existing LR ESolver, that usually indicates a framework design issue. And when problems aren't fixed at the root, they tend to grow bigger and harder to clean up over time. Conceptually, I'd strongly recommend merging this into the existing LR ESolver as much as possible. This is the biggest item on the list, so feel free to reach out if you'd like to discuss the design together!
    None of these are blockers in spirit — items 1–3 are quick fixes, and for items 4–7 I'm glad to help brainstorm the refactoring approach. Looking forward to the next iteration!

@maki49

maki49 commented Oct 7, 2026 •

Copy link
Copy Markdown
Collaborator Author

@mohanchen For the last point (7.), no new ESolver is added: esolver_lr_grad.cpp is just a separate .cpp file from esolver_lr_lcao_tddft.cpp, corresponding to esolver_lr_lcao_tddft.h.
All of your other suggestions have been followed by the last several refactor commits.

maki49 and others added 15 commits October 7, 2026 00:19
Port the LR-Grad series (44 commits, 7c6398d..6aea1fc, all by maki49)
onto current develop as a single squashed change, adapted to develop's
large-scale refactor.

Path/layout migration:
  module_lr        -> source_lcao/module_lr
  module_ri        -> source_lcao/module_ri (Exx_LRI.h -> exx_lri.h)
  esolver_lrtd_lcao -> source_esolver/esolver_lr_lcao_tddft
  ESolver_LR now lives in namespace ModuleESolver, not LR.

Interface adaptation:
  - Grid integration: the Gint_Gamma/Gint_k objects and TGint<T> selector
    are gone. All call sites now use the free functions in ModuleGint
    (cal_gint_rho / cal_gint_vl / cal_gint_fvl); the gint pointer has been
    removed from LR_Force, LR_Density, OperatorLRHxc, HamiltLR/HamiltULR,
    the Z-equation Hamiltonians and the multiplier helpers.
  - PulayForceStress::cal_pulay_fs (grid overload) keeps the explicit
    `nspin` argument the gradient code needs; develop's own callers pass
    PARAM.inp.nspin.
  - Parallel_Orbitals::get_row_size(iat)/get_col_size(iat)
    -> get_nrow_atom(iat)/get_ncol_atom(iat).
  - elecstate::Potential::init_pot(istep, chg) -> init_pot(chg);
    get_effective_v -> get_eff_v; Structure_Factor::setup_structure_factor
    -> setup; Charge::allocate now takes kin_den.
  - ModuleIO::output_single_R takes a SparseWriteOptions struct.
  - ModuleIO::read_wfc_nao uses develop's signature (ik2iktot/nkstot/nspin).
  - ModuleBase::timer::tick -> start/end pairs.
  - LR_Util::gather_2d_to_full keeps develop's explicit-dimension form;
    scatter_full_to_2d is added for the Z-equation solver.
  - get_DMR_real_imag_part / set_HR_real_imag_part are now the templated,
    nat-free versions; module_bse call sites updated, the now-dead
    lr_util_hcontainer.cpp is removed.
  - DensityMatrix::cal_DMR is const (_DMR is mutable).

Known gap: LR_Force::cal_force_exx_gs_dm_relaxed_diff still uses
RI::LR from LibRI, which no longer derives from RI::Exx in the current
LibRI master (commit 5c6c262 refactored it into a k-space CVCX/BSE helper).
That single function does not compile; everything else does.
…X force kernel

LibRI commit 5c6c262 repurposed RI::LR: it no longer derives from RI::Exx and
lost set_Ds/cal_Hs/cal_force/force. Reintroduce just the piece the LR-TDDFT
gradient needs as LR::ExxForceTwoDM, a thin RI::Exx subclass whose cal_force
overload overwrites post_2D.saves["Ds_<suffix>"] with the left density matrix
before delegating to Exx::cal_force. This reproduces the old RI::LR semantics
(D_IJ from Ds_left, D_KL from set_Ds) without pinning LibRI.

Also drop dm_trans/dmr_complex.cpp from the build: develop moved the
DensityMatrix<complex,complex>::cal_DMR specialization into density_matrix.cpp,
so keeping both gave a duplicate definition at link time.
…ose HSolverLR timer

Two crashes found running the H2/TDHF/SZ gradient case:

1. PulayForceStress::cal_pulay_fs (the LR overload in pulay_force_hcontainer.h)
   passed the density matrix's spin count as `nspin` to ModuleGint::cal_gint_fvl
   while supplying a single potential channel. Gint_fvl::cal_fvl_svl_ indexes
   vr_eff_[is] for is<nspin, so with a two-channel DM (the test_force path feeds
   in the ground-state DM) it read past the vector and segfaulted inside the
   OpenMP region. The LR kernel produces one potential channel, so the grid
   integrals now run with nspin=1 on the leading DM channel -- exactly what the
   old `Gint_inout(is=0, ...)` calls did.

2. HSolverLR::solve started a timer that was never ended, so the second
   (triplet) call threw "timer::start HSolverLR::solve". Added the matching
   ModuleBase::timer::end on the common exit path.

Verified on H2 / TDHF / SZ (nspin=2, gamma-only):
  - SCF reproduces the pre-port reference bit-for-bit (E_KohnSham -2.3110142036 Ry,
    TOTAL-FORCE -14.7806682608 eV/Ang).
  - test_force cross-check: GS Hxc force via potential -28.4711400896 vs via the
    LR kernel -28.4711414920 eV/Ang (7 significant digits).
  - The reproduced EXX force (31.3817781891) is exactly 2x the SCF EXX force
    (15.6908890935), confirming ExxForceTwoDM replaces RI::LR correctly.
  - Excitations 22.484 eV (singlet) / 13.813 eV (triplet); excited-state
    gradients -0.6313 / -2.01497 eV/Ang, antisymmetric with ~1e-14 transverse
    components.
solve_Z_lapack built a replicated Hessian via HamiltLR::matrix() but handed
LAPACK the pX-distributed right-hand side `R` directly. lapack_linear_solver
then copies n_global*nstates entries out of a buffer that only holds
ld*nstates = nk*pX[0].get_local_size()*nstates -- an out-of-bounds read as soon
as the local size is smaller than the global one, and a segfault on any rank
whose local size is 0.

Gather R into a global buffer first, mirroring the scatter of Z_full that
already follows the solve. Serial results are bit-for-bit unchanged (verified
on H2/TDHF/SZ: excited-state gradients -0.631306 / -2.01497 eV/Ang, all
per-term forces identical to the pre-fix run).
`PotHxcLR::cal_v_eff` dominates an LR run: 74% of the total for 07_LiF/pbe
(949 s of 1289 s over 541 calls). Roughly half of that was FFT; the rest was
serial element-wise work and allocator traffic. Five changes, no algorithmic
change:

1. OpenMP on the element-wise loops. Everything these loops call
   (`grad_rho`, `grad_dot`, the Hartree kernel) was already parallel; the
   integrands here were the only part still running on one thread.

2. Raw pointers instead of `.at()` (43 sites). The bounds check is a branch
   per access, ~10 per grid point, and it blocks vectorization. The
   outer-vector lookups (`drho_gs.at(0)`) are hoisted out of the loops, where
   they keep their bounds check for free.

3. A scratch pool instead of per-call temporaries. `drho`, `gdot_terms` and
   `vxc_tmp` were allocated and value-initialized on every call and then
   overwritten before being read -- ~450 MB of pointless memset per call at a
   200^3 grid, plus the page faults. `grad_dot` assigns rather than
   accumulates, so even `vxc_tmp`'s zero-fill was dead. The pool is shared by
   all instances (a run holds three `PotHxcLR`) so it does not multiply.

4. `add_v_hartree` replaces `H_Hartree_pw::v_hartree`, which (a) re-did the
   forward FFT of rho^X that the GGA branch needs anyway -- 1 of the 10 FFTs
   per call was pure duplication, (b) reduced a Hartree "energy" of the
   transition density through `Parallel_Reduce::reduce_pool` every call, a
   collective nobody reads that also clobbers the global
   `H_Hartree_pw::hartree_energy`, and (c) returned a `matrix` by value, to
   which `v_eff += 2 * (...)` added a second full-size temporary.

5. The nspin=2 singlet/triplet combinations `v2rho2_uu -+ v2rho2_ud` and
   `2*vsigma_uu -+ vsigma_ud` are pre-contracted once per potential instead
   of being rebuilt at every grid point of every call, which also replaces
   two strided reads with one contiguous one. Costs 8-16 B/point, and only
   for closed-shell nspin=2.

`PotGradXCLR::cal_v_eff` (the g^xc branch feeding the Z-vector RHS) gets 1-3
of the same treatment, for -20% to -28%. The shared pool matters more there:
a `PotGradXCLR` is constructed inside the loop over excited states, so
per-object buffers would never be reused at all.

Measured (single node, 16 threads), analytic gradients unchanged:

  01_Si/lda    93 s -> 51 s   (-45%)   cal_v_eff  66 ->   45 ms/call
  01_Si/pbe   229 s -> 125 s  (-46%)   cal_v_eff 356 ->  268 ms/call
  02_Li2/lda  165 s -> 86 s   (-48%)   cal_v_eff 559 ->  278 ms/call
  02_Li2/pbe  530 s -> 253 s  (-52%)   cal_v_eff 2651 -> 1553 ms/call

(The total also benefits from solving the Z-vector equation once instead of
three times, a separate fix; the per-call figures above isolate this commit.)

Excitation energies are identical to every printed digit; the nspin=1 cases
agree to 1e-14, and of 60 force components per nspin=2 case, 2-3 differ by
exactly one unit in the last printed digit (the Hartree prefactor is now
associated differently).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QMV7xFZ69hhc3HnVBHbD96
maki49 and others added 29 commits October 7, 2026 00:19
This PR's own global-dependency diff vs develop was net_delta=121 (136
added, 15 removed) -- the governance checker's "Global dependency
budget" rule unconditionally blocks on any positive delta, with no
automatic honoring of its own "exception allowed" tag, so the only way
to turn the CI job green is to actually reduce PARAM/GlobalV/GlobalC
touches in the diff, not just document them.

Threaded explicit dependencies through instead of reading PARAM/GlobalV
directly, in order of where the duplication was worst:

- ESolver_LR (esolver_lr_lcao_tddft.h/.cpp, esolver_lr_grad.cpp): added
  cached ofs_running_/ofs_warning_/my_rank_ members (bound once from
  the corresponding globals) and switched ~50 PARAM.inp.* reads to the
  already-existing `this->inp_->*` (the class already follows this
  convention, see esolver.h's `inp_` comment) or the new members.
- LR_Force (lr_force.h/.cpp, lr_force_test.cpp, force_funcs_lcao.h):
  same pattern -- cached ofs_running_/nspin_/test_force_/vl_in_h_/
  vh_in_h_ members; threaded nspin/test_force/ofs_running as explicit
  params into ForcePWTerms::operator() instead of it reading PARAM.
- cal_edm_from_multipliers.h: threaded `test_force` as an explicit
  param through cal_edm_terms_from_XZWK / cal_edm_from_XZ_istate(
  _openshell) instead of an inline PARAM.inp.test_force read.
- hamilt_zeq_left.h/right.h, zeq_solver.hpp: threaded `in_dir`/
  `out_dir`/`ks_solver` as explicit params through Z_vector_equation ->
  Z_vector_L/R/UR instead of each reading PARAM.globalv directly in
  its HamiltLR base-class init list.
- operator_gxc_ulr.h: added an explicit `nspin`/`ks_solver` constructor
  parameter instead of reading PARAM.inp.nspin/ks_solver internally.
- pot_grad_xc.cpp / pot_hxc_lrtd.cpp: both independently duplicated the
  same `nspin==1 || (nspin==4 && !domag && !domag_z)` expression;
  consolidated into one shared `LR_Util::kernel_nspin()` (lr_util_xc.hpp)
  so there is a single PARAM touch instead of two.
- lr_density.hpp: cached out_chg[1]/global_out_dir as members.
- lr_util_hcontainer.h: threaded `out_dir`/`nlocal` as explicit params
  through save_sparse -> save_HR -> save_DMR instead of each reading
  PARAM.globalv directly.

Left as documented exceptions (reading PARAM/GlobalV is the
least-disruptive option here):
- esolver_lr_grad.cpp's `PARAM.globalv.has_float_data` read mirrors the
  identical, pre-existing call site in esolver_fp.cpp.
- operator_lr_exx.h's `gs_is_hybrid()` and the `PARAM.inp.cal_force`
  check in `OperatorLREXX`'s constructor: both are called from 8-9
  sites spanning the ground-state Casida Hamiltonian (hamilt_ulr.hpp,
  hamilt_casida.h) as well as the Grad operators; widening either
  signature would touch core files well outside this PR's scope.
- lr_util_hcontainer.h's `GlobalV::DRANK` single-writer-rank check
  matches the universal, unthreaded ABACUS convention for gating file
  writes to one rank.
- exx_lri.hpp's one new `PARAM.inp.nspin`-keyed SPIN_multiple lookup
  duplicates an idiom already repeated several times, unchanged, in
  the same core (non-Grad) file.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…r equation

`solve_Z_lapack` gathers the whole orbital Hessian onto every rank and has
each of them factorize the same n x n system. Add two solvers that keep both
the matrix and the factorization distributed:

* `solve_Z_scalapack`: LU with p?gesv, the same system as `solve_Z_lapack`.
* `solve_Z_elpa`: ELPA Cholesky of (A+A^H)/2, then p?potrs. A+B is positive
  definite at a stable ground state, so this needs about half the LU flops.

Both go through `solve_Z_2d`, which builds the Hessian column by column from
`hPsi` -- the same operator CG applies -- straight into a 2D block-cyclic
layout on the BLACS grid of pX, so only one column is replicated per rank.
The pX <-> global gather/scatter is factored into `zvec_local_to_full` /
`zvec_full_to_local`, now also used by `solve_Z_lapack`.

ELPA 2022.11 caveats, documented at `elpa_linear_solver`:
* its C binding declares `error` intent(in), so a failed Cholesky still
  returns ELPA_OK; the failure is detected from U's pivots instead;
* on a non-positive-definite matrix only the rank owning the failing block
  returns, so a multi-rank run hangs rather than erroring.

Adds pdgesv and p?potrs to ScalapackConnector and Cigsum2d to the BLACS
declarations.

Tests: MODULE_LR_GRAD_zeq_linear_solver checks both kernels against
LR_Util::lapack_linear_solver (real/complex, non-symmetric LU input, ELPA's
Hermitian part, ELPA rejecting an indefinite matrix on one rank); passes on
1, 2 and 4 MPI ranks.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…tor solver

The Z-vector solver was hard-wired to the `zvec_solver = "cg"` default of
`Z_vector_equation`. Add `lr_grad_solver` -- the linear-equation counterpart of
`lr_solver` -- with the values cg (default, unchanged behaviour), lapack,
scalapack and elpa, and pass it explicitly from `solve_zvector_eqation`.
The default arguments of `Z_vector_equation` are dropped; its one call site
already passed `openshell`.

check_value rejects unknown values, scalapack/elpa without MPI, and elpa
without ELPA, at input time instead of deep inside the force calculation.

Docs: the generated lr_grad_solver entry is added to parameters.yaml and
input-main.md.

Tests: read_input_item_test covers accepted/rejected values for the MPI and
ELPA configurations; read_input_ptest checks the default.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ons, regenerate

docs/parameters.yaml and input-main.md had drifted from the Input_Item
definitions, so the CI sync check could not pass:
* parameters.yaml lacked lr_target_state and lr_target_spin entirely;
* 6741e24 documented the relax-target scoping (ignored outside
  calculation = relax; an open-shell run accepts any lr_target_spin and only
  reports `triplet` as ignored; a closed-shell run rejects `updown`) by
  editing input-main.md alone, leaving the C++ descriptions on the older text;
* lr_relax_degen_mode's md text said "the reported energy" where the source
  named the internal `cal_energy`.

That text now lives in read_inp_tddft.cpp (cross-references as backticks,
as everywhere else in the Input_Item descriptions, since `abacus -h` prints
them verbatim), and both files are regenerated from the source:
`abacus --generate-parameters-yaml` output and
`python docs/generate_input_main.py docs/parameters.yaml`.

The one md-only statement not carried over is xc_kernel's "exx_erfc_omega
($\omega$)": the operator in the same sentence is erfc($\mu$ r), so the
source's ($\mu$) is the consistent one.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…sky Z-vector solve

`scalapack_cholesky_linear_solver` Hermitizes the distributed orbital Hessian,
factorizes it with p?potrf and solves with p?potrs -- the ELPA path with
p?potrf in place of elpa_cholesky. The Hermitization is now the shared
`hermitize_2d`, used by both Cholesky solvers.

Unlike ELPA's, p?potrf reports a non-positive-definite matrix through a global
INFO, so it throws on every rank instead of hanging the non-owning ranks, and
needs no pivot check to work around ELPA 2022.11 returning ELPA_OK on failure.
It also needs no ELPA build.

Timing of the solve alone, 4 MPI ranks, single-threaded BLAS, one rhs:

      n   scalapack (LU)   elpa    scalapack_chol   (of which Hermitize)
   4800        0.56         0.52        0.54              0.13
   8000        2.40         2.04        1.95              0.30
  12000        8.07         6.42        6.08              0.83
  16000       17.93        13.96       14.23              1.52

so scalapack_chol is on par with elpa and 15-25% faster than the LU; its gain
is robustness, not speed. In the Z-vector solve all of them are dwarfed by
building the matrix (one A+B application per column).

Tests: MODULE_LR_GRAD_zeq_linear_solver gains real/complex, Hermitian-part
and indefinite-matrix cases for the new kernel; the indefinite case runs on
every rank and passes on 1, 2 and 4 MPI ranks. read_input_item_test covers
the new value. parameters.yaml and input-main.md regenerated.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
xc_gga_vxc's first-derivative vsigma used raw, unclamped sigma while
the fxc/kxc (2nd/3rd derivative) calls a few lines below already
clamped sigma to work around a known HSE06/wpbeh libxc closed-form
instability near symmetry-forced grad-rho->0 vacuum grid points. The
unclamped vsigma blew up to ~3.9e7 at such a point for 09_CH4/hse,
which propagated through carbon's diffuse 2nd-zeta DZP orbitals into
huge AO-basis V_Hxc elements and produced three ~-19404 Ry ghost
eigenvalues in the full Casida/TDA spectrum (lr_solver=lapack),
negative total oscillator strength, and force blow-ups in the
Z-vector/gradient code. Hoisted the sigma_clamped computation above
all three libxc calls so vxc/fxc/kxc share one clamped array.
…d_to_full

The Test and CUDA Test CI jobs failed at the build step with
BUILD_TESTING=ON:

- dm_trans_test, dm_diff_test and CVCX_test called
  LR_Util::gather_2d_to_full with 3 arguments, but develop's signature
  takes (pv, sub, full, row_major, global_nrow, global_ncol) with no
  defaults. dm_trans_test is restored to develop's version; the two new
  tests now pass the dimensions of each Parallel_2D explicitly.
- gather_2d_to_full scatters into fullmat and MPI_Allreduce-sums it, so
  the destination must start at zero. The new tests now zero X_full,
  c_full, V_full, dm_gather and AX_gather before each gather; without it
  the parallel cases compared against garbage (NaN) on more than one rank,
  and AX_gather accumulated the occ result into the virt one.
- test_grad_matrix_degenerate failed to link: base needs device
  (gemm_op, memory ops). Link container and device like its siblings.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…chars

mohanchen's PR review left two blocking comments: no new subdirectories
under source_lcao/module_lr beyond the ones that already existed (utils,
ao_to_mo_transformer, dm_trans, operator_casida, potentials, ri_benchmark),
and every filename capped at 15 characters (stem, excluding extension).

This is a pure reorganization, no functional changes:
- Grad/, Grad/CVCX/, Grad/degenerate/, Grad/dm_diff/, Grad/force/,
  Grad/multipliers/, Grad/xc/, and dm_band/ (all introduced by this
  branch) are gone. Their files land either in the pre-existing sibling
  directory they semantically match (CVCX -> ao_to_mo_transformer/,
  dm_diff -> dm_trans/, pot_grad_xc -> potentials/, operator_gxc_ulr ->
  operator_casida/) or flat in module_lr/ itself when no existing
  directory fits (the force/ and multipliers/ Z-vector/gradient code,
  grad_degen.{h,cpp}, dm_band.{h,cpp}).
- esolver_lr_grad.cpp moves to source_esolver/, next to its siblings
  esolver_lr_lcao_tddft.cpp/esolver_lr_lcao_bse.cpp; its own filename was
  already under the cap so only its location changed.
- A new module_lr/test/ holds the two genuine CTest unit tests that used
  to live under Grad/*/test/ (test_zeqlin.cpp, test_grad_degen.cpp).
  lr_force_test.cpp is NOT a unit test -- it's compiled straight into the
  production `lr` object library -- so it stays flat in module_lr/, not
  in test/.
- Renamed over-length files to fit the cap: cal_edm_from_multipliers.*
  (the file mohanchen named directly) -> cal_edm.*, plus
  zeq_linear_solver.* -> zeqlin_solv.*, hamilt_zeq_left/right.h ->
  hamilt_zeq_l/r.h, hamilt_zeq_ulr.h -> hamilt_zequlr.h,
  cal_multiplier_w_from_z.h -> cal_w_from_z.h, exx_force_two_dm.h ->
  exx_force_dm.h, force_funcs_lcao.h -> force_funcs.h,
  pulay_force_hcontainer.h -> pulay_hc.h, CVCX_parallel.cpp -> CVCX_par.cpp,
  grad_matrix_degenerate.* -> grad_degen.*.
- Updated every #include site referencing the old Grad/dm_band paths
  (CMake and Makefile.Objects builds both verified), module_lr/CMakeLists.txt,
  source_esolver/CMakeLists.txt, source/CMakeLists.txt (dropped the merged
  lr_grad target), and the test CMakeLists that gained a test file.

Verified: CMake build of abacus_std_para is clean with no warnings.
The legacy source/Makefile.Objects build path was checked but could not
be locally verified end-to-end on this machine: CXX=mpiicpc needs the
classic icpc, which this oneAPI 2026 install no longer ships (only
icpx); CXX=mpicxx's non-Intel branch needs standalone FFTW3 dev headers
and a standalone ScaLAPACK library, neither of which exists here
(FFTW3/ScaLAPACK are only available bundled inside MKL). The generated
compiler invocations for every moved/renamed object file were confirmed
to resolve to the correct new paths before compilation failed on the
missing toolchain pieces.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The Z-vector/Hxc diagnostic getenv switches (ABACUS_LR_L_SYM,
ABACUS_LR_ZSOLVER, ABACUS_LR_SKIP_HXC_L/EXX_L, ABACUS_LR_ONLY_OP1,
ABACUS_LR_SKIP_GXC_L, ABACUS_LR_PRINT_VHXC_AO/HXC_MO) were already
removed from hamilt_zeq_l.h/r.h, zeq_solver.hpp, and
operator_casida/operator_lr_hxc.cpp. This drops the now-unused
<cstdlib>/<iostream>/<cmath> includes those switches needed.

The remaining diagnostic switches (ABACUS_LR_EXX_SWAP in
esolver_lr_grad.cpp, ABACUS_LR_CXCO_T in operator_lr_exx.cpp) are
EXX-specific and intentionally left in place, uncommitted.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
CI's legacy Makefile build failed with "No rule to make target
'build/obj/zeq_lin_solv.o'": the OBJS_LR_GRAD list still named the
object zeq_lin_solv.o from before the file-flattening rename, but the
actual source is zeqlin_solv.cpp (zeq_linear_solver.* was renamed to
fit the 15-char filename cap). Fixed the one stale name.

Verified locally with CXX=mpiicpx (icpx 2026.0.0): the previously
missing zeqlin_solv.o now compiles cleanly, both with and without
LIBRI_DIR set. A full end-to-end link isn't reachable on this machine
beyond that -- two separate, pre-existing issues already present
verbatim on origin/develop block it (OBJS_BSE unconditionally lists
utils/lr_io_krlist.o while OBJS_LR's ifdef LIBRI_DIR block adds it
again, and rdmft_pot.cpp references RI_2D_Comm::get_ik_list
unconditionally though it's only compiled in under LIBRI_DIR) -- both
unrelated to this branch and out of scope here.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The agent governance check now rejects #pragma once in new/changed
headers. Converted all 19 headers this branch added under
source_lcao/module_lr/ to the project's #ifndef/#define/#endif
convention, matching the macro naming already used by neighboring
files (ABACUS_SOURCE_<PATH-SANS-source_/UPPERCASE>_H/HPP).

Verified: `agent_governance_check.py --staged` reports no pragma
findings, and abacus_std_para builds clean.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
…move

Commit 4491667 ("derive occupied windows from charged spin
populations") moved nocc's resolution from INPUT-parse-time (a
reset_value lambda clamping it to nbands) to post-SCF, inside
esolver_lr_lcao_tddft.cpp, using actual KS spin populations instead of
a naive nelec/2. It updated the matching assertion in
test_serial/read_input_test.cpp but missed two sibling test files that
still expected the old behavior:

- source_io/test/read_input_ptest.cpp: InputParaTest.ParaRead expected
  an unset `nocc` to already equal `nbands` right after parsing. It
  no longer does -- that resolution only happens later -- so it stays
  at its raw default (-1). Failed in CI as MODULE_IO_input_test_para
  and MODULE_IO_input_test_para_4.

- source_io/test_serial/read_input_item_test.cpp: InputTest.Item_test2
  called nocc's `reset_value` directly, which commit 4491667 removed
  (no longer assigned since that behavior is implemented elsewhere) --
  the call hit an empty std::function and threw bad_function_call.
  Removed the now-meaningless block; there is nothing left to exercise
  at the INPUT-parsing layer. Failed in CI as
  MODULE_IO_read_item_serial.

Verified locally: built an isolated BUILD_TESTING=ON tree (build_test/,
same ENABLE_LIBRI/ELPA/LIBXC config as the main build) so as not to
disturb the running build/ directory, and confirmed all three tests
now pass:
  320 - MODULE_IO_read_item_serial ....... Passed
  327 - MODULE_IO_input_test_para ........ Passed
  328 - MODULE_IO_input_test_para_4 ...... Passed

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Pass matrix storage pointers to the closed/open-shell Pulay grid
integrator instead of taking &(0,0), which asserts on valid MPI ranks
with no local FFT grid. Keep all ranks in the integration and pool
reduction. Both HF gamma force cases retain their existing references.

Disable force in the LDA/PBE gamma cases because the CI distribution
libxc lacks kxc. Regenerate their references with the current XC
cutoffs, retaining excitation-energy checks. Explain in the kernel
comment that commit 40e8e95's cutoff change also affects excitation
energies even when cal_force is disabled. No cutoff or solver threshold
is changed here.

Verification:
- cmake --build build --target abacus_std_para -j8: passed for the
  empty-grid fix; subsequent source change is comment-only.
- bash .diagnostics/ci-lr-4mpi/run-fixed.sh: four sequential MPI runs,
  4 ranks, OMP_NUM_THREADS=1, MKL_NUM_THREADS=1, system libxc: all rc=0.
- Autotest.sh -a /home/fortneu/lr-grad/abacus-develop/build/abacus_std_para
  -n 4 -j 1 in an isolated case copy: rc=0, all 17 numeric checks pass.
  Uses original lr_thr=1e-2 and default comparison thresholds.
- python3 tools/03_code_analysis/agent_governance_check.py --staged:
  no blocking findings; documentation-sync warning only.
- git diff --cached --check: passed.

No user-facing INPUT parameter behavior changes: only integration-case
configuration, an explanatory comment, and an MPI empty-grid fix.
No parameter metadata regeneration is required. Diagnostic source edits,
ScaLAPACK work, and scratch/log files are excluded.
Reuse the grid density for each input spin and vector in open-shell LR and
Z-vector actions, while evaluating each output-spin potential separately.
Keep separate entries for transition and difference density matrices, and
retain the density formulas and numbered grid-integration steps in comments.

Validation: five paired CH3 runs before/after the optimization covered real
Gamma average gradients with 1/4 MPI, real Gamma JT gradients with 4 MPI,
and complex Gamma/two-k spectra with 4 MPI (OMP/MKL threads=1). Printed
energies, transition dipoles, forces and available CG residual histories
were unchanged. The 4-MPI average rho integral count fell from 433 to 237;
vl remained at 442. Restoring comments leaves the tested code unchanged.
git diff --check and the staged governance check passed without blockers.

No INPUT behavior changed, so parameter documentation is unchanged.
The new includes supply the map value type and complete Hxc type needed
for dispatch. Existing unrelated diagnostic changes are excluded.
Replace per-pair density probes with matrix contractions and one MO reduction per k point. Retire the legacy probe implementation and preserve nonsymmetric and complex projection conventions. Validated with six unit tests, Gamma and multi-k HF regressions, and paired LiF single-cell, 222-supercell and 444-k-point CPU runs.
Remove mutable DMR buffers and record readiness only when the owning object is updated. Pass writable transition densities to Hxc and EDM assembly explicitly. No user-facing INPUT or numerical behavior changes; generated parameter documentation is unchanged.

Validation: cmake --build build --target abacus_std_para -j8; cmake --build build_test --target MODULE_ESTATE_dm_cal_DMR_test -j4; OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 ctest --test-dir build_test -V -R ^MODULE_ESTATE_dm_cal_DMR_test$ (4 tests passed); staged governance check (documentation rationale above).
Remove density-matrix getters for orbital layout and k vectors. Both transpose overloads receive the layout explicitly and all four LR callers pass their own dependency. No INPUT or numerical behavior change; no generated parameter documentation update is needed.

Validation: full abacus_std_para build; four LR gamma cases with OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 and 4 MPI ranks, 17 numerical comparisons passed, including both force totals; staged governance check has only evidence/documentation review warnings addressed here.
Keep ESolver responsible for solving Z, constructing dependencies and dispatching force evaluation. Move closed/open-shell density and force assembly, formatted output, amplitude padding and root following to LR free functions with explicit borrowed inputs. No INPUT or numerical behavior changes, so generated parameter documentation is unchanged. Template headers require complete tensor/layout types and reduction declarations; output declarations require matrix and standard stream/container types. Diagnostic code is excluded.

Validation: abacus_std_para build; four 4-MPI gamma cases, 17 numerical comparisons passed; MODULE_LR_grad_degen and MODULE_LR_gradient_amplitudes CTest targets, all 18 tests passed; staged governance and diff whitespace checks (header/documentation review warnings explained above).
Rename CVCX sources and their unit test to lowercase filenames, updating all includes and CMake sources. Split gradient/relaxation orchestration and separate Jahn-Teller algebra/tests so every new C++ source/header in this PR is below 500 lines. No INPUT or numerical behavior changes; parameter documentation needs no update.

Validation: full abacus_std_para build; MODULE_LR_CVCX_test, MODULE_LR_grad_degen, MODULE_LR_gradient_amplitudes via CTest, all 22 tests passed; staged governance and whitespace checks passed with the documentation rationale above.
Make open-shell Z applications and their density setters non-const. Build fixed XC spin combinations in potential constructors, and pass call-local exchange projection buffers explicitly. This removes all mutable members introduced by this PR without hiding writable state behind const methods. No INPUT change; generated parameter documentation is unchanged. Exchange projection now allocates its gradient scratch per call.

Validation: full abacus_std_para build; four 4-MPI gamma cases with OMP_NUM_THREADS=1 MKL_NUM_THREADS=1, all 17 numerical comparisons passed and both force totals unchanged; staged governance and whitespace checks passed with evidence/documentation review rationale here.
Shorten overlong new filenames, excluding extensions, and update includes and CMake source lists. Rename paired projection/Jahn-Teller modules with their tests to preserve discoverability. No INPUT or numerical behavior change; parameter documentation is unchanged. Existing diagnostic blocks remain unstaged.

Validation: full abacus_std_para build; exchange projection, degenerate gradients, amplitude tracking and Z linear-solver CTest targets, all 35 tests passed; all new C++ files relative to origin/develop have lowercase stems <=15 characters and fewer than 500 lines; staged governance and whitespace checks passed with documentation rationale here.
Rename lr_relax_degen_mode to lr_degen_mode and lr_grad_degen_thr to lr_degen_thr, keeping names within fifteen characters. Update internal fields, checks, diagnostics and references. Regenerate docs/parameters.yaml from abacus --generate-parameters-yaml and input-main.md with docs/generate_input_main.py. Old unreleased names are rejected; defaults and numerical behavior are unchanged.

Validation: full abacus_std_para build; --version and -h for both new names; --check-input accepts state/average/jt with valid thresholds and rejects invalid mode, zero threshold with average, and both old names with the expected messages (7 cases); four 4-MPI LR gamma cases, 17 numerical comparisons passed; staged governance and whitespace checks passed.
Include dm_diff.h directly in the extracted force evaluators and in
cal_w_from_z.h before its templates are defined. Previously the declarations
were supplied indirectly by EXX headers; default GNU, without ELPA and serial
builds with LibRI disabled could not resolve cal_dm_diff_*.

Synchronize Makefile.Objects with the lowercase CVCX filenames and the six
new gradient/relaxation translation units, including exx_proj.cpp.

Verification (all passed, including executable linking):
- cmake -S . -B .diagnostics/ci-build-fix/gnu-make -DENABLE_LIBRI=OFF -DENABLE_LIBXC=OFF
  cmake --build .diagnostics/ci-build-fix/gnu-make -j4
- cmake -S . -B .diagnostics/ci-build-fix/noelpa -DENABLE_ELPA=OFF -DENABLE_LIBRI=OFF -DENABLE_LIBXC=OFF
  cmake --build .diagnostics/ci-build-fix/noelpa -j3
- cmake -S . -B .diagnostics/ci-build-fix/serial -DENABLE_MPI=OFF -DENABLE_LIBXC=OFF
  cmake --build .diagnostics/ci-build-fix/serial -j3
- source /opt/intel/oneapi/setvars.sh; export I_MPI_CXX=icpx
  make -f "$PWD/source/Makefile" -j4 BUILD_DIR="$PWD/.diagnostics/ci-build-fix/legacy" CXX=mpiicpx ELPA_LIB_DIR=/usr/local/lib ELPA_INCLUDE_DIR=/usr/local/include CEREAL_DIR=/usr/include/cereal OPENMP=ON
- cmake --build build_test --target MODULE_LR_exx_projection MODULE_LR_gradient_amplitudes MODULE_LR_grad_degen MODULE_LR_zeqlin_solv -j3
  env OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 ctest --test-dir build_test -V -R 'MODULE_LR_(exx_projection|gradient_amplitudes|grad_degen|zeqlin_solv)$'
  Four CTest targets / 35 unit cases passed outside the sandbox.
- All four resulting executables: --version exits 0, v3.11.0-beta10.
- GNU incremental rebuild and staged governance/diff checks passed.

Local GNU is 13.3.0, C++11, using oneMKL and Intel MPI. CMake generator is
Unix Makefiles because Ninja is unavailable; CI containers were not run.

Header dependency rationale: cal_w_from_z.h directly calls cal_dm_diff_*
in template definitions, requiring declarations/implementation before definition
for two-phase lookup, independently of __EXX and caller include order.
No numerical or INPUT behavior changes; parameter docs need no update.
No reference changes or diagnostic code included.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

EXX and lr-TDDFT Related to EXX or lr-TDDFT Features Needed The features are indeed needed, and developers should have sophisticated knowledge

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants