From dd67e6b237498d635b54988f5897bac3015cf88e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Mon, 8 Jun 2026 02:20:11 +0800 Subject: [PATCH 1/4] dH (Pulay only) --- include/RI/physics/Exx.h | 5 +++ include/RI/physics/Exx.hpp | 71 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) diff --git a/include/RI/physics/Exx.h b/include/RI/physics/Exx.h index e7457d1..c69154d 100644 --- a/include/RI/physics/Exx.h +++ b/include/RI/physics/Exx.h @@ -72,10 +72,15 @@ class Exx const Tdata_real &threshold, const std::string &save_name_suffix=""); + void set_flag_save_dHs(const bool flag) { this->flag_save_result.dHs = flag; } + void set_flag_save_dHRs(const bool flag) { this->flag_save_result.dHRs = flag; } + void cal_Hs( const std::array &save_names_suffix={"","",""}); // "Cs","Vs","Ds" void cal_force( const std::array &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCs","dVs" + void cal_dHs( + const std::array& save_names_suffix = { "","","","","" }); // "Cs","Vs","Ds","dCs","dVs" void cal_stress( const std::array &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCRs","dVRs" diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp index 2bb8be3..32ff8db 100644 --- a/include/RI/physics/Exx.hpp +++ b/include/RI/physics/Exx.hpp @@ -385,6 +385,77 @@ void Exx::cal_force( } +template +void Exx::cal_dHs( + const std::array& save_names_suffix) // "Cs","Vs","Ds","dCs","dVs" +{ + assert(this->flag_finish.stru); + assert(this->flag_finish.Cs); + assert(this->flag_finish.Vs); + assert(this->flag_finish.Ds); + assert(this->flag_finish.dCs); + assert(this->flag_finish.dVs); + + for (const Label::ab label : {Label::ab::a1b1, Label::ab::a1b2, Label::ab::a2b1, Label::ab::a2b2}) + this->lri.data_ab_name[label] = "Ds_" + save_names_suffix[2]; + for (std::size_t ipos = 0; ipos < Npos; ++ipos) + { + { + this->dHs[ipos][0].clear(); + + this->lri.data_ab_name[Label::ab::a] = "dCs_" + std::to_string(ipos) + "_" + save_names_suffix[3]; + this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + save_names_suffix[1]; + this->lri.data_ab_name[Label::ab::b] = "Cs_" + save_names_suffix[0]; + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a1b1, + Label::ab_ab::a0b0_a1b2, }, + this->dHs[ipos][0], + -1.0); + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a2b1, + Label::ab_ab::a0b0_a2b2 }, + this->dHs[ipos][0], + 1.0); + + this->lri.data_ab_name[Label::ab::a] = "Cs_" + save_names_suffix[0]; + this->lri.data_ab_name[Label::ab::a0b0] = "dVs_" + std::to_string(ipos) + "_" + save_names_suffix[4]; + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a2b2, + Label::ab_ab::a0b0_a2b1 }, + this->dHs[ipos][0], + 1.0); + } + { + this->dHs[ipos][1].clear(); + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a2b2, + Label::ab_ab::a0b0_a1b2 }, + this->dHs[ipos][1], + -1.0); + + this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + save_names_suffix[1]; + this->lri.data_ab_name[Label::ab::b] = "dCs_" + std::to_string(ipos) + "_" + save_names_suffix[3]; + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a1b1, + Label::ab_ab::a0b0_a2b1 }, + this->dHs[ipos][1], + -1.0); + + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a1b2, + Label::ab_ab::a0b0_a2b2 }, + this->dHs[ipos][1], + 1.0); + } + } // end for(ipos) +} + + template void Exx::cal_stress( From b058950a375af4ab95bf9fc1f0ffa7be21aa73dd Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Mon, 8 Jun 2026 04:58:44 +0800 Subject: [PATCH 2/4] add Hellmann-Feynman term of dHs --- include/RI/physics/Exx.h | 4 +- include/RI/physics/Exx.hpp | 92 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 1 deletion(-) diff --git a/include/RI/physics/Exx.h b/include/RI/physics/Exx.h index c69154d..77e483d 100644 --- a/include/RI/physics/Exx.h +++ b/include/RI/physics/Exx.h @@ -85,7 +85,9 @@ class Exx const std::array &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCRs","dVRs" std::map>> Hs; - std::array>> ,2>,Npos> dHs; + std::array>>, 2>, Npos> dHs; // Pulay term only. [direction][diffed atom 0/1][I][JR] + std::array>>>, Npos> dHs_HF; // Hellmann-Feynman term only. [direction][diffed atom on {K,L}, 0-Natom][I][JR]. Filled by cal_dHs; like dHs it is left rank-PARTIAL -- the consumer must sum-reduce it on (diffed atom, I, JR). + std::array>> ,Npos>,Npos> dHRs; Tdata energy = 0; std::array,Ndim> force; diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp index 32ff8db..ec7de7d 100644 --- a/include/RI/physics/Exx.hpp +++ b/include/RI/physics/Exx.hpp @@ -452,6 +452,98 @@ void Exx::cal_dHs( this->dHs[ipos][1], 1.0); } + + // ---- Hellmann-Feynman term: dHs_HF[ipos][Apin(I0)][I][{J,R}] ---- + // Derivative w.r.t. the CONTRACTED (density-matrix) atoms K (a-side) and L (b-side), + // grouped by that atom and accumulated into dHs_HF[ipos][Apin(I0)]. + // + // Physics (each stored 2-center C and V obeys d/dtau_K = -d/dtau_I): + // a-side (K): the dC part for ALL 4 labels equals MINUS the Pulay dC part (so the signs + // are flipped vs cal_dHs block 0), plus the dV part on the a1 labels only + // (V's a-abf center is the contracted K only there). + // b-side (L): symmetric -- dC in the b-slot (signs flipped vs Pulay block 1), dV on b1. + // + // Pinning: the density matrix's OUTER key is the a-side contracted atom K in every label + // (a1->Aa01, a2->Aa2); its INNER key is the b-side contracted atom L. So we restrict K + // (resp. L) by slicing the shared Ds to a single outer (resp. inner) atom and pointing all + // four D-slots at the slice. We must NOT pin the C slot: its outer key is the EXTERNAL atom + // I/J for the a2/b2 labels, which would mis-attribute those terms to I/J instead of K/L. + // NOTE: like dHs, dHs_HF is left rank-PARTIAL; the consumer must sum-reduce it on (Apin(I0),I,{J,R}). + { + this->dHs_HF[ipos].clear(); + const std::string &sfx0 = save_names_suffix[0]; // Cs + const std::string &sfx1 = save_names_suffix[1]; // Vs + const std::string dC = "dCs_" + std::to_string(ipos) + "_" + save_names_suffix[3]; + const std::string dV = "dVs_" + std::to_string(ipos) + "_" + save_names_suffix[4]; + const std::string Ds_full = "Ds_" + save_names_suffix[2]; + const std::string tmp = "__Exx_HF_Ds_tmp"; + + const std::map>>& Ds = this->lri.data_pool.at(Ds_full).Ds_ab; + std::vector pin_atoms; + pin_atoms.reserve(Ds.size()); + for (const auto& kv : Ds) + pin_atoms.push_back(kv.first); + + auto set_D_slots = [&](const std::string& name) { + for (const Label::ab lab : { Label::ab::a1b1, Label::ab::a1b2, Label::ab::a2b1, Label::ab::a2b2 }) + this->lri.data_ab_name[lab] = name; + }; + auto install_slice = [&](std::map>> slice) { + Data_Pack pack; + pack.Ds_ab = std::move(slice); + pack.index_Ds_ab = RI_Tools::get_index(pack.Ds_ab); + this->lri.data_pool[tmp] = std::move(pack); + set_D_slots(tmp); + }; + + for (const TA& Apin : pin_atoms) + { + // ===== a-side (K = Apin): slice Ds by OUTER key (TA) ===== + { + std::map>> slice; + slice[Apin] = Ds.at(Apin); // shallow; tensors shared + install_slice(std::move(slice)); + } + // dC part, sign-flipped vs Pulay block 0: a=dCs, a0b0=Vs, b=Cs + this->lri.data_ab_name[Label::ab::a] = dC; + this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + sfx1; + this->lri.data_ab_name[Label::ab::b] = "Cs_" + sfx0; + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], 1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], -1.0); + // dV part, a1 labels (V a-abf = K): a=Cs, a0b0=dVs, b=Cs + this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; + this->lri.data_ab_name[Label::ab::a0b0] = dV; + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], 1.0); + + // ===== b-side (L = Apin): slice Ds by INNER key (TAC.first) ===== + { + std::map>> slice; + for (const auto& kv : Ds) + { + std::map> inner; + for (const auto& lc : kv.second) + if (lc.first.first == Apin) + inner.insert(lc); // shallow + if (!inner.empty()) + slice[kv.first] = std::move(inner); + } + install_slice(std::move(slice)); + } + // dC part, sign-flipped vs Pulay block 1, in b-slot: a=Cs, a0b0=Vs, b=dCs + this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; + this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + sfx1; + this->lri.data_ab_name[Label::ab::b] = dC; + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], 1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], -1.0); + // dV part, b1 labels (V b-abf = L): a=Cs, a0b0=dVs, b=Cs + this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; + this->lri.data_ab_name[Label::ab::a0b0] = dV; + this->lri.data_ab_name[Label::ab::b] = "Cs_" + sfx0; + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], -1.0); + } + set_D_slots(Ds_full); // restore the four density-matrix slots + this->lri.data_pool.erase(tmp); + } } // end for(ipos) } From 78a0decea9901bab2a9529c9e7c6274a4d6ca4ad Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Mon, 8 Jun 2026 18:13:24 +0800 Subject: [PATCH 3/4] reverse signs of dH to match the real physics --- include/RI/physics/Exx.hpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp index ec7de7d..fef5ba1 100644 --- a/include/RI/physics/Exx.hpp +++ b/include/RI/physics/Exx.hpp @@ -411,13 +411,13 @@ void Exx::cal_dHs( { Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2, }, this->dHs[ipos][0], - -1.0); + 1.0); this->lri.cal_loop3( { Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2 }, this->dHs[ipos][0], - 1.0); + -1.0); this->lri.data_ab_name[Label::ab::a] = "Cs_" + save_names_suffix[0]; this->lri.data_ab_name[Label::ab::a0b0] = "dVs_" + std::to_string(ipos) + "_" + save_names_suffix[4]; @@ -426,7 +426,7 @@ void Exx::cal_dHs( { Label::ab_ab::a0b0_a2b2, Label::ab_ab::a0b0_a2b1 }, this->dHs[ipos][0], - 1.0); + -1.0); } { this->dHs[ipos][1].clear(); @@ -435,7 +435,7 @@ void Exx::cal_dHs( { Label::ab_ab::a0b0_a2b2, Label::ab_ab::a0b0_a1b2 }, this->dHs[ipos][1], - -1.0); + 1.0); this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + save_names_suffix[1]; this->lri.data_ab_name[Label::ab::b] = "dCs_" + std::to_string(ipos) + "_" + save_names_suffix[3]; @@ -444,13 +444,13 @@ void Exx::cal_dHs( { Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs[ipos][1], - -1.0); + 1.0); this->lri.cal_loop3( { Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b2 }, this->dHs[ipos][1], - 1.0); + -1.0); } // ---- Hellmann-Feynman term: dHs_HF[ipos][Apin(I0)][I][{J,R}] ---- @@ -508,12 +508,12 @@ void Exx::cal_dHs( this->lri.data_ab_name[Label::ab::a] = dC; this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + sfx1; this->lri.data_ab_name[Label::ab::b] = "Cs_" + sfx0; - this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], 1.0); - this->lri.cal_loop3({ Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], -1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], -1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a2b1, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], 1.0); // dV part, a1 labels (V a-abf = K): a=Cs, a0b0=dVs, b=Cs this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; this->lri.data_ab_name[Label::ab::a0b0] = dV; - this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], 1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], -1.0); // ===== b-side (L = Apin): slice Ds by INNER key (TAC.first) ===== { @@ -533,13 +533,13 @@ void Exx::cal_dHs( this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; this->lri.data_ab_name[Label::ab::a0b0] = "Vs_" + sfx1; this->lri.data_ab_name[Label::ab::b] = dC; - this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], 1.0); - this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], -1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], -1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b2, Label::ab_ab::a0b0_a2b2 }, this->dHs_HF[ipos][Apin], 1.0); // dV part, b1 labels (V b-abf = L): a=Cs, a0b0=dVs, b=Cs this->lri.data_ab_name[Label::ab::a] = "Cs_" + sfx0; this->lri.data_ab_name[Label::ab::a0b0] = dV; this->lri.data_ab_name[Label::ab::b] = "Cs_" + sfx0; - this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], -1.0); + this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a2b1 }, this->dHs_HF[ipos][Apin], 1.0); } set_D_slots(Ds_full); // restore the four density-matrix slots this->lri.data_pool.erase(tmp); From 3c9df93ed8dd3c63a9448397312b0ba647a4a4c8 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Thu, 18 Jun 2026 22:20:21 +0800 Subject: [PATCH 4/4] optimize performance according to copilot's suggestion --- include/RI/physics/Exx.hpp | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp index fef5ba1..6139b79 100644 --- a/include/RI/physics/Exx.hpp +++ b/include/RI/physics/Exx.hpp @@ -496,6 +496,14 @@ void Exx::cal_dHs( set_D_slots(tmp); }; + // Precompute reverse index: inner atom L -> { outer atom K -> entries }. + // still D[K][LR] but stored as D[L][KR] where R is always L-K + // Avoids an O(NK) full scan per pinned atom in the b-side slice below. + std::map>>> Ds_by_inner; + for (const auto& kv : Ds) + for (const auto& lc : kv.second) + Ds_by_inner[lc.first.first][kv.first].insert(lc); // shallow; tensors shared + for (const TA& Apin : pin_atoms) { // ===== a-side (K = Apin): slice Ds by OUTER key (TA) ===== @@ -515,18 +523,12 @@ void Exx::cal_dHs( this->lri.data_ab_name[Label::ab::a0b0] = dV; this->lri.cal_loop3({ Label::ab_ab::a0b0_a1b1, Label::ab_ab::a0b0_a1b2 }, this->dHs_HF[ipos][Apin], -1.0); - // ===== b-side (L = Apin): slice Ds by INNER key (TAC.first) ===== + // ===== b-side (L = Apin): use precomputed reverse index ===== { std::map>> slice; - for (const auto& kv : Ds) - { - std::map> inner; - for (const auto& lc : kv.second) - if (lc.first.first == Apin) - inner.insert(lc); // shallow - if (!inner.empty()) - slice[kv.first] = std::move(inner); - } + const auto it = Ds_by_inner.find(Apin); + if (it != Ds_by_inner.end()) + slice = it->second; // shallow; tensors shared install_slice(std::move(slice)); } // dC part, sign-flipped vs Pulay block 1, in b-slot: a=Cs, a0b0=Vs, b=dCs