feat(dfpt): complete the nlcc core-response terms (QE dynmat_nlcc + dynmatcc) and fix the mirrored structure-factor phase

source/source_esolver/esolver_dfpt_pw.cpp

@@ -426,10 +426,59 @@ void ESolver_DFPT_PW::init_dfpt(UnitCell& ucell)
         }
     }
 
+    // ground-state XC potential on the rho grid, evaluated on the SAME
+    // density the first-order kernel linearizes around (valence + core
+    // charge): consumed by the nlcc diagonal 2n+1 term (QE dynmatcc).
+    // Only needed when some species carries a nonlinear core correction.
+    bool nlcc_any = false;
+    for (int it = 0; it < ucell.ntype; ++it)
+    {
+        if (ucell.atoms[it].ncpp.nlcc)
+        {
+            nlcc_any = true;
+            break;
+        }
+    }
+    std::vector<double> vxc_r;
+    if (nlcc_any)
+    {
+        const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
+            = XC_Functional::v_xc(this->pw_rho->nrxx, this->pelec->charge, &ucell, 1,
+                                  false, false, this->inp_->gga_grad,
+                                  XC_Functional::get_hybrid_alpha(),
+                                  XC_Functional::get_hse_omega());
+        const ModuleBase::matrix& v_xc_mat = std::get<2>(etxc_vtxc_v);
+        if (v_xc_mat.nr > 0 && v_xc_mat.nc == this->pw_rho->nrxx)
+        {
+            vxc_r.resize(this->pw_rho->nrxx);
+            for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
+            {
+                vxc_r[ir] = v_xc_mat(0, ir);
+            }
+        }
+        else
+        {
+            ModuleBase::WARNING("ESolver_DFPT_PW::init_dfpt",
+                                "failed to evaluate the ground-state XC potential; "
+                                "the nlcc diagonal dynamical-matrix term is skipped");
+        }
+    }
+
     this->stp.sync_cpu();
-    dfpt_->init(ucell, *this->stp.psi_cpu, this->pw_rho, this->pw_wfc, &this->sf, veff_r,
-                this->pelec->wg, this->pelec->ekb, xc_adapter_, nelec_, ecutwfc_,
-                dft_plus_u_ ? this->dftu_.get() : nullptr);
+    const ModuleDFPT::DFPT_PW_InitContext ctx{&ucell,
+                                               this->stp.psi_cpu,
+                                               this->pw_rho,
+                                               this->pw_wfc,
+                                               &this->sf,
+                                               &veff_r,
+                                               &this->pelec->wg,
+                                               &this->pelec->ekb,
+                                               xc_adapter_,
+                                               nelec_,
+                                               ecutwfc_,
+                                               dft_plus_u_ ? this->dftu_.get() : nullptr,
+                                               &vxc_r};
+    dfpt_->init(ctx);
 }
 
 void ESolver_DFPT_PW::run_post_process(UnitCell& ucell)

source/source_pw/module_dfpt/dfpt_phon.cpp

@@ -11,6 +11,7 @@
 #include <cmath>
 #include <complex>
 #include <cstdlib>
+#include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
@@ -64,6 +65,21 @@ void DFPT_Phon::init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, DFPT_Pert* per
     ModuleBase::timer::end("DFPT_Phon", "init");
 }
 
+void DFPT_Phon::set_rho(const DFPT_Rho* rho)
+{
+    rho_ = rho;
+}
+
+void DFPT_Phon::set_vxc_r(const std::vector<double>* vxc_r)
+{
+    vxc_r_ = vxc_r;
+}
+
+void DFPT_Phon::set_dvxc_core_table(const std::vector<std::vector<std::complex<double>>>* table)
+{
+    dvxc_core_ = table;
+}
+
 // ---------------------------------------------------------------------------
 // assemble / diagonalize / LO-TO / sum rule
 // ---------------------------------------------------------------------------

source/source_pw/module_dfpt/dfpt_phon.h

@@ -5,7 +5,9 @@
 #include "source_cell/unitcell.h"
 #include "source_psi/psi.h"
 
+#include <complex>
 #include <string>
