diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 79637856b4..9e5eb6db3e 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/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 xc_func_id = XC_Functional::get_func_id(); + const std::vector> xc_ext_params + = XC_Functional_Libxc::ext_params_for(xc_func_id, XC_Functional::get_hybrid_alpha(), XC_Functional::get_hse_omega()); + std::vector 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 " diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index b4541c8a9e..a6846c3c0b 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/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> ext_params_for(const std::vector& func_id, + const double hybrid_alpha, + const double hse_omega); + extern void finish_func(std::vector &funcs); //------------------- diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 9f91c9276e..498116f1c3 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -333,4 +333,23 @@ void XC_Functional_Libxc::finish_func(std::vector &funcs) } } +std::vector> XC_Functional_Libxc::ext_params_for(const std::vector& func_id, + const double hybrid_alpha, + const double hse_omega) +{ + std::vector> params; + params.reserve(func_id.size()); + for (int id : func_id) + { + const std::vector 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 diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 3465a46794..a8b57f0046 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -193,6 +193,15 @@ class XC_Functional public: static std::vector 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 //------------------- diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index a5f3590509..878ad251cd 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/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('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('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"); diff --git a/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp b/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp index b485515e77..3b83d58ad7 100644 --- a/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp +++ b/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp @@ -44,6 +44,27 @@ XC_First_Order_Analytic::XC_First_Order_Analytic(const std::vector& func_id this->init_kernel_fields_(rho0_r); } +XC_First_Order_Analytic::XC_First_Order_Analytic(const std::vector& func_id, + const int func_type, + const std::vector& rho0_r, + ModulePW::PW_Basis* pw_rho, + const std::vector>& ext_params_per_func, + const std::vector& 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& q_cart) const { q_cart_ = q_cart; @@ -105,7 +126,29 @@ void XC_First_Order_Analytic::init_kernel_fields_(const std::vector& 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(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& 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& 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& 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); diff --git a/source/source_pw/module_dfpt/dfpt_xc_analytic.h b/source/source_pw/module_dfpt/dfpt_xc_analytic.h index a2d9d4bb7a..c8b5170544 100644 --- a/source/source_pw/module_dfpt/dfpt_xc_analytic.h +++ b/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& 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& func_id, + const int func_type, + const std::vector& rho0_r, + ModulePW::PW_Basis* pw_rho, + const std::vector>& ext_params_per_func, + const std::vector& weights); + void set_context(const ModuleBase::Vector3& q_cart) const override; void apply(const std::vector>& 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 func_id_; ///< libxc ids, exchange first (XC_Functional order) + std::vector> ext_params_; ///< per-id set_ext_params lists + std::vector weights_; ///< signed per-id combination factors ModulePW::PW_Basis* pw_rho_ = nullptr; mutable ModuleBase::Vector3 q_cart_; ///< phonon q, ucell::G units (1/lat0) diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_xc_analytic_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_xc_analytic_serial_test.cpp index 40d5d18918..265eba9aff 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_xc_analytic_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_xc_analytic_serial_test.cpp @@ -380,4 +380,165 @@ TEST(DFPTXcAnalytic, CarrierTermAppearsAtFiniteQOnlyForGga) EXPECT_EQ(ldiff, 0.0); } +TEST(DFPTXcAnalytic, UniformDensityKeepsLocalGgaResponse) +{ + // regression guard of the flat-gradient limit: with sigma = 0 everywhere + // the GGA response must reduce to the local (rho-rho) kernel v2rho2 -- + // finite and NONZERO -- not vanish under the gradient-cutoff guard + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0(pw.nrxx, 0.1); + const std::vector amp = cosine_amplitude(pw, 0.02); + const std::vector func_id = {XC_GGA_X_PBE, XC_GGA_C_PBE}; + ModuleDFPT::XC_First_Order_Analytic kernel(func_id, 2, rho0, &pw); + const ModuleBase::Vector3 q0(0.0, 0.0, 0.0); + kernel.set_context(q0); + + std::vector> a(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + a[ir] = std::complex(amp[ir], 0.3 * amp[ir]); + } + std::vector> dvxc; + kernel.apply(a, dvxc); + + // expected local kernel: sum over functional components, Ha -> Ry (x2) + double k_rr = 0.0; + for (size_t ifunc = 0; ifunc < func_id.size(); ++ifunc) + { + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, func_id[ifunc], XC_UNPOLARIZED), 0); + double rho = 0.1; + double sig = 0.0; + double v2rho2 = 0.0; + double v2rhosigma = 0.0; + double v2sigma2 = 0.0; + xc_gga_fxc(&func, 1, &rho, &sig, &v2rho2, &v2rhosigma, &v2sigma2); + k_rr += 2.0 * v2rho2; + xc_func_end(&func); + } + ASSERT_GT(std::abs(k_rr), 1.0e-3); // a vanishing local kernel is the bug + + double max_err = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + const std::complex expect = k_rr * a[ir]; + max_err = std::max(max_err, std::abs(dvxc[ir] - expect)); + } + EXPECT_LT(max_err / std::abs(k_rr), 1.0e-10); +} + +TEST(DFPTXcAnalytic, ConfiguredExtParamsDifferentiateTheConfiguredFunctional) +{ + // nondefault external parameters must reach the kernel: a finite-T + // LDA (KSDT) differentiated at theta = 0.007 must equal the manually + // parameterized libxc second derivative, not the default-parameter one + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0(pw.nrxx, 0.1); + const std::vector func_id = {XC_LDA_XC_KSDT}; + const std::vector> ext_params = {{0.007}}; + const std::vector weights = {1.0}; + ModuleDFPT::XC_First_Order_Analytic kernel(func_id, 1, rho0, &pw, ext_params, weights); + + std::vector> a(pw.nrxx, std::complex(0.05, 0.0)); + std::vector> dvxc; + kernel.apply(a, dvxc); + + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, XC_LDA_XC_KSDT, XC_UNPOLARIZED), 0); + const double theta = 0.007; + xc_func_set_ext_params(&func, &theta); + double rho = 0.1; + double v2rho2 = 0.0; + xc_lda_fxc(&func, 1, &rho, &v2rho2); + xc_func_end(&func); + const double k_rr = 2.0 * v2rho2; + + double max_err = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + max_err = std::max(max_err, std::abs(dvxc[ir] - k_rr * a[ir])); + } + ASSERT_GT(std::abs(k_rr), 1.0e-6); + EXPECT_LT(max_err / std::abs(k_rr), 1.0e-12); +} + +TEST(DFPTXcAnalytic, SignedComponentWeightsEnterEveryKernelField) +{ + // signed combination factors (the BLYP_LR pattern subtracts components) + // must scale every derivative field; with a uniform density only the + // local field is visible, so check k_rr against the weighted sum + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0(pw.nrxx, 0.1); + const std::vector func_id = {XC_GGA_X_PBE, XC_LDA_C_PZ}; + const std::vector> ext_params = {{}, {}}; + const std::vector weights = {-2.0, 3.0}; + ModuleDFPT::XC_First_Order_Analytic kernel(func_id, 2, rho0, &pw, ext_params, weights); + + std::vector> a(pw.nrxx, std::complex(0.05, 0.0)); + std::vector> dvxc; + kernel.apply(a, dvxc); + + double k_rr = 0.0; + const double weights_ref[] = {-2.0, 3.0}; + const int ids_ref[] = {XC_GGA_X_PBE, XC_LDA_C_PZ}; + for (int i = 0; i < 2; ++i) + { + xc_func_type func; + ASSERT_EQ(xc_func_init(&func, ids_ref[i], XC_UNPOLARIZED), 0); + double rho = 0.1; + if (ids_ref[i] == XC_GGA_X_PBE) + { + double sig = 0.0; + double v2rho2 = 0.0; + double v2rhosigma = 0.0; + double v2sigma2 = 0.0; + xc_gga_fxc(&func, 1, &rho, &sig, &v2rho2, &v2rhosigma, &v2sigma2); + k_rr += weights_ref[i] * 2.0 * v2rho2; + } + else + { + double v2rho2 = 0.0; + xc_lda_fxc(&func, 1, &rho, &v2rho2); + k_rr += weights_ref[i] * 2.0 * v2rho2; + } + xc_func_end(&func); + } + ASSERT_GT(std::abs(k_rr), 1.0e-6); + double max_err = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + max_err = std::max(max_err, std::abs(dvxc[ir] - k_rr * a[ir])); + } + EXPECT_LT(max_err / std::abs(k_rr), 1.0e-12); + + // configured ctor with neutral configuration must reproduce the + // plain ctor bit-for-bit + const std::vector plain_id = {XC_GGA_X_PBE, XC_GGA_C_PBE}; + const std::vector rho_g = gaussian_rho(pw, 0.3); + const std::vector> no_params = {{}, {}}; + const std::vector unit_weights = {1.0, 1.0}; + ModuleDFPT::XC_First_Order_Analytic plain(plain_id, 2, rho_g, &pw); + ModuleDFPT::XC_First_Order_Analytic configured(plain_id, 2, rho_g, &pw, no_params, unit_weights); + const ModuleBase::Vector3 q0(0.0, 0.0, 0.0); + plain.set_context(q0); + configured.set_context(q0); + const std::vector amp = cosine_amplitude(pw, 0.02); + std::vector> aa(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + aa[ir] = std::complex(amp[ir], 0.3 * amp[ir]); + } + std::vector> d_plain; + std::vector> d_conf; + plain.apply(aa, d_plain); + configured.apply(aa, d_conf); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + EXPECT_EQ(d_plain[ir], d_conf[ir]) << "ir=" << ir; + } +} + #endif // __LIBXC