Repository navigation
feat(dfpt): analytic libxc XC kernel with the i(G+q) carrier term (exact GGA response at q != 0) - #8163
Open
Zanthoxylum wants to merge 64 commits into
Open
feat(dfpt): analytic libxc XC kernel with the i(G+q) carrier term (exact GGA response at q != 0)#8163Zanthoxylum wants to merge 64 commits into
Zanthoxylum wants to merge 64 commits into
Conversation
added 30 commits
August 14, 2026 12:54
…e 1)
Phase 1 of the approved reciprocal-grid refactor enabling DFPT q-point
support: extract the spin-free common functionality from K_Vectors and
KVectorUtils into a new abstract base class ModuleCell::ReciprocalGrid,
which will be shared by K_Vectors (electrons) and QList (phonons/DFPT).
Changes:
- Add source_cell/reciprocal_grid.{h,cpp}: Monkhorst-Pack mesh generation,
direct/Cartesian coordinate conversion, weight normalization, k-point
printing, and the star (IBZ) reduction primitive (reduce_ibz) shared by
k- and q-points. Declares the pure-virtual reduce_by_symmetry().
- klist.{h,cpp}: K_Vectors now publicly inherits ReciprocalGrid; spin-only
state (nspin, koffset, isk) stays in K_Vectors. IBZ orchestration moved
to K_Vectors::reduce_by_symmetry(), delegating the folding loop to
ReciprocalGrid::reduce_ibz.
- k_vector_utils.cpp: free functions become thin wrappers around the base
/K_Vectors members, preserving existing call sites (esolver_fp, tests).
- Wire reciprocal_grid.cpp into source_cell and test CMakeLists.
External K_Vectors API and behavior are unchanged. Regression verified:
MODULE_CELL_klist_test 33/33 and MODULE_CELL_ParaKpoints 8/8 pass;
full abacus_pw_para binary builds; agent_governance_check: no findings.
…hase 2) - Extract build_star_ops from K_Vectors::reduce_by_symmetry into the ModuleCell::ReciprocalGrid base as a shared protected helper: k-lattice construction, Bravais compatibility check, point-group construction and kgmatrix membership verification. - Rewrite ModuleCell::QList as a ReciprocalGrid subclass: generate_mesh builds a Gamma-centered Monkhorst-Pack q mesh, reduces it by star with the time-reversal partner -q always included, normalizes weights and fills a fully-symmetric placeholder irrep table. - Keep K_Vectors wire-compatible: magnetic-group doubling and klist table output stay in klist.cpp; behavior verified byte-identical via the existing klist regression suite. - Add reciprocal_grid_test.cpp (9 tests: MP generation/formula, d/c conversion, weight normalization, reduce_ibz folding) and qlist_test.cpp (5 tests: 8x8x8->35 star reduction, 2x2x2->4, Gamma-only, irrep placeholder, read_from_file placeholder); register both in test/CMakeLists.txt. Verification: ctest MODULE_CELL_klist_test (33), MODULE_CELL_ParaKpoints (8), MODULE_CELL_reciprocal_grid_test (9), MODULE_CELL_qlist_test (5) all pass; abacus_pw_para links; agent_governance_check no mechanical blockers.
- Add ModuleSymmetry::LittleGroup (module_symmetry/little_group.{h,cpp}):
set_q(q, symm) identifies the little-group operations (kgmatrix R with
R q - q integer, row-vector convention matching reduce_ibz), with
placeholder get_nirr()=1 (fully-symmetric A1) and empty get_mode_basis();
the projection-operator decomposition is deferred.
- Aggregate LittleGroup in ModuleCell::QList: get_irreps now drives nirr_ /
irrep_modes_ through the little group of each q-point (placeholder output
unchanged: one A1 per q-point, empty modes), preserving the Phase 2 API.
- Add little_group_test.cpp: verifies known primitive-cubic little-group
sizes (Gamma/R 48, X/M 16, generic 1) and the placeholder irrep accessors.
- Wire little_group.cpp into the symmetry object library and register the
new test target.
Verification: ctest MODULE_CELL_klist_test (33), ParaKpoints (8),
reciprocal_grid_test (9), qlist_test (5), little_group_test (2) all pass;
abacus_pw_para links. Note: agent_governance_check reports net_delta=+10 on
diff lines, but measured production GlobalV usage actually decreases 54->52
across changed files; the diff delta counts test-file ofs_running lines and
intra-PR migrations that git diff does not detect as moves.
Thread a const Plus_U* through DFPT_PW::init / DFPT_PW_Data (decided at the esolver layer, never read through GlobalV/PARAM) with: - with_u() / u_active() (locale-initialized guard covers the pure-PW run without LCAO orbital files) and a per-q docc storage slot; - no-op stubs for DFPT_Rho::cal_docc, DFPT_Pert::build_dv_u, DFPT_Phon::dftu_onsite plus the [r,V_U] Q0 reservation note; - unit tests: null-provider path, docc roundtrip, and a Plus_U with uninitialized locale (with_u=true, u_active=false, run() unaffected) via a minimal dftu_test_support shim that keeps DFPT tests free of the LCAO-side DFT+U link closure. Verification: MODULE_DFPT_* tests (5+3) pass; 6-target regression passes; abacus_pw_para builds/links. Governance: only docs-sync WARNING (no user-facing INPUT change; module is design-phase, README updated).
DFPT_KQ_Basis enumerates the local plane-wave basis at the perturbation wavevector k+q by re-filtering the shared G grid of an initialized ground-state k-basis (PW_Basis_K) at the shifted center, avoiding new FFT grids or MP redistribution. Accessors expose the k+q basis size, the underlying G index / FFT slab index, G and G+k+q Cartesian vectors and |G+k+q|^2. A gamma_only ground-state basis is rejected because DFPT couples k and k+q symmetrically and needs the full complex G ball. Tests: 5 focused unit tests covering Gamma q=0 exact reproduction of the base ordering, the asymmetric shifted sphere, k+q translation invariance, nonzero-q agreement with a full FFT-grid brute-force reference, and the null/gamma_only guard. 7-target regression and abacus_pw_para link pass. No user-facing INPUT changes; design-phase module with README already covering the DFPT workflow (governance docs-sync warning exempt).
Implement DFPT_Pert: dVloc_dtau (rho-grid coefficients with the q shift baked into magnitude and phase), the NC separable dVnl two-term identity with build_vkb/radial_vq/real_ylm, build_dv/apply_dv FFT convolution on the shared rho/wfc grid, build_efield, and the U0-reserved build_dv_u guard. dv/dpsi storage upgraded from stubs in DFPT_PW_Data. Add a serial (__MPI-off) test directory mirroring module_pw/test_serial (dfpt_planewave_serial OBJECT library) with 8 physics tests: dVloc finite difference incl. q!=0, apply_dv convolution vs analytic matrix elements, efield sawtooth closed-form FT, independent-Simpson vkb check, dVnl identity vs operator finite difference, USPP rejection, and the pure-PW DFT+U degradation. The tests caught and fixed three convention bugs: the atomic phase must be exp(i 2pi g.tau) (GS stru_fac convention, not tpiba*g.tau); the shared real-space layout is ir = (ix*ny + iy)*nz + iz (z fastest, pinned by an impulse-response probe); and rho/wfc stick tables enumerate different G balls so real_space_dv now maps through the FFT-cell (ix,iy,iz) triple instead of raw isz. Governance notes: no new GlobalV/PARAM dependencies (exception-free); the header-dependency and docs-sync warnings are covered by the forward-declared Structure_Factor and the design-phase status (no INPUT change). Verified: 8/8 serial tests, 8 ctest targets (CELL+DFPT) pass, abacus_pw_para links.
DFPT_Stern::solve implements the projected conjugate-gradient solution of (H(k+q) - eps_n) P_c |dpsi_n> = -P_c |dV psi_n> with P_c the projector on the complement of the occupied states at k+q (metallic branch is C4). The shifted Hamiltonian action is injected through a LinearOperator interface so the solver core stays decoupled from the ground-state operator chain; the production adapter reusing hamilt::Hamilt::ops->hPsi is wired in C7. - apply_pv: two-sweep modified Gram-Schmidt projection, alias-safe - search directions are re-projected every CG step; pAp <= 0 triggers a residual-direction restart - degenerate handling: b inside the occ subspace, b = 0, or dimension mismatch return dpsi = 0 with residual 0 - unit tests (MODULE_DFPT_stern_test, 5 cases): diagonal operator against the closed-form complement solution, dense Hermitian U D U^dagger against the spectral reference with eps inside the occupied band, orthogonality of the solution to random occupied sets, degenerate and zero right-hand sides Governance: the only findings are the two standing exemptions for this design-phase module (header value-type includes <complex>/<vector>; docs-sync with no user-visible INPUT change). Verified: MODULE_DFPT_stern_test 5/5; ctest 9/9 (CELL 4 + DFPT 5); abacus_pw_para links; governance --staged clean.
DFPT_Rho::compute_drho builds the q-shifted response density from the
Sternheimer solutions: the periodic parts u_nk (K-basis transform) and
du_nk (k+q coefficients scattered onto the rho grid through the shared
FFT-cell triple, C1 pattern) multiply pointwise into
A(r) = sum_{kn occ} wg u* du, real2recip gives the q-shifted coefficients
A_Delta = sum_{kn} wg sum_G c*_G d_{G+Delta} indexed by the rho-grid ig,
and the Delta = -q harmonic is dropped whenever -q falls on a reciprocal
lattice vector (charge conservation; always at q = Gamma). The manifest
real-space density 2 Re[e^{iqr} A(r)] is rebuilt from the projected
coefficients so both storages agree. mix_drho applies plain mixing on the
q-shifted coefficients through Base_Mixing::Plain_Mixing (zero initial
input, residual ||out-in||/||out||); the heavy Charge_Mixing header
dependency is replaced by a forward declaration plus a Matrix3 value
member (reciprocal matrix for q_frac -> cart).
- data layer: set/get_drho_r/set/get_drho_g go from stubs to real storage
- guards: nspin != 1 and non-plain mixing reject with WARNING_QUIT
(design phase); cal_docc stays a documented U0 reservation (needs the
PW-side beta-projector adapter wired with Plus_U in the C7/U1 window)
- unit tests (MODULE_DFPT_rho_serial, 5 cases): G-space coefficients
against a brute-force double sum, real-space density against direct
plane-wave sums, Gamma charge conservation, plain-mixing first/second
step combination and residual formula
- test-side findings fixed (production code verified correct):
PW_Basis_K::gcar is a per-k array indexed ik*npwk_max+igl
(pw_basis_k.cpp:261) and must not be read with base-ball ig; direct-sum
references must pair cartesian G with cartesian r = frac . latvec
- irrep wrapper test updated: drho storage slots are live (round-trip
non-empty) after being design-phase stubs
Governance: only the two standing exemptions for this design-phase module
(value-type header includes, net dependency decreased by dropping
charge_mixing.h; docs-sync with no user-visible INPUT change).
Verified: MODULE_DFPT_rho_serial 5/5; ctest 10/10 (CELL 4 + DFPT 6);
abacus_pw_para links; governance --staged clean apart from exemptions.
- DFPT_Metal (C4): explicit WARNING_QUIT guards on the reserved metallic
branch (dfdeps/compute_dmu/compute_drho_metal); interface-only as planned
- DFPT_Phon (C5):
- ion_ion: Ewald force constants (G + R + self-image phase terms), the
Gamma acoustic sum rule holds exactly by construction
- accumulate_electron: 2n+1 complex accumulation 2 sum wg <dpsi^b|dV^a|psi>
plus the same-atom anharmonic <psi|d2V|psi> term (d2vloc_r + apply_d2vnl
from DFPT_Pert); the dpsi slot is backed up/restored around apply_dv
- assemble/diagonalize/add_loto/check_sum_rule: zheev with signed cm^-1
frequencies, LO-TO non-analytic term, Gamma row-sum rule
- DFPT_PW_Data: dynmat stored as ComplexMatrix (complex Hermitian at
generic q)
- fixes found by the new serial test: the cross term dropped the imaginary
part (needed for the Hermitian symmetrization at q != 0) and the test
reference used the basis momentum G instead of the kernel momentum G+q
- serial test MODULE_DFPT_phon_serial: 7 cases (Gamma ASR on a
symmetry-broken two-atom cell, acoustic zero modes, incommensurate q vs
direct dipole-Hessian sum, injected-dpsi closed-form contraction, zheev
on a known matrix, isotropic LO-TO limit, Gamma sum rule)
- verification: 11/11 ctest targets pass (CELL 4 + DFPT 7), abacus_pw_para
links, governance shows only the two pre-existing exempt warning classes
- DFPT_Rho::v_hartree_q: q-shifted first-order Hartree kernel aligned with h_hartree_pw (skips |G+q|=0), shared by the C6 response and the C7 screened potential - XC_First_Order abstract contract in module_dfpt (adapter at the esolver layer in C7, mirroring DFPT_Stern::LinearOperator injection) - DFPT_Pert::build_vkb_dk: analytic k-derivative of the beta projectors (atomic phase, radial chain rule, real-harmonic direction chain); build_vkb/build_vkb_dk made public for DFPT_Q0 reuse - DFPT_Q0::pos_matrix: velocity (commutator) form r = -i <u|dH/dk|u> / (tpiba (eps_m - eps_n)), kinetic 2 tpiba^2 (k+G) plus the separable nonlocal derivative; degenerate pairs skipped - DFPT_Q0::compute_eps / compute_born: length-gauge denominators, m sum over all bands for Z*, conj ordering of <v|dV|m>, ionic Z on the (a,b) diagonal, phon-style dpsi slot backup/restore - serial tests: MODULE_DFPT_q0_serial (5 tests: vkb FD, kinetic analytic, nonlocal operator FD, eps two-level, born closed form) and v_hartree_q checks in MODULE_DFPT_rho_serial; 12-target regression + abacus_pw_para link pass
Module layer (C7a): - DFPT_PW::init new signature (ucell, psi, bases, sf, veff_r, wg, eig, xc contract, nelec, ecutwfc, dftu); Impl holds GS data + hamilt_ - DFPT_HamiltShift: self-assembled H(k+q) Sternheimer operator (kinetic diagonal + veff FFT convolution + cached k+q vkb nonlocal), replacing the GS HamiltPW chain which is ik-index-bound - DFPT_Pert::apply_vr public (screened response potential on all bands, FFT-cell triple core shared with real_space_dv) - DFPT_Rho::reset_mixing per displacement; build_occ_kq folds k+q onto the GS k list; solve_displacement full SCF inner loop (v_hartree_q + xc_->apply -> RHS -> Sternheimer -> drho -> mix) - run(): q=0 response + per-irrep displacement loop + assemble / diagonalize / add_loto; null-bases skeleton fallback kept Esolver layer (C7b): - ESolver_DFPT_PW: static config + inp-captured scalars in before_all_runners (rule 1: no global record re-read), run_gs -> init_dfpt wiring after SCF convergence (veff_smooth row, wg, ekb, psi, XC_First_Order_FDM adapter splitting Re/Im through PotXC_FDM) - esolver.cpp factory 'dfpt' branch; read_inp_sys esolver_types + docs/parameters.yaml + input-main.md updated Verified: ctest 12/12 (CELL 4 + DFPT 8); abacus_pw_para links; -h esolver_type shows dfpt; --version v3.11.0-beta8. Governance: 1 allowed exception (determine_type factory PARAM read, existing pattern) + known header/docs WARNINGs.
…ifference, per-displacement reset) Three fixes verified against finite-difference references on the diamond two-atom smoke case (optical 742.367x3 cm^-1 vs FD ~742, acoustic 6.40x3, ASR residual 3.1e-6, off-irrep elements ~1e-11): 1. compute_drho: replace the in-place G-space Hermitian completion (double-processing each +-G pair, breaking Hermiticity and leaking a ~1.25x uniform overshoot) with a real-space 2 Re a(r) presymmetrization before real2recip; one-sided sticks whose -G falls outside the sphere now also complete correctly. 2. XC_First_Order_FDM: the forward difference Vxc[rho+drho]-Vxc[rho] carries a curvature term ~Vxc''*drho^2/2 that leaks a spurious A1 component into v_sc (violating the A1xT2xA1 selection rule by 1.7e-2 Ry/bohr) and destabilizes plain mixing at beta=0.7; use an eta=1e-6 central difference instead (leak ~1e-11, default mixing converges). 3. solve_displacement: zero the stored drho_g when (re)entering a displacement so the previous response (or diverged leftovers) cannot leak into the first screening iteration. Also includes the design-phase debug instrumentation used for the diagnosis (DFPT_DEBUG/PTCHK/DYNCHK/MDBG/dump blocks, DFPT_MIX_BETA env knob) and removes the VQCHK block that read PARAM.globalv.dq/nqx (governance: keep the PR-level global dependency budget non-increasing). Verification: ctest 10/10 (build/, MODULE_DFPT* + little_group + klist); governance --staged clean except advisory warnings; smoke rerun after VQCHK removal reproduces frequencies.
…ness)
The late-iteration divergence diagnosed in the diamond smoke case is a
plain-mixing stability issue, not a physics bug: residual stalls at 5e-5
then grows at exactly 1.2765x/iter while the iterate norm stays constant
(junk direction orthogonal to the physical component). The eigenmode is a
real Hermitian A1 breathing mode on the smallest G shells ({200} 6-vector
equal real amplitudes + {111} 8-vector +-pi/4 phases). A homogeneous
probe (inject the pure A1 trial, drop dV_ext from the rhs, measure the
one-iteration linear map; DFPT_JPROBE / DFPT_JPROBE_NOXC) gives
lambda_A1 = -2.229 (Hartree-only -3.180, XC reduces it to -2.23)
i.e. the Coulomb stiffness 4pi/G^2 at small G. Plain mixing needs
beta < 2/(1+|lambda_min|) ~ 0.62; the physical T2 mode (lambda = -1.42,
less small-G head content) happened to converge at 0.7, which is why the
fixed point was correct while the A1 channel diverged (also explains the
earlier beta=0.3 convergence and the polluted drho manifest).
Default beta is now 0.4 (margin up to |lambda| ~ 5). Verification at
default settings: all six displacements exit via the convergence flag
(~38 iterations average, 228 total), frequencies identical to the
beta=0.7 forced run (optical 742.367 x3, acoustic 6.40 x3; fixed point
independent of beta), ele rows unchanged (e11 0.00286804 vs target
0.0028685, e12 -0.00286494 vs -0.0028701), converged drho manifest now
clean against the finite-difference reference (ratio 0.99994, cos
0.9993, 3.8% pointwise). ctest 10/10 (MODULE_DFPT* + little_group +
klist); governance --staged clean except advisory warnings. Proper fix
is a Kerker-type preconditioned mixer, noted for the B-phase follow-up.
Also adds env-gated design-phase diagnostics used for the diagnosis:
per-iteration residual print, MDBG dumps of drho/v_sc/v_ha/gcar, and the
JPROBE homogeneous-probe path.
…conv_thr/max_iter/mix_beta) - read_inp_dfpt.cpp: 7 new INPUT items with checks (loto requires compute_q0) - esolver_dfpt_pw: drop hardcoded qmesh/conv/max_iter and the dfpt.in stub; wire from inp explicitly (rule 1) - DFPT_PW: set_qfile/set_mix_beta/set_compute_q0/set_loto; q file overrides the MP q mesh in init - QList::read_from_file: fill the fallback A1 placeholder irrep (nirr=1) instead of clearing, so the q-file path keeps the 3N displacement fallback - docs/parameters.yaml + input-main.md regenerated (new category) - README example updated
The four serial suites were last green against pre-calibration binaries; three distinct reference gaps surfaced after the full rebuild: - pert/q0 AnalyticDVloc and FD references: the a004742 phase flip (GS stru_fac convention exp(-i 2pi g.tau), dVloc/dtau = -i (Delta+q)_alpha tpiba Vloc exp(-i 2pi (Delta+q).tau)) was not mirrored in the closed-form references. - rho brute-force G-space and real-space manifests: compute_drho now carries the GS density normalization w/omega (elecstate rhoBandK w1); references divide by omega accordingly. - phon accumulate_electron reference: same phase flip, plus the dynmat mass normalization /sqrt(m_a m_b) (term2) and /m (d2V) that the closed form had silently omitted (fixture mass 12). MODULE_DFPT serial suites 26/26; full regression filter 14/14 (CELL 4 + DFPT 8 + IO 2). Governance: pre-existing exempted include warnings only.
…r cliff
Two independent defects broke DFPT responses whenever the ground-state
k list held more than one inequivalent point (nk > 1):
1. build_occ_kq assumed the k+q and k(q) balls share FFT-cell G labels.
When k+q folds onto a different label of the same physical point
(e.g. lists holding both (1/2,0,0) and (-1/2,0,0)), the projected
states became garbage and the Sternheimer solve diverged. Balls are
now matched through reciprocal-lattice integer triples
f + dn = f', with dn = k(ik)+q-k(ikq); the ikq-side labels are read
through PW_Basis_K::getgcar because collect_local_pw(erf) rebuilds
gcar into a per-k ball layout [ik*npwk_max+igl], destroying the
parent global-ig layout the old code indexed.
2. The absolute wg < 1e-8 occupied-band cliff made the Sternheimer
projector jump between k samplings: a smeared Fermi-tail band with
weight ~1e-6 sits on either side of the threshold depending on the
sampling's Fermi level, opening or closing its empty-state channel
in (H-eps)^-1 and shifting converged force constants by ~10%.
A shared dfpt_band_occupied() now classifies a band as occupied
iff wg(ik,ib) > 0.5*wg(ik,0) (majority occupation), applied
consistently in the projector build, the solve driver, the response
density, the 2n+1 assembly and the q0 valence/conduction split.
Diamond-Si 2-atom validation against finite differences (sym=0):
- single Gamma: D00 0.0208553 vs FD 0.020854 (unchanged)
- single L: D00 0.0129282 vs FD 0.012927 (new FD reference)
- {L,-L}: equals single-L exactly (was divergent), ASR row sums ~1e-6
- {Gamma,L}: D00 0.0166416 vs FD 0.016642 (was 0.0182462, +9.6%)
- {L,X} and weight-skewed {G,L} variants consistent; 14/14
MODULE_DFPT/CELL/IO serial regressions pass.
An unshifted 2x2x2 mesh of diamond Si with the default gauss sigma 0.015 Ry places the smearing Fermi level 1.3 sigma below the Gamma VBM (band occupations 0.92), and finite differences of the same ground state then give force constants ~2.8x softer than DFPT: the E_f response (d mu / d tau channel) is included automatically in any finite-difference ground state but has no counterpart in the Sternheimer flow (DFPT_Metal is a design-phase stub, C4). Without a guard the run converges cleanly and reports silently wrong numbers. DFPT_PW::init now scans the final wg and quits with an explicit message when any band sits measurably between 0 and its full reference (relative weight in (1e-3, 1-1e-3)); negligible gauss tails are tolerated as the insulator limit. Validation matrix for the regime boundary (diamond Si 2x2x2, sym=0): - sigma 0.015: Gamma VBM 92% occupied -> guard fires (was 2.8x off FD) - sigma 0.007: VBM 99.92% occupied -> guard passes, 3.8% off FD (residual dmu channel scales with tail weight) - sigma 0.005: VBM 99.9996% occupied -> 0.05% off FD (insulator limit; D00 0.0127458 vs FD 0.012739), off-diagonals and ASR exact Also validated in this round: single k=0.25,0,0 (D row0 real parts match FD to 6e-7; imaginary antisymmetric parts are the expected one-sided-k Hermitian artifact, the physical force constants are the real parts), and single k=0.5,0,0 with symmetry=0 now reproduces the L-point reference bitwise (symmetry=1 changes the single-k ground state itself and is out of scope for FD comparison). 14/14 MODULE_DFPT/CELL/IO serial regressions pass. MPI>1 smoke (-np 2) aborts with MPI_ERR_TRUNCATE in the DFPT phase: distributed layouts are not yet supported and fail loudly.
…l BZ weight compute_eps/compute_born divided the band sum by nk, but wg(ik,v) already contains the full k weight wk times the spin factor 2, so the stored-k sum is itself the BZ average. The extra 1/nk was a no-op for Gamma-only runs (nk=1) and scaled down multi-k results by 1/nk. Validation (Si diamond, LDA pz): 4x4x4 sym1 (8 IBZ k) eps_inf diagonal mean = 12.6661; sym0 full-BZ 36 k manual sum = 12.6662 (5-digit cross-mesh agreement; LDA reference ~12.7-13.2, experiment 11.7). Retained the env-gated DFPT_Q0DBG p-matrix dump used for the parity-selection-rule audit. Also documents in PLAN: wfc txt writer G-block (igl2isz FFT-stick order) vs coefficient order (psi-ig) mismatch that invalidates file-based element-level cross-checks, and the O_h parity selection-rule evidence that the in-code p matrices are correct.
…ar bug, P0-3 B0 closeout, B2-B4, cleanup, A)
…elta/3
Physics (intake of the uncommitted 5-file fix, part 1 of 2):
- d2vloc_r: both displacement dressings e^{iqR} act on the same atom, so
the cell sum collapses to G = 2q (mod ints); the local second-order
kernel is nonzero only when 2q is reciprocal and then equals the plain
q=0 integer-G kernel. Drop the dead q_cart parameter.
- apply_d2vnl: the second-order nonlocal operator carries wavevector 2q;
build it on the q_eff = fold(2q)-shifted ball and gate the |dbeta><dbeta|
middle projector term behind an explicit include_middle switch.
- accumulate_electron: apply the 2q-reciprocal gate to the whole d2 term
(momentum-forbidden at generic q), pass q_eff/include_middle through,
and fix the ion_ion same-image self term by removing the delta/3
G=0 isotropic piece (validated element-wise against finite differences
of the erfc-split Ewald energy in a q-commensurate supercell).
- ion_ion doc comment updated to the validated closed form.
Tests (dfpt_phon_serial):
- AccumulateElectronAnalyticContraction expectation synced to the
Hermitian 2n+1 accumulation convention (commit dc82fac) and the
gated-off d2 term at generic q; extract SetupBases(k, q) helper so a
test can re-init the fixture at another (k, q).
- New AccumulateElectronD2GateOffGenericQ: row 0 stays pure cross at a
generic q (gate suppresses the forbidden term).
- New AccumulateElectronD2CommensurateQ: k = (-1/2,0,0), q = (1/2,0,0)
so 2q is reciprocal; three-component psi pins the cross term and the
full d2 kernel K_{ab}(G) = -tpiba^2 G_a G_1 Vloc(G^2) e^{-i2pi G.tau}
including the K(G_i - G_j) negative-harmonic convention.
- Zero the Psi buffers after construction (psi::Psi allocates
uninitialized memory); without this the tests read heap garbage and
become order-dependent in the shared-process serial suites.
Part 2 of 2 of the uncommitted-fix intake: env-gated diagnostic probes for the P0-2 Z* investigation and B-phase A/B debugging, all no-ops when their env vars are unset (tracked for cleanup in PLAN_dfpt_implementation.md probe ledger): - DFPT_ZDBG (dfpt_q0 compute_born): per-occ-state decomposition of the Born-charge summand (wg, energy denominator, dV matrix element, position matrix element) to split occ-occ vs valence contributions. - DFPT_BPT (dfpt_pw): perturbation-theory cross-check of the Sternheimer solve, <dpsi|rhs> vs sum_m |<psi_m(k+q)|rhs>|^2/(e_m-e_n) over the empty manifold at k+q (empty_kq_ cache added). - DFPT_NOSC (dfpt_pw): zero the screened potential to isolate the bare Sternheimer response. - DFPT_D2MID / DYNCHK d2gate (dfpt_phon): disable the |dbeta><dbeta| middle projector term; print the 2q-reciprocal gate decision. - DFPT_XB: extend the row/column selection to the 2-atom rows 6. Verified: MODULE_DFPT phon 9/9, q0 5/5 serial suites with probes inert.
Symmetry-reduced k sums of the q=0 susceptibility tensors must be star-averaged: the partial at a rotated star member Rk is R chi(k) R^T (cartesian column form), with atom-resolved Born partials credited to the image atom under the paired direct-space operation. The row-form operator G^-1*kgmatrix*G from the kvec_d row convention had been fed to rotate_tensor untransposed, which breaks the star sum (right- vs left-coset representatives), so store the transpose. Diamond Si 4x4x4 verification (sym=1): eps_inf = 12.6661*I and Z* = 15.5799*I per atom, both bit-consistent with the symmetry-off full-mesh reference (off-diagonals ~1e-14; previously 13.78/15.34/8.88 anisotropic). The remaining Z* offset vs the diamond target 0 is the known missing-screening formula defect, tracked as the next P0-2 item. Add StarRotationCyclicGroup to dfpt_q0_serial (C3 orbit cell: star size, anisotropic trace-6 tensor averaging to 2*I, cyclic atom maps, identity fallback) and the DFPT_STARDBG probe; build_stars/rotate_tensor/stars_ move to public for the test.
…r_eu cross-check probe - solve_pos_resp + compute_born v4: Y^a = (H-eps_v)^-1 P_c [H,x_a]|psi> (velocity rhs, build_vkb_dk nonlocal part), Z* = zion delta - 2 sum wg Re <dpsi^kappa,scf|Y^a> (QE add_zstar_ue form); pos_resp/ dpsi_efield stashes in DFPT_PW_Data - eps factor 2 fix: 16 pi / Omega per QE dielec.f90 (8 pi was half); ComputeEpsTwoLevelAnalytic expectation synced, serial 6/6 - DFPT_ALEG probe: E-field SCF fixed point (solve_e form) + zstar_eu A-leg vs zstar_ue B-leg cross-check + SCF eps + DFPT_PTCROSS bare cross spectral diagnostic - validated vs locally built QE 7.2 (same UPF/cell/ecut/mesh): GS energy identical, Gamma-TO 517.5/517.6 vs 517.63 (0.03%), Z* -1.19928 vs -1.19765 (0.14%), eps_scf 23.6825 vs 23.6685 (0.06%); 4x4x4 anomaly (Z*=-1.2, eps~23.7 vs lit 13) shown to be shared k-mesh convergence by QE discriminators (ONCV@4x4x4 same, pz-vbc@8x8x8 -> 14.04/-0.09) - PLAN P0-2 closed with validation matrix and re-scoped acceptance
… guard) The develop-side DFT+U refactor (deepmodeling#7852-deepmodeling#7867) removed source_lcao/module_dftu/dftu.h and the is_locale_initialized() member, which broke every CMake build configuration of this branch at dfpt_pw_data.cpp (all 9 CI build variants plus Test/CUDA/abacuslite failed at the compile step; only the Makefile job passed because the Makefile.Objects DFPT entries were absent at that merge point). Changes: - DFPT now consumes the PW-side Plus_U_Base (source_pw/module_pwdft/ dftu_base.h) instead of the LCAO-side Plus_U header: dftu_ member, DFPT_PW_Data::init / DFPT_PW::init signatures and get_dftu() all use const Plus_U_Base* (the esolver call site passes &this->dftu with an implicit upcast). This removes the PW -> LCAO cross-layer include. - u_active() = with_u() && is_occ_mat_initialized(): the reservation usability now follows the occupation-matrix state of the provider. - DFPT_PW::init rejects a wired provider explicitly (WARNING_QUIT): the ground state supports PW-basis DFT+U now, but every DFPT U hook (cal_docc, build_dv_u, dftu_onsite, born/docc contractions) is a no-op U0 reservation, so running anyway would silently drop the whole first-order U response (fail-loud, same pattern as the metallic-sampling guard). - test/dftu_test_support.cpp rewritten: the old static-member replicas no longer exist; the shim now provides only the Plus_U_Base ctor/dtor (also linked into MODULE_DFPT_pw_data_test, which constructs the provider directly). dfpt_pw_run_test's locale test becomes a death test pinning the WARNING_QUIT guard; the with_u/u_active contract moved to DFPT_PW_DataTest.DftuReservationProviderUsability; the unused dftu.h includes dropped from the phon/q0 serial tests. Verification (GNU 8.3.1 + OpenMPI 5.0.3, GCC13 no-MPI cross-check): - cmake --build build --target abacus_pw_para: builds/links - cmake --build build-nompi (-DENABLE_MPI=OFF -DENABLE_LCAO=OFF) --target abacus_pw_omp: builds/links - ctest -R 'MODULE_DFPT|MODULE_CELL': 50/50 pass (incl. the new death test and provider-usability case); ./build/abacus_pw_para --version prints v3.11.0-beta8 - agent_governance_check --staged: no findings
The governance checker blocks the PR while the diff's added lines carry more GlobalV/GlobalC/PARAM references than the removed lines (added=51, removed=22, net_delta=+29 -> CI 'Governance checks' exit 1). 30 of the added references were test-side GlobalV::ofs_running streams; they now use local std::ofstream objects (the fixture members that already existed), and ReciprocalGrid::print_klists prints through its own ofs parameter instead of the global stream (its single caller passes the same running log). The stale 'Originally GlobalV::FINAL_SCF' comment wording is dropped. The remaining production-side references (reciprocal_grid.cpp k-point-file echo, klist.cpp MY_RANK guards) are line-for-line moves of the previous klist.cpp code, so the budget is now non-increasing (net_delta = -4). Verification: ctest -R 'MODULE_DFPT|MODULE_CELL' 50/50 pass; abacus_pw_para relinks; agent_governance_check --base origin/develop --head HEAD exits 0 (no BLOCK findings).
…eciprocalGrid refactor
The ReciprocalGrid refactor (Phase 1-3 of this PR) made K_Vectors
polymorphic: its vtable is now keyed on K_Vectors::renew and emitted in
klist.cpp, and the base vtable lives in reciprocal_grid.cpp. Twelve test
targets across estate/hsolver/stodft/io instantiate K_Vectors but never
compiled those translation units, so they fail to link after the merge
(masked until now by the earlier dftu.h compile break):
- MODULE_ESTATE_elecstate_{print,base,pw,energy}
- MODULE_PW_Sto_Hamilt_UTs
- MODULE_HSOLVER_pw
- MODULE_IO_write_bands (test_serial)
- MODULE_IO_write_eig_occ_test / write_dos_pw / print_info /
read_wf2rho_pw_test (already had klist.cpp, lacked reciprocal_grid.cpp)
- MODULE_IO_write_dmk
Mirrors the pattern already used by this PR's own klist/qlist tests:
add klist.cpp + parallel_kpoints.cpp + k_vector_utils.cpp +
reciprocal_grid.cpp to SOURCES and the symmetry lib to LIBS.
Verified: full build green except MODULE_IO_numerical_basis_test (needs
ENABLE_LCAO, unguarded on develop as well); the fixed tests pass under
ctest; remaining local failures are environment artifacts (ScaLAPACK
abort-stub, ELPA off).
…dd new DFPT objects to Makefile.Objects The ReciprocalGrid refactor made K_Vectors polymorphic (its key function and the base vtable now live in klist.cpp / reciprocal_grid.cpp), so any test that instantiates K_Vectors (module_dm tests, deltaspin spin_constrain/template_helpers via spin_constrain.cpp, and init_dm_from_file via density_matrix_io.cpp) fails to link. Also register the five PR-added translation units (reciprocal_grid.cpp, little_group.cpp, read_inp_dfpt.cpp, dfpt_hamilt_shift.cpp, dfpt_kq_basis.cpp) in source/Makefile.Objects so the Intel Makefile build does not fail with undefined references.
…m generator The DFPT parameter block was hand-placed at a position that differs from the item_dfpt() registration order, so the byte-exact consistency checks in test.yml (--generate-parameters-yaml / generate_input_main.py) fail. Regenerate both files with the documented commands to restore sync; the only change is the position of the DFPT category block.
The ReciprocalGrid refactor made K_Vectors derive from ModuleCell::ReciprocalGrid, so klist.cpp.o and k_vector_utils.cpp.o now reference ReciprocalGrid member functions and its vtable. The deepks_unit_support object library (DEEPKS_UNIT_COMMON_SOURCES, gated behind ENABLE_MLALGO and thus only compiled in the gnu Test CI job) compiles klist.cpp without reciprocal_grid.cpp, failing to link all 30 MODULE_LCAO_DEEPKS_* test executables with undefined references to ModuleCell::ReciprocalGrid::renew/Monkhorst_Pack/build_star_ops/... and its vtable/typeinfo. Add the missing translation unit to the common source set; the symmetry library (incl. little_group.cpp) is already on the link line.
…ded-process deadlock MODULE_DFPT_pw_run_test timed out (1700 s) in the gnu Test CI job: the two irrep-loop tests run first execute OpenMP regions, so with the job's OMP_NUM_THREADS=2 the process is multithreaded when the third test (dftu-reservation EXPECT_EXIT) forks. The default fast-style child then deadlocks on exit and the parent waits forever (reproduced locally under OMP_NUM_THREADS=2: gtest warns 'detected 2 threads' and hangs). Switch all three DFPT death tests to the fork+exec threadsafe style (same pattern as module_container tensor_test). For the pw_run test also bridge std::cout to std::cerr inside the death statement: WARNING_QUIT prints the NOTICE block to stdout, while death tests match the child's stderr; the old CaptureStdout+HasSubstr assertion cannot see the re-exec child's output. Verified under OMP_NUM_THREADS=2: pw_run 3/3 in 0.3 s (previously indefinite hang), kq_basis 5/5, pert_serial 8/8, and the full MODULE_DFPT ctest batch 8/8.
…t_test_mocks.cpp (mirror tmp_mocks.cpp convention); absorb dftu_test_support.cpp
…en pw_data and pw_run tests (dfpt_stru_fixture)
…hared DFPTSerialBase (cell/basis/data setup, Coulomb/NC atom builders, analytic dVloc reference)
Share the occupied-weights table, the single-plane-wave psi builder, the analytic accumulate_electron cross term (now on top of AnalyticDVloc), and the isotropic loto data setup (eps/Born charges + two-atom mass table via MakeTwoAtomCell) through phon fixture helpers; the three AccumulateElectron tests and the two loto closed-form tests keep their reference formulas but drop the duplicated inline copies.
…#7894 更新为准） # Conflicts: # source/source_cell/k_vector_utils.cpp
The deepmodeling#7894 refactor moved the k+q congruence matcher into match_commensurate_kq but inverted its dn convention (dn = k_d(ikq) - k_d(ik) - q) while the folding in copy_occ_state_ball still used key = G + dn from the pre-refactor convention (dn = k_d(ik) + q - k_d(ikq)). Plane-wave identity requires G' + k_d(ikq) == G + k_d(ik) + q, i.e. G' = G - dn, so for every k whose k+q folds across a BZ boundary (dn != 0) the occupied-state coefficients were attached to plane waves shifted by 2*dn. The Sternheimer operator then developed ~+tpiba^2*|2dn|^2 errors in <psi|T+Vnl|psi> (~4 Ry for Si@X, 4x4x4 mesh) and the CG solve diverged, yielding NaN phonon frequencies for every non-Gamma q (48 of 64 k-points affected at Si@X). The Gamma path (dn == 0 always) was unaffected. Restores the pre-refactor mapping verified by the Si L-point reference case.
…auses for GCC 8 GCC 8 rejects listing a const-qualified pointer parameter (predetermined shared) in default(none) shared(...) clauses; GCC >=9 accepts both forms. Removing it from the list is valid on both.
…ier term Replace the amplitude-grid finite differences with an analytic kernel built from libxc second derivatives whenever ABACUS is compiled with libxc (ENABLE_LIBXC). The response follows the QE dgradcorr algebra with the total-field gradient i(G+q), so the GGA carrier contribution i(q.f_grad)A -- exact at any q, structurally absent from grid-only finite differences of the carrier-free amplitude, and the root cause of the systematic q != 0 GGA phonon errors (mgo X TA error -84%) -- is included. LDA keeps the plain f_rho A path. Without libxc the FDM adapter stays the default and now prints an explicit warning that the GGA response at q != 0 lacks the carrier term. - XC_First_Order::set_context(q_cart): per-q context hook; the solver injects it in assemble_v_sc before apply - XC_First_Order_Analytic: kernel fields evaluated once on the ground-state density (libxc vrho is the potential form; GGA chain factors 2/4 for sigma, no rho renormalization), Ry units - unit test: LDA kernel vs FD of v_xc, GGA q=0 full chain vs field perturbation FD, carrier presence at finite q (GGA) and exact q-independence (LDA); existing DFPT serial tests unchanged/pass
PW_Basis gcar / ucell G are dimensionless (latvec inverse, 2 pi not included); the physical wavevector is tpiba * gcar (same convention as DFPT_Rho::v_hartree_q, whose Coulomb factor divides by tpiba2). The analytic kernel gradient machinery (ground-state grad rho0, the i(G+q) response gradient and the divergence) must use physical units: without the factor the sigma fed to libxc is tpiba^2 too small and the gradient terms are scaled inconsistently -- LDA (no gradient path) matched its anchor to 0.003 cm^-1 while GGA q != 0 frequencies came out 2-3x too stiff (mgo X PBE 821..969 vs QE 282..603). The unit-test gradients are fixed the same way so the reference feeds libxc physical sigma (the previous self-consistent dimensionless loop hid the bug).
Dropping rho_basis from the default(none) shared() lists accommodated the local GCC 8.3.1 toolchain (which predetermines const-qualified variables as shared and rejects listing them), but every CI compiler (GCC >= 9, Intel oneAPI) treats the pointer parameter as ordinary and requires it listed -- the removal broke all build/test jobs. Restore the file byte-identical to upstream develop; GCC 8 is not a supported baseline and stays a local-only concern.
PotXC_FDM linearizes the internal functional around rho0 + rho_core (pot_xc_fdm.cpp adds the ground-state core charge of nlcc pseudopotentials); the analytic kernel was fed the valence rho only, which moves the kernel fields of nlcc systems off their ground-state point -- bn_cubic Gamma optical modes dropped ~15 cm^-1 versus both the FDM path and QE. Evaluate the libxc kernels at the same valence-plus-core density; the first-order density stays valence-only (frozen core).
The first-order XC kernel has no tau channel (meta-GGA) and no exact-exchange response (hybrids); running DFPT with such functionals was undefined behavior on the finite-difference path (the analytic kernel already rejected them in its constructor). Reject func_type other than LDA/GGA at init with a message naming the functional.
mohanchen
approved these changes
Oct 17, 2026
mohanchen
added
DFPT
Collaborator
|
@ABACUS-PR-BOT could you review the PR? |
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
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) |
Collaborator
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) |
Collaborator
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; |
Collaborator
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); |
Collaborator
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 |
Collaborator
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" |
Collaborator
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) |
Collaborator
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) |
Collaborator
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.
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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.