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

[maki49]

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.

2. 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.

Other Behavior changes

• KXC-enabled LibXC build in CI/CD: the build will take longer time.
• XC kernel cutoffs: the density cutoff changes from 1e-6 to 1e-14, and the gradient cutoff from 1e-10 to 1e-20. This can change excitation energies even when forces are disabled. A four-rank A/B run reproduces the LDA gamma reference changes. The periodic Si multi-k case has a minimum spin density of approximately 1.29e-3, above the old cutoff; its four excitation energies agree at 17-digit output precision between the two settings.
• Closed-shell singlet weighting: nspin=1 singlet response uses the factor-two Hartree/XC contribution, while ground-state S1_gs retains weight one. Existing single-spin excitation results can therefore change. BSE grid initialization uses S1_gs because HamiltBSE already applies the singlet factor; the existing RI BSE path is unaffected by this grid correction.
• HSE06 kernel regularization: the LR HSE06 implementation clamps sigma to at least 1e-6 when evaluating XC derivatives. This can change HSE06 response energies and gradients near vanishing density gradients; it does not apply this clamp to ordinary PBE.
• Composite XC third derivatives: contributions from multiple XC functionals are accumulated when evaluating kxc, rather than allowing the last functional to overwrite earlier contributions. This affects the XC contribution to analytic excited-state forces.
• Davidson MPI broadcast: eigenvector broadcasts use the actual scalar type through the guarded wrapper, replacing an unconditional complex-double broadcast. This fixes real-valued MPI buffer handling and can affect both ground-state and LR calculations using this solver.
• Gradient orbital window: the Z-vector calculation uses all available ground-state virtual bands, while lr_nvirt continues to select the response basis. The gradient can depend on the available ground-state band window even when the response basis is held fixed.
• Failure handling: incomplete full-window wavefunction reads and unconverged Z-vector CG solves now fail explicitly instead of continuing to produce a force. CG convergence is checked using the final residual; its iteration limit and tolerance remain internal constants.
• MPI and buffer robustness: empty local grid matrices, spin-resolved EXX density-matrix buffers, and the single-channel complex spectrum buffer are handled without invalid indexing or incorrect allocation assumptions.
• Force output and diagnostics: excited-state force output and test force summaries retain ten decimal places. Full Z-vector dumps are restricted to test_force runs.

[Copilot]

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

[mohanchen]

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]

@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.

[ABACUS-PR-BOT]

Reviewed base 65a324ba3e5156c31f2818ccd93cd564985305ad -> head ab811c9b506165c7b1f1f3410c9306568672c461. The inline findings distinguish reproduced failures, small mathematical probes and conditional risks.

Follow-ups belong to the existing HSE, spin-basis tracking and operator-construction threads.

Remaining governance items from the earlier review:

• Existing interfaces still gain default arguments: ao_to_mo.h:22,32,47 and operator_lr_hxc.h:42-44. Use explicit arguments or a clear overload.
• New comma declarations remain, e.g. lr_grad_os.cpp:46 and force_funcs.h:29; declare one variable per line.
• zeq_solver.h:9 propagates a new implementation .hpp; move implementation out or record the narrow exception. The global-reference budget decreases by 32, so it is not a net-increase blocker.
• git diff --check reports 33 whitespace errors; clean only the added/changed lines.
• Non-blocking: correct the PR body's lr_relax_degen_mode / lr_nvirt to lr_degen_mode / nvirt.

Verification: MPI/LibXC/ELPA build without LibRI succeeded; 15 selected CTests, five additional four-rank tests, eight Python tests and 44 production-TU C++11 syntax checks passed. Independent finite differences and failure probes exposed the issues above. Governance: 0 mechanical errors, 216 warnings; 103/120 changed C++ files score >=60. Generated INPUT docs and CLI checks passed.

Not verified: EXX/HF/HSE/LRC force runs and BSE runtime (LibRI configuration lacked cereal), production open-shell finite differences, GPU/serial full builds, supercell replication and cutoff-boundary derivatives. CI status was unavailable.

[mohanchen]

I will mention again here: the esolver directory currently contains too many implementation details specific to LR algorithms. For instance, cal_nupdown_form_occ—along with several other functions—should not reside here and needs to be relocated or refactored. Please streamline the esolver_lr code to significantly reduce its line count and clutter."

