feat(dfpt): analytic libxc XC kernel with the i(G+q) carrier term (exact GGA response at q != 0)

source/source_esolver/esolver_dfpt_pw.cpp

@@ -13,6 +13,7 @@
 #include "source_base/macros.h"
 #include "source_base/math_polyint.h"
 #include "source_base/tool_quit.h"
+#include "source_hamilt/module_xc/libxc_abacus.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_pw/module_dfpt/dfpt_pert.h"
@@ -335,10 +336,24 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
     {
         rho0_r[ir] = chg0->rho[0][ir] + (has_core ? chg0->rho_core[ir] : 0.0);
     }
-    xc_adapter_ = new ModuleDFPT::XC_First_Order_Analytic(XC_Functional::get_func_id(),
+    // the kernel must differentiate the SAME effective functional the
+    // ground state used: external/in-built libxc parameters (finite-T
+    // LDA, tuned fractions) and the signed per-component combination
+    // factors (e.g. BLYP_LR subtracts its short-range components)
+    const std::vector<int> xc_func_id = XC_Functional::get_func_id();
+    const std::vector<std::vector<double>> xc_ext_params
+        = XC_Functional_Libxc::ext_params_for(xc_func_id, XC_Functional::get_hybrid_alpha(), XC_Functional::get_hse_omega());
+    std::vector<double> xc_weights(xc_func_id.size(), 1.0);
+    for (size_t i = 0; i < xc_func_id.size(); ++i)
+    {
+        xc_weights[i] = XC_Functional::get_weight(xc_func_id[i]);
+    }
+    xc_adapter_ = new ModuleDFPT::XC_First_Order_Analytic(xc_func_id,
                                                           XC_Functional::get_func_type(),
                                                           rho0_r,
-                                                          this->pw_rho);
+                                                          this->pw_rho,
+                                                          xc_ext_params,
+                                                          xc_weights);
 #else
     ModuleBase::WARNING("ESolver_DFPT_PW::init_dfpt",
                         "built without libxc: the DFPT XC kernel falls back to the finite-difference "

source/source_hamilt/module_xc/libxc_abacus.h

@@ -69,6 +69,22 @@ namespace XC_Functional_Libxc
         const double hybrid_alpha,
         const double hse_omega);
 
+    /**
+     * @brief Per-component external parameters that init_func would apply.
+     *
+     * Returns, for every functional id in func_id, the parameter list the
+     * ground-state init_func would pass to xc_func_set_ext_params (external
+     * xc_exch_ext/xc_corr_ext parameters take precedence over the in-built
+     * ones; empty vector = libxc defaults). Consumers that instantiate
+     * libxc handles themselves -- e.g. the DFPT first-order kernel, which
+     * needs the SAME effective functional as the ground state -- use this
+     * to configure their own handles without duplicating the parameter
+     * tables.
+     */
+    extern std::vector<std::vector<double>> ext_params_for(const std::vector<int>& func_id,
+                                                           const double hybrid_alpha,
+                                                           const double hse_omega);
+
     extern void finish_func(std::vector<xc_func_type> &funcs);
 
 //-------------------

source/source_hamilt/module_xc/libxc_setup.cpp

@@ -333,4 +333,23 @@ void XC_Functional_Libxc::finish_func(std::vector<xc_func_type> &funcs)
     }
 }
 
+std::vector<std::vector<double>> XC_Functional_Libxc::ext_params_for(const std::vector<int>& func_id,
+                                                                     const double hybrid_alpha,
+                                                                     const double hse_omega)
+{
+    std::vector<std::vector<double>> params;
+    params.reserve(func_id.size());
+    for (int id : func_id)
+    {
+        const std::vector<double> external_ext_params = external_xc_func_ext_params(id);
+        if (!external_ext_params.empty())
+        {
+            params.push_back(external_ext_params);
+            continue;
+        }
+        params.push_back(in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega));
+    }
+    return params;
+}
+
 #endif

source/source_hamilt/module_xc/xc_functional.h

@@ -193,6 +193,15 @@ class XC_Functional
     public:
     static std::vector<int> get_func_id() { return func_id; }
 
+    /// signed combination weight of a functional component id (1.0 when
+    /// the functional is not in the scaling table), e.g. BLYP_LR
+    /// subtracts its short-range components
+    static double get_weight(const int id)
+    {
+        const auto it = scaling_factor_xc.find(id);
+        return it != scaling_factor_xc.end() ? it->second : 1.0;
+    }
+
 //-------------------
 //  xc_lda_wrap.cpp
 //-------------------

source/source_pw/module_dfpt/dfpt_pw_solve.cpp

@@ -257,11 +257,17 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir)
         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
