feat(dfpt): complete the nlcc core-response terms (QE dynmat_nlcc + dynmatcc) and fix the mirrored structure-factor phase by Zanthoxylum · Pull Request #8237 · deepmodeling/abacus-develop

Zanthoxylum · 2026-10-29T13:24:44Z

Summary

The DFPT phonon implementation was missing the nonlinear core correction (nlcc) contribution to the force constants: the first-order core-density response drhoc entered the SCF screening channel (merged with #8182) but the two direct 2n+1 terms were absent, and the core-response structure factor carried a mirrored phase. This PR completes the nlcc implementation, matching QE 6.8 phonon frequencies to ~1 cm⁻¹ on both validation systems (previously off by 30-300 cm⁻¹).

Three changes, all pinned by element-level dynamical-matrix comparisons against QE dyn.G files:

The two missing direct terms (QE dynmat_nlcc + dynmatcc equivalents):
- Kernel term: D_ab += (Ω/nrxx) Σ_r dvxc^a(r) conj(drhoc^b(r)) · (1/2)(X_ab + conj(X_ba)), with dvxc^a = f_xc(drho_val^a + drhoc^a) built from a per-displacement table covering ALL atoms (coreless rows need their valence-only dvaux; leaving them empty silently dropped the entire one-sided cross block)
- Diagonal term: D_aa[i][j] += ∫ dr w_ij(r) v_xc⁰(r) with w_ij = -tpiba² G_i G_j ρ̃_c(|G|) e^{-iGτ}, evaluated through recip2real as a real field; v_xc⁰ computed in the esolver wiring and delivered via DFPT_PW_InitContext::vxc_r
Structure-factor phase fix (compute_drhoc_g): the phase was exp(+i(G+q)τ), placing the core response at the mirror site -τ instead of the atom position. For wurtzite B at (1/3,2/3,1/4), the mirror lands exactly on an N site of the other layer, producing spurious same-layer overlaps that collapsed the interlayer E2low shear to -399 cm⁻¹ (QE: -29.8). Fixed to exp(-i(G+q)τ) (the QE addcore convention). Every historical validation was blind to the sign: bn_cubic's B sits at the origin where both conventions coincide.
Regression tests replacing the temporary matrix-dump diagnostic:
- DrhocLocalizesAtTheAtomNotItsMirror: nlcc atom at non-origin τ = (0.25, 0.31, 0.17); |drhoc|² centroid must sit at the atom, not the mirror (catches the phase-sign bug directly; includes a generic non-zero q variant)
- NlccDirectTermsFillSemantics: synthetic kernel table (no XC kernel needed); verifies the ½-symmetrized fill, the non-nlcc row contribution, and Hermiticity

Also: DFPT_Rho::compute_drhoc_g extracted as the pure core-response builder; assemble_v_sc takes the q index explicitly (the drhoc slot is per q; last_q_ can be stale on the first assembly of a new q); per-displacement density stash (set/get_drho_disp) mirrors the existing v_sc stash.

Validation (identical settings as QE 6.8: ONCV/PBE, ecut 20/80 Ry, 4×4×4 k-mesh, Γ-centered):

bn_cubic Γ TO frequencies:

	QE	ours	Δ
TO	1328.040	1328.033	0.007
(B,B)xx	0.814388	0.814599	2.1e-4
(B,N)xx	-0.814659	-0.814584	7.4e-5

bn_hex Γ optical modes (the C5 outlier, rms was 28.65 → now 1.20 cm⁻¹):

mode	QE	ours	Δ
E2low (×2)	-29.79	-32.18	-2.39
A1 (TO)	236.96	237.71	+0.76
E2 (TO)	531.07	530.33	-0.74
E1 (×2)	1666.44	1666.33	-0.11
LO	1783.83	1783.43	-0.40

Multi-q validation: bn_cubic X rms 0.04 cm⁻¹, bn_hex A rms 0.82, bn_hex M matrix blocks all ≤6.4e-4. The E2low residual (-2.39) traces to the non-nlcc baseline ((N,N) block -8.4e-4 = 0.065% relative) amplified ~50× by the soft-mode sensitivity — the same -2.31 appears in the no-nlcc pseudo comparison.

Non-nlcc invariance: Si Γ frequencies bit-identical before/after this PR's changes.

dfpt_mix_type (plain|kerker) and dfpt_kerker_a2 replace the design-phase env knobs DFPT_MIX_TYPE / DFPT_KERKER_A2 / DFPT_MIX_BETA; the plain 0.4 default behavior is unchanged. Also emit an explicit ModuleBase::WARNING when a displacement response exhausts dfpt_max_iter (warning.log channel, next to the stdout summary from the R4a fix in deepmodeling#8163). Rebased from the dfpt-mix branch onto current develop (resolves the logging-format conflict by combining both versions).

…den) New INPUT parameter dfpt_broyden_ndim (default 8). When dfpt_mix_type = broyden, the first-order density mixing uses the Modified Broyden method (Johnson 1988) from the existing Base_Mixing::Broyden_Mixing -- the same algorithm QE ph.x uses (nmix_ph = 4 default there). The Broyden history is per (q, displacement-cycle) and is dropped by reset_mixing. The three-step Broyden interface (push_data + cal_coef + mix_data) is wired inside mix_drho; the plain and kerker paths are unchanged.

…_sc) When dfpt_mix_type = v_sc, the SCF fixed-point variable changes from the first-order density drho to the screening potential v_sc: each iteration solves the Sternheimer equation with the mixed potential, obtains the output density, rebuilds the output potential from it, and mixes the two potentials. The variable change provides implicit K^-1 preconditioning (the Hartree+XC kernel's inverse), flattening the Coulomb-stiff smallest-G shells that destabilize plain density mixing (tio2 X response Jacobian eigenvalue ~ -2.2 on the {111}/{200} shells). The residual in v_sc mode measures the potential self-consistency ||v_out - v_in||/||v_out||; the 2n+1 accumulation path uses the converged output density from the converged potential. plain, kerker, and broyden paths are unchanged.

The synthetic Broyden convergence test crashes inside Broyden_Mixing::push_data with the mock SetupBases fixture (likely an internal interface mismatch with the null screen/combine callbacks); the real-system G5 validation (si/tio2 X with broyden) provides end-to-end coverage of the mixing path. Re-enable after the mock interface is understood.

…d DFPT_PW setter Two validation points silently rejected the new mixing types: DFPT_PW::set_mix_type only accepted plain/kerker (broyden/v_sc were dropped, leaving the default plain), and DFPT_Rho::init WARNING_QUIT on any mix_type outside plain/kerker. Both now accept all four.

Broyden_Mixing::tem_push_data calls mix(data, data_in, F) unconditionally (broyden_mixing.cpp:47 has a nullptr guard for screen but NOT for mix); passing an empty std::function threw bad_function_call -> SIGABRT. Provide the actual preliminary mix (data = data_in + beta * F) that the Broyden algorithm needs. tio2 X now converges at iteration ~18 (first time this system converges in DFPT response SCF; plain and v_sc both oscillate at residual ~0.8-1.7).

…ments reset_mixing destroys the per-q Mixing_Data object but the Broyden mixer's internal 'address' pointer still referenced the freed object; the next displacement's push_data then hit the 'One Broyden_Mixing object can only bind one Mixing_Data' guard and WARNING_QUIT. Call broyden_mixer_->reset() alongside the mdata reset to clear the dangling binding.

…t; add v_sc Broyden dfpt_mix_type values renamed to a consistent naming scheme: plain -> rho_plain (default) kerker -> rho_kerker broyden -> rho_broyden v_sc -> v_sc_plain (new) -> v_sc_broyden v_sc_broyden uses a separate Broyden_Mixing instance on the nrxx-length real-space screening potential (the density Broyden operates on npw-length G-space arrays). The v_sc Mixing_Data is initialized and reset per displacement cycle in solve_displacement. The v_sc_* modes bypass DFPT_Rho's density mixing entirely (a rho_plain config is passed to DFPT_Rho as a no-op) and do their own potential mixing in solve_displacement's v_sc branch.

G5 validation across si X and tio2 X: both rho_broyden and v_sc_broyden converge identically (same fixed point, ~13-17 iterations per displacement; tio2 X converges for the first time ever in DFPT while rho_plain diverges at 12+ unconverged displacements). v_sc_broyden is chosen as the default because it exactly matches QE's approach (Modified Broyden on the screening potential), provides better-conditioned Broyden inner products through the implicit K^-1 weighting, and has theoretical headroom for nonlinear extensions. All five mixing methods remain available via dfpt_mix_type: rho_plain | rho_kerker | rho_broyden | v_sc_plain | v_sc_broyden

- migrate the remaining old mix_type tokens: the rho_kerker branch in mix_drho, the serial-test fixtures, the Impl default (now the valid v_sc_broyden), and the read_inp_dfpt empty-input fallback (was the rejected "plain"); CI MODULE_DFPT_rho_serial died on the fixture token - reset the shared Broyden mixer on every reset_mixing call, not only when the q already has a history entry: entering a fresh q left the mixer address-bound to the previous q's Mixing_Data and the first push died in the one-object-one-data guard - reduce the Broyden cal_coef inner products and the residual norms over the pool (Parallel_Reduce::reduce_pool behind __MPI guards, xc_pot.cpp convention): each rank only holds its FFT slab / G-shell share, so unreduced Gram systems and stopping norms were partition dependent - guard the relative-defect convergence metric: a zero output norm now only reports convergence when the defect itself is zero, and non-finite norms never converge (v_out = 1 - c v_in counterexample) - enable the Broyden regression with a sourced model fixed point (the sourceless variant cannot converge under the relative residual) and add a two-q reset_mixing regression for the address-guard fix - INPUT metadata: use the accepted rho_*/v_sc_* names in the descriptions, document which residual dfpt_conv_thr measures per scheme, and reword the v_sc comment (plain-mixing spectrum is invariant under the variable change; Broyden benefits from the K-weighted metric) - forward-declare the mixing types in dfpt_rho.h and drop the concrete includes (dfpt_pw_impl.h already includes them directly)

…regression Review follow-up: the binding populated shell is d = 0.5 with factor |1 - beta (1 - D)| = 0.8 (>= 103 plain steps to 1e-10), not the unpopulated-shell bound 0.6. Assertions are unchanged.

- E-field solve (dfpt_pw_q0.cpp): report non-convergence explicitly with the same ModuleBase::WARNING as solve_displacement (a silently truncated E-field solve corrupts Z* / LO-TO); the v_sc_* -> rho_plain mapping note is added at the init site - DFPT_Rho::init: reject broyden_ndim < 1 for rho_broyden (Config is a C++11 aggregate, so brace-initialized callers get 0 and Broyden_Mixing would divide by zero in the % mixing_ndim wrap); fix the 'default 8' comment; accept only the rho_* tokens (the solver-side init already maps v_sc_*, so the nested branch was unreachable) - drop the no-op set_drho_g write-back in the v_sc branch (data_ still holds the last compute_drho output) and reword the comment - drop the out-of-repository MIXING_20261020.md reference - rename safe_relative_defect_ -> safe_relative_defect (free function in an anonymous namespace; the trailing underscore is the member convention); sync the plain-mixing contraction notes in the test (roughly 100 iterations at 0.8/step on the binding shell)

…der v_sc_broyden Review follow-up (Critsium-xy): with the v_sc_broyden default, the E-field solve (Z* / epsilon_inf / LO-TO) silently fell back to plain density mixing at beta 0.4 -- the weakest member of the family exactly on the path that matters most for polar systems. Map v_sc_broyden to rho_broyden for DFPT_Rho (v_sc_plain keeps rho_plain), so the E-field loop runs Modified Broyden on the density while the displacement loop keeps its potential-space Broyden. Verification (si, Gamma + dfpt_compute_q0/dfpt_loto, default INPUT so the path is exercised through the mapping): E-field 3/3 directions converge in 12 iterations (bit-identical residuals to an explicit rho_broyden run) versus 38 with the previous plain fallback -- displacements unchanged at 13 iterations (v_sc Broyden); Gamma frequencies and the LO-TO block reproduce the pre-change values exactly. INPUT metadata documents that dfpt_broyden_ndim now sizes both Broyden instances and that the E-field residual is always the density residual.

…dynmatcc) The nlcc response was previously wired only through the SCF screening channel (f_xc.(drho_val + drhoc) inside assemble_v_sc, QE dv_of_drho semantics, verified line by line). QE ph.x adds two direct terms to the force constants that have no wavefunction representation, and the bn_hex isolation experiment (identical pseudos with the PP_NLCC block zeroed: ABACUS vs QE optical rms 0.9 cm^-1 without nlcc, while the SCF-only nlcc overshoots QE's nlcc effect by 4-12x with wrong signs on the A1/E2 soft modes) shows they are the dominant cancellation partner of the screening contribution: - kernel term (QE dynmat_nlcc): D_ab += Omega/nrxx sum_r dvxc^a(r) conj(drhoc^b(r)) with dvxc^a = f_xc.(drho_val^a + drhoc^a) of the converged displacement a and the core response of every column b; dvxc^a is rebuilt in the two-pass accumulation loop from the newly stashed per-displacement density response plus a recomputed core response (the working slots hold the last solved displacement). - diagonal term (QE dynmatcc): D_aa[i][j] += Int dr w_ij(r) v_xc^0(r) with w_ij the second-order shifted-core amplitude -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau), evaluated through recip2real as a real field so the half-space G-sum convention never enters directly; v_xc^0 is the ground-state XC potential on the same valence+core density the kernel linearizes around, computed once in the esolver wiring (XC_Functional::v_xc) and delivered through the new DFPT_PW_InitContext::vxc_r provider (null disables the term). Also: DFPT_Rho::compute_drhoc_g extracted as the pure core-response builder (compute_drhoc keeps its data-slot semantics on top of it); assemble_v_sc now takes the q index explicitly instead of reading the drhoc slot through last_q_ (stale on the first assembly of a new q); per-displacement density stash (set/get_drho_disp) mirrors the existing v_sc stash. Verification: MODULE_DFPT serial suite 9/9 (the new terms are inert without nlcc pseudos); bn_cubic / bn_hex Gamma-vs-QE validation runs attached in the PR.

…ubic) Three fixes over the first cut, pinned by the element-level matrix comparison against QE dyn.G (bn_cubic: all blocks now agree to 2e-4, TO frequencies to 0.007 cm-1): - the per-q kernel table must fill EVERY displacement slot: rows of non-nlcc atoms need their dvaux too (f_xc applied to the valence response alone) -- leaving their slots empty silently dropped the entire one-sided cross-block X(N,B), the dominant mixed piece - the effective Hermitian element is the symmetrized (1/2)(X_ab + conj(X_ba)) (QE dynmat_nlcc single ordering + final Hermitization); the earlier both-order full fill doubled the same-atom blocks, and the lower-triangle-only variant overshot the one-sided cross blocks - the vxc provider must live in Impl-owned storage (the init context points at esolver stack state); the diagonal block is added once per atom (dir == 0)

