Repository navigation
feat(dfpt): analytic libxc XC kernel with the i(G+q) carrier term (exact GGA response at q != 0) - #8163
feat(dfpt): analytic libxc XC kernel with the i(G+q) carrier term (exact GGA response at q != 0)#8163Zanthoxylum wants to merge 3 commits into
Conversation
mohanchen
added
DFPT
|
@ABACUS-PR-BOT could you review the PR? |
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
There was a problem hiding this comment.
Reviewed statically at 7ac3caf8990cf93977584f2fa3b891658b870b93 with independent verification. No builds or runtime/numerical tests were run; CI diagnostics were unavailable. Low/signed-density regularization and per-Libxc FXC availability remain unresolved.
Non-blocking: please add the exact validation commands, Libxc/MPI configuration, and observed pass/fail/not-run outcomes to the PR body. The current test descriptions and frequency observations do not establish which checks actually ran.
| { | ||
| const double rho = rho0_r[ir]; | ||
| const double sig = sigma[ir]; | ||
| if (std::abs(rho) <= epsr || sig <= epsg) |
There was a problem hiding this comment.
[P1] This guard suppresses k_rr_ along with the gradient fields when sigma <= epsg, even at a dense symmetry point or in the uniform-density limit. The local XC derivative is still nonzero there; unlike QE's dgradcorr, this provider has no separately accumulated local term. Please retain that response according to the ground-state density policy and add a zero-gradient regression; the current GGA test skips those points.
| for (size_t ifunc = 0; ifunc < func_id_.size(); ++ifunc) | ||
| { | ||
| xc_func_type func; | ||
| if (xc_func_init(&func, func_id_[ifunc], XC_UNPOLARIZED) != 0) |
There was a problem hiding this comment.
[P2] xc_func_init uses Libxc defaults, while the ground-state init_func applies xc_exch_ext/xc_corr_ext and finite-temperature LDA parameters. Allowed configured LDA/GGA calculations therefore differentiate a different functional in DFPT. Please pass or reuse the same effective parameters and test a nondefault configuration.
| // e_sigma_sigma term the factor 4 (QE xc_wrapper_d_gga | ||
| // convention; no rho renormalization, the potential-form | ||
| // kernel cancels it) | ||
| k_rr_[ir] += ha_to_ry * v2rho2; |
There was a problem hiding this comment.
[P2] The unconditional sum loses the ground-state component weights. For example, the accepted type-2 functional BLYP_LR subtracts its ITYH and LYPR components, while these fields add both with weight +1. Please carry the same weights into every derivative field; preserving external parameters alone will not fix this.
| diff = std::max(diff, std::abs(dq[ir] - d0[ir])); | ||
| scale = std::max(scale, std::abs(d0[ir])); | ||
| } | ||
| EXPECT_GT(diff / scale, 1.0e-3); |
There was a problem hiding this comment.
Non-blocking: this finite-q assertion only establishes that the output changes, so a wrong q sign or scale could also pass. Please add a quantitative oracle, such as a commensurate-supercell potential finite difference, that checks the complex response. Keep the useful LDA q-independence test.
| AddTest( | ||
| TARGET MODULE_DFPT_xc_analytic_serial | ||
| LIBS parameter dfpt_planewave_serial device base symmetry Libxc::xc | ||
| SOURCES dfpt_xc_analytic_serial_test.cpp |
There was a problem hiding this comment.
Non-blocking: please rename the new source to test_dfpt_xc_analytic.cpp and update this entry, following the project's naming rule for new unit tests. The surrounding historical filenames are exempt.
|
|
||
| #include "dfpt_rho.h" | ||
|
|
||
| #include "source_base/matrix3.h" |
There was a problem hiding this comment.
Non-blocking: this header never uses Matrix3, and it already gets PW_Basis through dfpt_rho.h. Please remove the redundant includes while retaining the complete base and Vector3/vector definitions. This small cleanup can be handled separately.
| // exists yet, so reject them explicitly instead of producing silently | ||
| // wrong frequencies. | ||
| const int xc_func_type = XC_Functional::get_func_type(); | ||
| if (xc_func_type != 1 && xc_func_type != 2) |
There was a problem hiding this comment.
Non-blocking: please describe the LDA/GGA-only DFPT scope and the no-Libxc finite-q GGA limitation in the C++ Input_Item help. The YAML and INPUT-reference files are ignored generated artifacts on this branch; validate the output through the existing --generate-parameters-yaml/Sphinx workflow, or give a concrete no-documentation-update rationale.
| // kernel linearization point: the ground-state XC density, i.e. valence | ||
| // plus the nonlinear core correction density when the pseudopotentials | ||
| // carry one (PotXC_FDM adds chg->rho_core the same way); the first-order | ||
| // density itself stays valence-only (a frozen core does not respond) |
There was a problem hiding this comment.
Non-blocking: frozen core means its atomic shape is fixed; its density still translates with a displaced ion. Please qualify this as valence-only electronic screening at fixed ions, separate from the ionic NLCC derivative. That wider ionic-NLCC omission predates this PR, so this is only a wording/scope clarification.
Review finding of deepmodeling#8164 (P1): the sigma <= epsg guard zeroed k_rr_ together with the gradient fields, so flat-gradient regions (vacuum, uniform-density interstitials) lost the local (rho-rho) XC response entirely -- unlike QE, whose dv_of_drho accumulates the local term separately from the guarded dgradcorr gradient terms. The local kernel is finite in the zero-gradient limit (the LDA-like term of the GGA functional), so only the density guard applies to it; the gradient chain fields keep the epsg guard. Adds the zero-gradient regression test the review asked for: a uniform density (sigma == 0 everywhere) must respond with dvxc = k_rr A, finite and nonzero, exact to 1e-10.
…n the kernel Review findings of deepmodeling#8164 (two P2): xc_func_init used libxc defaults, so configured LDA/GGA runs (xc_exch_ext / xc_corr_ext fractions, finite-T LDA) differentiated a different functional in DFPT than the ground state used; and the unconditional +1 sum lost the signed combination factors (BLYP_LR subtracts its short-range ITYH/LYPR components). - XC_Functional_Libxc::ext_params_for: expose the per-component parameter lists init_func would apply (external override in-built), so the kernel instantiates exactly the GS effective functional without duplicating the parameter tables - XC_First_Order_Analytic gains a configured ctor taking those lists and the signed weights; the esolver factory feeds both from XC_Functional (new XC_Functional::get_weight accessor) - per-component family dispatch: mixed LDA+GGA combinations evaluate LDA members through the LDA entry points - tests: KSDT at nondefault theta matches the manually parameterized libxc second derivative; signed weights ({-2,+3}) scale every field; neutral configuration reproduces the plain ctor bit-for-bit
Review finding of deepmodeling#8164 (non-blocking): the trace/summary lines mutated std::cout flags and precision, and the changes leaked into the opt-in DFPT_DEBUG snapshots of the next iteration (setprecision(1) truncating them). Both blocks now format into a local ostringstream, leaving the shared stream state untouched.
Review response — fixes for the findings at 8d93aff (2026-10-19)
[P1] local GGA response below the gradient cutoff — fixed in 549a6d1Finding: the Fix: the local kernel is finite in the zero-gradient limit (the LDA-like term of the [P2] configured Libxc parameters ignored — fixed in db71dcaFinding: Fix: new [P2] signed component weights lost — fixed in db71dcaFinding: the unconditional +1 sum lost the ground-state combination factors — e.g. Fix: Follow-ups included in the same commits
|
This branch was successfully deployed
3 participants
Summary
PW-DFPT phonons evaluated with GGA functionals carried a systematic error at every q != 0, while LDA agreed with reference codes everywhere and all functionals were correct at q = Gamma. This PR replaces the finite-difference first-order XC path with an analytic kernel built from libxc second derivatives whenever ABACUS is compiled with libxc (
ENABLE_LIBXC), and adds an explicit warning on the FDM fallback path.Part 1 — LDA vs GGA: what the bug was
For a phonon perturbation at wavevector q the first-order density is a carrier wave
drho(r) = 2 Re[A(r) e^{i q r}]with a carrier-free amplitude A stored on the unit-cell FFT grid. The GGA first-order XC response isThe third term is the gradient carried by the phase factor
e^{i q r}itself; it has no counterpart in grid gradients of A. The previous FDM path finite-differenced the XC potential on +/- A around the ground state, so its gradients only ever sawgrad A— the carrier term was structurally missing. The term vanishes exactly in the two limits that matched the observed failure pattern:f_grad = 0— LDA phonons were correct at every q;q = 0— optical modes were correct for every functional.In between (GGA at q != 0) the omission interacted with the variational stiffness of the Sternheimer solve and the large-cancellation structure of polar materials, amplifying a percent-level kernel deficit into mode errors of tens of percent. QE evaluates the same response with the gradient operator acting on the physical field through
i(G + q)(fft_qgradientin gradutils.f90,dgradcorrin LR_Modules), which includes the carrier exactly.Part 2 — Why libxc for the kernel
The carrier term is analytic: restoring it requires the exact GGA kernel tensor (
f_rho_rho,f_rho_sigma,f_sigma_sigma,f_sigma) evaluated on the ground-state density. Hand-deriving and maintaining analytic second derivatives per functional is error-prone and does not scale beyond PBE, while libxc (xc_gga_fxc/xc_lda_fxc, requirement >= 5.1.7) provides them for every functional it ships. Design decisions:use_libxc = false), so ground-state baselines do not move.xc_wrapper_d_gga):k_rr = v2rho2,k_sr = 2 v2rhosigma,k_ss = 4 v2sigma2,k_s = 2 vsigma; LDA reduces tok_rr = v2rho2(libxcvrhois already the potential form, so no additional rho renormalization appears — the potential-form kernels cancel it).e2 = 2(Hartree to Ry) andtpiba(PW_Basisgcaris dimensionless — the physical wavevector istpiba * gcar, the same conventionDFPT_Rho::v_hartree_quses with itstpiba2divisor).rho0 = rho + rho_core, exactly whatPotXC_FDMfeeds the internal functional, seepot_xc_fdm.cpp), while the first-order density itself stays valence-only (a frozen core does not respond). Feeding the valence density alone moves nlcc systems off their ground-state point: bn_cubic Gamma optical modes dropped ~15 cm^-1 versus both the FDM path and QE. Non-nlcc systems are unaffected.XC_First_Order::set_context(q_cart)— a const per-q context hook on the existing provider base; the solver injects it inassemble_v_scbeforeapply. Without libxc the FDM adapter remains the default and now warns that the GGA response at q != 0 lacks the carrier term (LDA stays exact).MODULE_DFPT_xc_analytic_serial— LDA kernel vs finite difference ofv_xc; the full GGA chain at Gamma cross-checked against a field-perturbation finite difference (an algebra-level independent path); carrier-term presence at finite q for GGA and exact q-independence for LDA. Existing DFPT serial suites are unchanged.