-        std::cout << "  DFPT mix iter=" << std::setw(3) << std::right << iter + 1 << "/" << max_iter_
-                  << " residual=" << std::setprecision(3) << std::scientific << residual
-                  << " (thr " << std::setprecision(1) << std::scientific << conv_thr_ << ")"
-                  << std::endl;
+        // above stay opt-in diagnostics. Formatting goes through a local
+        // buffer so the stream flags/precision are never disturbed for
+        // the surrounding output (#8164 review: setprecision pollution
+        // degraded the DBG snapshots of the next iteration)
+        {
+            std::ostringstream trace;
+            trace << "  DFPT mix iter=" << std::setw(3) << std::right << iter + 1 << "/" << max_iter_ << " residual="
+                  << std::setprecision(3) << std::scientific << residual << " (thr " << std::setprecision(1)
+                  << std::scientific << conv_thr_ << ")" << std::endl;
+            std::cout << trace.str();
+        }
         if (dbg)
         {
             std::cout << "DBG iter=" << iter << " residual=" << residual << " conv_thr=" << conv_thr_ << std::endl;
@@ -270,13 +276,16 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir)
         ++n_iter_done;
     }
     // convergence summary in the efield-solve style; names the state
-    // explicitly instead of failing silently at max_iter
-    std::cout << std::defaultfloat << std::setprecision(4);
-    std::cout << "DFPT disp q=(" << q_frac.x << "," << q_frac.y << "," << q_frac.z << ") atom=" << iat + 1 << "("
-              << ucell_->atoms[ucell_->iat2it[iat]].label << ") dir=" << static_cast<char>('x' + idir) << " "
-              << (converged ? "converged" : "NOT converged (max_iter reached)") << ", residual=" << std::scientific
-              << std::setprecision(6) << residual << " (iter=" << n_iter_done << ")" << std::defaultfloat
-              << std::endl;
+    // explicitly instead of failing silently at max_iter; formatted in a
+    // local buffer for the same stream-state reason
+    {
+        std::ostringstream summary;
+        summary << "DFPT disp q=(" << q_frac.x << "," << q_frac.y << "," << q_frac.z << ") atom=" << iat + 1 << "("
+                << ucell_->atoms[ucell_->iat2it[iat]].label << ") dir=" << static_cast<char>('x' + idir) << " "
+                << (converged ? "converged" : "NOT converged (max_iter reached)") << ", residual=" << std::scientific
+                << std::setprecision(6) << residual << " (iter=" << n_iter_done << ")" << std::endl;
+        std::cout << summary.str();
+    }
     stash_converged_disp_response(q_idx, iat, idir, v_sc_r_last, nk, nbands);
 
     ModuleBase::timer::end("DFPT_PW", "solve_displacement");

source/source_pw/module_dfpt/dfpt_xc_analytic.cpp

@@ -44,6 +44,27 @@ XC_First_Order_Analytic::XC_First_Order_Analytic(const std::vector<int>& func_id
     this->init_kernel_fields_(rho0_r);
 }
 
