diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 4956896a05..29d0be3276 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -434,7 +434,8 @@ OBJS_DFPT=dfpt_hamilt_shift.o\ dfpt_q0.o\ dfpt_q0_pos.o\ dfpt_rho.o\ - dfpt_stern.o + dfpt_stern.o\ + dfpt_xc_analytic.o OBJS_HAMILT_LCAO=hamilt_lcao.o\ hamilt_lcao_factory.o\ diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index d46bc7709f..79637856b4 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -13,15 +13,20 @@ #include "source_base/macros.h" #include "source_base/math_polyint.h" #include "source_base/tool_quit.h" +#include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/parameter.h" #include "source_pw/module_dfpt/dfpt_pert.h" #include "source_pw/module_dfpt/dfpt_pw.h" #include "source_pw/module_dfpt/dfpt_rho.h" +#ifdef __LIBXC +#include "source_pw/module_dfpt/dfpt_xc_analytic.h" +#endif #include #include #include #include +#include #include namespace { @@ -296,9 +301,52 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell) veff_r[ir] = veff_smooth(0, ir); } - // first-order XC kernel adapter around the converged ground-state density + // input guard: the first-order XC kernel supports LDA (1) and GGA (2) + // only. Meta-GGA needs the tau channel (absent in both kernel paths), + // and hybrid functionals need the exact-exchange response -- neither + // exists yet, so reject them explicitly instead of producing silently + // wrong frequencies. + const int xc_func_type = XC_Functional::get_func_type(); + if (xc_func_type != 1 && xc_func_type != 2) + { + std::stringstream msg; + msg << "DFPT supports LDA/GGA functionals only; the ground-state functional '" + << ucell.atoms[0].ncpp.xc_func << "' (func_type " << xc_func_type + << ") needs response channels that are not implemented"; + ModuleBase::WARNING_QUIT("ESolver_DFPT_PW::init_dfpt", msg.str()); + } + + // first-order XC kernel around the converged ground-state density: + // with libxc use the analytic kernel (libxc second derivatives in the + // dgradcorr form; the gradient operator acts as i(G + q) on the + // physical carrier, so the GGA response is exact at any q); without + // libxc fall back to the amplitude-grid finite differences (exact for + // LDA; at q != 0 the GGA response lacks the i(q . f_grad) A carrier + // term). +#ifdef __LIBXC + // kernel linearization point: the ground-state XC density, i.e. valence + // plus the nonlinear core correction density when the pseudopotentials + // carry one (PotXC_FDM adds chg->rho_core the same way); the first-order + // density itself stays valence-only (a frozen core does not respond) + const Charge* chg0 = this->pelec->charge; + const bool has_core = chg0->rho_core != nullptr; + std::vector rho0_r(this->pw_rho->nrxx, 0.0); + for (int ir = 0; ir < this->pw_rho->nrxx; ++ir) + { + 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(), + XC_Functional::get_func_type(), + rho0_r, + this->pw_rho); +#else + ModuleBase::WARNING("ESolver_DFPT_PW::init_dfpt", + "built without libxc: the DFPT XC kernel falls back to the finite-difference " + "amplitude path (LDA-only support; the GGA response at q != 0 lacks the " + "i(q.f_grad)A carrier term)"); xc_adapter_ = new XC_First_Order_FDM(this->inp_->nspin, PARAM.globalv.domag, PARAM.globalv.domag_z, this->inp_->gga_grad, this->inp_->out_elf[0] > 0, this->inp_->test_charge, this->pw_rho, this->pelec->charge, &ucell); +#endif if (getenv("DFPT_VKB") != nullptr) { diff --git a/source/source_pw/module_dfpt/CMakeLists.txt b/source/source_pw/module_dfpt/CMakeLists.txt index 4c9398d746..5e0c96f2b2 100644 --- a/source/source_pw/module_dfpt/CMakeLists.txt +++ b/source/source_pw/module_dfpt/CMakeLists.txt @@ -33,10 +33,19 @@ set(HEADERS dfpt_q0.h dfpt_metal.h dfpt_hamilt_shift.h + dfpt_xc_analytic.h ) add_library(${MODULE_NAME} ${SOURCES} ${HEADERS}) +# analytic libxc kernel for the DFPT first-order XC response; the source is +# an empty translation unit without libxc, so an unconditional build keeps +# every consumer linkable either way +if(ENABLE_LIBXC) + target_sources(${MODULE_NAME} PRIVATE dfpt_xc_analytic.cpp) + target_link_libraries(${MODULE_NAME} PUBLIC Libxc::xc) +endif() + target_include_directories(${MODULE_NAME} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR} diff --git a/source/source_pw/module_dfpt/dfpt_pert.cpp b/source/source_pw/module_dfpt/dfpt_pert.cpp index ad3a2628ba..c2dc764089 100644 --- a/source/source_pw/module_dfpt/dfpt_pert.cpp +++ b/source/source_pw/module_dfpt/dfpt_pert.cpp @@ -414,45 +414,4 @@ void DFPT_Pert::d2vloc_r(int atom_idx, int da, int db, std::vector& field, DFPT_PW_Data& data) -{ - ModuleBase::TITLE("DFPT_Pert", "build_efield"); - ModuleBase::timer::start("DFPT_Pert", "build_efield"); - // first-order electric-field potential: delta V(r) = - r . E (q=0 limit, - // position operator in the periodic cell). Computed directly on the shared - // real-space grid. Only relevant for the Q0 dielectric response (C6). - if (pw_rho_ == nullptr) - { - ModuleBase::timer::end("DFPT_Pert", "build_efield"); - return; - } - if (pw_rho_->gamma_only) - { - ModuleBase::WARNING_QUIT("DFPT_Pert::build_efield", - "DFPT requires a complex (gamma_only=false) real-space basis."); - } - std::vector> dv_real(pw_rho_->nrxx, std::complex(0.0, 0.0)); - const ModuleBase::Matrix3& latvec = ucell_->latvec; - const double lat0 = ucell_->lat0; - // shared real-space grid layout (serial pool): ir = (ix*ny + iy)*nz + iz, - // i.e. z runs fastest (verified against the impulse response of the FFT). - for (int ir = 0; ir < pw_rho_->nrxx; ++ir) - { - const int iz = ir % pw_rho_->nz; - const int rem = ir / pw_rho_->nz; - const int iy = rem % pw_rho_->ny; - const int ix = rem / pw_rho_->ny; - const double fx = ix / static_cast(pw_rho_->nx); - const double fy = iy / static_cast(pw_rho_->ny); - const double fz = iz / static_cast(pw_rho_->nz); - ModuleBase::Vector3 r; - r.x = (fx * latvec.e11 + fy * latvec.e12 + fz * latvec.e13) * lat0; - r.y = (fx * latvec.e21 + fy * latvec.e22 + fz * latvec.e23) * lat0; - r.z = (fx * latvec.e31 + fy * latvec.e32 + fz * latvec.e33) * lat0; - dv_real[ir] = -(field * r); // -e r.E (e absorbed in field convention) - } - data.set_dv_rc(0, 0, dv_real); - ModuleBase::timer::end("DFPT_Pert", "build_efield"); -} - } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pert.h b/source/source_pw/module_dfpt/dfpt_pert.h index 8da6662696..0581810676 100644 --- a/source/source_pw/module_dfpt/dfpt_pert.h +++ b/source/source_pw/module_dfpt/dfpt_pert.h @@ -36,8 +36,6 @@ class DFPT_Pert void apply_dv(int q_idx, int k_idx, const psi::Psi>& psi, DFPT_PW_Data& data); - void build_efield(const ModuleBase::Vector3& field, DFPT_PW_Data& data); - /// C5: real-space kernel of the same-atom second-order LOCAL potential /// d2V_loc(r) = d^2 V_loc / d tau_{da} d tau_{db} (both displacement /// dressings e^{i q.R} multiply on the SAME atom, so the cell sum diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 0355722dc5..a5f3590509 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -50,6 +51,9 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, { return; } + // per-q context: analytic GGA kernels act with i(G + q) on the + // physical carrier of the q-shifted amplitude + xc_->set_context(q_cart); std::vector> a_r(nrxx); pw_rho_->recip2real(drho_in_g.data(), a_r.data()); std::vector> b_r; @@ -207,6 +211,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) bool converged = false; double residual = 0.0; + int n_iter_done = 0; std::vector> v_sc_r_last; for (int iter = 0; iter < max_iter_ && !converged; ++iter) { @@ -251,12 +256,27 @@ 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_); + // 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; if (dbg) { std::cout << "DBG iter=" << iter << " residual=" << residual << " conv_thr=" << conv_thr_ << std::endl; } converged = (residual < conv_thr_); + ++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; 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.h b/source/source_pw/module_dfpt/dfpt_rho.h index a21c26686a..2d4876f3ac 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -37,6 +37,18 @@ class XC_First_Order public: virtual ~XC_First_Order() = default; + /// Optional per-q context before the first apply of a new phonon q + /// (cartesian, ucell::G units). Analytic GGA kernels need the phonon q + /// so the gradient operator acts as i(G + q) on the physical carrier; + /// amplitude-grid implementations ignore it (their gradient of the + /// carrier-free amplitude is only exact at q = 0). Const like apply: + /// the kernel is held as pointer-to-const by the solver, and context + /// injection between solves is the sanctioned state change. + virtual void set_context(const ModuleBase::Vector3& q_cart) const + { + (void)q_cart; + } + /// dvxc_r(r) = delta V_xc[drho_r](r), complex q-shifted amplitude on /// the shared real-space grid. Implementations must not resize or /// alias drho_r; dvxc_r is resized to drho_r.size() and fully diff --git a/source/source_pw/module_dfpt/dfpt_stern.h b/source/source_pw/module_dfpt/dfpt_stern.h index 325e6583f0..bc48aa666d 100644 --- a/source/source_pw/module_dfpt/dfpt_stern.h +++ b/source/source_pw/module_dfpt/dfpt_stern.h @@ -18,8 +18,10 @@ namespace ModuleDFPT * * The shifted Hamiltonian action is injected through LinearOperator so the * solver core stays decoupled from the ground-state operator chain: the - * production adapter reuses hamilt::Hamilt::ops->hPsi at the k+q point - * (wired in C7), while unit tests supply analytic operators. + * production adapter (DFPT_HamiltShift, dfpt_hamilt_shift.h) builds the + * (H(k+q) - eps) action itself on the k+q basis rather than reusing + * hamilt::Hamilt::ops->hPsi (a potential future refactor to shrink the + * reimplementation surface); unit tests supply analytic operators. */ class DFPT_Stern { diff --git a/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp b/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp new file mode 100644 index 0000000000..b485515e77 --- /dev/null +++ b/source/source_pw/module_dfpt/dfpt_xc_analytic.cpp @@ -0,0 +1,241 @@ +#ifdef __LIBXC + +#include "dfpt_xc_analytic.h" + +#include "source_base/constants.h" +#include "source_base/tool_quit.h" + +#include + +#include +#include +#include + +namespace ModuleDFPT +{ + +namespace +{ +// QE dgradcorr thresholds: no gradient-kernel response where the +// ground-state density or its gradient is negligible +const double epsr = 1.0e-6; +const double epsg = 1.0e-10; + +// libxc works in Hartree; the DFPT potentials are in Ry (e2 = 2) +const double ha_to_ry = 2.0; +} // namespace + +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) + : func_type_(func_type), func_id_(func_id), pw_rho_(pw_rho) +{ + ModuleBase::TITLE("XC_First_Order_Analytic", "XC_First_Order_Analytic"); + if (func_type != 1 && func_type != 2) + { + ModuleBase::WARNING_QUIT("XC_First_Order_Analytic", "only LDA (1) and GGA (2) kernels are supported"); + } + if (func_id_.empty() || pw_rho_ == nullptr || static_cast(rho0_r.size()) != pw_rho_->nrxx) + { + ModuleBase::WARNING_QUIT("XC_First_Order_Analytic", "invalid func_id / basis / rho0 size"); + } + this->init_gradient_(rho0_r); + this->init_kernel_fields_(rho0_r); +} + +void XC_First_Order_Analytic::set_context(const ModuleBase::Vector3& q_cart) const +{ + q_cart_ = q_cart; + q_set_ = true; +} + +void XC_First_Order_Analytic::init_gradient_(const std::vector& rho0_r) +{ + const int nrxx = pw_rho_->nrxx; + const int npw = pw_rho_->npw; + std::vector> rho0_g(npw); + pw_rho_->real2recip(rho0_r.data(), rho0_g.data()); + grho_.assign(3, std::vector(nrxx, 0.0)); + std::vector> g_g(npw); + for (int ipol = 0; ipol < 3; ++ipol) + { + for (int ig = 0; ig < npw; ++ig) + { + // pw gcar is dimensionless (2 pi not folded in): the physical + // wavevector is tpiba * gcar, in 1/lat0 + const double g = pw_rho_->tpiba * pw_rho_->gcar[ig][ipol]; + g_g[ig] = std::complex(0.0, g) * rho0_g[ig]; + } + std::vector> g_r(nrxx); + pw_rho_->recip2real(g_g.data(), g_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + grho_[ipol][ir] = g_r[ir].real(); + } + } +} + +void XC_First_Order_Analytic::init_kernel_fields_(const std::vector& rho0_r) +{ + const int nrxx = pw_rho_->nrxx; + k_rr_.assign(nrxx, 0.0); + k_sr_.assign(nrxx, 0.0); + k_ss_.assign(nrxx, 0.0); + k_s_.assign(nrxx, 0.0); + + // sigma = |grad rho0|^2 for the GGA branch + std::vector sigma(nrxx, 0.0); + if (func_type_ == 2) + { + for (int ir = 0; ir < nrxx; ++ir) + { + for (int ipol = 0; ipol < 3; ++ipol) + { + sigma[ir] += grho_[ipol][ir] * grho_[ipol][ir]; + } + } + } + + // open every functional component once, accumulate the fields, close + for (size_t ifunc = 0; ifunc < func_id_.size(); ++ifunc) + { + xc_func_type func; + if (xc_func_init(&func, func_id_[ifunc], XC_UNPOLARIZED) != 0) + { + ModuleBase::WARNING_QUIT("XC_First_Order_Analytic", "xc_func_init failed"); + } + if (func_type_ == 1) + { + for (int ir = 0; ir < nrxx; ++ir) + { + double v2rho2 = 0.0; + const double rho = rho0_r[ir]; + 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; + } + } + else + { + for (int ir = 0; ir < nrxx; ++ir) + { + const double rho = rho0_r[ir]; + const double sig = sigma[ir]; + if (std::abs(rho) <= epsr || sig <= epsg) + { + continue; // QE dgradcorr guard: no kernel response here + } + double zk = 0.0; + double vrho = 0.0; + double vsigma = 0.0; + xc_gga_exc_vxc(&func, 1, &rho, &sig, &zk, &vrho, &vsigma); + double v2rho2 = 0.0; + double v2rhosigma = 0.0; + double v2sigma2 = 0.0; + xc_gga_fxc(&func, 1, &rho, &sig, &v2rho2, &v2rhosigma, &v2sigma2); + // 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); + } + } + xc_func_end(&func); + } +} + +void XC_First_Order_Analytic::apply(const std::vector>& drho_r, + std::vector>& dvxc_r) const +{ + const int nrxx = pw_rho_->nrxx; + const int npw = pw_rho_->npw; + if (static_cast(drho_r.size()) != nrxx) + { + ModuleBase::WARNING_QUIT("XC_First_Order_Analytic", "drho_r is not on the rho grid"); + } + dvxc_r.assign(nrxx, std::complex(0.0, 0.0)); + if (func_type_ == 1) + { + // LDA: dvxc = f_rho A, no q dependence + for (int ir = 0; ir < nrxx; ++ir) + { + dvxc_r[ir] = k_rr_[ir] * drho_r[ir]; + } + return; + } + + // GGA: dgradcorr with the total-field gradient i (G + q) + const double qx = q_set_ ? q_cart_.x : 0.0; + const double qy = q_set_ ? q_cart_.y : 0.0; + const double qz = q_set_ ? q_cart_.z : 0.0; + std::vector> a_g(npw); + pw_rho_->real2recip(drho_r.data(), a_g.data()); + + // gdrho(ipol, r) = F^-1[ i tpiba (G + q)_ipol A(G) ]: physical + // gradient of the carrier-carrying field, 1/lat0 units + std::vector>> gdrho(3, std::vector>(nrxx, std::complex(0.0, 0.0))); + std::vector> g_g(npw); + for (int ipol = 0; ipol < 3; ++ipol) + { + const double q_ipol = ipol == 0 ? qx : (ipol == 1 ? qy : qz); + for (int ig = 0; ig < npw; ++ig) + { + const double w = pw_rho_->tpiba * (pw_rho_->gcar[ig][ipol] + q_ipol); + g_g[ig] = std::complex(0.0, w) * a_g[ig]; + } + pw_rho_->recip2real(g_g.data(), gdrho[ipol].data()); + } + + // s1 = grad rho0 . gdrho + std::vector> s1(nrxx, std::complex(0.0, 0.0)); + for (int ir = 0; ir < nrxx; ++ir) + { + for (int ipol = 0; ipol < 3; ++ipol) + { + s1[ir] += grho_[ipol][ir] * gdrho[ipol][ir]; + } + } + + // local part and the h field whose q-gradient completes the response + std::vector>> h(3, std::vector>(nrxx, std::complex(0.0, 0.0))); + for (int ir = 0; ir < nrxx; ++ir) + { + dvxc_r[ir] = k_rr_[ir] * drho_r[ir] + k_sr_[ir] * s1[ir]; + const std::complex fac = k_sr_[ir] * drho_r[ir] + k_ss_[ir] * s1[ir]; + for (int ipol = 0; ipol < 3; ++ipol) + { + h[ipol][ir] = fac * grho_[ipol][ir] + k_s_[ir] * gdrho[ipol][ir]; + } + } + + // dvxc -= F^-1[ i tpiba (G + q) . h(G) ] + std::vector> dh_g(npw, std::complex(0.0, 0.0)); + std::vector> h_g(npw); + for (int ipol = 0; ipol < 3; ++ipol) + { + const double q_ipol = ipol == 0 ? qx : (ipol == 1 ? qy : qz); + pw_rho_->real2recip(h[ipol].data(), h_g.data()); + for (int ig = 0; ig < npw; ++ig) + { + const double w = pw_rho_->tpiba * (pw_rho_->gcar[ig][ipol] + q_ipol); + dh_g[ig] += std::complex(0.0, w) * h_g[ig]; + } + } + std::vector> dh_r(nrxx); + pw_rho_->recip2real(dh_g.data(), dh_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + dvxc_r[ir] -= dh_r[ir]; + } +} + +} // namespace ModuleDFPT + +#endif // __LIBXC diff --git a/source/source_pw/module_dfpt/dfpt_xc_analytic.h b/source/source_pw/module_dfpt/dfpt_xc_analytic.h new file mode 100644 index 0000000000..a2d9d4bb7a --- /dev/null +++ b/source/source_pw/module_dfpt/dfpt_xc_analytic.h @@ -0,0 +1,75 @@ +#ifndef DFPT_XC_ANALYTIC_H +#define DFPT_XC_ANALYTIC_H + +#ifdef __LIBXC + +#include "dfpt_rho.h" + +#include "source_base/matrix3.h" +#include "source_basis/module_pw/pw_basis.h" + +#include + +namespace ModuleDFPT +{ + +/** + * @brief Analytic first-order XC kernel built from libxc second derivatives + * (C6 provider; replaces the finite-difference adapter when ABACUS + * is compiled with libxc). + * + * The response to the q-shifted density amplitude A(r) of the physical + * first-order density drho(r) = 2 Re[A(r) e^{i q r}] is evaluated in the + * dgradcorr form (QE LR_Modules/dgradcorr.f90): the gradient operator acts + * on the total physical field through i G -> i (G + q), so the carrier + * contribution i (q . f_grad) A -- exact for any q, and structurally + * missing from grid-only finite differences of the carrier-free amplitude + * -- is included. LDA functionals reduce to dvxc = f_rho A (their gradient + * kernel vanishes identically, no q dependence). + * + * The libxc kernel fields are evaluated once at construction on the + * ground-state density rho0 and stored in Ry units (libxc Hartree values + * scaled by e2 = 2). Real-space points with |rho0| <= epsr or + * |grad rho0|^2 <= epsg follow the QE dgradcorr guard: the gradient-kernel + * fields vanish there. Construction is the only libxc entry point + * (handles are opened and closed inside the constructor), so this header + * does not include xc.h. + */ +class XC_First_Order_Analytic : public XC_First_Order +{ + public: + XC_First_Order_Analytic(const std::vector& func_id, + const int func_type, + const std::vector& rho0_r, + ModulePW::PW_Basis* pw_rho); + + void set_context(const ModuleBase::Vector3& q_cart) const override; + + void apply(const std::vector>& drho_r, + std::vector>& dvxc_r) const override; + + private: + /// grho_(ipol)[ir] = d rho0 / d r_ipol on the rho grid (real part) + void init_gradient_(const std::vector& rho0_r); + + /// kernel fields k_rr_/k_sr_/k_ss_/k_s_ (LDA fills k_rr_ only) + void init_kernel_fields_(const std::vector& rho0_r); + + int func_type_ = 1; ///< 1 = LDA, 2 = GGA (XC_Functional convention) + std::vector func_id_; ///< libxc ids, exchange first (XC_Functional order) + ModulePW::PW_Basis* pw_rho_ = nullptr; + + mutable ModuleBase::Vector3 q_cart_; ///< phonon q, ucell::G units (1/lat0) + mutable bool q_set_ = false; ///< false: apply() behaves as q = 0 (Gamma) + + std::vector> grho_; ///< [3][nrxx] ground-state gradient + std::vector k_rr_; ///< rho-rho kernel field, Ry units + std::vector k_sr_; ///< rho-sigma field (s1 = grho.gdrho coefficient) + std::vector k_ss_; ///< sigma-sigma field (s1^2 coefficient) + std::vector k_s_; ///< sigma field (gdrho coefficient) +}; + +} // namespace ModuleDFPT + +#endif // __LIBXC +#endif // DFPT_XC_ANALYTIC_H diff --git a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt index b5cfc3ef70..c91bc36d2f 100644 --- a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt +++ b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt @@ -106,3 +106,25 @@ AddTest( # Shared serial-side cell/basis/data fixture. dfpt_serial_fixture.cpp ) + +# Analytic libxc first-order XC kernel (empty TU without libxc, so the test +# only exists when the feature is compiled in) +if(ENABLE_LIBXC) + AddTest( + TARGET MODULE_DFPT_xc_analytic_serial + LIBS parameter dfpt_planewave_serial device base symmetry Libxc::xc + SOURCES dfpt_xc_analytic_serial_test.cpp + ../dfpt_xc_analytic.cpp + ../dfpt_rho.cpp + ../dfpt_pw_data.cpp + ../dfpt_kq_basis.cpp + ../../../source_cell/qlist.cpp + ../../../source_cell/reciprocal_grid.cpp + ../../../source_psi/psi.cpp + # Shared ctor/dtor stubs for the cell/spepot/Plus_U link closures + # (see test/dfpt_test_mocks.cpp). + ../test/dfpt_test_mocks.cpp + # Shared serial-side cell/basis/data fixture. + dfpt_serial_fixture.cpp + ) +endif() diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_pert_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_pert_serial_test.cpp index cff890dbb9..06a2f76b14 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_pert_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_pert_serial_test.cpp @@ -41,7 +41,7 @@ * - build_dv + apply_dv: the stored reciprocal coefficients and the FFT * convolution dpsi(G'') = sum_G' psi(G') c(G''-G') are checked against * the analytic matrix elements on the k+q basis (no extra q-phase). - * - build_efield: -E.r ramp layout and its closed-form discrete FT. + * - build_vkb: l=0 projector against an independent hand-rolled Simpson * transform; pure phase response under a tau shift. * - dVnl_dtau: the two-term NC identity is validated against a finite @@ -178,66 +178,6 @@ TEST_F(DFPTPertSerialTest, ApplyDvConvolutionMatchesAnalyticMatrixElement) } } -TEST_F(DFPTPertSerialTest, BuildEfieldRampMatchesClosedForm) -{ - const double E0 = 0.5; - const double Lx = 10.0 * lat0_; - pert_.build_efield(ModuleBase::Vector3(E0, 0.0, 0.0), data_); - const std::vector> dv = data_.get_dv_rc(0, 0); - ASSERT_EQ(dv.size(), static_cast(pw_rho_.nrxx)); - const int nx = pw_rho_.nx; - - // real-space values form a pure x ramp -E0 (ix/nx) Lx - for (int ir = 0; ir < pw_rho_.nrxx; ++ir) - { - const double j = std::llround(-dv[ir].real() / (E0 * Lx) * nx); - EXPECT_GE(j, 0.0); - EXPECT_LE(j, nx - 1.0); - EXPECT_NEAR(dv[ir].real(), -E0 * (j / nx) * Lx, 1.0e-10); - EXPECT_EQ(dv[ir].imag(), 0.0); - } - - // closed-form discrete FT of the sawtooth discriminates the (ix,iy,iz) - // layout: c(k) = E0 Lx / (nx (1 - w)), w = exp(-2 pi i k / nx), k != 0 - std::vector> recip(pw_rho_.npw); - pw_rho_.real2recip(dv.data(), recip.data()); - int ig_p100 = -1, ig_m100 = -1, ig_p110 = -1, ig_p010 = -1; - for (int ig = 0; ig < pw_rho_.npw; ++ig) - { - const long long key = GKey(pw_rho_.gcar[ig]); - if (key == FKey(1, 0, 0)) - { - ig_p100 = ig; - } - else if (key == FKey(-1, 0, 0)) - { - ig_m100 = ig; - } - else if (key == FKey(1, 1, 0)) - { - ig_p110 = ig; - } - else if (key == FKey(0, 1, 0)) - { - ig_p010 = ig; - } - } - ASSERT_GE(ig_p100, 0); - ASSERT_GE(ig_m100, 0); - ASSERT_GE(ig_p110, 0); - ASSERT_GE(ig_p010, 0); - const std::complex w1 = std::polar(1.0, -ModuleBase::TWO_PI / nx); - const std::complex nx_c(nx); - const std::complex c1 = E0 * Lx / (nx_c * (1.0 - w1)); - const std::complex cm1 = E0 * Lx / (nx_c * (1.0 - std::conj(w1))); - EXPECT_NEAR(recip[ig_p100].real(), c1.real(), 1.0e-8); - EXPECT_NEAR(recip[ig_p100].imag(), c1.imag(), 1.0e-8); - EXPECT_NEAR(recip[ig_m100].real(), cm1.real(), 1.0e-8); - EXPECT_NEAR(recip[ig_m100].imag(), cm1.imag(), 1.0e-8); - EXPECT_NEAR(std::abs(recip[ig_p110]), 0.0, 1.0e-10); - EXPECT_NEAR(std::abs(recip[ig_p010]), 0.0, 1.0e-10); -} - TEST_F(DFPTPertSerialTest, BuildVkbL0MatchesIndependentSimpson) { MakeNCAtom(); 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 new file mode 100644 index 0000000000..d95b489da6 --- /dev/null +++ b/source/source_pw/module_dfpt/test_serial/dfpt_xc_analytic_serial_test.cpp @@ -0,0 +1,382 @@ +#include + +#ifdef __LIBXC + +#include "source_pw/module_dfpt/dfpt_xc_analytic.h" + +#include "source_base/constants.h" +#include "source_basis/module_pw/pw_basis.h" + +#include + +#include +#include +#include + +namespace +{ + +// simple serial PW_Basis on a cubic cell (rho-grid only); PW_Basis owns raw +// FFT state, so it is initialized in place instead of returned by value +void make_basis(ModulePW::PW_Basis& pw, const double lat0, const double a, const double ecut) +{ + const ModuleBase::Matrix3 latvec(a, 0.0, 0.0, 0.0, a, 0.0, 0.0, 0.0, a); + pw.initgrids(lat0, latvec, ecut); + pw.initparameters(false, ecut); + pw.fft_bundle.initfftmode(0); + pw.setuptransform(); + pw.collect_local_pw(); +} + +// smooth positive density: gaussian bump at the cell center +std::vector gaussian_rho(const ModulePW::PW_Basis& pw, const double amp) +{ + std::vector rho(pw.nrxx, 0.0); + for (int ix = 0; ix < pw.nx; ++ix) + { + for (int iy = 0; iy < pw.ny; ++iy) + { + for (int iz = 0; iz < pw.nz; ++iz) + { + const int ir = (ix * pw.ny + iy) * pw.nz + iz; + const double x = static_cast(ix) / pw.nx - 0.5; + const double y = static_cast(iy) / pw.ny - 0.5; + const double z = static_cast(iz) / pw.nz - 0.5; + rho[ir] = amp * std::exp(-(x * x + y * y + z * z) / 0.08); + } + } + } + return rho; +} + +// smooth real amplitude with gradient content +std::vector cosine_amplitude(const ModulePW::PW_Basis& pw, const double amp) +{ + std::vector a(pw.nrxx, 0.0); + for (int ix = 0; ix < pw.nx; ++ix) + { + for (int iy = 0; iy < pw.ny; ++iy) + { + for (int iz = 0; iz < pw.nz; ++iz) + { + const int ir = (ix * pw.ny + iy) * pw.nz + iz; + const double x = ModuleBase::TWO_PI * ix / pw.nx; + const double y = ModuleBase::TWO_PI * iy / pw.ny; + a[ir] = amp * (std::cos(x) + 0.5 * std::cos(y) + 0.3 * std::cos(x + y) * std::cos(2.0 * x)); + } + } + } + return a; +} + +// FFT gradient of a real field on the rho grid (physical units: pw gcar is +// dimensionless, the physical wavevector is tpiba * gcar) +void fft_gradient(const ModulePW::PW_Basis& pw, const std::vector& f, std::vector>& grad) +{ + const int npw = pw.npw; + std::vector> f_g(npw); + pw.real2recip(f.data(), f_g.data()); + grad.assign(3, std::vector(pw.nrxx, 0.0)); + std::vector> g_g(npw); + for (int ipol = 0; ipol < 3; ++ipol) + { + for (int ig = 0; ig < npw; ++ig) + { + g_g[ig] = std::complex(0.0, pw.tpiba * pw.gcar[ig][ipol]) * f_g[ig]; + } + std::vector> g_r(pw.nrxx); + pw.recip2real(g_g.data(), g_r.data()); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + grad[ipol][ir] = g_r[ir].real(); + } + } +} + +// divergence of a real vector field on the rho grid (physical units) +void fft_divergence(const ModulePW::PW_Basis& pw, + const std::vector>& h, + std::vector& div) +{ + const int npw = pw.npw; + std::vector> dh_g(npw, std::complex(0.0, 0.0)); + std::vector> h_g(npw); + for (int ipol = 0; ipol < 3; ++ipol) + { + pw.real2recip(h[ipol].data(), h_g.data()); + for (int ig = 0; ig < npw; ++ig) + { + dh_g[ig] += std::complex(0.0, pw.tpiba * pw.gcar[ig][ipol]) * h_g[ig]; + } + } + div.resize(pw.nrxx); + std::vector> dh_r(pw.nrxx); + pw.recip2real(dh_g.data(), dh_r.data()); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + div[ir] = dh_r[ir].real(); + } +} + +// sum of libxc vrho / vsigma over functional components (Hartree) +void libxc_vxc(const std::vector& func_id, + const std::vector& rho, + const std::vector& sigma, + std::vector& vrho, + std::vector& vsigma) +{ + const int n = rho.size(); + vrho.assign(n, 0.0); + vsigma.assign(n, 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); + std::vector vr(n, 0.0); + std::vector vs(n, 0.0); + std::vector zk(n, 0.0); + xc_gga_exc_vxc(&func, n, rho.data(), sigma.data(), zk.data(), vr.data(), vs.data()); + for (int ir = 0; ir < n; ++ir) + { + vrho[ir] += vr[ir]; + vsigma[ir] += vs[ir]; + } + xc_func_end(&func); + } +} + +} // namespace + +TEST(DFPTXcAnalytic, LdaKernelMatchesVrhofiniteDifference) +{ + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0 = gaussian_rho(pw, 0.3); + const std::vector func_id = {XC_LDA_X, XC_LDA_C_PZ}; + ModuleDFPT::XC_First_Order_Analytic kernel(func_id, 1, rho0, &pw); + + // constant amplitude A = 1: dvxc[ir] must equal k_rr[ir] + std::vector> a(pw.nrxx, std::complex(1.0, 0.0)); + std::vector> dvxc; + kernel.apply(a, dvxc); + + // reference: central difference of the LDA potential. libxc vxc is + // already d(rho e)/d rho, so the kernel is its plain derivative + const double eta = 1.0e-5; + std::vector rho_p(pw.nrxx); + std::vector rho_m(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + rho_p[ir] = rho0[ir] + eta; + rho_m[ir] = rho0[ir] - eta; + } + std::vector v_p(pw.nrxx, 0.0); + std::vector v_m(pw.nrxx, 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); + std::vector zk(pw.nrxx, 0.0); + std::vector vxc(pw.nrxx, 0.0); + xc_lda_exc_vxc(&func, pw.nrxx, rho_p.data(), zk.data(), vxc.data()); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + v_p[ir] += vxc[ir]; + } + xc_lda_exc_vxc(&func, pw.nrxx, rho_m.data(), zk.data(), vxc.data()); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + v_m[ir] += vxc[ir]; + } + xc_func_end(&func); + } + // sigma = |grad rho0|^2: the FD reference has no kernel guards, so the + // comparison is restricted to the region where the analytic kernel is + // active (rho > epsr, sigma > epsg) and dense enough that the FD + // truncation error (grows like rho^-7/3 in the tail) stays controlled + std::vector> grho0_l; + fft_gradient(pw, rho0, grho0_l); + double max_abs = 0.0; + double max_err = 0.0; + int n_checked = 0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + if (rho0[ir] < 0.01) + { + continue; + } + double s2 = 0.0; + for (int ipol = 0; ipol < 3; ++ipol) + { + s2 += grho0_l[ipol][ir] * grho0_l[ipol][ir]; + } + if (s2 < 1.0e-8) + { + continue; + } + const double ref = 2.0 * (v_p[ir] - v_m[ir]) / (2.0 * eta); // Ry + max_abs = std::max(max_abs, std::abs(ref)); + max_err = std::max(max_err, std::abs(ref - dvxc[ir].real())); + ++n_checked; + } + EXPECT_GT(n_checked, pw.nrxx / 10); // the checked region is representative + EXPECT_GT(max_abs, 1.0); // nontrivial kernel scale (Ry bohr^3) + EXPECT_LT(max_err / max_abs, 1.0e-6); +} + +TEST(DFPTXcAnalytic, GgaQ0MatchesFieldPerturbationFiniteDifference) +{ + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0 = gaussian_rho(pw, 0.3); + 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.0); + } + std::vector> dvxc; + kernel.apply(a, dvxc); + + // reference: central difference of the full GGA potential + // v_xc = vrho - div(2 vsigma grad rho) + // between rho0 +- eta A (Hartree -> Ry by the 2.0 below) + std::vector> grho0; + fft_gradient(pw, rho0, grho0); + std::vector> gamp; + fft_gradient(pw, amp, gamp); + const double eta = 2.0e-6; + std::vector rho_p(pw.nrxx); + std::vector rho_m(pw.nrxx); + std::vector sig_p(pw.nrxx); + std::vector sig_m(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + rho_p[ir] = rho0[ir] + eta * amp[ir]; + rho_m[ir] = rho0[ir] - eta * amp[ir]; + double sp = 0.0; + double sm = 0.0; + for (int ipol = 0; ipol < 3; ++ipol) + { + const double gp = grho0[ipol][ir] + eta * gamp[ipol][ir]; + const double gm = grho0[ipol][ir] - eta * gamp[ipol][ir]; + sp += gp * gp; + sm += gm * gm; + } + sig_p[ir] = sp; + sig_m[ir] = sm; + } + std::vector vr_p; + std::vector vs_p; + std::vector vr_m; + std::vector vs_m; + libxc_vxc(func_id, rho_p, sig_p, vr_p, vs_p); + libxc_vxc(func_id, rho_m, sig_m, vr_m, vs_m); + std::vector> h_p(3, std::vector(pw.nrxx, 0.0)); + std::vector> h_m(3, std::vector(pw.nrxx, 0.0)); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + for (int ipol = 0; ipol < 3; ++ipol) + { + const double gp = grho0[ipol][ir] + eta * gamp[ipol][ir]; + const double gm = grho0[ipol][ir] - eta * gamp[ipol][ir]; + h_p[ipol][ir] = 2.0 * vs_p[ir] * gp; + h_m[ipol][ir] = 2.0 * vs_m[ir] * gm; + } + } + std::vector div_p; + std::vector div_m; + fft_divergence(pw, h_p, div_p); + fft_divergence(pw, h_m, div_m); + // same region restriction as the LDA test: compare where the kernel is + // active and the FD truncation is controlled + double l2_num = 0.0; + double l2_den = 0.0; + double max_abs = 0.0; + int n_checked = 0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + if (rho0[ir] < 0.01) + { + continue; + } + double s2 = 0.0; + for (int ipol = 0; ipol < 3; ++ipol) + { + s2 += grho0[ipol][ir] * grho0[ipol][ir]; + } + if (s2 < 1.0e-8) + { + continue; + } + const double v_p = 2.0 * (vr_p[ir] - div_p[ir]); + const double v_m = 2.0 * (vr_m[ir] - div_m[ir]); + const double ref = (v_p - v_m) / (2.0 * eta); + const double err = ref - dvxc[ir].real(); + l2_num += err * err; + l2_den += ref * ref; + max_abs = std::max(max_abs, std::abs(ref)); + ++n_checked; + } + EXPECT_GT(n_checked, pw.nrxx / 10); // the checked region is representative + EXPECT_GT(max_abs, 1.0e-3); // nontrivial response scale + // the eta^2 FD truncation of the reference bounds the agreement around + // 1e-6..1e-5; structural kernel/assembly errors show up at percent level + EXPECT_LT(std::sqrt(l2_num / l2_den), 1.0e-5); +} + +TEST(DFPTXcAnalytic, CarrierTermAppearsAtFiniteQOnlyForGga) +{ + ModulePW::PW_Basis pw; + make_basis(pw, 1.8897261255, 6.0, 20.0); + const std::vector rho0 = gaussian_rho(pw, 0.3); + const std::vector amp = cosine_amplitude(pw, 0.02); + std::vector> a(pw.nrxx); + for (int ir = 0; ir < pw.nrxx; ++ir) + { + a[ir] = std::complex(amp[ir], 0.3 * amp[ir]); + } + + // GGA: the finite-q response must differ from the Gamma response + // (the i(q.f_grad)A carrier term); LDA must stay q-independent + const std::vector gga_id = {XC_GGA_X_PBE, XC_GGA_C_PBE}; + ModuleDFPT::XC_First_Order_Analytic gga(gga_id, 2, rho0, &pw); + const ModuleBase::Vector3 q0(0.0, 0.0, 0.0); + gga.set_context(q0); + std::vector> d0; + gga.apply(a, d0); + const ModuleBase::Vector3 qx(0.25, 0.0, 0.0); + gga.set_context(qx); + std::vector> dq; + gga.apply(a, dq); + double diff = 0.0; + double scale = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + diff = std::max(diff, std::abs(dq[ir] - d0[ir])); + scale = std::max(scale, std::abs(d0[ir])); + } + EXPECT_GT(diff / scale, 1.0e-3); + + const std::vector lda_id = {XC_LDA_X, XC_LDA_C_PZ}; + ModuleDFPT::XC_First_Order_Analytic lda(lda_id, 1, rho0, &pw); + lda.set_context(q0); + std::vector> l0; + lda.apply(a, l0); + lda.set_context(qx); + std::vector> lq; + lda.apply(a, lq); + double ldiff = 0.0; + for (int ir = 0; ir < pw.nrxx; ++ir) + { + ldiff = std::max(ldiff, std::abs(lq[ir] - l0[ir])); + } + EXPECT_EQ(ldiff, 0.0); +} + +#endif // __LIBXC