feat(dfpt): acoustic sum rule enforcement (dfpt_asr = no | simple | crystal)

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);
 }

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

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

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<double> 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<double> 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<double> 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<double> 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<double> 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<double> 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<std::complex<double>> s(3 * nat3, std::complex<double>(0.0, 0.0));
+        std::vector<std::complex<double>> t(3 * nat3, std::complex<double>(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<double> 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<double> 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<double> 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<double> 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<double>(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<double> 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<double>(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<double>& qhat, DFPT_PW_Data& data)
 {
     ModuleBase::TITLE("DFPT_Phon", "add_loto");