Skip to content

EXX cal_dHs including Hellmann-Feynman term #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
maki49 wants to merge 4 commits into
base: master
Choose a base branch
from
Open

EXX cal_dHs including Hellmann-Feynman term #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
dd67e6b
dH (Pulay only)
maki49 Jun 7, 2026
b058950
add Hellmann-Feynman term of dHs
maki49 Jun 7, 2026
78a0dec
reverse signs of dH to match the real physics
maki49 Jun 8, 2026
3c9df93
optimize performance according to copilot's suggestion
maki49 Jun 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion include/RI/physics/Exx.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,15 +72,22 @@ 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<std::string,3> &save_names_suffix={"","",""}); // "Cs","Vs","Ds"
void cal_force(
const std::array<std::string,5> &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCs","dVs"
void cal_dHs(
const std::array<std::string, 5>& save_names_suffix = { "","","","","" }); // "Cs","Vs","Ds","dCs","dVs"
void cal_stress(
const std::array<std::string,5> &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCRs","dVRs"

std::map<TA, std::map<TAC, Tensor<Tdata>>> Hs;
std::array<std::array< std::map<TA, std::map<TAC, Tensor<Tdata>>> ,2>,Npos> dHs;
std::array<std::array<std::map<TA, std::map<TAC, Tensor<Tdata>>>, 2>, Npos> dHs; // Pulay term only. [direction][diffed atom 0/1][I][JR]
std::array<std::map<TA, std::map<TA, std::map<TAC, Tensor<Tdata>>>>, 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<std::array< std::map<TA, std::map<TAC, Tensor<Tdata>>> ,Npos>,Npos> dHRs;
Tdata energy = 0;
std::array<std::map<TA,Tdata>,Ndim> force;
Expand Down
165 changes: 165 additions & 0 deletions include/RI/physics/Exx.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -385,6 +385,171 @@ void Exx<TA,Tcell,Ndim,Tdata>::cal_force(
}


template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
void Exx<TA, Tcell, Ndim, Tdata>::cal_dHs(
const std::array<std::string, 5>& 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);
}

// ---- 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<TA, std::map<TAC, Tensor<Tdata>>>& Ds = this->lri.data_pool.at(Ds_full).Ds_ab;
std::vector<TA> 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<TA, std::map<TAC, Tensor<Tdata>>> slice) {
Data_Pack<TA, TC, Tdata> 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);
};

// 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<TA, std::map<TA, std::map<TAC, Tensor<Tdata>>>> 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)
{
Comment thread
maki49 marked this conversation as resolved.
// ===== a-side (K = Apin): slice Ds by OUTER key (TA) =====
{
std::map<TA, std::map<TAC, Tensor<Tdata>>> 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): use precomputed reverse index =====
{
std::map<TA, std::map<TAC, Tensor<Tdata>>> slice;
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
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);
Comment on lines +487 to +547
}
} // end for(ipos)
}



template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
void Exx<TA,Tcell,Ndim,Tdata>::cal_stress(
Expand Down