Brings the response-mixing machinery (v_sc_broyden default, INPUT parameterization, E-field density-Broyden mapping) under the nlcc direct terms: ~3x faster displacement solves for the nlcc validation matrix (13 vs 38 iterations per displacement). Conflict resolutions in dfpt_pw_solve.cpp: keep the mixing vsc_mode iteration structure and the safe_relative_defect helper, thread the explicit q index through every assemble_v_sc call (the drhoc slot is per q; last_q_ can be stale on the first assembly of a new q).

compute_drhoc_g built the phase exp(+i(G+q).tau), which places the first-order core density at the MIRROR site -tau instead of the atom position tau (QE addcore: exp(-iG.tau) * exp(-i q.tau)). Every validation so far was structurally blind to the sign: bn_cubic's B sits at the origin (both conventions identical) and N carries no nlcc, and the F2 numerical cross-check ran on the same blind system. The bn_hex 16-block matrix comparison against QE dyn.G exposed it element by element: with wurtzite B1=(1/3,2/3,1/4), the mirror of B1 lands exactly on N2=(2/3,1/3,3/4) and the mirror of B2 on N1, so the mirrored core response produced spurious same-layer overlaps ((B1,B2)xx +0.053 vs the true 0.0013 interlayer scale) and lost the real diagonal ones ((B1,B1)xx -0.105), collapsing the E2low interlayer shear to -399 cm^-1 against QE's -29.8. After the one-line conjugation, bn_hex stays bit-identical (B at the origin) and bn_hex matches QE on every block to <= 2.5e-4 (B-involving) with the optical rms dropping from 28.7 (C5, no core response) and ~60 (previous iterations) to 1.2 cm^-1; E2low recovers to -32.2 vs -29.8 and the stiff modes to within 0.8 cm^-1.

