Repository navigation
feat: add BSE module - #7718
Merged
Merged
feat: add BSE module#7718
Conversation
mohanchen
added
Feature Discussed
mohanchen
reviewed
Jul 31, 2026
mohanchen
added
the
Tests/Examples
Collaborator
|
Just a note: the current code is not well encapsulated:
|
Fisherd99
force-pushed
the
develop
branch
7 times, most recently
from
August 12, 2026 20:41
124ce0e to
1962112
Compare
2. feat: add ExcitonPlotter 3. add input parameter `rpa_outdir` and `rpa_out_vel`
…outputs, ensuring robust handling of parallel I/O. - Refactor exciton plotting utilities. - Delete unnecessary 'PARAM‘ and 'GlobalV'.
fix without MPI compile bug fix elpa interface
mohanchen
added
Features Needed
mohanchen
approved these changes
Aug 13, 2026
mohanchen
left a comment
mohanchen
left a comment
Collaborator
There was a problem hiding this comment.
This is an important update, thanks for your contributions! However, keep in mind that refactors are needed in near future.
dyzheng
pushed a commit
to dyzheng/abacus-develop
that referenced
this pull request
Aug 18, 2026
…r-sign Conflict resolutions (build fixes required by the merged tree): - module_operator_lcao/test/CMakeLists.txt: ORB_atomic_lm.cpp -> orb_atomic_lm.cpp, matching develop's filename normalization (deepmodeling#7797) in the dm_trace test target sources. - module_lr/utils/exciton_plotter.cpp: add missing 'template' keyword for the dependent-template call eval_wfc_bloch<T> (upstream BSE module, deepmodeling#7718) that fails to compile with the standard C++11 rule.
maki49
added a commit
to maki49/abacus-develop
that referenced
this pull request
Sep 4, 2026
The LR part was being subtracted. It should be added: what cal_force(int) returns already follows the ABACUS force convention F = -dE/dR, despite the "Gradients of each excited state" heading it prints under. The finite-difference reference it was validated against, abacus-fd lr-custom, goes through run_diff_custom_lr, where energies[0] is the +dx/2 point and the result is (energies[1] - energies[0])/dx, i.e. -dE/dx -- so that table compares force to force. (run_diff_all_kslr in the same file uses the opposite ordering; do not read the convention off that one.) With the wrong sign a relaxation still converged, because F_gs - F_Omega does reach zero somewhere -- just not at a stationary point of E_gs + Omega. On BH/lda it stopped at 1.2864 A while the energy minimum was near 1.25 A. What made that visible is the other half of this commit: one line per ionic step giving E_gs, Omega, E_exc and the two force magnitudes. Before it the log carried only a force, so a force that was self-consistently zero looked like convergence and nothing could be cross-checked against the energy. After the fix, BH/lda relaxes the first singlet to 1.25628 A, which is the lowest E_exc of every geometry sampled and where the two force halves cancel (0.29326 vs 0.29359). The ground state relaxes to 1.27430 A, so the excited bond is shorter by 0.018 A -- the direction and size seen for BH A1Pi (1.2195 vs 1.2324 A). Also restores the aims guard on set_atomic_trace. d4fe3fe ("Support different basis number from aims") deliberately replaced the unconditional call with a guarded one, because with aims_nbasis the per-atom orbital counts behind iat2iwt do not match nbasis. 1e4c1c6 (BSE, deepmodeling#7718) then re-added an unconditional call above it, silently undoing that. The previous commit here kept the wrong one of the two. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu
mohanchen
added a commit
that referenced
this pull request
Oct 18, 2026
…d excited-state geometry relaxation (#8069) * LR-TDDFT analytical gradients: port onto refactored develop Port the LR-Grad series (44 commits, 7c6398d3..6aea1fc0, 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. * fix(lr-grad): replace LibRI RI::LR with a local two-density-matrix EXX 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. * fix(lr-grad): correct nspin in the LR grid-force ModuleGint calls, close 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. * fix(lr-grad): gather the RHS before the LAPACK Z-vector solve 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). * fix: singlet gate for gxc * fix: Hxc dm-trans term (symmetrization and factor) * Fix some factor bugs in ceZc, K_ii terms (visible in DZP) * Fix: transpose dm_trans in Z-vector eq. RHS * Fix pot_hxc_gs: correct TDRPA@LDA (linear vxc=fxc[gs](T+Z) instead of vxc[T+Z]) * Feature: implement gxc * Perf: share gxc kernels to save memory * Perf: put drho out to cal_v_eff for gxc, to save memory * fix nspin=1 * small fixes * Perf: speed up the LR kernel-to-potential integrands by ~45% `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 * feat: lr-grad openshell implementation * Feat: enable LRC kernel * Remove the LR-Grad EXX debug scaffolding `lr_grad_debug.h` held temporary ablation switches (LRDBG_EXX_ZRK / _KEXX / _FDMT) read from the environment. Nothing includes it any more, and it was only tracked because it got swept into a83cbb062 during the rebase onto develop. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * fix(lr): make the KS-orbital window copy and paraX_ setup re-entrant Three defects that all only bite once ESolver_LR runs more than once, or runs without MPI: - The serial branch of the [nocc, nvirt] window extraction computed the source and destination pointers and then copied nothing, leaving psi_ks uninitialized in every non-MPI build. - setup_eigenvectors_X() only ever emplace_back'd into paraX_, so a second call would double its size. - psi_ks / psi_ks_all were raw new with no matching delete before reassignment; they are unique_ptr now, which also drops the two deletes from the destructor. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * refactor(lr): keep the KS solver as a member and alias its geometry objects ESolver_LR built the ground state in a stack-local ESolver_KS_LCAO and then moved its innards out. That is a dead end for anything that has to follow moving ions, and it also left the excited-state force unable to run at all on the ks-lr path: initialize_from_ks_ never calls ESolver_FP::before_all_runners and only moved pw_rho, so pw_rhod stayed null and Pgrid / sf / locpp stayed empty -- all four are required by LR_Force and by init_pot_groundstate. The solver is now a member (ks_), and the objects it already builds for the current geometry are aliased instead of moved or rebuilt: pw_rho/pw_rhod/pw_big directly, and Pgrid / sf / locpp / gd / two_center_bundle_ through accessors, since ESolver_FP holds the first three by value. psi_ks_all aliases ks_->psi rather than stealing it -- it is only read here, and it is the one array too large to copy per step. pw_rho_flag stays false so the borrowed basis is never freed here. kv, eig_ks_all, wg_ks_all and paraMat_.atom_begin_* are copied rather than aliased: they are small, and the KS solver overwrites its own pelec on the next step. Two bugs fall out of the aliasing: - two_center_bundle_ was only populated for abs_gauge = velocity, but LR_Force and cal_hs_grad use it unconditionally, so the length gauge computed forces from an empty bundle. - wg_ks_all was never filled on the ks-lr path, though cal_dm_gs and test_force both read it. init_pot_groundstate no longer repeats sf.setup / locpp.init_vloc when ks_ is present; ESolver_FP::before_scf has already done both. ks_->runner is still called from before_all_runners -- moving it into runner(istep) is the next step. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * feat(lr): re-run the ground state from runner(istep) before_all_runners is called once, outside the relaxation loop, so a ground state computed there can never follow moving atoms. The KS SCF moves into runner(istep), which ESolver_KS::runner already supports: its before_scf rebuilds the neighbour lists, grid tables and Hamiltonian for the current geometry. initialize_from_ks_ splits accordingly. What survives a change of geometry stays there -- dimensions, the 2D distributions, the psi_ks / eig_ks buffers, pelec, orb_cutoff_, and the Exx_LRI object itself. What does not moves into a new refresh_from_ks_, run once per ionic step: the psi_ks_all alias, the eig_ks_all / wg_ks_all copies, the [nocc, nvirt] window extraction, init_pot, cal_exx_ions, and the re-publication of the static gint pointer. Three things that only a second call exposes: - reset_dim_spin2 shifts nocc/nvirt between the spin channels, so it compounds if repeated. It is now guarded to the first step, and a later change of the ground-state spin population is a hard error rather than a silent dimension mismatch. - init_pot used pot.resize, which keeps existing elements, so the second geometry would have kept the first one's potentials. It assigns now. - move_exx_lri nulled the ground-state solver's exx_ptr after taking it. That was harmless while the KS solver was a temporary; now it would leave the KS solver without EXX from the second step on. It is share_exx_lri now and shares rather than steals. Known gap, left for the commit that wires the force into the ESolver interface: initialize_from_ks_ only builds an Exx_LRI when xc_kernel is a hybrid, while initialize_from_unitcell_ also builds one for cal_force with a hybrid ground state. The ks-lr path therefore has no exx_lri for a local kernel on top of a hybrid ground state. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * feat(lr): drive geometry relaxation from an excited-state surface ESolver_LR computed excited-state gradients but had no way to hand them to anything: the ESolver interface overrides were no-op stubs and cal_energy returned zero, so Relax_Driver::esolve came back with a zero force matrix. They are implemented now, for calculation = relax with esolver_type = ks-lr. Sign and composition, since neither is obvious from the code: - cal_force(int ispin) returns the *gradient* +d(Omega)/dR, matching the "Gradients of each excited state" heading and the finite-difference reference in abacus_fd (core.py:615). ESolver::cal_force must return a force, F = -dE/dR (ions_move_basic.cpp:45 negates it back), hence the minus sign. - The LR terms only carry Omega. The ground-state force is a separate part of the excited-state gradient and comes from ks_->cal_force, so the total is force_gs_ - lr_grad_. Likewise cal_energy returns etot_gs_ + Omega, which is what the line searches in cg / bfgs / lbfgs need. Outside a relaxation it still returns zero, so existing single-point output is untouched. A relaxation follows one state, and solving the Z-vector equation dominates the cost of a gradient, so solve_zvector_eqation, cal_force and cal_force_openshell take an istate_only argument. Z_vector_equation already treats X and Z as nstates independent blocks, so one state is that block with nstates = 1; X keeps the real state index and Z is indexed within the solved range. istate_only = -1 keeps the all-states behaviour that single-point runs use. New input: lr_target_state (0-based) and lr_target_spin (singlet/triplet/updown). Whether the calculation is open- or closed-shell is only known after the ground-state occupations are read, so that half of the validation lives in setup_relax_target_, together with the gamma-only requirement -- there is no complex Z-vector solver. Rejected at input parsing, each with the reason in the message: cell-relax (no excited-state stress), md (no non-adiabatic couplings, and a fixed state index would walk straight through a conical intersection), esolver_type=lr with relax (its ground state belongs to one fixed geometry), and lr_solver=spectrum with relax. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * fix(lr): scope the relax-target inputs, and clear out the re-entrancy leftovers lr_target_state / lr_target_spin only steer a relaxation, but their check_value ran for every calculation, so a stray lr_target_spin in a single-point INPUT could abort a run the parameter has no effect on. Both checks now return early unless calculation = relax. An open-shell calculation solves a single spin-conserving channel, so lr_target_spin has nothing to select there: it no longer refuses singlet/triplet, it relaxes that one channel and only reports an explicit `triplet` as ignored. Closed shell keeps the error for `updown`, where singlet and triplet really are different states with different gradients and the choice is ambiguous rather than redundant. Remaining re-entrancy and correctness leftovers, all exposed by running the solver more than once: - spin_types was assigned inside three separate solver branches, and the open-shell `spectrum` branch assigned it nowhere, leaving it empty for after_all_runners and the gradient to index. It is set once now, from `openshell`, before any branch runs. - The excitation energy and amplitude files carry a _step<n> suffix during a relaxation; otherwise every ionic step overwrote the last and the trajectory could not be inspected afterwards. Single-point file names are unchanged. - initialize_from_unitcell_ derived nupdown from pelec->wg, which read_ks_wfc never fills -- it fills wg_ks. A spin-polarised ground state read from file therefore always took the closed-shell branch. - pelec was newed twice there, leaking the first one. - set_parallel_orbitals_band called set_atomic_trace twice, the second time behind a condition the first had already made moot. - setup_eigenvectors_X built a local spin_types that nothing read. - zeq_solver.h declared a Z_vector_equation overload with an older signature and no definition anywhere; a call matching it would have failed at link time. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * fix(lr): the excited-state force adds the LR part, and log the energy The LR part was being subtracted. It should be added: what cal_force(int) returns already follows the ABACUS force convention F = -dE/dR, despite the "Gradients of each excited state" heading it prints under. The finite-difference reference it was validated against, abacus-fd lr-custom, goes through run_diff_custom_lr, where energies[0] is the +dx/2 point and the result is (energies[1] - energies[0])/dx, i.e. -dE/dx -- so that table compares force to force. (run_diff_all_kslr in the same file uses the opposite ordering; do not read the convention off that one.) With the wrong sign a relaxation still converged, because F_gs - F_Omega does reach zero somewhere -- just not at a stationary point of E_gs + Omega. On BH/lda it stopped at 1.2864 A while the energy minimum was near 1.25 A. What made that visible is the other half of this commit: one line per ionic step giving E_gs, Omega, E_exc and the two force magnitudes. Before it the log carried only a force, so a force that was self-consistently zero looked like convergence and nothing could be cross-checked against the energy. After the fix, BH/lda relaxes the first singlet to 1.25628 A, which is the lowest E_exc of every geometry sampled and where the two force halves cancel (0.29326 vs 0.29359). The ground state relaxes to 1.27430 A, so the excited bond is shorter by 0.018 A -- the direction and size seen for BH A1Pi (1.2195 vs 1.2324 A). Also restores the aims guard on set_atomic_trace. d4fe3fe84 ("Support different basis number from aims") deliberately replaced the unconditional call with a guarded one, because with aims_nbasis the per-atom orbital counts behind iat2iwt do not match nbasis. 1e4c1c6af (BSE, #7718) then re-added an unconditional call above it, silently undoing that. The previous commit here kept the wrong one of the two. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * fix(lr): build an Exx_LRI on the ks-lr path for a hybrid ground state An Exx_LRI is needed for two independent reasons: the LR exchange kernel when xc_kernel is a hybrid, and the ground-state EXX terms of the gradient (the H_gs[T+Z] and W multipliers, and the EXX force term, all gated on gs_is_hybrid()). initialize_from_unitcell_ has always tested both; initialize_from_ks_ tested only the first, so a local kernel on top of a hybrid ground state -- the rpa_at_hse / rpa_at_hf tiers -- had no exx_lri at all on the ks-lr path while those terms went looking for one. Only the from-file path was correct, which is why the benchmark table never showed it. The share test widens too. Both objects are constructed from the same info_ri.coulomb_param, which input_conv derives from dft_functional alone, so whenever the ground-state solver holds one of the right type it is the same object we would build -- and it is already current for this geometry, which also saves a cal_exx_ions per ionic step. The previous `xc_kernel == dft_functional` guard would have forced a redundant from-scratch build in exactly the case being fixed. warn_if_kernel_differs_from_gs is self-guarding, so it moves out of the else branch and covers both. No change where a local kernel sits on a local ground state: both halves of the condition are false, as before. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01GUvgXxW4C5E1sRzUXqegEu * fix: use an spin-index mixing in close-shell grad: affect triplet gradients with degenerate KS basis up ≠ down * fix: expand Z-window to nbands for complete orbital rotation space * fix: small fixes about nupdown and degeneracy information * fix openshell bug: parallel index in HamiltULR and add-on op-lr-diag * fix: adapt LR-gradient code to develop interface changes after rebase Four upstream changes that the merge algorithm could not see, all found by the compiler rather than by a conflict: - read_wfc_nao gained a 'binary' argument (#7863). Added at both LR call sites; the second one, in the cal_force branch, is our own code and was not flagged as a conflict. - BaseCell::Kind::unit_cell was renamed to 'unitcell'. - DensityMatrix::cal_DMR is now const, so the '_dmr_ready' cache flag it sets must be mutable -- the same treatment '_DMR' already documents. - HamiltLR's constructor now takes 'in_dir' as well as 'out_dir' (#7849). Both Z-vector call sites passed only the output directory. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * fix compile warning * fix: clamp sigma for HSE06's wpbeh kernel; fix complex CXC index bug HSE06's short-range exchange (libxc's gga_x_wpbeh) evaluates its 2nd/3rd sigma-derivatives through a closed form that loses all precision via catastrophic cancellation near the reduced gradient s->0 (a symmetry-forced grad-rho=0 node at otherwise ordinary density). The true derivative is finite there; libxc's closed-form evaluation is not. Clamp sigma from below (sigma_cut=1e-6) before both the xc_gga_fxc and xc_gga_kxc calls whenever the functional being evaluated is XC_HYB_GGA_XC_HSE06 -- clamping only one of the two orders is insufficient and can be worse, since it breaks the partial cancellation that existed between them at the raw point. Plain GGA functionals (e.g. PBE) and xc_kernel=rpa are unaffected: sigma_cut=-1 is a true no-op since sigma=|grad rho|^2 >= 0 always, and rpa never reaches this switch at all. Verified against archived finite-difference references for 06_N2/hse, 09_CH4/hse and 08_BeH2/rpa_at_hse (MAE/|F| 0.015%-0.556%, matching functionals that were never broken) and confirmed to remove the divergence in 08_BeH2/hse itself, at both its DZP and TZDP bases. OperatorLREXX<complex<double>>::cal_DM_onebase's CXC term1 and CXC_o branches had their two band indices swapped relative to the real specialization and to their own comments: term1 read Co's band index against Cv's full-space slot and Cv's virtual index against CvX's occupied-only row count (out of range once nvirt > nocc), and CXC_o read (io, iv) where it needs (iv, io) to match the coxt_full/psi_ks_full argument order used to construct it. Casida (Omega, X) for a k-point (gamma_only=0) run is otherwise untouched by this file; the complex Z-vector solver remains unimplemented and still raises on the gradient path. * chore: remove diagnostic env-var switches from the LR-Grad investigation Removes every ABACUS_LR_*/ABACUS_GXC_* getenv-gated diagnostic branch, scale factor, and dump/print accumulated while tracking down the HSE sigma-clamp fix and the complex-path CXC index bug: Z_SCALE, T_SCALE, K_FDMT, EDM_T4, EDM_W, EDM_Z, EDM_OMEGA, K_CXCO_LOC, K_CXCO, ZR_NOTRANS, K_CXC, R_HT, ZSOLVER, DUMP_L, A_SYM, DUMP_A, ANA_THR, KOO_SYM, CXC_T1, CXC_T2, KERNEL_ALPHA, DIAG_CONTRIB, VEFF_DIAG (plus its g_veff_diag_istate global), DIAG, RHO_THR, SIGMA_THR, RHOMIN_THR (and the mask_gxc_coef function that housed the last four, which masked the wrong variable -- see the governance log), CXCO_T, EXX_SWAP, and ROT. Each removal keeps whatever the default (unset) behavior already was: the env var always multiplied a production quantity by 1.0 or took an already-default branch when unset, so none of this changes numerics. Verified by rebuilding and rerunning 08_BeH2/hse (DZP, full window) after the cleanup: forces match the pre-cleanup run bit for bit. * feat(lr-grad): algebra for the degenerate-subspace gradient matrix At a d-fold degeneracy no single excited state has a gradient vector: the branch slopes along a displacement u are the eigenvalues of M(u) = sum_a u_a G^(a), and its eigenvectors depend on u. The full first-order object is the 3N matrices G^(Aalpha)_kl = <X_k|dA/dR|X_l>, of which the existing per-state gradient only ever evaluates the diagonal, in whatever basis the diagonalizer returned. Restricted to one multiplet the existing gradient pipeline is a quadratic form in X (every ingredient is built from X (x) X, the CPSCF operator L depends only on the ground state, and Omega is constant inside the multiplet), and every X in the degenerate subspace is itself a legitimate eigenvector with the same Omega. So F[X] = G[X,X] on the whole diagonal, and a symmetric bilinear form is determined by its diagonal: the off-diagonal elements follow from the polarization identity without any new physics code. This commit adds only the basis-independent bookkeeping that route needs, kept free of the parallel layout and the Z-vector solver so that it is unit-testable without a ground state: - group_degenerate_states: split states into multiplets, anchoring each group to its lowest member rather than chaining consecutive gaps, so the spread inside a group is bounded by the threshold however many members it collects. - degenerate_pairs: the d(d-1)/2 pairs, fixing the order the combinations are evaluated and assembled in. - combine_normalized: X_+ = (X_k + X_l)/sqrt(2), normalized so that it is a genuine eigenvector and the untouched per-state path can evaluate it. - assemble_grad_matrix: G_kl = F[X_+] - (G_kk + G_ll)/2, reusing the d diagonal gradients already computed, so the whole matrix costs d(d+1)/2 evaluations -- exactly its number of independent components. Derivation: LR-Grad-formulas/2026-09-简并激发态梯度-实测和讨论.md section 5.4. Verification: new test_grad_matrix_degenerate, 10 tests, all passing. Built and run standalone (g++ -std=c++14, gtest 1.x, openblas) since the repository build was in use: [==========] 10 tests from 1 test suite ran. (0 ms total) [ PASSED ] 10 tests. The two substantive cases do not test the arithmetic but the route: a synthetic quadratic form X^T B X stands in for the gradient pipeline, and the tests check that the combinations reproduce G_kl = X_k^T B X_l to 1e-12 and that the result is covariant under a rotation of the subspace basis, G' = U^T G U. That is the self-check the document lists as the strongest one available without finite differences. Governance check: 4 warnings, no blockers. The three header includes (<cstddef>, <utility>, <vector>) are required by value in the declarations (size_t, std::pair, std::vector), not reducible to forward declarations. No documentation update is required: this commit adds no INPUT parameter and changes no existing behaviour -- nothing calls these functions yet. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * refactor(lr-grad): let the excited-state gradient take X from the caller `cal_force` read the excitation vectors straight out of the `X` member and looked up each Omega by state index, so the only gradients it could produce were those of the stored eigenvectors. Computing the degenerate-subspace gradient matrix needs the gradient of linear COMBINATIONS of the members of a multiplet, which are themselves perfectly legitimate eigenvectors with the same Omega. Split the existing body out into cal_force_Xz(ispin, Xz, omega, label_begin) cal_force_openshell_Xz(Xz, omega, label_begin) which take blocks already widened into the Z window plus their Omega, and leave `cal_force(ispin, istate_only)` / `cal_force_openshell(istate_only)` as thin wrappers that widen the stored X and look up Omega. `solve_zvector_eqation` now takes the block count `nst` directly instead of re-deriving it from `istate_only`. `ist_begin`/`ist_end` were kept as the loop bounds, now labels only, so the 220-line body is untouched apart from the Omega lookup becoming `omega[istate - ist_begin]`. No behaviour change: this is a pure plumbing split. Verification: built the same source tree twice in the same build directory, once with these two files reverted to HEAD and once with them applied, then ran both binaries on fullwin/02_Li2/lda (nbands 14, nvirt 11, 5 states, gamma, LDA kernel, OMP_NUM_THREADS=1, serial, systemd-run MemoryMax=8G): diff forces_pre.txt forces_post.txt -> IDENTICAL diff eig_pre.txt eig_post.txt -> EIG IDENTICAL The printed force block agrees down to the 1e-15 components whose exact value is zero; those digits are sensitive to any change in operation order, so they are the evidence that nothing was reordered. Full build: exit 0, no warnings from the changed files. Governance check: 2 warnings, no blockers. No test is added here because the commit changes no behaviour and the numerical evidence above is the check that applies -- the unit-testable algebra it exists to serve was tested in the preceding commit. No documentation update is required: no INPUT parameter and no user-visible behaviour changes. Three comments forward-reference `cal_grad_matrix_degenerate`, the driver that lands in the next commit, to say why the seam exists. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat(lr-grad): compute the gradient matrix of a degenerate multiplet Adds `lr_grad_degen_thr` (Ry, default 0 = off). When set, states whose excitation energies lie within it are grouped into a multiplet and the full gradient matrix G^(Aalpha)_kl = <X_k|dA/dR_Aalpha|X_l> is assembled for it, on top of the per-state gradients. Route (A2) of LR-Grad-formulas/2026-09-简并激发态梯度-实测和讨论.md section 5.4: no new physics code. Inside a multiplet the existing gradient pipeline is a quadratic form in X and every X in the subspace is a legitimate eigenvector with the same Omega, so F[X] = G[X,X] on the whole diagonal and a symmetric bilinear form is fixed by its diagonal. The off-diagonal elements then follow from the polarization identity, G_kl = F[(X_k+X_l)/sqrt(2)] - (G_kk + G_ll)/2, which reuses the d diagonal gradients already computed. Cost: d(d-1)/2 further Z-vector solves per multiplet, all of them in one solve. The preconditions are reported rather than assumed: each multiplet's actual Omega spread and max |<X_k|X_l> - delta_kl| go to the running log, with a warning above 1e-6. A loose threshold cannot distinguish a true degeneracy from an accidental near-degeneracy -- where the states have genuinely different Omega and the construction does not apply -- and the reported spread is what settles it. Also printed is Tr(G)/d, the one smooth basis-independent 3N vector field the multiplet has, which is what a symmetry-constrained relaxation should follow. Verification (fullwin/02_Li2/lda, gamma, LDA kernel, nbands 14, nvirt 11, 5 states, OMP_NUM_THREADS=1, serial, systemd-run MemoryMax=8G): 1. Default path unchanged. With `lr_grad_degen_thr` unset, the printed force block is byte-identical to the binary built before this series: diff forces_pre.txt forces_step3.txt -> IDENTICAL and no degeneracy output is emitted (grep count 0). 2. Preconditions on the exactly degenerate pair (states 2, 3, the pi pair): Omega = 0.393945 Ry, spread = 2.22e-15 Ry, max |<X_k|X_l> - delta_kl| = 7.77e-16 3. The resulting G is symmetry-correct, which is the physics check: - along the bond axis (z): G = 0.0519746 * I, off-diagonal 1.2e-13. A displacement along the axis is totally symmetric, so within the pi irrep Schur's lemma forces G proportional to the identity -- no splitting. - perpendicular (x, y): every entry ~1e-14. For a linear molecule the linear term of a bending displacement vanishes by symmetry; that is exactly why the effect is second-order (Renner-Teller) rather than first-order Jahn-Teller. Both patterns are symmetry-forced and would be hard to produce by accident. The 1.2e-13 off-diagonal is the expected cancellation residual of (A2) at this force scale (~2e-12 relative) through a Z-vector solve. Note this also means Li2 cannot exercise a NONZERO off-diagonal; a nonlinear Jahn-Teller case (09_CH4/hf T2) is still owed. Full build: exit 0, no warnings from the changed files. The algebra underneath is covered by test_grad_matrix_degenerate (10 tests) from the first commit of this series, including the rotation-covariance check G' = U^T G U. Governance exception (rule 1, global dependency budget, net_delta = +1): - Reason: one `GlobalV::ofs_running` at the call site in `after_all_runners`. Any new log output costs at least one such reference, and the alternative -- not reporting the matrix -- would make the feature useless. - Scope: a single line in source_esolver; `Grad/` gains none, the driver takes `std::ofstream&` explicitly. The PARAM reference on that line replaces the one removed from the line it rewrites, so PARAM is net zero. - Risk: none functionally; it is the same log stream the surrounding code and `cal_force`'s existing `print_force` already use. - Cleanup: folds into any future change that threads a logger through the esolver layer instead of reaching for GlobalV. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat(lr-grad): store the degenerate gradient matrix, print its eigenvalues, add the average relax mode Three things the previous commit left undone. 1. The gradient matrix is kept instead of discarded. `cal_grad_matrix_degenerate` always returned it, but the only caller printed it and threw it away. It now lands in `multiplet_lvc_` as `MultipletLVC{ispin, states, omega0, omega_spread, g}`, rebuilt per spin channel on each call so a multiplet from the previous geometry cannot survive into the next one. It is held as one object rather than as loose force matrices because H(dR) = Omega_0 * I + sum_a dR_a G^(a) is the linear vibronic coupling Hamiltonian, and the relaxation and (later) non-adiabatic dynamics want the same parameter set. 2. Each coordinate's matrix is now followed by its eigenvalues. "Diagonalizing G" is not a well-defined operation -- G is 3N matrices of d x d and the branch slopes belong to a chosen direction u -- but the eigenvalues of a SINGLE G^(Aalpha) are meaningful: they are the branch slopes for displacing that one atom along that one axis. The printout says so, and says they must not be combined across axes, since G^(x), G^(y), G^(z) do not commute and picking one eigenvalue per axis describes no adiabatic state. 3. New `lr_relax_degen_mode` (String, default `state`) selects what a relaxation follows when `lr_target_state` sits inside a multiplet: - `state`: the historical per-state gradient. Kept as the default and as the only way to reproduce earlier results, although inside a multiplet it is the diagonal of G in whatever basis the eigensolver returned. - `average`: the multiplet average, gradient Tr(G)/d. This is a genuine smooth basis-independent surface and its gradient is totally symmetric, so following it stays on the symmetric configuration. A string rather than a boolean because the Jahn-Teller mode lands next and three meanings do not fit in two values -- and because flipping the meaning of `false` later would silently change existing results. `cal_energy()` switches to the same average through `target_omega_()`. This is the subtle half: `cg`/`bfgs`/`lbfgs` line-search on `cal_energy()`, so a gradient of one surface against the energy of another would not converge. At the degeneracy the two agree numerically, so the inconsistency would only appear once a displacement split the multiplet. `average` needs only the DIAGONAL of G, so it costs d solves, not d(d+1)/2. The off-diagonal part is what the Jahn-Teller direction needs. Also: `average_forces` now holds the one implementation of "average these force matrices", which had been written three times; `cal_force_Xz` gained the caller `cal_lr_force_relax_`, which is why the previous commit's plumbing split existed; and every `LR-Grad-formulas/*.md` pointer is removed from code comments (7 sites, 6 of them added by this series, 1 pre-existing in `hamilt_zeq_right.h`). Those notes are working documents that will not be committed, so the pointers would dangle. Each comment keeps the statement it was making; where the pointer carried the only hint of where a derivation lives, an in-tree pointer replaces it. Verification (build exit 0, no warnings from the changed files): - test_grad_matrix_degenerate: 10 tests, all passing. - Default path untouched: with `lr_grad_degen_thr` unset, the printed force block of fullwin/02_Li2/lda is byte-identical to the binary built before this series (diff -> IDENTICAL) and no degeneracy output is emitted. - fullwin/02_Li2/lda with the matrix enabled (OMP_NUM_THREADS=1, serial, systemd-run MemoryMax=8G): multiplet spread 2.2e-15 Ry, orthonormality 7.8e-16, and the new eigenvalues are symmetry-correct -- along the bond axis `eig 0.0519746 0.0519746`, exactly degenerate as Schur's lemma requires for a totally symmetric displacement inside the pi irrep, and ~1e-14 in both perpendicular directions, where a linear molecule has no first-order term. - INPUT: `abacus_std_para -h lr_relax_degen_mode` renders the description; `--check-input` rejects `average` without `lr_grad_degen_thr` with the intended message, rejects an out-of-range value, and accepts the legal combination. - Governance check: 1 ERROR, no warnings. NOT verified, and the first thing to check after the branch-tracking merge: `average` has not been executed end to end. It only runs under `calculation = relax`, and per the plan real relaxation cases wait for that merge. The code path compiles and its INPUT surface is validated, but an unexecuted path is exactly the failure mode this project has been bitten by before, so it should not be trusted until a relaxation has gone through it. Also not locally validated: `docs/parameters.yaml` could not be parsed here (no PyYAML in this environment). The new entry mirrors the neighbouring `lr_grad_degen_thr` entry's structure exactly; CI is what confirms it. Governance exception (rule 1, global dependency budget, net_delta = +1): - Reason: one `GlobalV::ofs_running`, passed into `cal_lr_force_relax_` at the relax call site so the note about switching surfaces reaches the log. Any new log output costs at least one such reference. - Scope: a single line in source_esolver. `Grad/` gains none -- both `cal_lr_force_relax_` and `cal_grad_matrix_degenerate` take `std::ofstream&` explicitly. - Risk: none functionally; it is the same stream the surrounding code and `cal_force`'s existing `print_force` already write to. - Cleanup: folds into any future change that threads a logger through the esolver layer instead of reaching for GlobalV. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat: state following by max overlap Cherry-picked from a50c6efab on lif-investigation, with the changes below. Original commit: a relaxation followed `lr_target_state` by INDEX every step, and an index names the n-th lowest root rather than a state. Where surfaces are close -- near-degenerate excitons being the normal case in a symmetric crystal -- the ordering swaps as the geometry moves, so a fixed index silently hops between diabatic states; the energy along the path is then not one smooth surface and its "gradient" is not conservative, which is what makes a CG relaxation stall with large erratic forces. `follow_target_state_` re-chooses the root at each step as argmax_j |<X_prev|X_j>| and `target_state_` replaces the input index everywhere downstream. Changes made while resolving: 1. Conflict in `runner`: the followed-root call `cal_force(target_is_, target_state_)` and this branch's `cal_lr_force_relax_` both replaced the same line. Resolved to `cal_lr_force_relax_`, which dispatches on `lr_relax_degen_mode` and already routes through `cal_force` for the single-state case, so the two features compose rather than exclude each other. 2. Integration, not just conflict resolution: `resolve_target_multiplet_` and `cal_lr_force_relax_` were keyed on `inp_->lr_target_state` and now key on `target_state_`. This matters -- `lr_target_state` is only the seed, so after any crossing the multiplet would have been resolved around a root the relaxation had already stopped following, and the average would have been taken over the wrong group. The degeneracy log line says "Followed state" rather than "Target state" for the same reason. 3. AGENTS.md rule 9: the overlap reduction used a direct `MPI_Allreduce(MPI_IN_PLACE, ...)` with a hand-computed `sizeof(T)/sizeof(double)` element count. Replaced with `Parallel_Reduce::reduce_all(ov.data(), nstates)`, the guarded wrapper, which also removes the `#ifdef __MPI` since it is a no-op in a serial build and removes the reinterpretation of a complex array as doubles. 4. `follow_target_state_` takes `std::ofstream&` instead of reaching for `GlobalV::ofs_running`, matching `cal_lr_force_relax_` and `cal_grad_matrix_degenerate` on this branch. This is what drops the global dependency budget from ERROR to WARNING: two references leave `Grad/` and one enters at the esolver call site, so the PR total is non-increasing and no rule-1 exception is needed for this commit. Verification (build exit 0, no warnings from the changed files): - Single-point path untouched: fullwin/02_Li2/lda printed force block is byte-identical to the binary built before this series (diff -> IDENTICAL), and no degeneracy output when `lr_grad_degen_thr` is unset. Expected, since the tracking only runs under `calculation = relax`, but it is the check that proves the cherry-pick did not disturb the shared code. - fullwin/02_Li2/lda with the gradient matrix enabled reproduces the previous commit's numbers exactly: spread 2.2e-15 Ry, orthonormality 7.8e-16. - test_grad_matrix_degenerate: 10 tests, all passing. - Governance check: 0 errors, 3 warnings. NOT verified: the tracking itself, and `lr_relax_degen_mode = average`, have not been executed. Both live on the `calculation = relax` path, which no case here exercises yet. This is now the one thing blocking confidence in the whole relaxation route, and it needs a real relaxation on a case with a crossing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat(lr-grad): find the Jahn-Teller direction from the degenerate gradient matrix The Jahn-Teller direction is a JOINT optimization over the displacement and the mixing inside the degenerate subspace -- they are determined together, which is why it cannot be done one Cartesian axis at a time: min_{|u|=1} lambda_min( sum_a u_a G^(a) ). Written that way it looks like a non-convex problem on the unit sphere in 3N dimensions. It is not. Since lambda_min(M) = min_{|v|=1} v^T M v, the two minimizations can be swapped, and the inner one over u has a closed form: min_{|u|=1} u . q(v) = -|q(v)|, q(v)_a = v^T G^(a) v, leaving min_{|u|=1} lambda_min(M(u)) = -max_{|v|=1} |q(v)|, u* = -q(v*)/|q(v*)|. So the search runs over the d-dimensional subspace (d = 2 or 3 in practice), not over 3N coordinates. It also reads plainly: q(v) is the gradient of the mixed state sum_k v_k |X_k>, so the Jahn-Teller direction is the steepest-descent direction of whichever state in the multiplet has the largest gradient. Solved by alternating exact minimization -- u <- -q(v)/|q(v)|, then v <- the lambda_min eigenvector of M(u) -- which decreases the joint objective monotonically. The objective is homogeneous of degree one and only piecewise smooth, so several DETERMINISTIC starting points are tried and the best kept; `restarts_agreeing` reports how many landed on it, which is the practical signal that the optimum is global. Deterministic because a relaxation must give the same answer twice. `split_symmetric_part` separates q(v) into the multiplet average Tr(G)/d, which is totally symmetric and common to every branch and only relaxes the geometry, and the remainder, which is the part that actually breaks the degeneracy. At a stationary point of the average surface -- where an `lr_relax_degen_mode = average` relaxation ends up -- the first term vanishes and the whole direction is Jahn-Teller, which is the intended two-stage workflow. The small symmetric eigensolver is written out as Jacobi rotations rather than taken from LAPACK, so this file stays free of the linear-algebra and parallel layers and remains unit-testable on its own; d is the dimension of an electronic multiplet and is never large. This is first order only. It gives the DIRECTION: the distortion amplitude needs the harmonic term too (Q ~ -g/k), and a norm that means anything physically should be taken in mass-weighted rather than plain Cartesian coordinates. Both caveats are in the header. Verification: test_grad_matrix_degenerate, 15 tests (5 new), all passing. Built and run standalone (g++ -std=c++14, gtest, openblas): [==========] 15 tests from 2 test suites ran. (96 ms total) [ PASSED ] 15 tests. The two substantive cases check the mathematics against independent brute force rather than checking the arithmetic against itself: - ReductionIdentityHolds: for d = 2 and d = 3, lambda_min along the returned direction is computed by a dense scan of the Rayleigh quotient and must equal -slope. A wrong reduction cannot pass this. - IsGlobalOverTheSubspace: |q(v)| is scanned over the whole subspace sphere and must not exceed the reported optimum, which is what distinguishes "converged to the global optimum" from "converged to a stationary point". Plus: the all-identity multiplet (no splitting available) must return the common gradient; scaling G scales the slope and leaves the direction fixed; and the symmetric/Jahn-Teller split reconstructs q(v) with a symmetric part independent of the branch. Not wired into the relaxation yet -- `lr_relax_degen_mode = jt` lands next. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat(lr-grad): lr_relax_degen_mode = jt, descending the Jahn-Teller branch Wires the search from the previous commit into the relaxation as a third value of `lr_relax_degen_mode`. `cal_jt_force_` assembles the off-diagonal part of the gradient matrix -- which `average` does not need, so this costs d(d-1)/2 further Z-vector solves per step on top of the d diagonal ones -- flattens it to one d x d block per nuclear coordinate, and hands it to `LR::find_jt_direction`. FORCES go in rather than gradients, so the direction that comes back already points downhill and no sign flip is needed anywhere. The force handed to the optimiser is that of the descending branch, q(v) with v the optimal mixing: the branch's own gradient, which is what a relaxation has to follow. The mode is self-limiting by construction -- after a step has split the multiplet, `resolve_target_multiplet_` finds no group and the ordinary single-state path takes over, which is also the point at which the overlap tracking from the previous commit becomes the thing that matters. The log reports, per step: the branch's force magnitude, its mixing coefficients, how many of the deterministic starting points agreed on the optimum, and the split of the force into the part common to the whole multiplet and the symmetry-breaking remainder. That split is the useful diagnostic -- at a stationary point of the average surface the common part vanishes and the whole force is Jahn-Teller, which makes `average` then `jt` the natural two-stage workflow. When the symmetry-breaking part vanishes the log says so explicitly and names the expected cause (a linear molecule, where the effect is second-order Renner-Teller and there is no first-order term at all). `lr_grad_degen_thr > 0` is now required by both non-default modes, not just `average`, since both need the grouping to know what the multiplet is. Verification (build exit 0, no warnings from the changed files): - INPUT surface: `--check-input` rejects `jt` without the threshold with "lr_relax_degen_mode=jt requires lr_grad_degen_thr > 0 to define the multiplet", rejects an out-of-range value with "lr_relax_degen_mode must be state, average or jt", accepts the legal combination, and `-h lr_relax_degen_mode` lists all three modes. - Single-point path untouched: fullwin/02_Li2/lda printed force block still byte-identical to the binary built before this series (diff -> IDENTICAL), no degeneracy output when disabled. - test_grad_matrix_degenerate: 15 tests, all passing. - Governance check: 0 errors, 1 warning. NOT verified: no relaxation has been run. `state` is exercised by every existing case, but `average` and `jt` both live on the `calculation = relax` path, and the Jahn-Teller code additionally needs a case with a genuine non-linear degeneracy (09_CH4/hf T2 is the obvious candidate -- Li2 is linear and by construction has no first-order Jahn-Teller term, so it can only exercise the warning branch). Until such a run exists these two modes are compiled and INPUT-validated but not demonstrated, and the algebra underneath is what carries the confidence: the reduction identity and the global optimality of the search are both checked against independent brute force in the unit tests. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * fix: adapt LR-gradient code to develop interface changes after rebase (part 2) Resolves the remaining build breakage left after the 41-commit `git pull --rebase origin develop`: Charge's copy constructor is now deleted (switch dm_to_charge to an out-parameter), Charge::allocate/cal_force_cc/ Structure_Factor::setup gained new required parameters, cal_dm_psi.h was renamed to dm_from_psi.h, and DensityMatrix::cal_dmr lost its default ik_in=-1 argument so call sites now pass it explicitly. Also reverts OperatorLRHxc/OperatorLREXX's constructor back to taking `const DensityMatrix<T,T>&` (not `unique_ptr<DensityMatrix<T,T>>&`): that unique_ptr signature was an incorrect choice made while resolving the e58237768 rebase conflict, and it broke ~20 pre-existing call sites in cal_edm_from_multipliers.h, hamilt_zeq_right/left.h, cal_multiplier_w_from_z.h, and lr_util.hpp that build these operators from a plain DensityMatrix reference, not a unique_ptr-owning member. Untrack LR-Grad-formulas/log/2026-08-最新解析&差分结果.md per request (kept on disk, no longer git-tracked). Verified: `make -j8 abacus_std_para` is clean (0 errors, 0 warnings). Smoke-tested 13_BeH/pbe (GS scf + LR-TDDFT gradient, 4 states, nspin=2 unrestricted): run completes and converged forces are finite with correct antisymmetry between the two atoms, consistent with the pre-rebase behavior. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * test(lr): add cal_force to the gamma-only LR-TDDFT integrate tests Enables `cal_force` on the four gamma-only LR-TDDFT tests in tests/08_RI (lr_tddft_pbe_gamma, lr_tddft_lda_gamma, lr_tddft_hf_gamma, lr_tddft_hf_ulr_gamma) to exercise the excited-state analytic gradient path. Multi-k cases (lr_tddft_multik) and the unsupported BSE kernel (bse_lr_multik) are left untouched, per request. catch_properties.sh gets a new `totallrforceref` extraction for LR runs: ESolver_LR prints its own "Forces (-gradients) of each excited state" table instead of the ground-state TOTAL-FORCE one, so it needs a separate parse (the existing TOTAL-FORCE force block is now gated to is_lr==0 to avoid silently recording a bogus zero for LR cases). Turning on cal_force exposed two real bugs in the gradient code, both fixed here: - `dft_functional == "default"` was compared as the literal string against `xc_kernel` to decide whether the ground-state potential could share its g^xc kernel with the LR one; "default" never resolves in PARAM, so this always disagreed even when the ground state and the LR kernel were the same functional, and `pot_hxc_gs` then threw the first time `cal_W_from_Z` asked it for g^xc. Fixed by reading the actual resolved functional from `ucell.atoms[0].ncpp.xc_func` when `dft_functional` is literally "default". - `PotGradXCLR::cal_v_eff(_openshell)` picked its LDA-vs-GGA branch from the global `XC_Functional::get_func_type()`, which reflects the ground state's functional, not the specific `KernelXC` it was constructed with. A cross-functional kernel (e.g. TDLDA on a PBE ground state, as in lr_tddft_lda_gamma) built an LDA KernelXC with an empty `drho_gs_` while the global func_type still said GGA, so the GGA branch indexed into an empty vector. Added `KernelXC::is_gga()` (true iff `drho_gs_` was actually filled) and branch on that instead. Also fixes two C++11-baseline violations (AGENTS.md rule 7) found while testing a CI config without LibRI/EXX (which otherwise bumps CMAKE_CXX_STANDARD to 14): a generic lambda parameter in lr_util_hcontainer.h::save_sparse, and another in zeq_solver.hpp's build_and_solve helper, extracted into a plain template function (build_and_solve_zeq) since C++11 lambdas can't take `auto` parameters. Verified: `make -j8 abacus_std_para` is clean (0 errors, 0 warnings). Ran all 4 updated test cases (GS scf + LR-TDDFT nscf, serial, OMP_NUM_THREADS=1) to completion; excitation energies match the prior result.ref to 5-6 decimals (ground-state physics unchanged), and the new totallrforceref values are finite with the expected Newton's-third-law antisymmetry between atoms. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * fix: guard MPI-only and EXX-only symbols for the no-ELPA/no-EXX CI builds Three pre-existing gaps in the LR-TDDFT gradient code, each only visible when a CI matrix leg compiles out MPI or LibRI/EXX: - `LR_Util::transpose_DMR` called `mattrans` (PBLAS pdtran_/pztranc_) unconditionally, but `mattrans`'s own declaration is correctly `#ifdef __MPI`-only (there is no serial equivalent for a 2D-block-cyclic dm(k) transpose). Guard the call sites too, with a clear runtime_error on the no-MPI side. - `cal_multiplier_w_from_z.h` and `cal_edm_from_multipliers.h` each had one unguarded call to an EXX-only symbol (`LR::gs_is_hybrid()` / `LR::exx_kernel_list()`, both declared under `#ifdef __EXX` in operator_lr_exx.h) alongside the corresponding `op_ht_exx`/`op_K_exx` operator object, which is itself correctly `#ifdef __EXX`-only. Both call sites are now guarded to match. - `operator_gxc_ulr.h` used `K_Vectors` without including its header; the type was only visible transitively through the `#ifdef __EXX` branch of another header's includes, so it disappeared once EXX was off. Added the direct include (source_cell/klist.h). Verified with `-fsyntax-only`: cal_edm_from_multipliers.cpp and esolver_lr_grad.cpp are clean under `-D__MPI` without `-D__EXX` (matching the "Build without ELPA"/"Build without MPI" CI legs' actual macro state -- ENABLE_LIBRI defaults OFF in all three CI variants, so the discriminator is really EXX, not MPI), and esolver_factory.cpp is clean with neither defined. Full `make -j8 abacus_std_para` under the normal MPI+EXX dev config is still clean (0 errors, 0 warnings). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * fix: more MPI-only guards, and wire the Grad sources into Makefile.Objects Continuing the no-MPI CI fixes: three more unconditional uses of MPI/PBLAS-only symbols, all with an existing serial fallback already used elsewhere in the same codebase, just not at these call sites: - `lr_density.hpp::cal_eh_density_single_state` called `cal_dm_diff_pblas` unconditionally; fall back to `cal_dm_diff_blas` without MPI, matching the `#ifdef __MPI` pattern already used in hamilt_zeq_right.h and cal_multiplier_w_from_z.h for the same two functions. - `zeq_solver.hpp::solve_Z_lapack` scattered the LAPACK-solved Z back into the 2D-block-cyclic layout via `scatter_full_to_2d` (PBLAS) unconditionally; the matching `gather_2d_to_full` step just above it already had a `#else std::copy(...)` fallback, so mirror it here. - `esolver_lr_lcao_tddft.cpp::before_scf`/`refresh_from_ks_` passed `Parallel_Orbitals::desc_wfc` (an `#ifdef __MPI`-only ScaLAPACK descriptor) to `fill_z_window_`; the callee already ignores that argument outside its own `#ifdef __MPI` branch, so guard the two call sites and pass `nullptr` without MPI. - `lr_util_hcontainer.h::build_dm_from_dmk` called the 3-argument `matsym(T*, int, const Parallel_2D&)` (PBLAS symmetrization) unconditionally; fall back to the always-available 2-argument serial overload. Also wires three Grad/ source files that were missing from Makefile.Objects' VPATH and object lists (the Makefile build never picked them up, even though the CMake build already builds them): `Grad/force/lr_force_test.cpp`, `Grad/degenerate/grad_matrix_degenerate.cpp` (unconditional, mirroring their CMakeLists.txt treatment), and `utils/lr_io_krlist.cpp` (BSE/LibRI-only, added under the existing `ifdef LIBRI_DIR` convention used elsewhere in the Makefile, since it unconditionally includes a LibRI-only header). The VPATH list was also missing `dm_band/` and every `Grad/*` subdirectory entirely -- added all of them. Verified: `-fsyntax-only` on esolver_lr_lcao_tddft.cpp and esolver_factory.cpp is clean with neither `__MPI` nor `__EXX` defined. Compiled `lr_force_test.cpp` and `grad_matrix_degenerate.cpp` for real with the Intel compiler (icpx via mpiicpc) through the Makefile build to confirm the VPATH fix resolves them. `make -j8 abacus_std_para` under the normal CMake/MPI/EXX dev config is still clean (0 errors, 0 warnings). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> * fix: make the gradient path build again under -DENABLE_MPI=OFF/Intel Makefile CI, and fix a parallel force bug Reproduced both failing CI legs locally against their exact flags (`cmake -DENABLE_MPI=OFF`; `make -f ../Makefile CXX=mpiicpx ELPA_LIB_DIR=... OPENMP=ON` with no LIBRI_DIR) instead of patching from partial pasted errors, and iterated until both link successfully. Build-fix root causes, in the order they surfaced: - `dm_band.cpp`'s header includes `<RI/global/Tensor.h>` unconditionally (not behind `#ifdef __EXX`), so it cannot be translated at all without LibRI. Its only consumer, `operator_lr_exx.cpp`, is itself entirely `#ifdef __EXX`. Moved `dm_band.cpp` into the `if(ENABLE_LIBRI)` block in CMakeLists.txt (next to `lr_io_krlist.cpp`, which was already there for the same reason) and the matching `ifdef LIBRI_DIR` block in Makefile.Objects. - `cal_edm_from_multipliers.cpp`/`.h` and `cal_multiplier_w_from_z.h` called several PBLAS routines (`pdgemm_`/`pzgemm_` via `ScalapackConnector::gemm`, `cal_dm_trans_pblas`, `cal_dm_diff_pblas`, the 3-argument `matsym`) and read `Parallel_2D::desc`/`blacs_ctxt` unconditionally. All of these are `#ifdef __MPI`-only (PBLAS needs an actual BLACS context), with serial fallbacks that already exist elsewhere in the codebase (`cal_dm_trans_blas`, `cal_dm_diff_blas`, the 2-argument `matsym`, raw `dgemm_`/`zgemm_` via `BlasConnector::gemm_cm`) but weren't wired in at these call sites. Added the missing `#ifdef __MPI ... #else ... #endif` branches throughout, mirroring the pattern already used elsewhere in the same files for the same functions. - Same pattern in `esolver_lr_grad.cpp`'s `cal_force_Xz`/`cal_force_openshell_Xz`: six more unguarded `cal_dm_trans_pblas`/`cal_dm_diff_pblas`/`matsym` call sites, fixed the same way. - `operator_lr_hxc.cpp::act`'s serial (`#else`) CXC branch referenced an undefined `psi_in_bfirst` identifier and dereferenced `this->psi_ks` (already a reference member, not a pointer) with `*psi_ks` -- this branch had apparently never compiled since CXC/CXC_o were added. Fixed to pass `psi_in` (the function's own parameter) and the spin-selected `psil_ks` (matching what the `#ifdef __MPI` branch passes), plus `this->nocc[sl]`/`this->nvirt[sl]` to match `CVCX_virt_blas`/`CVCX_occ_blas`'s scalar (not vector) parameters. - `lr_density.hpp`, `zeq_solver.hpp`'s `solve_Z_lapack`, and `esolver_lr_lcao_tddft.cpp`'s `fill_z_window_` call sites had the same unguarded-PBLAS/missing-serial-fallback pattern (`cal_dm_diff_pblas`/`scatter_full_to_2d`/`Parallel_Orbitals::desc_wfc`), fixed the same way (the callee for `desc_wfc` already ignores the argument outside its own `#ifdef __MPI` branch, so the fix there is just passing `nullptr` without MPI). - Two C++11-baseline violations (generic lambda parameters, which need C++14): one in `lr_util_hcontainer.h::save_sparse`, one in `zeq_solver.hpp`'s `build_and_solve` helper, extracted into a plain template function (`build_and_solve_zeq`). - Wired three Grad/ source files that were missing from Makefile.Objects' VPATH and object lists entirely (the Makefile build never picked them up, even though CMake already built them): `Grad/force/lr_force_test.cpp`, `Grad/degenerate/grad_matrix_degenerate.cpp`, and `utils/lr_io_krlist.cpp`. VPATH was also missing `dm_band/` and every `Grad/*` subdirectory entirely. Separately, a real parallel-correctness bug found while chasing a report that `17_LiH/39-e100`'s excited-state force differs between 1 and 4 MPI processes (Z-vector solve matches; the force calculated from it doesn't). Localized via the `test_force=1` per-term force breakdown: VL_dVL/EWALD/NLCC/NONLOCAL matched exactly between 1 and 4 processes (core ABACUS PW/Ewald/nonlocal machinery is fine under MPI), but KINETIC already differed by ~2x with no z-component, and everything after it (LOCAL-PP Pulay, HARTREE+XC Pulay/HF, HXC/GXC DMTRANS, the final force) diverged badly with a spurious z-component. Root cause: `pulay_force_hcontainer.h`'s 2-center-integration `cal_pulay_fs` overload (used for the KINETIC term, via `cal_hs_grad.h`'s `build_hcontainer_local_op`) and its two grid-integration overloads (`cal_pulay_fs`/`cal_pulay_fs_openshell`, which delegate to `ModuleGint::cal_gint_fvl`) all sum only the atom pairs/grid points this rank owns and return that partial sum directly, with no `Parallel_Reduce` before returning. Core ABACUS always follows the same grid-based `cal_pulay_fs`/`cal_gint_fvl` calls with an explicit `Parallel_Reduce::reduce_pool` at the call site (confirmed in `force_stress_lcao.cpp`); LR-Grad's own force code in `lr_force.cpp`/ `lr_force_test.cpp` never did, at any of the many call sites that use these functions directly. Correct by coincidence at nprocs=1 (every local share is the whole array there), silently wrong otherwise. Fixed by adding `Parallel_Reduce::reduce_all`/`reduce_pool` inside the three `pulay_force_hcontainer.h` overloads (fixing all of their callers at once) and after every direct `ModuleGint::cal_gint_fvl` / direct core-`cal_pulay_fs` call in `lr_force.cpp`/`lr_force_test.cpp`. Also fixed, while investigating (not the actual cause of the above, but a real bug): `cal_dm_gs()`/`test_force()`'s `edm_gs` construction read `psi_ks_all_` via `&this->paraMat_all_` unconditionally, but on the ks-lr path `psi_ks_all_` is aliased straight from the ground-state solver's own `psi` and is actually distributed per `this->ks_->pv` (`refresh_from_ks_`/`fill_z_window_` already make this same distinction for their `Cpxgemr2d` source descriptor; `cal_dm_gs()` never got it). Verified: - Both CI configs built to a linked execut…
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
ATTENTION: BSE module depends on LibRI#11 Implement CVC algorithm and enhance LR class abacusmodeling/LibRI#11
rpa_outdirandrpa_out_velfor Output files are not in the correct directories #7661Reminder
AGENTS.mdanddocs/developers_guide/agent_governance.md.Not yet. The BSE module depends on LibRPA output files for reading. I don't know whether we can upload these files.
source/changes.Linked Issue
Fix #7661
Unit Tests and/or Case Tests for my changes
What's changed?
Governance Notes