Skip to content

Commit 94c4bc5

Browse files
ZanthoxylumZanthoxylum
andauthored
Initial implementation of a density-functional perturbation theory code using a plane-wave basis set (deepmodeling#7875)
* Refactor: extract ModuleCell::ReciprocalGrid base for k/q grids (Phase 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. * Refactor: QList on ReciprocalGrid base with star reduction + tests (Phase 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. * Feat: LittleGroup interface for q-point irreps + wire QList (Phase 3) - 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. * Feat: DFPT per-irrep SCF loop via DFPT_IrrepData adapter + tests (Phase 4) * Feat: complete QList q-point management (Cartesian, file read, print, use_irreps) * Feat: reserve DFT+U interface for DFPT (U0) 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). * Feat: k+q plane-wave basis enumeration for DFPT (C0) 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). * Feat: first-order perturbation potentials for DFPT (C1) 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. * Feat: projected CG Sternheimer solver for DFPT (C2) 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. * Feat: first-order density response for DFPT (C3) 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. * Feat: dynamical matrix for DFPT (C4 guard + C5) - 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 * Feat: q->0 response for DFPT (C6: eps, Born, v_hartree_q, XC contract) - 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 * Feat: wire DFPT driver and esolver factory (C7) 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. * Fix: DFPT screening-channel calibration (q=0 completion, XC central difference, 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. * Fix: DFPT plain-mixing default beta 0.7 -> 0.4 (small-G Coulomb stiffness) 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. * Docs: record DFPT stage-B gap audit and revised execution plan * Feat: INPUT-driven DFPT parameters (dfpt_qmesh/qfile/compute_q0/loto/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 * Test: sync DFPT serial references to production conventions 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. * Fix: multi-k DFPT ball-label matching and smeared-occupation projector 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. * Docs: record multi-k DFPT root causes and FD validation matrix in PLAN * Fix: reject metallic smearing occupations in DFPT with an explicit guard 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. * Docs: record validation-ladder extension, metallic-regime boundary, MPI smoke in PLAN * Docs: record non-Gamma q smoke results (dfpt_qfile end-to-end, q<->-q consistency) * Fix: drop spurious 1/nk in DFPT eps/born sums; wg already carries full 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. * Docs: record continuation plan (P0-1 uncommitted-fix intake, P0-2 Zstar bug, P0-3 B0 closeout, B2-B4, cleanup, A) * Fix: gate same-atom d2V_ext on 2q reciprocal; drop spurious ion_ion delta/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. * Debug: DFPT design-phase probes (ZDBG/BPT/NOSC/D2MID/DYNCHK/XB) 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. * Docs: P0-1 done (2q-reciprocal d2 gate intake, order-dependence fix, 28/28 serial) * DFPT q0: star-rotate the symmetry-reduced eps/Z* tensor sums 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. * DFPT q0: Sternheimer screened Z* (v4), QE-anchored eps 16pi fix, zstar_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 * DFPT q0: promote the E-field SCF solve, compute_eps to the dielec.f90 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 * DFPT: build_occ_kq diagnostic detail in the commensurability error; PLAN P0-3 intake (non-Gamma-q chain defect, eps SCF promotion record) * DFPT: fix q!=Gamma phonon frequencies (missing spin factor 2 in drho), 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. * DFPT PLAN: P0-3 non-Gamma-q defect root-caused and fixed (drho spin factor 2, a915352) * DFPT B2: formalize the phonon output (multi-q report, LO-TO corrected 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 B3: Kerker-preconditioned density mixing in DFPT_Rho - 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 B4: sink the (q,irrep) SCF ledger into DFPT_PW_Data, retire the 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. * DFPT: retire the B-phase validation instrumentation (net -977 lines) - 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). * delete PLAN * Fix: adapt DFPT to the refactored Plus_U interface (compile break + U guard) The develop-side DFT+U refactor (deepmodeling#7852-deepmodeling#7867) removed source_lcao/module_dftu/dftu.h and the is_locale_initialized() member, which broke every CMake build configuration of this branch at dfpt_pw_data.cpp (all 9 CI build variants plus Test/CUDA/abacuslite failed at the compile step; only the Makefile job passed because the Makefile.Objects DFPT entries were absent at that merge point). Changes: - DFPT now consumes the PW-side Plus_U_Base (source_pw/module_pwdft/ dftu_base.h) instead of the LCAO-side Plus_U header: dftu_ member, DFPT_PW_Data::init / DFPT_PW::init signatures and get_dftu() all use const Plus_U_Base* (the esolver call site passes &this->dftu with an implicit upcast). This removes the PW -> LCAO cross-layer include. - u_active() = with_u() && is_occ_mat_initialized(): the reservation usability now follows the occupation-matrix state of the provider. - DFPT_PW::init rejects a wired provider explicitly (WARNING_QUIT): the ground state supports PW-basis DFT+U now, but every DFPT U hook (cal_docc, build_dv_u, dftu_onsite, born/docc contractions) is a no-op U0 reservation, so running anyway would silently drop the whole first-order U response (fail-loud, same pattern as the metallic-sampling guard). - test/dftu_test_support.cpp rewritten: the old static-member replicas no longer exist; the shim now provides only the Plus_U_Base ctor/dtor (also linked into MODULE_DFPT_pw_data_test, which constructs the provider directly). dfpt_pw_run_test's locale test becomes a death test pinning the WARNING_QUIT guard; the with_u/u_active contract moved to DFPT_PW_DataTest.DftuReservationProviderUsability; the unused dftu.h includes dropped from the phon/q0 serial tests. Verification (GNU 8.3.1 + OpenMPI 5.0.3, GCC13 no-MPI cross-check): - cmake --build build --target abacus_pw_para: builds/links - cmake --build build-nompi (-DENABLE_MPI=OFF -DENABLE_LCAO=OFF) --target abacus_pw_omp: builds/links - ctest -R 'MODULE_DFPT|MODULE_CELL': 50/50 pass (incl. the new death test and provider-usability case); ./build/abacus_pw_para --version prints v3.11.0-beta8 - agent_governance_check --staged: no findings * Fix: reduce PR global dependency budget to non-increasing The governance checker blocks the PR while the diff's added lines carry more GlobalV/GlobalC/PARAM references than the removed lines (added=51, removed=22, net_delta=+29 -> CI 'Governance checks' exit 1). 30 of the added references were test-side GlobalV::ofs_running streams; they now use local std::ofstream objects (the fixture members that already existed), and ReciprocalGrid::print_klists prints through its own ofs parameter instead of the global stream (its single caller passes the same running log). The stale 'Originally GlobalV::FINAL_SCF' comment wording is dropped. The remaining production-side references (reciprocal_grid.cpp k-point-file echo, klist.cpp MY_RANK guards) are line-for-line moves of the previous klist.cpp code, so the budget is now non-increasing (net_delta = -4). Verification: ctest -R 'MODULE_DFPT|MODULE_CELL' 50/50 pass; abacus_pw_para relinks; agent_governance_check --base origin/develop --head HEAD exits 0 (no BLOCK findings). * Fix: link K_Vectors/ReciprocalGrid sources into tests broken by the ReciprocalGrid refactor The ReciprocalGrid refactor (Phase 1-3 of this PR) made K_Vectors polymorphic: its vtable is now keyed on K_Vectors::renew and emitted in klist.cpp, and the base vtable lives in reciprocal_grid.cpp. Twelve test targets across estate/hsolver/stodft/io instantiate K_Vectors but never compiled those translation units, so they fail to link after the merge (masked until now by the earlier dftu.h compile break): - MODULE_ESTATE_elecstate_{print,base,pw,energy} - MODULE_PW_Sto_Hamilt_UTs - MODULE_HSOLVER_pw - MODULE_IO_write_bands (test_serial) - MODULE_IO_write_eig_occ_test / write_dos_pw / print_info / read_wf2rho_pw_test (already had klist.cpp, lacked reciprocal_grid.cpp) - MODULE_IO_write_dmk Mirrors the pattern already used by this PR's own klist/qlist tests: add klist.cpp + parallel_kpoints.cpp + k_vector_utils.cpp + reciprocal_grid.cpp to SOURCES and the symmetry lib to LIBS. Verified: full build green except MODULE_IO_numerical_basis_test (needs ENABLE_LCAO, unguarded on develop as well); the fixed tests pass under ctest; remaining local failures are environment artifacts (ScaLAPACK abort-stub, ELPA off). * Fix: link K_Vectors/ReciprocalGrid sources into LCAO-side tests and add new DFPT objects to Makefile.Objects The ReciprocalGrid refactor made K_Vectors polymorphic (its key function and the base vtable now live in klist.cpp / reciprocal_grid.cpp), so any test that instantiates K_Vectors (module_dm tests, deltaspin spin_constrain/template_helpers via spin_constrain.cpp, and init_dm_from_file via density_matrix_io.cpp) fails to link. Also register the five PR-added translation units (reciprocal_grid.cpp, little_group.cpp, read_inp_dfpt.cpp, dfpt_hamilt_shift.cpp, dfpt_kq_basis.cpp) in source/Makefile.Objects so the Intel Makefile build does not fail with undefined references. * Docs: resync parameters.yaml and input-main.md with the C++ Input_Item generator The DFPT parameter block was hand-placed at a position that differs from the item_dfpt() registration order, so the byte-exact consistency checks in test.yml (--generate-parameters-yaml / generate_input_main.py) fail. Regenerate both files with the documented commands to restore sync; the only change is the position of the DFPT category block. * Fix: compile reciprocal_grid.cpp in deepks unit tests The ReciprocalGrid refactor made K_Vectors derive from ModuleCell::ReciprocalGrid, so klist.cpp.o and k_vector_utils.cpp.o now reference ReciprocalGrid member functions and its vtable. The deepks_unit_support object library (DEEPKS_UNIT_COMMON_SOURCES, gated behind ENABLE_MLALGO and thus only compiled in the gnu Test CI job) compiles klist.cpp without reciprocal_grid.cpp, failing to link all 30 MODULE_LCAO_DEEPKS_* test executables with undefined references to ModuleCell::ReciprocalGrid::renew/Monkhorst_Pack/build_star_ops/... and its vtable/typeinfo. Add the missing translation unit to the common source set; the symmetry library (incl. little_group.cpp) is already on the link line. * Fix: use threadsafe death tests in DFPT suites to avoid fork-in-threaded-process deadlock MODULE_DFPT_pw_run_test timed out (1700 s) in the gnu Test CI job: the two irrep-loop tests run first execute OpenMP regions, so with the job's OMP_NUM_THREADS=2 the process is multithreaded when the third test (dftu-reservation EXPECT_EXIT) forks. The default fast-style child then deadlocks on exit and the parent waits forever (reproduced locally under OMP_NUM_THREADS=2: gtest warns 'detected 2 threads' and hangs). Switch all three DFPT death tests to the fork+exec threadsafe style (same pattern as module_container tensor_test). For the pw_run test also bridge std::cout to std::cerr inside the death statement: WARNING_QUIT prints the NOTICE block to stdout, while death tests match the child's stderr; the old CaptureStdout+HasSubstr assertion cannot see the re-exec child's output. Verified under OMP_NUM_THREADS=2: pw_run 3/3 in 0.3 s (previously indefinite hang), kq_basis 5/5, pert_serial 8/8, and the full MODULE_DFPT ctest batch 8/8. * Refactor DFPT unit tests: consolidate ctor/dtor stubs into shared dfpt_test_mocks.cpp (mirror tmp_mocks.cpp convention); absorb dftu_test_support.cpp * Refactor DFPT unit tests: share the cubic-cell/stru_lib fixture between pw_data and pw_run tests (dfpt_stru_fixture) * Refactor DFPT serial tests: derive pert/rho/phon/q0 fixtures from a shared DFPTSerialBase (cell/basis/data setup, Coulomb/NC atom builders, analytic dVloc reference) * test(dfpt): dedupe repeated analytic blocks in the phon serial test Share the occupied-weights table, the single-plane-wave psi builder, the analytic accumulate_electron cross term (now on top of AnalyticDVloc), and the isotropic loto data setup (eps/Born charges + two-atom mass table via MakeTwoAtomCell) through phon fixture helpers; the three AccumulateElectron tests and the two loto closed-form tests keep their reference formulas but drop the duplicated inline copies. --------- Co-authored-by: Zanthoxylum <chenshengjun@localhost.localdomain>
1 parent 3df3bb4 commit 94c4bc5

72 files changed

Lines changed: 12306 additions & 1307 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

docs/advanced/input_files/input-main.md

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -589,6 +589,14 @@
589589
- [Reduced Density Matrix Functional Theory](#reduced-density-matrix-functional-theory)
590590
- [rdmft](#rdmft)
591591
- [rdmft\_power\_alpha](#rdmft_power_alpha)
592+
- [Density functional perturbation theory](#density-functional-perturbation-theory)
593+
- [dfpt\_qmesh](#dfpt_qmesh)
594+
- [dfpt\_qfile](#dfpt_qfile)
595+
- [dfpt\_compute\_q0](#dfpt_compute_q0)
596+
- [dfpt\_loto](#dfpt_loto)
597+
- [dfpt\_conv\_thr](#dfpt_conv_thr)
598+
- [dfpt\_max\_iter](#dfpt_max_iter)
599+
- [dfpt\_mix\_beta](#dfpt_mix_beta)
592600

593601
## System variables
594602

@@ -643,6 +651,7 @@
643651
- nep: Neuroevolution Potential
644652
- ks-lr: Kohn-Sham density functional theory + LR-TDDFT (Under Development Feature)
645653
- lr: LR-TDDFT with given KS orbitals (Under Development Feature)
654+
- dfpt: density functional perturbation theory (Under Development Feature)
646655
- **Default**: ksdft
647656

648657
### symmetry
@@ -5192,3 +5201,49 @@
51925201
- **Default**: 0.656
51935202

51945203
[back to top](#full-list-of-input-keywords)
5204+
5205+
## Density functional perturbation theory
5206+
5207+
### dfpt_qmesh
5208+
5209+
- **Type**: Vector of Int (1 or 3 values)
5210+
- **Description**: Set the Monkhorst-Pack q mesh (gamma-centered) for DFPT phonon calculations. The q mesh must be commensurate with the ground-state k mesh: k + q must be a point of the k list (modulo a reciprocal lattice vector). For example, a 4x4x4 KPT mesh is commensurate with dfpt_qmesh values of 1, 2, or 4 along each direction. This parameter is ignored when dfpt_qfile is set.
5211+
- **Default**: 1 1 1
5212+
5213+
### dfpt_qfile
5214+
5215+
- **Type**: String
5216+
- **Description**: Set the file containing the q points for DFPT, in the same format as the KPT file (Q_POINTS card: Gamma/Monkhorst-Pack mesh, or an explicit Direct/Cartesian list; symmetry reduction is not applied to file q lists). When set, it overrides dfpt_qmesh. Each q point must still be commensurate with the ground-state k mesh.
5217+
- **Default**: ""
5218+
5219+
### dfpt_compute_q0
5220+
5221+
- **Type**: Boolean
5222+
- **Description**: Whether to compute the macroscopic dielectric tensor (epsilon_inf) and the Born effective charges at q = 0 within the same DFPT run. Requires a q point at Gamma (the default dfpt_qmesh 1 1 1).
5223+
- **Default**: false
5224+
5225+
### dfpt_loto
5226+
5227+
- **Type**: Boolean
5228+
- **Description**: Whether to apply the Lyddane-Sachs-Teller non-analytic correction to the Gamma-point dynamical matrix, which splits the longitudinal and transverse optical modes. Requires dfpt_compute_q0 to be true, since the correction is built from epsilon_inf and the Born effective charges.
5229+
- **Default**: false
5230+
5231+
### dfpt_conv_thr
5232+
5233+
- **Type**: Real
5234+
- **Description**: Set the convergence threshold of the self-consistent DFPT cycle: the iteration stops when the relative residual of the first-order density ||drho_out - drho_in|| / ||drho_out|| drops below this value for every displacement.
5235+
- **Default**: 1.0e-8
5236+
5237+
### dfpt_max_iter
5238+
5239+
- **Type**: Integer
5240+
- **Description**: Set the maximum number of self-consistent DFPT iterations for each atomic displacement.
5241+
- **Default**: 100
5242+
5243+
### dfpt_mix_beta
5244+
5245+
- **Type**: Real
5246+
- **Description**: Set the plain-mixing coefficient of the first-order density in the self-consistent DFPT cycle. The response Jacobian has strongly negative eigenvalues on the smallest-G shells (Coulomb stiffness), so beta must stay below 2 / (1 + |lambda_min|); the default 0.4 keeps margin up to |lambda_min| ~ 3. A larger value accelerates convergence for weakly screened systems but may diverge.
5247+
- **Default**: 0.4
5248+
5249+
[back to top](#full-list-of-input-keywords)

docs/parameters.yaml

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@ parameters:
6161
* nep: Neuroevolution Potential
6262
* ks-lr: Kohn-Sham density functional theory + LR-TDDFT (Under Development Feature)
6363
* lr: LR-TDDFT with given KS orbitals (Under Development Feature)
64+
* dfpt: density functional perturbation theory (Under Development Feature)
6465
default_value: ksdft
6566
unit: ""
6667
availability: ""
@@ -3046,6 +3047,62 @@ parameters:
30463047
default_value: "-1 2 -1 2"
30473048
unit: primitive cells
30483049
availability: "lr_solver==plot and exciton_plot_format in [slice, both]"
3050+
- name: dfpt_qmesh
3051+
category: Density functional perturbation theory
3052+
type: Vector of Int (1 or 3 values)
3053+
description: |
3054+
Set the Monkhorst-Pack q mesh (gamma-centered) for DFPT phonon calculations. The q mesh must be commensurate with the ground-state k mesh: k + q must be a point of the k list (modulo a reciprocal lattice vector). For example, a 4x4x4 KPT mesh is commensurate with dfpt_qmesh values of 1, 2, or 4 along each direction. This parameter is ignored when dfpt_qfile is set.
3055+
default_value: 1 1 1
3056+
unit: ""
3057+
availability: ""
3058+
- name: dfpt_qfile
3059+
category: Density functional perturbation theory
3060+
type: String
3061+
description: |
3062+
Set the file containing the q points for DFPT, in the same format as the KPT file (Q_POINTS card: Gamma/Monkhorst-Pack mesh, or an explicit Direct/Cartesian list; symmetry reduction is not applied to file q lists). When set, it overrides dfpt_qmesh. Each q point must still be commensurate with the ground-state k mesh.
3063+
default_value: "\"\""
3064+
unit: ""
3065+
availability: ""
3066+
- name: dfpt_compute_q0
3067+
category: Density functional perturbation theory
3068+
type: Boolean
3069+
description: |
3070+
Whether to compute the macroscopic dielectric tensor (epsilon_inf) and the Born effective charges at q = 0 within the same DFPT run. Requires a q point at Gamma (the default dfpt_qmesh 1 1 1).
3071+
default_value: "false"
3072+
unit: ""
3073+
availability: ""
3074+
- name: dfpt_loto
3075+
category: Density functional perturbation theory
3076+
type: Boolean
3077+
description: |
3078+
Whether to apply the Lyddane-Sachs-Teller non-analytic correction to the Gamma-point dynamical matrix, which splits the longitudinal and transverse optical modes. Requires dfpt_compute_q0 to be true, since the correction is built from epsilon_inf and the Born effective charges.
3079+
default_value: "false"
3080+
unit: ""
3081+
availability: ""
3082+
- name: dfpt_conv_thr
3083+
category: Density functional perturbation theory
3084+
type: Real
3085+
description: |
3086+
Set the convergence threshold of the self-consistent DFPT cycle: the iteration stops when the relative residual of the first-order density ||drho_out - drho_in|| / ||drho_out|| drops below this value for every displacement.
3087+
default_value: "1.0e-8"
3088+
unit: ""
3089+
availability: ""
3090+
- name: dfpt_max_iter
3091+
category: Density functional perturbation theory
3092+
type: Integer
3093+
description: |
3094+
Set the maximum number of self-consistent DFPT iterations for each atomic displacement.
3095+
default_value: "100"
3096+
unit: ""
3097+
availability: ""
3098+
- name: dfpt_mix_beta
3099+
category: Density functional perturbation theory
3100+
type: Real
3101+
description: |
3102+
Set the plain-mixing coefficient of the first-order density in the self-consistent DFPT cycle. The response Jacobian has strongly negative eigenvalues on the smallest-G shells (Coulomb stiffness), so beta must stay below 2 / (1 + |lambda_min|); the default 0.4 keeps margin up to |lambda_min| ~ 3. A larger value accelerates convergence for weakly screened systems but may diverge.
3103+
default_value: "0.4"
3104+
unit: ""
3105+
availability: ""
30493106
- name: out_freq_ion
30503107
category: Output information
30513108
type: Integer

source/Makefile.Objects

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -212,6 +212,7 @@ OBJS_CELL=atom_pseudo.o\
212212
setup_nonlocal.o\
213213
klist.o\
214214
k_vector_utils.o\
215+
reciprocal_grid.o\
215216
cell_index.o\
216217
cell_tools.o\
217218
check_atomic_stru.o\
@@ -374,6 +375,8 @@ OBJS_HAMILT_OF=kedf_tf.o\
374375
evolve_ofdft.o\
375376

376377
OBJS_DFPT=dfpt_metal.o\
378+
dfpt_hamilt_shift.o\
379+
dfpt_kq_basis.o\
377380
dfpt_pert.o\
378381
dfpt_phon.o\
379382
dfpt_pw.o\
@@ -564,6 +567,7 @@ OBJS_SYMMETRY=symm_other.o\
564567
symm_magnetic.o\
565568
symm_pricell.o\
566569
symm_rho.o\
570+
little_group.o\
567571
symmetry.o\
568572

569573
OBJS_XC=xc_functional.o\
@@ -674,6 +678,7 @@ OBJS_IO=module_parameter/input_conv.o\
674678
module_parameter/read_inp_model.o\
675679
module_parameter/read_inp_postproc.o\
676680
module_parameter/read_inp_exx_dftu.o\
681+
module_parameter/read_inp_dfpt.o\
677682
module_parameter/read_inp_other.o\
678683
module_parameter/read_inp_out.o\
679684
module_parameter/read_set_globalv.o\

source/source_cell/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ add_library(
2222
read_atoms_helper.cpp
2323
read_orb.cpp
2424
klist.cpp
25+
reciprocal_grid.cpp
2526
parallel_kpoints.cpp
2627
cell_index.cpp
2728
cell_tools.cpp

0 commit comments

Comments
 (0)