+XC_First_Order_Analytic::XC_First_Order_Analytic(const std::vector<int>& func_id,
+                                                 const int func_type,
+                                                 const std::vector<double>& rho0_r,
+                                                 ModulePW::PW_Basis* pw_rho,
+                                                 const std::vector<std::vector<double>>& ext_params_per_func,
+                                                 const std::vector<double>& weights)
+    : XC_First_Order_Analytic(func_id, func_type, rho0_r, pw_rho)
+{
+    // delegating ctor already built the default kernel fields; the
+    // configured factors require rebuilding them with the parameters
+    // applied, so validate first and reconstruct
+    if (ext_params_per_func.size() != func_id.size() || weights.size() != func_id.size())
+    {
+        ModuleBase::WARNING_QUIT("XC_First_Order_Analytic",
+                                 "ext_params_per_func / weights must match func_id elementwise");
+    }
+    ext_params_ = ext_params_per_func;
+    weights_ = weights;
+    this->init_kernel_fields_(rho0_r);
+}
+
 void XC_First_Order_Analytic::set_context(const ModuleBase::Vector3<double>& q_cart) const
 {
     q_cart_ = q_cart;
@@ -105,7 +126,29 @@ void XC_First_Order_Analytic::init_kernel_fields_(const std::vector<double>& rho
         {
             ModuleBase::WARNING_QUIT("XC_First_Order_Analytic", "xc_func_init failed");
         }
-        if (func_type_ == 1)
+        // apply the ground-state effective parameters (external override
+        // in-built; empty = libxc defaults) so the kernel differentiates
+        // the SAME functional the GS run used
+        if (ifunc < ext_params_.size() && !ext_params_[ifunc].empty())
+        {
+            const xc_func_info_type* info = xc_func_get_info(&func);
+            const int nref = xc_func_info_get_n_ext_params(info);
+            if (static_cast<int>(ext_params_[ifunc].size()) != nref)
+            {
+                ModuleBase::WARNING_QUIT("XC_First_Order_Analytic",
+                                         "invalid number of external parameters for libxc functional id "
+                                             + std::to_string(func_id_[ifunc]));
+            }
+            xc_func_set_ext_params(&func, ext_params_[ifunc].data());
+        }
+        // signed per-component combination factor (e.g. BLYP_LR subtracts
+        // its short-range components)
+        const double w = ifunc < weights_.size() ? weights_[ifunc] : 1.0;
+        // dispatch per COMPONENT family: a functional combination may mix
+        // LDA and GGA members (func_type_ only reflects the overall type),
+        // and the LDA entry points must be used for LDA-family components
+        const int family = xc_func_info_get_family(xc_func_get_info(&func));
+        if (family == XC_FAMILY_LDA)
         {
             for (int ir = 0; ir < nrxx; ++ir)
             {
@@ -114,7 +157,7 @@ void XC_First_Order_Analytic::init_kernel_fields_(const std::vector<double>& rho
                 xc_lda_fxc(&func, 1, &rho, &v2rho2);
                 // libxc vrho is already the potential d(rho e)/d rho, so
                 // the kernel is its plain derivative (QE dmuxc convention)
-                k_rr_[ir] += ha_to_ry * v2rho2;
+                k_rr_[ir] += ha_to_ry * w * v2rho2;
             }
         }
         else
@@ -123,9 +166,9 @@ void XC_First_Order_Analytic::init_kernel_fields_(const std::vector<double>& rho
             {
                 const double rho = rho0_r[ir];
                 const double sig = sigma[ir];
-                if (std::abs(rho) <= epsr || sig <= epsg)
+                if (std::abs(rho) <= epsr)
                 {
-                    continue; // QE dgradcorr guard: no kernel response here
+                    continue; // no XC response where the density vanishes
                 }
                 double zk = 0.0;
                 double vrho = 0.0;
@@ -135,16 +178,24 @@ void XC_First_Order_Analytic::init_kernel_fields_(const std::vector<double>& rho
                 double v2rhosigma = 0.0;
                 double v2sigma2 = 0.0;
                 xc_gga_fxc(&func, 1, &rho, &sig, &v2rho2, &v2rhosigma, &v2sigma2);
+                // the local (rho-rho) response is finite in the zero-gradient
+                // limit (the LDA-like term of the GGA functional), so only
+                // the density guard applies; QE retains it through the
+                // separately accumulated local term of dv_of_drho
+                k_rr_[ir] += ha_to_ry * w * v2rho2;
+                if (sig <= epsg)
+                {
+                    continue; // gradient kernel vanishes in the flat limit
+                }
                 // sigma = |grad rho0|^2 chain factors of the dgradcorr
                 // algebra: s1 = grho.gdrho = d sigma / 2, so the response
                 // e_rho_sigma d(sigma) carries the factor 2 and the
                 // e_sigma_sigma term the factor 4 (QE xc_wrapper_d_gga
                 // convention; no rho renormalization, the potential-form
                 // kernel cancels it)
-                k_rr_[ir] += ha_to_ry * v2rho2;
-                k_sr_[ir] += ha_to_ry * (2.0 * v2rhosigma);
-                k_ss_[ir] += ha_to_ry * (4.0 * v2sigma2);
-                k_s_[ir] += ha_to_ry * (2.0 * vsigma);
+                k_sr_[ir] += ha_to_ry * w * (2.0 * v2rhosigma);
+                k_ss_[ir] += ha_to_ry * w * (4.0 * v2sigma2);
+                k_s_[ir] += ha_to_ry * w * (2.0 * vsigma);
             }
         }
         xc_func_end(&func);

source/source_pw/module_dfpt/dfpt_xc_analytic.h

@@ -43,6 +43,19 @@ class XC_First_Order_Analytic : public XC_First_Order
                             const std::vector<double>& rho0_r,
                             ModulePW::PW_Basis* pw_rho);
 
+    /// Configured instantiation: ext_params_per_func carries the same
+    /// external parameter lists the ground-state init_func would apply
+    /// (XC_Functional_Libxc::ext_params_for; empty entry = libxc
+    /// defaults), weights the signed per-component combination factors
+    /// (e.g. BLYP_LR subtracts its short-range components; empty or
+    /// 1.0 = plain addition). Both must match func_id elementwise.
+    XC_First_Order_Analytic(const std::vector<int>& func_id,
+                            const int func_type,
+                            const std::vector<double>& rho0_r,
+                            ModulePW::PW_Basis* pw_rho,
+                            const std::vector<std::vector<double>>& ext_params_per_func,
+                            const std::vector<double>& weights);
+
     void set_context(const ModuleBase::Vector3<double>& q_cart) const override;
 
     void apply(const std::vector<std::complex<double>>& drho_r,
@@ -57,6 +70,8 @@ class XC_First_Order_Analytic : public XC_First_Order
 
     int func_type_ = 1; ///< 1 = LDA, 2 = GGA (XC_Functional convention)
     std::vector<int> func_id_; ///< libxc ids, exchange first (XC_Functional order)
+    std::vector<std::vector<double>> ext_params_; ///< per-id set_ext_params lists
+    std::vector<double> weights_; ///< signed per-id combination factors
     ModulePW::PW_Basis* pw_rho_ = nullptr;
 
     mutable ModuleBase::Vector3<double> q_cart_; ///< phonon q, ucell::G units (1/lat0)