Repository navigation
Conversation
mohanchen
added
Features Needed
e5a8e55 to
cd8346b
Compare
|
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.
|
|
@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. |
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
Keep the existing dimensionless 1e-6 diagnostic threshold and explain its meaning in the warning. Rename state tracking sources and checker without changing numerical behavior. Test basenames use grp/chk abbreviations to retain the requested length limit; no INPUT or documentation behavior changes.
Separate orbital distributions and ground-state input, KS initialization and refresh, and dimension/kernel configuration into esolver_lr_orb.cpp, esolver_lr_init.cpp, and esolver_lr_cfg.cpp. Keep member definitions under source_esolver; update CMake and Make. Preserve numerical behavior and global dependency counts. No INPUT changes. Verified: full GNU MPI builds after each extraction; Intel Make object compilation; C++11 checks without MPI and ELPA; four-rank BeH2 average, CH4 JT, and unrestricted HF regressions.
|
@mohanchen The files named |
Apply the flattened UnitCell symmetry-analysis interface to esolver_lr_init.cpp while preserving the ESolver split and original feature SHAs. Verified full GNU MPI build, C++11 initialization checks without MPI and without ELPA, and four-rank two-step BeH2 state tracking. No INPUT behavior changes.
Use vector initialization for iat2it and iat2ia after the develop UnitCell refactor. Fix both the overlap block fixture and the cross-geometry AO overlap fixture. Verified both MPI-linked test targets, the four-rank state-overlap suite, and the serial state-overlap suite. No production behavior or reference changes.
Redistribute occupied and virtual MO-overlap subblocks to descriptors starting at (1,1) on the amplitude BLACS grid before PBLAS contractions. Preserve the existing block size and local distribution without gathering matrices. Remove the unit-block restriction while retaining dimension and layout checks. Regression tests cover block sizes 2, 4 and 32, misaligned occupied windows, multiple blocks, complex adjoints and empty ranks against independent contractions. Verified full GNU MPI build; 19 overlap tests in serial, MPI single-rank and four-rank runs; projection tests on a six-rank non-square grid; two-step BeH2 tracking with nb2d=32 on one rank and nb2d=2 on four ranks. No INPUT defaults or reference changes.
Redistribute selected old/new MO coefficient columns before AO-to-MO contraction, constructing only oo and vv overlaps once per spin. Reuse these aligned blocks for every tracked reference; projection now performs only the two matrix multiplications. Reuse the coefficient work buffer and avoid storing unused occupied-virtual overlaps. Preserve block sizes, numerical formulas and INPUT behavior. Verified full GNU MPI build; 21 overlap tests in serial and MPI single/four-rank runs; distributed tests on six ranks; two-step BeH2 tracking with nb2d=32 on one rank and nb2d=2 on four ranks. Added independent occupied/virtual construction checks including complex and misaligned large blocks. No reference changes or user-facing documentation changes required.
Merge cd2c7ce, including the upstream consistent AO radial spline values/derivatives correction (deepmodeling#8119). Resolve the unrestricted HF reference conflict by retaining the LR force check and adopting the verified new spectrum. Refresh force sums in the five Gamma LR cases; multi-k spectra and BeH2/CH4 state-following references remain unchanged. Keep existing thresholds and timing references. Verification: cmake --build build -j8 passed. Run all eight tests/08_RI/lr_* integration cases with OMP_NUM_THREADS=1, MKL_NUM_THREADS=1 and mpirun -np 4, using LibXC 5.1.7 with KXC as in CI; collect with catch_properties.sh. All 28 reference checks pass. Repeat unrestricted HF on four ranks (same force sum) and one rank (3.12e-7 deviation, below 1e-6). ABACUS v3.11.0-beta10.
|
If this PR is ready for review again, let me know. @maki49 |
Yes, ready. |
There was a problem hiding this comment.
Let's keep this file as simple as possible, we don't have exx, ofdft, sdft or rt-tddft information here, so I think the codes related to LR-TDDFT can be moved to somewhere else?
There was a problem hiding this comment.
how about 'exited_state_check' instead of state_check?
There was a problem hiding this comment.
this file name still contains 'root'?
There was a problem hiding this comment.
this file name contains 'root' ?
|
Did you introduce KXC and LXC in this PR? If so, please make sure to mention it. I noticed that compiling KXC increases the overall CI/CD build time. |
Yes, kxc is necessary for lda/pbe LR force and relax tests. lxc is not introduced. I've mentioned it in the description. |
# Conflicts: # source/source_lcao/module_lr/utils/test/CMakeLists.txt # source/source_lcao/module_ri/test/CMakeLists.txt # source/source_lcao/module_ri/unittests/test_exx_lr_ws.cpp
# Conflicts: # docs/advanced/input_files/input-main.md # docs/parameters.yaml
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
There was a problem hiding this comment.
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,47andoperator_lr_hxc.h:42-44. Use explicit arguments or a clear overload. - New comma declarations remain, e.g.
lr_grad_os.cpp:46andforce_funcs.h:29; declare one variable per line. zeq_solver.h:9propagates 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 --checkreports 33 whitespace errors; clean only the added/changed lines.- Non-blocking: correct the PR body's
lr_relax_degen_mode/lr_nvirttolr_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.
| // the Brillouin condition in every occupied-virtual rotation. X is zero-padded into it, so | ||
| // D^X and T come out bit-identical -- and Omega is untouched, so an existing finite-difference | ||
| // reference stays valid. See `fill_z_window_`. | ||
| const ct::Tensor Xz = this->pad_X_to_z_(ispin, ist_begin, nst); |
There was a problem hiding this comment.
[P1] Differentiate the actual truncated response problem.
Zero-padding X and widening only the virtual Z space does not produce the derivative of the selected excitation energy. For Gamma BeH2/RPA@LDA, fixed occupations, nbands=NLOCAL=17, nocc=3, nvirt=1 and LAPACK, H1-z is -1.7630594252 analytically versus -1.7791264016 eV/Angstrom by central differences (h=1e-3 Bohr); the complete nvirt=14 window agrees within 6e-6. With nocc=1, the triplet also differs: -1.4320207370 versus -1.4067535000. Include the omitted occupied and active/inactive orbital-response terms, or reject unsupported truncated force windows until implemented; do not silently change the response energy.
There was a problem hiding this comment.
Documented in 1ccd41e0d. Thank you. This is a known limitation of the current force formulation for truncated orbital windows, also seen in our earlier window-convergence tests. We retain the available-Z-space expansion and zero-padding of X as an approximation, without changing the selected truncated excitation energy. We do not claim that padding supplies all active/inactive orbital-response terms or makes the force the exact derivative of that truncated eigenproblem.
The INPUT reference now recommends all occupied orbitals, nbands=NLOCAL, and the full spin-resolved virtual window for reliable forces and relaxation. Setting nvirt=nbands-nocc alone is insufficient when nbands<NLOCAL. Initialization now warns once when the resolved occupied X window, virtual X window, or ground-state band space is incomplete, explicitly stating that the force may differ from finite differences of the truncated excitation energy and that window convergence must be checked before relaxation. The previous message suggesting that nvirt could remain unchanged for a converged gradient has been removed.
For context, our archived N2/HF comparisons report relative force discrepancies of 7.18% for (nvirt, nbands)=(8,26), 0.648% for (15,20), and 0.00% at the reported precision for the full (21,26) window. These comparisons vary X and/or the available KS band space and do not isolate a universal benefit from widening Z alone. In a separate BeH2/LDA Z-window comparison, the two reported singlet-group discrepancies decreased from +15.1%/+11.0% to -1.3%/-1.5% by padding Z window to nbands. These are archived finite-difference results, not newly rerun measurements or accuracy guarantees. An exact gradient for arbitrary truncated active spaces remains follow-up work; the truncated-force counterexample is not claimed fixed by this warning/documentation change.
| // On a non-positive-definite A only the rank owning the failing diagonal block returns | ||
| // (elpa_cholesky_template.F90) while the others wait in a broadcast, so with several | ||
| // ranks the check below is never reached and the run hangs. | ||
| elpa_cholesky(handle, A, &status); |
There was a problem hiding this comment.
[P1] Avoid the hanging ELPA failure path.
With ELPA 2022.11.001 and four ranks, the indefinite-matrix test hangs inside this call (15-second timeout, exit 124), while the ScaLAPACK Cholesky test returns the expected error. The pivot check cannot help because it is never reached; the current test skips this MPI case. Disable this backend/version combination or use a backend with collective failure handling.
There was a problem hiding this comment.
Fixed in 5b27615eb. Removed the ELPA-backed Z-vector solver option and its implementation. It differed from scalapack_chol only in the Cholesky factorization; both already used ScaLAPACK for the triangular solve. lr_grad_solver=elpa is now rejected during INPUT validation. The retained scalapack_chol reports a non-positive-definite factorization on every rank and recommends the LU-based lr_grad_solver=scalapack, without automatically changing the solver. Tests cover real and complex indefinite/singular matrices failing beyond the first distributed block, plus successful LU solutions for the nonsingular indefinite cases. The four-rank failure regression has a 60-second timeout. The tests pass with and without ELPA enabled; KS ELPA solvers are unchanged.
| // the ground-state SCF. `ESolver_KS::runner` is re-entrant: its `before_scf` rebuilds the | ||
| // neighbour lists, the grid tables and the Hamiltonian for the geometry of this step. | ||
| XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); | ||
| this->ks_->runner(ucell, istep); |
There was a problem hiding this comment.
[P2] Propagate the embedded KS convergence status.
The LR wrapper keeps its inherited conv_esolver=true after ks_->runner fails to converge. A scf_nmax=1 relaxation with a loose force threshold prints both 'SCF IS NOT CONVERGED' and 'Relaxation is converged!', then exits 0. Propagate the KS status and reject the optimization step/success result when the ground state is unconverged.
There was a problem hiding this comment.
Fixed in 324cacf9c. The LR wrapper now propagates the embedded KS convergence status and rejects an excited-state relaxation step before LR initialization or force evaluation when the ground-state SCF is unconverged.
| const int ekb_off = this->openshell ? 0 : this->target_is_ * this->nstates; | ||
| double sum = 0.0; | ||
| for (const int ist : this->target_group_) { sum += this->pelec->ekb.c[ekb_off + ist]; } | ||
| return sum / static_cast<double>(this->target_group_.size()); |
There was a problem hiding this comment.
[P2] Return the same JT branch energy and force.
This averages every nonempty group, including JT mode, while cal_jt_force_ returns a selected branch force. For H(q)=diag(1-q,1+q), q=1e-6 and threshold=1e-4, the actual JT helper returns force 1 but the returned mean energy has derivative 0. Use a consistent branch energy/force for line searches; equality at the exact degeneracy does not fix the finite grouping interval.
There was a problem hiding this comment.
Documented in 0bf7c327c. The JT construction is intended for a numerically degenerate subspace. At an exact degeneracy, the selected mixture gives the first-order branch slope through degenerate perturbation theory; no additional derivative of its mixing coefficients is required, and splitting after displacement is expected. Your example does identify a limitation of the finite grouping tolerance: while split roots remain in one lr_degen_thr group, the code reports the group mean energy but returns the selected JT force, so an energy-based line search need not be consistent. The existing single-root path returns that root's own energy and force once threshold-based splitting leaves the tracked group with one member. I have documented this approximation explicitly under lr_degen_thr, including that a small energy spread does not bound the force error, and clarified the switching condition under lr_degen_mode. The JT algorithm is unchanged; this documentation update does not claim to restore energy–force consistency throughout the finite grouping interval.
| const bool& openshell = openshell_; | ||
| const std::vector<std::vector<ModuleBase::Vector3<double>>>& drho_gs = drho_gs_; | ||
| /// Whether THIS kernel (built for its own functional name, which may differ from the | ||
| /// ground state's `dft_functional` in a cross-functional run such as TDLDA@PBE) needs |
There was a problem hiding this comment.
[P2] Reject file-based kernels before requesting unsupported forces.
The lr_init_xc_kernel=file branch in xc_kernel.cpp:39-56 returns with fxc populated but no gxc, so force evaluation reaches this throw. The same input succeeds with cal_force=0; cal_force=1 finishes SCF/LR and then aborts with 'singlet g^xc coefficients were not built' (exit 134). Reject this combination during input/initialization, or define and implement the external kernel's geometry derivatives.
There was a problem hiding this comment.
Fixed in 324cacf9c. INPUT validation now rejects lr_init_xc_kernel=file for LR forces and relaxation, before SCF.
| AddTest( | ||
| TARGET MODULE_LR_state_group | ||
| LIBS base parameter ${math_libs} container device | ||
| SOURCES test_state_grp.cpp ../state_group.cpp |
There was a problem hiding this comment.
Governance: align new test filenames with their modules.
Rule 11 requires test_.cpp: rename test_state_grp.cpp to test_state_group.cpp, test_zeqlin.cpp to test_zeqlin_solv.cpp, and dm_trans/test/dm_diff_tst.cpp to test_dm_diff.cpp. Update their build registrations; historical untouched tests need no renaming.
There was a problem hiding this comment.
Fixed in 324cacf9c. Renamed the three new test files to match the source modules exactly and updated their CMake registrations.
| :PotLRBase(rho_basis, LR_Util::kernel_nspin(), chg_gs.nrxx, ucell.tpiba), | ||
| xc_kernel_(xc_kernel), spin_type_(st), | ||
| pot_hartree_(LR_Util::make_unique<elecstate::PotHartree>(&rho_basis)), | ||
| xc_kernel_components_(make_kernel(xc_kernel, rho_basis, ucell, chg_gs, pgrid, (st == SpinType::S2_updown), gxc_spin, lr_init_xc_kernel)), |
There was a problem hiding this comment.
Governance: move nontrivial construction work into the body.
The initializer lists call kernel_nspin(), make_kernel(), allocation helpers and XC queries, contrary to rule 13. Move computed preparation into inspectable constructor-body steps, or record a narrow exception where the member design requires restructuring.
There was a problem hiding this comment.
Fixed in 34ee45493. Kernel preparation is now an explicit build_xc_kernel free-function call. PotHxcLR only binds an already-built shared kernel; it no longer computes XC coefficients in its initializer list or reads the current global XC type or spin count. Each kernel retains its own functional metadata, so binding an earlier LDA kernel after preparing PBE, or vice versa, selects the correct integration path. Existing same-functional sharing is preserved, and LR/BSE callers have been migrated. Four focused tests passed, along with the GNU/MPI build and four-rank before/after comparisons: all 162 signed force components and 28 transition data files were unchanged.
| #endif | ||
|
|
||
| // use lapack to solve the linear equation | ||
| ModuleBase::timer::start("Z_vector", "lapack_solver"); |
There was a problem hiding this comment.
Governance: time complete functions.
Here and at the Z_vector_R block below, start/end wrap isolated calls inside another function. Rule 12 requires the enclosing function's entry/exit and name; move these scopes to the corresponding functions.
There was a problem hiding this comment.
Fixed in 324cacf9c. The timer scopes now cover the complete solve_Z_lapack and build_and_solve_zeq functions, using their function names. The latter explicitly covers both RHS construction and the solve.
| const ModulePW::PW_Basis& rho_basis_; | ||
| const int nspin_ = 1; | ||
| const int nrxx_ = 1; | ||
| std::unique_ptr<elecstate::PotHartree> pot_hartree_; |
There was a problem hiding this comment.
Non-blocking: remove the unused Hartree object.
pot_hartree_ is allocated by both constructors but never read; add_v_hartree uses its own implementation. Remove the member, allocations and unnecessary include. This cleanup can follow separately.
There was a problem hiding this comment.
Fixed in 324cacf9c. Removed the unused Hartree member, its allocations and the corresponding includes. The dedicated Hartree application remains unchanged.
| #include <gtest/gtest.h> | ||
| #include "../zeq_solver.hpp" | ||
| #include "../cal_edm.h" | ||
| #include "../cal_edm.h" |
There was a problem hiding this comment.
Non-blocking: remove the duplicate include.
cal_edm.h is included twice. Drop the duplicate and check whether this linear-solver test needs the EDM header at all; this can be a follow-up cleanup.
There was a problem hiding this comment.
Fixed in 324cacf9c. Removed both EDM includes; this linear-solver test does not require the EDM header.
Add a protected HamiltLR storage-only constructor and delegate the regular Casida constructor to it. Closed-shell Z operators initialize their own chains directly, keeping the shared transition density instead of replacing it after deleting an unused Casida chain. RI benchmark DMR initialization and the generic Operator implementation remain unchanged. Verification (OMP_NUM_THREADS=1, GNU/MPI, LibXC 5.2.3-kxc): - cmake --build build -j4 --target abacus_std_para MODULE_LR_zeqlin_solv - ctest --test-dir build --output-on-failure -R '^MODULE_LR_zeqlin_solv$' - Four-rank before/after runs of lr_tddft_lda_gamma and lr_tddft_hf_gamma: all 36 signed force components per case identical. - Four-rank lr_tddft_multik spectra and the two Gamma cases: five transition data files per case byte-identical before and after. - git diff --check; staged agent governance check: no errors. The focused existing runtime cases exercise both Z chains and ordinary Casida construction. No user documentation or INPUT changes are needed: this only changes internal construction and preserves numerical behavior.
Move kernel creation from PotHxcLR to the build_xc_kernel free function. Potentials now bind shared, already-evaluated coefficients and obtain the functional type from kernel-owned metadata rather than the current global XC state. Pass nspin explicitly and reject missing or mismatched kernels. Initialize spin metadata on every KernelXC path, including file and HF/RPA kernels, while preserving same-functional coefficient sharing. Migrate LR and BSE initialization to the explicit preparation interface. Kernel evaluation still uses the existing global XC setup internally; this change removes the ordering dependency from subsequent potential binding. Verification (GNU/MPI, LibXC 5.2.3-kxc, OMP/BLAS threads=1): - cmake --build build -j3 --target abacus_std_para MODULE_LR_pot_hxc_lrtd (targets built separately): passed. - ctest --test-dir build --output-on-failure -R '^MODULE_LR_pot_hxc_lrtd$': passed; four tests cover interleaved LDA/PBE binding, pointer sharing, file/HF/RPA spin metadata, and invalid kernel bindings. - The same four tests also passed using isolated staged sources, without the unrelated, uncommitted HSE component experiment or an ELPA link. - Four-MPI before/after LDA (nspin=1 and 2), PBE, HF, multi-k spectra, and LDA@PBE runs: all 162 signed force components identical and all 28 transition data files byte-identical. - git diff --cached --check and staged agent governance: no errors. The memory header is required by the shared_ptr factory declaration. No user documentation or INPUT update is needed for this internal refactor.
Remove lr_grad_solver=elpa, its implementation, dispatch, and test linkage. Keep ScaLAPACK Cholesky for positive-definite systems; report factorization failure on all ranks and recommend the LU-based scalapack option instead. Stop function timers before throwing so repeated calls report the actual matrix error rather than a timer-state error. KS ELPA solvers are unchanged. Parameter metadata documents the supported choices and failure behavior; user-facing YAML/Markdown are generated from this metadata. Verification (OMP/MKL/OpenBLAS threads=1): - GNU/MPI build: abacus_std_para, MODULE_LR_zeqlin_solv, MODULE_IO_input_test_para passed. - ctest --test-dir build --output-on-failure -R '^MODULE_LR_zeqlin_(solv|failure_mpi)$': 2/2 passed, including the four-rank indefinite/singular regression with a 60-second timeout. - All 14 solver tests passed on four ranks, including real/complex LU solutions for nonsingular indefinite matrices rejected by Cholesky. - Without ELPA: solver target rebuilt and the same two CTests passed. - InputParaTest.RemovedElpaZSolverIsRejected passed; --check-input rejects elpa with exit 1 and an explicit removal diagnostic. - CLI help, generated parameter YAML/Markdown, staged governance and git diff --cached --check verified; no governance errors.
Retain available-Z-space padding while diagnosing incomplete occupied, virtual, or KS band windows once during initialization. State explicitly that padding does not restore every orbital-response term and the force may differ from finite differences of the truncated excitation energy. Recommend full spin-resolved windows and force/window convergence checks. Remove the previous message suggesting nvirt could remain unchanged for a converged gradient. Update the generated INPUT reference metadata for nocc/nvirt; generated YAML/Markdown are not committed. Verification: - GNU/MPI executable built successfully with this change. - cmake --build build -j3 --target MODULE_LR_lr_util_pp and OMP_NUM_THREADS=1 ctest --test-dir build --output-on-failure -R '^MODULE_LR_lr_util_pp$': passed again before this commit. - Earlier four-rank, two-step BeH2 checks of this implementation passed for full, virtual-truncated, occupied-truncated, and band-truncated windows; full windows emit no warning, each truncated case warns once. - CLI help and generated parameter YAML/Markdown contain the caveats. - Staged governance: no errors; git diff --cached --check passed. This documents an approximation, not an exact truncated-window gradient fix. Archived FD measurements were not rerun for this documentation change.
Explain under lr_degen_thr that a threshold-defined group uses the mean excitation energy with a selected JT force. At exact degeneracy this gives the first-order branch slope without a mixing-coefficient response; with finite splitting inside the tolerance, energy-based line searches need not be consistent. A small spread does not bound the force error. Clarify that single-root energy and force resume only when threshold-based splitting leaves the tracked group with one member. The algorithm is unchanged and finite-interval energy-force consistency is not claimed. Verification: abacus_std_para rebuilt successfully; -h lr_degen_thr, --generate-parameters-yaml and docs/generate_input_main.py output all contain the limitation. git diff --cached --check and staged governance passed. User docs are generated from the updated parameter metadata.
Keep the historical absolute-force totals and add a separate signed comparison keyed by spin channel, state, atom, and Cartesian component. Require exact entry sets, reject malformed/duplicate/nonfinite data, and use each case's existing force_threshold (1e-6 eV/Angstrom). Only the separate LR force tables are read; GS and excited-state total forces and the existing state-tracking validators are unaffected. Add component references for the five existing LDA/PBE/HF force cases from pre-change logs. Their absolute sums match the existing scalar references; no existing numerical values or tolerances are changed. These are behavior baselines, not independent physical validation. No finite-difference calculations are added to CI. Verification (OMP/MKL/OpenBLAS threads=1): - python3 -m unittest discover -s tests/integrate/validation_tools -p test_lr_force_check.py -v: 7 tests passed, including sign reversal, atom/axis/state permutations, missing data, duplicates/nonfinite values, tolerance handling and shell collector failure propagation. - ctest --test-dir build --output-on-failure -R '^lr_force_collector$': passed; the lightweight unit tests are registered with CTest. - Real Autotest.sh runs of all five cases, sequentially with 1 and 4 MPI processes: all 153 signed components and existing properties passed, without warnings. No FD calculations were performed. - bash -n props_lr.sh, staged governance and git diff --cached --check: passed with no findings.
The serial read-item tests still expected the old lr_grad_solver option list and build-dependent ELPA diagnostic. Update these assertions for the removed option and require ELPA rejection independently of build flags. Keep valid solver and serial ScaLAPACK rejection checks unchanged. Verification (OMP/MKL/OpenBLAS threads=1): - Reproduced InputTest.Item_test2 locally with the same two assertion failures reported by MODULE_IO_read_item_serial in CI. - cmake --build build -j3 --target MODULE_IO_read_item_serial: passed. - ctest --test-dir build --output-on-failure -R '^MODULE_IO_read_item_serial$': passed; all seven tests passed. - Staged governance and git diff --cached --check passed. Only test expectations change; no production behavior changes.
Features
Refactors
Other Behavior changes
1e-6to1e-14, and the gradient cutoff from1e-10to1e-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 approximately1.29e-3, above the old cutoff; its four excitation energies agree at 17-digit output precision between the two settings.nspin=1singlet response uses the factor-two Hartree/XC contribution, while ground-stateS1_gsretains weight one. Existing single-spin excitation results can therefore change. BSE grid initialization usesS1_gsbecauseHamiltBSEalready applies the singlet factor; the existing RI BSE path is unaffected by this grid correction.1e-6when evaluating XC derivatives. This can change HSE06 response energies and gradients near vanishing density gradients; it does not apply this clamp to ordinary PBE.kxc, rather than allowing the last functional to overwrite earlier contributions. This affects the XC contribution to analytic excited-state forces.lr_nvirtcontinues to select the response basis. The gradient can depend on the available ground-state band window even when the response basis is held fixed.test_forceruns.