Skip to content

Commit 91dc587

Browse files
maki49maki49
committed
restrict to irreducible sector
loop1-filter and if cal_energy refactor: class Symmetry_Filter
1 parent 3883b69 commit 91dc587

File tree

5 files changed

+113
-14
lines changed

5 files changed

+113
-14
lines changed

‎include/RI/physics/Exx.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,9 @@ class Exx
6868
const std::string &save_name_suffix="");
6969

7070
void cal_Hs(
71-
const std::array<std::string,3> &save_names_suffix={"","",""}); // "Cs","Vs","Ds"
71+
const std::array<std::string, 3>& save_names_suffix = { "","","" }, // "Cs","Vs","Ds"
72+
const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector = {},
73+
const bool if_cal_energy = true);
7274
void cal_force(
7375
const std::array<std::string,5> &save_names_suffix={"","","","",""}); // "Cs","Vs","Ds","dCs","dVs"
7476
void cal_stress(

‎include/RI/physics/Exx.hpp‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -177,7 +177,9 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_dVRs(
177177

178178
template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
179179
void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
180-
const std::array<std::string,3> &save_names_suffix) // "Cs","Vs","Ds"
180+
const std::array<std::string, 3>& save_names_suffix, // "Cs","Vs","Ds"
181+
const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector,
182+
const bool if_cal_energy)
181183
{
182184
using namespace Map_Operator;
183185

@@ -208,9 +210,11 @@ void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
208210
Label::ab_ab::a0b0_a1b2,
209211
Label::ab_ab::a0b0_a2b1,
210212
Label::ab_ab::a0b0_a2b2},
211-
this->Hs);
213+
this->Hs,
214+
1.0/*fac_add_Ds*/,
215+
irreducible_sector);
212216

213-
//if()
217+
if (if_cal_energy)
214218
this->energy = this->post_2D.cal_energy(
215219
this->post_2D.saves["Ds_"+save_names_suffix[2]],
216220
this->post_2D.set_tensors_map2(this->Hs) );

‎include/RI/ri/LRI-cal_loop3.hpp‎

