A beta version of the DFPT program - #7875
Open
Zanthoxylum wants to merge 40 commits into
Open
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
added 10 commits
August 20, 2026 20:06
… screened form - solve_efield_resp is now production (QE solve_e order): runs after solve_pos_resp, before the displacement solves; converged dpsi^E,a stashed through DFPT_PW_Data (dpsi_efield) - compute_eps consumes pos_resp + dpsi_efield: eps = 1 - (16 pi/Omega) sum_k wg sum_occ Re<Y^a|dpsi^E,b>, star-rotated on symmetry-reduced meshes; the PT r-matrix path is retired (pos_matrix kept as the design-phase analytic reference for its serial tests) - serial test ComputeEpsScfSyntheticStash replaces the PT two-level case (prefactor, wg, occupied sum, conj/index pinning, empty-row skip); 6/6 - end-to-end sym 4x4x4: eps = 23.35 delta (was IPA 12.67), consistent with the nosym ALEG value 23.68 and QE dielec.f90 anchor 23.67
…LAN P0-3 intake (non-Gamma-q chain defect, eps SCF promotion record)
…, KQ dual-reservoir completeness, term3 d2 ungating - compute_drho: include the spin factor 2 at every q (QE incdrhoscf wgt = 2*weight/omega); the q=0 Hermitian completion now keeps Re only instead of 2 Re. Previously the screening was half strength away from Gamma, which collapsed the L-point Si frequencies to -948/-148/182/199 cm^-1. After the fix: 100.49/100.49/380.41/402.11/485.93/485.93 cm^-1 vs QE 101.61x2/380.54/402.24/486.28x2 (Si NC 4x4x4, 0.1-1.1%); Gamma stays 517.491 cm^-1 (QE 517.633). - dfpt_kq_basis: dual-reservoir G assembly so the k and k+q balls share the same igl2ig maps (fixes silent truncation when one ball exhausts the rho-grid reservoir). - dfpt_phon: drop the 2q-reciprocal gate on the same-atom d2 term (it is q-independent by construction; the old gate silently dropped it and produced imaginary branches). - Verification: ctest 12/12 (MODULE_CELL x4 + MODULE_DFPT x8); serial 4/4 (pert/phon/q0/rho); bare-response L run matches QE niter_ph=1 to 0.008-0.4% (-2281.83 vs -2282.01 etc.). - No docs change: module_dfpt is design-phase, no INPUT parameter touched.
… frequencies, data-layer loto direction)
- DFPT_PW_Data: loto_dir_ (unit-normalized setter, isotropic (1,1,1)/sqrt(3)
default) and phon_freq_loto_ storage.
- DFPT_Phon: diagonalize_loto re-diagonalizes the Gamma matrix after
add_loto and stores signed frequencies separately (plain phon_freq(0)
stays intact); format_q_report/format_loto_report provide deterministic
fixed-precision blocks (header with direct q coordinates and the
correction direction).
- DFPT_PW::run uses data_.get_loto_dir() instead of the hardcoded
(1,1,1)/sqrt(3); new accessors get_nq/get_qvec/get_loto_dir/
get_phon_freq_loto/set_loto_dir plus the format forwarders.
- esolver run_post_process prints one block per q of the list plus the
LO-TO Gamma block when enabled; tensor blocks only print when computed.
- Serial regression: 3 new cases (direction normalization, closed-form
LO-TO spectrum {0, 13/12*pref}, char-exact format strings); phon 12/12,
ctest 12/12, all 4 DFPT serial tests pass.
- End-to-end smoke (Gamma, compute_q0+loto, 4x4x4): TO 517.490709
unchanged, LO-TO block along (0.577350 0.577350 0.577350), eps_inf
23.6825 and Z*=-1.19928d for both atoms vs QE 23.6685/-1.19765 (0.13%).
QE itself prints same-sign Z* with asr Sum=-2.395 for this setup; the
acoustic-branch lift is the faithful consequence, not a defect.
- No docs change: module_dfpt is design-phase, no INPUT parameter touched.
- DFPT_Rho::init gains mix_type (plain/kerker) and kerker_a2 (1/lat0^2); no charge_mixing.h dependency, screen f_g = |G+q|^2/(|G+q|^2+a^2) built with the v_hartree_q convention (gcar + q_frac*G). Screen both inputs, plain_mix, add the screened part back: mixed = rin + beta*f*(out-rin) (QE semantics, stored density stays physical; |G+q|=0 harmonic frozen, consistent with its drop in compute_drho). Init signature extended with an explicit kerker_a2 argument (no default arg; both call sites updated). - Wiring: env DFPT_MIX_TYPE / DFPT_KERKER_A2 design-phase knobs mirroring the DFPT_MIX_BETA precedent; default plain keeps behavior identical and the beta=0.4 default (and its stability rationale) stays documented in the init comment. No INPUT parameter change: no docs update required (env knobs are internal calibration aids, same category as DFPT_MIX_BETA). - Tests (dfpt_rho_serial, 6 -> 8): analytic first Kerker step; lambda=-2.2 stiff-shell model problem where plain beta=0.7 diverges (residual > 1) and kerker converges (< 1e-8) to the target. - Fixed latent breaks masked by a stale test binary since a915352: kq0.init not updated to the 4-arg DFPT_KQ_Basis::init signature, and the brute-force references missing the band-weight spin factor 2. - End-to-end (L point, 4x4x4, abacus_pw_para v3.11.0-beta8): plain beta=0.7 diverges (|drho| -> 1e20); kerker beta=0.7 converges in 1393 s (vs 2332 s plain beta=0.4); frequencies identical across plain 0.4 / kerker 0.4 / kerker 0.7 to 8-9 digits (100.487828 x2 / 380.41385 / 402.10912 / 485.93199 x2 cm^-1). - Verification: OMP_NUM_THREADS=1 ctest -R 'MODULE_CELL_klist_test$| MODULE_CELL_reciprocal_grid_test|MODULE_CELL_qlist_test| MODULE_CELL_little_group_test|MODULE_DFPT' -> 12/12; serial suites pert 8 / phon 12 / q0 6 / rho 8 all pass; governance --staged clean except the expected no-docs-needed WARNING recorded here.
…DFPT_IrrepData adapter - DFPT_PW_Data: the write-only single-slot ledger (set_current_iter(int)/ set_converged(bool)/add_residual(double)) is replaced by the (q,irrep)-keyed six-accessor ledger sunk from DFPT_IrrepData (std::map value members, missing keys read as not-converged / empty history / iteration 0, clean() drops the ledger). The irrep dimension stays as the stage-A slot: the fallback irrep 0 carries the full 3N displacement basis. The new <map>/ <utility> includes are required by the map value members the header owns. - DFPT_IrrepData adapter deleted (git rm): its irrep==0 forwarding of dpsi/drho/dv duplicated the existing per-q data API, and its own keyed maps moved to the data layer. get_dpsi_obj (static dummy, zero callers) removed. Both CMakeLists updated, including the pw_run_test source list. - run() outer-while accounting made honest: current_iter now increments per pass and convergence is worst-final-displacement-residual < conv_thr instead of an unconditional single pass. An unconverged pass re-runs the full solve (solve_displacement restarts from a zero input), bounded by max_iter outer passes, with the residual history keeping a record. Behavior on converged runs is bit-identical. - solve_displacement / solve_efield_resp: write-only inner ledger writes removed; per-displacement state stays local and the final residual returns to run() for aggregation. - Tests: dfpt_irrep_data_test.cpp renamed/rewritten as dfpt_pw_data_test.cpp (target MODULE_DFPT_pw_data_test, 5 cases: QList delegation, bound-safe accessors with the (q,spin) signature, setter round trip, keyed-ledger independence + clean() reset, U0 reservation). - Verification: OMP_NUM_THREADS=1 ctest -R 'MODULE_CELL_klist_test$| MODULE_CELL_reciprocal_grid_test|MODULE_CELL_qlist_test| MODULE_CELL_little_group_test|MODULE_DFPT' -> 12/12 (pw_data_test fills the retired irrep_data_test slot); serial suites pert 8 / phon 12 / q0 6 / rho 8 all pass; end-to-end L-point default-config smoke (abacus_pw_para v3.11.0-beta8) reproduces the reference frequencies bit-consistently (100.487828/100.487829/380.413847/402.109158/ 485.931988/485.931988 cm^-1, TOTAL 2332 s, same as the pre-B4 reference). Governance --staged: header-include warning justified by map value members; no INPUT behavior change so no docs update required.
- Deleted (acceptance complete): PTCHK gauge/term2/HF-channel probes and the drho_dfpt.dat dump; the DYNCHK family (term2/d2gate/d2k/d2/ion/ele/elei and the DYNCHK4 double-zheev comparison); MDBG binary dumps (x2); JPROBE + JPROBE_NOXC (B3 acceptance done, delete as planned); OCCCHK incl. the dbg_miss label analysis and the empty_kq_/empty_kq_eig_ companion storage; XB; BPT incl. the want_empty projector expansion; NOSC; XCS/NOXC (v_sc assembly simplified to the knob-free path); DKCHK; YCHK; D2MID (include_middle sunk to literal true, q-independence settled); ALEG + PTCROSS (the whole aleg_crosscheck method); STARDBG; Q0DBG. Dead accumulators (d2sum_loc/nl, cross_k) and the now-purposeless <fstream>/ <set> includes removed with them. - Kept: DFPT_DEBUG (SCF residual tracing + posresp tracking, the B3/B4 acceptance instrument and routine convergence diagnostics) and the B3 calibration knobs DFPT_MIX_BETA / DFPT_MIX_TYPE / DFPT_KERKER_A2 (documented in the DFPT_Rho::init comment). - Behavior-preserving: every deleted probe was env-gated off by default; include_middle and want_empty defaults equal the sunk values. - Verification: OMP_NUM_THREADS=1 ctest -R 'MODULE_CELL_klist_test$| MODULE_CELL_reciprocal_grid_test|MODULE_CELL_qlist_test| MODULE_CELL_little_group_test|MODULE_DFPT' -> 12/12; serial suites pert 8 / phon 12 / q0 6 / rho 8 pass; end-to-end L-point default-config smoke (abacus_pw_para v3.11.0-beta8) reproduces the reference frequencies bit-consistently (100.487828/100.487829/380.413847/402.109158/ 485.931988/485.931988 cm^-1). Governance --staged clean except the expected no-docs-needed WARNING (internal env probes, no INPUT change).
mohanchen
self-requested a review
August 28, 2026 13:56
mohanchen
added
Features Needed
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
2 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.
Reminder
AGENTS.mdanddocs/developers_guide/agent_governance.md.source/changes.Linked Issue
Fix #
Unit Tests and/or Case Tests for my changes
What's changed?
Governance Notes