+#include <vector>
 
 namespace ModulePW
 {
@@ -16,6 +18,7 @@ namespace ModuleDFPT
 {
 
 class DFPT_Pert;
+class DFPT_Rho;
 
 /**
  * @brief Dynamical matrix of DFPT (C5).
@@ -35,6 +38,23 @@ class DFPT_Pert;
  * row D[b][*] is filled right after displacement b converges, so the dpsi
  * storage never needs a direction dimension (data-layer refactor reserved
  * for phase B).
+ *
+ * nlcc systems add two direct terms (QE ph.x ground truth: dynmat_nlcc
+ * and dynmatcc, both absent from the wavefunction response). Kernel
+ * term (row-conjugated, pool-summed):
+ *   X_ab = (Omega/nxyz) sum_pool conj(dvxc^a(r)) drhoc^b(r),
+ *   D_ab += (1/2)(X_ab + conj(X_ba)) / sqrt(m_a m_b)
+ * where dvxc^a = f_xc . (drho_val^a + drhoc^a) is the first-order XC
+ * potential of the TOTAL response of displacement a, supplied by the
+ * caller through the per-q kernel table (covering ALL atoms: rows of
+ * non-nlcc atoms carry their valence-only dvxc; leaving their slots
+ * empty silently drops the one-sided cross blocks). Diagonal
+ * correction (q-independent, pool-summed):
+ *   D_aa[i][j] += (Omega/nxyz) sum_pool Re[w_ij(r) vxc^0(r)] / m_a
+ * where w_ij = -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau_a), evaluated
+ * through recip2real as a real field. The ground-state XC potential
+ * vxc^0 is the valence + core charge density's XC potential, supplied
+ * through the InitContext by the esolver wiring.
  */
 class DFPT_Phon
 {
@@ -44,18 +64,38 @@ class DFPT_Phon
 
     void init(UnitCell& ucell, ModulePW::PW_Basis* pw_rho, DFPT_Pert* pert);
 
+    /// providers for the nlcc direct terms: the first-order core response
+    /// (DFPT_Rho::compute_drhoc_g) and the ground-state XC potential on
+    /// the rho grid (for the diagonal term). Null/empty providers leave
+    /// the nlcc terms inert.
+    void set_rho(const DFPT_Rho* rho);
+    void set_vxc_r(const std::vector<double>* vxc_r);
+
+    /// per-q table of dvxc^a = f_xc . (drho_val^a + drhoc^a) for every
+    /// displacement slot [3*nat]: the first-order XC potential of the
+    /// TOTAL response (valence + core) supplied by the solver's
+    /// accumulation loop; coreless atoms carry their valence-only dvxc
+    /// (the kernel term contracts it against the stashed valence
+    /// responses of the mixed piece). Set by the solver once per q,
+    /// before the accumulation pass; owned by the caller.
+    void set_dvxc_core_table(const std::vector<std::vector<std::complex<double>>>* table);
+
     void assemble(int q_idx, DFPT_PW_Data& data);
 
     /// Fill the D[b][*] row of the electronic dynamical-matrix contribution
     /// for the converged displacement (atom_idx, dir); psi/wg are the
-    /// ground-state wavefunctions and occupations. Requires a wired
+    /// ground-state wavefunctions and occupations. dvxc_r is the first-order
+    /// XC potential of this displacement's TOTAL density response (valence
+    /// + core), computed by the caller from the analytic kernel; it gates
+    /// the nlcc kernel term (empty = no nlcc). Requires a wired
     /// DFPT_Pert (init); a null pert leaves the row untouched.
     void accumulate_electron(int q_idx,
                              int atom_idx,
                              int dir,
                              const psi::Psi<std::complex<double>>& psi,
                              const ModuleBase::matrix& wg,
-                             DFPT_PW_Data& data);
+                             DFPT_PW_Data& data,
+                             const std::vector<std::complex<double>>& dvxc_r);
 
     void diagonalize(int q_idx, DFPT_PW_Data& data);
 
@@ -138,6 +178,9 @@ class DFPT_Phon
     UnitCell* ucell_ = nullptr;
     ModulePW::PW_Basis* pw_rho_ = nullptr;
     DFPT_Pert* pert_ = nullptr;
+    const DFPT_Rho* rho_ = nullptr;
+    const std::vector<double>* vxc_r_ = nullptr;
+    const std::vector<std::vector<std::complex<double>>>* dvxc_core_ = nullptr;
 
     double ewald_alpha_ = 0.0;
     double ewald_rcut_ = 0.0;
@@ -153,6 +196,23 @@ class DFPT_Phon
                             const psi::Psi<std::complex<double>>& psi,
                             const ModuleBase::matrix& wg);
 
+    /// nlcc kernel term (QE dynmat_nlcc): for the converged displacement
+    /// (atom_idx, dir) with first-order XC potential dvxc_r, contract with
+    /// the core response of every displacement (iat2, idir2) and fill the
+    /// Hermitian pair (row, col) / (col, row). Gated on dvxc_r and the
+    /// providers set by set_rho.
+    void accum_nlcc_kernel(int atom_idx,
+                           int dir,
+                           int q_idx,
+                           DFPT_PW_Data& data,
+                           const std::vector<std::complex<double>>& dvxc_r);
+
+    /// nlcc diagonal term (QE dynmatcc): q-independent same-atom block
+    /// -Omega tpiba^2 sum_G G_i G_j rhoc~(|G|) exp(-i G tau_a)
+    /// conj(vxc^0(G)) added into the (3a+i, 3a+j) entries. Gated on the
+    /// ground-state XC potential provider (set_vxc_r).
+    void accum_nlcc_diag(int atom_idx, int q_idx);
+
     /// DFT+U contribution to the dynamical matrix (U0 reservation).
     void dftu_onsite(int q_idx, DFPT_PW_Data& data);