diff --git a/source/source_esolver/esolver_dfpt_pw.cpp b/source/source_esolver/esolver_dfpt_pw.cpp index 79637856b4..a9c1bfca8c 100644 --- a/source/source_esolver/esolver_dfpt_pw.cpp +++ b/source/source_esolver/esolver_dfpt_pw.cpp @@ -223,6 +223,7 @@ void ESolver_DFPT_PW::before_all_runners(BaseCell& basecell, const Input_para& i dfpt_->set_conv_thr(inp.dfpt_conv_thr); dfpt_->set_max_iter(inp.dfpt_max_iter); dfpt_->set_mix_beta(inp.dfpt_mix_beta); + dfpt_->set_asr(inp.dfpt_asr); dfpt_->set_compute_q0(inp.dfpt_compute_q0); dfpt_->set_loto(inp.dfpt_loto); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 5fba78e9b2..3edfddb418 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -439,6 +439,7 @@ struct Input_para double dfpt_conv_thr = 1.0e-8; ///< convergence threshold of the DFPT first-order density int dfpt_max_iter = 100; ///< max iterations of the DFPT first-order density mixing double dfpt_mix_beta = 0.4; ///< mixing coefficient of the DFPT first-order density + std::string dfpt_asr = "no"; ///< acoustic sum rule on the Gamma dynamical matrix: no | simple | crystal // ============== #Parameters (11.Output) =========================== int out_stru = 1; ///< output stru file each ion step diff --git a/source/source_io/module_parameter/read_inp_dfpt.cpp b/source/source_io/module_parameter/read_inp_dfpt.cpp index f4eceef9d5..3c9052bbbb 100644 --- a/source/source_io/module_parameter/read_inp_dfpt.cpp +++ b/source/source_io/module_parameter/read_inp_dfpt.cpp @@ -168,5 +168,41 @@ void ReadInput::item_dfpt() }; this->add_item(item); } + { + Input_Item item("dfpt_asr"); + item.annotation = "acoustic sum rule on the Gamma dynamical matrix: no, simple or crystal"; + item.category = category; + item.type = "String"; + item.description = "Enforce the translational acoustic sum rule (sum over atoms of the " + "Gamma force-constant rows = 0) on the RAW Gamma dynamical matrix, " + "before the LO-TO nonanalytic term, together with the Z* sum rule " + "sum_a Z*_a = 0 on the Born effective charges. Unconstrained " + "linear-response self-consistency violates both at the percent to " + "tens-of-percent level (QE ph.x shows the same), and the reference " + "workflow post-processes with asr='crystal'. simple: subtract each " + "row sum from its diagonal element. crystal: the optimized " + "minimal-norm symmetric projection (keeps Hermiticity and " + "symmetry-equivalent atoms equivalent). no: leave the matrix as " + "computed."; + item.default_value = "no"; + item.read_value = [](const Input_Item& item, Parameter& para) { + if (item.get_size() == 0) + { + para.input.dfpt_asr = "no"; + } + else + { + para.input.dfpt_asr = strvalue; + } + }; + sync_string(input.dfpt_asr); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.dfpt_asr != "no" && para.input.dfpt_asr != "simple" && para.input.dfpt_asr != "crystal") + { + ModuleBase::WARNING_QUIT("ReadInput", "dfpt_asr must be no, simple or crystal."); + } + }; + this->add_item(item); + } } } // namespace ModuleIO diff --git a/source/source_pw/module_dfpt/dfpt_phon.cpp b/source/source_pw/module_dfpt/dfpt_phon.cpp index 678f5071a5..de9fa67306 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.cpp +++ b/source/source_pw/module_dfpt/dfpt_phon.cpp @@ -163,6 +163,315 @@ void DFPT_Phon::diagonalize(int q_idx, DFPT_PW_Data& data) ModuleBase::timer::end("DFPT_Phon", "diagonalize"); } +void DFPT_Phon::audit_asr(const std::string& label, int q_idx, const DFPT_PW_Data& data) const +{ + ModuleBase::TITLE("DFPT_Phon", "audit_asr"); + ModuleBase::timer::start("DFPT_Phon", "audit_asr"); + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const bool at_gamma = std::abs(q_frac.x) <= gamma_tol && std::abs(q_frac.y) <= gamma_tol + && std::abs(q_frac.z) <= gamma_tol; + const int nat = ucell_->nat; + const int nat3 = 3 * nat; + const ModuleBase::ComplexMatrix dyn = data.get_dynmat(q_idx); + if (!at_gamma || dyn.nr != nat3) + { + ModuleBase::timer::end("DFPT_Phon", "audit_asr"); + return; + } + double max_elem = 0.0; + for (int i = 0; i < nat3; ++i) + { + for (int j = 0; j < nat3; ++j) + { + max_elem = std::max(max_elem, std::abs(dyn(i, j))); + } + } + // per-atom PHYSICAL translational row sums: the stored matrix is + // mass-reduced (D = Phi/sqrt(m_a m_b)), so the Phi row sums read + // s(ia)[da] = sum_jb sqrt(m_jb) D(3ia+da, jb) -- the same quantity + // QE's unreduced dyn rows sum to + std::vector sqm(nat, 0.0); + for (int ia = 0; ia < nat; ++ia) + { + sqm[ia] = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass); + } + std::cout << "DFPT ASR audit (" << label << "):"; + double worst = 0.0; + int worst_atom = -1; + for (int ia = 0; ia < nat; ++ia) + { + // per-DIRECTION residuals R(a,alpha,i) = sum_b sqrt(m_b) D(3a+alpha, + // 3b+i): summing over atoms at a FIXED translation direction i + // before combining norms, so opposite x/y residuals cannot cancel + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int ib = 0; ib < nat; ++ib) + { + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); + } + s2 += std::norm(s); + } + } + const double snorm = std::sqrt(s2); + std::cout << " atom" << ia + 1 << "=" << std::setprecision(3) << std::scientific << snorm; + if (snorm > worst) + { + worst = snorm; + worst_atom = ia; + } + } + std::cout << std::defaultfloat << std::setprecision(6); + if (max_elem > 0.0 && worst_atom >= 0) + { + std::cout << " | worst rel = " << std::setprecision(3) << std::scientific << worst / max_elem + << std::defaultfloat << std::setprecision(6); + } + std::cout << std::endl; + ModuleBase::timer::end("DFPT_Phon", "audit_asr"); +} + +void DFPT_Phon::apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data) +{ + ModuleBase::TITLE("DFPT_Phon", "apply_asr"); + ModuleBase::timer::start("DFPT_Phon", "apply_asr"); + const double gamma_tol = 1.0e-8; ///< empirical parameter: fractional-q Gamma tolerance + const ModuleBase::Vector3 q_frac = data.get_qvec(q_idx); + const bool at_gamma = std::abs(q_frac.x) <= gamma_tol && std::abs(q_frac.y) <= gamma_tol + && std::abs(q_frac.z) <= gamma_tol; + const int nat = ucell_->nat; + const int nat3 = 3 * nat; + ModuleBase::ComplexMatrix dyn = data.get_dynmat(q_idx); + if (!at_gamma || dyn.nr != nat3) + { + std::cout << "DFPT ASR (" << asr_type << ") skipped: not a Gamma matrix" << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); + return; + } + + // The stored matrix is mass-reduced, D = Phi / sqrt(m_a m_b) (amu masses + // from the cell); the PHYSICAL sum rule acts on Phi rows, + // sum_b sqrt(m_b) D(3a+da, 3b+db) = 0, and the acoustic eigenvectors of + // D are v(3a+da) = sqrt(m_a) e_da (uniform physical displacement). QE + // dynmat works on the unreduced Phi, so its unweighted rows are this + // same rule; here the constraint vectors carry the sqrt(m) weights. + std::vector sqm(nat, 0.0); + double mtot = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + sqm[ia] = std::sqrt(ucell_->atoms[ucell_->iat2it[ia]].mass); + mtot += ucell_->atoms[ucell_->iat2it[ia]].mass; + } + // per-DIRECTION physical residual R(a,alpha,i) = sum_b sqrt(m_b) + // D(3a+alpha, 3b+i): summing over atoms at a fixed translation + // direction i so opposite x/y residuals cannot cancel + const auto physical_residual_norm = [&](int ia) { + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int ib = 0; ib < nat; ++ib) + { + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); + } + s2 += std::norm(s); + } + } + return std::sqrt(s2); + }; + + // violation BEFORE the correction (physical residuals) + double worst_before = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + worst_before = std::max(worst_before, physical_residual_norm(ia)); + } + + // the crystal projection D -> (1-P) D (1-P) with the orthonormal + // physical translations w_i(3a+da) = sqrt(m_a) delta(da,i)/sqrt(sum m): + // zeroes every per-direction physical residual exactly, keeps + // Hermiticity and symmetry-equivalent atoms equivalent, and leaves the + // acoustic eigenvectors in the null space so the LO-TO term (which + // decouples from them through sum_a Z*_a = 0) cannot lift the acoustic + // frequencies. Both 'simple' and 'crystal' finish with this projection: + // 'simple' first dumps the residual onto the onsite blocks (the QE + // legacy flavor), then the projection restores Hermiticity and removes + // whatever onsite-only correction could not absorb (cross-Cartesian + // antisymmetric residuals have no exact Hermitian onsite-only fix) + const auto apply_crystal_projection = [&]() { + std::vector> s(3 * nat3, std::complex(0.0, 0.0)); + std::vector> t(3 * nat3, std::complex(0.0, 0.0)); + ModuleBase::ComplexMatrix u(3, 3, true); + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int nb = 0; nb < nat; ++nb) + { + for (int db = 0; db < 3; ++db) + { + const std::complex d = dyn(3 * na + da, 3 * nb + db); + s[da * nat3 + 3 * nb + db] += sqm[na] * d; + t[(3 * na + da) * 3 + db] += sqm[nb] * d; + u(da, db) += sqm[na] * sqm[nb] * d; + } + } + } + } + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int nb = 0; nb < nat; ++nb) + { + for (int db = 0; db < 3; ++db) + { + std::complex d = dyn(3 * na + da, 3 * nb + db); + d -= sqm[na] / mtot * s[da * nat3 + 3 * nb + db]; + d -= sqm[nb] / mtot * t[(3 * na + da) * 3 + db]; + d += sqm[na] * sqm[nb] / (mtot * mtot) * u(da, db); + dyn(3 * na + da, 3 * nb + db) = d; + } + } + } + } + }; + + if (asr_type == "simple") + { + // legacy QE 'simple' translated to the reduced matrix: dump the + // per-direction residual onto the onsite 3x3 block, + // D(a,a)[da][db] = -sum_{b!=a} sqrt(m_b/m_a) D(a,da; b,db); + // the onsite dump may break Hermiticity (cross-Cartesian parts), + // so symmetrize, then the crystal projection below removes the + // residual the symmetrization reintroduced (at most half the + // anti-Hermitian part) and guarantees the final matrix is + // Hermitian and exactly satisfies all three translation constraints + for (int na = 0; na < nat; ++na) + { + for (int da = 0; da < 3; ++da) + { + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int nb = 0; nb < nat; ++nb) + { + if (nb != na) + { + s += sqm[nb] * dyn(3 * na + da, 3 * nb + db); + } + } + dyn(3 * na + da, 3 * na + db) = -s / sqm[na]; + } + } + } + // Hermitian re-symmetrization of the onsite dump (the projection + // is Hermitian-preserving only on Hermitian input); the j > i loop + // alone leaves the diagonal untouched, and the onsite dump may + // have written complex values there -- take the real part so the + // matrix is genuinely Hermitian before the projection + for (int i = 0; i < nat3; ++i) + { + for (int j = i + 1; j < nat3; ++j) + { + const std::complex avg = 0.5 * (dyn(i, j) + std::conj(dyn(j, i))); + dyn(i, j) = avg; + dyn(j, i) = std::conj(avg); + } + dyn(i, i) = std::complex(dyn(i, i).real(), 0.0); + } + apply_crystal_projection(); + } + else if (asr_type == "crystal") + { + apply_crystal_projection(); + } + else + { + std::cout << "DFPT ASR: unknown type '" << asr_type << "', nothing applied" << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); + return; + } + + // Hermitian re-symmetrization (rows are corrected independently) + for (int i = 0; i < nat3; ++i) + { + for (int j = i + 1; j < nat3; ++j) + { + const std::complex avg = 0.5 * (dyn(i, j) + std::conj(dyn(j, i))); + dyn(i, j) = avg; + dyn(j, i) = std::conj(avg); + } + } + + // Z* sum rule (both enforcing modes, QE set_asr convention): shift by + // the per-component mean over atoms, sum_a Z*_a -> 0 + double zworst_before = 0.0; + ModuleBase::ComplexMatrix zsum(3, 3, true); + int n_born = 0; + for (int ia = 0; ia < nat; ++ia) + { + const ModuleBase::matrix za = data.get_born(ia); + if (za.nr != 3 || za.nc != 3) + { + continue; + } + ++n_born; + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + zsum(a, b) += std::complex(za(a, b), 0.0); + } + } + } + if (n_born == nat) + { + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + zworst_before = std::max(zworst_before, std::abs(zsum(a, b))); + } + } + for (int ia = 0; ia < nat; ++ia) + { + ModuleBase::matrix za = data.get_born(ia); + for (int a = 0; a < 3; ++a) + { + for (int b = 0; b < 3; ++b) + { + za(a, b) -= zsum(a, b).real() / nat; + } + } + data.set_born(ia, za); + } + } + + data.set_dynmat(q_idx, dyn); + + // violation AFTER the correction (physical residuals) + double worst_after = 0.0; + for (int ia = 0; ia < nat; ++ia) + { + worst_after = std::max(worst_after, physical_residual_norm(ia)); + } + std::cout << "DFPT ASR (" << asr_type << ") applied: worst physical residual " << std::setprecision(3) + << std::scientific << worst_before << " -> " << worst_after; + if (n_born == nat) + { + std::cout << " | max |sum_a Z*_a| " << zworst_before << " -> 0"; + } + std::cout << std::defaultfloat << std::setprecision(6) << std::endl; + ModuleBase::timer::end("DFPT_Phon", "apply_asr"); +} + void DFPT_Phon::add_loto(const ModuleBase::Vector3& qhat, DFPT_PW_Data& data) { ModuleBase::TITLE("DFPT_Phon", "add_loto"); diff --git a/source/source_pw/module_dfpt/dfpt_phon.h b/source/source_pw/module_dfpt/dfpt_phon.h index 9e8c18cdce..9bcd8dbaec 100644 --- a/source/source_pw/module_dfpt/dfpt_phon.h +++ b/source/source_pw/module_dfpt/dfpt_phon.h @@ -87,6 +87,40 @@ class DFPT_Phon /// from Gamma, where the rule does not apply). bool check_sum_rule(int q_idx, DFPT_PW_Data& data) const; + /// translational sum-rule audit of the Gamma dynamical matrix + /// (per-atom 3-vector row sums, relative to the matrix scale); + /// diagnostics for the acoustic-branch investigation: the raw matrix + /// and the LO-TO-corrected one are audited separately because the + /// nonanalytic term carries its own Z* sum rule + void audit_asr(const std::string& label, int q_idx, const DFPT_PW_Data& data) const; + + /// enforce the translational acoustic sum rule on the Gamma dynamical + /// matrix (and, in every enforcing mode, the Z* sum rule + /// sum_a Z*_a = 0 on the stored Born charges, like QE dynmat's + /// set_asr). Both modes operate on the mass-reduced matrix D = + /// Phi/sqrt(m_a m_b); the physical rule reads sum_b sqrt(m_b) + /// D(3a+alpha, 3b+beta) = 0 for each translation direction beta + /// independently. + /// + /// Both modes finish with the exact Hermitian projection D -> + /// (1-P) D (1-P), P = sum_i w_i w_i^T, with the orthonormal physical + /// translations w_i(3a+alpha) = sqrt(m_a) delta(alpha,i) / + /// sqrt(sum m): the final matrix satisfies all three mass-weighted + /// translation constraints exactly and is Hermitian. + /// + /// "simple": the QE legacy flavor -- the residual is first dumped + /// onto the onsite 3x3 blocks, D(a,a)[alpha][beta] = + /// -sum_{b!=a} sqrt(m_b/m_a) D(a,alpha; b,beta), then the projection + /// absorbs whatever onsite-only correction could not (a purely + /// onsite Hermitian fix does not exist when the residual has a + /// cross-Cartesian antisymmetric part). + /// + /// "crystal": the optimized minimal-norm fix -- the projection alone, + /// which also keeps symmetry-equivalent atoms equivalent. + /// + /// Away from Gamma the call is a no-op. + void apply_asr(const std::string& asr_type, int q_idx, DFPT_PW_Data& data); + /// Ewald ion-ion force constants C^ewald_ab(q) (G-space + real-space + /// Gaussian self term), mass-reduced by 1/sqrt(M_a M_b). A stateless /// building block exposed publicly so the serial analytic tests can diff --git a/source/source_pw/module_dfpt/dfpt_pw.cpp b/source/source_pw/module_dfpt/dfpt_pw.cpp index 2d325c9005..de3b7a955a 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw.cpp @@ -138,6 +138,14 @@ void DFPT_PW::set_mix_beta(double beta) } } +void DFPT_PW::set_asr(const std::string& asr_type) +{ + if (asr_type == "no" || asr_type == "simple" || asr_type == "crystal") + { + pimpl_->asr_type_ = asr_type; + } +} + void DFPT_PW::set_compute_q0(bool flag) { pimpl_->data_.set_compute_q0(flag); diff --git a/source/source_pw/module_dfpt/dfpt_pw.h b/source/source_pw/module_dfpt/dfpt_pw.h index 83560769a1..ea611743b9 100644 --- a/source/source_pw/module_dfpt/dfpt_pw.h +++ b/source/source_pw/module_dfpt/dfpt_pw.h @@ -128,6 +128,10 @@ class DFPT_PW void set_mix_beta(double beta); + /// acoustic sum rule on the Gamma dynamical matrix ("no" | "simple" | + /// "crystal"); applied to the raw matrix before the LO-TO term + void set_asr(const std::string& asr_type); + /// q = 0 response switches (epsilon_inf / Born charges / LO-TO) void set_compute_q0(bool flag); diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.cpp b/source/source_pw/module_dfpt/dfpt_pw_data.cpp index 754f4e063d..529517c621 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_data.cpp @@ -646,4 +646,34 @@ void DFPT_PW_Data::deallocate_memory() ModuleBase::timer::end("DFPT_PW_Data", "deallocate_memory"); } +void DFPT_PW_Data::set_drhoc_g(int q_idx, const std::vector>& rho_c) +{ + ModuleBase::TITLE("DFPT_PW_Data", "set_drhoc_g"); + ModuleBase::timer::start("DFPT_PW_Data", "set_drhoc_g"); + if (q_idx < 0) + { + ModuleBase::timer::end("DFPT_PW_Data", "set_drhoc_g"); + return; + } + if (q_idx >= static_cast(drhoc_g_.size())) + { + drhoc_g_.resize(q_idx + 1); + } + drhoc_g_[q_idx] = rho_c; + ModuleBase::timer::end("DFPT_PW_Data", "set_drhoc_g"); +} + +std::vector> DFPT_PW_Data::get_drhoc_g(int q_idx) const +{ + ModuleBase::TITLE("DFPT_PW_Data", "get_drhoc_g"); + ModuleBase::timer::start("DFPT_PW_Data", "get_drhoc_g"); + if (q_idx >= 0 && q_idx < static_cast(drhoc_g_.size())) + { + ModuleBase::timer::end("DFPT_PW_Data", "get_drhoc_g"); + return drhoc_g_[q_idx]; + } + ModuleBase::timer::end("DFPT_PW_Data", "get_drhoc_g"); + return {}; +} + } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_data.h b/source/source_pw/module_dfpt/dfpt_pw_data.h index d1ea8dd6f1..239c77fa4a 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_data.h +++ b/source/source_pw/module_dfpt/dfpt_pw_data.h @@ -62,6 +62,11 @@ class DFPT_PW_Data void set_drho_r(int q_idx, int spin, const std::vector& rho); std::vector get_drho_r(int q_idx, int spin) const; void set_drho_g(int q_idx, int spin, const std::vector>& rho); + + /// q-shifted first-order core-density amplitude (G space, npw); + /// zero-size means no nlcc atoms / not yet computed + void set_drhoc_g(int q_idx, const std::vector>& rho_c); + std::vector> get_drhoc_g(int q_idx) const; std::vector> get_drho_g(int q_idx, int spin) const; void set_dv_r(int q_idx, int spin, const std::vector& v); @@ -267,6 +272,8 @@ class DFPT_PW_Data std::vector>> drho_r_; std::vector>>> drho_g_; + ///< q-shifted first-order core density (G space), [q] + std::vector>> drhoc_g_; std::vector>> dv_r_; diff --git a/source/source_pw/module_dfpt/dfpt_pw_impl.h b/source/source_pw/module_dfpt/dfpt_pw_impl.h index 7448034d70..ed32d6ad79 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_impl.h +++ b/source/source_pw/module_dfpt/dfpt_pw_impl.h @@ -177,6 +177,7 @@ class DFPT_PW::Impl double conv_thr_ = 1e-8; int max_iter_ = 100; double mix_beta_ = 0.4; + std::string asr_type_ = "no"; ///< acoustic sum rule mode on the Gamma matrix }; } // namespace ModuleDFPT diff --git a/source/source_pw/module_dfpt/dfpt_pw_run.cpp b/source/source_pw/module_dfpt/dfpt_pw_run.cpp index 16aec15df2..9db6f17e94 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_run.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_run.cpp @@ -108,7 +108,25 @@ void DFPT_PW::Impl::run_q0_post(int q_idx) void DFPT_PW::Impl::run_assemble(int q_idx) { phon_.assemble(q_idx, data_); - phon_.diagonalize(q_idx, data_); + // translational sum-rule audit of the raw Gamma matrix (acoustic-branch + // diagnostics; the loto-corrected matrix is audited separately below). + // The raw audit does not consume frequencies, so the ordinary + // eigensolve runs once after any applicable ASR correction + if (asr_type_ != "no") + { + phon_.audit_asr("raw", q_idx, data_); + // enforce the sum rule on the raw matrix (with the Z* sum rule), + // diagonalize once and audit; the LO-TO term below then consumes + // the corrected Born charges (QE dynmat set_asr order) + phon_.apply_asr(asr_type_, q_idx, data_); + phon_.diagonalize(q_idx, data_); + phon_.audit_asr("asr", q_idx, data_); + } + else + { + phon_.diagonalize(q_idx, data_); + phon_.audit_asr("raw", q_idx, data_); + } if (q_idx != 0 || !data_.get_loto()) { return; @@ -118,6 +136,7 @@ void DFPT_PW::Impl::run_assemble(int q_idx) // set_loto_dir overrides, e.g. per irrep direction in stage A) phon_.add_loto(data_.get_loto_dir(), data_); phon_.diagonalize_loto(data_); + phon_.audit_asr("loto", q_idx, data_); } void DFPT_PW::run() diff --git a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp index a5f3590509..df06df88dc 100644 --- a/source/source_pw/module_dfpt/dfpt_pw_solve.cpp +++ b/source/source_pw/module_dfpt/dfpt_pw_solve.cpp @@ -56,6 +56,19 @@ void DFPT_PW::Impl::assemble_v_sc(const ModuleBase::Vector3& q_cart, xc_->set_context(q_cart); std::vector> a_r(nrxx); 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_); + if (static_cast(drhoc_g.size()) == pw_rho_->npw) + { + std::vector> drhoc_r(nrxx); + pw_rho_->recip2real(drhoc_g.data(), drhoc_r.data()); + for (int ir = 0; ir < nrxx; ++ir) + { + a_r[ir] += drhoc_r[ir]; + } + } std::vector> b_r; xc_->apply(a_r, b_r); if (static_cast(b_r.size()) != nrxx) @@ -203,6 +216,9 @@ double DFPT_PW::Impl::solve_displacement(int q_idx, int iat, int idir) const bool dbg = (getenv("DFPT_DEBUG") != nullptr); pert_.build_dv(q_idx, iat, idir, data_); + // nlcc: q-shifted core-density derivative for the displaced atom + // (fixed per displacement, consumed by assemble_v_sc's XC channel) + rho_.compute_drhoc(*ucell_, iat, idir, q_idx, data_); rho_.reset_mixing(q_idx); data_.set_drho_g(q_idx, 0, std::vector>(npw, std::complex(0.0, 0.0))); diff --git a/source/source_pw/module_dfpt/dfpt_rho.cpp b/source/source_pw/module_dfpt/dfpt_rho.cpp index 2a271c9716..239c7804ab 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.cpp +++ b/source/source_pw/module_dfpt/dfpt_rho.cpp @@ -6,6 +6,9 @@ #include "source_base/timer.h" #include "source_base/tool_title.h" #include "source_base/module_mixing/plain_mixing.h" +#include "source_base/math_integral.h" +#include "source_base/math_sphbes.h" +#include "source_estate/module_charge/chg_tools.h" #include #include @@ -242,6 +245,103 @@ void DFPT_Rho::compute_drho(const psi::Psi>& psi, ModuleBase::timer::end("DFPT_Rho", "compute_drho"); } +void DFPT_Rho::compute_drhoc(const UnitCell& ucell, + int iat, + int dir, + int q_idx, + DFPT_PW_Data& data) +{ + ModuleBase::TITLE("DFPT_Rho", "compute_drhoc"); + ModuleBase::timer::start("DFPT_Rho", "compute_drhoc"); + const int npw = pw_rho_ ? pw_rho_->npw : 0; + 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)); + std::vector aux(mesh); + + for (int ig = 0; ig < npw; ++ig) + { + const ModuleBase::Vector3 w = pw_rho_->gcar[ig] + q_cart; + const double w2 = w * w; + // skip the |G+q| = 0 shell (dropped by compute_drho, zero density there) + if (w2 < 1.0e-12) + { + continue; + } + const double g = std::sqrt(w2) * ucell.tpiba; // physical 1/bohr + // radial FT: rho_core~(|G+q|) = 4pi/Omega * Int r^2 rho_atc(r) j0(gr) dr + ModuleBase::Sphbes::Spherical_Bessel(mesh, p.r.data(), g, 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; + + // structure factor phase: exp(i (G+q) . tau_iat) + const double phase = ModuleBase::TWO_PI * (w * tau_cart) / ucell.lat0; + 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, const ModuleBase::matrix& wg, int q_idx, diff --git a/source/source_pw/module_dfpt/dfpt_rho.h b/source/source_pw/module_dfpt/dfpt_rho.h index 2d4876f3ac..83ae0cda25 100644 --- a/source/source_pw/module_dfpt/dfpt_rho.h +++ b/source/source_pw/module_dfpt/dfpt_rho.h @@ -5,6 +5,7 @@ #include "source_base/matrix3.h" #include "source_basis/module_pw/pw_basis.h" #include "source_basis/module_pw/pw_basis_k.h" +#include "source_cell/unitcell.h" #include "source_psi/psi.h" #include @@ -114,6 +115,22 @@ class DFPT_Rho int q_idx, DFPT_PW_Data& data); + /// q-shifted first-order core-density amplitude for an nlcc + /// pseudopotential atom displaced along (iat, dir), in G space: + /// d rho_core/d tau_dir (G+q) = -i tpiba (G+q)_dir + /// * rho_core~(|G+q|) * exp(i (G+q) tau_iat) + /// where rho_core~ is the radial spherical-Bessel FT of the core + /// charge (the same non_linear_core_correction used by the ground + /// state). Non-nlcc atoms / types give zero. The result enters the + /// XC channel of assemble_v_sc only (the core Coulomb is already in + /// the ionic pseudopotential, so no Hartree contribution). + /// ucell is passed explicitly (DFPT_Rho does not own it). + void compute_drhoc(const UnitCell& ucell, + int iat, + int dir, + int q_idx, + DFPT_PW_Data& data); + /// 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/dfpt_stru_fixture.cpp b/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp index 4df392da1c..4192fe161b 100644 --- a/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp +++ b/source/source_pw/module_dfpt/test/dfpt_stru_fixture.cpp @@ -64,6 +64,21 @@ void DFPTStruTestFixture::construct_ucell(stru_& stru) } ucell.nat += ucell.atoms[i].na; } + // atom-index tables: real cells always carry them (iat2it/iat2ia map the + // running atom index to type and in-type index); populate them so any + // consumer (e.g. DFPT_Phon::apply_asr mass lookup) works on this fixture + ucell.iat2it.resize(ucell.nat); + ucell.iat2ia.resize(ucell.nat); + int iat = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + ucell.iat2it[iat] = it; + ucell.iat2ia[iat] = ia; + ++iat; + } + } } void DFPTStruTestFixture::ClearUcell() 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 2511020066..085de1a93d 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 @@ -691,3 +691,276 @@ TEST_F(DFPTPhonSerialTest, FormatReportsRegression) " mode 2 : 520.123457 cm^-1\n"; EXPECT_EQ(lrep, expect_l); } + +namespace +{ +/// per-atom PHYSICAL residual norm of a 2-atom (nat3=6) mass-reduced +/// dynamical matrix: R(a,alpha,beta) = sum_b sqrt(m_b) D(3a+alpha, 3b+beta) +/// at a FIXED direction beta, then combined -- opposite x/y residuals +/// cannot cancel +double row_sum_norm(const UnitCell& ucell, const ModuleBase::ComplexMatrix& dyn, int ia) +{ + double sqm[2] = {std::sqrt(ucell.atoms[ucell.iat2it[0]].mass), + std::sqrt(ucell.atoms[ucell.iat2it[1]].mass)}; + double s2 = 0.0; + for (int da = 0; da < 3; ++da) + { + for (int db = 0; db < 3; ++db) + { + std::complex s(0.0, 0.0); + for (int ib = 0; ib < 2; ++ib) + { + s += sqm[ib] * dyn(3 * ia + da, 3 * ib + db); + } + s2 += std::norm(s); + } + } + return std::sqrt(s2); +} + +/// Hermitian 2-atom matrix with an asymmetric translational violation: +/// self blocks A1 != A2, off block B, plus an extra bump on atom 0 +ModuleBase::ComplexMatrix make_violating_dyn() +{ + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double self1 = 0.05; + const double self2 = 0.04; + const double off = 0.02; + const double bump = 0.012; // asymmetry: only on atom 0, x row + for (int d = 0; d < 3; ++d) + { + dyn(d, d) = std::complex(self1 + (d == 0 ? bump : 0.0), 0.0); + dyn(3 + d, 3 + d) = std::complex(self2, 0.0); + dyn(d, 3 + d) = std::complex(off, 0.01 * d); + dyn(3 + d, d) = std::conj(dyn(d, 3 + d)); + } + return dyn; +} +} // namespace + +TEST_F(DFPTPhonSerialTest, AsrSimpleZeroesRowSums) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + const ModuleBase::ComplexMatrix dyn0 = make_violating_dyn(); + data_.set_dynmat(0, dyn0); + phon_.apply_asr("simple", 0, data_); + const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); + // legacy fix leaves only the tiny imaginary Hermitization residual + const double scale = 0.06; + EXPECT_LT(row_sum_norm(ucell_, dyn, 0) / scale, 1.0e-12); + EXPECT_LT(row_sum_norm(ucell_, dyn, 1) / scale, 1.0e-12); + // Hermiticity preserved by the re-symmetrization pass + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(dyn(i, j), std::conj(dyn(j, i))); + } + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalZeroesRowSumsExactlyAndProjectsAcousticToZero) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix dyn = data_.get_dynmat(0); + EXPECT_LT(row_sum_norm(ucell_, dyn, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, dyn, 1), 1.0e-14); + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(dyn(i, j), std::conj(dyn(j, i))); + } + } + // the three translation eigenvalues are projected out exactly: the + // uniform-displacement subspace maps to zero under (1-P) D (1-P). + // The matrix residual is ~1e-17; in frequency units the square root + // rescaling lifts that to ~1e-5 cm^-1, still an exact zero physically + phon_.diagonalize(0, data_); + const std::vector freq = data_.get_phon_freq(0); + std::vector byabs = freq; + std::sort(byabs.begin(), byabs.end(), [](double a, double b) { return std::abs(a) < std::abs(b); }); + for (int m = 0; m < 3; ++m) + { + EXPECT_NEAR(byabs[m], 0.0, 1.0e-4); + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalKeepsEquivalentAtomsEquivalent) +{ + MakeTwoAtomCell(); + // symmetry-equivalent atoms carry equal masses; make the two fixture + // species mass-equal so the equivalence assertion is meaningful + ucell_.atoms[1].mass = ucell_.atoms[0].mass; + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + // symmetric input: identical self blocks, common off block, plus a + // violation spread equally (row sums nonzero but atom-symmetric) + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double self = 0.05; + const double off = 0.02; + for (int d = 0; d < 3; ++d) + { + dyn(d, d) = std::complex(self, 0.0); + dyn(3 + d, 3 + d) = std::complex(self, 0.0); + dyn(d, 3 + d) = std::complex(off, 0.0); + dyn(3 + d, d) = std::complex(off, 0.0); + } + data_.set_dynmat(0, dyn); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + // crystal redistributes the mean uniformly: both self blocks stay equal + // and both row sums vanish (simple would dump them on the diagonals, + // also staying equal here -- the discriminating case is the asymmetric + // one above, where only crystal shares the correction) + for (int d = 0; d < 3; ++d) + { + EXPECT_EQ(out(d, d), out(3 + d, 3 + d)); + } +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalCorrectsBornChargeSumRule) +{ + MakeTwoAtomCell(); + // rewire the q list to Gamma (the fixture default q is not 0) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + ModuleBase::matrix z1(3, 3, true); + ModuleBase::matrix z2(3, 3, true); + z1(0, 0) = 1.3; + z1(1, 1) = 1.1; + z1(2, 2) = 1.2; + z2(0, 0) = -1.1; + z2(1, 1) = -1.0; + z2(2, 2) = -1.05; // sums: 0.2 / 0.1 / 0.15 != 0 + data_.set_born(0, z1); + data_.set_born(1, z2); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::matrix o1 = data_.get_born(0); + const ModuleBase::matrix o2 = data_.get_born(1); + for (int a = 0; a < 3; ++a) + { + EXPECT_NEAR(o1(a, a) + o2(a, a), 0.0, 1.0e-14); + } + // shifted by the per-component mean: 0.1 / 0.05 / 0.075 + EXPECT_NEAR(o1(0, 0), 1.3 - 0.1, 1.0e-14); + EXPECT_NEAR(o2(1, 1), -1.1 + 0.05, 1.0e-14); +} + +TEST_F(DFPTPhonSerialTest, AsrCrystalPlusLotoKeepsAcousticExactlyZero) +{ + // regression of the mass-weighting bug: with UNequal masses the + // translation constraints must live on the physical (sqrt(m)-weighted) + // vectors, otherwise the LO-TO nonanalytic term -- which decouples from + // the physical acoustic subspace only through sum_a Z*_a = 0 -- mixes + // into the projected null space and lifts one acoustic frequency (the + // bn_cubic Gamma run showed +36 cm^-1 before the fix) + MakeTwoAtomCell(); // masses 12 / 4: unequal, the trigger + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + data_.set_dynmat(0, make_violating_dyn()); + ModuleBase::matrix eps(3, 3, true); + for (int d = 0; d < 3; ++d) + { + eps(d, d) = 3.0; + } + data_.set_dielectric(eps); + ModuleBase::matrix z1(3, 3, true); + ModuleBase::matrix z2(3, 3, true); + z1(0, 0) = z1(1, 1) = z1(2, 2) = 1.3; + z2(0, 0) = z2(1, 1) = z2(2, 2) = -1.1; // sums 0.2 / 0.2 / 0.2 != 0 + data_.set_born(0, z1); + data_.set_born(1, z2); + phon_.apply_asr("crystal", 0, data_); + // LO-TO along (1,1,1)/sqrt(3): the acoustic triplet must stay at zero + phon_.add_loto(ModuleBase::Vector3(1.0, 1.0, 1.0), data_); + phon_.diagonalize_loto(data_); + const std::vector freq = data_.get_phon_freq_loto(); + ASSERT_EQ(freq.size(), 6u); + std::vector byabs = freq; + std::sort(byabs.begin(), byabs.end(), [](double a, double b) { return std::abs(a) < std::abs(b); }); + for (int m = 0; m < 3; ++m) + { + EXPECT_NEAR(byabs[m], 0.0, 1.0e-4); + } +} + +TEST_F(DFPTPhonSerialTest, AsrDirectionalCancellationDetected) +{ + // regression of the audit-metric fix: the block [[1,-1,0],[-1,1,0]] + // has zero row sums (sum over ALL columns) but a NONZERO per-direction + // residual D*t_x; the audit and apply_asr report must not miss it + MakeTwoAtomCell(); + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + ModuleBase::ComplexMatrix dyn(6, 6, true); + // atom0 self-block with the cancellation pattern + dyn(0, 0) = std::complex(1.0, 0.0); + dyn(0, 1) = std::complex(-1.0, 0.0); + dyn(1, 0) = std::complex(-1.0, 0.0); + dyn(1, 1) = std::complex(1.0, 0.0); + // atom1 self-block mirrors it (keeps total charge-conservation flavor) + dyn(3, 3) = std::complex(1.0, 0.0); + dyn(3, 4) = std::complex(-1.0, 0.0); + dyn(4, 3) = std::complex(-1.0, 0.0); + dyn(4, 4) = std::complex(1.0, 0.0); + data_.set_dynmat(0, dyn); + + // per-direction residual BEFORE any correction is nonzero + const double residual_before = row_sum_norm(ucell_, dyn, 0); + EXPECT_GT(residual_before, 0.1); + + // crystal ASR zeroes it + ModuleBase::ComplexMatrix dyn_copy = dyn; + data_.set_dynmat(0, dyn_copy); + phon_.apply_asr("crystal", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + EXPECT_LT(row_sum_norm(ucell_, out, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, out, 1), 1.0e-14); +} + +TEST_F(DFPTPhonSerialTest, AsrSimpleCrossCartesianRegression) +{ + // regression of the simple-mode Hermitian fix: the cross-Cartesian + // coupling D(0,4) = D(4,0) = delta (atom0-x <-> atom1-y) has no exact + // Hermitian onsite-only fix; simple must still produce a final matrix + // that is Hermitian, satisfies all three translation constraints + // exactly, and has a real diagonal + MakeTwoAtomCell(); // masses 12/4 (unequal) + SetupPhon(k_d_, ModuleBase::Vector3(0.0, 0.0, 0.0)); + ModuleBase::ComplexMatrix dyn(6, 6, true); + const double delta = 0.07; + dyn(0, 4) = std::complex(delta, 0.0); + dyn(4, 0) = std::complex(delta, 0.0); + // also a diagonal self term so the matrix is not purely offsite + dyn(0, 0) = std::complex(0.05, 0.0); + dyn(4, 4) = std::complex(0.04, 0.0); + dyn(2, 2) = std::complex(0.03, 0.0); + dyn(5, 5) = std::complex(0.06, 0.0); + data_.set_dynmat(0, dyn); + phon_.apply_asr("simple", 0, data_); + const ModuleBase::ComplexMatrix out = data_.get_dynmat(0); + + // translation constraints hold exactly (per-direction residuals) + EXPECT_LT(row_sum_norm(ucell_, out, 0), 1.0e-14); + EXPECT_LT(row_sum_norm(ucell_, out, 1), 1.0e-14); + + // Hermiticity + for (int i = 0; i < 6; ++i) + { + for (int j = i + 1; j < 6; ++j) + { + EXPECT_EQ(out(i, j), std::conj(out(j, i))) << "i=" << i << " j=" << j; + } + } + + // diagonal is real + for (int i = 0; i < 6; ++i) + { + EXPECT_EQ(out(i, i).imag(), 0.0) << "i=" << i; + } +}