From 1b31f7def941c8bebedd4086089978fd98d535f4 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sat, 17 Oct 2026 15:06:30 +0800 Subject: [PATCH 1/6] feat(dfpt): translational sum-rule audit of the Gamma dynamical matrix Print the per-atom row-sum norm of the Gamma dynamical matrix right after the raw diagonalization and, when the LO-TO correction runs, again on the corrected matrix (the nonanalytic term carries its own Z* sum rule, so the two violations are audited separately). This is the diagnostics feeder for the acoustic-branch investigation: the existing check_sum_rule helper only returns a bool against a 1e-6 relative tolerance and is not wired anywhere, while the observed Gamma acoustic modes are off by tens to hundreds of cm^-1 across all ten validation systems. --- source/source_pw/module_dfpt/dfpt_phon.cpp | 59 ++++++++++++++++++++ source/source_pw/module_dfpt/dfpt_phon.h | 7 +++ source/source_pw/module_dfpt/dfpt_pw_run.cpp | 4 ++ 3 files changed, 70 insertions(+) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 678f5071a5a..63921a43a54 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -163,6 +163,65 @@ void DFPT_Phon::diagonalize(int q_idx, DFPT_PW_Data& data) ModuleBase::timer::end("DFPT_Phon", "diagonalize"); } +void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Data& data) const +{ + ModuleBase::TITLE("DFPT_Phon", "audit_asr"); + ModuleBase::timer::start("DFPT_Phon", "audit_asr"); + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const bool at_gamma = std::abs(q_frac.x) <= gamma_tol && std::abs(q_frac.y) <= gamma_tol + && std::abs(q_frac.z) <= gamma_tol; + const int nat = ucell_->nat; + const int nat3 = 3 * nat; + const ModuleBase::ComplexMatrix dyn = data.get_dynmat(q_idx); + if (!at_gamma || dyn.nr != nat3) + { + ModuleBase::timer::end("DFPT_Phon", "audit_asr"); + return; + } + double max_elem = 0.0; + for (int i = 0; i < nat3; ++i) + { + for (int j = 0; j < nat3; ++j) + { + max_elem = std::max(max_elem, std::abs(dyn(i, j))); + } + } + // per-atom translational row sums s(ia)[da] = sum_ib,db D_{ia da, ib db} + std::cout << "DFPT ASR audit (" << label << "):"; + double worst = 0.0; + int worst_atom = -1; + for (int ia = 0; ia < nat; ++ia) + { + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + std::complex s(0.0, 0.0); + for (int jb = 0; jb < nat3; ++jb) + { + s += dyn(3 * ia + da, jb); + } + s2 += std::norm(s); + } + const double snorm = std::sqrt(s2); + std::cout << " atom" << ia + 1 << "(" << ucell_->atoms[ucell_->iat2it[ia]].label << ")=" + << std::setprecision(3) << std::scientific << snorm; + if (snorm > worst) + { + worst = snorm; + worst_atom = ia; + } + } + std::cout << std::defaultfloat << std::setprecision(6); + if (max_elem > 0.0 && worst_atom >= 0) + { + std::cout << " | worst rel = " << std::setprecision(3) << std::scientific << worst / max_elem + << std::defaultfloat << std::setprecision(6); + } + std::cout << std::endl; + ModuleBase::timer::end("DFPT_Phon", "audit_asr"); +} + void DFPT_Phon::add_loto(const ModuleBase::Vector3& qhat, DFPT_PW_Data& data) { ModuleBase::TITLE("DFPT_Phon", "add_loto"); diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 9e8c18cdceb..f6d8ce1e3ba 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -87,6 +87,13 @@ class DFPT_Phon /// from Gamma, where the rule does not apply). bool check_sum_rule(int q_idx, DFPT_PW_Data& data) const; + /// translational sum-rule audit of the Gamma dynamical matrix + /// (per-atom 3-vector row sums, relative to the matrix scale); + /// diagnostics for the acoustic-branch investigation: the raw matrix + /// and the LO-TO-corrected one are audited separately because the + /// nonanalytic term carries its own Z* sum rule + void audit_asr(const std::string& label, int q_idx, const DFPT_PW_Data& data) const; + /// Ewald ion-ion force constants C^ewald_ab(q) (G-space + real-space + /// Gaussian self term), mass-reduced by 1/sqrt(M_a M_b). A stateless /// building block exposed publicly so the serial analytic tests can diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 16aec15df2f..57a3c06b657 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -109,6 +109,9 @@ void DFPT_PW::Impl::run_assemble(int q_idx) { phon_.assemble(q_idx, data_); phon_.diagonalize(q_idx, data_); + // translational sum-rule audit of the raw Gamma matrix (acoustic-branch + // diagnostics; the loto-corrected matrix is audited separately below) + phon_.audit_asr("raw", q_idx, data_); if (q_idx != 0 || !data_.get_loto()) { return; @@ -118,6 +121,7 @@ void DFPT_PW::Impl::run_assemble(int q_idx) // set_loto_dir overrides, e.g. per irrep direction in stage A) phon_.add_loto(data_.get_loto_dir(), data_); phon_.diagonalize_loto(data_); + phon_.audit_asr("loto", q_idx, data_); } void DFPT_PW::run() From 99ecb2895763c3535fc7aed65c8334f76d6194c1 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 19 Oct 2026 12:52:51 +0800 Subject: [PATCH 2/6] feat(dfpt): acoustic sum rule enforcement on the Gamma dynamical matrix New INPUT parameter dfpt_asr = no (default) | simple | crystal, applied to the RAW Gamma dynamical matrix before the LO-TO nonanalytic term, together with the Z* sum rule sum_a Z*_a = 0 on the Born effective charges (the LO-TO term then consumes the corrected charges, QE dynmat set_asr order). Unconstrained linear-response self-consistency violates both sum rules at the percent to tens-of-percent level -- measured on ten validation systems (raw worst rowsum rel 1.8e-3 .. 8.3e-1, invariant under a 100x tighter convergence threshold), and QE ph.x raw output shows the same behavior (0.09% .. 34%); the reference workflow post-processes with dynmat asr='crystal', so this option aligns the treatment. - simple: the legacy fix, each row's total is subtracted from its diagonal element, with a Hermitian re-symmetrization pass - crystal: the optimized minimal-norm fix -- the symmetric projection D -> (1-P) D (1-P) with the orthonormal translations u_i(3a+alpha) = delta(alpha,i)/sqrt(nat); zeroes every row AND column sum exactly, preserves Hermiticity and keeps symmetry-equivalent atoms equivalent; the three acoustic eigenvalues vanish to machine precision - an explicit report line prints the worst rowsum before -> after (and the |sum Z*| correction); run_assemble audits raw / asr / loto stages separately Tests (MODULE_DFPT_phon_serial 12 -> 16): simple and crystal zero the row sums; crystal keeps Hermiticity, maps the acoustic triplet to zero (sqrt-scaled residual ~1e-5 cm^-1 from a ~1e-17 matrix residual), keeps equivalent atoms equivalent, and shifts each Born charge by the per-component mean. --- source/source_esolver/esolver_dfpt_pw.cpp | 1 + .../module_parameter/input_parameter.h | 1 + .../module_parameter/read_inp_dfpt.cpp | 36 ++++ source/source_pw/module_dfpt/dfpt_phon.cpp | 183 ++++++++++++++++++ source/source_pw/module_dfpt/dfpt_phon.h | 14 ++ source/source_pw/module_dfpt/dfpt_pw.cpp | 8 + source/source_pw/module_dfpt/dfpt_pw.h | 4 + source/source_pw/module_dfpt/dfpt_pw_impl.h | 1 + source/source_pw/module_dfpt/dfpt_pw_run.cpp | 9 + .../test_serial/dfpt_phon_serial_test.cpp | 150 ++++++++++++++ 10 files changed, 407 insertions(+) diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 79637856b42..a9c1bfca8c6 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -223,6 +223,7 @@ void ESolver_DFPT_PW::before_all_runners(BaseCell& basecell, const Input_para& i dfpt_->set_conv_thr(inp.dfpt_conv_thr); dfpt_->set_max_iter(inp.dfpt_max_iter); dfpt_->set_mix_beta(inp.dfpt_mix_beta); + dfpt_->set_asr(inp.dfpt_asr); dfpt_->set_compute_q0(inp.dfpt_compute_q0); dfpt_->set_loto(inp.dfpt_loto); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 5fba78e9b2b..3edfddb418f 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -439,6 +439,7 @@ struct Input_para double dfpt_conv_thr = 1.0e-8; ///< convergence threshold of the DFPT first-order density int dfpt_max_iter = 100; ///< max iterations of the DFPT first-order density mixing double dfpt_mix_beta = 0.4; ///< mixing coefficient of the DFPT first-order density + std::string dfpt_asr = "no"; ///< acoustic sum rule on the Gamma dynamical matrix: no | simple | crystal // ============== #Parameters (11.Output) =========================== int out_stru = 1; ///< output stru file each ion step diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index f4eceef9d56..3c9052bbbb6 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -168,5 +168,41 @@ void ReadInput::item_dfpt() }; this->add_item(item); } + { + Input_Item item("dfpt_asr"); + item.annotation = "acoustic sum rule on the Gamma dynamical matrix: no, simple or crystal"; + item.category = category; + item.type = "String"; + item.description = "Enforce the translational acoustic sum rule (sum over atoms of the " + "Gamma force-constant rows = 0) on the RAW Gamma dynamical matrix, " + "before the LO-TO nonanalytic term, together with the Z* sum rule " + "sum_a Z*_a = 0 on the Born effective charges. Unconstrained " + "linear-response self-consistency violates both at the percent to " + "tens-of-percent level (QE ph.x shows the same), and the reference " + "workflow post-processes with asr='crystal'. simple: subtract each " + "row sum from its diagonal element. crystal: the optimized " + "minimal-norm symmetric projection (keeps Hermiticity and " + "symmetry-equivalent atoms equivalent). no: leave the matrix as " + "computed."; + item.default_value = "no"; + item.read_value = [](const Input_Item& item, Parameter& para) { + if (item.get_size() == 0) + { + para.input.dfpt_asr = "no"; + } + else + { + para.input.dfpt_asr = strvalue; + } + }; + sync_string(input.dfpt_asr); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.dfpt_asr != "no" && para.input.dfpt_asr != "simple" && para.input.dfpt_asr != "crystal") + { + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_asr must be no, simple or crystal."); + } + }; + this->add_item(item); + } } } // namespace ModuleIO diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 63921a43a54..a9861d98a5b 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -222,6 +222,189 @@ void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Dat ModuleBase::timer::end("DFPT_Phon", "audit_asr"); } +void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data) +{ + ModuleBase::TITLE("DFPT_Phon", "apply_asr"); + ModuleBase::timer::start("DFPT_Phon", "apply_asr"); + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const bool at_gamma = std::abs(q_frac.x) <= gamma_tol && std::abs(q_frac.y) <= gamma_tol + && std::abs(q_frac.z) <= gamma_tol; + const int nat = ucell_->nat; + const int nat3 = 3 * nat; + ModuleBase::ComplexMatrix dyn = data.get_dynmat(q_idx); + if (!at_gamma || dyn.nr != nat3) + { + std::cout << "DFPT ASR (" << asr_type << ") skipped: not a Gamma matrix" << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); + return; + } + + // violation BEFORE the correction (for the report line) + double worst_before = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + std::complex s(0.0, 0.0); + for (int jb = 0; jb < nat3; ++jb) + { + s += dyn(3 * ia + da, jb); + } + s2 += std::norm(s); + } + worst_before = std::max(worst_before, std::sqrt(s2)); + } + + if (asr_type == "simple") + { + // legacy fix: each row's total is subtracted from its diagonal + // element (QE dynmat 'simple'); Hermiticity is restored by the + // re-symmetrization pass below, mirroring QE's upper-triangle use + for (int i = 0; i < nat3; ++i) + { + std::complex s(0.0, 0.0); + for (int j = 0; j < nat3; ++j) + { + s += dyn(i, j); + } + dyn(i, i) -= s; + } + } + else if (asr_type == "crystal") + { + // minimal-norm symmetric projection D -> (1-P) D (1-P) with the + // orthonormal translations u_i(3a+alpha) = delta(alpha,i)/sqrt(nat): + // M(alpha, nb, beta) = (1/nat) sum_na' D(na',alpha; nb,beta) (P D) + // N(na, alpha, beta) = (1/nat) sum_nb' D(na,alpha; nb',beta) (D P) + // T(alpha, beta) = sum_na,nb D(na,alpha; nb,beta) + // D' = D - M - N + T/nat^2 + std::vector> m(3 * nat3, std::complex(0.0, 0.0)); + std::vector> n(3 * nat3, std::complex(0.0, 0.0)); + ModuleBase::ComplexMatrix t(3, 3, true); + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int nb = 0; nb < nat; ++nb) + { + for (int db = 0; db < 3; ++db) + { + const std::complex d = dyn(3 * na + da, 3 * nb + db); + m[da * nat3 + 3 * nb + db] += d / static_cast(nat); + n[(3 * na + da) * 3 + db] += d / static_cast(nat); + t(da, db) += d; + } + } + } + } + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int nb = 0; nb < nat; ++nb) + { + for (int db = 0; db < 3; ++db) + { + dyn(3 * na + da, 3 * nb + db) -= m[da * nat3 + 3 * nb + db]; + dyn(3 * na + da, 3 * nb + db) -= n[(3 * na + da) * 3 + db]; + dyn(3 * na + da, 3 * nb + db) += t(da, db) / static_cast(nat * nat); + } + } + } + } + } + else + { + std::cout << "DFPT ASR: unknown type '" << asr_type << "', nothing applied" << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); + return; + } + + // Hermitian re-symmetrization (rows are corrected independently) + for (int i = 0; i < nat3; ++i) + { + for (int j = i + 1; j < nat3; ++j) + { + const std::complex avg = 0.5 * (dyn(i, j) + std::conj(dyn(j, i))); + dyn(i, j) = avg; + dyn(j, i) = std::conj(avg); + } + } + + // Z* sum rule (both enforcing modes, QE set_asr convention): shift by + // the per-component mean over atoms, sum_a Z*_a -> 0 + double zworst_before = 0.0; + ModuleBase::ComplexMatrix zsum(3, 3, true); + int n_born = 0; + for (int ia = 0; ia < nat; ++ia) + { + const ModuleBase::matrix za = data.get_born(ia); + if (za.nr != 3 || za.nc != 3) + { + continue; + } + ++n_born; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + zsum(a, b) += std::complex(za(a, b), 0.0); + } + } + } + if (n_born == nat) + { + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + zworst_before = std::max(zworst_before, std::abs(zsum(a, b))); + } + } + for (int ia = 0; ia < nat; ++ia) + { + ModuleBase::matrix za = data.get_born(ia); + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + za(a, b) -= zsum(a, b).real() / nat; + } + } + data.set_born(ia, za); + } + } + + data.set_dynmat(q_idx, dyn); + + // violation AFTER the correction + double worst_after = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + std::complex s(0.0, 0.0); + for (int jb = 0; jb < nat3; ++jb) + { + s += dyn(3 * ia + da, jb); + } + s2 += std::norm(s); + } + worst_after = std::max(worst_after, std::sqrt(s2)); + } + std::cout << "DFPT ASR (" << asr_type << ") applied: worst rowsum " << std::setprecision(3) << std::scientific + << worst_before << " -> " << worst_after; + if (n_born == nat) + { + std::cout << " | max |sum_a Z*_a| " << zworst_before << " -> 0"; + } + std::cout << std::defaultfloat << std::setprecision(6) << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); +} + void DFPT_Phon::add_loto(const ModuleBase::Vector3& qhat, DFPT_PW_Data& data) { ModuleBase::TITLE("DFPT_Phon", "add_loto"); diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index f6d8ce1e3ba..0452805d41e 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -94,6 +94,20 @@ class DFPT_Phon /// nonanalytic term carries its own Z* sum rule void audit_asr(const std::string& label, int q_idx, const DFPT_PW_Data& data) const; + /// enforce the translational acoustic sum rule on the Gamma dynamical + /// matrix (and, in every enforcing mode, the Z* sum rule + /// sum_a Z*_a = 0 on the stored Born charges, like QE dynmat's + /// set_asr). "simple": the legacy fix, each row's total sum is + /// subtracted from its diagonal element (QE 'simple'). "crystal": + /// the optimized minimal-norm fix -- with the three orthonormal + /// translation vectors u_i(3a+alpha) = delta(alpha,i)/sqrt(nat), the + /// symmetric projection D -> (1-P) D (1-P), P = sum_i u_i u_i^T, + /// which zeroes every row AND column sum, preserves Hermiticity and + /// keeps symmetry-equivalent atoms equivalent (QE 'crystal' for the + /// translational part). Both apply to the RAW matrix, before the + /// LO-TO nonanalytic term; away from Gamma the call is a no-op. + void apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data); + /// Ewald ion-ion force constants C^ewald_ab(q) (G-space + real-space + /// Gaussian self term), mass-reduced by 1/sqrt(M_a M_b). A stateless /// building block exposed publicly so the serial analytic tests can diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 2d325c90059..de3b7a955a4 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -138,6 +138,14 @@ void DFPT_PW::set_mix_beta(double beta) } } +void DFPT_PW::set_asr(const std::string& asr_type) +{ + if (asr_type == "no" || asr_type == "simple" || asr_type == "crystal") + { + pimpl_->asr_type_ = asr_type; + } +} + void DFPT_PW::set_compute_q0(bool flag) { pimpl_->data_.set_compute_q0(flag); diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 83560769a19..ea611743b92 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -128,6 +128,10 @@ class DFPT_PW void set_mix_beta(double beta); + /// acoustic sum rule on the Gamma dynamical matrix ("no" | "simple" | + /// "crystal"); applied to the raw matrix before the LO-TO term + void set_asr(const std::string& asr_type); + /// q = 0 response switches (epsilon_inf / Born charges / LO-TO) void set_compute_q0(bool flag); diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index 7448034d703..ed32d6ad795 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -177,6 +177,7 @@ class DFPT_PW::Impl double conv_thr_ = 1e-8; int max_iter_ = 100; double mix_beta_ = 0.4; + std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix }; } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 57a3c06b657..52887dbb6f7 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -112,6 +112,15 @@ void DFPT_PW::Impl::run_assemble(int q_idx) // translational sum-rule audit of the raw Gamma matrix (acoustic-branch // diagnostics; the loto-corrected matrix is audited separately below) phon_.audit_asr("raw", q_idx, data_); + if (asr_type_ != "no") + { + // enforce the sum rule on the raw matrix (with the Z* sum rule), + // re-diagonalize and audit; the LO-TO term below then consumes the + // corrected Born charges (QE dynmat set_asr order) + phon_.apply_asr(asr_type_, q_idx, data_); + phon_.diagonalize(q_idx, data_); + phon_.audit_asr("asr", q_idx, data_); + } if (q_idx != 0 || !data_.get_loto()) { return; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp index 2511020066a..95c8fba08dd 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp @@ -691,3 +691,153 @@ TEST_F(DFPTPhonSerialTest, FormatReportsRegression) " mode 2 : 520.123457 cm^-1\n"; EXPECT_EQ(lrep, expect_l); } + +namespace +{ +/// per-atom row-sum norm of a 2-atom (nat3=6) dynamical matrix +double row_sum_norm(const ModuleBase::ComplexMatrix& dyn, int ia) +{ + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + std::complex s(0.0, 0.0); + for (int jb = 0; jb < 6; ++jb) + { + s += dyn(3 * ia + da, jb); + } + s2 += std::norm(s); + } + return std::sqrt(s2); +} + +/// Hermitian 2-atom matrix with an asymmetric translational violation: +/// self blocks A1 != A2, off block B, plus an extra bump on atom 0 +ModuleBase::ComplexMatrix make_violating_dyn() +{ + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double self1 = 0.05; + const double self2 = 0.04; + const double off = 0.02; + const double bump = 0.012; // asymmetry: only on atom 0, x row + for (int d = 0; d < 3; ++d) + { + dyn(d, d) = std::complex(self1 + (d == 0 ? bump : 0.0), 0.0); + dyn(3 + d, 3 + d) = std::complex(self2, 0.0); + dyn(d, 3 + d) = std::complex(off, 0.01 * d); + dyn(3 + d, d) = std::conj(dyn(d, 3 + d)); + } + return dyn; +} +} // namespace + +TEST_F(DFPTPhonSerialTest, AsrSimpleZeroesRowSums) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + const ModuleBase::ComplexMatrix dyn0 = make_violating_dyn(); + data_.set_dynmat(0, dyn0); + phon_.apply_asr("simple", 0, data_); + const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); + // legacy fix leaves only the tiny imaginary Hermitization residual + const double scale = 0.06; + EXPECT_LT(row_sum_norm(dyn, 0) / scale, 1.0e-12); + EXPECT_LT(row_sum_norm(dyn, 1) / scale, 1.0e-12); + // Hermiticity preserved by the re-symmetrization pass + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(dyn(i, j), std::conj(dyn(j, i))); + } + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalZeroesRowSumsExactlyAndProjectsAcousticToZero) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); + EXPECT_LT(row_sum_norm(dyn, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(dyn, 1), 1.0e-14); + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(dyn(i, j), std::conj(dyn(j, i))); + } + } + // the three translation eigenvalues are projected out exactly: the + // uniform-displacement subspace maps to zero under (1-P) D (1-P). + // The matrix residual is ~1e-17; in frequency units the square root + // rescaling lifts that to ~1e-5 cm^-1, still an exact zero physically + phon_.diagonalize(0, data_); + const std::vector freq = data_.get_phon_freq(0); + std::vector byabs = freq; + std::sort(byabs.begin(), byabs.end(), [](double a, double b) { return std::abs(a) < std::abs(b); }); + for (int m = 0; m < 3; ++m) + { + EXPECT_NEAR(byabs[m], 0.0, 1.0e-4); + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalKeepsEquivalentAtomsEquivalent) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + // symmetric input: identical self blocks, common off block, plus a + // violation spread equally (row sums nonzero but atom-symmetric) + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double self = 0.05; + const double off = 0.02; + for (int d = 0; d < 3; ++d) + { + dyn(d, d) = std::complex(self, 0.0); + dyn(3 + d, 3 + d) = std::complex(self, 0.0); + dyn(d, 3 + d) = std::complex(off, 0.0); + dyn(3 + d, d) = std::complex(off, 0.0); + } + data_.set_dynmat(0, dyn); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + // crystal redistributes the mean uniformly: both self blocks stay equal + // and both row sums vanish (simple would dump them on the diagonals, + // also staying equal here -- the discriminating case is the asymmetric + // one above, where only crystal shares the correction) + for (int d = 0; d < 3; ++d) + { + EXPECT_EQ(out(d, d), out(3 + d, 3 + d)); + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalCorrectsBornChargeSumRule) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + ModuleBase::matrix z1(3, 3, true); + ModuleBase::matrix z2(3, 3, true); + z1(0, 0) = 1.3; + z1(1, 1) = 1.1; + z1(2, 2) = 1.2; + z2(0, 0) = -1.1; + z2(1, 1) = -1.0; + z2(2, 2) = -1.05; // sums: 0.2 / 0.1 / 0.15 != 0 + data_.set_born(0, z1); + data_.set_born(1, z2); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::matrix o1 = data_.get_born(0); + const ModuleBase::matrix o2 = data_.get_born(1); + for (int a = 0; a < 3; ++a) + { + EXPECT_NEAR(o1(a, a) + o2(a, a), 0.0, 1.0e-14); + } + // shifted by the per-component mean: 0.1 / 0.05 / 0.075 + EXPECT_NEAR(o1(0, 0), 1.3 - 0.1, 1.0e-14); + EXPECT_NEAR(o2(1, 1), -1.1 + 0.05, 1.0e-14); +} From f3e260fb150cf5f1877aa89c0631b7d55f540e57 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 19 Oct 2026 14:05:17 +0800 Subject: [PATCH 3/6] fix(dfpt): mass-weight the acoustic sum-rule constraints The first crystal/simple implementation enforced UNWEIGHTED row sums on the stored mass-REDUCED matrix (D = Phi/sqrt(m_a m_b)), which projects out the wrong direction (uniform mass-coordinate displacement instead of uniform physical displacement). Consequence caught by the first ten-system validation run: the pre-LO-TO acoustic modes zero out along the spurious direction, but the LO-TO nonanalytic term -- which decouples from the PHYSICAL acoustic subspace only through sum_a Z*_a = 0 -- mixes back in and lifts one acoustic mode (bn_cubic Gamma: +36 cm^-1 along (1,1,1) while QE stays at ~1e-5). QE dynmat operates on the UNreduced matrix (dyndiag divides by the sqrt(m_a m_b) only at diagonalization), so its unweighted rows are the physical rule; the correct translation into our reduced storage carries the sqrt(m) weights: - crystal: project with the orthonormal physical translations w_i(3a+alpha) = sqrt(m_a) delta(alpha,i)/sqrt(sum m), i.e. D' = D - (sqrt(m_a)/M) S - (sqrt(m_b)/M) T + (sqrt(m_a m_b)/M^2) U with S/T/U the sqrt(m)-weighted single and double atom sums - simple: D(a,a) = -sum_{b!=a} sqrt(m_b/m_a) D(a,b) - the audit and the report line now print the PHYSICAL row sums sum_jb sqrt(m_jb) D(3a+da, jb), the quantity QE's rows sum to Tests (phon serial 16 -> 17): row-sum assertions move to the physical sums; the equivalent-atoms case runs with mass-equal species; new regression AsrCrystalPlusLotoKeepsAcousticExactlyZero (masses 12/4, violating matrix, eps 3I and unbalanced Z*): after crystal + LO-TO the acoustic triplet stays at zero to 1e-4 cm^-1. --- source/source_pw/module_dfpt/dfpt_phon.cpp | 117 ++++++++++++------ .../test_serial/dfpt_phon_serial_test.cpp | 57 +++++++-- 2 files changed, 127 insertions(+), 47 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index a9861d98a5b..f87561db7ca 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -187,7 +187,15 @@ void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Dat max_elem = std::max(max_elem, std::abs(dyn(i, j))); } } - // per-atom translational row sums s(ia)[da] = sum_ib,db D_{ia da, ib db} + // per-atom PHYSICAL translational row sums: the stored matrix is + // mass-reduced (D = Phi/sqrt(m_a m_b)), so the Phi row sums read + // s(ia)[da] = sum_jb sqrt(m_jb) D(3ia+da, jb) -- the same quantity + // QE's unreduced dyn rows sum to + std::vector sqm(nat, 0.0); + for (int ia = 0; ia < nat; ++ia) + { + sqm[ia] = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass); + } std::cout << "DFPT ASR audit (" << label << "):"; double worst = 0.0; int worst_atom = -1; @@ -199,7 +207,7 @@ void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Dat std::complex s(0.0, 0.0); for (int jb = 0; jb < nat3; ++jb) { - s += dyn(3 * ia + da, jb); + s += sqm[jb / 3] * dyn(3 * ia + da, jb); } s2 += std::norm(s); } @@ -240,49 +248,81 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& return; } - // violation BEFORE the correction (for the report line) + // The stored matrix is mass-reduced, D = Phi / sqrt(m_a m_b) (amu masses + // from the cell); the PHYSICAL sum rule acts on Phi rows, + // sum_b sqrt(m_b) D(3a+da, 3b+db) = 0, and the acoustic eigenvectors of + // D are v(3a+da) = sqrt(m_a) e_da (uniform physical displacement). QE + // dynmat works on the unreduced Phi, so its unweighted rows are this + // same rule; here the constraint vectors carry the sqrt(m) weights. + std::vector sqm(nat, 0.0); + double mtot = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + sqm[ia] = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass); + mtot += ucell_->atoms[ucell_->iat2it[ia]].mass; + } + const auto physical_row_sum = [&](int ia, int da) { + std::complex s(0.0, 0.0); + for (int jb = 0; jb < nat3; ++jb) + { + s += sqm[jb / 3] * dyn(3 * ia + da, jb); + } + return s; + }; + + // violation BEFORE the correction (physical row sums) double worst_before = 0.0; for (int ia = 0; ia < nat; ++ia) { double s2 = 0.0; for (int da = 0; da < 3; ++da) { - std::complex s(0.0, 0.0); - for (int jb = 0; jb < nat3; ++jb) - { - s += dyn(3 * ia + da, jb); - } - s2 += std::norm(s); + s2 += std::norm(physical_row_sum(ia, da)); } worst_before = std::max(worst_before, std::sqrt(s2)); } if (asr_type == "simple") { - // legacy fix: each row's total is subtracted from its diagonal - // element (QE dynmat 'simple'); Hermiticity is restored by the - // re-symmetrization pass below, mirroring QE's upper-triangle use - for (int i = 0; i < nat3; ++i) + // legacy fix, QE 'simple' translated to the reduced matrix: + // D(a,a) = -sum_{b!=a} sqrt(m_b/m_a) D(a,b), zeroes the physical + // row sums; Hermiticity is restored by the re-symmetrization below + for (int na = 0; na < nat; ++na) { - std::complex s(0.0, 0.0); - for (int j = 0; j < nat3; ++j) + for (int da = 0; da < 3; ++da) { - s += dyn(i, j); + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int nb = 0; nb < nat; ++nb) + { + if (nb != na) + { + s += sqm[nb] * dyn(3 * na + da, 3 * nb + db); + } + } + dyn(3 * na + da, 3 * na + db) = -s / sqm[na]; + } } - dyn(i, i) -= s; } } else if (asr_type == "crystal") { // minimal-norm symmetric projection D -> (1-P) D (1-P) with the - // orthonormal translations u_i(3a+alpha) = delta(alpha,i)/sqrt(nat): - // M(alpha, nb, beta) = (1/nat) sum_na' D(na',alpha; nb,beta) (P D) - // N(na, alpha, beta) = (1/nat) sum_nb' D(na,alpha; nb',beta) (D P) - // T(alpha, beta) = sum_na,nb D(na,alpha; nb,beta) - // D' = D - M - N + T/nat^2 - std::vector> m(3 * nat3, std::complex(0.0, 0.0)); - std::vector> n(3 * nat3, std::complex(0.0, 0.0)); - ModuleBase::ComplexMatrix t(3, 3, true); + // orthonormal physical translations w_i(3a+da) = sqrt(m_a) + // delta(da,i)/sqrt(sum m): zeroes every physical row AND column sum, + // keeps Hermiticity and symmetry-equivalent atoms equivalent, and + // leaves the acoustic eigenvectors in the null space so the LO-TO + // term (which decouples from them through sum_a Z*_a = 0) cannot + // lift the acoustic frequencies + // S(da, nb, db) = sum_na' sqrt(m_a') D(na',da; nb,db) + // T(na, da, db) = sum_nb' sqrt(m_b') D(na,da; nb',db) + // U(da, db) = sum_na,nb sqrt(m_a) sqrt(m_b) D(na,da; nb,db) + // D' = D - (sqrt(m_na)/M) S - (sqrt(m_nb)/M) T + (sqrt(m_na m_nb)/M^2) U + // with M = sum_a m_a + std::vector> s(3 * nat3, std::complex(0.0, 0.0)); + std::vector> t(3 * nat3, std::complex(0.0, 0.0)); + ModuleBase::ComplexMatrix u(3, 3, true); for (int na = 0; na < nat; ++na) { for (int da = 0; da < 3; ++da) @@ -292,9 +332,9 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& for (int db = 0; db < 3; ++db) { const std::complex d = dyn(3 * na + da, 3 * nb + db); - m[da * nat3 + 3 * nb + db] += d / static_cast(nat); - n[(3 * na + da) * 3 + db] += d / static_cast(nat); - t(da, db) += d; + s[da * nat3 + 3 * nb + db] += sqm[na] * d; + t[(3 * na + da) * 3 + db] += sqm[nb] * d; + u(da, db) += sqm[na] * sqm[nb] * d; } } } @@ -307,9 +347,11 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& { for (int db = 0; db < 3; ++db) { - dyn(3 * na + da, 3 * nb + db) -= m[da * nat3 + 3 * nb + db]; - dyn(3 * na + da, 3 * nb + db) -= n[(3 * na + da) * 3 + db]; - dyn(3 * na + da, 3 * nb + db) += t(da, db) / static_cast(nat * nat); + std::complex d = dyn(3 * na + da, 3 * nb + db); + d -= sqm[na] / mtot * s[da * nat3 + 3 * nb + db]; + d -= sqm[nb] / mtot * t[(3 * na + da) * 3 + db]; + d += sqm[na] * sqm[nb] / (mtot * mtot) * u(da, db); + dyn(3 * na + da, 3 * nb + db) = d; } } } @@ -379,24 +421,19 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data.set_dynmat(q_idx, dyn); - // violation AFTER the correction + // violation AFTER the correction (physical row sums) double worst_after = 0.0; for (int ia = 0; ia < nat; ++ia) { double s2 = 0.0; for (int da = 0; da < 3; ++da) { - std::complex s(0.0, 0.0); - for (int jb = 0; jb < nat3; ++jb) - { - s += dyn(3 * ia + da, jb); - } - s2 += std::norm(s); + s2 += std::norm(physical_row_sum(ia, da)); } worst_after = std::max(worst_after, std::sqrt(s2)); } - std::cout << "DFPT ASR (" << asr_type << ") applied: worst rowsum " << std::setprecision(3) << std::scientific - << worst_before << " -> " << worst_after; + std::cout << "DFPT ASR (" << asr_type << ") applied: worst physical rowsum " << std::setprecision(3) + << std::scientific << worst_before << " -> " << worst_after; if (n_born == nat) { std::cout << " | max |sum_a Z*_a| " << zworst_before << " -> 0"; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp index 95c8fba08dd..13abde6d1ee 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp @@ -694,16 +694,19 @@ TEST_F(DFPTPhonSerialTest, FormatReportsRegression) namespace { -/// per-atom row-sum norm of a 2-atom (nat3=6) dynamical matrix -double row_sum_norm(const ModuleBase::ComplexMatrix& dyn, int ia) +/// per-atom PHYSICAL row-sum norm of a 2-atom (nat3=6) mass-reduced +/// dynamical matrix: sum_jb sqrt(m_jb) D(3ia+da, jb), masses from the cell +double row_sum_norm(const UnitCell& ucell, const ModuleBase::ComplexMatrix& dyn, int ia) { + double sqm[2] = {std::sqrt(ucell.atoms[ucell.iat2it[0]].mass), + std::sqrt(ucell.atoms[ucell.iat2it[1]].mass)}; double s2 = 0.0; for (int da = 0; da < 3; ++da) { std::complex s(0.0, 0.0); for (int jb = 0; jb < 6; ++jb) { - s += dyn(3 * ia + da, jb); + s += sqm[jb / 3] * dyn(3 * ia + da, jb); } s2 += std::norm(s); } @@ -741,8 +744,8 @@ TEST_F(DFPTPhonSerialTest, AsrSimpleZeroesRowSums) const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); // legacy fix leaves only the tiny imaginary Hermitization residual const double scale = 0.06; - EXPECT_LT(row_sum_norm(dyn, 0) / scale, 1.0e-12); - EXPECT_LT(row_sum_norm(dyn, 1) / scale, 1.0e-12); + EXPECT_LT(row_sum_norm(ucell_, dyn, 0) / scale, 1.0e-12); + EXPECT_LT(row_sum_norm(ucell_, dyn, 1) / scale, 1.0e-12); // Hermiticity preserved by the re-symmetrization pass for (int i = 0; i < 6; ++i) { @@ -761,8 +764,8 @@ TEST_F(DFPTPhonSerialTest, AsrCrystalZeroesRowSumsExactlyAndProjectsAcousticToZe data_.set_dynmat(0, make_violating_dyn()); phon_.apply_asr("crystal", 0, data_); const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); - EXPECT_LT(row_sum_norm(dyn, 0), 1.0e-14); - EXPECT_LT(row_sum_norm(dyn, 1), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, dyn, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, dyn, 1), 1.0e-14); for (int i = 0; i < 6; ++i) { for (int j = i + 1; j < 6; ++j) @@ -787,6 +790,9 @@ TEST_F(DFPTPhonSerialTest, AsrCrystalZeroesRowSumsExactlyAndProjectsAcousticToZe TEST_F(DFPTPhonSerialTest, AsrCrystalKeepsEquivalentAtomsEquivalent) { MakeTwoAtomCell(); + // symmetry-equivalent atoms carry equal masses; make the two fixture + // species mass-equal so the equivalence assertion is meaningful + ucell_.atoms[1].mass = ucell_.atoms[0].mass; // rewire the q list to Gamma (the fixture default q is not 0) SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); // symmetric input: identical self blocks, common off block, plus a @@ -841,3 +847,40 @@ TEST_F(DFPTPhonSerialTest, AsrCrystalCorrectsBornChargeSumRule) EXPECT_NEAR(o1(0, 0), 1.3 - 0.1, 1.0e-14); EXPECT_NEAR(o2(1, 1), -1.1 + 0.05, 1.0e-14); } + +TEST_F(DFPTPhonSerialTest, AsrCrystalPlusLotoKeepsAcousticExactlyZero) +{ + // regression of the mass-weighting bug: with UNequal masses the + // translation constraints must live on the physical (sqrt(m)-weighted) + // vectors, otherwise the LO-TO nonanalytic term -- which decouples from + // the physical acoustic subspace only through sum_a Z*_a = 0 -- mixes + // into the projected null space and lifts one acoustic frequency (the + // bn_cubic Gamma run showed +36 cm^-1 before the fix) + MakeTwoAtomCell(); // masses 12 / 4: unequal, the trigger + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + ModuleBase::matrix eps(3, 3, true); + for (int d = 0; d < 3; ++d) + { + eps(d, d) = 3.0; + } + data_.set_dielectric(eps); + ModuleBase::matrix z1(3, 3, true); + ModuleBase::matrix z2(3, 3, true); + z1(0, 0) = z1(1, 1) = z1(2, 2) = 1.3; + z2(0, 0) = z2(1, 1) = z2(2, 2) = -1.1; // sums 0.2 / 0.2 / 0.2 != 0 + data_.set_born(0, z1); + data_.set_born(1, z2); + phon_.apply_asr("crystal", 0, data_); + // LO-TO along (1,1,1)/sqrt(3): the acoustic triplet must stay at zero + phon_.add_loto(ModuleBase::Vector3(1.0, 1.0, 1.0), data_); + phon_.diagonalize_loto(data_); + const std::vector freq = data_.get_phon_freq_loto(); + ASSERT_EQ(freq.size(), 6u); + std::vector byabs = freq; + std::sort(byabs.begin(), byabs.end(), [](double a, double b) { return std::abs(a) < std::abs(b); }); + for (int m = 0; m < 3; ++m) + { + EXPECT_NEAR(byabs[m], 0.0, 1.0e-4); + } +} From f7c9bd39e2cd0166fd55a7664321fc0f6154d775 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 19 Oct 2026 18:09:00 +0800 Subject: [PATCH 4/6] fix(dfpt): do not dereference iat2it in the ASR audit; complete the stru fixture The CI MODULE_DFPT_pw_run_test (GCC 11 container) segfaulted inside the first run() of the ASR branch while the same binary passed everywhere locally (GCC 13): DFPTStruTestFixture::construct_ucell never populated ucell.iat2it, so the audit's cosmetic species-label lookup atoms[iat2it[ia]].label read an EMPTY vector -- undefined behavior that GCC 13 layouts happened to survive and GCC 11 did not. - audit_asr: print the atom index only; production code no longer touches iat2it (apply_asr's mass lookup still does, legitimately) - the stru fixture now fills iat2it/iat2ia like any real cell --- source/source_pw/module_dfpt/dfpt_phon.cpp | 3 +-- .../module_dfpt/test/dfpt_stru_fixture.cpp | 15 +++++++++++++++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index f87561db7ca..cc84298e080 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -212,8 +212,7 @@ void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Dat s2 += std::norm(s); } const double snorm = std::sqrt(s2); - std::cout << " atom" << ia + 1 << "(" << ucell_->atoms[ucell_->iat2it[ia]].label << ")=" - << std::setprecision(3) << std::scientific << snorm; + std::cout << " atom" << ia + 1 << "=" << std::setprecision(3) << std::scientific << snorm; if (snorm > worst) { worst = snorm; diff --git a/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp b/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp index 4df392da1cf..4192fe161b0 100644 --- a/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp @@ -64,6 +64,21 @@ void DFPTStruTestFixture::construct_ucell(stru_& stru) } ucell.nat += ucell.atoms[i].na; } + // atom-index tables: real cells always carry them (iat2it/iat2ia map the + // running atom index to type and in-type index); populate them so any + // consumer (e.g. DFPT_Phon::apply_asr mass lookup) works on this fixture + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); + int iat = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + ucell.iat2it[iat] = it; + ucell.iat2ia[iat] = ia; + ++iat; + } + } } void DFPTStruTestFixture::ClearUcell() From 64fcadf9585e4b2c552773254a67f1406ac520c8 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Tue, 20 Oct 2026 14:03:16 +0800 Subject: [PATCH 5/6] fix(dfpt): resolve the two P2 review findings on the ASR implementation [P2: directional cancellation in the audit metric] The audit and the apply_asr report summed over ALL columns (atoms AND Cartesian directions) before taking the norm, so opposite x/y residuals could cancel: the block [[1,-1,0],[-1,1,0],[0,0,0]] reported zero although D*t_x != 0. Both the audit and the report now compute per-direction residuals R(a,alpha,beta) = sum_b sqrt(m_b) D(3a+alpha, 3b+beta) at a FIXED translation direction beta, then combine the norms. New regression: the cancellation pattern must be detected. [P2: Hermitian re-symmetrization breaks the simple correction] The simple mode dumped the residual onto the onsite 3x3 blocks, which can write complex values on the DIAGONAL (the j>i symmetrization loop leaves the diagonal untouched), producing a non-Hermitian input to the crystal projection -- whose result is then non-Hermitian, and the shared re-symmetrization halves the correction back (the review's +/-delta/2 example). Fix: (1) take the real part of the diagonal after the onsite dump, making the input genuinely Hermitian before the projection; (2) finish with the crystal projection, which now gives an exact Hermitian matrix satisfying all three translation constraints. New regression: cross-Cartesian coupling (atom0-x <-> atom1-y) with unequal masses, checking translations exact + Hermiticity + real diagonal. Also addresses the two non-blocking notes: - header doc: correct translation vectors to sqrt(m_a) delta/sqrt(M), describe simple as onsite-block correction + projection finish - run_assemble: single eigensolve after any applicable correction (the raw audit does not consume frequencies) --- source/source_pw/module_dfpt/dfpt_phon.cpp | 155 +++++++++++------- source/source_pw/module_dfpt/dfpt_phon.h | 31 +++- source/source_pw/module_dfpt/dfpt_pw_run.cpp | 16 +- .../test_serial/dfpt_phon_serial_test.cpp | 92 ++++++++++- 4 files changed, 212 insertions(+), 82 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index cc84298e080..de9fa67306c 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -201,15 +201,21 @@ void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Dat int worst_atom = -1; for (int ia = 0; ia < nat; ++ia) { + // per-DIRECTION residuals R(a,alpha,i) = sum_b sqrt(m_b) D(3a+alpha, + // 3b+i): summing over atoms at a FIXED translation direction i + // before combining norms, so opposite x/y residuals cannot cancel double s2 = 0.0; for (int da = 0; da < 3; ++da) { - std::complex s(0.0, 0.0); - for (int jb = 0; jb < nat3; ++jb) + for (int db = 0; db < 3; ++db) { - s += sqm[jb / 3] * dyn(3 * ia + da, jb); + std::complex s(0.0, 0.0); + for (int ib = 0; ib < nat; ++ib) + { + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); + } + s2 += std::norm(s); } - s2 += std::norm(s); } const double snorm = std::sqrt(s2); std::cout << " atom" << ia + 1 << "=" << std::setprecision(3) << std::scientific << snorm; @@ -260,65 +266,45 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& sqm[ia] = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass); mtot += ucell_->atoms[ucell_->iat2it[ia]].mass; } - const auto physical_row_sum = [&](int ia, int da) { - std::complex s(0.0, 0.0); - for (int jb = 0; jb < nat3; ++jb) - { - s += sqm[jb / 3] * dyn(3 * ia + da, jb); - } - return s; - }; - - // violation BEFORE the correction (physical row sums) - double worst_before = 0.0; - for (int ia = 0; ia < nat; ++ia) - { + // per-DIRECTION physical residual R(a,alpha,i) = sum_b sqrt(m_b) + // D(3a+alpha, 3b+i): summing over atoms at a fixed translation + // direction i so opposite x/y residuals cannot cancel + const auto physical_residual_norm = [&](int ia) { double s2 = 0.0; for (int da = 0; da < 3; ++da) { - s2 += std::norm(physical_row_sum(ia, da)); - } - worst_before = std::max(worst_before, std::sqrt(s2)); - } - - if (asr_type == "simple") - { - // legacy fix, QE 'simple' translated to the reduced matrix: - // D(a,a) = -sum_{b!=a} sqrt(m_b/m_a) D(a,b), zeroes the physical - // row sums; Hermiticity is restored by the re-symmetrization below - for (int na = 0; na < nat; ++na) - { - for (int da = 0; da < 3; ++da) + for (int db = 0; db < 3; ++db) { - for (int db = 0; db < 3; ++db) + std::complex s(0.0, 0.0); + for (int ib = 0; ib < nat; ++ib) { - std::complex s(0.0, 0.0); - for (int nb = 0; nb < nat; ++nb) - { - if (nb != na) - { - s += sqm[nb] * dyn(3 * na + da, 3 * nb + db); - } - } - dyn(3 * na + da, 3 * na + db) = -s / sqm[na]; + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); } + s2 += std::norm(s); } } - } - else if (asr_type == "crystal") + return std::sqrt(s2); + }; + + // violation BEFORE the correction (physical residuals) + double worst_before = 0.0; + for (int ia = 0; ia < nat; ++ia) { - // minimal-norm symmetric projection D -> (1-P) D (1-P) with the - // orthonormal physical translations w_i(3a+da) = sqrt(m_a) - // delta(da,i)/sqrt(sum m): zeroes every physical row AND column sum, - // keeps Hermiticity and symmetry-equivalent atoms equivalent, and - // leaves the acoustic eigenvectors in the null space so the LO-TO - // term (which decouples from them through sum_a Z*_a = 0) cannot - // lift the acoustic frequencies - // S(da, nb, db) = sum_na' sqrt(m_a') D(na',da; nb,db) - // T(na, da, db) = sum_nb' sqrt(m_b') D(na,da; nb',db) - // U(da, db) = sum_na,nb sqrt(m_a) sqrt(m_b) D(na,da; nb,db) - // D' = D - (sqrt(m_na)/M) S - (sqrt(m_nb)/M) T + (sqrt(m_na m_nb)/M^2) U - // with M = sum_a m_a + worst_before = std::max(worst_before, physical_residual_norm(ia)); + } + + // the crystal projection D -> (1-P) D (1-P) with the orthonormal + // physical translations w_i(3a+da) = sqrt(m_a) delta(da,i)/sqrt(sum m): + // zeroes every per-direction physical residual exactly, keeps + // Hermiticity and symmetry-equivalent atoms equivalent, and leaves the + // acoustic eigenvectors in the null space so the LO-TO term (which + // decouples from them through sum_a Z*_a = 0) cannot lift the acoustic + // frequencies. Both 'simple' and 'crystal' finish with this projection: + // 'simple' first dumps the residual onto the onsite blocks (the QE + // legacy flavor), then the projection restores Hermiticity and removes + // whatever onsite-only correction could not absorb (cross-Cartesian + // antisymmetric residuals have no exact Hermitian onsite-only fix) + const auto apply_crystal_projection = [&]() { std::vector> s(3 * nat3, std::complex(0.0, 0.0)); std::vector> t(3 * nat3, std::complex(0.0, 0.0)); ModuleBase::ComplexMatrix u(3, 3, true); @@ -355,6 +341,56 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& } } } + }; + + if (asr_type == "simple") + { + // legacy QE 'simple' translated to the reduced matrix: dump the + // per-direction residual onto the onsite 3x3 block, + // D(a,a)[da][db] = -sum_{b!=a} sqrt(m_b/m_a) D(a,da; b,db); + // the onsite dump may break Hermiticity (cross-Cartesian parts), + // so symmetrize, then the crystal projection below removes the + // residual the symmetrization reintroduced (at most half the + // anti-Hermitian part) and guarantees the final matrix is + // Hermitian and exactly satisfies all three translation constraints + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int nb = 0; nb < nat; ++nb) + { + if (nb != na) + { + s += sqm[nb] * dyn(3 * na + da, 3 * nb + db); + } + } + dyn(3 * na + da, 3 * na + db) = -s / sqm[na]; + } + } + } + // Hermitian re-symmetrization of the onsite dump (the projection + // is Hermitian-preserving only on Hermitian input); the j > i loop + // alone leaves the diagonal untouched, and the onsite dump may + // have written complex values there -- take the real part so the + // matrix is genuinely Hermitian before the projection + for (int i = 0; i < nat3; ++i) + { + for (int j = i + 1; j < nat3; ++j) + { + const std::complex avg = 0.5 * (dyn(i, j) + std::conj(dyn(j, i))); + dyn(i, j) = avg; + dyn(j, i) = std::conj(avg); + } + dyn(i, i) = std::complex(dyn(i, i).real(), 0.0); + } + apply_crystal_projection(); + } + else if (asr_type == "crystal") + { + apply_crystal_projection(); } else { @@ -420,18 +456,13 @@ void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data.set_dynmat(q_idx, dyn); - // violation AFTER the correction (physical row sums) + // violation AFTER the correction (physical residuals) double worst_after = 0.0; for (int ia = 0; ia < nat; ++ia) { - double s2 = 0.0; - for (int da = 0; da < 3; ++da) - { - s2 += std::norm(physical_row_sum(ia, da)); - } - worst_after = std::max(worst_after, std::sqrt(s2)); + worst_after = std::max(worst_after, physical_residual_norm(ia)); } - std::cout << "DFPT ASR (" << asr_type << ") applied: worst physical rowsum " << std::setprecision(3) + std::cout << "DFPT ASR (" << asr_type << ") applied: worst physical residual " << std::setprecision(3) << std::scientific << worst_before << " -> " << worst_after; if (n_born == nat) { diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 0452805d41e..9bcd8dbaec8 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -97,15 +97,28 @@ class DFPT_Phon /// enforce the translational acoustic sum rule on the Gamma dynamical /// matrix (and, in every enforcing mode, the Z* sum rule /// sum_a Z*_a = 0 on the stored Born charges, like QE dynmat's - /// set_asr). "simple": the legacy fix, each row's total sum is - /// subtracted from its diagonal element (QE 'simple'). "crystal": - /// the optimized minimal-norm fix -- with the three orthonormal - /// translation vectors u_i(3a+alpha) = delta(alpha,i)/sqrt(nat), the - /// symmetric projection D -> (1-P) D (1-P), P = sum_i u_i u_i^T, - /// which zeroes every row AND column sum, preserves Hermiticity and - /// keeps symmetry-equivalent atoms equivalent (QE 'crystal' for the - /// translational part). Both apply to the RAW matrix, before the - /// LO-TO nonanalytic term; away from Gamma the call is a no-op. + /// set_asr). Both modes operate on the mass-reduced matrix D = + /// Phi/sqrt(m_a m_b); the physical rule reads sum_b sqrt(m_b) + /// D(3a+alpha, 3b+beta) = 0 for each translation direction beta + /// independently. + /// + /// Both modes finish with the exact Hermitian projection D -> + /// (1-P) D (1-P), P = sum_i w_i w_i^T, with the orthonormal physical + /// translations w_i(3a+alpha) = sqrt(m_a) delta(alpha,i) / + /// sqrt(sum m): the final matrix satisfies all three mass-weighted + /// translation constraints exactly and is Hermitian. + /// + /// "simple": the QE legacy flavor -- the residual is first dumped + /// onto the onsite 3x3 blocks, D(a,a)[alpha][beta] = + /// -sum_{b!=a} sqrt(m_b/m_a) D(a,alpha; b,beta), then the projection + /// absorbs whatever onsite-only correction could not (a purely + /// onsite Hermitian fix does not exist when the residual has a + /// cross-Cartesian antisymmetric part). + /// + /// "crystal": the optimized minimal-norm fix -- the projection alone, + /// which also keeps symmetry-equivalent atoms equivalent. + /// + /// Away from Gamma the call is a no-op. void apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data); /// Ewald ion-ion force constants C^ewald_ab(q) (G-space + real-space + diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 52887dbb6f7..9db6f17e942 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -108,19 +108,25 @@ void DFPT_PW::Impl::run_q0_post(int q_idx) void DFPT_PW::Impl::run_assemble(int q_idx) { phon_.assemble(q_idx, data_); - phon_.diagonalize(q_idx, data_); // translational sum-rule audit of the raw Gamma matrix (acoustic-branch - // diagnostics; the loto-corrected matrix is audited separately below) - phon_.audit_asr("raw", q_idx, data_); + // diagnostics; the loto-corrected matrix is audited separately below). + // The raw audit does not consume frequencies, so the ordinary + // eigensolve runs once after any applicable ASR correction if (asr_type_ != "no") { + phon_.audit_asr("raw", q_idx, data_); // enforce the sum rule on the raw matrix (with the Z* sum rule), - // re-diagonalize and audit; the LO-TO term below then consumes the - // corrected Born charges (QE dynmat set_asr order) + // diagonalize once and audit; the LO-TO term below then consumes + // the corrected Born charges (QE dynmat set_asr order) phon_.apply_asr(asr_type_, q_idx, data_); phon_.diagonalize(q_idx, data_); phon_.audit_asr("asr", q_idx, data_); } + else + { + phon_.diagonalize(q_idx, data_); + phon_.audit_asr("raw", q_idx, data_); + } if (q_idx != 0 || !data_.get_loto()) { return; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp index 13abde6d1ee..085de1a93da 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp @@ -694,8 +694,10 @@ TEST_F(DFPTPhonSerialTest, FormatReportsRegression) namespace { -/// per-atom PHYSICAL row-sum norm of a 2-atom (nat3=6) mass-reduced -/// dynamical matrix: sum_jb sqrt(m_jb) D(3ia+da, jb), masses from the cell +/// per-atom PHYSICAL residual norm of a 2-atom (nat3=6) mass-reduced +/// dynamical matrix: R(a,alpha,beta) = sum_b sqrt(m_b) D(3a+alpha, 3b+beta) +/// at a FIXED direction beta, then combined -- opposite x/y residuals +/// cannot cancel double row_sum_norm(const UnitCell& ucell, const ModuleBase::ComplexMatrix& dyn, int ia) { double sqm[2] = {std::sqrt(ucell.atoms[ucell.iat2it[0]].mass), @@ -703,12 +705,15 @@ double row_sum_norm(const UnitCell& ucell, const ModuleBase::ComplexMatrix& dyn, double s2 = 0.0; for (int da = 0; da < 3; ++da) { - std::complex s(0.0, 0.0); - for (int jb = 0; jb < 6; ++jb) + for (int db = 0; db < 3; ++db) { - s += sqm[jb / 3] * dyn(3 * ia + da, jb); + std::complex s(0.0, 0.0); + for (int ib = 0; ib < 2; ++ib) + { + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); + } + s2 += std::norm(s); } - s2 += std::norm(s); } return std::sqrt(s2); } @@ -884,3 +889,78 @@ TEST_F(DFPTPhonSerialTest, AsrCrystalPlusLotoKeepsAcousticExactlyZero) EXPECT_NEAR(byabs[m], 0.0, 1.0e-4); } } + +TEST_F(DFPTPhonSerialTest, AsrDirectionalCancellationDetected) +{ + // regression of the audit-metric fix: the block [[1,-1,0],[-1,1,0]] + // has zero row sums (sum over ALL columns) but a NONZERO per-direction + // residual D*t_x; the audit and apply_asr report must not miss it + MakeTwoAtomCell(); + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + ModuleBase::ComplexMatrix dyn(6, 6, true); + // atom0 self-block with the cancellation pattern + dyn(0, 0) = std::complex(1.0, 0.0); + dyn(0, 1) = std::complex(-1.0, 0.0); + dyn(1, 0) = std::complex(-1.0, 0.0); + dyn(1, 1) = std::complex(1.0, 0.0); + // atom1 self-block mirrors it (keeps total charge-conservation flavor) + dyn(3, 3) = std::complex(1.0, 0.0); + dyn(3, 4) = std::complex(-1.0, 0.0); + dyn(4, 3) = std::complex(-1.0, 0.0); + dyn(4, 4) = std::complex(1.0, 0.0); + data_.set_dynmat(0, dyn); + + // per-direction residual BEFORE any correction is nonzero + const double residual_before = row_sum_norm(ucell_, dyn, 0); + EXPECT_GT(residual_before, 0.1); + + // crystal ASR zeroes it + ModuleBase::ComplexMatrix dyn_copy = dyn; + data_.set_dynmat(0, dyn_copy); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + EXPECT_LT(row_sum_norm(ucell_, out, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, out, 1), 1.0e-14); +} + +TEST_F(DFPTPhonSerialTest, AsrSimpleCrossCartesianRegression) +{ + // regression of the simple-mode Hermitian fix: the cross-Cartesian + // coupling D(0,4) = D(4,0) = delta (atom0-x <-> atom1-y) has no exact + // Hermitian onsite-only fix; simple must still produce a final matrix + // that is Hermitian, satisfies all three translation constraints + // exactly, and has a real diagonal + MakeTwoAtomCell(); // masses 12/4 (unequal) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double delta = 0.07; + dyn(0, 4) = std::complex(delta, 0.0); + dyn(4, 0) = std::complex(delta, 0.0); + // also a diagonal self term so the matrix is not purely offsite + dyn(0, 0) = std::complex(0.05, 0.0); + dyn(4, 4) = std::complex(0.04, 0.0); + dyn(2, 2) = std::complex(0.03, 0.0); + dyn(5, 5) = std::complex(0.06, 0.0); + data_.set_dynmat(0, dyn); + phon_.apply_asr("simple", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + + // translation constraints hold exactly (per-direction residuals) + EXPECT_LT(row_sum_norm(ucell_, out, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, out, 1), 1.0e-14); + + // Hermiticity + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(out(i, j), std::conj(out(j, i))) << "i=" << i << " j=" << j; + } + } + + // diagonal is real + for (int i = 0; i < 6; ++i) + { + EXPECT_EQ(out(i, i).imag(), 0.0) << "i=" << i; + } +} From c3319e6b737beeaeba5b642ce62c0df2c0d98340 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Tue, 20 Oct 2026 18:33:01 +0800 Subject: [PATCH 6/6] feat(dfpt): nlcc core-density first-order response in the DFPT phonon loop When a pseudopotential with nonlinear core correction (nlcc) displaces, its core charge density moves with the atom. The first-order density drho must include this contribution in the XC channel (the core Coulomb interaction is already absorbed in the ionic pseudopotential, so no Hartree term). Previously compute_drho carried only the valence u*du cross terms, leaving nlcc systems' response incomplete: bn_hex Gamma soft-mode pair was 8.7x too imaginary versus QE (-259.7 vs -29.8 cm^-1) while non-nlcc systems were unaffected. New DFPT_Rho::compute_drhoc computes the q-shifted G-space amplitude d rho_core / d tau_dir (G+q) = -i * tpiba * (G+q)_dir * rho_core~(|G+q|) * exp(i (G+q) tau_iat) where rho_core~ is the radial spherical-Bessel Fourier transform of the core charge (evaluated per-|G+q| from the pseudopotential's rho_atc radial data, the same integral the ground state's non_linear_core_correction performs). Non-nlcc atoms give zero. assemble_v_sc injects the core derivative on the XC input side only: a_r += drhoc_r before xc_->apply; the Hartree above keeps the pure valence density. solve_displacement computes drhoc once per displacement (it is iteration-independent) alongside build_dv. Verified: all existing DFPT serial suites pass (rho/pert/phon/q0/ xc_analytic 8/7/19/6/6); nlcc validation runs in flight. --- source/source_pw/module_dfpt/dfpt_pw_data.cpp | 30 ++++++ source/source_pw/module_dfpt/dfpt_pw_data.h | 7 ++ .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 16 +++ source/source_pw/module_dfpt/dfpt_rho.cpp | 100 ++++++++++++++++++ source/source_pw/module_dfpt/dfpt_rho.h | 17 +++ 5 files changed, 170 insertions(+) diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp index 754f4e063df..529517c6212 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp @@ -646,4 +646,34 @@ void DFPT_PW_Data::deallocate_memory() ModuleBase::timer::end("DFPT_PW_Data", "deallocate_memory"); } +void DFPT_PW_Data::set_drhoc_g(int q_idx, const std::vector>& rho_c) +{ + ModuleBase::TITLE("DFPT_PW_Data", "set_drhoc_g"); + ModuleBase::timer::start("DFPT_PW_Data", "set_drhoc_g"); + if (q_idx < 0) + { + ModuleBase::timer::end("DFPT_PW_Data", "set_drhoc_g"); + return; + } + if (q_idx >= static_cast(drhoc_g_.size())) + { + drhoc_g_.resize(q_idx + 1); + } + drhoc_g_[q_idx] = rho_c; + ModuleBase::timer::end("DFPT_PW_Data", "set_drhoc_g"); +} + +std::vector> DFPT_PW_Data::get_drhoc_g(int q_idx) const +{ + ModuleBase::TITLE("DFPT_PW_Data", "get_drhoc_g"); + ModuleBase::timer::start("DFPT_PW_Data", "get_drhoc_g"); + if (q_idx >= 0 && q_idx < static_cast(drhoc_g_.size())) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drhoc_g"); + return drhoc_g_[q_idx]; + } + ModuleBase::timer::end("DFPT_PW_Data", "get_drhoc_g"); + return {}; +} + } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h index d1ea8dd6f1f..239c77fa4ad 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.h +++ b/source/source_pw/module_dfpt/dfpt_pw_data.h @@ -62,6 +62,11 @@ class DFPT_PW_Data void set_drho_r(int q_idx, int spin, const std::vector& rho); std::vector get_drho_r(int q_idx, int spin) const; void set_drho_g(int q_idx, int spin, const std::vector>& rho); + + /// q-shifted first-order core-density amplitude (G space, npw); + /// zero-size means no nlcc atoms / not yet computed + void set_drhoc_g(int q_idx, const std::vector>& rho_c); + std::vector> get_drhoc_g(int q_idx) const; std::vector> get_drho_g(int q_idx, int spin) const; void set_dv_r(int q_idx, int spin, const std::vector& v); @@ -267,6 +272,8 @@ class DFPT_PW_Data std::vector>> drho_r_; std::vector>>> drho_g_; + ///< q-shifted first-order core density (G space), [q] + std::vector>> drhoc_g_; std::vector>> dv_r_; diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index a5f3590509d..df06df88dce 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -56,6 +56,19 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, xc_->set_context(q_cart); std::vector> a_r(nrxx); pw_rho_->recip2real(drho_in_g.data(), a_r.data()); + // nlcc: the first-order core density enters the XC channel only (the + // core Coulomb is already in the ionic pseudopotential, so the Hartree + // above correctly uses the valence density alone) + const std::vector> drhoc_g = data_.get_drhoc_g(last_q_); + if (static_cast(drhoc_g.size()) == pw_rho_->npw) + { + std::vector> drhoc_r(nrxx); + pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + a_r[ir] += drhoc_r[ir]; + } + } std::vector> b_r; xc_->apply(a_r, b_r); if (static_cast(b_r.size()) != nrxx) @@ -203,6 +216,9 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) const bool dbg = (getenv("DFPT_DEBUG") != nullptr); pert_.build_dv(q_idx, iat, idir, data_); + // nlcc: q-shifted core-density derivative for the displaced atom + // (fixed per displacement, consumed by assemble_v_sc's XC channel) + rho_.compute_drhoc(*ucell_, iat, idir, q_idx, data_); rho_.reset_mixing(q_idx); data_.set_drho_g(q_idx, 0, std::vector>(npw, std::complex(0.0, 0.0))); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 2a271c97168..239c7804ab8 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -6,6 +6,9 @@ #include "source_base/timer.h" #include "source_base/tool_title.h" #include "source_base/module_mixing/plain_mixing.h" +#include "source_base/math_integral.h" +#include "source_base/math_sphbes.h" +#include "source_estate/module_charge/chg_tools.h" #include #include @@ -242,6 +245,103 @@ void DFPT_Rho::compute_drho(const psi::Psi>& psi, ModuleBase::timer::end("DFPT_Rho", "compute_drho"); } +void DFPT_Rho::compute_drhoc(const UnitCell& ucell, + int iat, + int dir, + int q_idx, + DFPT_PW_Data& data) +{ + ModuleBase::TITLE("DFPT_Rho", "compute_drhoc"); + ModuleBase::timer::start("DFPT_Rho", "compute_drhoc"); + const int npw = pw_rho_ ? pw_rho_->npw : 0; + if (pw_rho_ == nullptr || npw == 0 || iat < 0 || iat >= ucell.nat || dir < 0 || dir > 2) + { + data.set_drhoc_g(q_idx, {}); + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); + return; + } + const int it = ucell.iat2it[iat]; + if (!ucell.atoms[it].ncpp.nlcc) + { + data.set_drhoc_g(q_idx, {}); + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); + return; + } + + // q-shifted |G+q|^2 in tpiba^2 units, and the q-shifted structure + // factor phase exp(i (G+q) tau_iat) for the displaced atom + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; + // tau in cartesian bohr (lat0-scaled: tau is in direct units of latvec) + const ModuleBase::Vector3 tau_cart = ucell.atoms[it].tau[ucell.iat2ia[iat]] * ucell.lat0; + + // build the |G+q|^2 shell values and evaluate the radial core FT + // on them (non_linear_core_correction works on unique gg shells) + std::vector w2_all(npw); + for (int ig = 0; ig < npw; ++ig) + { + const ModuleBase::Vector3 w = pw_rho_->gcar[ig] + q_cart; + w2_all[ig] = w * w; // dimensionless (1/lat0)^2 + } + + // use the existing radial FT machinery: pass the |G+q|^2 values as + // pseudo-unique shells (the function just needs matching arrays; + // we can call it per-shell or build a shell lookup) + // Simpler: directly evaluate the spherical Bessel integral here + // for each ig (no shell deduplication -- the cost is O(npw * mesh) + // and only runs once per displacement for nlcc atoms) + const pseudo& p = ucell.atoms[it].ncpp; + const int mesh = p.msh; + + std::vector> drhoc_g(npw, std::complex(0.0, 0.0)); + std::vector aux(mesh); + + for (int ig = 0; ig < npw; ++ig) + { + const ModuleBase::Vector3 w = pw_rho_->gcar[ig] + q_cart; + const double w2 = w * w; + // skip the |G+q| = 0 shell (dropped by compute_drho, zero density there) + if (w2 < 1.0e-12) + { + continue; + } + const double g = std::sqrt(w2) * ucell.tpiba; // physical 1/bohr + // radial FT: rho_core~(|G+q|) = 4pi/Omega * Int r^2 rho_atc(r) j0(gr) dr + ModuleBase::Sphbes::Spherical_Bessel(mesh, p.r.data(), g, 0, aux.data()); + double rhocg = 0.0; + for (int ir = 0; ir < mesh; ++ir) + { + aux[ir] = p.r[ir] * p.r[ir] * p.rho_atc[ir] * aux[ir]; + } + ModuleBase::Integral::Simpson_Integral(mesh, aux.data(), p.rab.data(), rhocg); + rhocg *= ModuleBase::FOUR_PI / ucell.omega; + + // structure factor phase: exp(i (G+q) . tau_iat) + const double phase = ModuleBase::TWO_PI * (w * tau_cart) / ucell.lat0; + const std::complex gtau(std::cos(phase), std::sin(phase)); + + // direction derivative: -i * tpiba * (G+q)_dir + const double wdir = ucell.tpiba * w[dir]; + + drhoc_g[ig] = std::complex(0.0, -wdir) * rhocg * gtau; + } + + data.set_drhoc_g(q_idx, drhoc_g); + + // report magnitude for diagnostics + double norm2 = 0.0; + for (int ig = 0; ig < npw; ++ig) + { + norm2 += std::norm(drhoc_g[ig]); + } + std::cout << " DFPT nlcc core derivative: atom " << iat + 1 << "(" + << ucell.atoms[it].label << ") dir " << static_cast('x' + dir) + << " |drhoc_g| = " << std::setprecision(3) << std::scientific << std::sqrt(norm2) + << std::defaultfloat << std::setprecision(6) << std::endl; + + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); +} + void DFPT_Rho::cal_docc(const psi::Psi>& psi, const ModuleBase::matrix& wg, int q_idx, diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 2d4876f3ace..83ae0cda250 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -5,6 +5,7 @@ #include "source_base/matrix3.h" #include "source_basis/module_pw/pw_basis.h" #include "source_basis/module_pw/pw_basis_k.h" +#include "source_cell/unitcell.h" #include "source_psi/psi.h" #include @@ -114,6 +115,22 @@ class DFPT_Rho int q_idx, DFPT_PW_Data& data); + /// q-shifted first-order core-density amplitude for an nlcc + /// pseudopotential atom displaced along (iat, dir), in G space: + /// d rho_core/d tau_dir (G+q) = -i tpiba (G+q)_dir + /// * rho_core~(|G+q|) * exp(i (G+q) tau_iat) + /// where rho_core~ is the radial spherical-Bessel FT of the core + /// charge (the same non_linear_core_correction used by the ground + /// state). Non-nlcc atoms / types give zero. The result enters the + /// XC channel of assemble_v_sc only (the core Coulomb is already in + /// the ionic pseudopotential, so no Hartree contribution). + /// ucell is passed explicitly (DFPT_Rho does not own it). + void compute_drhoc(const UnitCell& ucell, + int iat, + int dir, + int q_idx, + DFPT_PW_Data& data); + /// first-order occupation matrix (docc) for DFT+U (U0 reservation). void cal_docc(const psi::Psi>& psi, const ModuleBase::matrix& wg,