Symmetry: Irreducible sector selection based on cal_loop3

include/RI/global/Tensor.h

@@ -45,6 +45,7 @@ class Tensor
 	inline T& operator() (const std::size_t i0, const std::size_t i1, const std::size_t i2, const std::size_t i3) const;
 
 	Tensor transpose() const;
+	Tensor dagger() const;
 
 	// ||d||_p = (|d_1|^p+|d_2|^p+...)^{1/p}
 	// if(p==std::numeric_limits<double>::max())    ||d||_max = max_i |d_i|

include/RI/global/Tensor.hpp

@@ -178,6 +178,17 @@ Tensor<T> Tensor<T>::transpose() const
 	return t;
 }
 
+template<typename T>
+Tensor<T> Tensor<T>::dagger() const
+{
+	assert(this->shape.size() == 2);
+	Tensor<T> t({ this->shape[1], this->shape[0] });
+	for (std::size_t i0 = 0; i0 < this->shape[0]; ++i0)
+		for (std::size_t i1 = 0; i1 < this->shape[1]; ++i1)
+			t(i1, i0) = std::conj((*this)(i0, i1));
+	return t;
+}
+
 template<typename T>
 Global_Func::To_Real_t<T> Tensor<T>::norm(const double p) const
 {

include/RI/physics/Exx.h

@@ -68,7 +68,9 @@ class Exx
 		const std::string &save_name_suffix="");
 
 	void cal_Hs(
-		const std::array<std::string,3> &save_names_suffix={"","",""});		// "Cs","Vs","Ds"
+		const std::array<std::string, 3>& save_names_suffix = { "","","" },		// "Cs","Vs","Ds"
+		const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector = {},
+		const bool if_cal_energy = true);
 	void cal_force(
 		const std::array<std::string,5> &save_names_suffix={"","","","",""});	// "Cs","Vs","Ds","dCs","dVs"
 	void cal_stress(

include/RI/physics/Exx.hpp

@@ -177,7 +177,9 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_dVRs(
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
-	const std::array<std::string,3> &save_names_suffix)						// "Cs","Vs","Ds"
+	const std::array<std::string, 3>& save_names_suffix,						// "Cs","Vs","Ds"
+	const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector,
+	const bool if_cal_energy)
 {
 	using namespace Map_Operator;
 
@@ -208,9 +210,11 @@ void Exx<TA,Tcell,Ndim,Tdata>::cal_Hs(
 		 Label::ab_ab::a0b0_a1b2,
 		 Label::ab_ab::a0b0_a2b1,
 		 Label::ab_ab::a0b0_a2b2},
-		this->Hs);
+		this->Hs,
+		1.0/*fac_add_Ds*/,
+		irreducible_sector);
 
-	//if()
+	if (if_cal_energy)
 		this->energy = this->post_2D.cal_energy(
 			this->post_2D.saves["Ds_"+save_names_suffix[2]],
 			this->post_2D.set_tensors_map2(this->Hs) );

include/RI/ri/LRI-cal_loop3.hpp

@@ -9,7 +9,7 @@
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Tensor_Multiply.h"
-
+#include "../symmetry/Symmetry_Filter.h"
 #include <omp.h>
 #ifdef __MKL_RI
 #include <mkl_service.h>