Lines changed: 33 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
#include "LRI_Cal_Aux.h"
1010
#include "../global/Array_Operator.h"
1111
#include "../global/Tensor_Multiply.h"
12-
12+
#include "../symmetry/Symmetry_Filter.h"
1313
#include <omp.h>
1414
#ifdef __MKL_RI
1515
#include <mkl_service.h>
@@ -22,8 +22,10 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
2222
void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
2323
const std::vector<Label::ab_ab> &labels,
2424
std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
25-
const double fac_add_Ds)
25+
const double fac_add_Ds,
26+
const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector)
2627
{
28+
std::cout << "using loop3 LibRI" << std::endl;
2729
using namespace Array_Operator;
2830

2931
const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
@@ -70,6 +72,8 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
7072
? LRI_Cal_Aux::cal_Ds_transpose( data_wrapper(Label::ab::b).Ds_ab )
7173
: std::map<TA, std::map<TAC, Tensor<Tdata>>>{};
7274

75+
Symmetry_Filter<TA, Tcell, Ndim, Tdata> symmetry_filter(this->period, irreducible_sector);
76+
7377
omp_lock_t lock_Ds_result_add;
7478
omp_init_lock(&lock_Ds_result_add);
7579

@@ -110,7 +114,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
110114

111115
for(const TAC &Aa2 : list_Aa2)
112116
{
113-
std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
117+
if (!symmetry_filter.is_I_in_irreducible_sector(Aa2.first)) continue;
118+
119+
std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
114120

115121
#pragma omp for schedule(dynamic) nowait
116122
for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
@@ -138,6 +144,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
138144
// D_result = D_mul * D_b
139145
for(const TAC &Ab2 : list_Ab2)
140146
{
147+
// if ((Ab2-Aa2) exceeds the irreducible sector) continue (known Ab01)
148+
if (!symmetry_filter.in_irreducible_sector(Aa2, Ab2)) continue;
149+
141150
const Tensor<Tdata> D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
142151
if(D_b.empty()) continue;
143152
// a2b2 = a2b0b1 * b0b1b2
@@ -232,7 +241,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
232241

233242
for(const TAC &Ab01 : list_Ab01)
234243
{
235-
std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
244+
if (!symmetry_filter.is_J_in_irreducible_sector(Ab01.first)) continue;
245+
246+
std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
236247

237248
#pragma omp for schedule(dynamic) nowait
238249
for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
@@ -258,7 +269,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
258269
// D_result = D_mul * D_a * D_a0b0
259270
for(const TAC &Aa2 : list_Aa2)
260271
{
261-
const Tensor<Tdata> &D_a_transpose = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
272+
// if ((Ab01-Aa2) exceeds the irreducible sector) continue
273+
if (!symmetry_filter.in_irreducible_sector(Aa2, Ab01)) continue;
274+
275+
const Tensor<Tdata>& D_a_transpose = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
262276
if(D_a_transpose.empty()) continue;
263277
// b1a1a0 = b0b1a1 * a0b0
264278
const Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1x2y0_ax1x2_y0a(D_mul, D_a0b0);
@@ -474,7 +488,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
474488

475489
for(const TA &Aa01 : list_Aa01)
476490
{
477-
std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
491+
if (!symmetry_filter.is_I_in_irreducible_sector(Aa01)) continue;
492+
493+
std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
478494

479495
#pragma omp for schedule(dynamic) nowait
480496
for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
@@ -500,8 +516,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
500516
// D_result = D_mul * D_a0b0 * D_b
501517
for(const TAC &Ab2 : list_Ab2)
502518
{
519+
// if ((Ab2-Aa01) exceeds the irreducible sector) continue
520+
if (!symmetry_filter.in_irreducible_sector(Aa01, Ab2)) continue;
521+
503522
const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
504-
if(D_b.empty()) continue;
523+
if (D_b.empty()) continue;
505524

506525
// b0b1a1 = a0b0 * b1a1a0
507526
const Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1y0y1_ax1_y0y1a(D_a0b0, D_mul);
@@ -720,7 +739,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
720739

721740
for(const TA &Aa01 : list_Aa01)
722741
{
723-
std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
742+
if (!symmetry_filter.is_I_in_irreducible_sector(Aa01)) continue;
743+
744+
std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
724745

725746
#pragma omp for schedule(dynamic) nowait
726747
for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
@@ -744,7 +765,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
744765
// D_result = D_mul * D_a0b0 * D_b
745766
for(const TAC &Ab01 : list_Ab01)
746767
{
747-
const Tensor<Tdata> &D_b_transpose = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
768+
// if ((Ab01-Aa01) exceeds the irreducible sector) continu
769+
if (!symmetry_filter.in_irreducible_sector(Aa01, Ab01)) continue;
770+
771+
const Tensor<Tdata>& D_b_transpose = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{ Ab2.first, (Ab2.second - Ab01.second) % this->period });
748772
if(D_b_transpose.empty()) continue;
749773
const Tensor<Tdata> D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
750774
if(D_a0b0.empty()) continue;

‎include/RI/ri/LRI.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,8 @@ class LRI
5959
void cal_loop3(
6060
const std::vector<Label::ab_ab> &labels,
6161
std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
62-
const double fac_add_Ds = 1.0);
62+
const double fac_add_Ds = 1.0,
63+
const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector = {});
6364

6465
public:
6566
std::shared_ptr<Parallel_LRI<TA,Tcell,Ndim,Tdata>>

‎include/RI/symmetry/Symmetry_Filter.h‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
#include <array>
2+
#include <map>
3+
#include <set>
4+
#include <tuple>
5+
#define NO_SEC_RETURN_TRUE if(this->irreducible_sector_.empty()) return true;
6+
#include "../global/Array_Operator.h"
7+
namespace RI
8+
{
9+
using namespace Array_Operator;
10+
template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
11+
class Symmetry_Filter
12+
{
13+
using TC = std::array<Tcell, Ndim>;
14+
using TAC = std::pair<TA, TC>;
15+
16+
using TIJ = std::pair<TA, TA>;
17+
using TIJR = std::pair<TIJ, TC>;
18+
using Tsec = std::map<TIJ, std::set<TC>>;
19+
public:
20+
Symmetry_Filter(const TC& period_in, const Tsec& irsec)
21+
:period(period_in), irreducible_sector_(irsec) {}
22+
bool in_irreducible_sector(const TA& Aa, const TAC& Ab) const
23+
{
24+
NO_SEC_RETURN_TRUE;
25+
const TIJ& ap = { Aa, Ab.first };
26+
if (irreducible_sector_.find(ap) != irreducible_sector_.end())
27+
if (irreducible_sector_.at(ap).find(Ab.second % this->period) != irreducible_sector_.at(ap).end())
28+
return true;
29+
return false;
30+
}
31+
bool in_irreducible_sector(const TAC& Aa, const TAC& Ab) const
32+
{
33+
NO_SEC_RETURN_TRUE;
34+
const TC dR = (Ab.second - Aa.second) % this->period;
35+
const std::pair<TA, TA> ap = { Aa.first, Ab.first };
36+
if (irreducible_sector_.find(ap) != irreducible_sector_.end())
37+
if (irreducible_sector_.at(ap).find(dR) != irreducible_sector_.at(ap).end())
38+
return true;
39+
return false;
40+
}
41+
bool is_I_in_irreducible_sector(const TA& Aa) const
42+
{
43+
NO_SEC_RETURN_TRUE;
44+
for (const auto& apRs : irreducible_sector_)
45+
if (apRs.first.first == Aa)return true;
46+
return false;
47+
}
48+
bool is_J_in_irreducible_sector(const TA& Ab) const
49+
{
50+
NO_SEC_RETURN_TRUE;
51+
for (const auto& apRs : irreducible_sector_)
52+
if (apRs.first.second == Ab)return true;
53+
return false;
54+
}
55+
TIJR get_IJR(const TA& I, const TAC& J) const
56+
{
57+
return { {I,J.first}, J.second % this->period };
58+
}
59+
TIJR get_IJR(const TAC& I, const TAC& J) const
60+
{
61+
return { {I.first,J.first}, (J.second - I.second) % this->period };
62+
}
63+
private:
64+
const Tsec& irreducible_sector_;
65+
const TC& period;
66+
};
67+
68+
}

0 commit comments

Comments
 (0)