Skip to content

Fix LRI force with a spectrally truncated Coulomb metric - #8224

Open
maki49 wants to merge 5 commits into
deepmodeling:developfrom
maki49:fix-lrc-force-dC
Open

maki49 wants to merge 5 commits into
deepmodeling:developfrom
maki49:fix-lrc-force-dC

Conversation

@maki49

@maki49 maki49 commented Oct 27, 2026 •

Copy link
Copy Markdown
Collaborator

Correct the derivative of LRI coefficients with a spectrally truncated Coulomb metric.

In N₂/LC-PBE, the local auxiliary metric has a condition number of approximately $1.4\times10^{12}$, producing unphysical forces without regularization. Spectral truncation (setting non-zero exx_cs_inv_thr) stabilizes the inverse, but the previous derivative still treated the truncated pseudoinverse as an ordinary inverse, omitting rotations between retained and discarded auxiliary subspaces.

For a Coulomb matrix $V=U\Lambda U^\dagger$, define $f_i=1/\lambda_i$ for $\lambda_i>\tau\lambda_{\max}$ and zero otherwise. With $F(V)=U\mathrm{diag}(f_i)U^\dagger$ and $C=F(V)A$ (where $A$ is the LRI-fitted 3-center integration $(P\vert\phi_1\phi_2)$ and $F(V)=V^{-1}$ when $V$ is not spectrally truncated), use

$$ C^x=F(V)A^x+ U\left[\Gamma\circ(U^\dagger V^xU)\right]U^\dagger A, \qquad \Gamma_{ij}=\begin{cases} (f_i-f_j)/(\lambda_i-\lambda_j),&i\ne j,\\ f'(\lambda_i),&i=j. \end{cases} $$

to replace the dCs formula with a ordinary Coulomb inverse: $C^x=V^{-1}(A^x-V^xC)$.
See the full derivation in LC-PBE-force-anomaly-2026-10-27.md

The cross-subspace terms track how the retained fitting space rotates with atomic motion. Reuse one eigendecomposition for all three Cartesian derivatives and retain the existing untruncated path. The derivative assumes no eigenvalue crosses the cutoff.

N₂ validation (forces in eV/Å):

Quantity Before fix After fix Shared FD
Ground-state total force -16.19786482 -16.19792180 -16.19794562
Lowest singlet excitation force +14.56047593 +14.56051600 +14.56050000
Lowest triplet excitation force +14.52237009 +14.52241051 +14.52240000

Both versions use exx_cs_inv_thr=1e-5 and the same converged N₂ reference; the FD values are shared. Printed ground-state and excitation energies are identical. Both force implementations agree with FD at approximately 1e-4 eV/Å. This N₂ force check is a regression test, not evidence that the derivative correction alone removes the unregularized force failure. Separately, at the same cutoff, the offline local-coefficient derivative (dCs) error decreases from 36.13% with the old formula to 0.00365% with the full derivative.

Multi-k regression (Si DZP 6 bohr orbitals): distorted two-atom Si, Γ-centered 2×2×2, LC-PBE, exx_cs_inv_thr=1e-5 (ground-state forces only). The controlled multi-k comparison below reduces the absolute force error from 0.03990 to 0.000451 eV/Å (approximately 88-fold). A finite residual remains.

Force (eV/Å) Before fix After fix Shared FD
Si atom 2, z -0.46568500 -0.42533206 -0.42578338

Validation: full ABACUS build and four real/complex finite-difference unit tests passed. No new INPUT parameters; existing parameter metadata now correctly describes Cholesky inversion and the relative cutoff, from which user documentation is generated.

Copilot AI balanced review requested due to automatic review settings October 27, 2026 18:32

Copilot AI left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Copilot was unable to review this pull request because the user who requested the review has insufficient remaining quota.

An organization owner can choose the repository organization as the billing source for reviews requested by members with a paid Copilot license. Learn about Copilot code review billing. Review the Copilot code review billing settings.

@maki49 maki49 changed the title Fix LRI force derivatives with a spectrally truncated Coulomb metric Fix LRI force with a spectrally truncated Coulomb metric Oct 27, 2026

@ABACUS-PR-BOT ABACUS-PR-BOT left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Static review of feff584; the suggestions below are nonblocking. No builds or runtime tests were run, and CI verification was unavailable.

Non-blocking: The reported validation is useful; please add the exact configure/build/CTest and N2/Si run commands, plus input locations, versions and checksums (or a shareable fixture). The posted evidence points to local checks/... paths, so reviewers cannot reproduce the candidate configuration and finite differences from the PR alone.

item.category = "Exact Exchange (LCAO)";
item.type = "Real";
item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using LU decomposition. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues smaller than exx_cs_inv_thr are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5.";
item.description = "By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using Cholesky factorization. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues smaller than exx_cs_inv_thr times the largest eigenvalue are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5.";

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Non-blocking: Please carry the fixed-membership limitation into this generated INPUT description. The derivative assumes modes remain on the same side of the relative cutoff; crossings can make the truncated model nonsmooth during relaxation or MD. The PR already explains this, and a short note here would make it visible to users; follow-up is fine.

