Repository navigation
feat(dfpt): complete the nlcc core-response terms (QE dynmat_nlcc + dynmatcc) and fix the mirrored structure-factor phase - #8237
Merged
Conversation
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.
Zanthoxylum
added
Refactor
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
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
There was a problem hiding this comment.
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.
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
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
There was a problem hiding this comment.
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.
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
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
There was a problem hiding this comment.
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.
- 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)
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
There was a problem hiding this comment.
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.
mohanchen
approved these changes
Oct 31, 2026
Zanthoxylum
added a commit
to Zanthoxylum/abacus-develop
that referenced
this pull request
Oct 31, 2026
…ut-dyn Conflict resolutions: keep both the myrank field (dyn output, guarded CommunicationDomain wrapper) and the vxc_r field (nlcc diagonal term) in DFPT_PW_InitContext; the wrapper passes myrank + nullptr for vxc_r; the esolver uses the full context with both fields.
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
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
The DFPT phonon implementation was missing the nonlinear core correction (nlcc) contribution to the force constants: the first-order core-density response
drhocentered 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.Gfiles:The two missing direct terms (QE
dynmat_nlcc+dynmatccequivalents):D_ab += (Ω/nrxx) Σ_r dvxc^a(r) conj(drhoc^b(r)) · (1/2)(X_ab + conj(X_ba)), withdvxc^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)D_aa[i][j] += ∫ dr w_ij(r) v_xc⁰(r)withw_ij = -tpiba² G_i G_j ρ̃_c(|G|) e^{-iGτ}, evaluated throughrecip2realas a real field;v_xc⁰computed in the esolver wiring and delivered viaDFPT_PW_InitContext::vxc_rStructure-factor phase fix (
compute_drhoc_g): the phase wasexp(+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 toexp(-i(G+q)τ)(the QEaddcoreconvention). 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 HermiticityAlso:
DFPT_Rho::compute_drhoc_gextracted as the pure core-response builder;assemble_v_sctakes 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:
bn_hex Γ optical modes (the C5 outlier, rms was 28.65 → now 1.20 cm⁻¹):
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.