Skip to content

Fix(module_xc): compute real density laplacian for meta-GGA functionals#7286

Draft
dyzheng wants to merge 6 commits intodeepmodeling:LTSfrom
dyzheng:fix/scanl-laplacian
Draft

Fix(module_xc): compute real density laplacian for meta-GGA functionals#7286
dyzheng wants to merge 6 commits intodeepmodeling:LTSfrom
dyzheng:fix/scanl-laplacian

Conversation

@dyzheng
Copy link
Copy Markdown
Collaborator

@dyzheng dyzheng commented Apr 26, 2026

SCANL, R2SCANL, and REVSCAN require the density laplacian (∇²ρ) but v_xc_meta was passing |∇ρ|² instead. This caused ~150 kbar stress error for Na with SCANL vs SCAN.

  • Add cal_lapl() to compute ∇²ρ via reciprocal space
  • Fix v_xc_meta to pass real laplacian and handle vlapl output
  • Fix tau_xc and tau_xc_spin wrapper laplacian handling
  • Fix xc_functional_gradcorr to precompute ∇²ρ for stress/force paths

Reminder

  • Have you linked an issue with this pull request?
  • Have you added adequate unit tests and/or case tests for your pull request?
  • Have you noticed possible changes of behavior below or in the linked issue?
  • Have you explained the changes of codes in core modules of ESolver, HSolver, ElecState, Hamilt, Operator or Psi? (ignore if not applicable)

Linked Issue

Fix #...

Unit Tests and/or Case Tests for my changes

  • A unit test is added for each new feature or bug fix.

What's changed?

  • Example: My changes might affect the performance of the application under certain conditions, and I have tested the impact on various scenarios...

Any changes of core modules? (ignore if not applicable)

  • Example: I have added a new virtual function in the esolver base class in order to ...

@dyzheng dyzheng force-pushed the fix/scanl-laplacian branch from 32e29f4 to 0d7fbe2 Compare April 26, 2026 16:59
…eta-GGA functionals

SCANL, R2SCANL, and REVSCAN require the density laplacian (∇²ρ) but
v_xc_meta was passing |∇ρ|² instead. This caused ~150 kbar stress error
for Na with SCANL vs SCAN. Additionally, the vlapl stress contribution
was completely missing.

Laplacian computation method:
  ∇²ρ(r) = -∑_G |G|² ρ(G) e^{iGr}
  - ρ(r) is Fourier-transformed to ρ(G) via real2recip()
  - Each G-component is multiplied by -|G|² (i.e. -gg[ig] * tpiba²)
  - Inverse FFT via recip2real() yields ∇²ρ(r) on the real-space grid
  - Only the real part is retained (imaginary part should be ~0 for real ρ)
  This is exact within the plane-wave basis, not numerical differentiation.

Spin handling:
  - cal_lapl(nspin, nrxx, rho, tpiba2, chr) computes ∇²ρ for each spin
    channel: rho is interleaved [rho_up[0], rho_dw[0], rho_up[1], ...]
    and output lapl is [lapl_up[0], lapl_dw[0], lapl_up[1], ...]
  - nspin=1: single channel (spin-unpolarized), loop runs once for is=0
  - nspin=2: separate laplacian for spin-up and spin-down channels,
    each computed independently from its own real-space density
  - nspin=4 (noncollinear): blocked by set_xc_type for meta-GGA

Hessian computation for vlapl stress:
  - cal_rho_hessian(nspin, nrxx, rho, chr) computes ∂²ρ/∂r_α∂r_β
    via reciprocal space: FFT[ -G_α G_β ρ(G) ]
  - Returns 6 independent components per spin (xx, yy, zz, xy, yz, zx)
  - Stress contribution: σ_αβ^vlapl = -2 ∫ vlapl(r) ∂²ρ/∂r_α∂r_β dr
  - Applied in gradcorr is_stress path for both nspin=1 and nspin=2

Changes:
  - Add cal_lapl() and cal_rho_hessian() as reusable functions
  - Fix v_xc_meta to call cal_lapl() and pass real laplacian to LibXC;
    handle vlapl output by computing ∇²(vlapl) via same G-space method
  - Fix tau_xc/tau_xc_spin wrappers to return vlapl derivative
  - Fix gradcorr stress/force paths: replace duplicated inline Laplacian
    code with cal_lapl() calls; add Hessian computation and vlapl stress
    contribution for meta-GGA functionals
  - Update test_xc4.cpp and test CMakeLists for new tau_xc signature
@dyzheng dyzheng force-pushed the fix/scanl-laplacian branch from 0d7fbe2 to 812061e Compare April 27, 2026 15:09
@dyzheng dyzheng marked this pull request as draft April 28, 2026 00:11
@mohanchen mohanchen added Bugs Bugs that only solvable with sufficient knowledge of DFT Refactor Refactor ABACUS codes labels Apr 28, 2026
@jeanwsr
Copy link
Copy Markdown

jeanwsr commented Apr 30, 2026

should this be removed

bool not_supported_xc_with_laplacian(const std::string& xc_func_in)

dyzheng added 4 commits April 30, 2026 21:47
…arning

For meta-GGA functionals that depend on the Laplacian of the density
(e.g., SCANL), the vlapl contribution to V_xc (e2*∇²(vlapl)) is omitted
because it causes |G|² amplification in reciprocal space leading to
SCF divergence. vtxc also omits vlapl for consistency with V. The vlapl
energy is still included in etxc via exc, and the vlapl stress term
(2*vlapl*Hess*e2) is computed in gradcorr. A runtime warning is printed
when a Laplacian-dependent functional is detected, recommending r2SCAN
or SCAN instead. See docs/scanl_laplacian_implementation_report.md for
the full analysis including FD stress validation results.
Add ecutwfc 60/80/100 Ry PW and LCAO FD test results for scheme 3.
All schemes that omit vlapl from V show 5-10% stress error due to
etxc-vtxc mismatch. Update summary table accordingly.
…ential

Replace spectral Laplacian (|G|² amplification → SCF divergence) with
finite-difference Laplacian kernel in G-space. The FD kernel matches
the spectral kernel at low G but is bounded at high G, enabling SCF
convergence while preserving density-potential consistency.

Key changes:
- xc_functional_libxc_vxc.cpp: compute ∇²(vlapl·sgn) using FD kernel
  in G-space with metric tensor GGT for non-orthogonal cells
- xc_functional.cpp: register SCANL functional name (XC_MGGA_X_SCANL +
  XC_MGGA_C_SCANL)
- Update vtxc to include vlapl contribution for stress consistency

FD stress validation (Si2 FCC, Γ-point):
  SCAN:  FD/AB = 1.0000 (all ecutwfc)
  SCANL: FD/AB = 0.995  (ecutwfc=80, -0.51% error)
  vs scheme 3 (no vlapl in V): FD/AB = 0.906 (-9.5% error)
…st data

Add complete derivation of FD Laplacian kernel in G-space for general
(non-orthogonal) cells, including metric tensor cross terms. Document
numerical verification of gg_FD vs gg_spectral at multiple FFT grid
sizes, full FD stress validation results with energy data, and error
analysis.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Bugs Bugs that only solvable with sufficient knowledge of DFT Refactor Refactor ABACUS codes

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants