diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 777a678d4a..c22017f783 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -225,6 +225,9 @@ 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_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 1c624aa5a0..8bf7d91c91 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -440,6 +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 = "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 cfea717beb..a673bf2548 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -113,13 +113,20 @@ 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, 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) { @@ -147,12 +154,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) { @@ -199,5 +208,88 @@ void ReadInput::item_dfpt() }; 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; + item.type = "String"; + 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. 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 = "v_sc_broyden"; + } + else + { + para.input.dfpt_mix_type = strvalue; + } + }; + sync_string(input.dfpt_mix_type); + item.check_value = [](const Input_Item& item, const Parameter& para) { + 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 rho_plain, rho_kerker, rho_broyden, v_sc_plain or v_sc_broyden."); + } + }; + this->add_item(item); + } + { + Input_Item item("dfpt_kerker_a2"); + 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 = 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 other schemes."; + 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); + } + { + Input_Item item("dfpt_broyden_ndim"); + 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 = 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) { + 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 de3b7a955a..4a00e89b80 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -146,6 +146,31 @@ void DFPT_PW::set_asr(const std::string& 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" + || mix_type == "v_sc_plain" || mix_type == "v_sc_broyden") + { + pimpl_->mix_type_ = mix_type; + } +} + +void DFPT_PW::set_kerker_a2(double a2) +{ + if (a2 > 0.0) + { + pimpl_->kerker_a2_ = 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 ea611743b9..69ec3892d1 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -131,6 +131,17 @@ 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); + /// 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 (rho_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 ed32d6ad79..f7ec543731 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" @@ -178,6 +180,13 @@ 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_ = "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_; + 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 2b702514f6..47240b6618 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -87,38 +87,30 @@ 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")) + // |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 / + // 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_; + // 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") { - 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; - } + rho_mix_type = "rho_plain"; } - double kerker_a2 = 1.0; - if (const char* env_a2 = getenv("DFPT_KERKER_A2")) + else if (mix_type == "v_sc_broyden") { - const double parsed = atof(env_a2); - if (parsed > 0.0) - { - kerker_a2 = parsed; - } + rho_mix_type = "rho_broyden"; } - 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, 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_q0.cpp b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp index 39dba1d3f8..16758a7231 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 d467aa3622..876519826f 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -11,19 +11,43 @@ #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" #include #include #include #include +#include +#include #include #include 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 @@ -229,16 +253,50 @@ 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. 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 + 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; - 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) { @@ -270,8 +328,98 @@ 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); + 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(); + } +#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; 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) + { + dn2 += std::norm(v_sc_out_r[ir] - v_sc_in_r[ir]); + o2 += std::norm(v_sc_out_r[ir]); + } +#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; 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) + { + 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; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(dn2); + Parallel_Reduce::reduce_pool(o2); +#endif + residual = safe_relative_defect(dn2, o2); + } + // 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 + { + 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 @@ -302,6 +450,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"); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 239c7804ab..326c4c7f04 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) @@ -37,14 +60,27 @@ 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") + // 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") { - ModuleBase::WARNING_QUIT("DFPT_Rho", "unsupported mix_type, expected plain or kerker"); + 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"); } @@ -386,6 +422,21 @@ 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(); + } + // 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"); } @@ -422,7 +473,64 @@ 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_ == "rho_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*)>(); + // 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(), + screen, + combine, + true); + // 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) { + double s = 0.0; + for (int i = 0; i < npw; ++i) + { + 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_ == "rho_kerker") { const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; std::vector> rin_s(npw); @@ -455,7 +563,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) @@ -463,7 +572,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 83ae0cda25..89cb83a95a 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -15,6 +15,8 @@ namespace Base_Mixing { class Plain_Mixing; +class Broyden_Mixing; +class Mixing_Data; } namespace ModuleDFPT @@ -76,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 @@ -100,9 +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" or "kerker" + 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 (kerker only) + double kerker_a2; // Kerker screen a^2, 1/lat0^2 (rho_kerker only) + int broyden_ndim; // Broyden history depth (rho_broyden only, must be >= 1; no C++11 aggregate default) }; DFPT_Rho(); @@ -186,12 +189,18 @@ 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; 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_; 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 4e4a328f14..66c0543f4a 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) { @@ -529,3 +529,134 @@ TEST_F(DFPTRhoSerialTest, MixDrhoKerkerStabilizesStiffModelProblem) EXPECT_NEAR(final_in[ig].imag(), target[ig].imag(), 1.0e-10); } } + +TEST_F(DFPTRhoSerialTest, BroydenConvergesDiagonalSystemFast) +{ + // 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 populated shell by |1 - beta (1 - D)| = 0.8 (d = 0.5, + // 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; + 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> 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; + 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_); + ++n_iter; + const double res = rho.get_residual(0, data_); + if (res < 1.0e-10) + { + break; + } + 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 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); + // 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); + } +}