@@ -22,7 +22,8 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 	const std::vector<Label::ab_ab> &labels,
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-	const double fac_add_Ds)
+	const double fac_add_Ds,
+	const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector)
 {
 	using namespace Array_Operator;
 
@@ -70,6 +71,8 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 		? LRI_Cal_Aux::cal_Ds_transpose( data_wrapper(Label::ab::b).Ds_ab )
 		: std::map<TA, std::map<TAC, Tensor<Tdata>>>{};
 
+	Symmetry_Filter<TA, Tcell, Ndim, Tdata> symmetry_filter(this->period, irreducible_sector);
+
 	omp_lock_t lock_Ds_result_add;
 	omp_init_lock(&lock_Ds_result_add);
 
@@ -110,7 +113,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 
 					for(const TAC &Aa2 : list_Aa2)
 					{
-						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
+						if (!symmetry_filter.is_I_in_irreducible_sector(Aa2.first)) continue;
+
+						std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
@@ -138,6 +143,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 							// D_result = D_mul * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
+								// if ((Ab2-Aa2) exceeds the irreducible sector) continue (known Ab01)
+								if (!symmetry_filter.in_irreducible_sector(Aa2, Ab2))	continue;
+
 								const Tensor<Tdata> D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								// a2b2 = a2b0b1 * b0b1b2
@@ -232,7 +240,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
+						if (!symmetry_filter.is_J_in_irreducible_sector(Ab01.first)) continue;
+
+						std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
@@ -258,7 +268,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 							// D_result = D_mul * D_a * D_a0b0
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								const Tensor<Tdata> &D_a_transpose = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
+								// if ((Ab01-Aa2) exceeds the irreducible sector) continue
+								if (!symmetry_filter.in_irreducible_sector(Aa2, Ab01))	continue;
+
+								const Tensor<Tdata>& D_a_transpose = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a_transpose.empty())	continue;
 								// b1a1a0 = b0b1a1 * a0b0
 								const Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1x2y0_ax1x2_y0a(D_mul, D_a0b0);
@@ -474,7 +487,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
+						if (!symmetry_filter.is_I_in_irreducible_sector(Aa01)) continue;
+
+						std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
@@ -500,8 +515,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
+								// if ((Ab2-Aa01) exceeds the irreducible sector) continue
+								if (!symmetry_filter.in_irreducible_sector(Aa01, Ab2))	continue;
+
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
-								if(D_b.empty())	continue;
+								if (D_b.empty())	continue;
 
 								// b0b1a1 = a0b0 * b1a1a0
 								const Tensor<Tdata> D_tmp2 = Tensor_Multiply::x1y0y1_ax1_y0y1a(D_a0b0, D_mul);
@@ -720,7 +738,9 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
+						if (!symmetry_filter.is_I_in_irreducible_sector(Aa01))	continue;
+
+						std::map<TAC, Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
@@ -744,7 +764,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								const Tensor<Tdata> &D_b_transpose = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
+								// if ((Ab01-Aa01) exceeds the irreducible sector) continu
+								if (!symmetry_filter.in_irreducible_sector(Aa01, Ab01))	continue;
+
+								const Tensor<Tdata>& D_b_transpose = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{ Ab2.first, (Ab2.second - Ab01.second) % this->period });
 								if(D_b_transpose.empty())	continue;
 								const Tensor<Tdata> D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 								if(D_a0b0.empty())	continue;

include/RI/ri/LRI.h

@@ -59,7 +59,8 @@ class LRI
 	void cal_loop3(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const double fac_add_Ds = 1.0);
+		const double fac_add_Ds = 1.0,
+		const std::map<std::pair<TA, TA>, std::set<TC>>& irreducible_sector = {});
 
 public:
 	std::shared_ptr<Parallel_LRI<TA,Tcell,Ndim,Tdata>>

include/RI/symmetry/Symmetry_Filter.h

@@ -0,0 +1,68 @@
+#include <array>
+#include <map>
+#include <set>
+#include <tuple>
+#define NO_SEC_RETURN_TRUE if(this->irreducible_sector_.empty()) return true;
+#include "../global/Array_Operator.h"
+namespace RI
+{
+	using namespace Array_Operator;
+	template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
+	class Symmetry_Filter
+	{
+		using TC = std::array<Tcell, Ndim>;
+		using TAC = std::pair<TA, TC>;
+
+		using TIJ = std::pair<TA, TA>;
+		using TIJR = std::pair<TIJ, TC>;
+		using Tsec = std::map<TIJ, std::set<TC>>;
+	public:
+		Symmetry_Filter(const TC& period_in, const Tsec& irsec)
+			:period(period_in), irreducible_sector_(irsec) {}
+		bool in_irreducible_sector(const TA& Aa, const TAC& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			const TIJ& ap = { Aa, Ab.first };
+			if (irreducible_sector_.find(ap) != irreducible_sector_.end())
+				if (irreducible_sector_.at(ap).find(Ab.second % this->period) != irreducible_sector_.at(ap).end())
+					return true;
+			return false;
+		}
+		bool in_irreducible_sector(const TAC& Aa, const TAC& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			const TC dR = (Ab.second - Aa.second) % this->period;
+			const std::pair<TA, TA> ap = { Aa.first, Ab.first };
+			if (irreducible_sector_.find(ap) != irreducible_sector_.end())
+				if (irreducible_sector_.at(ap).find(dR) != irreducible_sector_.at(ap).end())
+					return true;
+			return false;
+		}
+		bool is_I_in_irreducible_sector(const TA& Aa) const
+		{
+			NO_SEC_RETURN_TRUE;
+			for (const auto& apRs : irreducible_sector_)
+				if (apRs.first.first == Aa)return true;
+			return false;
+		}
+		bool is_J_in_irreducible_sector(const TA& Ab) const
+		{
+			NO_SEC_RETURN_TRUE;
+			for (const auto& apRs : irreducible_sector_)
+				if (apRs.first.second == Ab)return true;
+			return false;
+		}
+		TIJR get_IJR(const TA& I, const TAC& J) const
+		{
+			return { {I,J.first}, J.second % this->period };
+		}
+		TIJR get_IJR(const TAC& I, const TAC& J) const
+		{
+			return { {I.first,J.first}, (J.second - I.second) % this->period };
+		}
+	private:
+		const Tsec& irreducible_sector_;
+		const TC& period;
+	};
+
+}