Two-level regression suite replacing the temporary matrix dump (whose diagnostic role -- the element-level QE comparison that pinned every formula bug -- now lives as executable invariants): Level 1 (dfpt_rho_serial_test): DrhocLocalizesAtTheAtomNotItsMirror. An nlcc atom at the non-origin tau = (0.25, 0.31, 0.17) (the anti-blind-spot position: every structure-factor phase is active) must have its first-order core response |drhoc|^2 centroid AT the atom. The mirrored phase exp(+i(G+q)tau) -- the bug every historical validation missed because bn_cubic's B sits at the origin -- places it at the mirror site instead. Includes a generic non-zero q variant. Level 2 (dfpt_phon_serial_test): NlccDirectTermsFillSemantics. With a synthetic all-slots kernel table and no dpsi (isolating the kernel contribution from the bare cross term), the one-sided cross block must carry exactly (1/2) conj(X_ba) (the symmetrized fill pinned by the element-level QE comparison), and the accumulated matrix must stay Hermitian. Together these guard against the three fill bugs fixed in the parent commits (dead coreless rows, doubled same-atom blocks, and the halved one-sided cross blocks). The temporary DFPT_DUMP_DYN production block is removed (its commit 045f050 in the branch history retains the identical code for local re-use during Project E matrix comparisons).

