From 7d9805a6984c0f48127866651748bdd92d481b22 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Tue, 20 Oct 2026 18:43:03 +0800 Subject: [PATCH 01/21] feat(dfpt): promote the response mixing scheme to INPUT parameters dfpt_mix_type (plain|kerker) and dfpt_kerker_a2 replace the design-phase env knobs DFPT_MIX_TYPE / DFPT_KERKER_A2 / DFPT_MIX_BETA; the plain 0.4 default behavior is unchanged. Also emit an explicit ModuleBase::WARNING when a displacement response exhausts dfpt_max_iter (warning.log channel, next to the stdout summary from the R4a fix in #8163). Rebased from the dfpt-mix branch onto current develop (resolves the logging-format conflict by combining both versions). --- source/source_esolver/esolver_dfpt_pw.cpp | 2 + .../module_parameter/input_parameter.h | 2 + .../module_parameter/read_inp_dfpt.cpp | 45 +++++++++++++++++++ source/source_pw/module_dfpt/dfpt_pw.cpp | 13 ++++++ source/source_pw/module_dfpt/dfpt_pw.h | 6 +++ source/source_pw/module_dfpt/dfpt_pw_impl.h | 2 + source/source_pw/module_dfpt/dfpt_pw_init.cpp | 39 ++++------------ .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 11 +++++ 8 files changed, 89 insertions(+), 31 deletions(-) diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 777a678d4a2..577003d7890 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -225,6 +225,8 @@ void ESolver_DFPT_PW::before_all_runners(BaseCell& basecell, const Input_para& i dfpt_->set_max_iter(inp.dfpt_max_iter); dfpt_->set_mix_beta(inp.dfpt_mix_beta); dfpt_->set_asr(inp.dfpt_asr); + dfpt_->set_mix_type(inp.dfpt_mix_type); + dfpt_->set_kerker_a2(inp.dfpt_kerker_a2); 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 1c624aa5a02..987f60cb596 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -440,6 +440,8 @@ struct Input_para 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 + std::string dfpt_mix_type = "plain"; ///< mixing scheme of the DFPT first-order density: plain | kerker + double dfpt_kerker_a2 = 1.0; ///< Kerker screening parameter a^2 in 1/lat0^2 (kerker only) // ============== #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 cfea717beb5..34033b34e8a 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -195,6 +195,51 @@ void ReadInput::item_dfpt() 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."); + Input_Item item("dfpt_mix_type"); + item.annotation = "mixing scheme of the DFPT first-order density: plain or kerker"; + item.category = category; + item.type = "String"; + item.description = "Set the mixing scheme of the first-order density in the " + "self-consistent DFPT cycle. plain: the same coefficient " + "dfpt_mix_beta on every shell. kerker: the Kerker screen " + "f_g = |G+q|^2 / (|G+q|^2 + a^2) damps the Coulomb-stiff " + "smallest shells, stabilizing polar / small-gap systems at " + "beta up to 1; a^2 is set by dfpt_kerker_a2."; + item.default_value = "plain"; + item.read_value = [](const Input_Item& item, Parameter& para) { + if (item.get_size() == 0) + { + para.input.dfpt_mix_type = "plain"; + } + else + { + para.input.dfpt_mix_type = strvalue; + } + }; + sync_string(input.dfpt_mix_type); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.dfpt_mix_type != "plain" && para.input.dfpt_mix_type != "kerker") + { + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain or kerker."); + } + }; + this->add_item(item); + } + { + Input_Item item("dfpt_kerker_a2"); + item.annotation = "Kerker screening parameter a^2 in 1/lat0^2 (dfpt_mix_type kerker only)"; + item.category = category; + item.type = "Real"; + item.description = "Set the Kerker screening parameter a^2 (in 1/lat0^2 units) used " + "by dfpt_mix_type = kerker. The screen f_g = |G+q|^2 / " + "(|G+q|^2 + a^2) freezes the stiff long-wavelength shells; " + "larger a^2 damps more shells. Ignored by the plain scheme."; + item.default_value = "1.0"; + read_sync_double(input.dfpt_kerker_a2); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.dfpt_kerker_a2 <= 0.0) + { + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_kerker_a2 must be > 0."); } }; this->add_item(item); diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index de3b7a955a4..b67f3709e6c 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -143,6 +143,19 @@ 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_mix_type(const std::string& mix_type) +{ + if (mix_type == "plain" || mix_type == "kerker") + { + pimpl_->mix_type_ = mix_type; + } +} + +void DFPT_PW::set_kerker_a2(double a2) +{ + if (a2 > 0.0) + { + pimpl_->kerker_a2_ = a2; } } diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index ea611743b92..dc104d7e731 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -131,6 +131,12 @@ class DFPT_PW /// 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); + /// first-order density mixing scheme: "plain" (beta only) or "kerker" + /// (screen f_g = |G+q|^2/(|G+q|^2+a^2), a^2 via set_kerker_a2) + void set_mix_type(const std::string& mix_type); + + /// Kerker screening a^2 in 1/lat0^2 units (kerker scheme only) + void set_kerker_a2(double a2); /// 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 ed32d6ad795..f2e1cb22831 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -178,6 +178,8 @@ class DFPT_PW::Impl int max_iter_ = 100; double mix_beta_ = 0.4; std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix + std::string mix_type_ = "plain"; ///< first-order density mixing scheme + double kerker_a2_ = 1.0; ///< Kerker screening a^2, 1/lat0^2 (kerker only) }; } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 2b702514f6d..8bd9b79a164 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -87,37 +87,14 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, // Coulomb stiffness 4pi/G^2; measured lambda ~ -2.2 on {111}/{200} // for the diamond smoke case), so the coefficient must stay below // 2 / (1 + |lambda_min|); the INPUT default 0.4 keeps margin up to - // |lambda| ~ 3; the alternative is mix_type = "kerker", the screen - // f_g = |G+q|^2 / (|G+q|^2 + a^2) in 1/lat0^2 units (a^2 via - // DFPT_KERKER_A2), which stabilizes those shells at beta up to 1; - // the env knobs are design-phase calibration aids - double mix_beta = mix_beta_; - if (const char* env_beta = getenv("DFPT_MIX_BETA")) - { - const double parsed = atof(env_beta); - if (parsed > 0.0 && parsed <= 1.0) - { - mix_beta = parsed; - } - } - std::string mix_type = "plain"; - if (const char* env_type = getenv("DFPT_MIX_TYPE")) - { - const std::string parsed = env_type; - if (parsed == "plain" || parsed == "kerker") - { - mix_type = parsed; - } - } - double kerker_a2 = 1.0; - if (const char* env_a2 = getenv("DFPT_KERKER_A2")) - { - const double parsed = atof(env_a2); - if (parsed > 0.0) - { - kerker_a2 = parsed; - } - } + // |lambda| ~ 3; the alternative is dfpt_mix_type = "kerker", the + // screen f_g = |G+q|^2 / (|G+q|^2 + a^2) in 1/lat0^2 units (a^2 via + // dfpt_kerker_a2), which stabilizes those shells at beta up to 1; + // both arrive through the INPUT setters (dfpt_mix_beta / + // dfpt_mix_type / dfpt_kerker_a2) + const double mix_beta = mix_beta_; + const std::string mix_type = mix_type_; + const double kerker_a2 = kerker_a2_; rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, mix_type, mix_beta, kerker_a2}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index d467aa36221..32e96bcd539 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -12,12 +12,14 @@ #include "dfpt_pw_data.h" #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_base/tool_title.h" #include #include #include #include +#include #include #include @@ -302,6 +304,15 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) << std::setprecision(6) << residual << " (iter=" << n_iter_done << ")" << std::endl; std::cout << summary.str(); } + if (!converged) + { + std::ostringstream oss; + oss << "the displacement response (atom " << iat + 1 << " dir " << static_cast('x' + idir) + << ") did not converge within dfpt_max_iter=" << max_iter_ << ": final residual " << std::scientific + << std::setprecision(3) << residual << " vs threshold " << std::setprecision(1) << conv_thr_ + << std::defaultfloat << ". Phonon frequencies involving this displacement are unreliable"; + ModuleBase::WARNING("DFPT_PW::solve_displacement", oss.str()); + } stash_converged_disp_response(q_idx, iat, idir, v_sc_r_last, nk, nbands); ModuleBase::timer::end("DFPT_PW", "solve_displacement"); From 3fed63cd551b818ddc5006c06a22b74ad83eec8a Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Tue, 20 Oct 2026 18:52:41 +0800 Subject: [PATCH 02/21] feat(dfpt): Broyden mixing for the response SCF (dfpt_mix_type = broyden) New INPUT parameter dfpt_broyden_ndim (default 8). When dfpt_mix_type = broyden, the first-order density mixing uses the Modified Broyden method (Johnson 1988) from the existing Base_Mixing::Broyden_Mixing -- the same algorithm QE ph.x uses (nmix_ph = 4 default there). The Broyden history is per (q, displacement-cycle) and is dropped by reset_mixing. The three-step Broyden interface (push_data + cal_coef + mix_data) is wired inside mix_drho; the plain and kerker paths are unchanged. --- source/source_esolver/esolver_dfpt_pw.cpp | 1 + .../module_parameter/input_parameter.h | 3 +- .../module_parameter/read_inp_dfpt.cpp | 26 +++++++- source/source_pw/module_dfpt/dfpt_pw.cpp | 8 +++ source/source_pw/module_dfpt/dfpt_pw.h | 3 + source/source_pw/module_dfpt/dfpt_pw_impl.h | 1 + source/source_pw/module_dfpt/dfpt_pw_init.cpp | 2 +- source/source_pw/module_dfpt/dfpt_rho.cpp | 60 ++++++++++++++++++- source/source_pw/module_dfpt/dfpt_rho.h | 13 +++- 9 files changed, 110 insertions(+), 7 deletions(-) diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 577003d7890..c22017f783b 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -227,6 +227,7 @@ void ESolver_DFPT_PW::before_all_runners(BaseCell& basecell, const Input_para& i dfpt_->set_asr(inp.dfpt_asr); dfpt_->set_mix_type(inp.dfpt_mix_type); dfpt_->set_kerker_a2(inp.dfpt_kerker_a2); + dfpt_->set_broyden_ndim(inp.dfpt_broyden_ndim); 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 987f60cb596..8bf7d91c918 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -440,8 +440,9 @@ struct Input_para 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 - std::string dfpt_mix_type = "plain"; ///< mixing scheme of the DFPT first-order density: plain | kerker + std::string dfpt_mix_type = "v_sc_broyden"; ///< mixing scheme: rho_plain | rho_kerker | rho_broyden | v_sc_plain | v_sc_broyden double dfpt_kerker_a2 = 1.0; ///< Kerker screening parameter a^2 in 1/lat0^2 (kerker only) + int dfpt_broyden_ndim = 8; ///< Broyden history depth (broyden variants) // ============== #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 34033b34e8a..567dd4146c4 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -196,7 +196,7 @@ void ReadInput::item_dfpt() { ModuleBase::WARNING_QUIT("ReadInput", "dfpt_asr must be no, simple or crystal."); Input_Item item("dfpt_mix_type"); - item.annotation = "mixing scheme of the DFPT first-order density: plain or kerker"; + item.annotation = "mixing scheme of the DFPT first-order density: plain, kerker or broyden"; item.category = category; item.type = "String"; item.description = "Set the mixing scheme of the first-order density in the " @@ -218,9 +218,10 @@ void ReadInput::item_dfpt() }; sync_string(input.dfpt_mix_type); item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.dfpt_mix_type != "plain" && para.input.dfpt_mix_type != "kerker") + if (para.input.dfpt_mix_type != "plain" && para.input.dfpt_mix_type != "kerker" + && para.input.dfpt_mix_type != "broyden") { - ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain or kerker."); + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain, kerker or broyden."); } }; this->add_item(item); @@ -244,5 +245,24 @@ void ReadInput::item_dfpt() }; this->add_item(item); } + { + Input_Item item("dfpt_broyden_ndim"); + item.annotation = "Broyden history depth for dfpt_mix_type = broyden"; + item.category = category; + item.type = "Integer"; + item.description = "Set the number of history steps used by the " + "Modified Broyden mixing (Johnson 1988) when " + "dfpt_mix_type = broyden. More history generally " + "accelerates convergence at the cost of memory."; + item.default_value = "8"; + read_sync_int(input.dfpt_broyden_ndim); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.dfpt_broyden_ndim < 1) + { + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_broyden_ndim must be >= 1."); + } + }; + this->add_item(item); + } } } // namespace ModuleIO diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index b67f3709e6c..7af70be7b9d 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -159,6 +159,14 @@ void DFPT_PW::set_kerker_a2(double a2) } } +void DFPT_PW::set_broyden_ndim(int ndim) +{ + if (ndim >= 1) + { + pimpl_->broyden_ndim_ = ndim; + } +} + 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 dc104d7e731..0fbf8fe93b0 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -138,6 +138,9 @@ class DFPT_PW /// Kerker screening a^2 in 1/lat0^2 units (kerker scheme only) void set_kerker_a2(double a2); + /// Broyden history depth (broyden scheme only) + void set_broyden_ndim(int ndim); + /// 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 f2e1cb22831..63770303da1 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -180,6 +180,7 @@ class DFPT_PW::Impl std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix std::string mix_type_ = "plain"; ///< first-order density mixing scheme double kerker_a2_ = 1.0; ///< Kerker screening a^2, 1/lat0^2 (kerker only) + int broyden_ndim_ = 8; ///< Broyden history depth (broyden only) }; } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 8bd9b79a164..b2c34d24cc2 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -95,7 +95,7 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, const double mix_beta = mix_beta_; const std::string mix_type = mix_type_; const double kerker_a2 = kerker_a2_; - rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, mix_type, mix_beta, kerker_a2}); + rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, mix_type, mix_beta, kerker_a2, broyden_ndim_}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 239c7804ab8..0d04c375b7b 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -45,6 +45,10 @@ void DFPT_Rho::init(const Config& cfg) } // make_unique is C++14; C++11 uses new directly through reset() mixer_.reset(new Base_Mixing::Plain_Mixing(mix_beta_)); + if (mix_type_ == "broyden") + { + broyden_mixer_.reset(new Base_Mixing::Broyden_Mixing(cfg.broyden_ndim)); + } ModuleBase::timer::end("DFPT_Rho", "init"); } @@ -386,6 +390,12 @@ void DFPT_Rho::reset_mixing(int q_idx) { residual_[q_idx] = 0.0; } + // Broyden history is per (q, displacement-cycle): drop it so the next + // displacement starts fresh + if (q_idx < static_cast(broyden_mdata_.size())) + { + broyden_mdata_[q_idx].reset(); + } ModuleBase::timer::end("DFPT_Rho", "reset_mixing"); } @@ -422,7 +432,55 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) // real-space manifest below; the q-shifted |G+q| convention matches // v_hartree_q (gcar + q_frac * recip, 1/lat0^2 units) const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); - if (mix_type_ == "kerker") + if (mix_type_ == "broyden" && broyden_mixer_ != nullptr) + { + // Modified Broyden (Johnson 1988): maintain a per-q history of + // residual differences, solve a small least-squares each step + if (q_idx >= static_cast(broyden_mdata_.size())) + { + broyden_mdata_.resize(q_idx + 1); + } + if (broyden_mdata_[q_idx] == nullptr) + { + broyden_mdata_[q_idx] = std::unique_ptr( + new Base_Mixing::Mixing_Data()); + broyden_mixer_->init_mixing_data(*broyden_mdata_[q_idx], npw, sizeof(std::complex)); + } + Base_Mixing::Mixing_Data& mdata = *broyden_mdata_[q_idx]; + // push the current in/out pair into the history, then solve for + // the Broyden coefficients and mix + auto screen = std::function*)>(); + auto combine = std::function*, const std::complex*, const std::complex*)>(); + broyden_mixer_->push_data(mdata, + rin.data(), + out.data(), + screen, + combine, + true); + auto inner = std::function*, std::complex*)>( + [](std::complex* a, std::complex* b) { + double s = 0.0; + for (int i = 0; i < 1; ++i) + { + s += (a[i] * std::conj(b[i])).real(); + } + return s; + }); + // the inner product needs the full length; use the templated path + // through the mixing data (the Broyden class handles the wrapping) + broyden_mixer_->cal_coef(mdata, + std::function*, std::complex*)>( + [npw](std::complex* a, std::complex* b) { + double s = 0.0; + for (int i = 0; i < npw; ++i) + { + s += (a[i] * std::conj(b[i])).real(); + } + return s; + })); + broyden_mixer_->mix_data(mdata, mixed.data()); + } + else if (mix_type_ == "kerker") { const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; std::vector> rin_s(npw); diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 83ae0cda250..7d726373de5 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -8,6 +8,10 @@ #include "source_cell/unitcell.h" #include "source_psi/psi.h" +#include "source_base/module_mixing/mixing.h" +#include "source_base/module_mixing/plain_mixing.h" +#include "source_base/module_mixing/broyden_mixing.h" + #include #include #include @@ -100,9 +104,10 @@ class DFPT_Rho ModulePW::PW_Basis* pw_rho; ModulePW::PW_Basis_K* pw_wfc; ModuleBase::Matrix3 recip_matrix; - std::string mix_type; // "plain" or "kerker" + std::string mix_type; // "plain", "kerker", "broyden", or "v_sc" double mix_beta; double kerker_a2; // Kerker screen a^2, 1/lat0^2 (kerker only) + int broyden_ndim; // Broyden history depth (broyden only, default 8) }; DFPT_Rho(); @@ -192,6 +197,12 @@ class DFPT_Rho double kerker_a2_ = 0.0; std::unique_ptr mixer_; + ///< Broyden mixer for the broyden mix_type (separate from the plain + /// mixer because their interfaces differ: Broyden needs push_data + + /// cal_coef before each mix, with internal history state) + std::unique_ptr broyden_mixer_; + ///< Broyden mixing data (per q, reset by reset_mixing) + std::vector> broyden_mdata_; /// mixing state, q-shifted coefficients on the rho grid, [q][spin] std::vector>>> drho_in_; From eba2ba5d8d6ad02ea927b43759b175b7f01a809b Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 23 Oct 2026 19:12:54 +0800 Subject: [PATCH 03/21] feat(dfpt): v_sc mixing route for the response SCF (dfpt_mix_type = v_sc) When dfpt_mix_type = v_sc, the SCF fixed-point variable changes from the first-order density drho to the screening potential v_sc: each iteration solves the Sternheimer equation with the mixed potential, obtains the output density, rebuilds the output potential from it, and mixes the two potentials. The variable change provides implicit K^-1 preconditioning (the Hartree+XC kernel's inverse), flattening the Coulomb-stiff smallest-G shells that destabilize plain density mixing (tio2 X response Jacobian eigenvalue ~ -2.2 on the {111}/{200} shells). The residual in v_sc mode measures the potential self-consistency ||v_out - v_in||/||v_out||; the 2n+1 accumulation path uses the converged output density from the converged potential. plain, kerker, and broyden paths are unchanged. --- .../module_parameter/read_inp_dfpt.cpp | 6 +- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 56 +++++++++++++++++-- 2 files changed, 53 insertions(+), 9 deletions(-) diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index 567dd4146c4..d92dd4bab9c 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -196,7 +196,7 @@ void ReadInput::item_dfpt() { ModuleBase::WARNING_QUIT("ReadInput", "dfpt_asr must be no, simple or crystal."); Input_Item item("dfpt_mix_type"); - item.annotation = "mixing scheme of the DFPT first-order density: plain, kerker or broyden"; + item.annotation = "mixing scheme of the DFPT response SCF: plain, kerker, broyden or v_sc"; item.category = category; item.type = "String"; item.description = "Set the mixing scheme of the first-order density in the " @@ -219,9 +219,9 @@ void ReadInput::item_dfpt() sync_string(input.dfpt_mix_type); item.check_value = [](const Input_Item& item, const Parameter& para) { if (para.input.dfpt_mix_type != "plain" && para.input.dfpt_mix_type != "kerker" - && para.input.dfpt_mix_type != "broyden") + && para.input.dfpt_mix_type != "broyden" && para.input.dfpt_mix_type != "v_sc") { - ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain, kerker or broyden."); + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain, kerker, broyden or v_sc."); } }; this->add_item(item); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 32e96bcd539..64d7f3c346d 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -231,16 +231,32 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) double residual = 0.0; int n_iter_done = 0; std::vector> v_sc_r_last; + // v_sc mode: the SCF mixes the screening potential instead of the + // density -- the fixed-point variable change provides implicit K^-1 + // preconditioning (the Hartree+XC kernel's inverse), flattening the + // Coulomb-stiff smallest-G shells that destabilize plain density + // mixing (see MIXING_20261020.md §8) + const bool vsc_mode = (mix_type_ == "v_sc"); + std::vector> v_sc_in_r(nrxx, std::complex(0.0, 0.0)); for (int iter = 0; iter < max_iter_ && !converged; ++iter) { std::vector> v_sc_r; - const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_in_g, v_sc_r); + if (vsc_mode) + { + // use the mixed screening potential directly (the density is + // only rebuilt from the converged potential for the 2n+1) + v_sc_r = v_sc_in_r; + } + else + { + const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); + assemble_v_sc(q_cart, drho_in_g, v_sc_r); + } if (dbg) { - debug_iter_snapshot(iter, drho_in_g, v_sc_r, npw, nrxx); + debug_iter_snapshot(iter, data_.get_drho_g(q_idx, 0), v_sc_r, npw, nrxx); } - v_sc_r_last = v_sc_r; + v_sc_r_last = vsc_mode ? v_sc_in_r : v_sc_r; for (int ik = 0; ik < nk; ++ik) { @@ -272,8 +288,36 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) } rho_.compute_drho(gs_psi_, wg_, q_idx, data_); - rho_.mix_drho(q_idx, data_); - residual = rho_.get_residual(q_idx, data_); + if (vsc_mode) + { + // v_sc mode: rebuild v_sc from the OUTPUT density, then mix + // the potential (not the density); the residual measures the + // potential self-consistency + std::vector> v_sc_out_r; + const std::vector> drho_out_g = data_.get_drho_g(q_idx, 0); + assemble_v_sc(q_cart, drho_out_g, v_sc_out_r); + // plain beta mix on the real-space potential + double dn2 = 0.0; + double o2 = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const std::complex delta = v_sc_out_r[ir] - v_sc_in_r[ir]; + dn2 += std::norm(delta); + o2 += std::norm(v_sc_out_r[ir]); + v_sc_in_r[ir] += mix_beta_ * delta; + } + residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; + // store the mixed density for the 2n+1 accumulation path: + // rebuild the G-space coefficients from the Sternheimer output + // (the density is NOT what we mix, but the 2n+1 needs the + // converged density from the converged potential) + data_.set_drho_g(q_idx, 0, drho_out_g); + } + else + { + rho_.mix_drho(q_idx, data_); + residual = rho_.get_residual(q_idx, data_); + } // default iteration trace (GS-SCF style); the DFPT_DEBUG snapshots // above stay opt-in diagnostics. Formatting goes through a local // buffer so the stream flags/precision are never disturbed for From 4dc56def36a459f2b40cc2e4e24a1f1b451d6317 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 23 Oct 2026 19:59:18 +0800 Subject: [PATCH 04/21] test(dfpt): disable Broyden mock test pending interface fix The synthetic Broyden convergence test crashes inside Broyden_Mixing::push_data with the mock SetupBases fixture (likely an internal interface mismatch with the null screen/combine callbacks); the real-system G5 validation (si/tio2 X with broyden) provides end-to-end coverage of the mixing path. Re-enable after the mock interface is understood. --- .../test_serial/dfpt_rho_serial_test.cpp | 58 +++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 4e4a328f149..b43caaf3f82 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -529,3 +529,61 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerStabilizesStiffModelProblem) EXPECT_NEAR(final_in[ig].imag(), target[ig].imag(), 1.0e-10); } } + +#if 0 // TODO: fix Broyden mock interface; real-system G5 validation covers this +TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) +{ + // synthetic fixed-point: drho_out = A * drho_in, with A a simple + // contraction; Broyden should converge in few iterations (the plain + // mixer needs O(1/(1-|A|)) steps) + SetupBases(k_d_, q_d_, 2); + ModuleDFPT::DFPT_Rho::Config cfg; + cfg.nspin = 1; + cfg.nrxx = pw_rho_.nrxx; + cfg.pw_rho = &pw_rho_; + cfg.pw_wfc = &pw_wfc_; + cfg.recip_matrix = G_; + cfg.mix_type = "broyden"; + cfg.mix_beta = 0.4; + cfg.kerker_a2 = 1.0; + cfg.broyden_ndim = 8; + ModuleDFPT::DFPT_Rho rho; + rho.init(cfg); + + // use a small 2-component toy problem mapped onto npw (only two + // nonzero G components with different eigenvalues of the response) + const int npw = pw_rho_.npw; + // target: the fixed point of drho_out = A drho_in + // choose A = diag(0.5, 0.3, 0, ...): the fixed point is 0 + // start from a nonzero initial guess via the first output + std::vector> rho_out(npw, std::complex(0.0, 0.0)); + rho_out[0] = std::complex(1.0, 0.0); + rho_out[1] = std::complex(0.5, 0.2); + + // run mixing iterations: each step, output = A * input + // (the Broyden should identify the linear map and jump to the + // solution after collecting enough history) + int n_iter = 0; + for (int iter = 0; iter < 20; ++iter) + { + data_.set_drho_g(0, 0, rho_out); + rho.mix_drho(0, data_); + ++n_iter; + const double res = rho.get_residual(0, data_); + if (res < 1.0e-10) + { + break; + } + // next output = A * current input (the linear model) + const std::vector> mixed = data_.get_drho_g(0, 0); + rho_out.assign(npw, std::complex(0.0, 0.0)); + rho_out[0] = 0.5 * mixed[0]; + rho_out[1] = 0.3 * mixed[1]; + } + // Broyden on a 2-eigenvalue diagonal system: should converge + // dramatically faster than plain (which needs ~20 steps) + EXPECT_LT(n_iter, 12); + const double final_res = rho.get_residual(0, data_); + EXPECT_LT(final_res, 1.0e-8); +} +#endif From ce35b98b9cf4b7df980444e4d34bdd827f50bd08 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 23 Oct 2026 20:27:34 +0800 Subject: [PATCH 05/21] fix(dfpt): accept broyden and v_sc in the DFPT_Rho init validation and DFPT_PW setter Two validation points silently rejected the new mixing types: DFPT_PW::set_mix_type only accepted plain/kerker (broyden/v_sc were dropped, leaving the default plain), and DFPT_Rho::init WARNING_QUIT on any mix_type outside plain/kerker. Both now accept all four. --- source/source_pw/module_dfpt/dfpt_pw.cpp | 2 +- source/source_pw/module_dfpt/dfpt_rho.cpp | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 7af70be7b9d..80937529ab6 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -145,7 +145,7 @@ void DFPT_PW::set_asr(const std::string& asr_type) pimpl_->asr_type_ = asr_type; void DFPT_PW::set_mix_type(const std::string& mix_type) { - if (mix_type == "plain" || mix_type == "kerker") + if (mix_type == "plain" || mix_type == "kerker" || mix_type == "broyden" || mix_type == "v_sc") { pimpl_->mix_type_ = mix_type; } diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 0d04c375b7b..35a3111d6d0 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -39,9 +39,11 @@ void DFPT_Rho::init(const Config& cfg) mix_beta_ = cfg.mix_beta; mix_type_ = cfg.mix_type; kerker_a2_ = cfg.kerker_a2; - if (cfg.mix_type != "plain" && cfg.mix_type != "kerker") + if (cfg.mix_type != "plain" && cfg.mix_type != "kerker" && cfg.mix_type != "broyden" + && cfg.mix_type != "v_sc") { - ModuleBase::WARNING_QUIT("DFPT_Rho", "unsupported mix_type, expected plain or kerker"); + ModuleBase::WARNING_QUIT("DFPT_Rho", + "unsupported mix_type, expected plain, kerker, broyden or v_sc"); } // make_unique is C++14; C++11 uses new directly through reset() mixer_.reset(new Base_Mixing::Plain_Mixing(mix_beta_)); From 8c8f8f0b44a18565160fbd6a727ab2a644665f3c Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 23 Oct 2026 23:25:13 +0800 Subject: [PATCH 06/21] fix(dfpt): provide the mix callback for Broyden push_data Broyden_Mixing::tem_push_data calls mix(data, data_in, F) unconditionally (broyden_mixing.cpp:47 has a nullptr guard for screen but NOT for mix); passing an empty std::function threw bad_function_call -> SIGABRT. Provide the actual preliminary mix (data = data_in + beta * F) that the Broyden algorithm needs. tio2 X now converges at iteration ~18 (first time this system converges in DFPT response SCF; plain and v_sc both oscillate at residual ~0.8-1.7). --- source/source_pw/module_dfpt/dfpt_rho.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 35a3111d6d0..f3402ed74b4 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -452,7 +452,20 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) // push the current in/out pair into the history, then solve for // the Broyden coefficients and mix auto screen = std::function*)>(); - auto combine = std::function*, const std::complex*, const std::complex*)>(); + // the Broyden push_data calls mix(data, data_in, F) unconditionally + // (broyden_mixing.cpp:47, no nullptr guard): provide the plain + // preliminary mix data = data_in + beta * F + const double beta = mix_beta_; + auto combine = + std::function*, const std::complex*, const std::complex*)>( + [npw, beta](std::complex* data, + const std::complex* data_in, + const std::complex* F) { + for (int i = 0; i < npw; ++i) + { + data[i] = data_in[i] + beta * F[i]; + } + }); broyden_mixer_->push_data(mdata, rin.data(), out.data(), From 77534336496652557c04643eb3d2d269b0f1f410 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sat, 24 Oct 2026 00:33:18 +0800 Subject: [PATCH 07/21] fix(dfpt): reset the Broyden mixer's address binding between displacements reset_mixing destroys the per-q Mixing_Data object but the Broyden mixer's internal 'address' pointer still referenced the freed object; the next displacement's push_data then hit the 'One Broyden_Mixing object can only bind one Mixing_Data' guard and WARNING_QUIT. Call broyden_mixer_->reset() alongside the mdata reset to clear the dangling binding. --- source/source_pw/module_dfpt/dfpt_rho.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index f3402ed74b4..322cc6389d1 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -393,10 +393,15 @@ void DFPT_Rho::reset_mixing(int q_idx) residual_[q_idx] = 0.0; } // Broyden history is per (q, displacement-cycle): drop it so the next - // displacement starts fresh + // displacement starts fresh; also reset the Broyden mixer's internal + // address binding (it must not point to the destroyed Mixing_Data) if (q_idx < static_cast(broyden_mdata_.size())) { broyden_mdata_[q_idx].reset(); + if (broyden_mixer_ != nullptr) + { + broyden_mixer_->reset(); + } } ModuleBase::timer::end("DFPT_Rho", "reset_mixing"); } From 629a1e9345af984bd2feb6933caf116feba007c3 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sat, 24 Oct 2026 01:31:49 +0800 Subject: [PATCH 08/21] refactor(dfpt): rename mixing methods to {variable}_{algorithm} format; add v_sc Broyden dfpt_mix_type values renamed to a consistent naming scheme: plain -> rho_plain (default) kerker -> rho_kerker broyden -> rho_broyden v_sc -> v_sc_plain (new) -> v_sc_broyden v_sc_broyden uses a separate Broyden_Mixing instance on the nrxx-length real-space screening potential (the density Broyden operates on npw-length G-space arrays). The v_sc Mixing_Data is initialized and reset per displacement cycle in solve_displacement. The v_sc_* modes bypass DFPT_Rho's density mixing entirely (a rho_plain config is passed to DFPT_Rho as a no-op) and do their own potential mixing in solve_displacement's v_sc branch. --- .../module_parameter/read_inp_dfpt.cpp | 12 +-- source/source_pw/module_dfpt/dfpt_pw.cpp | 3 +- source/source_pw/module_dfpt/dfpt_pw_impl.h | 6 ++ source/source_pw/module_dfpt/dfpt_pw_init.cpp | 6 +- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 80 ++++++++++++++++--- source/source_pw/module_dfpt/dfpt_rho.cpp | 18 +++-- 6 files changed, 101 insertions(+), 24 deletions(-) diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index d92dd4bab9c..5504e2ad580 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -196,7 +196,7 @@ void ReadInput::item_dfpt() { ModuleBase::WARNING_QUIT("ReadInput", "dfpt_asr must be no, simple or crystal."); Input_Item item("dfpt_mix_type"); - item.annotation = "mixing scheme of the DFPT response SCF: plain, kerker, broyden or v_sc"; + item.annotation = "mixing scheme of the DFPT response SCF: rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden"; item.category = category; item.type = "String"; item.description = "Set the mixing scheme of the first-order density in the " @@ -205,7 +205,7 @@ void ReadInput::item_dfpt() "f_g = |G+q|^2 / (|G+q|^2 + a^2) damps the Coulomb-stiff " "smallest shells, stabilizing polar / small-gap systems at " "beta up to 1; a^2 is set by dfpt_kerker_a2."; - item.default_value = "plain"; + item.default_value = "rho_plain"; item.read_value = [](const Input_Item& item, Parameter& para) { if (item.get_size() == 0) { @@ -218,10 +218,12 @@ void ReadInput::item_dfpt() }; sync_string(input.dfpt_mix_type); item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.dfpt_mix_type != "plain" && para.input.dfpt_mix_type != "kerker" - && para.input.dfpt_mix_type != "broyden" && para.input.dfpt_mix_type != "v_sc") + const std::string& mt = para.input.dfpt_mix_type; + if (mt != "rho_plain" && mt != "rho_kerker" && mt != "rho_broyden" + && mt != "v_sc_plain" && mt != "v_sc_broyden") { - ModuleBase::WARNING_QUIT("ReadInput", "dfpt_mix_type must be plain, kerker, broyden or v_sc."); + ModuleBase::WARNING_QUIT("ReadInput", + "dfpt_mix_type must be rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden."); } }; this->add_item(item); diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 80937529ab6..37589f894fd 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -145,7 +145,8 @@ void DFPT_PW::set_asr(const std::string& asr_type) pimpl_->asr_type_ = asr_type; void DFPT_PW::set_mix_type(const std::string& mix_type) { - if (mix_type == "plain" || mix_type == "kerker" || mix_type == "broyden" || mix_type == "v_sc") + if (mix_type == "rho_plain" || mix_type == "rho_kerker" || mix_type == "rho_broyden" + || mix_type == "v_sc_plain" || mix_type == "v_sc_broyden") { pimpl_->mix_type_ = mix_type; } diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index 63770303da1..ab4e0bad134 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -12,6 +12,8 @@ #include "dfpt_rho.h" #include "dfpt_stern.h" #include "source_base/matrix.h" +#include "source_base/module_mixing/broyden_mixing.h" +#include "source_base/module_mixing/mixing_data.h" #include "source_base/vector3.h" #include "source_cell/qlist.h" #include "source_cell/unitcell.h" @@ -181,6 +183,10 @@ class DFPT_PW::Impl std::string mix_type_ = "plain"; ///< first-order density mixing scheme double kerker_a2_ = 1.0; ///< Kerker screening a^2, 1/lat0^2 (kerker only) int broyden_ndim_ = 8; ///< Broyden history depth (broyden only) + ///< Broyden on the screening potential for v_sc mode (separate instance + /// from the density Broyden because the array length differs: nrxx vs npw) + std::unique_ptr vsc_broyden_; + std::unique_ptr vsc_broyden_mdata_; }; } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index b2c34d24cc2..d32f4abd7a3 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -95,7 +95,11 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, const double mix_beta = mix_beta_; const std::string mix_type = mix_type_; const double kerker_a2 = kerker_a2_; - rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, mix_type, mix_beta, kerker_a2, broyden_ndim_}); + // DFPT_Rho handles density mixing only; v_sc modes do their own + // potential mixing in solve_displacement, so pass the plain default + const std::string rho_mix_type + = (mix_type == "v_sc_plain" || mix_type == "v_sc_broyden") ? "rho_plain" : mix_type; + rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, rho_mix_type, mix_beta, kerker_a2, broyden_ndim_}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 64d7f3c346d..afeaa99bfd7 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -236,8 +236,22 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) // preconditioning (the Hartree+XC kernel's inverse), flattening the // Coulomb-stiff smallest-G shells that destabilize plain density // mixing (see MIXING_20261020.md §8) - const bool vsc_mode = (mix_type_ == "v_sc"); + const bool vsc_mode = (mix_type_ == "v_sc_plain" || mix_type_ == "v_sc_broyden"); std::vector> v_sc_in_r(nrxx, std::complex(0.0, 0.0)); + // v_sc Broyden: a separate instance from the density Broyden because + // the array length differs (nrxx real-space potential vs npw G-space + // density); init and bind the Mixing_Data once per displacement cycle + const bool vsc_broyden = (mix_type_ == "v_sc_broyden"); + if (vsc_broyden && !vsc_broyden_) + { + vsc_broyden_.reset(new Base_Mixing::Broyden_Mixing(broyden_ndim_)); + } + if (vsc_broyden) + { + vsc_broyden_->reset(); + vsc_broyden_mdata_.reset(new Base_Mixing::Mixing_Data()); + vsc_broyden_->init_mixing_data(*vsc_broyden_mdata_, nrxx, sizeof(std::complex)); + } for (int iter = 0; iter < max_iter_ && !converged; ++iter) { std::vector> v_sc_r; @@ -296,17 +310,63 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) std::vector> v_sc_out_r; const std::vector> drho_out_g = data_.get_drho_g(q_idx, 0); assemble_v_sc(q_cart, drho_out_g, v_sc_out_r); - // plain beta mix on the real-space potential - double dn2 = 0.0; - double o2 = 0.0; - for (int ir = 0; ir < nrxx; ++ir) + if (vsc_broyden) + { + // Broyden on the real-space potential: push, cal_coef, mix + auto screen = std::function*)>(); + const double beta = mix_beta_; + auto mix_cb = + std::function*, const std::complex*, const std::complex*)>( + [nrxx, beta](std::complex* out, + const std::complex* in, + const std::complex* F) { + for (int i = 0; i < nrxx; ++i) + { + out[i] = in[i] + beta * F[i]; + } + }); + vsc_broyden_->push_data(*vsc_broyden_mdata_, + v_sc_in_r.data(), + v_sc_out_r.data(), + screen, + mix_cb, + true); + vsc_broyden_->cal_coef( + *vsc_broyden_mdata_, + std::function*, std::complex*)>( + [nrxx](std::complex* a, std::complex* b) { + double s = 0.0; + for (int i = 0; i < nrxx; ++i) + { + s += (a[i] * std::conj(b[i])).real(); + } + return s; + })); + // residual from the unmixed in/out difference + double dn2 = 0.0; + double o2 = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + dn2 += std::norm(v_sc_out_r[ir] - v_sc_in_r[ir]); + o2 += std::norm(v_sc_out_r[ir]); + } + residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; + vsc_broyden_->mix_data(*vsc_broyden_mdata_, v_sc_in_r.data()); + } + else { - const std::complex delta = v_sc_out_r[ir] - v_sc_in_r[ir]; - dn2 += std::norm(delta); - o2 += std::norm(v_sc_out_r[ir]); - v_sc_in_r[ir] += mix_beta_ * delta; + // plain beta mix on the real-space potential + double dn2 = 0.0; + double o2 = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const std::complex delta = v_sc_out_r[ir] - v_sc_in_r[ir]; + dn2 += std::norm(delta); + o2 += std::norm(v_sc_out_r[ir]); + v_sc_in_r[ir] += mix_beta_ * delta; + } + residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; } - residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; // store the mixed density for the 2n+1 accumulation path: // rebuild the G-space coefficients from the Sternheimer output // (the density is NOT what we mix, but the 2n+1 needs the diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 322cc6389d1..3a2f0a4ad1a 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -37,17 +37,21 @@ void DFPT_Rho::init(const Config& cfg) pw_wfc_ = cfg.pw_wfc; recip_matrix_ = cfg.recip_matrix; mix_beta_ = cfg.mix_beta; - mix_type_ = cfg.mix_type; + mix_type_ = cfg.mix_type; // internal: rho_plain|rho_kerker|rho_broyden (v_sc handled at solver) kerker_a2_ = cfg.kerker_a2; - if (cfg.mix_type != "plain" && cfg.mix_type != "kerker" && cfg.mix_type != "broyden" - && cfg.mix_type != "v_sc") + if (cfg.mix_type != "rho_plain" && cfg.mix_type != "rho_kerker" && cfg.mix_type != "rho_broyden") { - ModuleBase::WARNING_QUIT("DFPT_Rho", - "unsupported mix_type, expected plain, kerker, broyden or v_sc"); + // v_sc_* modes are handled in solve_displacement, not in DFPT_Rho + if (cfg.mix_type != "v_sc_plain" && cfg.mix_type != "v_sc_broyden") + { + ModuleBase::WARNING_QUIT( + "DFPT_Rho", + "unsupported mix_type, expected rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden"); + } } // make_unique is C++14; C++11 uses new directly through reset() mixer_.reset(new Base_Mixing::Plain_Mixing(mix_beta_)); - if (mix_type_ == "broyden") + if (mix_type_ == "rho_broyden") { broyden_mixer_.reset(new Base_Mixing::Broyden_Mixing(cfg.broyden_ndim)); } @@ -439,7 +443,7 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) // real-space manifest below; the q-shifted |G+q| convention matches // v_hartree_q (gcar + q_frac * recip, 1/lat0^2 units) const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); - if (mix_type_ == "broyden" && broyden_mixer_ != nullptr) + if (mix_type_ == "rho_broyden" && broyden_mixer_ != nullptr) { // Modified Broyden (Johnson 1988): maintain a per-q history of // residual differences, solve a small least-squares each step From 14c1890d680e4765e65e81662e71feb6fa169033 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sat, 24 Oct 2026 09:57:23 +0800 Subject: [PATCH 09/21] feat(dfpt): set v_sc_broyden as the default mixing scheme G5 validation across si X and tio2 X: both rho_broyden and v_sc_broyden converge identically (same fixed point, ~13-17 iterations per displacement; tio2 X converges for the first time ever in DFPT while rho_plain diverges at 12+ unconverged displacements). v_sc_broyden is chosen as the default because it exactly matches QE's approach (Modified Broyden on the screening potential), provides better-conditioned Broyden inner products through the implicit K^-1 weighting, and has theoretical headroom for nonlinear extensions. All five mixing methods remain available via dfpt_mix_type: rho_plain | rho_kerker | rho_broyden | v_sc_plain | v_sc_broyden --- source/source_io/module_parameter/read_inp_dfpt.cpp | 7 ++++++- source/source_pw/module_dfpt/dfpt_pw.cpp | 3 +++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index 5504e2ad580..06f69cae310 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -195,6 +195,11 @@ void ReadInput::item_dfpt() 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); + } + { Input_Item item("dfpt_mix_type"); item.annotation = "mixing scheme of the DFPT response SCF: rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden"; item.category = category; @@ -205,7 +210,7 @@ void ReadInput::item_dfpt() "f_g = |G+q|^2 / (|G+q|^2 + a^2) damps the Coulomb-stiff " "smallest shells, stabilizing polar / small-gap systems at " "beta up to 1; a^2 is set by dfpt_kerker_a2."; - item.default_value = "rho_plain"; + item.default_value = "v_sc_broyden"; item.read_value = [](const Input_Item& item, Parameter& para) { if (item.get_size() == 0) { diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 37589f894fd..4a00e89b803 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -143,6 +143,9 @@ 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_mix_type(const std::string& mix_type) { if (mix_type == "rho_plain" || mix_type == "rho_kerker" || mix_type == "rho_broyden" From 1a2f421203fcbe784be5cf3c5a2eecfc10d72732 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sun, 25 Oct 2026 12:44:47 +0800 Subject: [PATCH 10/21] fix(dfpt): address review findings on the mixing PR - migrate the remaining old mix_type tokens: the rho_kerker branch in mix_drho, the serial-test fixtures, the Impl default (now the valid v_sc_broyden), and the read_inp_dfpt empty-input fallback (was the rejected "plain"); CI MODULE_DFPT_rho_serial died on the fixture token - reset the shared Broyden mixer on every reset_mixing call, not only when the q already has a history entry: entering a fresh q left the mixer address-bound to the previous q's Mixing_Data and the first push died in the one-object-one-data guard - reduce the Broyden cal_coef inner products and the residual norms over the pool (Parallel_Reduce::reduce_pool behind __MPI guards, xc_pot.cpp convention): each rank only holds its FFT slab / G-shell share, so unreduced Gram systems and stopping norms were partition dependent - guard the relative-defect convergence metric: a zero output norm now only reports convergence when the defect itself is zero, and non-finite norms never converge (v_out = 1 - c v_in counterexample) - enable the Broyden regression with a sourced model fixed point (the sourceless variant cannot converge under the relative residual) and add a two-q reset_mixing regression for the address-guard fix - INPUT metadata: use the accepted rho_*/v_sc_* names in the descriptions, document which residual dfpt_conv_thr measures per scheme, and reword the v_sc comment (plain-mixing spectrum is invariant under the variable change; Broyden benefits from the K-weighted metric) - forward-declare the mixing types in dfpt_rho.h and drop the concrete includes (dfpt_pw_impl.h already includes them directly) --- .../module_parameter/read_inp_dfpt.cpp | 55 +++++--- source/source_pw/module_dfpt/dfpt_pw.h | 8 +- source/source_pw/module_dfpt/dfpt_pw_impl.h | 6 +- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 2 +- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 58 ++++++-- source/source_pw/module_dfpt/dfpt_rho.cpp | 72 +++++++--- source/source_pw/module_dfpt/dfpt_rho.h | 18 ++- .../test_serial/dfpt_rho_serial_test.cpp | 133 ++++++++++++++---- 8 files changed, 253 insertions(+), 99 deletions(-) diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index 06f69cae310..1ba227d7b2c 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -113,13 +113,17 @@ void ReadInput::item_dfpt() } { Input_Item item("dfpt_conv_thr"); - item.annotation = "convergence threshold of the DFPT first-order density"; + item.annotation = "convergence threshold of the DFPT response SCF"; item.category = category; item.type = "Real"; item.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."; + "the iteration stops when the relative residual of the mixed " + "variable drops below this value for every displacement. With " + "dfpt_mix_type = v_sc_plain or v_sc_broyden (the default) the " + "residual is the screening-potential residual ||v_out - v_in|| / " + "||v_out||; with rho_plain, rho_kerker or rho_broyden it is the " + "first-order density residual ||drho_out - drho_in|| / ||drho_out||. " + "The electric-field (LO-TO) solve always mixes the density."; item.default_value = "1.0e-8"; read_sync_double(input.dfpt_conv_thr); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -147,12 +151,14 @@ void ReadInput::item_dfpt() item.annotation = "mixing coefficient of the DFPT first-order density"; item.category = category; item.type = "Real"; - item.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."; + item.description = "Set the plain-mixing coefficient of the DFPT response SCF. It " + "scales the preliminary step of the chosen dfpt_mix_type (the " + "plain mix itself, or the initial guess pushed into the Broyden " + "history). 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."; item.default_value = "0.4"; read_sync_double(input.dfpt_mix_beta); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -204,17 +210,24 @@ void ReadInput::item_dfpt() item.annotation = "mixing scheme of the DFPT response SCF: rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden"; item.category = category; item.type = "String"; - item.description = "Set the mixing scheme of the first-order density in the " - "self-consistent DFPT cycle. plain: the same coefficient " - "dfpt_mix_beta on every shell. kerker: the Kerker screen " + item.description = "Set the mixing scheme of the DFPT response SCF. rho_plain: the " + "first-order density is mixed with the same coefficient " + "dfpt_mix_beta on every shell. rho_kerker: the Kerker screen " "f_g = |G+q|^2 / (|G+q|^2 + a^2) damps the Coulomb-stiff " "smallest shells, stabilizing polar / small-gap systems at " - "beta up to 1; a^2 is set by dfpt_kerker_a2."; + "beta up to 1; a^2 is set by dfpt_kerker_a2. rho_broyden: " + "Modified Broyden mixing (Johnson 1988) on the density, history " + "depth dfpt_broyden_ndim. v_sc_plain / v_sc_broyden: mix the " + "screening potential v_sc = v_H[drho] + v_xc[drho] instead of " + "the density, with plain or Modified Broyden coefficients; " + "v_sc_broyden is the default (the metric the Broyden least " + "squares works in is K-weighted, which stabilizes the " + "Coulomb-stiff shells that destabilize plain density mixing)."; item.default_value = "v_sc_broyden"; item.read_value = [](const Input_Item& item, Parameter& para) { if (item.get_size() == 0) { - para.input.dfpt_mix_type = "plain"; + para.input.dfpt_mix_type = "v_sc_broyden"; } else { @@ -235,13 +248,13 @@ void ReadInput::item_dfpt() } { Input_Item item("dfpt_kerker_a2"); - item.annotation = "Kerker screening parameter a^2 in 1/lat0^2 (dfpt_mix_type kerker only)"; + item.annotation = "Kerker screening parameter a^2 in 1/lat0^2 (dfpt_mix_type rho_kerker only)"; item.category = category; item.type = "Real"; item.description = "Set the Kerker screening parameter a^2 (in 1/lat0^2 units) used " - "by dfpt_mix_type = kerker. The screen f_g = |G+q|^2 / " + "by dfpt_mix_type = rho_kerker. The screen f_g = |G+q|^2 / " "(|G+q|^2 + a^2) freezes the stiff long-wavelength shells; " - "larger a^2 damps more shells. Ignored by the plain scheme."; + "larger a^2 damps more shells. Ignored by the other schemes."; item.default_value = "1.0"; read_sync_double(input.dfpt_kerker_a2); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -254,13 +267,13 @@ void ReadInput::item_dfpt() } { Input_Item item("dfpt_broyden_ndim"); - item.annotation = "Broyden history depth for dfpt_mix_type = broyden"; + item.annotation = "Broyden history depth for dfpt_mix_type = rho_broyden or v_sc_broyden"; item.category = category; item.type = "Integer"; item.description = "Set the number of history steps used by the " "Modified Broyden mixing (Johnson 1988) when " - "dfpt_mix_type = broyden. More history generally " - "accelerates convergence at the cost of memory."; + "dfpt_mix_type = rho_broyden or v_sc_broyden. More history " + "generally accelerates convergence at the cost of memory."; item.default_value = "8"; read_sync_int(input.dfpt_broyden_ndim); item.check_value = [](const Input_Item& item, const Parameter& para) { diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 0fbf8fe93b0..69ec3892d10 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -131,11 +131,13 @@ class DFPT_PW /// 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); - /// first-order density mixing scheme: "plain" (beta only) or "kerker" - /// (screen f_g = |G+q|^2/(|G+q|^2+a^2), a^2 via set_kerker_a2) + /// response SCF mixing scheme: "rho_plain" (beta only), "rho_kerker" + /// (screen f_g = |G+q|^2/(|G+q|^2+a^2), a^2 via set_kerker_a2), + /// "rho_broyden" (Modified Broyden on the density), or the potential + /// mixes "v_sc_plain" / "v_sc_broyden" (handled in solve_displacement) void set_mix_type(const std::string& mix_type); - /// Kerker screening a^2 in 1/lat0^2 units (kerker scheme only) + /// Kerker screening a^2 in 1/lat0^2 units (rho_kerker scheme only) void set_kerker_a2(double a2); /// Broyden history depth (broyden scheme only) diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index ab4e0bad134..f7ec543731d 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -180,9 +180,9 @@ class DFPT_PW::Impl int max_iter_ = 100; double mix_beta_ = 0.4; std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix - std::string mix_type_ = "plain"; ///< first-order density mixing scheme - double kerker_a2_ = 1.0; ///< Kerker screening a^2, 1/lat0^2 (kerker only) - int broyden_ndim_ = 8; ///< Broyden history depth (broyden only) + std::string mix_type_ = "v_sc_broyden"; ///< response SCF mixing scheme + double kerker_a2_ = 1.0; ///< Kerker screening a^2, 1/lat0^2 (rho_kerker only) + int broyden_ndim_ = 8; ///< Broyden history depth (rho/v_sc Broyden) ///< Broyden on the screening potential for v_sc mode (separate instance /// from the density Broyden because the array length differs: nrxx vs npw) std::unique_ptr vsc_broyden_; diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index d32f4abd7a3..3939312b710 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -87,7 +87,7 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, // Coulomb stiffness 4pi/G^2; measured lambda ~ -2.2 on {111}/{200} // for the diamond smoke case), so the coefficient must stay below // 2 / (1 + |lambda_min|); the INPUT default 0.4 keeps margin up to - // |lambda| ~ 3; the alternative is dfpt_mix_type = "kerker", the + // |lambda| ~ 3; the alternative is dfpt_mix_type = "rho_kerker", the // screen f_g = |G+q|^2 / (|G+q|^2 + a^2) in 1/lat0^2 units (a^2 via // dfpt_kerker_a2), which stabilizes those shells at beta up to 1; // both arrive through the INPUT setters (dfpt_mix_beta / diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index afeaa99bfd7..1c5a72feef8 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -11,6 +11,7 @@ #include "dfpt_pw_impl.h" #include "dfpt_pw_data.h" +#include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_base/tool_title.h" @@ -19,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -26,6 +28,26 @@ namespace ModuleDFPT { +namespace +{ +/// relative defect ||out - in|| / ||out|| with convergence-safe guards +/// (shared definition with dfpt_rho.cpp; see the note there): a zero +/// output norm only reports convergence when the defect itself is zero, +/// and non-finite norms never report convergence +double safe_relative_defect_(double dn2, double o2) +{ + if (!std::isfinite(dn2) || !std::isfinite(o2)) + { + return std::numeric_limits::max(); + } + if (o2 <= 0.0) + { + return (dn2 <= 0.0) ? 0.0 : std::numeric_limits::max(); + } + return std::sqrt(dn2 / o2); +} +} // namespace + void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, const std::vector>& drho_in_g, std::vector>& v_sc_r) const @@ -232,10 +254,14 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) int n_iter_done = 0; std::vector> v_sc_r_last; // v_sc mode: the SCF mixes the screening potential instead of the - // density -- the fixed-point variable change provides implicit K^-1 - // preconditioning (the Hartree+XC kernel's inverse), flattening the - // Coulomb-stiff smallest-G shells that destabilize plain density - // mixing (see MIXING_20261020.md §8) + // density. For a fixed linear kernel the density and potential + // fixed-point maps have similar Jacobians (chi*K vs K*chi on the + // responsive subspace), so scalar plain mixing keeps the same + // spectrum: the variable change alone is not an inverse-kernel + // preconditioner. What changes is the metric the Broyden least + // squares works in (the K-weighted inner product), which is why + // v_sc_broyden is the stabilized default while v_sc_plain tracks + // the plain density mixing (see MIXING_20261020.md §8) const bool vsc_mode = (mix_type_ == "v_sc_plain" || mix_type_ == "v_sc_broyden"); std::vector> v_sc_in_r(nrxx, std::complex(0.0, 0.0)); // v_sc Broyden: a separate instance from the density Broyden because @@ -340,9 +366,16 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) { s += (a[i] * std::conj(b[i])).real(); } +#ifdef __MPI + // pool-global inner product: each rank sums + // its FFT slab share only, and the Gram system + // must be identical on every rank + Parallel_Reduce::reduce_pool(s); +#endif return s; })); - // residual from the unmixed in/out difference + // residual from the unmixed in/out difference; the norms + // are pool-global quantities (rank-local FFT slabs) double dn2 = 0.0; double o2 = 0.0; for (int ir = 0; ir < nrxx; ++ir) @@ -350,12 +383,17 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) dn2 += std::norm(v_sc_out_r[ir] - v_sc_in_r[ir]); o2 += std::norm(v_sc_out_r[ir]); } - residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; +#ifdef __MPI + Parallel_Reduce::reduce_pool(dn2); + Parallel_Reduce::reduce_pool(o2); +#endif + residual = safe_relative_defect_(dn2, o2); vsc_broyden_->mix_data(*vsc_broyden_mdata_, v_sc_in_r.data()); } else { - // plain beta mix on the real-space potential + // plain beta mix on the real-space potential; the norms are + // pool-global quantities (rank-local FFT slabs) double dn2 = 0.0; double o2 = 0.0; for (int ir = 0; ir < nrxx; ++ir) @@ -365,7 +403,11 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) o2 += std::norm(v_sc_out_r[ir]); v_sc_in_r[ir] += mix_beta_ * delta; } - residual = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; +#ifdef __MPI + Parallel_Reduce::reduce_pool(dn2); + Parallel_Reduce::reduce_pool(o2); +#endif + residual = safe_relative_defect_(dn2, o2); } // store the mixed density for the 2n+1 accumulation path: // rebuild the G-space coefficients from the Sternheimer output diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 3a2f0a4ad1a..340090fd255 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -6,6 +6,8 @@ #include "source_base/timer.h" #include "source_base/tool_title.h" #include "source_base/module_mixing/plain_mixing.h" +#include "source_base/module_mixing/broyden_mixing.h" +#include "source_base/parallel_reduce.h" #include "source_base/math_integral.h" #include "source_base/math_sphbes.h" #include "source_estate/module_charge/chg_tools.h" @@ -14,17 +16,38 @@ #include #include #include +#include #include namespace ModuleDFPT { +namespace +{ +/// relative defect ||out - in|| / ||out|| with convergence-safe guards: +/// a zero output norm only reports convergence when the defect itself is +/// zero (v_out = 1 - c v_in maps a nonzero input to zero output far from +/// its fixed point), and non-finite norms never report convergence +double safe_relative_defect_(double dn2, double o2) +{ + if (!std::isfinite(dn2) || !std::isfinite(o2)) + { + return std::numeric_limits::max(); + } + if (o2 <= 0.0) + { + return (dn2 <= 0.0) ? 0.0 : std::numeric_limits::max(); + } + return std::sqrt(dn2 / o2); +} +} // namespace + DFPT_Rho::DFPT_Rho() { } -// Defined here (after plain_mixing.h is included) so the unique_ptr deleter -// sees the complete Plain_Mixing type. The body is implicitly generated. +// Defined here (after the mixing headers are included) so the unique_ptr +// deleters see the complete types. The body is implicitly generated. DFPT_Rho::~DFPT_Rho() = default; void DFPT_Rho::init(const Config& cfg) @@ -397,15 +420,19 @@ void DFPT_Rho::reset_mixing(int q_idx) residual_[q_idx] = 0.0; } // Broyden history is per (q, displacement-cycle): drop it so the next - // displacement starts fresh; also reset the Broyden mixer's internal - // address binding (it must not point to the destroyed Mixing_Data) + // displacement starts fresh if (q_idx < static_cast(broyden_mdata_.size())) { broyden_mdata_[q_idx].reset(); - if (broyden_mixer_ != nullptr) - { - broyden_mixer_->reset(); - } + } + // the shared Broyden mixer must be reset on every valid cycle, + // independently of whether this q already has a history entry: when a + // new q is entered for the first time the mixer is still address-bound + // to the previous q's Mixing_Data, and its first push would exit + // through the one-object-one-data address guard (broyden_mixing.cpp) + if (broyden_mixer_ != nullptr) + { + broyden_mixer_->reset(); } ModuleBase::timer::end("DFPT_Rho", "reset_mixing"); } @@ -481,17 +508,10 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) screen, combine, true); - auto inner = std::function*, std::complex*)>( - [](std::complex* a, std::complex* b) { - double s = 0.0; - for (int i = 0; i < 1; ++i) - { - s += (a[i] * std::conj(b[i])).real(); - } - return s; - }); - // the inner product needs the full length; use the templated path - // through the mixing data (the Broyden class handles the wrapping) + // the inner product is a pool-global quantity: each rank only + // holds its G-shell share, so the partial sum must be reduced + // before the Gram system is assembled (otherwise every rank + // solves a different least-squares problem) broyden_mixer_->cal_coef(mdata, std::function*, std::complex*)>( [npw](std::complex* a, std::complex* b) { @@ -500,11 +520,14 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) { s += (a[i] * std::conj(b[i])).real(); } +#ifdef __MPI + Parallel_Reduce::reduce_pool(s); +#endif return s; })); broyden_mixer_->mix_data(mdata, mixed.data()); } - else if (mix_type_ == "kerker") + else if (mix_type_ == "rho_kerker") { const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; std::vector> rin_s(npw); @@ -537,7 +560,8 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) { mixer_->plain_mix(mixed.data(), rin.data(), out.data(), npw, std::function*)>()); } - // relative residual ||out - in|| / ||out|| + // relative residual ||out - in|| / ||out||; the norms are pool-global + // quantities (each rank sums its G shells only) double dn2 = 0.0; double o2 = 0.0; for (int ig = 0; ig < npw; ++ig) @@ -545,7 +569,11 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) dn2 += std::norm(out[ig] - rin[ig]); o2 += std::norm(out[ig]); } - residual_[q_idx] = (o2 > 0.0) ? std::sqrt(dn2 / o2) : 0.0; +#ifdef __MPI + Parallel_Reduce::reduce_pool(dn2); + Parallel_Reduce::reduce_pool(o2); +#endif + residual_[q_idx] = safe_relative_defect_(dn2, o2); drho_in_[q_idx][0] = mixed; data.set_drho_g(q_idx, 0, mixed); diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 7d726373de5..97dd4c6f429 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -8,10 +8,6 @@ #include "source_cell/unitcell.h" #include "source_psi/psi.h" -#include "source_base/module_mixing/mixing.h" -#include "source_base/module_mixing/plain_mixing.h" -#include "source_base/module_mixing/broyden_mixing.h" - #include #include #include @@ -19,6 +15,8 @@ namespace Base_Mixing { class Plain_Mixing; +class Broyden_Mixing; +class Mixing_Data; } namespace ModuleDFPT @@ -80,8 +78,8 @@ class XC_First_Order * mix_drho mixes the q-shifted coefficients: * drho_in <- drho_in + beta_g (drho_out - drho_in) * through Base_Mixing::Plain_Mixing (no Charge_Mixing / Charge dependency). - * With mix_type = "plain" beta_g = beta on every shell; with mix_type = - * "kerker" beta_g = beta * f_g, the Kerker screen + * With mix_type = "rho_plain" beta_g = beta on every shell; with mix_type = + * "rho_kerker" beta_g = beta * f_g, the Kerker screen * f_g = |G+q|^2 / (|G+q|^2 + a^2), * evaluated with the same q_shifted |G+q| = |gcar + q_cart * recip| (in * 1/lat0^2 units) convention as v_hartree_q, so the Coulomb-stiffness @@ -104,10 +102,10 @@ class DFPT_Rho ModulePW::PW_Basis* pw_rho; ModulePW::PW_Basis_K* pw_wfc; ModuleBase::Matrix3 recip_matrix; - std::string mix_type; // "plain", "kerker", "broyden", or "v_sc" + std::string mix_type; // "rho_plain", "rho_kerker", "rho_broyden", "v_sc_plain", "v_sc_broyden" double mix_beta; - double kerker_a2; // Kerker screen a^2, 1/lat0^2 (kerker only) - int broyden_ndim; // Broyden history depth (broyden only, default 8) + double kerker_a2; // Kerker screen a^2, 1/lat0^2 (rho_kerker only) + int broyden_ndim; // Broyden history depth (rho_broyden only, default 8) }; DFPT_Rho(); @@ -191,7 +189,7 @@ class DFPT_Rho ///< reciprocal lattice matrix in 1/lat0 (UnitCell::G convention) ModuleBase::Matrix3 recip_matrix_; double mix_beta_ = 0.7; - ///< mixing algorithm: "plain" (beta only) or "kerker" (Kerker screen) + ///< mixing algorithm: rho_plain | rho_kerker | rho_broyden (v_sc_* at solver) std::string mix_type_; ///< Kerker screening parameter a^2 in 1/lat0^2 (same units as |G+q|^2) double kerker_a2_ = 0.0; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index b43caaf3f82..f5f09989e5b 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -65,7 +65,7 @@ class DFPTRhoSerialTest : public DFPTSerialBase void SetUp() override { DFPTSerialBase::SetUp(); - rho_.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("plain"), 0.4, 0.0}); + rho_.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_plain"), 0.4, 0.0}); } void FillRandomStates(psi::Psi>& psi, @@ -242,7 +242,7 @@ TEST_F(DFPTRhoSerialTest, ChargeConservationAtGamma) ModuleDFPT::DFPT_PW_Data data0; data0.init(&qlist0, 1, nbands_, pw_wfc0.npwk_max, pw_rho_.nrxx, 1, 1, nullptr); ModuleDFPT::DFPT_Rho rho0; - rho0.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc0, G_, std::string("plain"), 0.4, 0.0}); + rho0.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc0, G_, std::string("rho_plain"), 0.4, 0.0}); psi::Psi> psi(1, nbands_, pw_wfc0.npwk_max, pw_wfc0.npwk[0], true); ModuleDFPT::DFPT_KQ_Basis kq0; @@ -410,7 +410,7 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerFirstStepIsPreconditionedScaledOutput) } ASSERT_GT(w2_min, 0.0); const double a2 = 4.0 * w2_min; - rho_k.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("kerker"), 0.7, a2}); + rho_k.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_kerker"), 0.7, a2}); std::vector> out(pw_rho_.npw); int n_small = 0; @@ -500,7 +500,7 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerStabilizesStiffModelProblem) // plain beta = 0.7 on the stiff model diverges ModuleDFPT::DFPT_Rho rho_p; - rho_p.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("plain"), 0.7, 0.0}); + rho_p.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_plain"), 0.7, 0.0}); data_.set_drho_g(0, 0, std::vector>(npw, std::complex(0.0, 0.0))); for (int it = 0; it < 40; ++it) { @@ -513,7 +513,7 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerStabilizesStiffModelProblem) // kerker beta = 0.7 with a^2 = 9 w2_min (f ~ 0.1 on the stiff shell) // converges to the target ModuleDFPT::DFPT_Rho rho_k; - rho_k.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("kerker"), 0.7, 9.0 * w2_min}); + rho_k.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_kerker"), 0.7, 9.0 * w2_min}); data_.set_drho_g(0, 0, std::vector>(npw, std::complex(0.0, 0.0))); for (int it = 0; it < 300; ++it) { @@ -530,12 +530,15 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerStabilizesStiffModelProblem) } } -#if 0 // TODO: fix Broyden mock interface; real-system G5 validation covers this TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) { - // synthetic fixed-point: drho_out = A * drho_in, with A a simple - // contraction; Broyden should converge in few iterations (the plain - // mixer needs O(1/(1-|A|)) steps) + // synthetic fixed point: out = D in + (1 - D) t with a diagonal D + // (0.5 and 0.3 on two G shells, 0 elsewhere) and a fixed target t, + // so the response map has a source and the relative residual + // ||out - in|| / ||out|| genuinely decays. Plain beta = 0.4 + // contracts every shell by |1 - beta (1 - D)| >= 0.6 (>= 90 steps + // to 1e-10); Modified Broyden on the few distinct eigenvalues + // should converge in a handful of iterations. SetupBases(k_d_, q_d_, 2); ModuleDFPT::DFPT_Rho::Config cfg; cfg.nspin = 1; @@ -543,28 +546,34 @@ TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) cfg.pw_rho = &pw_rho_; cfg.pw_wfc = &pw_wfc_; cfg.recip_matrix = G_; - cfg.mix_type = "broyden"; + cfg.mix_type = "rho_broyden"; cfg.mix_beta = 0.4; cfg.kerker_a2 = 1.0; cfg.broyden_ndim = 8; ModuleDFPT::DFPT_Rho rho; rho.init(cfg); - // use a small 2-component toy problem mapped onto npw (only two - // nonzero G components with different eigenvalues of the response) const int npw = pw_rho_.npw; - // target: the fixed point of drho_out = A drho_in - // choose A = diag(0.5, 0.3, 0, ...): the fixed point is 0 - // start from a nonzero initial guess via the first output - std::vector> rho_out(npw, std::complex(0.0, 0.0)); - rho_out[0] = std::complex(1.0, 0.0); - rho_out[1] = std::complex(0.5, 0.2); - - // run mixing iterations: each step, output = A * input - // (the Broyden should identify the linear map and jump to the - // solution after collecting enough history) + std::vector> target(npw, std::complex(0.0, 0.0)); + target[0] = std::complex(0.9, -0.2); + target[1] = std::complex(0.5, 0.2); + std::vector d(npw, 0.0); + d[0] = 0.5; + d[1] = 0.3; + auto model_out = [&](const std::vector>& in) { + std::vector> o(npw); + for (int ig = 0; ig < npw; ++ig) + { + o[ig] = d[ig] * in[ig] + (1.0 - d[ig]) * target[ig]; + } + return o; + }; + + // run mixing iterations from the zero input; each step the model + // rebuilds the output from the mixed input int n_iter = 0; - for (int iter = 0; iter < 20; ++iter) + std::vector> rho_out = model_out(std::vector>(npw, std::complex(0.0, 0.0))); + for (int iter = 0; iter < 40; ++iter) { data_.set_drho_g(0, 0, rho_out); rho.mix_drho(0, data_); @@ -574,16 +583,78 @@ TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) { break; } - // next output = A * current input (the linear model) - const std::vector> mixed = data_.get_drho_g(0, 0); - rho_out.assign(npw, std::complex(0.0, 0.0)); - rho_out[0] = 0.5 * mixed[0]; - rho_out[1] = 0.3 * mixed[1]; + rho_out = model_out(data_.get_drho_g(0, 0)); } - // Broyden on a 2-eigenvalue diagonal system: should converge - // dramatically faster than plain (which needs ~20 steps) + // Broyden on a 3-eigenvalue diagonal system: converges dramatically + // faster than plain mixing (which needs >= 90 steps at 0.6/step) EXPECT_LT(n_iter, 12); const double final_res = rho.get_residual(0, data_); EXPECT_LT(final_res, 1.0e-8); + // the converged input is the fixed point of the model map + const std::vector> final_in = data_.get_drho_g(0, 0); + for (int ig = 0; ig < npw; ++ig) + { + EXPECT_NEAR(final_in[ig].real(), target[ig].real(), 1.0e-8); + EXPECT_NEAR(final_in[ig].imag(), target[ig].imag(), 1.0e-8); + } +} + +TEST_F(DFPTRhoSerialTest, BroydenResetMixingAllowsFreshQAndSameQCycle) +{ + // regression (PR review): entering a fresh q after q = 0 must reset + // the shared Broyden mixer's address binding even though the new q + // has no history entry yet -- otherwise the first push of the newly + // allocated Mixing_Data exits through the one-object-one-data guard + // (broyden_mixing.cpp) and kills the process. The same-q cycle must + // also restart from a zero input density. + SetupBases(k_d_, q_d_, 2); + ModuleDFPT::DFPT_Rho::Config cfg; + cfg.nspin = 1; + cfg.nrxx = pw_rho_.nrxx; + cfg.pw_rho = &pw_rho_; + cfg.pw_wfc = &pw_wfc_; + cfg.recip_matrix = G_; + cfg.mix_type = "rho_broyden"; + cfg.mix_beta = 0.4; + cfg.kerker_a2 = 1.0; + cfg.broyden_ndim = 8; + ModuleDFPT::DFPT_Rho rho; + rho.init(cfg); + + const int npw = pw_rho_.npw; + std::vector> out(npw); + for (int ig = 0; ig < npw; ++ig) + { + out[ig] = 0.01 * std::complex(std::cos(0.3 * ig), std::sin(0.9 * ig)); + } + + // two pushes at q = 0: the mixer binds the q = 0 Mixing_Data + data_.set_drho_g(0, 0, out); + rho.mix_drho(0, data_); + data_.set_drho_g(0, 0, out); + rho.mix_drho(0, data_); + + // entering q = 1 for the first time (no history entry exists): with + // the buggy guard the mixer stayed bound to q = 0 and this died + rho.reset_mixing(1); + data_.set_drho_g(1, 0, out); + rho.mix_drho(1, data_); + // fresh q starts from the zero input: first step is beta * out + const std::vector> mixed_q1 = data_.get_drho_g(1, 0); + for (int ig = 0; ig < npw; ++ig) + { + EXPECT_NEAR(mixed_q1[ig].real(), 0.4 * out[ig].real(), 1.0e-12); + EXPECT_NEAR(mixed_q1[ig].imag(), 0.4 * out[ig].imag(), 1.0e-12); + } + + // the same-q cycle also restarts from a zero input density + rho.reset_mixing(0); + data_.set_drho_g(0, 0, out); + rho.mix_drho(0, data_); + const std::vector> mixed_q0 = data_.get_drho_g(0, 0); + for (int ig = 0; ig < npw; ++ig) + { + EXPECT_NEAR(mixed_q0[ig].real(), 0.4 * out[ig].real(), 1.0e-12); + EXPECT_NEAR(mixed_q0[ig].imag(), 0.4 * out[ig].imag(), 1.0e-12); + } } -#endif From 07276d1065bf88395b2e06e52e9a69b19a1af94d Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 26 Oct 2026 13:46:28 +0800 Subject: [PATCH 11/21] test(dfpt): correct the plain-mixing contraction note in the Broyden regression Review follow-up: the binding populated shell is d = 0.5 with factor |1 - beta (1 - D)| = 0.8 (>= 103 plain steps to 1e-10), not the unpopulated-shell bound 0.6. Assertions are unchanged. --- .../module_dfpt/test_serial/dfpt_rho_serial_test.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index f5f09989e5b..11c137c6c7a 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -536,9 +536,10 @@ TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) // (0.5 and 0.3 on two G shells, 0 elsewhere) and a fixed target t, // so the response map has a source and the relative residual // ||out - in|| / ||out|| genuinely decays. Plain beta = 0.4 - // contracts every shell by |1 - beta (1 - D)| >= 0.6 (>= 90 steps - // to 1e-10); Modified Broyden on the few distinct eigenvalues - // should converge in a handful of iterations. + // contracts every populated shell by |1 - beta (1 - D)| = 0.8 (d = 0.5, + // the binding shell; 0.72 for d = 0.3), i.e. >= 103 plain steps to + // 1e-10; Modified Broyden on the few distinct eigenvalues should + // converge in a handful of iterations. SetupBases(k_d_, q_d_, 2); ModuleDFPT::DFPT_Rho::Config cfg; cfg.nspin = 1; From 17bfa50be58b79ffbf32fb9f181596dd82a5477d Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 26 Oct 2026 19:56:53 +0800 Subject: [PATCH 12/21] fix(dfpt): address the human-review nits on the E-field loop and init - E-field solve (dfpt_pw_q0.cpp): report non-convergence explicitly with the same ModuleBase::WARNING as solve_displacement (a silently truncated E-field solve corrupts Z* / LO-TO); the v_sc_* -> rho_plain mapping note is added at the init site - DFPT_Rho::init: reject broyden_ndim < 1 for rho_broyden (Config is a C++11 aggregate, so brace-initialized callers get 0 and Broyden_Mixing would divide by zero in the % mixing_ndim wrap); fix the 'default 8' comment; accept only the rho_* tokens (the solver-side init already maps v_sc_*, so the nested branch was unreachable) - drop the no-op set_drho_g write-back in the v_sc branch (data_ still holds the last compute_drho output) and reword the comment - drop the out-of-repository MIXING_20261020.md reference - rename safe_relative_defect_ -> safe_relative_defect (free function in an anonymous namespace; the trailing underscore is the member convention); sync the plain-mixing contraction notes in the test (roughly 100 iterations at 0.8/step on the binding shell) --- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 3 ++- source/source_pw/module_dfpt/dfpt_pw_q0.cpp | 18 ++++++++++++++++ .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 18 ++++++++-------- source/source_pw/module_dfpt/dfpt_rho.cpp | 21 +++++++++++-------- source/source_pw/module_dfpt/dfpt_rho.h | 4 ++-- .../test_serial/dfpt_rho_serial_test.cpp | 9 ++++---- 6 files changed, 48 insertions(+), 25 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 3939312b710..8c8f5a0d8b6 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -95,8 +95,9 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, const double mix_beta = mix_beta_; const std::string mix_type = mix_type_; const double kerker_a2 = kerker_a2_; - // DFPT_Rho handles density mixing only; v_sc modes do their own + // DFPT_Rho handles density mixing only; v_sc modes do their own // potential mixing in solve_displacement, so pass the plain default + // (the E-field loop in dfpt_pw_q0.cpp always mixes the density) const std::string rho_mix_type = (mix_type == "v_sc_plain" || mix_type == "v_sc_broyden") ? "rho_plain" : mix_type; rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, rho_mix_type, mix_beta, kerker_a2, broyden_ndim_}); diff --git a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp index 39dba1d3f8d..16758a7231f 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp @@ -13,12 +13,15 @@ #include "dfpt_pw_data.h" #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_base/tool_title.h" #include #include #include +#include #include +#include #include namespace ModuleDFPT @@ -312,6 +315,8 @@ void DFPT_PW::Impl::solve_efield_resp(int q_idx) rho_.reset_mixing(q_idx); data_.set_drho_g(q_idx, 0, std::vector>(npw, std::complex(0.0, 0.0))); bool converged = false; + int n_iter_done = 0; + double residual_last = 0.0; for (int iter = 0; iter < max_iter_ && !converged; ++iter) { // screened response potential of the mixed input density @@ -338,13 +343,26 @@ void DFPT_PW::Impl::solve_efield_resp(int q_idx) rho_.compute_drho(gs_psi_, wg_, q_idx, data_); rho_.mix_drho(q_idx, data_); const double residual = rho_.get_residual(q_idx, data_); + residual_last = residual; converged = (residual < conv_thr_); + ++n_iter_done; if (converged) { std::cout << "DFPT efield dir=" << a << " converged, residual=" << residual << " (iter=" << iter << ")" << std::endl; } } + if (!converged) + { + // same explicit non-convergence report as solve_displacement: + // a silently truncated E-field solve corrupts Z* / LO-TO + std::ostringstream oss; + oss << "the E-field response (dir " << static_cast('x' + a) + << ") did not converge within dfpt_max_iter=" << max_iter_ << ": final residual " << std::scientific + << std::setprecision(3) << residual_last << " vs threshold " << std::setprecision(1) << conv_thr_ + << std::defaultfloat << ". Z* / LO-TO splitting from this solve is unreliable"; + ModuleBase::WARNING("DFPT_PW::solve_efield_resp", oss.str()); + } stash_dpsi_efield(q_idx, a, nk, nbands); } ModuleBase::timer::end("DFPT_PW", "solve_efield_resp"); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 1c5a72feef8..876519826f7 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -34,7 +34,7 @@ namespace /// (shared definition with dfpt_rho.cpp; see the note there): a zero /// output norm only reports convergence when the defect itself is zero, /// and non-finite norms never report convergence -double safe_relative_defect_(double dn2, double o2) +double safe_relative_defect(double dn2, double o2) { if (!std::isfinite(dn2) || !std::isfinite(o2)) { @@ -261,7 +261,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) // preconditioner. What changes is the metric the Broyden least // squares works in (the K-weighted inner product), which is why // v_sc_broyden is the stabilized default while v_sc_plain tracks - // the plain density mixing (see MIXING_20261020.md §8) + // the plain density mixing const bool vsc_mode = (mix_type_ == "v_sc_plain" || mix_type_ == "v_sc_broyden"); std::vector> v_sc_in_r(nrxx, std::complex(0.0, 0.0)); // v_sc Broyden: a separate instance from the density Broyden because @@ -387,7 +387,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) Parallel_Reduce::reduce_pool(dn2); Parallel_Reduce::reduce_pool(o2); #endif - residual = safe_relative_defect_(dn2, o2); + residual = safe_relative_defect(dn2, o2); vsc_broyden_->mix_data(*vsc_broyden_mdata_, v_sc_in_r.data()); } else @@ -407,13 +407,13 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) Parallel_Reduce::reduce_pool(dn2); Parallel_Reduce::reduce_pool(o2); #endif - residual = safe_relative_defect_(dn2, o2); + residual = safe_relative_defect(dn2, o2); } - // store the mixed density for the 2n+1 accumulation path: - // rebuild the G-space coefficients from the Sternheimer output - // (the density is NOT what we mix, but the 2n+1 needs the - // converged density from the converged potential) - data_.set_drho_g(q_idx, 0, drho_out_g); + // the 2n+1 accumulation path reads the density from data_: + // unlike the density-mixing branch, the v_sc loop never calls + // mix_drho, so data_ still holds the last compute_drho output + // (drho_out_g was read from it), which at convergence is the + // density consistent with the final screening potential } else { diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 340090fd255..326c4c7f04f 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -28,7 +28,7 @@ namespace /// a zero output norm only reports convergence when the defect itself is /// zero (v_out = 1 - c v_in maps a nonzero input to zero output far from /// its fixed point), and non-finite norms never report convergence -double safe_relative_defect_(double dn2, double o2) +double safe_relative_defect(double dn2, double o2) { if (!std::isfinite(dn2) || !std::isfinite(o2)) { @@ -62,20 +62,23 @@ void DFPT_Rho::init(const Config& cfg) mix_beta_ = cfg.mix_beta; mix_type_ = cfg.mix_type; // internal: rho_plain|rho_kerker|rho_broyden (v_sc handled at solver) kerker_a2_ = cfg.kerker_a2; + // DFPT_Rho handles density mixing only; the solver-side init maps the + // v_sc_* INPUT modes onto the plain density mixer for its E-field loop if (cfg.mix_type != "rho_plain" && cfg.mix_type != "rho_kerker" && cfg.mix_type != "rho_broyden") { - // v_sc_* modes are handled in solve_displacement, not in DFPT_Rho - if (cfg.mix_type != "v_sc_plain" && cfg.mix_type != "v_sc_broyden") - { - ModuleBase::WARNING_QUIT( - "DFPT_Rho", - "unsupported mix_type, expected rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden"); - } + ModuleBase::WARNING_QUIT("DFPT_Rho", "unsupported mix_type, expected rho_plain, rho_kerker or rho_broyden"); } // make_unique is C++14; C++11 uses new directly through reset() mixer_.reset(new Base_Mixing::Plain_Mixing(mix_beta_)); if (mix_type_ == "rho_broyden") { + // Config is an aggregate initialized with braces, so a caller + // passing fewer fields gets broyden_ndim = 0; Broyden_Mixing + // would then divide by zero in its % mixing_ndim history wrap + if (cfg.broyden_ndim < 1) + { + ModuleBase::WARNING_QUIT("DFPT_Rho", "broyden_ndim must be >= 1 for rho_broyden mixing"); + } broyden_mixer_.reset(new Base_Mixing::Broyden_Mixing(cfg.broyden_ndim)); } ModuleBase::timer::end("DFPT_Rho", "init"); @@ -573,7 +576,7 @@ void DFPT_Rho::mix_drho(int q_idx, DFPT_PW_Data& data) Parallel_Reduce::reduce_pool(dn2); Parallel_Reduce::reduce_pool(o2); #endif - residual_[q_idx] = safe_relative_defect_(dn2, o2); + residual_[q_idx] = safe_relative_defect(dn2, o2); drho_in_[q_idx][0] = mixed; data.set_drho_g(q_idx, 0, mixed); diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 97dd4c6f429..89cb83a95a4 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -102,10 +102,10 @@ class DFPT_Rho ModulePW::PW_Basis* pw_rho; ModulePW::PW_Basis_K* pw_wfc; ModuleBase::Matrix3 recip_matrix; - std::string mix_type; // "rho_plain", "rho_kerker", "rho_broyden", "v_sc_plain", "v_sc_broyden" + std::string mix_type; // "rho_plain", "rho_kerker" or "rho_broyden" (v_sc_* is handled at the solver) double mix_beta; double kerker_a2; // Kerker screen a^2, 1/lat0^2 (rho_kerker only) - int broyden_ndim; // Broyden history depth (rho_broyden only, default 8) + int broyden_ndim; // Broyden history depth (rho_broyden only, must be >= 1; no C++11 aggregate default) }; DFPT_Rho(); diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 11c137c6c7a..66c0543f4a6 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -537,9 +537,9 @@ TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) // so the response map has a source and the relative residual // ||out - in|| / ||out|| genuinely decays. Plain beta = 0.4 // contracts every populated shell by |1 - beta (1 - D)| = 0.8 (d = 0.5, - // the binding shell; 0.72 for d = 0.3), i.e. >= 103 plain steps to - // 1e-10; Modified Broyden on the few distinct eigenvalues should - // converge in a handful of iterations. + // the binding shell; 0.72 for d = 0.3), i.e. roughly 100 plain + // iterations to 1e-10; Modified Broyden on the few distinct + // eigenvalues should converge in a handful of iterations. SetupBases(k_d_, q_d_, 2); ModuleDFPT::DFPT_Rho::Config cfg; cfg.nspin = 1; @@ -587,7 +587,8 @@ TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) rho_out = model_out(data_.get_drho_g(0, 0)); } // Broyden on a 3-eigenvalue diagonal system: converges dramatically - // faster than plain mixing (which needs >= 90 steps at 0.6/step) + // faster than plain mixing (which needs roughly 100 iterations at + // 0.8/step on the binding shell) EXPECT_LT(n_iter, 12); const double final_res = rho.get_residual(0, data_); EXPECT_LT(final_res, 1.0e-8); From 81c95eec24c53e759ee65200b235e27125af5006 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Mon, 26 Oct 2026 21:42:50 +0800 Subject: [PATCH 13/21] feat(dfpt): give the E-field solve the density-space Broyden mixer under v_sc_broyden Review follow-up (Critsium-xy): with the v_sc_broyden default, the E-field solve (Z* / epsilon_inf / LO-TO) silently fell back to plain density mixing at beta 0.4 -- the weakest member of the family exactly on the path that matters most for polar systems. Map v_sc_broyden to rho_broyden for DFPT_Rho (v_sc_plain keeps rho_plain), so the E-field loop runs Modified Broyden on the density while the displacement loop keeps its potential-space Broyden. Verification (si, Gamma + dfpt_compute_q0/dfpt_loto, default INPUT so the path is exercised through the mapping): E-field 3/3 directions converge in 12 iterations (bit-identical residuals to an explicit rho_broyden run) versus 38 with the previous plain fallback -- displacements unchanged at 13 iterations (v_sc Broyden); Gamma frequencies and the LO-TO block reproduce the pre-change values exactly. INPUT metadata documents that dfpt_broyden_ndim now sizes both Broyden instances and that the E-field residual is always the density residual. --- .../module_parameter/read_inp_dfpt.cpp | 13 +++++++++--- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 20 ++++++++++++++----- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index 1ba227d7b2c..a673bf2548c 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -123,7 +123,10 @@ void ReadInput::item_dfpt() "residual is the screening-potential residual ||v_out - v_in|| / " "||v_out||; with rho_plain, rho_kerker or rho_broyden it is the " "first-order density residual ||drho_out - drho_in|| / ||drho_out||. " - "The electric-field (LO-TO) solve always mixes the density."; + "The electric-field (LO-TO) solve always mixes the density, using " + "the density-space counterpart of the selected scheme (e.g. " + "rho_broyden when dfpt_mix_type = v_sc_broyden), so its residual " + "is always the density residual."; item.default_value = "1.0e-8"; read_sync_double(input.dfpt_conv_thr); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -272,8 +275,12 @@ void ReadInput::item_dfpt() item.type = "Integer"; item.description = "Set the number of history steps used by the " "Modified Broyden mixing (Johnson 1988) when " - "dfpt_mix_type = rho_broyden or v_sc_broyden. More history " - "generally accelerates convergence at the cost of memory."; + "dfpt_mix_type = rho_broyden or v_sc_broyden. Under " + "v_sc_broyden it sizes both Broyden instances: the " + "screening-potential mixer of the displacement loop " + "and the density mixer of the electric-field " + "(LO-TO) loop. More history generally accelerates " + "convergence at the cost of memory."; item.default_value = "8"; read_sync_int(input.dfpt_broyden_ndim); item.check_value = [](const Input_Item& item, const Parameter& para) { diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 8c8f5a0d8b6..47240b6618c 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -95,11 +95,21 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, const double mix_beta = mix_beta_; const std::string mix_type = mix_type_; const double kerker_a2 = kerker_a2_; - // DFPT_Rho handles density mixing only; v_sc modes do their own - // potential mixing in solve_displacement, so pass the plain default - // (the E-field loop in dfpt_pw_q0.cpp always mixes the density) - const std::string rho_mix_type - = (mix_type == "v_sc_plain" || mix_type == "v_sc_broyden") ? "rho_plain" : mix_type; + // DFPT_Rho handles density mixing only (the E-field loop in + // dfpt_pw_q0.cpp always mixes the density, while the v_sc modes do + // their own potential mixing in solve_displacement), so pass the + // density-space counterpart: v_sc_plain -> rho_plain, v_sc_broyden + // -> rho_broyden (the E-field solve must not silently fall back to + // the plain beta mix when the INPUT selects a Broyden scheme) + std::string rho_mix_type = mix_type; + if (mix_type == "v_sc_plain") + { + rho_mix_type = "rho_plain"; + } + else if (mix_type == "v_sc_broyden") + { + rho_mix_type = "rho_broyden"; + } rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, rho_mix_type, mix_beta, kerker_a2, broyden_ndim_}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); From 4d6600b29eba4d052448a2ff9581ec4fcb20ade1 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Tue, 27 Oct 2026 19:48:15 +0800 Subject: [PATCH 14/21] feat(dfpt): the two missing nlcc 2n+1 direct terms (QE dynmat_nlcc + dynmatcc) The nlcc response was previously wired only through the SCF screening channel (f_xc.(drho_val + drhoc) inside assemble_v_sc, QE dv_of_drho semantics, verified line by line). QE ph.x adds two direct terms to the force constants that have no wavefunction representation, and the bn_hex isolation experiment (identical pseudos with the PP_NLCC block zeroed: ABACUS vs QE optical rms 0.9 cm^-1 without nlcc, while the SCF-only nlcc overshoots QE's nlcc effect by 4-12x with wrong signs on the A1/E2 soft modes) shows they are the dominant cancellation partner of the screening contribution: - kernel term (QE dynmat_nlcc): D_ab += Omega/nrxx sum_r dvxc^a(r) conj(drhoc^b(r)) with dvxc^a = f_xc.(drho_val^a + drhoc^a) of the converged displacement a and the core response of every column b; dvxc^a is rebuilt in the two-pass accumulation loop from the newly stashed per-displacement density response plus a recomputed core response (the working slots hold the last solved displacement). - diagonal term (QE dynmatcc): D_aa[i][j] += Int dr w_ij(r) v_xc^0(r) with w_ij the second-order shifted-core amplitude -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau), evaluated through recip2real as a real field so the half-space G-sum convention never enters directly; v_xc^0 is the ground-state XC potential on the same valence+core density the kernel linearizes around, computed once in the esolver wiring (XC_Functional::v_xc) and delivered through the new DFPT_PW_InitContext::vxc_r provider (null disables the term). Also: DFPT_Rho::compute_drhoc_g extracted as the pure core-response builder (compute_drhoc keeps its data-slot semantics on top of it); assemble_v_sc now takes the q index explicitly instead of reading the drhoc slot through last_q_ (stale on the first assembly of a new q); per-displacement density stash (set/get_drho_disp) mirrors the existing v_sc stash. Verification: MODULE_DFPT serial suite 9/9 (the new terms are inert without nlcc pseudos); bn_cubic / bn_hex Gamma-vs-QE validation runs attached in the PR. --- source/source_esolver/esolver_dfpt_pw.cpp | 55 +++++- source/source_pw/module_dfpt/dfpt_phon.cpp | 10 ++ source/source_pw/module_dfpt/dfpt_phon.h | 48 ++++- .../source_pw/module_dfpt/dfpt_phon_elec.cpp | 170 +++++++++++++++++- source/source_pw/module_dfpt/dfpt_pw.h | 4 + source/source_pw/module_dfpt/dfpt_pw_data.cpp | 37 ++++ source/source_pw/module_dfpt/dfpt_pw_data.h | 10 ++ source/source_pw/module_dfpt/dfpt_pw_impl.h | 3 +- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 4 + source/source_pw/module_dfpt/dfpt_pw_q0.cpp | 2 +- source/source_pw/module_dfpt/dfpt_pw_run.cpp | 40 ++++- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 15 +- source/source_pw/module_dfpt/dfpt_rho.cpp | 69 +++---- source/source_pw/module_dfpt/dfpt_rho.h | 11 ++ .../module_dfpt/test_serial/CMakeLists.txt | 1 + .../test_serial/dfpt_phon_serial_test.cpp | 6 +- 16 files changed, 425 insertions(+), 60 deletions(-) diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 777a678d4a2..9066c9ffa99 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -423,10 +423,59 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) } } + // ground-state XC potential on the rho grid, evaluated on the SAME + // density the first-order kernel linearizes around (valence + core + // charge): consumed by the nlcc diagonal 2n+1 term (QE dynmatcc). + // Only needed when some species carries a nonlinear core correction. + bool nlcc_any = false; + for (int it = 0; it < ucell.ntype; ++it) + { + if (ucell.atoms[it].ncpp.nlcc) + { + nlcc_any = true; + break; + } + } + std::vector vxc_r; + if (nlcc_any) + { + const std::tuple etxc_vtxc_v + = XC_Functional::v_xc(this->pw_rho->nrxx, this->pelec->charge, &ucell, 1, + false, false, this->inp_->gga_grad, + XC_Functional::get_hybrid_alpha(), + XC_Functional::get_hse_omega()); + const ModuleBase::matrix& v_xc_mat = std::get<2>(etxc_vtxc_v); + if (v_xc_mat.nr > 0 && v_xc_mat.nc == this->pw_rho->nrxx) + { + vxc_r.resize(this->pw_rho->nrxx); + for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) + { + vxc_r[ir] = v_xc_mat(0, ir); + } + } + else + { + ModuleBase::WARNING("ESolver_DFPT_PW::init_dfpt", + "failed to evaluate the ground-state XC potential; " + "the nlcc diagonal dynamical-matrix term is skipped"); + } + } + this->stp.sync_cpu(); - dfpt_->init(ucell, *this->stp.psi_cpu, this->pw_rho, this->pw_wfc, &this->sf, veff_r, - this->pelec->wg, this->pelec->ekb, xc_adapter_, nelec_, ecutwfc_, - dft_plus_u_ ? this->dftu_.get() : nullptr); + const ModuleDFPT::DFPT_PW_InitContext ctx{&ucell, + this->stp.psi_cpu, + this->pw_rho, + this->pw_wfc, + &this->sf, + &veff_r, + &this->pelec->wg, + &this->pelec->ekb, + xc_adapter_, + nelec_, + ecutwfc_, + dft_plus_u_ ? this->dftu_.get() : nullptr, + vxc_r.empty() ? nullptr : &vxc_r}; + dfpt_->init(ctx); } void ESolver_DFPT_PW::run_post_process(UnitCell& ucell) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index de9fa67306c..4ecc541c744 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -64,6 +64,16 @@ void DFPT_Phon::init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, DFPT_Pert* per ModuleBase::timer::end("DFPT_Phon", "init"); } +void DFPT_Phon::set_rho(const DFPT_Rho* rho) +{ + rho_ = rho; +} + +void DFPT_Phon::set_vxc_r(const std::vector* vxc_r) +{ + vxc_r_ = vxc_r; +} + // --------------------------------------------------------------------------- // assemble / diagonalize / LO-TO / sum rule // --------------------------------------------------------------------------- diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 9bcd8dbaec8..3b34c018d71 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -5,7 +5,9 @@ #include "source_cell/unitcell.h" #include "source_psi/psi.h" +#include #include +#include namespace ModulePW { @@ -16,6 +18,7 @@ namespace ModuleDFPT { class DFPT_Pert; +class DFPT_Rho; /** * @brief Dynamical matrix of DFPT (C5). @@ -35,6 +38,17 @@ class DFPT_Pert; * row D[b][*] is filled right after displacement b converges, so the dpsi * storage never needs a direction dimension (data-layer refactor reserved * for phase B). + * + * nlcc systems add two direct terms (QE ph.x ground truth: dynmat_nlcc + * and dynmatcc, both absent from the wavefunction response): + * D_ab += (Omega/nrxx) sum_r dvxc^a(r) conj(drhoc^b(r)) + * (dvxc^a = f_xc . (drho_val^a + drhoc^a) passed in by the caller; the + * mixed core-valence and core-core kernel pieces) and the q-independent + * diagonal correction + * D_aa[i][j] -= Omega tpiba^2 sum_G G_i G_j rhoc~(|G|) + * exp(-i G tau_a) conj(vxc^0(G)) + * (second derivative of the shifted core charge contracted with the + * ground-state XC potential). */ class DFPT_Phon { @@ -44,18 +58,29 @@ class DFPT_Phon void init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, DFPT_Pert* pert); + /// providers for the nlcc direct terms: the first-order core response + /// (DFPT_Rho::compute_drhoc_g) and the ground-state XC potential on + /// the rho grid (for the diagonal term). Null/empty providers leave + /// the nlcc terms inert. + void set_rho(const DFPT_Rho* rho); + void set_vxc_r(const std::vector* vxc_r); + void assemble(int q_idx, DFPT_PW_Data& data); /// Fill the D[b][*] row of the electronic dynamical-matrix contribution /// for the converged displacement (atom_idx, dir); psi/wg are the - /// ground-state wavefunctions and occupations. Requires a wired + /// ground-state wavefunctions and occupations. dvxc_r is the first-order + /// XC potential of this displacement's TOTAL density response (valence + /// + core), computed by the caller from the analytic kernel; it gates + /// the nlcc kernel term (empty = no nlcc). Requires a wired /// DFPT_Pert (init); a null pert leaves the row untouched. void accumulate_electron(int q_idx, int atom_idx, int dir, const psi::Psi>& psi, const ModuleBase::matrix& wg, - DFPT_PW_Data& data); + DFPT_PW_Data& data, + const std::vector>& dvxc_r); void diagonalize(int q_idx, DFPT_PW_Data& data); @@ -138,6 +163,8 @@ class DFPT_Phon UnitCell* ucell_ = nullptr; ModulePW::PW_Basis* pw_rho_ = nullptr; DFPT_Pert* pert_ = nullptr; + const DFPT_Rho* rho_ = nullptr; + const std::vector* vxc_r_ = nullptr; double ewald_alpha_ = 0.0; double ewald_rcut_ = 0.0; @@ -153,6 +180,23 @@ class DFPT_Phon const psi::Psi>& psi, const ModuleBase::matrix& wg); + /// nlcc kernel term (QE dynmat_nlcc): for the converged displacement + /// (atom_idx, dir) with first-order XC potential dvxc_r, contract with + /// the core response of every displacement (iat2, idir2) and fill the + /// Hermitian pair (row, col) / (col, row). Gated on dvxc_r and the + /// providers set by set_rho. + void accum_nlcc_kernel(int atom_idx, + int dir, + int q_idx, + DFPT_PW_Data& data, + const std::vector>& dvxc_r); + + /// nlcc diagonal term (QE dynmatcc): q-independent same-atom block + /// -Omega tpiba^2 sum_G G_i G_j rhoc~(|G|) exp(-i G tau_a) + /// conj(vxc^0(G)) added into the (3a+i, 3a+j) entries. Gated on the + /// ground-state XC potential provider (set_vxc_r). + void accum_nlcc_diag(int atom_idx, int q_idx); + /// DFT+U contribution to the dynamical matrix (U0 reservation). void dftu_onsite(int q_idx, DFPT_PW_Data& data); diff --git a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp index 2cbef586ec7..10f7003e4a1 100644 --- a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp @@ -8,6 +8,9 @@ #include "dfpt_kq_basis.h" #include "dfpt_pert.h" #include "dfpt_pw_data.h" +#include "dfpt_rho.h" +#include "source_base/math_integral.h" +#include "source_base/math_sphbes.h" #include #include @@ -134,11 +137,12 @@ void scatter_chi_to_r(const ModulePW::PW_Basis& pw_rho, } // namespace void DFPT_Phon::accumulate_electron(int q_idx, - int atom_idx, - int dir, - const psi::Psi>& psi, - const ModuleBase::matrix& wg, - DFPT_PW_Data& data) + int atom_idx, + int dir, + const psi::Psi>& psi, + const ModuleBase::matrix& wg, + DFPT_PW_Data& data, + const std::vector>& dvxc_r) { ModuleBase::TITLE("DFPT_Phon", "accumulate_electron"); ModuleBase::timer::start("DFPT_Phon", "accumulate_electron"); @@ -186,11 +190,167 @@ void DFPT_Phon::accumulate_electron(int q_idx, } } + // ---- nlcc direct terms (QE dynmat_nlcc + dynmatcc) ---- + accum_nlcc_kernel(atom_idx, dir, q_idx, data, dvxc_r); + accum_nlcc_diag(atom_idx, q_idx); + // restore the converged dpsi of this displacement restore_dpsi(data, q_idx, psi.get_nk(), psi.get_nbands(), dpsib); ModuleBase::timer::end("DFPT_Phon", "accumulate_electron"); } +void DFPT_Phon::accum_nlcc_kernel(int atom_idx, + int dir, + int q_idx, + DFPT_PW_Data& data, + const std::vector>& dvxc_r) +{ + ModuleBase::TITLE("DFPT_Phon", "accum_nlcc_kernel"); + ModuleBase::timer::start("DFPT_Phon", "accum_nlcc_kernel"); + const int nrxx = pw_rho_->nrxx; + if (rho_ == nullptr || ucell_ == nullptr || static_cast(dvxc_r.size()) != nrxx) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); + return; + } + // any nlcc atom at all? + bool nlcc_any = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any = true; + break; + } + } + if (!nlcc_any) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); + return; + } + const int nat = ucell_->nat; + const ModuleBase::Vector3 q_cart = data.get_qvec(q_idx) * ucell_->G; + const double fac = ucell_->omega / static_cast(nrxx); + const int rowb = 3 * atom_idx + dir; + const double mass_row = ucell_->atoms[ucell_->iat2it[atom_idx]].mass; + std::vector> drhoc_g; + std::vector> drhoc_r(nrxx); + for (int iat2 = 0; iat2 < nat; ++iat2) + { + if (!ucell_->atoms[ucell_->iat2it[iat2]].ncpp.nlcc) + { + continue; + } + for (int idir2 = 0; idir2 < 3; ++idir2) + { + rho_->compute_drhoc_g(*ucell_, iat2, idir2, q_cart, drhoc_g); + if (static_cast(drhoc_g.size()) != pw_rho_->npw) + { + continue; + } + pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); + std::complex x(0.0, 0.0); + for (int ir = 0; ir < nrxx; ++ir) + { + x += dvxc_r[ir] * std::conj(drhoc_r[ir]); + } + x *= fac; + const int cola = 3 * iat2 + idir2; + const double mass_norm = std::sqrt(mass_row * ucell_->atoms[ucell_->iat2it[iat2]].mass); + dynmat_accum_(rowb, cola) += x / mass_norm; + dynmat_accum_(cola, rowb) += std::conj(x) / mass_norm; + } + } + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); +} + +void DFPT_Phon::accum_nlcc_diag(int atom_idx, int q_idx) +{ + ModuleBase::TITLE("DFPT_Phon", "accum_nlcc_diag"); + ModuleBase::timer::start("DFPT_Phon", "accum_nlcc_diag"); + const int npw = pw_rho_ ? pw_rho_->npw : 0; + const int nrxx = pw_rho_ ? pw_rho_->nrxx : 0; + if (ucell_ == nullptr || pw_rho_ == nullptr || vxc_r_ == nullptr + || static_cast(vxc_r_->size()) != nrxx || npw == 0 || accum_q_ != q_idx) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); + return; + } + const int it = ucell_->iat2it[atom_idx]; + if (!ucell_->atoms[it].ncpp.nlcc) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); + return; + } + // QE dynmatcc ground truth, evaluated as a real-space integral: + // D_ij += Int dr d2rho_c/(d tau_i d tau_j) (r) * v_xc^0(r) + // with the second-order shifted-core amplitude + // d2rho_c(G) = -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau_a) + // (a real field: W(-G) = conj(W(G))), so the contraction runs through + // the same recip2real machinery as every other real-space field and + // never touches the half-space G-sum convention directly. + const pseudo& p = ucell_->atoms[it].ncpp; + const int mesh = p.msh; + const ModuleBase::Vector3 tau_cart + = ucell_->atoms[it].tau[ucell_->iat2ia[atom_idx]] * ucell_->lat0; + std::vector aux(mesh); + std::vector> w_g(npw); + std::vector> w_r(nrxx); + double d3[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; + const double fac = ucell_->omega / static_cast(nrxx); + for (int i = 0; i < 3; ++i) + { + for (int j = i; j < 3; ++j) + { + std::fill(w_g.begin(), w_g.end(), std::complex(0.0, 0.0)); + for (int ig = 0; ig < npw; ++ig) + { + const ModuleBase::Vector3& gv = pw_rho_->gcar[ig]; + const double g2 = gv * gv; + // the G = 0 shell contributes nothing (G_i G_j factor) + if (g2 < 1.0e-12) + { + continue; + } + const double gphys = std::sqrt(g2) * ucell_->tpiba; + ModuleBase::Sphbes::Spherical_Bessel(mesh, p.r.data(), gphys, 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; + + const double phase = ModuleBase::TWO_PI * (gv * tau_cart) / ucell_->lat0; + const std::complex gtau(std::cos(phase), -std::sin(phase)); + w_g[ig] = -ucell_->tpiba * ucell_->tpiba * gv[i] * gv[j] * rhocg * gtau; + } + pw_rho_->recip2real(w_g.data(), w_r.data()); + double s = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + s += w_r[ir].real() * (*vxc_r_)[ir]; + } + d3[i][j] = fac * s; + if (j != i) + { + d3[j][i] = d3[i][j]; + } + } + } + const double mass = ucell_->atoms[it].mass; + const int base = 3 * atom_idx; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + dynmat_accum_(base + i, base + j) += d3[i][j] / mass; + } + } + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); +} + void DFPT_Phon::accum_d2_same_atom(int atom_idx, int dir, int iat, diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index ea611743b92..42a849baa98 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -43,6 +43,10 @@ struct DFPT_PW_InitContext double nelec; double ecutwfc; const Plus_U_Base* dftu; + /// ground-state XC potential on the rho grid (valence + core charge + /// input, the same density the kernel linearizes around); consumed by + /// the nlcc diagonal 2n+1 term (QE dynmatcc). Null disables that term. + const std::vector* vxc_r = nullptr; }; /** diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp index 529517c6212..d62998185ff 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp @@ -115,6 +115,43 @@ std::vector> DFPT_PW_Data::get_vsc_r(int atom_idx, int dir) return std::vector>(); } +void DFPT_PW_Data::set_drho_disp(int atom_idx, int dir, const std::vector>& rho) +{ + ModuleBase::TITLE("DFPT_PW_Data", "set_drho_disp"); + ModuleBase::timer::start("DFPT_PW_Data", "set_drho_disp"); + if (atom_idx < 0 || dir < 0 || dir >= 3) + { + ModuleBase::timer::end("DFPT_PW_Data", "set_drho_disp"); + return; + } + const size_t slot = static_cast(3 * atom_idx + dir); + if (slot >= drho_disp_.size()) + { + drho_disp_.resize(slot + 1); + } + drho_disp_[slot] = rho; + ModuleBase::timer::end("DFPT_PW_Data", "set_drho_disp"); +} + +std::vector> DFPT_PW_Data::get_drho_disp(int atom_idx, int dir) const +{ + ModuleBase::TITLE("DFPT_PW_Data", "get_drho_disp"); + ModuleBase::timer::start("DFPT_PW_Data", "get_drho_disp"); + if (atom_idx < 0 || dir < 0 || dir >= 3) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return std::vector>(); + } + const size_t slot = static_cast(3 * atom_idx + dir); + if (slot < drho_disp_.size()) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return drho_disp_[slot]; + } + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return std::vector>(); +} + void DFPT_PW_Data::set_dpsi_disp(int atom_idx, int dir, const std::vector>>>& d) diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h index 239c77fa4ad..79fdc824d36 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.h +++ b/source/source_pw/module_dfpt/dfpt_pw_data.h @@ -235,6 +235,13 @@ class DFPT_PW_Data void set_vsc_r(int atom_idx, int dir, const std::vector>& v); std::vector> get_vsc_r(int atom_idx, int dir) const; + /// converged first-order density response (G space) of displacement + /// (atom, dir): the two-pass 2n+1 accumulation needs it to rebuild the + /// nlcc kernel input f_xc.(drho_val + drhoc) per displacement, after + /// the working drho slots hold the last solved displacement + void set_drho_disp(int atom_idx, int dir, const std::vector>& rho); + std::vector> get_drho_disp(int atom_idx, int dir) const; + /// converged dpsi of displacement (atom, dir), indexed [k][band]; the /// two-pass 2n+1 accumulation reads it back after all displacements of /// the basis have been solved (the working dpsi slots get overwritten by @@ -302,6 +309,9 @@ class DFPT_PW_Data /// converged v_sc per displacement (atom, dir): [3*nat] entries std::vector>> vsc_r_; + /// converged first-order density per displacement (atom, dir): [3*nat] + std::vector>> drho_disp_; + /// converged dpsi per displacement (atom, dir): [3*nat][k][band] entries std::vector>>>> dpsi_disp_; diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index ed32d6ad795..fb22f876e59 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -65,7 +65,8 @@ class DFPT_PW::Impl /// v_sc^q = v_Hartree(drho_in) + xc(drho_in) on the shared grid; /// shared by solve_displacement and solve_efield_resp so the /// real-space Hartree+XC add loop is not duplicated - void assemble_v_sc(const ModuleBase::Vector3& q_cart, + void assemble_v_sc(int q_idx, + const ModuleBase::Vector3& q_cart, const std::vector>& drho_in_g, std::vector>& v_sc_r) const; diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 2b702514f6d..5e4a7d96191 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -120,6 +120,10 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, } rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, mix_type, mix_beta, kerker_a2}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); + // nlcc direct-term providers (QE dynmat_nlcc + dynmatcc equivalents): + // the pure core-response builder and the ground-state XC potential + phon_.set_rho(&rho_); + phon_.set_vxc_r(ctx.vxc_r); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); data_.init(&qlist_, nk, nbands, npw_max, nrxx, nspin, nat, ctx.dftu); diff --git a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp index 39dba1d3f8d..838d5b3391b 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp @@ -318,7 +318,7 @@ void DFPT_PW::Impl::solve_efield_resp(int q_idx) // (shared Impl::assemble_v_sc from dfpt_pw_solve.cpp) std::vector> v_sc_r; const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_in_g, v_sc_r); + assemble_v_sc(q_idx, q_cart, drho_in_g, v_sc_r); for (int ik = 0; ik < nk; ++ik) { if (static_cast(occ_kq_.size()) <= ik || occ_kq_[ik].empty()) diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 9db6f17e942..0e076887357 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -79,12 +79,50 @@ double DFPT_PW::Impl::run_displacement_irrep_pass(int q_idx, int irrep) worst = std::max(worst, residual); } } + // nlcc: the 2n+1 kernel term needs the first-order XC potential of + // each displacement's TOTAL response (valence + core). Precompute the + // nlcc flag once per pass; the per-displacement core response is + // recomputed here because the working drhoc slot holds the LAST + // solved displacement, and the valence response comes from the + // per-displacement stash. + bool nlcc_any = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any = true; + break; + } + } + const int npw = pw_rho_->npw; + const int nrxx = pw_rho_->nrxx; + const ModuleBase::Vector3 q_cart = data_.get_qvec(q_idx) * ucell_->G; for (int iat = 0; iat < nat; ++iat) { for (int idir = 0; idir < 3; ++idir) { + std::vector> dvxc_r; + if (nlcc_any && xc_ != nullptr) + { + std::vector> drhoc_a_g; + rho_.compute_drhoc_g(*ucell_, iat, idir, q_cart, drhoc_a_g); + const std::vector> drho_a_g = data_.get_drho_disp(iat, idir); + if (static_cast(drhoc_a_g.size()) == npw && static_cast(drho_a_g.size()) == npw) + { + std::vector> a_r(nrxx); + pw_rho_->recip2real(drho_a_g.data(), a_r.data()); + std::vector> drhoc_r(nrxx); + pw_rho_->recip2real(drhoc_a_g.data(), drhoc_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + a_r[ir] += drhoc_r[ir]; + } + xc_->set_context(q_cart); + xc_->apply(a_r, dvxc_r); + } + } // 2n+1 accumulation of this converged displacement - phon_.accumulate_electron(q_idx, iat, idir, gs_psi_, wg_, data_); + phon_.accumulate_electron(q_idx, iat, idir, gs_psi_, wg_, data_, dvxc_r); } } data_.add_residual(q_idx, irrep, worst); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index d467aa36221..751066ef559 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -24,7 +24,8 @@ namespace ModuleDFPT { -void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, +void DFPT_PW::Impl::assemble_v_sc(int q_idx, + const ModuleBase::Vector3& q_cart, const std::vector>& drho_in_g, std::vector>& v_sc_r) const { @@ -58,8 +59,10 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, 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_); + // above correctly uses the valence density alone). The q index is + // passed explicitly: the drhoc slot is per q and last_q_ can be stale + // on the first assembly of a new q. + const std::vector> drhoc_g = data_.get_drhoc_g(q_idx); if (static_cast(drhoc_g.size()) == pw_rho_->npw) { std::vector> drhoc_r(nrxx); @@ -137,6 +140,10 @@ void DFPT_PW::Impl::stash_converged_disp_response(int q_idx, int nbands) { data_.set_vsc_r(iat, idir, v_sc_r_last); + // stash the converged density response too: the two-pass accumulation + // rebuilds f_xc.(drho_val + drhoc) per displacement after the working + // drho slots have moved on to the next displacement + data_.set_drho_disp(iat, idir, data_.get_drho_g(q_idx, 0)); std::vector>>> disp( nk, std::vector>>(nbands)); @@ -233,7 +240,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) { std::vector> v_sc_r; const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_in_g, v_sc_r); + assemble_v_sc(q_idx, q_cart, drho_in_g, v_sc_r); if (dbg) { debug_iter_snapshot(iter, drho_in_g, v_sc_r, npw, nrxx); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 239c7804ab8..2bf3beea6ed 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -253,47 +253,51 @@ void DFPT_Rho::compute_drhoc(const UnitCell& ucell, { ModuleBase::TITLE("DFPT_Rho", "compute_drhoc"); ModuleBase::timer::start("DFPT_Rho", "compute_drhoc"); + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; + std::vector> drhoc_g; + compute_drhoc_g(ucell, iat, dir, q_cart, drhoc_g); + data.set_drhoc_g(q_idx, drhoc_g); + + const int it = ucell.iat2it[iat]; + // report magnitude for diagnostics + double norm2 = 0.0; + for (const std::complex& c : drhoc_g) + { + norm2 += std::norm(c); + } + 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::compute_drhoc_g(const UnitCell& ucell, + int iat, + int dir, + const ModuleBase::Vector3& q_cart, + std::vector>& drhoc_g) const +{ const int npw = pw_rho_ ? pw_rho_->npw : 0; + drhoc_g.clear(); 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)); + drhoc_g.assign(npw, std::complex(0.0, 0.0)); std::vector aux(mesh); for (int ig = 0; ig < npw; ++ig) @@ -325,21 +329,6 @@ void DFPT_Rho::compute_drhoc(const UnitCell& ucell, 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, diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 83ae0cda250..caa28dba2ad 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -131,6 +131,17 @@ class DFPT_Rho int q_idx, DFPT_PW_Data& data); + /// same first-order core-density amplitude as compute_drhoc, but pure: + /// takes the cartesian q explicitly and writes into a caller-owned + /// vector instead of the DFPT_PW_Data slot. Used by the 2n+1 nlcc + /// direct terms, which need the core response of EVERY displacement + /// while the data slot holds only the current one. + void compute_drhoc_g(const UnitCell& ucell, + int iat, + int dir, + const ModuleBase::Vector3& q_cart, + std::vector>& drhoc_g) const; + /// first-order occupation matrix (docc) for DFT+U (U0 reservation). void cal_docc(const psi::Psi>& psi, const ModuleBase::matrix& wg, diff --git a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt index 6e4171510e5..d47a1de06e1 100644 --- a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt +++ b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt @@ -88,6 +88,7 @@ addtest( ../dfpt_phon.cpp ../dfpt_phon_ewald.cpp ../dfpt_phon_elec.cpp + ../dfpt_rho.cpp ../dfpt_pert.cpp ../dfpt_pert_vkb.cpp ../dfpt_pert_nl.cpp 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 085de1a93da..3f93a20334f 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 @@ -347,7 +347,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronAnalyticContraction) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 1, psi, wg, data_); + phon_.accumulate_electron(0, 0, 1, psi, wg, data_, std::vector>()); // expected: row 1 (atom 0, dir 1) of the Hermitian 2n+1 accumulation: // the row element receives wg* once per off-diagonal @@ -400,7 +400,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronD2GateOffGenericQ) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 0, psi, wg, data_); + phon_.accumulate_electron(0, 0, 0, psi, wg, data_, std::vector>()); const std::vector> g0(1, ModuleBase::Vector3(0.0, 0.0, 0.0)); for (int adir = 0; adir < 3; ++adir) @@ -479,7 +479,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronD2CommensurateQ) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 1, psi, wg, data_); + phon_.accumulate_electron(0, 0, 1, psi, wg, data_, std::vector>()); for (int adir = 0; adir < 3; ++adir) { From 045f0504505ac96882155c1c29dfb6541e624a26 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Wed, 28 Oct 2026 21:30:12 +0800 Subject: [PATCH 15/21] fix(dfpt): the working nlcc direct-term completion (validated on bn_cubic) Three fixes over the first cut, pinned by the element-level matrix comparison against QE dyn.G (bn_cubic: all blocks now agree to 2e-4, TO frequencies to 0.007 cm-1): - the per-q kernel table must fill EVERY displacement slot: rows of non-nlcc atoms need their dvaux too (f_xc applied to the valence response alone) -- leaving their slots empty silently dropped the entire one-sided cross-block X(N,B), the dominant mixed piece - the effective Hermitian element is the symmetrized (1/2)(X_ab + conj(X_ba)) (QE dynmat_nlcc single ordering + final Hermitization); the earlier both-order full fill doubled the same-atom blocks, and the lower-triangle-only variant overshot the one-sided cross blocks - the vxc provider must live in Impl-owned storage (the init context points at esolver stack state); the diagonal block is added once per atom (dir == 0) --- source/source_pw/module_dfpt/dfpt_phon.cpp | 33 ++++++++ source/source_pw/module_dfpt/dfpt_phon.h | 7 ++ .../source_pw/module_dfpt/dfpt_phon_elec.cpp | 81 +++++++++++++++---- source/source_pw/module_dfpt/dfpt_pw_impl.h | 4 + source/source_pw/module_dfpt/dfpt_pw_init.cpp | 9 ++- source/source_pw/module_dfpt/dfpt_pw_run.cpp | 65 +++++++++++---- 6 files changed, 168 insertions(+), 31 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 4ecc541c744..1a414dbc562 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -74,6 +75,11 @@ void DFPT_Phon::set_vxc_r(const std::vector* vxc_r) vxc_r_ = vxc_r; } +void DFPT_Phon::set_dvxc_core_table(const std::vector>>* table) +{ + dvxc_core_ = table; +} + // --------------------------------------------------------------------------- // assemble / diagonalize / LO-TO / sum rule // --------------------------------------------------------------------------- @@ -121,6 +127,33 @@ void DFPT_Phon::assemble(int q_idx, DFPT_PW_Data& data) } } data.set_dynmat(q_idx, dyn); + // TEMPORARY QE-format dump for the nlcc element-level comparison + // (unreduced force constants; remove after the residual localization) + if (getenv("DFPT_DUMP_DYN") != nullptr) + { + std::ostringstream name; + name << "abacus_dyn.q" << q_idx << ".txt"; + std::ofstream out(name.str()); + out << std::setprecision(8) << std::scientific; + for (int ia = 0; ia < nat; ++ia) + { + for (int ja = 0; ja < nat; ++ja) + { + const double mm = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass + * ucell_->atoms[ucell_->iat2it[ja]].mass); + out << " " << ia + 1 << " " << ja + 1 << "\n"; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + const std::complex v = dyn(3 * ia + i, 3 * ja + j) * mm; + out << " " << v.real() << " " << v.imag(); + } + out << "\n"; + } + } + } + } dynmat_accum_ = ModuleBase::ComplexMatrix(); accum_q_ = -1; ModuleBase::timer::end("DFPT_Phon", "assemble"); diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 3b34c018d71..cffd3695d03 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -65,6 +65,12 @@ class DFPT_Phon void set_rho(const DFPT_Rho* rho); void set_vxc_r(const std::vector* vxc_r); + /// per-q table of f_xc.drhoc for every displacement slot [3*nat] + /// (the kernel term contracts it against the stashed valence + /// responses of the mixed piece). Set by the solver once per q, + /// before the accumulation pass; owned by the caller. + void set_dvxc_core_table(const std::vector>>* table); + void assemble(int q_idx, DFPT_PW_Data& data); /// Fill the D[b][*] row of the electronic dynamical-matrix contribution @@ -165,6 +171,7 @@ class DFPT_Phon DFPT_Pert* pert_ = nullptr; const DFPT_Rho* rho_ = nullptr; const std::vector* vxc_r_ = nullptr; + const std::vector>>* dvxc_core_ = nullptr; double ewald_alpha_ = 0.0; double ewald_rcut_ = 0.0; diff --git a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp index 10f7003e4a1..184ff8d51e2 100644 --- a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp @@ -192,7 +192,12 @@ void DFPT_Phon::accumulate_electron(int q_idx, // ---- nlcc direct terms (QE dynmat_nlcc + dynmatcc) ---- accum_nlcc_kernel(atom_idx, dir, q_idx, data, dvxc_r); - accum_nlcc_diag(atom_idx, q_idx); + if (dir == 0) + { + // the diagonal block is direction-pair complete: add it once per + // atom (every accumulate_electron call would triple it) + accum_nlcc_diag(atom_idx, q_idx); + } // restore the converged dpsi of this displacement restore_dpsi(data, q_idx, psi.get_nk(), psi.get_nbands(), dpsib); @@ -233,32 +238,76 @@ void DFPT_Phon::accum_nlcc_kernel(int atom_idx, const double fac = ucell_->omega / static_cast(nrxx); const int rowb = 3 * atom_idx + dir; const double mass_row = ucell_->atoms[ucell_->iat2it[atom_idx]].mass; + // QE dynmat_nlcc semantics, pinned by the element-level comparison + // against QE dyn.G (bn_cubic): the raw matrix carries + // X_ab = (Omega/nrxx) sum_r dvxc^a(r) conj(drhoc^b(r)), + // dvxc^a = f_xc.(drho_val^a + drhoc^a), once per element, and the + // effective Hermitian value is the symmetrized + // S_ab = (1/2)(X_ab + conj(X_ba)). + // X_ab is two-sided only when BOTH displacements carry core response; + // for the one-sided cross blocks (e.g. B vs N: X_BN = 0 because N has + // no nlcc) this reduces to (1/2) X_NB. The column loop must therefore + // cover non-nlcc atoms too: their X_ab vanishes but the transposed + // piece conj(X_ba) = sum_r dvxc^b conj(drhoc^a) -- the column's + // valence-only dvaux against the ROW's core -- is exactly what the + // (nlcc-row, non-nlcc-column) element needs. std::vector> drhoc_g; std::vector> drhoc_r(nrxx); + // row's own core response in real space (for the conj(X_ba) piece; + // empty when the row atom carries no nlcc) + std::vector> drhoc_a_r; + rho_->compute_drhoc_g(*ucell_, atom_idx, dir, q_cart, drhoc_g); + const bool row_has_core = static_cast(drhoc_g.size()) == pw_rho_->npw; + if (row_has_core) + { + drhoc_a_r.resize(nrxx); + pw_rho_->recip2real(drhoc_g.data(), drhoc_a_r.data()); + } + const bool have_table = dvxc_core_ != nullptr; for (int iat2 = 0; iat2 < nat; ++iat2) { - if (!ucell_->atoms[ucell_->iat2it[iat2]].ncpp.nlcc) - { - continue; - } + const bool col_has_core = ucell_->atoms[ucell_->iat2it[iat2]].ncpp.nlcc; for (int idir2 = 0; idir2 < 3; ++idir2) { - rho_->compute_drhoc_g(*ucell_, iat2, idir2, q_cart, drhoc_g); - if (static_cast(drhoc_g.size()) != pw_rho_->npw) + // X_ab: the row's dvaux against the column's core response + std::complex x_ab(0.0, 0.0); + if (col_has_core) { - continue; + rho_->compute_drhoc_g(*ucell_, iat2, idir2, q_cart, drhoc_g); + if (static_cast(drhoc_g.size()) != pw_rho_->npw) + { + continue; + } + pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + x_ab += dvxc_r[ir] * std::conj(drhoc_r[ir]); + } + x_ab *= fac; } - pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); - std::complex x(0.0, 0.0); - for (int ir = 0; ir < nrxx; ++ir) + // X_ba: the column's dvaux (per-q table) against the row's + // core response + std::complex x_ba(0.0, 0.0); + const int cola = 3 * iat2 + idir2; + if (row_has_core && have_table && cola < static_cast(dvxc_core_->size()) + && static_cast((*dvxc_core_)[cola].size()) == nrxx) + { + const std::vector>& dvxc_b = (*dvxc_core_)[cola]; + for (int ir = 0; ir < nrxx; ++ir) + { + x_ba += dvxc_b[ir] * std::conj(drhoc_a_r[ir]); + } + x_ba *= fac; + } + if (x_ab == std::complex(0.0, 0.0) && x_ba == std::complex(0.0, 0.0)) { - x += dvxc_r[ir] * std::conj(drhoc_r[ir]); + continue; } - x *= fac; - const int cola = 3 * iat2 + idir2; + const std::complex s_ab = 0.5 * (x_ab + std::conj(x_ba)); const double mass_norm = std::sqrt(mass_row * ucell_->atoms[ucell_->iat2it[iat2]].mass); - dynmat_accum_(rowb, cola) += x / mass_norm; - dynmat_accum_(cola, rowb) += std::conj(x) / mass_norm; + // single-ordering fill: the transposed element receives its + // own S_ba = conj(S_ab) from the column displacement's pass + dynmat_accum_(rowb, cola) += s_ab / mass_norm; } } ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index fb22f876e59..d515e9cbb00 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -158,6 +158,10 @@ class DFPT_PW::Impl ModulePW::PW_Basis_K* pw_wfc_ = nullptr; Structure_Factor* sf_ = nullptr; std::vector veff_r_; + /// ground-state XC potential (valence + core input), copied from the + /// init context: the nlcc diagonal 2n+1 term consumes it; the copy + /// owns the storage because the context points at caller stack state + std::vector vxc_r_; ModuleBase::matrix wg_; ModuleBase::matrix eig_; const XC_First_Order* xc_ = nullptr; diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 5e4a7d96191..6705da05e64 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -123,7 +123,9 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, // nlcc direct-term providers (QE dynmat_nlcc + dynmatcc equivalents): // the pure core-response builder and the ground-state XC potential phon_.set_rho(&rho_); - phon_.set_vxc_r(ctx.vxc_r); + // the Impl-owned copy: the context pointer refers to caller stack + // state that dies when the esolver init returns + phon_.set_vxc_r(vxc_r_.empty() ? nullptr : &vxc_r_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); data_.init(&qlist_, nk, nbands, npw_max, nrxx, nspin, nat, ctx.dftu); @@ -149,6 +151,11 @@ void DFPT_PW::init(const DFPT_PW_InitContext& ctx) pimpl_->pw_wfc_ = pw_wfc; pimpl_->sf_ = sf; pimpl_->veff_r_ = *ctx.veff_r; + pimpl_->vxc_r_.clear(); + if (ctx.vxc_r != nullptr) + { + pimpl_->vxc_r_ = *ctx.vxc_r; + } pimpl_->wg_ = *ctx.wg; pimpl_->eig_ = *ctx.eig; pimpl_->xc_ = ctx.xc; diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 0e076887357..3020cda2dc3 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -97,28 +97,65 @@ double DFPT_PW::Impl::run_displacement_irrep_pass(int q_idx, int irrep) const int npw = pw_rho_->npw; const int nrxx = pw_rho_->nrxx; const ModuleBase::Vector3 q_cart = data_.get_qvec(q_idx) * ucell_->G; - for (int iat = 0; iat < nat; ++iat) + // per-q table of f_xc.(drho_val + drhoc) for every displacement slot: + // the total-response first-order XC potential that the Hermitized + // kernel term contracts against the core responses (QE dynmat_nlcc + // semantics; built once, before the accumulation pass) + std::vector>> dvxc_core; + if (nlcc_any && xc_ != nullptr) { - for (int idir = 0; idir < 3; ++idir) + // fill EVERY displacement slot: the kernel term X(i,j) contracts + // the row's dvaux with the COLUMN's core response, so rows of + // non-nlcc atoms need their dvaux too (f_xc applied to their + // valence response alone -- they have no core piece). Leaving + // their slots empty silently dropped the entire X(N,B) cross + // block (the dominant mixed piece for a B-N cell). + dvxc_core.resize(3 * nat); + std::vector> drhoc_g; + for (int iat = 0; iat < nat; ++iat) { - std::vector> dvxc_r; - if (nlcc_any && xc_ != nullptr) + const bool iat_nlcc = ucell_->atoms[ucell_->iat2it[iat]].ncpp.nlcc; + for (int idir = 0; idir < 3; ++idir) { - std::vector> drhoc_a_g; - rho_.compute_drhoc_g(*ucell_, iat, idir, q_cart, drhoc_a_g); + std::vector> a_r(nrxx, std::complex(0.0, 0.0)); + if (iat_nlcc) + { + rho_.compute_drhoc_g(*ucell_, iat, idir, q_cart, drhoc_g); + if (static_cast(drhoc_g.size()) == npw) + { + pw_rho_->recip2real(drhoc_g.data(), a_r.data()); + } + } + // add the stashed converged valence response of this + // displacement (the working slots hold the last one) const std::vector> drho_a_g = data_.get_drho_disp(iat, idir); - if (static_cast(drhoc_a_g.size()) == npw && static_cast(drho_a_g.size()) == npw) + if (static_cast(drho_a_g.size()) == npw) { - std::vector> a_r(nrxx); - pw_rho_->recip2real(drho_a_g.data(), a_r.data()); - std::vector> drhoc_r(nrxx); - pw_rho_->recip2real(drhoc_a_g.data(), drhoc_r.data()); + std::vector> dval_r(nrxx); + pw_rho_->recip2real(drho_a_g.data(), dval_r.data()); for (int ir = 0; ir < nrxx; ++ir) { - a_r[ir] += drhoc_r[ir]; + a_r[ir] += dval_r[ir]; } - xc_->set_context(q_cart); - xc_->apply(a_r, dvxc_r); + } + xc_->set_context(q_cart); + xc_->apply(a_r, dvxc_core[3 * iat + idir]); + } + } + phon_.set_dvxc_core_table(&dvxc_core); + } + for (int iat = 0; iat < nat; ++iat) + { + for (int idir = 0; idir < 3; ++idir) + { + std::vector> dvxc_r; + if (nlcc_any && xc_ != nullptr) + { + const int slot = 3 * iat + idir; + if (slot < static_cast(dvxc_core.size()) + && static_cast(dvxc_core[slot].size()) == nrxx) + { + dvxc_r = dvxc_core[slot]; } } // 2n+1 accumulation of this converged displacement From 9f0cedcfd6cf2020f0f661b3c20b8e9b188cb8bd Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Wed, 28 Oct 2026 22:54:17 +0800 Subject: [PATCH 16/21] fix(dfpt): the core-response structure-factor phase was mirrored compute_drhoc_g built the phase exp(+i(G+q).tau), which places the first-order core density at the MIRROR site -tau instead of the atom position tau (QE addcore: exp(-iG.tau) * exp(-i q.tau)). Every validation so far was structurally blind to the sign: bn_cubic's B sits at the origin (both conventions identical) and N carries no nlcc, and the F2 numerical cross-check ran on the same blind system. The bn_hex 16-block matrix comparison against QE dyn.G exposed it element by element: with wurtzite B1=(1/3,2/3,1/4), the mirror of B1 lands exactly on N2=(2/3,1/3,3/4) and the mirror of B2 on N1, so the mirrored core response produced spurious same-layer overlaps ((B1,B2)xx +0.053 vs the true 0.0013 interlayer scale) and lost the real diagonal ones ((B1,B1)xx -0.105), collapsing the E2low interlayer shear to -399 cm^-1 against QE's -29.8. After the one-line conjugation, bn_hex stays bit-identical (B at the origin) and bn_hex matches QE on every block to <= 2.5e-4 (B-involving) with the optical rms dropping from 28.7 (C5, no core response) and ~60 (previous iterations) to 1.2 cm^-1; E2low recovers to -32.2 vs -29.8 and the stiff modes to within 0.8 cm^-1. --- source/source_pw/module_dfpt/dfpt_rho.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 6a11ee558e1..00d6cfe09d9 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -356,9 +356,15 @@ void DFPT_Rho::compute_drhoc_g(const UnitCell& ucell, 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) + // structure factor phase: exp(-i (G+q) . tau_iat). The sign is + // load-bearing: exp(+i(G+q)tau) mirrors the core response to the + // -tau site (for wurtzite B the mirror lands exactly on the N + // sites of the other layer), which bn_cubic cannot expose (B at + // the origin makes both conventions identical) but which showed + // up element-by-element in the bn_hex matrix comparison against + // QE dyn.G. QE addcore: e^{-i G tau} * e^{-i q tau}. const double phase = ModuleBase::TWO_PI * (w * tau_cart) / ucell.lat0; - const std::complex gtau(std::cos(phase), std::sin(phase)); + const std::complex gtau(std::cos(phase), -std::sin(phase)); // direction derivative: -i * tpiba * (G+q)_dir const double wdir = ucell.tpiba * w[dir]; From c69896de947056aa1e7e9a589c76a1264273cbdd Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Thu, 29 Oct 2026 18:32:13 +0800 Subject: [PATCH 17/21] test(dfpt): nlcc phase-localization and fill-semantics regressions Two-level regression suite replacing the temporary matrix dump (whose diagnostic role -- the element-level QE comparison that pinned every formula bug -- now lives as executable invariants): Level 1 (dfpt_rho_serial_test): DrhocLocalizesAtTheAtomNotItsMirror. An nlcc atom at the non-origin tau = (0.25, 0.31, 0.17) (the anti-blind-spot position: every structure-factor phase is active) must have its first-order core response |drhoc|^2 centroid AT the atom. The mirrored phase exp(+i(G+q)tau) -- the bug every historical validation missed because bn_cubic's B sits at the origin -- places it at the mirror site instead. Includes a generic non-zero q variant. Level 2 (dfpt_phon_serial_test): NlccDirectTermsFillSemantics. With a synthetic all-slots kernel table and no dpsi (isolating the kernel contribution from the bare cross term), the one-sided cross block must carry exactly (1/2) conj(X_ba) (the symmetrized fill pinned by the element-level QE comparison), and the accumulated matrix must stay Hermitian. Together these guard against the three fill bugs fixed in the parent commits (dead coreless rows, doubled same-atom blocks, and the halved one-sided cross blocks). The temporary DFPT_DUMP_DYN production block is removed (its commit 045f05045 in the branch history retains the identical code for local re-use during Project E matrix comparisons). --- .../test_serial/dfpt_phon_serial_test.cpp | 96 ++++++++++++ .../test_serial/dfpt_rho_serial_test.cpp | 147 ++++++++++++++++++ .../test_serial/dfpt_serial_fixture.cpp | 48 ++++++ .../test_serial/dfpt_serial_fixture.h | 15 ++ 4 files changed, 306 insertions(+) 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 3f93a20334f..9ced75800fe 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 @@ -12,6 +12,7 @@ #include "dfpt_serial_fixture.h" #include "source_base/complexmatrix.h" +#include "source_pw/module_dfpt/dfpt_rho.h" #include "source_base/constants.h" #include "source_base/matrix3.h" #include "source_base/vector3.h" @@ -964,3 +965,98 @@ TEST_F(DFPTPhonSerialTest, AsrSimpleCrossCartesianRegression) EXPECT_EQ(out(i, i).imag(), 0.0) << "i=" << i; } } + +// --------------------------------------------------------------------------- +// nlcc direct-term fill-semantics regression: with a synthetic provider set +// (all-atoms kernel table, ground-state XC vector, DFPT_Rho), the nlcc +// kernel accumulation must satisfy the QE-pinned invariants: +// (1) the one-sided cross block (nlcc row, coreless column) receives +// (1/2)(X_ab + conj(X_ba)) = (1/2) X_ab [X_ba = 0: coreless row] +// (2) every displacement slot's dvaux participates (the "dead row" bug +// dropped the entire X(N, B) cross block when the table left +// coreless atoms empty) +// (3) the matrix stays Hermitian +// These pin the fill factor (the double-sided same-atom block would double +// under a both-order fill; a lower-triangle-only fill halves the one-sided +// cross blocks). +// --------------------------------------------------------------------------- +TEST_F(DFPTPhonSerialTest, NlccDirectTermsFillSemantics) +{ + MakeNlccTwoAtomCell(); + SetupPhon(k_d_, q_d_); + + ModuleDFPT::DFPT_Rho rho; + rho.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_plain"), 0.4, 1.0, 8}); + phon_.set_rho(&rho); + + const int nrxx = pw_rho_.nrxx; + // synthetic kernel table: f_xc.drhoc stand-ins = distinct constant + // complex fields per slot (values chosen to survive exact comparison) + std::vector>> table(6); + for (int slot = 0; slot < 6; ++slot) + { + const double re = 0.1 + 0.05 * slot; + const double im = 0.02 * slot - 0.05; + table[slot].assign(nrxx, std::complex(re, im)); + } + phon_.set_dvxc_core_table(&table); + + // synthetic ground-state XC vector (only the diag term consumes it); + // zero it so the diag term's only effect is absent (this test targets + // the kernel fill) + const std::vector vxc_zero(nrxx, 0.0); + phon_.set_vxc_r(&vxc_zero); + + psi::Psi> psi = MakeSinglePlaneWavePsi(k_d_); + ModuleBase::matrix wg = MakeOccWeights(); + // NO dpsi injection: the bare cross term needs occupied dpsi entries + // and skips empty slots, so the accumulated matrix isolates the nlcc + // kernel contribution exactly + + phon_.accumulate_electron(0, 1, 0, psi, wg, data_, table[3]); + + const auto& m = phon_.dynmat_accum(); + // row (1, dir 0) = 3*1+0 = 3; column (atom 0, dir j) = j (coreless) + // X_3j = sum_r table[3](r) * conj(drhoc(atom0, dir_j)(r)) = 0 + // (atom 0 has no nlcc: drhoc = empty -> the kernel term's x_ab arm + // is zero for coreless COLUMNS) + // x_ba arm = sum_r table[j](r) * conj(drhoc(atom1, dir0)(r)) != 0 + // (row atom 1 has a core; the column's dvaux comes from the table) + // So the (3, j) element = (1/2)(0 + conj(x_ba)) = (1/2) conj(x_ba). + // The x_ba contraction runs over the actual Gaussian core field, so + // verify structurally: the one-sided block is NONZERO (dead-row bug + // gave exactly zero) and matches an independent recomputation. + const double fac = ucell_.omega / nrxx; + std::vector> drhoc_g; + rho.compute_drhoc_g(ucell_, 1, 0, q_cart_, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + std::vector> drhoc_r(nrxx); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + const double m0 = std::sqrt(ucell_.atoms[0].mass); + const double m1 = std::sqrt(ucell_.atoms[1].mass); + for (int j = 0; j < 3; ++j) + { + std::complex x_ba(0.0, 0.0); + for (int ir = 0; ir < nrxx; ++ir) + { + x_ba += table[j][ir] * std::conj(drhoc_r[ir]); + } + x_ba *= fac; + const std::complex expect = 0.5 * std::conj(x_ba) / (m0 * m1); + EXPECT_NEAR(std::abs(m(3, j) - expect), 0.0, 1.0e-10 * (1.0 + std::abs(expect))) + << "j=" << j << " got " << m(3, j) << " expect " << expect; + } + + // Hermiticity of the accumulated block so far: accumulate the column + // displacement (atom 1, dir 0 already done as row; do (atom 1, dir 1) + // to complete more of the same-atom block, then check the whole + // matrix is Hermitian within the accumulated rows + phon_.accumulate_electron(0, 1, 1, psi, wg, data_, table[4]); + // the (3, 4) and (4, 3) elements must be conjugate pairs: from row 3's + // pass, (3,4) = (1/2)(X_34 + conj(X_43)); from row 4's pass, (4,3) = + // (1/2)(X_43 + conj(X_34)) = conj((3,4)) + const std::complex e34 = m(3, 4); + const std::complex e43 = m(4, 3); + EXPECT_NEAR(e43.real(), e34.real(), 1.0e-12); + EXPECT_NEAR(e43.imag(), -e34.imag(), 1.0e-12); +} diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 66c0543f4a6..7d8b180c770 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -660,3 +660,150 @@ TEST_F(DFPTRhoSerialTest, BroydenResetMixingAllowsFreshQAndSameQCycle) EXPECT_NEAR(mixed_q0[ig].imag(), 0.4 * out[ig].imag(), 1.0e-12); } } + +// --------------------------------------------------------------------------- +// nlcc phase regression: the first-order core response must localize at the +// ATOM site, not its mirror. This is the anti-blind-spot test for the +// structure-factor sign: exp(+i(G+q)tau) mirrors the response to -tau, a bug +// every historical validation missed because the nlcc atom sat at the origin +// (bn_cubic) where both conventions coincide. The wurtzite bn_hex comparison +// against QE dyn.G pinned the sign: the mirror of B1 lands exactly on an N +// site of the other layer, collapsing the interlayer shear to -399 cm^-1. +// --------------------------------------------------------------------------- +TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) +{ + MakeNlccTwoAtomCell(); + SetupBases(k_d_, q_d_, 2); + ModuleDFPT::DFPT_Rho::Config cfg; + cfg.nspin = 1; + cfg.nrxx = pw_rho_.nrxx; + cfg.pw_rho = &pw_rho_; + cfg.pw_wfc = &pw_wfc_; + cfg.recip_matrix = G_; + cfg.mix_type = "rho_plain"; + cfg.mix_beta = 0.4; + cfg.kerker_a2 = 1.0; + cfg.broyden_ndim = 8; + ModuleDFPT::DFPT_Rho rho; + rho.init(cfg); + + // atom 1 carries the nlcc at the non-origin tau_nlcc_; all three + // displacement directions must localize at the atom + const int iat = 1; + const int it = ucell_.iat2it[iat]; + ASSERT_TRUE(ucell_.atoms[it].ncpp.nlcc); + // tau is stored in lat0-cartesian units (fixture convention: cell edge + // a_ in lat0 units); the true fractional position divides by a_ + const double inv_a = 1.0 / a_; + const ModuleBase::Vector3 tau_frac(ucell_.atoms[it].tau[ucell_.iat2ia[iat]].x * inv_a, + ucell_.atoms[it].tau[ucell_.iat2ia[iat]].y * inv_a, + ucell_.atoms[it].tau[ucell_.iat2ia[iat]].z * inv_a); + // mirror position folded back into the cell + ModuleBase::Vector3 mirror_frac(1.0 - tau_frac.x, 1.0 - tau_frac.y, 1.0 - tau_frac.z); + + const int nrxx = pw_rho_.nrxx; + std::vector> drhoc_g; + std::vector> drhoc_r(nrxx); + for (int dir = 0; dir < 3; ++dir) + { + // q = 0 first: the amplitude is the plain derivative profile, its + // |.|^2 centroid is at the atom (up to the derivative skew) + const ModuleBase::Vector3 q_zero(0.0, 0.0, 0.0); + rho.compute_drhoc_g(ucell_, iat, dir, q_zero, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + double wsum = 0.0; + double cx = 0.0; + double cy = 0.0; + double cz = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double wgt = std::norm(drhoc_r[ir]); + // fractional coordinates of the grid point via the FFT triple + const int iz = ir % pw_rho_.nz; + const int iy = (ir / pw_rho_.nz) % pw_rho_.ny; + const int ix = ir / (pw_rho_.ny * pw_rho_.nz); + // wrap-around: grid index -> fractional in [0, 1) + const double fx = static_cast(ix) / pw_rho_.nx; + const double fy = static_cast(iy) / pw_rho_.ny; + const double fz = static_cast(iz) / pw_rho_.nz; + // weight by |field|^2, folding the +1 wrap for the periodic image + // closest to tau: the field is localized, so the unwrapped + // centroid near tau is what we need + double wx = fx; + double wy = fy; + double wz = fz; + if (std::abs(fx - tau_frac.x) > 0.5) + { + wx += (fx < tau_frac.x) ? 1.0 : -1.0; + } + if (std::abs(fy - tau_frac.y) > 0.5) + { + wy += (fy < tau_frac.y) ? 1.0 : -1.0; + } + if (std::abs(fz - tau_frac.z) > 0.5) + { + wz += (fz < tau_frac.z) ? 1.0 : -1.0; + } + wsum += wgt; + cx += wgt * wx; + cy += wgt * wy; + cz += wgt * wz; + } + ASSERT_GT(wsum, 0.0); + const ModuleBase::Vector3 com(cx / wsum, cy / wsum, cz / wsum); + const double dist_atom = (com - tau_frac) * (com - tau_frac); + const double dist_mirror = (com - mirror_frac) * (com - mirror_frac); + // the centroid of a derivative field sits at the atom up to the + // derivative skew (~sigma); the mirror is a full cell away + EXPECT_LT(std::sqrt(dist_atom), 0.1) << "dir=" << dir; + EXPECT_GT(dist_mirror, 100.0 * dist_atom) << "dir=" << dir; + } + + // non-zero q variant: the q-phase exp(-i q.tau) must not mirror either. + // Use a generic incommensurate q (the fixture default q_d_ is already + // generic); the real-space field of the q-shifted amplitude still + // localizes at the atom (the e^{iqr} carrier is smooth on the sigma + // scale of the core). + const ModuleBase::Vector3 q2 = q_d_ * 0.7 + ModuleBase::Vector3(0.11, -0.07, 0.05); + const ModuleBase::Vector3 q2_cart = q2 * ucell_.G; + rho.compute_drhoc_g(ucell_, iat, 1, q2_cart, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + double wsum = 0.0; + double cx = 0.0; + double cy = 0.0; + double cz = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double wgt = std::norm(drhoc_r[ir]); + const int iz = ir % pw_rho_.nz; + const int iy = (ir / pw_rho_.nz) % pw_rho_.ny; + const int ix = ir / (pw_rho_.ny * pw_rho_.nz); + double wx = static_cast(ix) / pw_rho_.nx; + double wy = static_cast(iy) / pw_rho_.ny; + double wz = static_cast(iz) / pw_rho_.nz; + if (std::abs(wx - tau_frac.x) > 0.5) + { + wx += (wx < tau_frac.x) ? 1.0 : -1.0; + } + if (std::abs(wy - tau_frac.y) > 0.5) + { + wy += (wy < tau_frac.y) ? 1.0 : -1.0; + } + if (std::abs(wz - tau_frac.z) > 0.5) + { + wz += (wz < tau_frac.z) ? 1.0 : -1.0; + } + wsum += wgt; + cx += wgt * wx; + cy += wgt * wy; + cz += wgt * wz; + } + ASSERT_GT(wsum, 0.0); + const ModuleBase::Vector3 com2(cx / wsum, cy / wsum, cz / wsum); + const double d2a = (com2 - tau_frac) * (com2 - tau_frac); + const double d2m = (com2 - mirror_frac) * (com2 - mirror_frac); + EXPECT_LT(std::sqrt(d2a), 0.15); + EXPECT_GT(d2m, 100.0 * d2a); +} diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp index dfbc4f70923..aa8ed639268 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp @@ -2,6 +2,8 @@ #include "source_base/constants.h" +#include + void DFPTSerialBase::SetUp() { SetUpCell(); @@ -147,6 +149,52 @@ void DFPTSerialBase::MakeTwoAtomCell() ucell_.iat2ia[1] = 0; } +void DFPTSerialBase::MakeNlccTwoAtomCell() +{ + // start from the symmetric-broken two-atom cell, then upgrade the + // second species (the one at the non-origin tau_nlcc_) to an NC + // pseudo with a nonlinear core correction + MakeTwoAtomCell(); + Atom& at = ucell_.atoms[1]; + at.label = "Bc"; + at.coulomb_potential = false; + pseudo& p = at.ncpp; + p.zv = 2.0; + p.tvanp = false; + p.has_so = false; + p.nbeta = 0; + p.lll = {}; + p.nh = 0; + p.msh = 121; + p.kkbeta = 121; + p.r.resize(121); + p.rab.resize(121); + p.vloc_at.assign(121, 0.0); + const double dx = 0.025; + for (int i = 0; i < 121; ++i) + { + p.r[i] = i * dx; + p.rab[i] = dx; + } + // Gaussian core charge normalized to the integrated nlcc_core_charge_ + // (rho_atc carries the 4 pi r^2 measure of the radial grid) + const double s2 = nlcc_sigma_ * nlcc_sigma_; + double norm = 0.0; + for (int i = 0; i < 121; ++i) + { + const double r = p.r[i]; + norm += r * r * std::exp(-r * r / (2.0 * s2)) * p.rab[i]; + } + const double amp = nlcc_core_charge_ / (norm * ModuleBase::FOUR_PI); + p.rho_atc.resize(121); + for (int i = 0; i < 121; ++i) + { + const double r = p.r[i]; + p.rho_atc[i] = amp * std::exp(-r * r / (2.0 * s2)); + } + p.nlcc = true; +} + long long DFPTSerialBase::FKey(int ix, int iy, int iz) const { return (static_cast(ix + 64) * 128 + (iy + 64)) * 128 + (iz + 64); diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h index 77285d8c7ce..065ee223fd2 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h @@ -63,6 +63,21 @@ class DFPTSerialBase : public testing::Test void MakeCoulombAtom(); void MakeNCAtom(); + /// extend the two-atom cell of MakeTwoAtomCell with a nonlinear core + /// correction on the second species: a Gaussian rho_atc of total charge + /// nlcc_core_charge_ on the existing radial grid (msh points). The + /// atom sits at the non-origin tau_nlcc_ = (0.25, 0.31, 0.17) -- the + /// anti-blind-spot position: every structure-factor phase is active, + /// so a mirrored phase sign (exp(+iG.tau) instead of the QE addcore + /// exp(-iG.tau)) places the response at the wrong site and is caught + /// directly by the drhoc-localization test. + void MakeNlccTwoAtomCell(); + + /// total integrated core charge (in e) of the nlcc Gaussian above + static constexpr double nlcc_core_charge_ = 0.24; + const ModuleBase::Vector3 tau_nlcc_{0.25, 0.31, 0.17}; // lat0 units + double nlcc_sigma_ = 1.5; // bohr, grid-resolvable (spacing ~0.63 bohr) + // reconfigure the cell as a two-atom Z=4/Z=2 crystal breaking all symmetry void MakeTwoAtomCell(); From 2580552557cbaf3408b039e8bd02b664e9c10415 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Thu, 29 Oct 2026 21:35:03 +0800 Subject: [PATCH 18/21] fix(dfpt): C++11 aggregate: no default member initializer on InitContext The vxc_r field carried a '= nullptr' default, which under C++11 makes DFPT_PW_InitContext a non-aggregate; the 12-argument thin wrapper's brace-init then fails to compile in the non-MPI / non-LCAO CI configurations. Remove the default and pass nullptr explicitly as the 13th initializer in the wrapper. --- source/source_pw/module_dfpt/dfpt_pw.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 65d12dc9dd9..0b014cc8b22 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -46,7 +46,8 @@ struct DFPT_PW_InitContext /// ground-state XC potential on the rho grid (valence + core charge /// input, the same density the kernel linearizes around); consumed by /// the nlcc diagonal 2n+1 term (QE dynmatcc). Null disables that term. - const std::vector* vxc_r = nullptr; + /// No default member initializer: C++11 aggregate rules require it. + const std::vector* vxc_r; }; /** @@ -95,7 +96,7 @@ class DFPT_PW double ecutwfc, const Plus_U_Base* dftu) { - const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu}; + const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu, nullptr}; init(ctx); } From a2889fe6d9ad91d295b005e44c2df709f536f0a8 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 30 Oct 2026 09:39:28 +0800 Subject: [PATCH 19/21] fix(dfpt): address review findings on the nlcc PR P1 conjugation: the kernel contractions now conjugate the ROW response (sum(conj(row)*col)), matching the wavefunction and Ewald terms; the old unconjugated form flips the imaginary cross terms at nonzero q. The fill-semantics regression test's independent recomputation is updated to the same convention. P1 pool reduction: the kernel and diagonal terms now normalize by Omega/nxyz (the full grid, not the rank-local slab nrxx) and reduce their products over the pool with Parallel_Reduce::reduce_pool behind #ifdef __MPI guards (the mixing PR's convention). Non-blocking: remove the leftover DFPT_DUMP_DYN production block; extract the conditional setter arguments to named locals (AGENTS rule 14); tighten the mirror-distance assertions to minimum-image distances that reject the reflected profile at any periodic image. --- source/source_pw/module_dfpt/dfpt_phon.cpp | 27 -------------- .../source_pw/module_dfpt/dfpt_phon_elec.cpp | 28 +++++++++++---- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 4 ++- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 3 +- .../test_serial/dfpt_phon_serial_test.cpp | 2 +- .../test_serial/dfpt_rho_serial_test.cpp | 36 ++++++++++++++----- 6 files changed, 55 insertions(+), 45 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 1a414dbc562..4b04c7e24c2 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -127,33 +127,6 @@ void DFPT_Phon::assemble(int q_idx, DFPT_PW_Data& data) } } data.set_dynmat(q_idx, dyn); - // TEMPORARY QE-format dump for the nlcc element-level comparison - // (unreduced force constants; remove after the residual localization) - if (getenv("DFPT_DUMP_DYN") != nullptr) - { - std::ostringstream name; - name << "abacus_dyn.q" << q_idx << ".txt"; - std::ofstream out(name.str()); - out << std::setprecision(8) << std::scientific; - for (int ia = 0; ia < nat; ++ia) - { - for (int ja = 0; ja < nat; ++ja) - { - const double mm = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass - * ucell_->atoms[ucell_->iat2it[ja]].mass); - out << " " << ia + 1 << " " << ja + 1 << "\n"; - for (int i = 0; i < 3; ++i) - { - for (int j = 0; j < 3; ++j) - { - const std::complex v = dyn(3 * ia + i, 3 * ja + j) * mm; - out << " " << v.real() << " " << v.imag(); - } - out << "\n"; - } - } - } - } dynmat_accum_ = ModuleBase::ComplexMatrix(); accum_q_ = -1; ModuleBase::timer::end("DFPT_Phon", "assemble"); diff --git a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp index 184ff8d51e2..ce1a780d4fa 100644 --- a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp @@ -11,6 +11,7 @@ #include "dfpt_rho.h" #include "source_base/math_integral.h" #include "source_base/math_sphbes.h" +#include "source_base/parallel_reduce.h" #include #include @@ -235,7 +236,9 @@ void DFPT_Phon::accum_nlcc_kernel(int atom_idx, } const int nat = ucell_->nat; const ModuleBase::Vector3 q_cart = data.get_qvec(q_idx) * ucell_->G; - const double fac = ucell_->omega / static_cast(nrxx); + // pool-global normalization: nrxx is the rank-local FFT slab size, + // so the raw sum needs a pool reduction (the mixing PR's convention) + const double fac = ucell_->omega / static_cast(pw_rho_->nxyz); const int rowb = 3 * atom_idx + dir; const double mass_row = ucell_->atoms[ucell_->iat2it[atom_idx]].mass; // QE dynmat_nlcc semantics, pinned by the element-level comparison @@ -269,7 +272,11 @@ void DFPT_Phon::accum_nlcc_kernel(int atom_idx, const bool col_has_core = ucell_->atoms[ucell_->iat2it[iat2]].ncpp.nlcc; for (int idir2 = 0; idir2 < 3; ++idir2) { - // X_ab: the row's dvaux against the column's core response + // X_ab: the row's dvaux against the column's core response. + // Conjugation convention: the ROW response carries the + // conjugate (sum(conj(row)*col), matching the wavefunction + // and Ewald terms; at Gamma with real amplitudes this is a + // no-op, at nonzero q it fixes the imaginary cross terms) std::complex x_ab(0.0, 0.0); if (col_has_core) { @@ -281,12 +288,15 @@ void DFPT_Phon::accum_nlcc_kernel(int atom_idx, pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); for (int ir = 0; ir < nrxx; ++ir) { - x_ab += dvxc_r[ir] * std::conj(drhoc_r[ir]); + x_ab += std::conj(dvxc_r[ir]) * drhoc_r[ir]; } +#ifdef __MPI + Parallel_Reduce::reduce_pool(x_ab); +#endif x_ab *= fac; } // X_ba: the column's dvaux (per-q table) against the row's - // core response + // core response; same conjugation swap std::complex x_ba(0.0, 0.0); const int cola = 3 * iat2 + idir2; if (row_has_core && have_table && cola < static_cast(dvxc_core_->size()) @@ -295,8 +305,11 @@ void DFPT_Phon::accum_nlcc_kernel(int atom_idx, const std::vector>& dvxc_b = (*dvxc_core_)[cola]; for (int ir = 0; ir < nrxx; ++ir) { - x_ba += dvxc_b[ir] * std::conj(drhoc_a_r[ir]); + x_ba += std::conj(dvxc_b[ir]) * drhoc_a_r[ir]; } +#ifdef __MPI + Parallel_Reduce::reduce_pool(x_ba); +#endif x_ba *= fac; } if (x_ab == std::complex(0.0, 0.0) && x_ba == std::complex(0.0, 0.0)) @@ -346,7 +359,7 @@ void DFPT_Phon::accum_nlcc_diag(int atom_idx, int q_idx) std::vector> w_g(npw); std::vector> w_r(nrxx); double d3[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - const double fac = ucell_->omega / static_cast(nrxx); + const double fac = ucell_->omega / static_cast(pw_rho_->nxyz); for (int i = 0; i < 3; ++i) { for (int j = i; j < 3; ++j) @@ -381,6 +394,9 @@ void DFPT_Phon::accum_nlcc_diag(int atom_idx, int q_idx) { s += w_r[ir].real() * (*vxc_r_)[ir]; } +#ifdef __MPI + Parallel_Reduce::reduce_pool(s); +#endif d3[i][j] = fac * s; if (j != i) { diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 2606bb7fb29..c334fc8f795 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -117,7 +117,9 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, phon_.set_rho(&rho_); // the Impl-owned copy: the context pointer refers to caller stack // state that dies when the esolver init returns - phon_.set_vxc_r(vxc_r_.empty() ? nullptr : &vxc_r_); + // AGENTS rule 14: named local instead of a conditional in the call + const std::vector* vxc_provider = vxc_r_.empty() ? nullptr : &vxc_r_; + phon_.set_vxc_r(vxc_provider); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); data_.init(&qlist_, nk, nbands, npw_max, nrxx, nspin, nat, ctx.dftu); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 25c3784a268..b4bce1b02a1 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -167,7 +167,8 @@ void DFPT_PW::Impl::stash_converged_disp_response(int q_idx, // stash the converged density response too: the two-pass accumulation // rebuilds f_xc.(drho_val + drhoc) per displacement after the working // drho slots have moved on to the next displacement - data_.set_drho_disp(iat, idir, data_.get_drho_g(q_idx, 0)); + const std::vector> drho_conv = data_.get_drho_g(q_idx, 0); + data_.set_drho_disp(iat, idir, drho_conv); std::vector>>> disp( nk, std::vector>>(nbands)); 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 9ced75800fe..9ede7275d94 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 @@ -1039,7 +1039,7 @@ TEST_F(DFPTPhonSerialTest, NlccDirectTermsFillSemantics) std::complex x_ba(0.0, 0.0); for (int ir = 0; ir < nrxx; ++ir) { - x_ba += table[j][ir] * std::conj(drhoc_r[ir]); + x_ba += std::conj(table[j][ir]) * drhoc_r[ir]; } x_ba *= fac; const std::complex expect = 0.5 * std::conj(x_ba) / (m0 * m1); diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 7d8b180c770..5ccc8c574df 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -752,12 +752,22 @@ TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) } ASSERT_GT(wsum, 0.0); const ModuleBase::Vector3 com(cx / wsum, cy / wsum, cz / wsum); - const double dist_atom = (com - tau_frac) * (com - tau_frac); - const double dist_mirror = (com - mirror_frac) * (com - mirror_frac); + // minimum-image distance: wrap each component into [-0.5, 0.5) + auto mindist2 = [](const ModuleBase::Vector3& a, + const ModuleBase::Vector3& b) { + double dx = a.x - b.x; dx -= std::round(dx); + double dy = a.y - b.y; dy -= std::round(dy); + double dz = a.z - b.z; dz -= std::round(dz); + return dx * dx + dy * dy + dz * dz; + }; + const double dist_atom = mindist2(com, tau_frac); + const double dist_mirror = mindist2(com, mirror_frac); // the centroid of a derivative field sits at the atom up to the - // derivative skew (~sigma); the mirror is a full cell away - EXPECT_LT(std::sqrt(dist_atom), 0.1) << "dir=" << dir; - EXPECT_GT(dist_mirror, 100.0 * dist_atom) << "dir=" << dir; + // derivative skew; the mirror is a full half-cell away in the + // minimum-image sense, far beyond any tolerance that would admit + // the reflected profile + EXPECT_LT(dist_atom, 0.01) << "dir=" << dir; + EXPECT_GT(dist_mirror, 25.0 * dist_atom) << "dir=" << dir; } // non-zero q variant: the q-phase exp(-i q.tau) must not mirror either. @@ -802,8 +812,16 @@ TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) } ASSERT_GT(wsum, 0.0); const ModuleBase::Vector3 com2(cx / wsum, cy / wsum, cz / wsum); - const double d2a = (com2 - tau_frac) * (com2 - tau_frac); - const double d2m = (com2 - mirror_frac) * (com2 - mirror_frac); - EXPECT_LT(std::sqrt(d2a), 0.15); - EXPECT_GT(d2m, 100.0 * d2a); + // same minimum-image comparison as above + auto mindist2 = [](const ModuleBase::Vector3& a, + const ModuleBase::Vector3& b) { + double dx = a.x - b.x; dx -= std::round(dx); + double dy = a.y - b.y; dy -= std::round(dy); + double dz = a.z - b.z; dz -= std::round(dz); + return dx * dx + dy * dy + dz * dz; + }; + const double d2a = mindist2(com2, tau_frac); + const double d2m = mindist2(com2, mirror_frac); + EXPECT_LT(d2a, 0.02); + EXPECT_GT(d2m, 25.0 * d2a); } From 824be4fe2cd1d646a272f8a4a13d8f747683e4b1 Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Fri, 30 Oct 2026 19:07:39 +0800 Subject: [PATCH 20/21] fix(dfpt): address second-round review on the nlcc PR P1 MPI participation (bot): the vxc provider now always passes a valid pointer (never nullptr on an empty local slab), so every rank enters accum_nlcc_diag's collective FFT/reduction path; the size checks inside gate only the local work. Rule 14 (bot): the gtau construction extracts cos(phase) and -sin(phase) to named locals. Formula comment (bot): dfpt_phon.h's class doc now states the row-conjugated, pool-summed Omega/nxyz kernel convention, the 1/2-symmetrized fill, and the full-table f_xc.(drho_val + drhoc) contract. Phase-sensitive q test (bot): the generic-q drhoc check now also verifies the G-space coefficient phase against the expected exp(-i q_cart . tau_cart) via a dominant-G ratio between q and q=0, catching a sign flip that the |field|^2 centroid is blind to. Nonconstant diag oracle (bot): tracked as follow-up (the current fill-semantics test targets the kernel path; a dedicated diag test with an analytic vxc field is the next test addition). --- source/source_esolver/esolver_dfpt_pw.cpp | 2 +- source/source_pw/module_dfpt/dfpt_phon.h | 24 +++++++---- .../source_pw/module_dfpt/dfpt_phon_elec.cpp | 5 ++- source/source_pw/module_dfpt/dfpt_pw_init.cpp | 7 ++-- .../test_serial/dfpt_rho_serial_test.cpp | 41 +++++++++++++++++++ 5 files changed, 65 insertions(+), 14 deletions(-) diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 3c753354180..e887601e8a9 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -477,7 +477,7 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) nelec_, ecutwfc_, dft_plus_u_ ? this->dftu_.get() : nullptr, - vxc_r.empty() ? nullptr : &vxc_r}; + &vxc_r}; dfpt_->init(ctx); } diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index cffd3695d03..11e6a52099a 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -40,15 +40,21 @@ class DFPT_Rho; * for phase B). * * nlcc systems add two direct terms (QE ph.x ground truth: dynmat_nlcc - * and dynmatcc, both absent from the wavefunction response): - * D_ab += (Omega/nrxx) sum_r dvxc^a(r) conj(drhoc^b(r)) - * (dvxc^a = f_xc . (drho_val^a + drhoc^a) passed in by the caller; the - * mixed core-valence and core-core kernel pieces) and the q-independent - * diagonal correction - * D_aa[i][j] -= Omega tpiba^2 sum_G G_i G_j rhoc~(|G|) - * exp(-i G tau_a) conj(vxc^0(G)) - * (second derivative of the shifted core charge contracted with the - * ground-state XC potential). + * and dynmatcc, both absent from the wavefunction response). Kernel + * term (row-conjugated, pool-summed): + * X_ab = (Omega/nxyz) sum_pool conj(dvxc^a(r)) drhoc^b(r), + * D_ab += (1/2)(X_ab + conj(X_ba)) / sqrt(m_a m_b) + * where dvxc^a = f_xc . (drho_val^a + drhoc^a) is the first-order XC + * potential of the TOTAL response of displacement a, supplied by the + * caller through the per-q kernel table (covering ALL atoms: rows of + * non-nlcc atoms carry their valence-only dvxc; leaving their slots + * empty silently drops the one-sided cross blocks). Diagonal + * correction (q-independent, pool-summed): + * D_aa[i][j] += (Omega/nxyz) sum_pool Re[w_ij(r) vxc^0(r)] / m_a + * where w_ij = -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau_a), evaluated + * through recip2real as a real field. The ground-state XC potential + * vxc^0 is the valence + core charge density's XC potential, supplied + * through the InitContext by the esolver wiring. */ class DFPT_Phon { diff --git a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp index ce1a780d4fa..38b86915ef4 100644 --- a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp @@ -385,7 +385,10 @@ void DFPT_Phon::accum_nlcc_diag(int atom_idx, int q_idx) rhocg *= ModuleBase::FOUR_PI / ucell_->omega; const double phase = ModuleBase::TWO_PI * (gv * tau_cart) / ucell_->lat0; - const std::complex gtau(std::cos(phase), -std::sin(phase)); + // AGENTS rule 14: named locals instead of computed ctor args + const double gtau_re = std::cos(phase); + const double gtau_im = -std::sin(phase); + const std::complex gtau(gtau_re, gtau_im); w_g[ig] = -ucell_->tpiba * ucell_->tpiba * gv[i] * gv[j] * rhocg * gtau; } pw_rho_->recip2real(w_g.data(), w_r.data()); diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index c334fc8f795..179c001a665 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -117,9 +117,10 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, phon_.set_rho(&rho_); // the Impl-owned copy: the context pointer refers to caller stack // state that dies when the esolver init returns - // AGENTS rule 14: named local instead of a conditional in the call - const std::vector* vxc_provider = vxc_r_.empty() ? nullptr : &vxc_r_; - phon_.set_vxc_r(vxc_provider); + // always pass the pointer (never nullptr): a rank with an empty + // local field must still enter accum_nlcc_diag's collective + // participation; the size checks inside gate only the local work + phon_.set_vxc_r(&vxc_r_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); data_.init(&qlist_, nk, nbands, npw_max, nrxx, nspin, nat, ctx.dftu); diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 5ccc8c574df..cc0e1a1b8d4 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -824,4 +824,45 @@ TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) const double d2m = mindist2(com2, mirror_frac); EXPECT_LT(d2a, 0.02); EXPECT_GT(d2m, 25.0 * d2a); + + // phase-sensitive check at generic q: the centroid via |field|^2 is + // invariant under a global phase rotation, so flipping the sign of + // exp(i q.tau) alone would leave it unchanged. Instead, verify the + // G-space coefficient phase directly: drhoc_g(G) must carry + // exp(-i(G+q).tau) (the QE addcore convention), not the conjugate. + // Compare the q=0 and q=q2 coefficients: their ratio at each G must + // be exp(-i q.cart . tau) * (ratio of the radial/Bessel factor at + // |G| vs |G+q|). For a single G with a dominant coefficient this + // isolates the structure-factor phase. + const ModuleBase::Vector3 q_zero_cart(0.0, 0.0, 0.0); + std::vector> drhoc_g0; + rho.compute_drhoc_g(ucell_, iat, 1, q_zero_cart, drhoc_g0); + ASSERT_EQ(static_cast(drhoc_g0.size()), pw_rho_.npw); + // find the G with the largest |G|=0-adjacent coefficient (where the + // radial factor is comparable) and verify the phase relation + int ig_best = -1; + double best_mag = 0.0; + for (int ig = 0; ig < pw_rho_.npw; ++ig) + { + const double m = std::norm(drhoc_g0[ig]); + if (m > best_mag) + { + best_mag = m; + ig_best = ig; + } + } + ASSERT_GE(ig_best, 0); + // the expected extra phase from q is exp(-i q_cart . tau_cart): + // drhoc_q(G) / drhoc_0(G) ~= [radial(|G+q|)/radial(|G|)] * exp(-i q.tau) + // for the dominant G (where |G+q| ~= |G| for small q), the ratio's + // phase is approximately -q_cart . tau_cart + const std::complex ratio = drhoc_g[ig_best] / drhoc_g0[ig_best]; + const double expected_phase_arg = -(q2_cart * (ucell_.atoms[it].tau[ucell_.iat2ia[iat]] * ucell_.lat0)); + const std::complex expected_phase(std::cos(expected_phase_arg), std::sin(expected_phase_arg)); + // the phase of the ratio should match the expected exp(-i q.tau) + // (within the radial-factor ratio, which is positive real) + const double phase_diff = std::arg(ratio / expected_phase); + // wrap to [-pi, pi] and allow the radial-ratio approximation + const double phase_diff_wrapped = std::atan2(std::sin(phase_diff), std::cos(phase_diff)); + EXPECT_NEAR(phase_diff_wrapped, 0.0, 0.5) << "phase of drhoc_q/drhoc_0 ratio vs expected exp(-i q.tau)"; } From 03f0d2d7bf63a354ffc6dd636ee9492311e2a32a Mon Sep 17 00:00:00 2001 From: Zanthoxylum <105619360+Zanthoxylum@users.noreply.github.com> Date: Sat, 31 Oct 2026 00:56:45 +0800 Subject: [PATCH 21/21] fix(dfpt): address third-round review on the nlcc PR - compute_drhoc: early return for invalid/coreless atoms before the diagnostic iat2it indexing (public API safety) - drho_disp stash: gated on nlcc_any, so coreless runs skip the 3*nat*npw complex copies - phase test: full TWO_PI-scaled phase reference with a tight 0.05-rad tolerance (both omission and reversal of exp(-i q.tau) are now rejected); ratio and trig arguments extracted to named locals (rule 14) - set_dvxc_core_table comment: updated to describe the total-response f_xc.(drho_val + drhoc) table (not core-only) --- source/source_pw/module_dfpt/dfpt_phon.h | 5 ++- .../source_pw/module_dfpt/dfpt_pw_solve.cpp | 20 +++++++++-- source/source_pw/module_dfpt/dfpt_rho.cpp | 8 +++++ .../test_serial/dfpt_rho_serial_test.cpp | 34 +++++++++++++------ 4 files changed, 52 insertions(+), 15 deletions(-) diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 11e6a52099a..7e460315ac1 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -71,7 +71,10 @@ class DFPT_Phon void set_rho(const DFPT_Rho* rho); void set_vxc_r(const std::vector* vxc_r); - /// per-q table of f_xc.drhoc for every displacement slot [3*nat] + /// per-q table of dvxc^a = f_xc . (drho_val^a + drhoc^a) for every + /// displacement slot [3*nat]: the first-order XC potential of the + /// TOTAL response (valence + core) supplied by the solver's + /// accumulation loop; coreless atoms carry their valence-only dvxc /// (the kernel term contracts it against the stashed valence /// responses of the mixed piece). Set by the solver once per q, /// before the accumulation pass; owned by the caller. diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index b4bce1b02a1..1a30c81fe88 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -166,9 +166,23 @@ void DFPT_PW::Impl::stash_converged_disp_response(int q_idx, data_.set_vsc_r(iat, idir, v_sc_r_last); // stash the converged density response too: the two-pass accumulation // rebuilds f_xc.(drho_val + drhoc) per displacement after the working - // drho slots have moved on to the next displacement - const std::vector> drho_conv = data_.get_drho_g(q_idx, 0); - data_.set_drho_disp(iat, idir, drho_conv); + // drho slots have moved on; only needed when some species carries + // a nonlinear core correction (skip the 3*nat*npw complex copies + // for coreless runs) + bool nlcc_any_solve = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any_solve = true; + break; + } + } + if (nlcc_any_solve) + { + const std::vector> drho_conv = data_.get_drho_g(q_idx, 0); + data_.set_drho_disp(iat, idir, drho_conv); + } std::vector>>> disp( nk, std::vector>>(nbands)); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 00d6cfe09d9..5dd8abe3349 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -295,6 +295,14 @@ void DFPT_Rho::compute_drhoc(const UnitCell& ucell, compute_drhoc_g(ucell, iat, dir, q_cart, drhoc_g); data.set_drhoc_g(q_idx, drhoc_g); + // early return for invalid/coreless atoms: compute_drhoc_g already + // returned empty, so the iat2it indexing below would be out of + // bounds for a bad index — skip the diagnostic in that case + if (drhoc_g.empty() || iat < 0 || iat >= ucell.nat) + { + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); + return; + } const int it = ucell.iat2it[iat]; // report magnitude for diagnostics double norm2 = 0.0; diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index cc0e1a1b8d4..3e7b268d4c7 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -853,16 +853,28 @@ TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) } ASSERT_GE(ig_best, 0); // the expected extra phase from q is exp(-i q_cart . tau_cart): - // drhoc_q(G) / drhoc_0(G) ~= [radial(|G+q|)/radial(|G|)] * exp(-i q.tau) - // for the dominant G (where |G+q| ~= |G| for small q), the ratio's - // phase is approximately -q_cart . tau_cart + // drhoc_q(G) = -i|G+q| tpiba rho~(|G+q|) exp(-i(G+q).tau) + // drhoc_0(G) = -i|G| tpiba rho~(|G|) exp(-i G. tau) + // ratio = [|G+q|/|G| · rho~(|G+q|)/rho~(|G|)] · exp(-i q. tau) + // The bracket is positive real; the phase is exactly -q_cart . tau_cart. + // Use the full TWO_PI-scaled phase (the tau and q_cart are already in + // consistent units: tau in lat0-cartesian, q_cart in 2pi/lat0 units, + // so the dot product directly gives the phase in radians). const std::complex ratio = drhoc_g[ig_best] / drhoc_g0[ig_best]; - const double expected_phase_arg = -(q2_cart * (ucell_.atoms[it].tau[ucell_.iat2ia[iat]] * ucell_.lat0)); - const std::complex expected_phase(std::cos(expected_phase_arg), std::sin(expected_phase_arg)); - // the phase of the ratio should match the expected exp(-i q.tau) - // (within the radial-factor ratio, which is positive real) - const double phase_diff = std::arg(ratio / expected_phase); - // wrap to [-pi, pi] and allow the radial-ratio approximation - const double phase_diff_wrapped = std::atan2(std::sin(phase_diff), std::cos(phase_diff)); - EXPECT_NEAR(phase_diff_wrapped, 0.0, 0.5) << "phase of drhoc_q/drhoc_0 ratio vs expected exp(-i q.tau)"; + const ModuleBase::Vector3 tau_cart_alat = ucell_.atoms[it].tau[ucell_.iat2ia[iat]]; + const double expected_phase_arg = -ModuleBase::TWO_PI * (q2_cart * tau_cart_alat) / ucell_.lat0; + const double cos_arg = std::cos(expected_phase_arg); + const double sin_arg = std::sin(expected_phase_arg); + const std::complex expected_phase(cos_arg, sin_arg); + // tight tolerance: the exact phase relation is exp(-i q.tau), so both + // omission (arg=0) and reversal (arg=+q.tau) must be rejected; the + // only residual is the positive-real radial-ratio bracket (which has + // zero phase contribution) and floating-point noise + // rule 14: named locals for computed function arguments + const std::complex unit_phase_ratio = ratio / expected_phase; + const double phase_diff = std::arg(unit_phase_ratio); + const double sin_diff = std::sin(phase_diff); + const double cos_diff = std::cos(phase_diff); + const double phase_diff_wrapped = std::atan2(sin_diff, cos_diff); + EXPECT_NEAR(phase_diff_wrapped, 0.0, 0.05) << "phase of drhoc_q/drhoc_0 ratio vs expected exp(-i q.tau)"; }