Repository navigation
Feature: LR-TDDFT analytical forces (gradients) by Z-vector method and excited-state geometry relaxation - #8069
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
|
487: [ RUN ] lr_tddft_hf_ulr_gamma 487: [ RUN ] lr_relax_lda_beh2 487: [ RUN ] lr_relax_lda_ch4 487: [ RUN ] lr_tddft_hf_gamma The example above takes a significant amount of testing time in CI/CD. Could you further reduce its testing time? If there are other examples where the testing time could also be reduced, that would be much appreciated as well. |
|
I hope the above two will be the last two comments of this PR |
|
Let me see if the bot have further comments. |
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
There was a problem hiding this comment.
Reviewed ac4bf0d7b2007160dcd4e5238493a58892e8d6b5 against 27bbd6cfa6ec40ef4c5699b0c46f2e7412b840b7 with independent static source/algebra reviews. New findings are inline; the signed-HSE issue remains in its existing discussion and is not repeated.
Existing governance blocker: the new-default-argument feedback remains unresolved. Ten new defaults across six existing interfaces (OperatorLRDiag, the three AO-to-MO functions, OperatorLRHxc, and OperatorLREXX) require explicit call-site updates or a maintainer-approved scoped exception under the governance matrix.
Non-blocking metadata: please use the registered names lr_degen_mode / nvirt in the PR description and add the exact validation commands and outcomes. Explain that YAML/Markdown INPUT references are build-generated from the updated Input_Item metadata; no generated files need to be committed.
Verification this round: trusted governance checker exit 0 (200 warnings, no mechanical errors); the manually verified default-argument rule above still applies. git diff --check failed with 32 trailing-whitespace diagnostics; all four changed Python files parsed. No builds, test suites, numerical/MPI jobs or CI verification were performed; CI diagnostics were unavailable. Full force accuracy and the recorded HSE/window/JT/occupation/parallel/BSE limitations remain unverified.
Prepare unrestricted Z dimensions at the call site, align the EXX workspace test name, and pass former default arguments explicitly throughout LR/BSE. Document parameter roles and reject external kernels/densities and fixed electric fields for LR forces or relaxation in INPUT metadata and validation.
Declare each force matrix on its own line without changing initialization. GNU/MPI main-program build, no-MPI LR module build and both focused Z-solver CTests pass. No behavior or INPUT changes; no new tests or generated documentation are needed for this declaration-only cleanup.
Fixed in 7b9c4eb. Removed all ten newly added default arguments from the six existing interfaces, together with preceding defaults where required by C++ parameter ordering. All call sites now pass the values explicitly, including LR/BSE and serial/MPI AO-to-MO tests. Added concise comments for scaling, accumulation, contraction selection, benchmark probe indices and operator type. The explicit values preserve the previous numerical behavior; no replacement defaulted overloads were introduced. The comma declarations are also fixed in 2e13b88. |
Keep save_HR/save_DMR as declaration-only adapters and put their common-writer calls in lr_sparse.cpp. Synchronize structural R coordinates before conversion, use explicit spin suffixes, and keep the LR gradient call site high-level. Diagnostic headers follow the common CSR format; force formulas are unchanged. Wire the new implementation through CMake and GNU Makefile. Pyabacus has no independent LR source list. Move the disjoint-R regression into the LR DM test and exercise both wrapper layers with two spin channels, a zero rank and a globally zero R block. Common IO tests retain no dependency on LR. GNU/MPI main build and serial LR build pass. Serial and four-rank wrapper tests pass; actual one/four-rank LDA test_force CSR files are byte-identical to the pre-wrapper version. Unit-test global references decrease by three; the retained DRANK assignment adapts the common writer's existing output-rank interface.
Temporarily select the resolved GS functional while Potential initializes its cached Hxc values, then restore all incoming XC selection and runtime settings with an exception-safe scoped adapter. Potential, PotXC and XC evaluation still use process-wide state; removing that dependency requires a separate three-layer refactor. Handle HF's empty functional-ID list safely. Four focused scope tests cover restoration, exception unwinding, nesting and LibXC backend selection. MPI and no-MPI builds/tests pass. HF GS reconstruction is invariant under RPA/LDA/PBE LR kernel selection (previous deviations 3.35/3.48 eV/Angstrom); one/four-rank probes use bz=1 to avoid an independent empty-grid diagnostic limitation. All six existing LR cases pass at one and four ranks with unchanged references.
Adapt the response-density screening fix from GPU branch commit e82898e. Retain density blocks in Casida/Z exchange and both sides of the LR exchange-force contractions; leave the ground-state SCF threshold unchanged and document its scope. Always evaluate the true residual of the returned CG vector, including zero/small RHS exits. Preserve the normal iteration sequence and the existing 1e-6 absolute tolerance. Add false-convergence and zero-RHS operator checks, also exercised on four MPI ranks. Update only the two HF force baselines after checking operator linearity, CG versus LU, MPI agreement, and full-window H2 finite differences. Excitation-energy and LDA/PBE references are unchanged. Temporary operator probes and extra-precision FD output are excluded from this commit. Validation: GNU/MPI executable build; no-MPI/no-ELPA lr/esolver builds; four focused CTests pass; six original LR cases pass on one/four ranks. Full-window H2 force FD relative errors are below 3e-7 at 0.001 Angstrom.
Add abs_gauge=none for LR-TDDFT and BSE while preserving excitation-energy checks, charge output, forces and relaxation. Initialize required projector metadata independently of the optical gauge. Test the energy collector, including stale optical files. Use coarser tracking-regression inputs and a direct solver for the small unrestricted HF Z equation. Keep the closed-shell HF CG coverage and all existing references and tolerances. Validation: GNU MPI executable and INPUT test builds passed; serial no-MPI/no-ELPA lr and esolver module builds passed. Eight integration cases passed with one and four MPI ranks. BSE TDA/full none-versus-velocity energies, standalone LR, INPUT/collector tests, CLI checks and generated parameter documentation passed. Governance: moved BSE rank guards retain the existing global dependency count. The print helper includes provide its inline numerical output. INPUT documentation is generated from the updated metadata.
@mohanchen Fixed in 05e814650. I reduced the cost of the four regression cases by:
A local comparison using four MPI ranks, one thread per rank, and identical CPU pinning gives the following wall-clock medians of two runs per configuration, including MPI startup, using the same executable:
The sum of these medians decreased from 65.73 to 26.46 s (59.7%). These are local measurements rather than predicted CI timings; the closed-shell HF improvement is small. EXX thresholds, numerical references, and comparison tolerances are unchanged. Other LR cases continue to exercise the optical paths. |
Include parallel_reduce.h directly for the collective used by the empty-block regression, fixing the undeclared Parallel_Reduce error in CUDA CI builds. No user-facing behavior or documentation changes. Validation: MODULE_LR_zeqlin_solv built with GNU/MPI; ctest -R ^MODULE_LR_zeqlin_(solv|failure_mpi)$ passed both tests, including the four-rank regression. Full CUDA build was not tested locally because nvcc is unavailable.
|
There are some further places that can be improved. However, since the PR is already large. I will accept it in its current form. Please keep in mind that further refactors are still needed, I do hope that you can help in future @maki49 |
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.