…-nlcc-direct Conflict resolutions: keep the nlcc-side q_idx parameter in every assemble_v_sc call and signature (the drhoc slot is per q; the merged develop code from deepmodeling#8214 lacks the explicit index), keep the nlcc phase-localization and fill-semantics regressions in the serial tests.

The vxc_r field carried a '= nullptr' default, which under C++11 makes DFPT_PW_InitContext a non-aggregate; the 12-argument thin wrapper's brace-init then fails to compile in the non-MPI / non-LCAO CI configurations. Remove the default and pass nullptr explicitly as the 13th initializer in the wrapper.

ABACUS-PR-BOT

Two NLCC correctness issues remain: finite-q matrix orientation and MPI cell-integral normalization. Please address the inline findings; three additional notes are nonblocking.

This review covers 2580552557cbaf3408b039e8bd02b664e9c10415 through static source inspection and analytic derivation. git diff --check and source-integrity/LF checks passed; no builds, tests or numerical/MPI runs were performed. The reported QE 6.8 comparisons were not independently verified.

ABACUS-PR-BOT · 2026-10-29T16:29:56Z

+            {
+                continue;
+            }
+            const std::complex<double> s_ab = 0.5 * (x_ab + std::conj(x_ba));


[P1] Match the NLCC block to the existing complex matrix convention

At nonzero q, these contractions leave the row response unconjugated, whereas the wavefunction and Ewald terms use the conjugated row. Even a constant-real-kernel core-only Hessian requires sum(conj(c_row) * f * c_col); the current expression inserts its conjugate and reverses imaginary cross terms, which Hermitian averaging preserves. Please correct the conjugation while retaining the 1/2 average, and add an independent finite-q complex cross-element regression.

Fixed in a2889fe: both kernel contractions now conjugate the ROW field (sum(conj(row)*col)), matching the wavefunction and Ewald terms. The fill-semantics regression test's independent recomputation is updated to the same convention (the imaginary cross terms flip sign, which the old test caught at the element level). At Gamma with real amplitudes this is a no-op (all bn_cubic/bn_hex Gamma validations unchanged).

ABACUS-PR-BOT · 2026-10-29T16:29:56Z

+    }
+    const int nat = ucell_->nat;
+    const ModuleBase::Vector3<double> q_cart = data.get_qvec(q_idx) * ucell_->G;
+    const double fac = ucell_->omega / static_cast<double>(nrxx);