include/RI/symmetry/Symmetry_Rotation.h

@@ -0,0 +1,55 @@
+#include <RI/global/Tensor.h>
+namespace RI
+{
+	namespace Sym
+	{
+		template <typename T>
+		inline void T1_HR(T* TA, const T* A, const Tensor<T>& T1, const int& n2)
+		{
+			// C' = T1^\dagger * C
+			const int& n1 = T1.shape[0];
+			Blas_Interface::gemm('C', 'N', n1, n2, n1,
+				T(1), T1.ptr(), n1, A, n2, T(0), TA, n2);
+			// zgemm_(&notrans, &dagger, &n12, &nabf, &nabf,
+			// 	&alpha, A, &n12, T1.ptr(), &nabf, &beta, TA, &n12);
+		}
+
+		template <typename T>
+		inline void T1_HR_T2(T* TAT, const T* A, const Tensor<T>& T1, const Tensor<T>& T2)
+		{
+			// H' = T1^\dagger * H * T2
+			const int& n2 = T2.shape[0], & n1 = T1.shape[0];
+			const RI::Shape_Vector& shape = { static_cast<size_t>(n1),static_cast<size_t>(n2) };
+			RI::Tensor<T> AT2(shape);
+			Blas_Interface::gemm('N', 'N', n1, n2, n2,
+				T(1), A, n2, T2.ptr(), n2, T(0), AT2.ptr(), n2);
+			Blas_Interface::gemm('C', 'N', n1, n2, n1,
+				T(1), T1.ptr(), n1, AT2.ptr(), n2, T(0), TAT, n2);
+			// col-major version
+			// zgemm_(&notrans, &notrans, &n2, &n1, &n2,
+			// 	&alpha, T2.ptr(), &n2, A, &n2, &beta, AT2.ptr(), &n2);
+			// zgemm_(&notrans, &dagger, &n2, &n1, &n1,
+			// 	&alpha, AT2.ptr(), &n2, T1.ptr(), &n1, &beta, TAT, &n2);
+		}
+
+		template<typename T>
+		inline void T1_DR_T2(T* TAT, const T* A, const Tensor<T>& T1, const Tensor<T>& T2)
+		{
+			// D' = T1^T * D * T2^*  = T1^T * [T2^\dagger * D^T]^T
+			const int& n2 = T2.shape[0], & n1 = T1.shape[0];
+			const RI::Shape_Vector& shape = { static_cast<size_t>(n1),static_cast<size_t>(n2) };
+			RI::Tensor<T> AT2(shape);
+			BlasConnector::gemm('C', 'T', n2, n1, n2,
+				T(1), T2.ptr(), n2, A, n2, T(0), AT2.ptr(), n1);
+			BlasConnector::gemm('T', 'T', n1, n2, n1,
+				T(1), T1.ptr(), n1, AT2.ptr(), n1, T(0), TAT, n2);
+			// col-major version
+			// zgemm_(&transpose, &dagger, &nw1, &nw2, &nw2,
+			// 	&alpha, A, &nw2, T2.ptr(), &nw2, &beta, AT2.ptr(), &nw1);
+			// zgemm_(&transpose, &transpose, &nw2, &nw1, &nw1,
+			// 	&alpha, AT2.ptr(), &nw1, T1.ptr(), &nw1, &beta, TAT, &nw2);
+		}
+
+	}
+
+}