diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index c22017f783b..e887601e8a9 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/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 vxc_r; + if (nlcc_any) + { + const std::tuple 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) diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index de9fa67306c..4b04c7e24c2 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -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* vxc_r) +{ + vxc_r_ = vxc_r; +} + +void DFPT_Phon::set_dvxc_core_table(const std::vector>>* table) +{ + dvxc_core_ = table; +} + // --------------------------------------------------------------------------- // assemble / diagonalize / LO-TO / sum rule // --------------------------------------------------------------------------- diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 9bcd8dbaec8..7e460315ac1 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -5,7 +5,9 @@ #include "source_cell/unitcell.h" #include "source_psi/psi.h" +#include #include +#include 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* 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>>* 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>& psi, const ModuleBase::matrix& wg, - DFPT_PW_Data& data); + DFPT_PW_Data& data, + const std::vector>& 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* vxc_r_ = nullptr; + const std::vector>>* dvxc_core_ = nullptr; double ewald_alpha_ = 0.0; double ewald_rcut_ = 0.0; @@ -153,6 +196,23 @@ class DFPT_Phon const psi::Psi>& 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>& 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); diff --git a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp index 2cbef586ec7..38b86915ef4 100644 --- a/source/source_pw/module_dfpt/dfpt_phon_elec.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon_elec.cpp @@ -8,6 +8,10 @@ #include "dfpt_kq_basis.h" #include "dfpt_pert.h" #include "dfpt_pw_data.h" +#include "dfpt_rho.h" +#include "source_base/math_integral.h" +#include "source_base/math_sphbes.h" +#include "source_base/parallel_reduce.h" #include #include @@ -134,11 +138,12 @@ void scatter_chi_to_r(const ModulePW::PW_Basis& pw_rho, } // namespace void DFPT_Phon::accumulate_electron(int q_idx, - int atom_idx, - int dir, - const psi::Psi>& psi, - const ModuleBase::matrix& wg, - DFPT_PW_Data& data) + int atom_idx, + int dir, + const psi::Psi>& psi, + const ModuleBase::matrix& wg, + DFPT_PW_Data& data, + const std::vector>& dvxc_r) { ModuleBase::TITLE("DFPT_Phon", "accumulate_electron"); ModuleBase::timer::start("DFPT_Phon", "accumulate_electron"); @@ -186,11 +191,234 @@ void DFPT_Phon::accumulate_electron(int q_idx, } } + // ---- nlcc direct terms (QE dynmat_nlcc + dynmatcc) ---- + accum_nlcc_kernel(atom_idx, dir, q_idx, data, dvxc_r); + if (dir == 0) + { + // the diagonal block is direction-pair complete: add it once per + // atom (every accumulate_electron call would triple it) + accum_nlcc_diag(atom_idx, q_idx); + } + // restore the converged dpsi of this displacement restore_dpsi(data, q_idx, psi.get_nk(), psi.get_nbands(), dpsib); ModuleBase::timer::end("DFPT_Phon", "accumulate_electron"); } +void DFPT_Phon::accum_nlcc_kernel(int atom_idx, + int dir, + int q_idx, + DFPT_PW_Data& data, + const std::vector>& dvxc_r) +{ + ModuleBase::TITLE("DFPT_Phon", "accum_nlcc_kernel"); + ModuleBase::timer::start("DFPT_Phon", "accum_nlcc_kernel"); + const int nrxx = pw_rho_->nrxx; + if (rho_ == nullptr || ucell_ == nullptr || static_cast(dvxc_r.size()) != nrxx) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); + return; + } + // any nlcc atom at all? + bool nlcc_any = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any = true; + break; + } + } + if (!nlcc_any) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); + return; + } + const int nat = ucell_->nat; + const ModuleBase::Vector3 q_cart = data.get_qvec(q_idx) * ucell_->G; + // pool-global normalization: nrxx is the rank-local FFT slab size, + // so the raw sum needs a pool reduction (the mixing PR's convention) + const double fac = ucell_->omega / static_cast(pw_rho_->nxyz); + const int rowb = 3 * atom_idx + dir; + const double mass_row = ucell_->atoms[ucell_->iat2it[atom_idx]].mass; + // QE dynmat_nlcc semantics, pinned by the element-level comparison + // against QE dyn.G (bn_cubic): the raw matrix carries + // X_ab = (Omega/nrxx) sum_r dvxc^a(r) conj(drhoc^b(r)), + // dvxc^a = f_xc.(drho_val^a + drhoc^a), once per element, and the + // effective Hermitian value is the symmetrized + // S_ab = (1/2)(X_ab + conj(X_ba)). + // X_ab is two-sided only when BOTH displacements carry core response; + // for the one-sided cross blocks (e.g. B vs N: X_BN = 0 because N has + // no nlcc) this reduces to (1/2) X_NB. The column loop must therefore + // cover non-nlcc atoms too: their X_ab vanishes but the transposed + // piece conj(X_ba) = sum_r dvxc^b conj(drhoc^a) -- the column's + // valence-only dvaux against the ROW's core -- is exactly what the + // (nlcc-row, non-nlcc-column) element needs. + std::vector> drhoc_g; + std::vector> drhoc_r(nrxx); + // row's own core response in real space (for the conj(X_ba) piece; + // empty when the row atom carries no nlcc) + std::vector> drhoc_a_r; + rho_->compute_drhoc_g(*ucell_, atom_idx, dir, q_cart, drhoc_g); + const bool row_has_core = static_cast(drhoc_g.size()) == pw_rho_->npw; + if (row_has_core) + { + drhoc_a_r.resize(nrxx); + pw_rho_->recip2real(drhoc_g.data(), drhoc_a_r.data()); + } + const bool have_table = dvxc_core_ != nullptr; + for (int iat2 = 0; iat2 < nat; ++iat2) + { + const bool col_has_core = ucell_->atoms[ucell_->iat2it[iat2]].ncpp.nlcc; + for (int idir2 = 0; idir2 < 3; ++idir2) + { + // X_ab: the row's dvaux against the column's core response. + // Conjugation convention: the ROW response carries the + // conjugate (sum(conj(row)*col), matching the wavefunction + // and Ewald terms; at Gamma with real amplitudes this is a + // no-op, at nonzero q it fixes the imaginary cross terms) + std::complex x_ab(0.0, 0.0); + if (col_has_core) + { + rho_->compute_drhoc_g(*ucell_, iat2, idir2, q_cart, drhoc_g); + if (static_cast(drhoc_g.size()) != pw_rho_->npw) + { + continue; + } + pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + x_ab += std::conj(dvxc_r[ir]) * drhoc_r[ir]; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(x_ab); +#endif + x_ab *= fac; + } + // X_ba: the column's dvaux (per-q table) against the row's + // core response; same conjugation swap + std::complex x_ba(0.0, 0.0); + const int cola = 3 * iat2 + idir2; + if (row_has_core && have_table && cola < static_cast(dvxc_core_->size()) + && static_cast((*dvxc_core_)[cola].size()) == nrxx) + { + const std::vector>& dvxc_b = (*dvxc_core_)[cola]; + for (int ir = 0; ir < nrxx; ++ir) + { + x_ba += std::conj(dvxc_b[ir]) * drhoc_a_r[ir]; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(x_ba); +#endif + x_ba *= fac; + } + if (x_ab == std::complex(0.0, 0.0) && x_ba == std::complex(0.0, 0.0)) + { + continue; + } + const std::complex s_ab = 0.5 * (x_ab + std::conj(x_ba)); + const double mass_norm = std::sqrt(mass_row * ucell_->atoms[ucell_->iat2it[iat2]].mass); + // single-ordering fill: the transposed element receives its + // own S_ba = conj(S_ab) from the column displacement's pass + dynmat_accum_(rowb, cola) += s_ab / mass_norm; + } + } + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_kernel"); +} + +void DFPT_Phon::accum_nlcc_diag(int atom_idx, int q_idx) +{ + ModuleBase::TITLE("DFPT_Phon", "accum_nlcc_diag"); + ModuleBase::timer::start("DFPT_Phon", "accum_nlcc_diag"); + const int npw = pw_rho_ ? pw_rho_->npw : 0; + const int nrxx = pw_rho_ ? pw_rho_->nrxx : 0; + if (ucell_ == nullptr || pw_rho_ == nullptr || vxc_r_ == nullptr + || static_cast(vxc_r_->size()) != nrxx || npw == 0 || accum_q_ != q_idx) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); + return; + } + const int it = ucell_->iat2it[atom_idx]; + if (!ucell_->atoms[it].ncpp.nlcc) + { + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); + return; + } + // QE dynmatcc ground truth, evaluated as a real-space integral: + // D_ij += Int dr d2rho_c/(d tau_i d tau_j) (r) * v_xc^0(r) + // with the second-order shifted-core amplitude + // d2rho_c(G) = -tpiba^2 G_i G_j rhoc~(|G|) exp(-i G tau_a) + // (a real field: W(-G) = conj(W(G))), so the contraction runs through + // the same recip2real machinery as every other real-space field and + // never touches the half-space G-sum convention directly. + const pseudo& p = ucell_->atoms[it].ncpp; + const int mesh = p.msh; + const ModuleBase::Vector3 tau_cart + = ucell_->atoms[it].tau[ucell_->iat2ia[atom_idx]] * ucell_->lat0; + std::vector aux(mesh); + std::vector> w_g(npw); + std::vector> w_r(nrxx); + double d3[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; + const double fac = ucell_->omega / static_cast(pw_rho_->nxyz); + for (int i = 0; i < 3; ++i) + { + for (int j = i; j < 3; ++j) + { + std::fill(w_g.begin(), w_g.end(), std::complex(0.0, 0.0)); + for (int ig = 0; ig < npw; ++ig) + { + const ModuleBase::Vector3& gv = pw_rho_->gcar[ig]; + const double g2 = gv * gv; + // the G = 0 shell contributes nothing (G_i G_j factor) + if (g2 < 1.0e-12) + { + continue; + } + const double gphys = std::sqrt(g2) * ucell_->tpiba; + ModuleBase::Sphbes::Spherical_Bessel(mesh, p.r.data(), gphys, 0, aux.data()); + double rhocg = 0.0; + for (int ir = 0; ir < mesh; ++ir) + { + aux[ir] = p.r[ir] * p.r[ir] * p.rho_atc[ir] * aux[ir]; + } + ModuleBase::Integral::Simpson_Integral(mesh, aux.data(), p.rab.data(), rhocg); + rhocg *= ModuleBase::FOUR_PI / ucell_->omega; + + const double phase = ModuleBase::TWO_PI * (gv * tau_cart) / ucell_->lat0; + // AGENTS rule 14: named locals instead of computed ctor args + const double gtau_re = std::cos(phase); + const double gtau_im = -std::sin(phase); + const std::complex gtau(gtau_re, gtau_im); + w_g[ig] = -ucell_->tpiba * ucell_->tpiba * gv[i] * gv[j] * rhocg * gtau; + } + pw_rho_->recip2real(w_g.data(), w_r.data()); + double s = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + s += w_r[ir].real() * (*vxc_r_)[ir]; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(s); +#endif + d3[i][j] = fac * s; + if (j != i) + { + d3[j][i] = d3[i][j]; + } + } + } + const double mass = ucell_->atoms[it].mass; + const int base = 3 * atom_idx; + for (int i = 0; i < 3; ++i) + { + for (int j = 0; j < 3; ++j) + { + dynmat_accum_(base + i, base + j) += d3[i][j] / mass; + } + } + ModuleBase::timer::end("DFPT_Phon", "accum_nlcc_diag"); +} + void DFPT_Phon::accum_d2_same_atom(int atom_idx, int dir, int iat, diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 69ec3892d10..0b014cc8b22 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -43,6 +43,11 @@ struct DFPT_PW_InitContext double nelec; double ecutwfc; const Plus_U_Base* dftu; + /// ground-state XC potential on the rho grid (valence + core charge + /// input, the same density the kernel linearizes around); consumed by + /// the nlcc diagonal 2n+1 term (QE dynmatcc). Null disables that term. + /// No default member initializer: C++11 aggregate rules require it. + const std::vector* vxc_r; }; /** @@ -91,7 +96,7 @@ class DFPT_PW double ecutwfc, const Plus_U_Base* dftu) { - const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu}; + const DFPT_PW_InitContext ctx{&ucell, &psi, pw_rho, pw_wfc, sf, &veff_r, &wg, &eig, xc, nelec, ecutwfc, dftu, nullptr}; init(ctx); } diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp index 529517c6212..d62998185ff 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp @@ -115,6 +115,43 @@ std::vector> DFPT_PW_Data::get_vsc_r(int atom_idx, int dir) return std::vector>(); } +void DFPT_PW_Data::set_drho_disp(int atom_idx, int dir, const std::vector>& rho) +{ + ModuleBase::TITLE("DFPT_PW_Data", "set_drho_disp"); + ModuleBase::timer::start("DFPT_PW_Data", "set_drho_disp"); + if (atom_idx < 0 || dir < 0 || dir >= 3) + { + ModuleBase::timer::end("DFPT_PW_Data", "set_drho_disp"); + return; + } + const size_t slot = static_cast(3 * atom_idx + dir); + if (slot >= drho_disp_.size()) + { + drho_disp_.resize(slot + 1); + } + drho_disp_[slot] = rho; + ModuleBase::timer::end("DFPT_PW_Data", "set_drho_disp"); +} + +std::vector> DFPT_PW_Data::get_drho_disp(int atom_idx, int dir) const +{ + ModuleBase::TITLE("DFPT_PW_Data", "get_drho_disp"); + ModuleBase::timer::start("DFPT_PW_Data", "get_drho_disp"); + if (atom_idx < 0 || dir < 0 || dir >= 3) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return std::vector>(); + } + const size_t slot = static_cast(3 * atom_idx + dir); + if (slot < drho_disp_.size()) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return drho_disp_[slot]; + } + ModuleBase::timer::end("DFPT_PW_Data", "get_drho_disp"); + return std::vector>(); +} + void DFPT_PW_Data::set_dpsi_disp(int atom_idx, int dir, const std::vector>>>& d) diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h index 239c77fa4ad..79fdc824d36 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.h +++ b/source/source_pw/module_dfpt/dfpt_pw_data.h @@ -235,6 +235,13 @@ class DFPT_PW_Data void set_vsc_r(int atom_idx, int dir, const std::vector>& v); std::vector> get_vsc_r(int atom_idx, int dir) const; + /// converged first-order density response (G space) of displacement + /// (atom, dir): the two-pass 2n+1 accumulation needs it to rebuild the + /// nlcc kernel input f_xc.(drho_val + drhoc) per displacement, after + /// the working drho slots hold the last solved displacement + void set_drho_disp(int atom_idx, int dir, const std::vector>& rho); + std::vector> get_drho_disp(int atom_idx, int dir) const; + /// converged dpsi of displacement (atom, dir), indexed [k][band]; the /// two-pass 2n+1 accumulation reads it back after all displacements of /// the basis have been solved (the working dpsi slots get overwritten by @@ -302,6 +309,9 @@ class DFPT_PW_Data /// converged v_sc per displacement (atom, dir): [3*nat] entries std::vector>> vsc_r_; + /// converged first-order density per displacement (atom, dir): [3*nat] + std::vector>> drho_disp_; + /// converged dpsi per displacement (atom, dir): [3*nat][k][band] entries std::vector>>>> dpsi_disp_; diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index f7ec543731d..9bb03fbf1fc 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -67,7 +67,8 @@ class DFPT_PW::Impl /// v_sc^q = v_Hartree(drho_in) + xc(drho_in) on the shared grid; /// shared by solve_displacement and solve_efield_resp so the /// real-space Hartree+XC add loop is not duplicated - void assemble_v_sc(const ModuleBase::Vector3& q_cart, + void assemble_v_sc(int q_idx, + const ModuleBase::Vector3& q_cart, const std::vector>& drho_in_g, std::vector>& v_sc_r) const; @@ -159,6 +160,10 @@ class DFPT_PW::Impl ModulePW::PW_Basis_K* pw_wfc_ = nullptr; Structure_Factor* sf_ = nullptr; std::vector veff_r_; + /// ground-state XC potential (valence + core input), copied from the + /// init context: the nlcc diagonal 2n+1 term consumes it; the copy + /// owns the storage because the context points at caller stack state + std::vector vxc_r_; ModuleBase::matrix wg_; ModuleBase::matrix eig_; const XC_First_Order* xc_ = nullptr; diff --git a/source/source_pw/module_dfpt/dfpt_pw_init.cpp b/source/source_pw/module_dfpt/dfpt_pw_init.cpp index 47240b6618c..179c001a665 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_init.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_init.cpp @@ -112,6 +112,15 @@ void DFPT_PW::Impl::init_submodules(const DFPT_PW_InitContext& ctx, } rho_.init({nspin, nrxx, ctx.pw_rho, ctx.pw_wfc, ctx.ucell->G, rho_mix_type, mix_beta, kerker_a2, broyden_ndim_}); phon_.init(*ctx.ucell, ctx.pw_rho, &pert_); + // nlcc direct-term providers (QE dynmat_nlcc + dynmatcc equivalents): + // the pure core-response builder and the ground-state XC potential + phon_.set_rho(&rho_); + // the Impl-owned copy: the context pointer refers to caller stack + // state that dies when the esolver init returns + // always pass the pointer (never nullptr): a rank with an empty + // local field must still enter accum_nlcc_diag's collective + // participation; the size checks inside gate only the local work + phon_.set_vxc_r(&vxc_r_); q0_.init(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, &pert_); hamilt_.reset(new DFPT_HamiltShift(*ctx.ucell, ctx.pw_rho, ctx.pw_wfc, *ctx.veff_r, &pert_)); data_.init(&qlist_, nk, nbands, npw_max, nrxx, nspin, nat, ctx.dftu); @@ -137,6 +146,11 @@ void DFPT_PW::init(const DFPT_PW_InitContext& ctx) pimpl_->pw_wfc_ = pw_wfc; pimpl_->sf_ = sf; pimpl_->veff_r_ = *ctx.veff_r; + pimpl_->vxc_r_.clear(); + if (ctx.vxc_r != nullptr) + { + pimpl_->vxc_r_ = *ctx.vxc_r; + } pimpl_->wg_ = *ctx.wg; pimpl_->eig_ = *ctx.eig; pimpl_->xc_ = ctx.xc; diff --git a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp index 16758a7231f..5238592ae12 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_q0.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_q0.cpp @@ -323,7 +323,7 @@ void DFPT_PW::Impl::solve_efield_resp(int q_idx) // (shared Impl::assemble_v_sc from dfpt_pw_solve.cpp) std::vector> v_sc_r; const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_in_g, v_sc_r); + assemble_v_sc(q_idx, q_cart, drho_in_g, v_sc_r); for (int ik = 0; ik < nk; ++ik) { if (static_cast(occ_kq_.size()) <= ik || occ_kq_[ik].empty()) diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 9db6f17e942..3020cda2dc3 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -79,12 +79,87 @@ double DFPT_PW::Impl::run_displacement_irrep_pass(int q_idx, int irrep) worst = std::max(worst, residual); } } + // nlcc: the 2n+1 kernel term needs the first-order XC potential of + // each displacement's TOTAL response (valence + core). Precompute the + // nlcc flag once per pass; the per-displacement core response is + // recomputed here because the working drhoc slot holds the LAST + // solved displacement, and the valence response comes from the + // per-displacement stash. + bool nlcc_any = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any = true; + break; + } + } + const int npw = pw_rho_->npw; + const int nrxx = pw_rho_->nrxx; + const ModuleBase::Vector3 q_cart = data_.get_qvec(q_idx) * ucell_->G; + // per-q table of f_xc.(drho_val + drhoc) for every displacement slot: + // the total-response first-order XC potential that the Hermitized + // kernel term contracts against the core responses (QE dynmat_nlcc + // semantics; built once, before the accumulation pass) + std::vector>> dvxc_core; + if (nlcc_any && xc_ != nullptr) + { + // fill EVERY displacement slot: the kernel term X(i,j) contracts + // the row's dvaux with the COLUMN's core response, so rows of + // non-nlcc atoms need their dvaux too (f_xc applied to their + // valence response alone -- they have no core piece). Leaving + // their slots empty silently dropped the entire X(N,B) cross + // block (the dominant mixed piece for a B-N cell). + dvxc_core.resize(3 * nat); + std::vector> drhoc_g; + for (int iat = 0; iat < nat; ++iat) + { + const bool iat_nlcc = ucell_->atoms[ucell_->iat2it[iat]].ncpp.nlcc; + for (int idir = 0; idir < 3; ++idir) + { + std::vector> a_r(nrxx, std::complex(0.0, 0.0)); + if (iat_nlcc) + { + rho_.compute_drhoc_g(*ucell_, iat, idir, q_cart, drhoc_g); + if (static_cast(drhoc_g.size()) == npw) + { + pw_rho_->recip2real(drhoc_g.data(), a_r.data()); + } + } + // add the stashed converged valence response of this + // displacement (the working slots hold the last one) + const std::vector> drho_a_g = data_.get_drho_disp(iat, idir); + if (static_cast(drho_a_g.size()) == npw) + { + std::vector> dval_r(nrxx); + pw_rho_->recip2real(drho_a_g.data(), dval_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + a_r[ir] += dval_r[ir]; + } + } + xc_->set_context(q_cart); + xc_->apply(a_r, dvxc_core[3 * iat + idir]); + } + } + phon_.set_dvxc_core_table(&dvxc_core); + } for (int iat = 0; iat < nat; ++iat) { for (int idir = 0; idir < 3; ++idir) { + std::vector> dvxc_r; + if (nlcc_any && xc_ != nullptr) + { + const int slot = 3 * iat + idir; + if (slot < static_cast(dvxc_core.size()) + && static_cast(dvxc_core[slot].size()) == nrxx) + { + dvxc_r = dvxc_core[slot]; + } + } // 2n+1 accumulation of this converged displacement - phon_.accumulate_electron(q_idx, iat, idir, gs_psi_, wg_, data_); + phon_.accumulate_electron(q_idx, iat, idir, gs_psi_, wg_, data_, dvxc_r); } } data_.add_residual(q_idx, irrep, worst); diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index 876519826f7..1a30c81fe88 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -48,7 +48,8 @@ double safe_relative_defect(double dn2, double o2) } } // namespace -void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, +void DFPT_PW::Impl::assemble_v_sc(int q_idx, + const ModuleBase::Vector3& q_cart, const std::vector>& drho_in_g, std::vector>& v_sc_r) const { @@ -82,8 +83,10 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, pw_rho_->recip2real(drho_in_g.data(), a_r.data()); // nlcc: the first-order core density enters the XC channel only (the // core Coulomb is already in the ionic pseudopotential, so the Hartree - // above correctly uses the valence density alone) - const std::vector> drhoc_g = data_.get_drhoc_g(last_q_); + // above correctly uses the valence density alone). The q index is + // passed explicitly: the drhoc slot is per q and last_q_ can be stale + // on the first assembly of a new q. + const std::vector> drhoc_g = data_.get_drhoc_g(q_idx); if (static_cast(drhoc_g.size()) == pw_rho_->npw) { std::vector> drhoc_r(nrxx); @@ -161,6 +164,25 @@ void DFPT_PW::Impl::stash_converged_disp_response(int q_idx, int nbands) { data_.set_vsc_r(iat, idir, v_sc_r_last); + // stash the converged density response too: the two-pass accumulation + // rebuilds f_xc.(drho_val + drhoc) per displacement after the working + // drho slots have moved on; only needed when some species carries + // a nonlinear core correction (skip the 3*nat*npw complex copies + // for coreless runs) + bool nlcc_any_solve = false; + for (int it = 0; it < ucell_->ntype; ++it) + { + if (ucell_->atoms[it].ncpp.nlcc) + { + nlcc_any_solve = true; + break; + } + } + if (nlcc_any_solve) + { + const std::vector> drho_conv = data_.get_drho_g(q_idx, 0); + data_.set_drho_disp(iat, idir, drho_conv); + } std::vector>>> disp( nk, std::vector>>(nbands)); @@ -290,7 +312,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) else { const std::vector> drho_in_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_in_g, v_sc_r); + assemble_v_sc(q_idx, q_cart, drho_in_g, v_sc_r); } if (dbg) { @@ -335,7 +357,7 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) // potential self-consistency std::vector> v_sc_out_r; const std::vector> drho_out_g = data_.get_drho_g(q_idx, 0); - assemble_v_sc(q_cart, drho_out_g, v_sc_out_r); + assemble_v_sc(q_idx, q_cart, drho_out_g, v_sc_out_r); if (vsc_broyden) { // Broyden on the real-space potential: push, cal_coef, mix diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 326c4c7f04f..5dd8abe3349 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -289,47 +289,59 @@ void DFPT_Rho::compute_drhoc(const UnitCell& ucell, { ModuleBase::TITLE("DFPT_Rho", "compute_drhoc"); ModuleBase::timer::start("DFPT_Rho", "compute_drhoc"); + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; + std::vector> drhoc_g; + compute_drhoc_g(ucell, iat, dir, q_cart, drhoc_g); + data.set_drhoc_g(q_idx, drhoc_g); + + // early return for invalid/coreless atoms: compute_drhoc_g already + // returned empty, so the iat2it indexing below would be out of + // bounds for a bad index — skip the diagnostic in that case + if (drhoc_g.empty() || iat < 0 || iat >= ucell.nat) + { + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); + return; + } + const int it = ucell.iat2it[iat]; + // report magnitude for diagnostics + double norm2 = 0.0; + for (const std::complex& c : drhoc_g) + { + norm2 += std::norm(c); + } + std::cout << " DFPT nlcc core derivative: atom " << iat + 1 << "(" + << ucell.atoms[it].label << ") dir " << static_cast('x' + dir) + << " |drhoc_g| = " << std::setprecision(3) << std::scientific << std::sqrt(norm2) + << std::defaultfloat << std::setprecision(6) << std::endl; + + ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); +} + +void DFPT_Rho::compute_drhoc_g(const UnitCell& ucell, + int iat, + int dir, + const ModuleBase::Vector3& q_cart, + std::vector>& drhoc_g) const +{ const int npw = pw_rho_ ? pw_rho_->npw : 0; + drhoc_g.clear(); if (pw_rho_ == nullptr || npw == 0 || iat < 0 || iat >= ucell.nat || dir < 0 || dir > 2) { - data.set_drhoc_g(q_idx, {}); - ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); return; } const int it = ucell.iat2it[iat]; if (!ucell.atoms[it].ncpp.nlcc) { - data.set_drhoc_g(q_idx, {}); - ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); return; } - // q-shifted |G+q|^2 in tpiba^2 units, and the q-shifted structure - // factor phase exp(i (G+q) tau_iat) for the displaced atom - const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); - const ModuleBase::Vector3 q_cart = q_frac * recip_matrix_; // tau in cartesian bohr (lat0-scaled: tau is in direct units of latvec) const ModuleBase::Vector3 tau_cart = ucell.atoms[it].tau[ucell.iat2ia[iat]] * ucell.lat0; - - // build the |G+q|^2 shell values and evaluate the radial core FT - // on them (non_linear_core_correction works on unique gg shells) - std::vector w2_all(npw); - for (int ig = 0; ig < npw; ++ig) - { - const ModuleBase::Vector3 w = pw_rho_->gcar[ig] + q_cart; - w2_all[ig] = w * w; // dimensionless (1/lat0)^2 - } - - // use the existing radial FT machinery: pass the |G+q|^2 values as - // pseudo-unique shells (the function just needs matching arrays; - // we can call it per-shell or build a shell lookup) - // Simpler: directly evaluate the spherical Bessel integral here - // for each ig (no shell deduplication -- the cost is O(npw * mesh) - // and only runs once per displacement for nlcc atoms) const pseudo& p = ucell.atoms[it].ncpp; const int mesh = p.msh; - std::vector> drhoc_g(npw, std::complex(0.0, 0.0)); + drhoc_g.assign(npw, std::complex(0.0, 0.0)); std::vector aux(mesh); for (int ig = 0; ig < npw; ++ig) @@ -352,30 +364,21 @@ void DFPT_Rho::compute_drhoc(const UnitCell& ucell, ModuleBase::Integral::Simpson_Integral(mesh, aux.data(), p.rab.data(), rhocg); rhocg *= ModuleBase::FOUR_PI / ucell.omega; - // structure factor phase: exp(i (G+q) . tau_iat) + // structure factor phase: exp(-i (G+q) . tau_iat). The sign is + // load-bearing: exp(+i(G+q)tau) mirrors the core response to the + // -tau site (for wurtzite B the mirror lands exactly on the N + // sites of the other layer), which bn_cubic cannot expose (B at + // the origin makes both conventions identical) but which showed + // up element-by-element in the bn_hex matrix comparison against + // QE dyn.G. QE addcore: e^{-i G tau} * e^{-i q tau}. const double phase = ModuleBase::TWO_PI * (w * tau_cart) / ucell.lat0; - const std::complex gtau(std::cos(phase), std::sin(phase)); + const std::complex gtau(std::cos(phase), -std::sin(phase)); // direction derivative: -i * tpiba * (G+q)_dir const double wdir = ucell.tpiba * w[dir]; drhoc_g[ig] = std::complex(0.0, -wdir) * rhocg * gtau; } - - data.set_drhoc_g(q_idx, drhoc_g); - - // report magnitude for diagnostics - double norm2 = 0.0; - for (int ig = 0; ig < npw; ++ig) - { - norm2 += std::norm(drhoc_g[ig]); - } - std::cout << " DFPT nlcc core derivative: atom " << iat + 1 << "(" - << ucell.atoms[it].label << ") dir " << static_cast('x' + dir) - << " |drhoc_g| = " << std::setprecision(3) << std::scientific << std::sqrt(norm2) - << std::defaultfloat << std::setprecision(6) << std::endl; - - ModuleBase::timer::end("DFPT_Rho", "compute_drhoc"); } void DFPT_Rho::cal_docc(const psi::Psi>& psi, diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 89cb83a95a4..f1083693c38 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -134,6 +134,17 @@ class DFPT_Rho int q_idx, DFPT_PW_Data& data); + /// same first-order core-density amplitude as compute_drhoc, but pure: + /// takes the cartesian q explicitly and writes into a caller-owned + /// vector instead of the DFPT_PW_Data slot. Used by the 2n+1 nlcc + /// direct terms, which need the core response of EVERY displacement + /// while the data slot holds only the current one. + void compute_drhoc_g(const UnitCell& ucell, + int iat, + int dir, + const ModuleBase::Vector3& q_cart, + std::vector>& drhoc_g) const; + /// first-order occupation matrix (docc) for DFT+U (U0 reservation). void cal_docc(const psi::Psi>& psi, const ModuleBase::matrix& wg, diff --git a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt index 6e4171510e5..d47a1de06e1 100644 --- a/source/source_pw/module_dfpt/test_serial/CMakeLists.txt +++ b/source/source_pw/module_dfpt/test_serial/CMakeLists.txt @@ -88,6 +88,7 @@ addtest( ../dfpt_phon.cpp ../dfpt_phon_ewald.cpp ../dfpt_phon_elec.cpp + ../dfpt_rho.cpp ../dfpt_pert.cpp ../dfpt_pert_vkb.cpp ../dfpt_pert_nl.cpp diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp index aadaa8e91e0..d0874671e66 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_phon_serial_test.cpp @@ -12,6 +12,7 @@ #include "dfpt_serial_fixture.h" #include "source_base/complexmatrix.h" +#include "source_pw/module_dfpt/dfpt_rho.h" #include "source_base/constants.h" #include "source_base/matrix3.h" #include "source_base/vector3.h" @@ -347,7 +348,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronAnalyticContraction) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 1, psi, wg, data_); + phon_.accumulate_electron(0, 0, 1, psi, wg, data_, std::vector>()); // expected: row 1 (atom 0, dir 1) of the Hermitian 2n+1 accumulation: // the row element receives wg* once per off-diagonal @@ -400,7 +401,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronD2GateOffGenericQ) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 0, psi, wg, data_); + phon_.accumulate_electron(0, 0, 0, psi, wg, data_, std::vector>()); const std::vector> g0(1, ModuleBase::Vector3(0.0, 0.0, 0.0)); for (int adir = 0; adir < 3; ++adir) @@ -479,7 +480,7 @@ TEST_F(DFPTPhonSerialTest, AccumulateElectronD2CommensurateQ) } data_.set_dpsi(0, 0, 0, dpsi_inj); - phon_.accumulate_electron(0, 0, 1, psi, wg, data_); + phon_.accumulate_electron(0, 0, 1, psi, wg, data_, std::vector>()); for (int adir = 0; adir < 3; ++adir) { @@ -964,3 +965,98 @@ TEST_F(DFPTPhonSerialTest, AsrSimpleCrossCartesianRegression) EXPECT_EQ(out(i, i).imag(), 0.0) << "i=" << i; } } + +// --------------------------------------------------------------------------- +// nlcc direct-term fill-semantics regression: with a synthetic provider set +// (all-atoms kernel table, ground-state XC vector, DFPT_Rho), the nlcc +// kernel accumulation must satisfy the QE-pinned invariants: +// (1) the one-sided cross block (nlcc row, coreless column) receives +// (1/2)(X_ab + conj(X_ba)) = (1/2) X_ab [X_ba = 0: coreless row] +// (2) every displacement slot's dvaux participates (the "dead row" bug +// dropped the entire X(N, B) cross block when the table left +// coreless atoms empty) +// (3) the matrix stays Hermitian +// These pin the fill factor (the double-sided same-atom block would double +// under a both-order fill; a lower-triangle-only fill halves the one-sided +// cross blocks). +// --------------------------------------------------------------------------- +TEST_F(DFPTPhonSerialTest, NlccDirectTermsFillSemantics) +{ + MakeNlccTwoAtomCell(); + SetupPhon(k_d_, q_d_); + + ModuleDFPT::DFPT_Rho rho; + rho.init({1, pw_rho_.nrxx, &pw_rho_, &pw_wfc_, G_, std::string("rho_plain"), 0.4, 1.0, 8}); + phon_.set_rho(&rho); + + const int nrxx = pw_rho_.nrxx; + // synthetic kernel table: f_xc.drhoc stand-ins = distinct constant + // complex fields per slot (values chosen to survive exact comparison) + std::vector>> table(6); + for (int slot = 0; slot < 6; ++slot) + { + const double re = 0.1 + 0.05 * slot; + const double im = 0.02 * slot - 0.05; + table[slot].assign(nrxx, std::complex(re, im)); + } + phon_.set_dvxc_core_table(&table); + + // synthetic ground-state XC vector (only the diag term consumes it); + // zero it so the diag term's only effect is absent (this test targets + // the kernel fill) + const std::vector vxc_zero(nrxx, 0.0); + phon_.set_vxc_r(&vxc_zero); + + psi::Psi> psi = MakeSinglePlaneWavePsi(k_d_); + ModuleBase::matrix wg = MakeOccWeights(); + // NO dpsi injection: the bare cross term needs occupied dpsi entries + // and skips empty slots, so the accumulated matrix isolates the nlcc + // kernel contribution exactly + + phon_.accumulate_electron(0, 1, 0, psi, wg, data_, table[3]); + + const auto& m = phon_.dynmat_accum(); + // row (1, dir 0) = 3*1+0 = 3; column (atom 0, dir j) = j (coreless) + // X_3j = sum_r table[3](r) * conj(drhoc(atom0, dir_j)(r)) = 0 + // (atom 0 has no nlcc: drhoc = empty -> the kernel term's x_ab arm + // is zero for coreless COLUMNS) + // x_ba arm = sum_r table[j](r) * conj(drhoc(atom1, dir0)(r)) != 0 + // (row atom 1 has a core; the column's dvaux comes from the table) + // So the (3, j) element = (1/2)(0 + conj(x_ba)) = (1/2) conj(x_ba). + // The x_ba contraction runs over the actual Gaussian core field, so + // verify structurally: the one-sided block is NONZERO (dead-row bug + // gave exactly zero) and matches an independent recomputation. + const double fac = ucell_.omega / nrxx; + std::vector> drhoc_g; + rho.compute_drhoc_g(ucell_, 1, 0, q_cart_, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + std::vector> drhoc_r(nrxx); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + const double m0 = std::sqrt(ucell_.atoms[0].mass); + const double m1 = std::sqrt(ucell_.atoms[1].mass); + for (int j = 0; j < 3; ++j) + { + std::complex x_ba(0.0, 0.0); + for (int ir = 0; ir < nrxx; ++ir) + { + x_ba += std::conj(table[j][ir]) * drhoc_r[ir]; + } + x_ba *= fac; + const std::complex expect = 0.5 * std::conj(x_ba) / (m0 * m1); + EXPECT_NEAR(std::abs(m(3, j) - expect), 0.0, 1.0e-10 * (1.0 + std::abs(expect))) + << "j=" << j << " got " << m(3, j) << " expect " << expect; + } + + // Hermiticity of the accumulated block so far: accumulate the column + // displacement (atom 1, dir 0 already done as row; do (atom 1, dir 1) + // to complete more of the same-atom block, then check the whole + // matrix is Hermitian within the accumulated rows + phon_.accumulate_electron(0, 1, 1, psi, wg, data_, table[4]); + // the (3, 4) and (4, 3) elements must be conjugate pairs: from row 3's + // pass, (3,4) = (1/2)(X_34 + conj(X_43)); from row 4's pass, (4,3) = + // (1/2)(X_43 + conj(X_34)) = conj((3,4)) + const std::complex e34 = m(3, 4); + const std::complex e43 = m(4, 3); + EXPECT_NEAR(e43.real(), e34.real(), 1.0e-12); + EXPECT_NEAR(e43.imag(), -e34.imag(), 1.0e-12); +} diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp index 66c0543f4a6..3e7b268d4c7 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_rho_serial_test.cpp @@ -660,3 +660,221 @@ TEST_F(DFPTRhoSerialTest, BroydenResetMixingAllowsFreshQAndSameQCycle) EXPECT_NEAR(mixed_q0[ig].imag(), 0.4 * out[ig].imag(), 1.0e-12); } } + +// --------------------------------------------------------------------------- +// nlcc phase regression: the first-order core response must localize at the +// ATOM site, not its mirror. This is the anti-blind-spot test for the +// structure-factor sign: exp(+i(G+q)tau) mirrors the response to -tau, a bug +// every historical validation missed because the nlcc atom sat at the origin +// (bn_cubic) where both conventions coincide. The wurtzite bn_hex comparison +// against QE dyn.G pinned the sign: the mirror of B1 lands exactly on an N +// site of the other layer, collapsing the interlayer shear to -399 cm^-1. +// --------------------------------------------------------------------------- +TEST_F(DFPTRhoSerialTest, DrhocLocalizesAtTheAtomNotItsMirror) +{ + MakeNlccTwoAtomCell(); + SetupBases(k_d_, q_d_, 2); + ModuleDFPT::DFPT_Rho::Config cfg; + cfg.nspin = 1; + cfg.nrxx = pw_rho_.nrxx; + cfg.pw_rho = &pw_rho_; + cfg.pw_wfc = &pw_wfc_; + cfg.recip_matrix = G_; + cfg.mix_type = "rho_plain"; + cfg.mix_beta = 0.4; + cfg.kerker_a2 = 1.0; + cfg.broyden_ndim = 8; + ModuleDFPT::DFPT_Rho rho; + rho.init(cfg); + + // atom 1 carries the nlcc at the non-origin tau_nlcc_; all three + // displacement directions must localize at the atom + const int iat = 1; + const int it = ucell_.iat2it[iat]; + ASSERT_TRUE(ucell_.atoms[it].ncpp.nlcc); + // tau is stored in lat0-cartesian units (fixture convention: cell edge + // a_ in lat0 units); the true fractional position divides by a_ + const double inv_a = 1.0 / a_; + const ModuleBase::Vector3 tau_frac(ucell_.atoms[it].tau[ucell_.iat2ia[iat]].x * inv_a, + ucell_.atoms[it].tau[ucell_.iat2ia[iat]].y * inv_a, + ucell_.atoms[it].tau[ucell_.iat2ia[iat]].z * inv_a); + // mirror position folded back into the cell + ModuleBase::Vector3 mirror_frac(1.0 - tau_frac.x, 1.0 - tau_frac.y, 1.0 - tau_frac.z); + + const int nrxx = pw_rho_.nrxx; + std::vector> drhoc_g; + std::vector> drhoc_r(nrxx); + for (int dir = 0; dir < 3; ++dir) + { + // q = 0 first: the amplitude is the plain derivative profile, its + // |.|^2 centroid is at the atom (up to the derivative skew) + const ModuleBase::Vector3 q_zero(0.0, 0.0, 0.0); + rho.compute_drhoc_g(ucell_, iat, dir, q_zero, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + double wsum = 0.0; + double cx = 0.0; + double cy = 0.0; + double cz = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double wgt = std::norm(drhoc_r[ir]); + // fractional coordinates of the grid point via the FFT triple + const int iz = ir % pw_rho_.nz; + const int iy = (ir / pw_rho_.nz) % pw_rho_.ny; + const int ix = ir / (pw_rho_.ny * pw_rho_.nz); + // wrap-around: grid index -> fractional in [0, 1) + const double fx = static_cast(ix) / pw_rho_.nx; + const double fy = static_cast(iy) / pw_rho_.ny; + const double fz = static_cast(iz) / pw_rho_.nz; + // weight by |field|^2, folding the +1 wrap for the periodic image + // closest to tau: the field is localized, so the unwrapped + // centroid near tau is what we need + double wx = fx; + double wy = fy; + double wz = fz; + if (std::abs(fx - tau_frac.x) > 0.5) + { + wx += (fx < tau_frac.x) ? 1.0 : -1.0; + } + if (std::abs(fy - tau_frac.y) > 0.5) + { + wy += (fy < tau_frac.y) ? 1.0 : -1.0; + } + if (std::abs(fz - tau_frac.z) > 0.5) + { + wz += (fz < tau_frac.z) ? 1.0 : -1.0; + } + wsum += wgt; + cx += wgt * wx; + cy += wgt * wy; + cz += wgt * wz; + } + ASSERT_GT(wsum, 0.0); + const ModuleBase::Vector3 com(cx / wsum, cy / wsum, cz / wsum); + // minimum-image distance: wrap each component into [-0.5, 0.5) + auto mindist2 = [](const ModuleBase::Vector3& a, + const ModuleBase::Vector3& b) { + double dx = a.x - b.x; dx -= std::round(dx); + double dy = a.y - b.y; dy -= std::round(dy); + double dz = a.z - b.z; dz -= std::round(dz); + return dx * dx + dy * dy + dz * dz; + }; + const double dist_atom = mindist2(com, tau_frac); + const double dist_mirror = mindist2(com, mirror_frac); + // the centroid of a derivative field sits at the atom up to the + // derivative skew; the mirror is a full half-cell away in the + // minimum-image sense, far beyond any tolerance that would admit + // the reflected profile + EXPECT_LT(dist_atom, 0.01) << "dir=" << dir; + EXPECT_GT(dist_mirror, 25.0 * dist_atom) << "dir=" << dir; + } + + // non-zero q variant: the q-phase exp(-i q.tau) must not mirror either. + // Use a generic incommensurate q (the fixture default q_d_ is already + // generic); the real-space field of the q-shifted amplitude still + // localizes at the atom (the e^{iqr} carrier is smooth on the sigma + // scale of the core). + const ModuleBase::Vector3 q2 = q_d_ * 0.7 + ModuleBase::Vector3(0.11, -0.07, 0.05); + const ModuleBase::Vector3 q2_cart = q2 * ucell_.G; + rho.compute_drhoc_g(ucell_, iat, 1, q2_cart, drhoc_g); + ASSERT_EQ(static_cast(drhoc_g.size()), pw_rho_.npw); + pw_rho_.recip2real(drhoc_g.data(), drhoc_r.data()); + double wsum = 0.0; + double cx = 0.0; + double cy = 0.0; + double cz = 0.0; + for (int ir = 0; ir < nrxx; ++ir) + { + const double wgt = std::norm(drhoc_r[ir]); + const int iz = ir % pw_rho_.nz; + const int iy = (ir / pw_rho_.nz) % pw_rho_.ny; + const int ix = ir / (pw_rho_.ny * pw_rho_.nz); + double wx = static_cast(ix) / pw_rho_.nx; + double wy = static_cast(iy) / pw_rho_.ny; + double wz = static_cast(iz) / pw_rho_.nz; + if (std::abs(wx - tau_frac.x) > 0.5) + { + wx += (wx < tau_frac.x) ? 1.0 : -1.0; + } + if (std::abs(wy - tau_frac.y) > 0.5) + { + wy += (wy < tau_frac.y) ? 1.0 : -1.0; + } + if (std::abs(wz - tau_frac.z) > 0.5) + { + wz += (wz < tau_frac.z) ? 1.0 : -1.0; + } + wsum += wgt; + cx += wgt * wx; + cy += wgt * wy; + cz += wgt * wz; + } + ASSERT_GT(wsum, 0.0); + const ModuleBase::Vector3 com2(cx / wsum, cy / wsum, cz / wsum); + // same minimum-image comparison as above + auto mindist2 = [](const ModuleBase::Vector3& a, + const ModuleBase::Vector3& b) { + double dx = a.x - b.x; dx -= std::round(dx); + double dy = a.y - b.y; dy -= std::round(dy); + double dz = a.z - b.z; dz -= std::round(dz); + return dx * dx + dy * dy + dz * dz; + }; + const double d2a = mindist2(com2, tau_frac); + const double d2m = mindist2(com2, mirror_frac); + EXPECT_LT(d2a, 0.02); + EXPECT_GT(d2m, 25.0 * d2a); + + // phase-sensitive check at generic q: the centroid via |field|^2 is + // invariant under a global phase rotation, so flipping the sign of + // exp(i q.tau) alone would leave it unchanged. Instead, verify the + // G-space coefficient phase directly: drhoc_g(G) must carry + // exp(-i(G+q).tau) (the QE addcore convention), not the conjugate. + // Compare the q=0 and q=q2 coefficients: their ratio at each G must + // be exp(-i q.cart . tau) * (ratio of the radial/Bessel factor at + // |G| vs |G+q|). For a single G with a dominant coefficient this + // isolates the structure-factor phase. + const ModuleBase::Vector3 q_zero_cart(0.0, 0.0, 0.0); + std::vector> drhoc_g0; + rho.compute_drhoc_g(ucell_, iat, 1, q_zero_cart, drhoc_g0); + ASSERT_EQ(static_cast(drhoc_g0.size()), pw_rho_.npw); + // find the G with the largest |G|=0-adjacent coefficient (where the + // radial factor is comparable) and verify the phase relation + int ig_best = -1; + double best_mag = 0.0; + for (int ig = 0; ig < pw_rho_.npw; ++ig) + { + const double m = std::norm(drhoc_g0[ig]); + if (m > best_mag) + { + best_mag = m; + ig_best = ig; + } + } + ASSERT_GE(ig_best, 0); + // the expected extra phase from q is exp(-i q_cart . tau_cart): + // drhoc_q(G) = -i|G+q| tpiba rho~(|G+q|) exp(-i(G+q).tau) + // drhoc_0(G) = -i|G| tpiba rho~(|G|) exp(-i G. tau) + // ratio = [|G+q|/|G| · rho~(|G+q|)/rho~(|G|)] · exp(-i q. tau) + // The bracket is positive real; the phase is exactly -q_cart . tau_cart. + // Use the full TWO_PI-scaled phase (the tau and q_cart are already in + // consistent units: tau in lat0-cartesian, q_cart in 2pi/lat0 units, + // so the dot product directly gives the phase in radians). + const std::complex ratio = drhoc_g[ig_best] / drhoc_g0[ig_best]; + const ModuleBase::Vector3 tau_cart_alat = ucell_.atoms[it].tau[ucell_.iat2ia[iat]]; + const double expected_phase_arg = -ModuleBase::TWO_PI * (q2_cart * tau_cart_alat) / ucell_.lat0; + const double cos_arg = std::cos(expected_phase_arg); + const double sin_arg = std::sin(expected_phase_arg); + const std::complex expected_phase(cos_arg, sin_arg); + // tight tolerance: the exact phase relation is exp(-i q.tau), so both + // omission (arg=0) and reversal (arg=+q.tau) must be rejected; the + // only residual is the positive-real radial-ratio bracket (which has + // zero phase contribution) and floating-point noise + // rule 14: named locals for computed function arguments + const std::complex unit_phase_ratio = ratio / expected_phase; + const double phase_diff = std::arg(unit_phase_ratio); + const double sin_diff = std::sin(phase_diff); + const double cos_diff = std::cos(phase_diff); + const double phase_diff_wrapped = std::atan2(sin_diff, cos_diff); + EXPECT_NEAR(phase_diff_wrapped, 0.0, 0.05) << "phase of drhoc_q/drhoc_0 ratio vs expected exp(-i q.tau)"; +} diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp index dfbc4f70923..aa8ed639268 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.cpp @@ -2,6 +2,8 @@ #include "source_base/constants.h" +#include + void DFPTSerialBase::SetUp() { SetUpCell(); @@ -147,6 +149,52 @@ void DFPTSerialBase::MakeTwoAtomCell() ucell_.iat2ia[1] = 0; } +void DFPTSerialBase::MakeNlccTwoAtomCell() +{ + // start from the symmetric-broken two-atom cell, then upgrade the + // second species (the one at the non-origin tau_nlcc_) to an NC + // pseudo with a nonlinear core correction + MakeTwoAtomCell(); + Atom& at = ucell_.atoms[1]; + at.label = "Bc"; + at.coulomb_potential = false; + pseudo& p = at.ncpp; + p.zv = 2.0; + p.tvanp = false; + p.has_so = false; + p.nbeta = 0; + p.lll = {}; + p.nh = 0; + p.msh = 121; + p.kkbeta = 121; + p.r.resize(121); + p.rab.resize(121); + p.vloc_at.assign(121, 0.0); + const double dx = 0.025; + for (int i = 0; i < 121; ++i) + { + p.r[i] = i * dx; + p.rab[i] = dx; + } + // Gaussian core charge normalized to the integrated nlcc_core_charge_ + // (rho_atc carries the 4 pi r^2 measure of the radial grid) + const double s2 = nlcc_sigma_ * nlcc_sigma_; + double norm = 0.0; + for (int i = 0; i < 121; ++i) + { + const double r = p.r[i]; + norm += r * r * std::exp(-r * r / (2.0 * s2)) * p.rab[i]; + } + const double amp = nlcc_core_charge_ / (norm * ModuleBase::FOUR_PI); + p.rho_atc.resize(121); + for (int i = 0; i < 121; ++i) + { + const double r = p.r[i]; + p.rho_atc[i] = amp * std::exp(-r * r / (2.0 * s2)); + } + p.nlcc = true; +} + long long DFPTSerialBase::FKey(int ix, int iy, int iz) const { return (static_cast(ix + 64) * 128 + (iy + 64)) * 128 + (iz + 64); diff --git a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h index 77285d8c7ce..065ee223fd2 100644 --- a/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h +++ b/source/source_pw/module_dfpt/test_serial/dfpt_serial_fixture.h @@ -63,6 +63,21 @@ class DFPTSerialBase : public testing::Test void MakeCoulombAtom(); void MakeNCAtom(); + /// extend the two-atom cell of MakeTwoAtomCell with a nonlinear core + /// correction on the second species: a Gaussian rho_atc of total charge + /// nlcc_core_charge_ on the existing radial grid (msh points). The + /// atom sits at the non-origin tau_nlcc_ = (0.25, 0.31, 0.17) -- the + /// anti-blind-spot position: every structure-factor phase is active, + /// so a mirrored phase sign (exp(+iG.tau) instead of the QE addcore + /// exp(-iG.tau)) places the response at the wrong site and is caught + /// directly by the drhoc-localization test. + void MakeNlccTwoAtomCell(); + + /// total integrated core charge (in e) of the nlcc Gaussian above + static constexpr double nlcc_core_charge_ = 0.24; + const ModuleBase::Vector3 tau_nlcc_{0.25, 0.31, 0.17}; // lat0 units + double nlcc_sigma_ = 1.5; // bohr, grid-resolvable (spacing ~0.63 bohr) + // reconfigure the cell as a two-atom Z=4/Z=2 crystal breaking all symmetry void MakeTwoAtomCell();