Furthermore, while there are still some legacy LCAO-related codes in esolver that shouldn't be there, that is not a justification to follow suit, especially since I don't have the time to clean up the LCAO parts myself. Leaving so many LR details here only increases my workload, and realistically, you are the one best positioned to help keep this clean."

[mohanchen]

To follow our design specifications: esolver should serve as the cockpit control panel of a car—it needs to remain as clean and simple as possible. While a car naturally relies on countless intricate parts and wiring to function, these internals should be encapsulated under the hood, not exposed on the dashboard.

The purpose of esolver is to present a high-level physical workflow, not the granular mechanics of how every individual physical quantity is computed. It is sufficient for the reader to know that a quantity has been computed.

Furthermore, we enforce strict naming conventions for functions inside esolver—you cannot simply introduce arbitrary helper functions. Just like modern vehicle control panels maintain standard layouts so drivers can transition seamlessly between different cars, esolver maintains a standardized interface across solvers. Unsanctioned additions degrade the readability and consistent experience that esolver is designed to provide.

[maki49]

@mohanchen Refactored in 398a482. All the newly added .hpp files are changed into .h+.cpp. exx_lri.hpp is an existing file and I only change a few lines, so refactor it is beyond the feature's scope. Some member functions are extracted outside from ESolverLR. Ready for review again.

[mohanchen]

Not good to use PARAM here— new code should avoid this. The reasoning comes from the project governance rules and testing practice:

1. Violates the cross-layer control rule

Governance rule 1 explicitly forbids increasing cross-layer control throughGlobalV ,GlobalC , orPARAM . ReadingPARAM.inp.xxx in a member initializer list implicitly injects the global input configuration into the class — nothing in the constructor signature reveals which INPUT parameters the class depends on.

2. Poor testability

To drive such a class in a unit test, you must set up global PARAM first, or resort to a#define private public hack — and the governance checker now blocks any net increase of such hacks. The fix pattern named in the governance doc is exactly this: pass INPUT values as explicit arguments (seeRelax_Criteria andK_Vectors::read_kpoints ).

3. Hidden initialization-order risk

Member initializer lists execute before the constructor body, so the code silently requires "PARAM must be filled before construction" — with no mechanism to enforce or even see that precondition.

Recommended pattern : read PARAM once at the wiring point (e.g., Driver / ESolver construction), assign to named local variables, then pass them explicitly to the constructor; the initializer list only does direct parameter passthrough (this also satisfies governance rule 13: no non-trivial function calls in the initializer list).

[mohanchen]

Do not callpgemm ordgemm directly; you must use the math library connectors provided insource_base .

[mohanchen]

Do not callpgemm ordgemm directly; you must use the math library connectors provided insource_base .

[mohanchen]

Do not callpgemm ordgemm directly; you must use the math library connectors provided insource_base .

[mohanchen]

Do not callpgemm ordgemm directly; you must use the math library connectors provided insource_base .

[mohanchen]

use connector

[mohanchen]

use connector

[mohanchen]

use std::vector instead of new/delete?

[mohanchen]

use std::vector instead of new/delete?

[mohanchen]

. Harder to Debug (the most practical issue)

• Breakpoints are set per line: the entire line gets a single breakpoint — you cannot stop at the loop body or theo += ... accumulation statement to inspect intermediate values
• Step-through granularity is per line: multiple statements on one line cannot be stepped individually
• Crash/out-of-bounds stacks only report line numbers: this line contains index arithmeticis * nk + ik , container accesswg(...) , and loop control — an error message cannot immediately tell you which step failed

2. Version-Control Unfriendly

• Coarse diff granularity : changing any single character (e.g., renamingnk ) marks the whole line red in the diff, forcing reviewers to hunt character-by-character for the actual change
• Polluted git blame : any later small edit overwrites blame ownership for the entire line, losing the "who wrote which part" history
• Higher merge-conflict probability : parallel edits to different parts of the same line always conflict and cannot be auto-merged by git

3. Readability and Inspectability