[P1] Integrate the NLCC terms over the whole PW pool

nrxx is the current rank's FFT slab size, so this forms a slab average when the pool has multiple ranks. accum_nlcc_diag has the same issue, and assemble does not reduce either contribution. Please use Omega/nxyz and guarded pool reductions of the local products, with the same collective participation on every rank, then compare raw matrices across rank counts.

Fixed in a2889fe: both the kernel and diagonal terms now normalize by Omega/nxyz (the full grid) and reduce their products over the pool with Parallel_Reduce::reduce_pool behind #ifdef __MPI guards — the same convention as the mixing PR (#8214). The serial test suite (9/9) is unaffected (no __MPI).

The Omega/nxyz scale and pool sums are now in place. A remaining participation case is a decomposition reaching accumulation with nrxx == 0: esolver_dfpt_pw.cpp:480 and dfpt_pw_init.cpp:121 turn that valid empty local XC field into nullptr, so the rank skips accum_nlcc_diag while peers enter its FFT and reduction. Please preserve provider availability independently of local size in both places and cover an empty-slab decomposition. This is a conditional source-level mismatch, not a reproduced hang.

Fixed in 824be4f: the vxc provider now always passes a valid pointer (&vxc_r_, never nullptr on empty local slabs) from both the esolver InitContext and the Impl init_submodules, so every rank enters accum_nlcc_diag's collective path. The size checks inside gate only the local work, not the participation. No empty-slab reproduction available (single-rank DFPT only), but the source-level mismatch is resolved.

ABACUS-PR-BOT · 2026-10-29T16:29:56Z

+        EXPECT_LT(std::sqrt(dist_atom), 0.1) << "dir=" << dir;
+        EXPECT_GT(dist_mirror, 100.0 * dist_atom) << "dir=" << dir;


Non-blocking: these assertions admit an ideal reflected profile. The fixture's fractional center is (0.025,0.031,0.017); its mirror unwraps to -tau, giving a distance of about 0.0866 < 0.1, while comparison to the still-wrapped mirror also passes. Please compare both centers in the same periodic image and choose a position/tolerance that rejects the old sign; the actual mutated test has not been run.

Fixed in a2889fe: the distance computation now uses the minimum-image convention (each component wrapped to [-0.5, 0.5) via round), which correctly distinguishes the atom from the mirror in every periodic image. The thresholds are tightened to dist_atom < 0.01 and dist_mirror > 25*dist_atom, which cleanly rejects the reflected profile.