Comment thread source/source_lcao/module_ri/unittests/test_invmat_ri.cpp
pack.input(responses[axis]);
block_responses[axis] = pack.output(sizes, sizes);
}
return std::make_pair(inverse.output(sizes, sizes), block_responses);

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Non-blocking governance: Root rule 14 asks for computed call arguments to be named first. Please store these inverse blocks before make_pair and apply the same pattern to the new displaced-metric inputs and reverse DPcal_dV call. This is an inspectability cleanup.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

try not using .hpp file, use .h and .cpp files instead

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

same as above, try not using hpp file

@mohanchen mohanchen added Refactor Refactor ABACUS codes EXX/lr-TDDFT/RPA/BSE/GW Related to beyond DFT methods labels Oct 29, 2026
…iles

Replace invmat_ri.hpp and lri_cv.hpp with compiled implementations and
explicit double/complex<double> instantiations. Move the new cal_I_dI
adapter into lri_cv_tools.cpp, keeping the existing generic tools intact.
Make Ewald/EXX tool dependencies explicit and update CMake, the focused
unit-test target, and the legacy Makefile object list.

No numerical formulas or INPUT behavior change; no user documentation
update is needed. The direct tools includes are required for the existing
Ewald/EXX template implementation calls.

Validation:
- cmake --build build --target MODULE_RI_invmat_ri ri -j 4
  (RI initially exposed missing direct includes; passed after correction)
- cmake --build build --target abacus_std_para -j 4: passed
- OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_RI_invmat_ri$'
  --output-on-failure: passed (8 cases)
- git diff --cached --check: passed
- agent_governance_check.py --staged: no errors; advisory warnings only
After the implementation-header split, the standalone workspace target
still compiled exx_workspace.cpp without lri_cv.cpp. Its instantiated
workspace factories require real and complex LRI_CV destructors, which
made the CI all-target build fail at link time although the main binary
and inverse-derivative test built successfully.

Include the compiled LRI implementation and use function/data sections
with linker garbage collection, as in the existing focused inverse test,
to keep unrelated integral assembly outside this geometry-sharing test.
No numerical formulas, INPUT behavior, or user documentation change.

Validation:
- Reproduced the missing LRI_CV destructors with
  cmake --build build --target MODULE_RI_exx_workspace -j 4 before fix.
- Rebuilt all six MODULE_RI targets successfully after fix.
- OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_RI_'
  --output-on-failure: 6/6 passed.
- cmake --build build --target abacus_std_para -j 4: passed.
- agent_governance_check.py --staged and git diff --cached --check: passed.

@ABACUS-PR-BOT ABACUS-PR-BOT left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Static incremental review at 8c51f32f4281; no builds or runtime tests were run.

RI::Tensor<Tdata> A;
};

#include "invmat_ri.hpp"

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Non-blocking: with the inverse definitions now in invmat_ri.cpp, move the remaining abfs_construct_pca.h include there too. Only the implementation traits use PCA::tensor_syev; the public declarations do not need the orbital and UnitCell headers it brings in.

The CUDA CI build succeeded, but scf_out_hk_spin2 aborted in host
ScaLAPACK with info=1/ifail=5 (eigenvector nonconvergence). The CUDA job
has LibRI disabled, so this is independent of the RI implementation.

Let this matrix-output fixture use the existing device-default solver:
cuSOLVER on GPU and the available backend on CPU. Preserve all energy,
H(k)/S(k), H(R)/S(R) references and tolerances. ScaLAPACK remains covered
by its dedicated solver tests and other integration fixtures. This does
not change production solver policy or claim to repair ScaLAPACK itself.

Also apply the review suggestion to move abfs_construct_pca.h from the
inverse public header into its cpp file. Give the existing tools template
bodies their own required UnitCell include rather than inheriting it via
PCA/orbital headers. A complete UnitCell type is required for member access
in those bodies, so a forward declaration alone is insufficient.

Validation:
- Built RI targets and abacus_std_para successfully.
- OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_RI_'
  --output-on-failure: 6/6 targets passed.
- Two-rank CPU output fixture: original ScaLAPACK and ELPA each pass all
  11 original comparisons; modified device-default input also passes all
  11 comparisons with the rebuilt binary.
- Governance check: no errors; include/doc advisory warnings explained.
- git diff --cached --check: passed.
- No local CUDA compiler/GPU; cuSOLVER verification requires CUDA CI.

No public INPUT default changes. The fixture README documents the
regression setup; no user manual update is needed.

This branch has not been deployed

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

Labels

EXX/lr-TDDFT/RPA/BSE/GW Related to beyond DFT methods Refactor Refactor ABACUS codes

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants