feat(dfpt): iteration logging and dead-code cleanup

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\

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 <cmath>
 #include <complex>
 #include <cstdlib>
 #include <iostream>
+#include <sstream>
 #include <vector>
 
 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<double> 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)
     {

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}

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<std::complex<
     ModuleBase::timer::end("DFPT_Pert", "d2vloc_r");
 }
 
-void DFPT_Pert::build_efield(const ModuleBase::Vector3<double>& 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<std::complex<double>> dv_real(pw_rho_->nrxx, std::complex<double>(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<double>(pw_rho_->nx);
-        const double fy = iy / static_cast<double>(pw_rho_->ny);
-        const double fz = iz / static_cast<double>(pw_rho_->nz);
-        ModuleBase::Vector3<double> 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

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<std::complex<double>>& psi, DFPT_PW_Data& data);
 
-    void build_efield(const ModuleBase::Vector3<double>& 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

source/source_pw/module_dfpt/dfpt_pw_solve.cpp

@@ -17,6 +17,7 @@
 #include <cmath>
 #include <complex>
 #include <cstdlib>
+#include <iomanip>
 #include <iostream>
 #include <vector>
 
@@ -50,6 +51,9 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3<double>& 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<std::complex<double>> a_r(nrxx);
     pw_rho_->recip2real(drho_in_g.data(), a_r.data());
     std::vector<std::complex<double>> 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<std::complex<double>> 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<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;
     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_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<double>& 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

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
 {