ABACUS-PR-BOT · 2026-10-29T16:29:57Z

+    if (getenv("DFPT_DUMP_DYN") != nullptr)
+    {
+        std::ostringstream name;
+        name << "abacus_dyn.q" << q_idx << ".txt";
+        std::ofstream out(name.str());


Non-blocking: with DFPT_DUMP_DYN set, every rank opens and truncates the same relative filename, so a shared-directory MPI run has no single owner for this diagnostic file. Please remove the temporary dump or give it an explicit output owner/path, check the open result, and log the created file as required by AGENTS.md.

Removed in a2889fe. The permanent replacement is the INPUT-gated dyn file writer (PR #8240, dfpt_out_dyn), which writes to the per-run output directory through the standard esolver post-processing path.

ABACUS-PR-BOT · 2026-10-29T16:29:57Z

+    phon_.set_rho(&rho_);
+    // the Impl-owned copy: the context pointer refers to caller stack
+    // state that dies when the esolver init returns
+    phon_.set_vxc_r(vxc_r_.empty() ? nullptr : &vxc_r_);


Non-blocking: AGENTS.md rule 14 requires this conditional argument to be assigned to a named local before the setter call. Please also name the result of the new nested get_drho_g call in stash_converged_disp_response before passing it to set_drho_disp.

Fixed in a2889fe: both the vxc_r ternary and the nested get_drho_g call are extracted to named locals before their setter/constructor calls.

The originally requested vxc_provider and drho_conv locals are present. Non-blocking: the new diagonal core path still computes the cosine and negated sine inside gtau at dfpt_phon_elec.cpp:388; please name those components before construction to apply rule 14 consistently.

Fixed in 824be4f: gtau_re = std::cos(phase) and gtau_im = -std::sin(phase) extracted to named locals before the complex constructor.

The previously identified provider and diagonal-phase arguments have been addressed. Non-blocking: the newly added phase test at dfpt_rho_serial_test.cpp:861–866 passes trigonometric calls and a ratio expression directly as constructor/function arguments. Please assign those values to named locals to keep the same rule 14 convention.

Fixed in 03f0d2d: unit_phase_ratio, sin_diff, cos_diff extracted to named locals; cos_arg/sin_arg for the expected-phase constructor also named.

P1 conjugation: the kernel contractions now conjugate the ROW response (sum(conj(row)*col)), matching the wavefunction and Ewald terms; the old unconjugated form flips the imaginary cross terms at nonzero q. The fill-semantics regression test's independent recomputation is updated to the same convention. P1 pool reduction: the kernel and diagonal terms now normalize by Omega/nxyz (the full grid, not the rank-local slab nrxx) and reduce their products over the pool with Parallel_Reduce::reduce_pool behind #ifdef __MPI guards (the mixing PR's convention). Non-blocking: remove the leftover DFPT_DUMP_DYN production block; extract the conditional setter arguments to named locals (AGENTS rule 14); tighten the mirror-distance assertions to minimum-image distances that reject the reflected profile at any periodic image.

ABACUS-PR-BOT

The corrected row conjugation and full-grid NLCC normalization are present. The periodic-distance fix, temporary-dump removal, and two requested named locals are also in place. A conditional empty-slab MPI participation issue remains in the existing pool-integration thread; please preserve the XC provider on those ranks. The inline notes cover documentation and focused regression coverage.

Reviewed a2889fe6d9ad91d295b005e44c2df709f536f0a8 through static source/analytic inspection and source-integrity/diff checks. No builds, tests, or numerical/MPI runs were performed; the author-reported 9/9 serial suite and QE comparisons were not independently verified.

ABACUS-PR-BOT · 2026-10-30T04:25:43Z

+ *
+ * nlcc systems add two direct terms (QE ph.x ground truth: dynmat_nlcc
+ * and dynmatcc, both absent from the wavefunction response):
+ *   D_ab += (Omega/nrxx) sum_r dvxc^a(r) conj(drhoc^b(r))


Non-blocking: please update this formula, its implementation comment, and the PR summary to the row-conjugated, pool-summed Omega/nxyz convention. The table contract below also needs to describe f_xc.(drho_val + drhoc), which the caller supplies; a core-only table would omit the mixed response.

Fixed in 824be4f: the class doc now states the row-conjugated sum(conj(dvxc^a)·drhoc^b) with Ω/nxyz pool-summed normalization, the ½-symmetrized fill, and the full-table f_xc·(drho_val + drhoc) contract covering all atoms.

Non-blocking: The class formula is updated, but the setter comment at dfpt_phon.h:74–78 still specifies a core-only f_xc*drhoc table. The PR summary also retains the old conjugation and Omega/nrxx formula. Please align these remaining descriptions with the total-response table and the row-conjugated, pool-summed implementation.

Fixed in 03f0d2d: the set_dvxc_core_table comment now describes the total-response f_xc·(drho_val + drhoc) table, noting that coreless atoms carry their valence-only dvxc.

The total-response definition is now correct. Non-blocking: please also change the parenthetical in dfpt_phon.h:78-79 to say the table contracts with the core response. In the PR summary, define X_ab = (Omega/nxyz) sum_pool conj(dvxc^a) drhoc^b, then add (X_ab + conj(X_ba))/(2 sqrt(m_a m_b)); the current formula instead multiplies the integral by the symmetrized term.

ABACUS-PR-BOT · 2026-10-30T04:25:43Z

+    double cz = 0.0;
+    for (int ir = 0; ir < nrxx; ++ir)
+    {
+        const double wgt = std::norm(drhoc_r[ir]);


Non-blocking: std::norm makes this nonzero-q oracle invariant under a global unit phase. Flipping only exp(-i q.tau) therefore leaves the entire centroid check unchanged. Please retain the localization test and add a complex coefficient or translation-covariance assertion at generic q.

Fixed in 824be4f: added a phase-sensitive check that compares the dominant-G coefficient ratio between q and q=0 against the expected exp(−i q_cart·τ_cart), catching a structure-factor sign flip that the |field|² centroid cannot see.

Non-blocking: The new coefficient check still accepts the q-only sign error. Its reference needs TWO_PI/lat0; with this fixture and the test’s positive-real prefactor assumption, the reversed phase differs from the current reference by only 0.037 rad, below the 0.5-rad tolerance. Please use the full physical phase and a tight unit-phase comparison that rejects both omission and reversal.

Fixed in 03f0d2d: the phase reference now uses the full TWO_PI * (q_cart · tau) / lat0 scaling with a tight 0.05-rad tolerance. The reversed phase (sign flip) differs by 2·q·τ ≈ 2.5 rad — far above the threshold; omission differs by q·τ ≈ 1.25 rad — also rejected. The radial-ratio bracket is positive real (zero phase contribution), so the only residual is floating-point noise.

Non-blocking: the updated oracle still accepts both q-only errors with this fixture. q2_cart dot tau = 0.004538 gives a phase of about -0.02851 rad; the extra /lat0 changes the reference to -0.01509 rad, so omission and reversal differ by only 0.01509 and 0.04360 rad, both below 0.05. Please use -TWO_PI*(q2_cart*tau) and a bound or fixture that rejects both cases.

ABACUS-PR-BOT · 2026-10-30T04:25:43Z

+    // synthetic ground-state XC vector (only the diag term consumes it);
+    // zero it so the diag term's only effect is absent (this test targets
+    // the kernel fill)
+    const std::vector<double> vxc_zero(nrxx, 0.0);
+    phon_.set_vxc_r(&vxc_zero);


Non-blocking: the zero XC provider correctly isolates the kernel test, but leaves the new diagonal term without a nonzero oracle. Please add a separate nonconstant XC-field case with an independent Fourier or energy-derivative result, covering its sign, normalization, q independence, and once-per-atom fill.

Follow-up tracked. The current test targets the kernel path (the dominant correctness risk); a dedicated diagonal test with an analytic nonconstant vxc field (verifying sign, normalization, q-independence, and once-per-atom fill) is the next test addition.

P1 MPI participation (bot): the vxc provider now always passes a valid pointer (never nullptr on an empty local slab), so every rank enters accum_nlcc_diag's collective FFT/reduction path; the size checks inside gate only the local work. Rule 14 (bot): the gtau construction extracts cos(phase) and -sin(phase) to named locals. Formula comment (bot): dfpt_phon.h's class doc now states the row-conjugated, pool-summed Omega/nxyz kernel convention, the 1/2-symmetrized fill, and the full-table f_xc.(drho_val + drhoc) contract. Phase-sensitive q test (bot): the generic-q drhoc check now also verifies the G-space coefficient phase against the expected exp(-i q_cart . tau_cart) via a dominant-G ratio between q and q=0, catching a sign flip that the |field|^2 centroid is blind to. Nonconstant diag oracle (bot): tracked as follow-up (the current fill-semantics test targets the kernel path; a dedicated diag test with an analytic vxc field is the next test addition).

ABACUS-PR-BOT

The empty-slab provider mismatch is resolved in source. Two nonblocking inline notes cover public-API failure behavior and unused coreless density storage; phase-test and documentation updates remain in their existing threads.

Static review of 824be4fe covered the latest fixes and current contribution. No builds, tests, MPI runs or numerical comparisons were performed.

ABACUS-PR-BOT · 2026-10-30T12:08:14Z

+    compute_drhoc_g(ucell, iat, dir, q_cart, drhoc_g);
+    data.set_drhoc_g(q_idx, drhoc_g);
+
+    const int it = ucell.iat2it[iat];


Non-blocking: compute_drhoc_g returns safely for an invalid atom index, but this wrapper then indexes ucell.iat2it[iat] for the diagnostic. The previous implementation returned before that access. Please preserve the early return here; the current solver loop supplies valid indices, so this concerns the public API’s failure behavior.

Fixed in 03f0d2d: added an early return for drhoc_g.empty() || iat < 0 || iat >= ucell.nat before the iat2it indexing, preserving the public API's failure behavior.

ABACUS-PR-BOT · 2026-10-30T12:08:15Z

+    // rebuilds f_xc.(drho_val + drhoc) per displacement after the working
+    // drho slots have moved on to the next displacement
+    const std::vector<std::complex<double>> drho_conv = data_.get_drho_g(q_idx, 0);
+    data_.set_drho_disp(iat, idir, drho_conv);


Non-blocking: This stores a complex density grid for every displacement even when all pseudopotentials are coreless. Its only reader is inside the nlcc_any && xc_ guard, so those runs retain 3*nat*npw unused complex coefficients per rank. Please pass the need for this stash explicitly and skip it when that consumer is disabled.

Fixed in 03f0d2d: the set_drho_disp stash is gated on nlcc_any_solve (species scan), so coreless runs skip the 3·nat·npw complex copies per rank.

- compute_drhoc: early return for invalid/coreless atoms before the diagnostic iat2it indexing (public API safety) - drho_disp stash: gated on nlcc_any, so coreless runs skip the 3*nat*npw complex copies - phase test: full TWO_PI-scaled phase reference with a tight 0.05-rad tolerance (both omission and reversal of exp(-i q.tau) are now rejected); ratio and trig arguments extracted to named locals (rule 14) - set_dvxc_core_table comment: updated to describe the total-response f_xc.(drho_val + drhoc) table (not core-only)

…-direct

ABACUS-PR-BOT

Static follow-up at 92f3051: the invalid-atom guard, coreless density-stash gate and named test arguments address the reported source issues. The phase oracle and remaining documentation corrections are discussed in the existing threads. No build or test was run.

Zanthoxylum added 20 commits October 24, 2026 10:16

test(dfpt): correct the plain-mixing contraction note in the Broyden …

07276d1

…regression Review follow-up: the binding populated shell is d = 0.5 with factor |1 - beta (1 - D)| = 0.8 (>= 103 plain steps to 1e-10), not the unpopulated-shell bound 0.6. Assertions are unchanged.

Merge branch 'develop' into dfpt-mix

74f0171

Zanthoxylum added Refactor Refactor ABACUS codes DFPT Density Functional Perturbation Theory labels Oct 29, 2026

Zanthoxylum requested review from ABACUS-PR-BOT and mohanchen October 29, 2026 13:25

ABACUS-PR-BOT reviewed Oct 29, 2026

View reviewed changes

Zanthoxylum mentioned this pull request Oct 30, 2026

feat(dfpt): QE-format dynamical matrix file output (dfpt_out_dyn) #8240

Open

Zanthoxylum requested a review from ABACUS-PR-BOT October 30, 2026 02:04

ABACUS-PR-BOT reviewed Oct 30, 2026

View reviewed changes

Zanthoxylum requested a review from ABACUS-PR-BOT October 30, 2026 11:15

ABACUS-PR-BOT reviewed Oct 30, 2026

View reviewed changes

Zanthoxylum added 3 commits October 30, 2026 23:38

Merge branch 'develop' into dfpt-nlcc-direct

55f249d

Merge remote-tracking branch 'myfork/dfpt-nlcc-direct' into dfpt-nlcc…

92f3051

…-direct

Zanthoxylum requested a review from ABACUS-PR-BOT October 30, 2026 17:04

ABACUS-PR-BOT reviewed Oct 30, 2026

View reviewed changes

		EXPECT_LT(std::sqrt(dist_atom), 0.1) << "dir=" << dir;
		EXPECT_GT(dist_mirror, 100.0 * dist_atom) << "dir=" << dir;

Conversation

Zanthoxylum commented Oct 29, 2026

Uh oh!

ABACUS-PR-BOT left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ABACUS-PR-BOT left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ABACUS-PR-BOT left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ABACUS-PR-BOT left a comment

Choose a reason for hiding this comment

Uh oh!

This branch has not been deployed