• Readers must parse statement boundaries themselves: the lambda signature, loop header, loop body, and return value — four layers of logic crammed together
• Expressions likeik < nk andis * nk + ik lack spacing and named intermediate variables, so operator precedence must be mentally verified — the same spirit as project rule 1. modify some variable name in LCAO_gen_fixedH::build_Nonlocal_beta() #14 (no conditional expressions inline in argument lists; assign to a named local first): every value should be inspectable
• Static analysis / coverage tools (gcov, llvm-cov) count per line; multiple statements on one line distort coverage data

4. Higher Cost of Modification and Refactoring

• Adding a log statement, an assertion, or extracting the loop body requires manually re-wrapping and reformatting first — a "one-line addition" becomes a "rewrite of a whole block"
• Compact writing hides small errors (e.g., the missing space before++ik here makes the loop structure easy to miss when skimming)

5. An Independent Issue in This Line: [&] Capture-All

The capture list says nothing about what the lambda actually depends on (in realitywg andnk ); the reader must read the entire enclosing function to confirm it does not accidentally capturethis or other mutable state. Writing[&wg, nk] would be far clearer.

[maki49]

@mohanchen All fixed in 17d5a66:

• Expanded the lambda and population loop into separate statements, named the spin/k-point row index, and replaced capture-all with [&wg, nk].
• The real and complex eigenproblem workspaces now use std::vector and pass .data() to LAPACK. No raw new[]/delete[] remains in lr_util.cpp.
• The EDM contractions now use ScalapackConnector::gemm and BlasConnector::gemm_cm, with shared real/complex implementations. Both real and complex LU solves now call container::lapackConnector::getrf/getrs, the existing column-major wrappers recommended by lapack_connector.h.
• Removed the global-reading kernel_nspin() helper. PotGradXCLR now receives the spin count explicitly from its LR callers and checks it against the supplied kernel metadata. Its base constructor only receives direct parameters.

Ready for review again.

[mohanchen]

For issues that a reviewer has pointed out before, please learn to draw inferences from one instance — in other words, fix all similar problems throughout the code as well, because the reviewer may have missed some of them too. Reviewing more than ten thousand lines of code is already an enormous burden on the reviewer, so please don't count on anything slipping through by luck. The less thorough your fixes are, the more they erode the reviewer's trust. Merging poor-quality code into the community does not reduce the submitter's workload, nor is the reviewer the only one who suffers — it harms the interests of everyone in the community who contributes to the ABACUS codebase. Please keep this sense of responsibility in mind.

[mohanchen]

For example, when I say that functions likepzgemm_ should be wrapped with a connector, you should check your own submitted code for similar issues, instead of only fixing the exact line pointed out in the review — doing only that is mere formalism. The more you do this, the more the reviewer becomes convinced that this code will never be improved in the future, because your goal is to get the code accepted rather than to make it better. Likewise, one-liner-style code appears in many places — please fix all similar occurrences as well. Please hold onto the aspiration of writing higher-quality code, not just getting your code accepted faster.

[mohanchen]

For issues that a reviewer has pointed out before, please learn to draw inferences from one instance — in other words, fix all similar problems throughout the code as well, because the reviewer may have missed some of them too. Reviewing more than ten thousand lines of code is already an enormous burden on the reviewer, so please don't count on anything slipping through by luck. The less thorough your fixes are, the more they erode the reviewer's trust. Merging poor-quality code into the community does not reduce the submitter's workload, nor is the reviewer the only one who suffers — it harms the interests of everyone in the community who contributes to the ABACUS codebase. Please keep this sense of responsibility in mind.

[mohanchen]

please use connector

[mohanchen]

use connector

[mohanchen]

use connector

[mohanchen]

use connector

[mohanchen]

use connector

[mohanchen]

Please break one-liner code into multiple lines.

[mohanchen]

Please break one-liner code into multiple lines.

[mohanchen]

Please break one-liner code into multiple lines.

[mohanchen]

Please break one-liner code into multiple lines.

[mohanchen]

Please break one-liner code into multiple lines.

[maki49]

All fixed throughout the additions:

• e74163267: use connector.
• 9bc3be36d: explicit lambda capture and static_cast.
• 0a212fda7: expand the one-line codes.

For static_cast, the capacity comparisons do not require narrowing the container size to int: the cast is redundant when the size fits in int and can produce an implementation-defined result when it does not. The comparisons now use a std::size_t local initialized from the nonnegative grid size int nrxx, without an explicit cast. Casts at interfaces that actually require int dimensions are treated separately.

Ready for review again.