feat(dfpt): v_sc Broyden mixing for the response SCF (default), plus rho_broyden, rho_kerker, v_sc_plain

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);
 }

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

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

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);

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);

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<Base_Mixing::Broyden_Mixing> vsc_broyden_;
+    std::unique_ptr<Base_Mixing::Mixing_Data> vsc_broyden_mdata_;
 };
 
 } // namespace ModuleDFPT

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_));

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 <cmath>
 #include <complex>
 #include <cstdlib>
+#include <iomanip>
 #include <iostream>
+#include <sstream>
 #include <vector>
 
 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<std::complex<double>>(npw, std::complex<double>(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<char>('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");