diff --git a/include/RI/global/Shared_Vector.h b/include/RI/global/Shared_Vector.h index 2041638..19b469f 100644 --- a/include/RI/global/Shared_Vector.h +++ b/include/RI/global/Shared_Vector.h @@ -28,6 +28,18 @@ class Shape_Vector for(auto ptr_in=v_in.begin(); ptr_in& v_in) + :size_(v_in.size()) + { + assert(v_in.size() <= sizeof(v) / sizeof(*v)); + for (std::size_t i = 0;i < size_;++i) this->v[i] = v_in[i]; + } + Shape_Vector(std::vector&& v_in) + :size_(v_in.size()) + { + assert(v_in.size() <= sizeof(v) / sizeof(*v)); + for (std::size_t i = 0;i < size_;++i) this->v[i] = v_in[i]; + } const std::size_t* begin() const noexcept { return this->v; } const std::size_t* end() const noexcept { return this->v+size_; } diff --git a/include/RI/global/Tensor.h b/include/RI/global/Tensor.h index bd840c9..3eeef09 100644 --- a/include/RI/global/Tensor.h +++ b/include/RI/global/Tensor.h @@ -47,6 +47,11 @@ class Tensor Tensor transpose() const; Tensor dagger() const; + /// permute from the input index order to {0, 1, 2, ..., N} + Tensor permute_from(const std::vector& order) const; + /// permute from {0, 1, 2, ..., N } to the input new index order + // Tensor permute_to(const std::vector) const; + // ||d||_p = (|d_1|^p+|d_2|^p+...)^{1/p} // if(p==std::numeric_limits::max()) ||d||_max = max_i |d_i| Global_Func::To_Real_t norm(const double p) const; diff --git a/include/RI/global/Tensor.hpp b/include/RI/global/Tensor.hpp index a0adae2..2508720 100644 --- a/include/RI/global/Tensor.hpp +++ b/include/RI/global/Tensor.hpp @@ -324,4 +324,74 @@ std::array,N2>,N1>,N0> to_array(const Ten return a; } +// check whether an array is a permutation of {from, from+1, ..., to} +inline bool is_permutation(const std::vector& arr, const std::size_t from, const std::size_t to) +{ + std::vector seen(arr.size(), false); + if (arr.size() != to - from + 1) return false; + for (std::size_t i = 0; i < arr.size(); ++i) { + if (arr[i] < from || arr[i] > to) return false; + const std::size_t idx = arr[i] - from; + if (seen[idx]) return false; + seen[idx] = true; + } + return true; +} + +template +Tensor Tensor::permute_from(const std::vector& order) const +{ + const std::size_t N = order.size(); + assert(this->shape.size() == N); + assert(is_permutation(order, 0, N-1)); + Shape_Vector new_shape(this->shape); + for (std::size_t d = 0; d < N; ++d) + new_shape[d] = this->shape[order[d]]; + + Tensor t(new_shape); + + /* N-permutation + for (std::size_t i0 = 0; i0 < new_shape[0], ++i0) + for (std::size_t i1 = 0, i1 < new_shape[1], ++i1) + ... + for (std::size_t iN = 0, iN < new_shape[N], ++iN) + t(i0, i1, ..., iN) = this->operator()(order(i0), order(i1), ..., order(N)) + */ + std::vector strides(N); + strides[N - 1] = 1; + for (int d = N - 2; d >= 0; --d) + strides[d] = strides[d + 1] * this->shape[d + 1]; + + std::vector new_strides(N); + new_strides[N - 1] = 1; + for (int d = N - 2; d >= 0; --d) + new_strides[d] = new_strides[d + 1] * new_shape[d + 1]; + + + std::vector new_idx(N, 0); + std::vector idx(N, 0); + for (std::size_t i = 0; i < t.get_shape_all(); ++i) // new 1D-index + { + // new ND-index + std::size_t tmp = i; + for (std::size_t d = 0; d < N; ++d) + { + new_idx[d] = tmp / new_strides[d]; + tmp %= new_strides[d]; + } + // old ND-index + for (std::size_t d = 0; d < N; ++d) + idx[order[d]] = new_idx[d]; + + // old 1D-index + std::size_t i0 = 0; + for (std::size_t d = 0; d < N; ++d) + i0 += idx[d] * strides[d]; + + (*t.data)[i] = (*this->data)[i0]; + } + + return t; +} + } \ No newline at end of file diff --git a/include/RI/global/Tensor_Multiply-32.hpp b/include/RI/global/Tensor_Multiply-32.hpp index a46c3f9..005ed0f 100644 --- a/include/RI/global/Tensor_Multiply-32.hpp +++ b/include/RI/global/Tensor_Multiply-32.hpp @@ -80,6 +80,28 @@ namespace Tensor_Multiply return Txy; } + // Txy(x0,x1,..., xM) = Tx(x0, x1,x2, ..., xM, y0, y1, ..., yN) * Vy(y0, y1, ..., yN) + template + Tensor gemv(const Tensor& Tx, const Tensor& Vy) + { + assert(Tx.shape.size() >= Vy.shape.size()); + const std::size_t ny = Vy.get_shape_all(); + assert(Tx.get_shape_all() % ny == 0); + const std::size_t dim = Tx.shape.size() - Vy.shape.size(); + std::vector shape_vector; + if (dim == 0) + shape_vector.push_back(1); + else + for (int d = 0; d < dim;++d) + shape_vector.push_back(Tx.shape[d]); + Tensor Txy(shape_vector); + Blas_Interface::gemv( + 'N', Txy.get_shape_all(), ny, + Tdata(1.0), Tx.ptr(), ny, Vy.ptr(), 1, + Tdata(0.0), Txy.ptr(), 1); + return Txy; + } + } } diff --git a/include/RI/physics/Exx.h b/include/RI/physics/Exx.h index e7457d1..bdb8ee9 100644 --- a/include/RI/physics/Exx.h +++ b/include/RI/physics/Exx.h @@ -29,6 +29,10 @@ class Exx constexpr static std::size_t Npos = Ndim; // tmp using Tatom_pos = std::array; // tmp + Exx(const std::string method = "loop3") { this->method = method; } + Exx(Exx&&) noexcept = default; + Exx& operator=(Exx&&) noexcept = default; + void set_parallel( const MPI_Comm &mpi_comm, const std::map &atoms_pos, @@ -96,6 +100,7 @@ class Exx void free_dVs(const std::string &save_name_suffix=""); void free_dCRs(const std::string &save_name_suffix=""); void free_dVRs(const std::string &save_name_suffix=""); + void free_cvc() { this->cvc_.clear(); } public: LRI lri; @@ -122,6 +127,9 @@ class Exx }; Flag_Save_Result flag_save_result; + std::string method = "loop3"; // "loop3", "cvc" + std::map>>>> cvc_; + std::map>>>>& get_cvc() const { return this->cvc_; } }; } diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp index 2bb8be3..a1073eb 100644 --- a/include/RI/physics/Exx.hpp +++ b/include/RI/physics/Exx.hpp @@ -266,12 +266,21 @@ void Exx::cal_Hs( if(!this->flag_finish.Ds_delta) this->Hs.clear(); - this->lri.cal_loop3( - {Label::ab_ab::a0b0_a1b1, - Label::ab_ab::a0b0_a1b2, - Label::ab_ab::a0b0_a2b1, - Label::ab_ab::a0b0_a2b2}, - this->Hs); + if (this->method == "loop3") + this->lri.cal_loop3( + { Label::ab_ab::a0b0_a1b1, + Label::ab_ab::a0b0_a1b2, + Label::ab_ab::a0b0_a2b1, + Label::ab_ab::a0b0_a2b2 }, + this->Hs); + else if (this->method == "cvc") + { + if (this->cvc_.empty()) + this->cvc_ = this->lri.cal_cvc(); + this->Hs = this->lri.constract_cvc_ds(this->cvc_); + } + else + assert(false && "Unknown method in Exx::cal_Hs"); //if() this->energy = this->post_2D.cal_energy( diff --git a/include/RI/physics/LR.h b/include/RI/physics/LR.h new file mode 100644 index 0000000..3ea19c9 --- /dev/null +++ b/include/RI/physics/LR.h @@ -0,0 +1,26 @@ +#pragma once +#include "./Exx.h" + +namespace RI +{ +// Nothing different from Exx, +// Except the two density matrices can be different +template +class LR : public Exx +{ +public: + LR(const std::string method = "cvc") { this->method = method; } // Exx default mehtod: "loop3" + LR(Exx&& exx) : Exx::Exx(std::move(exx)) {} + using Exx::cal_force; + // New function: cal_force with two different density matrices + void cal_force(const std::map>, Tensor>>& Ds_left, + const std::array& save_names_suffix = { "","","","","" }) // "Cs","Vs","Ds","dCs","dVs" + { + // The only difference from Exx::cal_force: save Ds_left if not empty + if (!Ds_left.empty()) + this->post_2D.saves["Ds_" + save_names_suffix[2]] = this->post_2D.set_tensors_map2(Ds_left); + this->cal_force(save_names_suffix); + } +}; + +} \ No newline at end of file diff --git a/include/RI/ri/LRI-cal_cvc.hpp b/include/RI/ri/LRI-cal_cvc.hpp new file mode 100644 index 0000000..9b17fa7 --- /dev/null +++ b/include/RI/ri/LRI-cal_cvc.hpp @@ -0,0 +1,315 @@ +// =================== +// Author: Peize Lin +// date: 2023.08.02 +// =================== + +#pragma once + +#include "LRI.h" +#include "LRI_Cal_Aux.h" +#include "../global/Array_Operator.h" +#include "../global/Tensor_Multiply.h" +#include +#include +#ifdef __MKL_RI +#include +#endif + +namespace RI +{ + +template +std::map>, std::map>, Tensor>>>> +LRI::cal_cvc() +{ + using namespace Array_Operator; + + const Data_Pack_Wrapper data_wrapper(this->data_pool, this->data_ab_name); + const LRI_Cal_Tools tools(this->period, this->data_pool, this->data_ab_name); + + #ifdef __MKL_RI + const std::size_t mkl_threads = mkl_get_max_threads(); + mkl_set_num_threads(1); + #endif + + std::map>>>> cvc; + /// ??? + std::map lock_cvc_result_add_map = LRI_Cal_Aux::init_lock_result({ Label::ab_ab::a0b0_a2b2 }, this->parallel->list_A, cvc); + + #pragma omp parallel + { + std::map>>>> cvc_thread; +/* pseudo code +for I + for J + for K + [CV]_{IK,J}=C^I_{IK}V_{IJ} + for L + [CVC]_{IJKL}+=[CV]_{IK,J} C^J_{JL} #1 + for L + for K + [CV]_{IK,L}=C^I_{IK}V_{IL} + for J + [CVC]_{IJKL}+=[CV]_{IK,L} C^L_{LJ} #2 +for K + for J + for I + [CV]_{KI,J}=C^K_{KI}V_{KJ} + for L + [CVC]_{IJKL}+=[CV]_{KI,J}C^J_{JL} #3 + for L + for I + [CV]_{KI,L}=C^K_{KI}V_{KL} + for J + [CVC]_{IJKL}+=[CV]_{KI,L}C^L_{LJ} #4 +*/ + // Main Problems: + // How to filter list_I, list_Aa02, list_K, list_L only using the range of Vs? + // How to deal with thread lock? + // 转置无法通过改变乘法顺序避免： 输出KL连续，输入KL分别在两个C上 + // for debug + auto print_a = [](const std::vector& vec, const std::string name) -> void + { + std::cout << name << ": "; + for (auto& v : vec) { std::cout << v << " "; } + std::cout << std::endl; + }; + auto print_ac = [](const std::vector& vec, const std::string name) -> void + { + std::cout << name << ": "; + for (auto& v : vec) { std::cout << v.first << " "; } + std::cout << std::endl; + }; + const std::vector list_I0 = LRI_Cal_Aux::filter_list_map(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).a01, data_wrapper(Label::ab::a).Ds_ab); + const std::vector list_J0 = LRI_Cal_Aux::filter_list_map(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).b01, data_wrapper(Label::ab::b).Ds_ab); + const std::vector list_K = LRI_Cal_Aux::filter_list_set(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).a2, data_wrapper(Label::ab::a).index_Ds_ab[0]); + const std::vector list_L = LRI_Cal_Aux::filter_list_set(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).b2, data_wrapper(Label::ab::b).index_Ds_ab[0]); + // filter Cs by the range of Vs + const std::vector list_I = LRI_Cal_Aux::filter_list_map(list_I0, data_wrapper(Label::ab::a0b0).Ds_ab); + const std::vector list_J = LRI_Cal_Aux::filter_list_set(list_J0, data_wrapper(Label::ab::a0b0).index_Ds_ab[0]); + // print_a(list_I0, "list_I0"); + // print_ac(list_J0, "list_J0"); + // print_ac(list_K, "list_K"); + // print_ac(list_L, "list_L"); + // print_a(list_I, "list_I"); + // print_ac(list_J, "list_J"); + + // allocate + // for (TA I : list_I) + // for (TAC J : list_J) + // for (TAC K : list_K) + // { + // const Tensor& C_I_IK = tools.get_Ds_ab(Label::ab::a, I, K); + // const int& ni = C_I_IK.shape[1]; + // const int& nk = C_I_IK.shape[2]; + // for (TAC L : list_L) + // { + // const TC R_KL = (L.second - K.second) % period; + // const Tensor& C_J_JL = tools.get_Ds_ab(Label::ab::b, J, L); + // const int& nj = C_J_JL.shape[1]; + // const int& nl = C_J_JL.shape[2]; + // cvc_thread[I][J][K.first][{L.first, R_KL}] = Tensor({ ni, nj, nk, nl }); + // } + // } + + // calculate +#pragma omp for schedule(static) collapse(2) nowait + for (TA I : list_I) + { + if (this->filter_atom->filter_for1(Label::ab_ab::a0b0_a2b2, I)) continue; // restrict I in the irreducible sector + auto& cvc_thread_I = cvc_thread[I]; + for (TAC J : list_J) // term 1 + { + if (this->filter_atom->filter_for32(Label::ab_ab::a0b0_a2b2, I, J, J)) continue; // restrict (I, J) in the irreducible sector + const Tensor& V_IJ = tools.get_Ds_ab(Label::ab::a0b0, I, J); + if (V_IJ.empty()) continue; + // symmetry: check if (I, J) in irreducible sector + auto& cvc_thread_IJ = cvc_thread[I][J]; + for (TAC K : list_K) + { + const Tensor& C_I_IK = tools.get_Ds_ab(Label::ab::a, I, K); + if (C_I_IK.empty()) continue; + // CV_{IK,J}=C^I_{IK}V_{IJ} + const Tensor CV_IK_J = Tensor_Multiply::x1x2y1_ax1x2_ay1(C_I_IK, V_IJ); + auto& cvc_thread_IJK = cvc_thread_IJ[K.first]; + for (TAC L : list_L) + { + const TC R_KL = (L.second - K.second) % period; + // std::array key = { I, J, K, L }; + // [CVC]_{IJKL}+=[CV]_{IK,J}C^J_{JL} + const Tensor& C_J_JL = tools.get_Ds_ab(Label::ab::b, J, L); + if (C_J_JL.empty()) continue; + //(ika) * (ajl) = (ikjl) -> (ijkl) + // cvc_thread[{I, J, K, L}][R_KL] += einsum("ika, ajl -> ijkl", CV_IK_J, C_J_JL); + // cvc_thread[I][J][K.first][{L.first, R_KL}] += Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_IK_J, C_J_JL).permute_from({ 0,2,1,3 }); + LRI_Cal_Aux::add_Ds( + Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_IK_J, C_J_JL).permute_from({ 0,2,1,3 }), + cvc_thread_IJK[{L.first, R_KL}]); + } + } + LRI_Cal_Aux::add_Ds_omp_try_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + } + } + LRI_Cal_Aux::add_Ds_omp_wait_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + +#pragma omp for schedule(static) collapse(2) nowait + for (TA I : list_I) + { + if (this->filter_atom->filter_for1(Label::ab_ab::a0b0_a2b2, I)) continue; + auto& cvc_thread_I = cvc_thread[I]; + for (TAC L : list_L) // term 2 + { + const Tensor& V_IL = tools.get_Ds_ab(Label::ab::a0b0, I, L); + if (V_IL.empty()) continue; + for (TAC K : list_K) + { + const TC& R_KL = (L.second - K.second) % period; + const Tensor& C_I_IK = tools.get_Ds_ab(Label::ab::a, I, K); + if (C_I_IK.empty()) continue; + // CV_{IK,L}: a1a2b0 = a0a1a2 * a0b0 + const Tensor CV_IK_L = Tensor_Multiply::x1x2y1_ax1x2_ay1(C_I_IK, V_IL); + for (TAC J : list_J) + { + if (this->filter_atom->filter_for32(Label::ab_ab::a0b0_a2b2, I, J, J)) continue; // restrict (I, J) in the irreducible sector + // [CVC]_{IJKL}+=[CV]_{IK,L}C^L_{LJ} + const Tensor& C_L_LJ = tools.get_Ds_ab(Label::ab::b, L, J); + if (C_L_LJ.empty()) continue; + // (ika) * (alj) = (iklj) -> (ijkl) + // cvc_thread[{I, J, K, L}][R_KL] += einsum("ika, alj -> ijkl", CV_IK_L, C_L_LJ); + // cvc_thread[I][J][K.first][{L.first, R_KL}] += Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_IK_L, C_L_LJ).permute_from({ 0,3,1,2 }); + LRI_Cal_Aux::add_Ds( + Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_IK_L, C_L_LJ).permute_from({ 0,3,1,2 }), + cvc_thread_I[J][K.first][{L.first, R_KL}]); + } + } + LRI_Cal_Aux::add_Ds_omp_try_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + } + } + LRI_Cal_Aux::add_Ds_omp_wait_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + +#pragma omp for schedule(static) collapse(2) nowait + for (TAC K : list_K) + { + for (TAC J : list_J) //term 3 + { + if (this->filter_atom->filter_for1(Label::ab_ab::a0b0_a1b2, J)) continue; // restrict J in the irreducible sector + const Tensor& V_KJ = tools.get_Ds_ab(Label::ab::a0b0, K, J); + if (V_KJ.empty()) continue; + for (TA I : list_I) + { + if (this->filter_atom->filter_for32(Label::ab_ab::a0b0_a2b2, I, J, J)) continue; // restrict (I, J) in the irreducible sector + const Tensor& C_K_KI = tools.get_Ds_ab(Label::ab::a, K, I); + if (C_K_KI.empty()) continue; + //[CV]_{KI,J}=C^K_{KI}V_{KJ} + const Tensor CV_KI_J = Tensor_Multiply::x1x2y1_ax1x2_ay1(C_K_KI, V_KJ); + auto& cvc_thread_IJK = cvc_thread[I][J][K.first]; + for (TAC L : list_L) + { + const TC& R_KL = (L.second - K.second) % period; + // [CVC]_{IJKL}+=[CV]_{KI,J}C^J_{JL} #(3) + const Tensor& C_J_JL = tools.get_Ds_ab(Label::ab::b, J, L); + if (C_J_JL.empty()) continue; + // (kia) * (ajl) = (kijl)->(ijkl) + // cvc_thread[{I, J, K, L}][R_KL] += einsum("kia, ajl -> ijkl", CV_KI_J, C_J_JL); + // cvc_thread[I][J][K.first][{L.first, R_KL}] += Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_KI_J, C_J_JL).permute_from({ 1,2,0,3 }); + LRI_Cal_Aux::add_Ds( + Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_KI_J, C_J_JL).permute_from({ 1,2,0,3 }), + cvc_thread_IJK[{L.first, R_KL}]); + } + } + LRI_Cal_Aux::add_Ds_omp_try_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + } + } + LRI_Cal_Aux::add_Ds_omp_wait_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + +#pragma omp for schedule(static) collapse(2) nowait + for (TAC K : list_K) + { + for (TAC L : list_L) // term 4 + { + const Tensor& V_KL = tools.get_Ds_ab(Label::ab::a0b0, K, L); + if (V_KL.empty()) continue; + const TC& R_KL = (L.second - K.second) % period; + for (TA I : list_I) + { + if (this->filter_atom->filter_for1(Label::ab_ab::a0b0_a2b2, I)) continue; // restrict I in the irreducible sector + const Tensor& C_K_KI = tools.get_Ds_ab(Label::ab::a, K, I); + if (C_K_KI.empty()) continue; + // CV_{IK,L}: a1a2b0 = a0a1a2 * a0b0 + const Tensor CV_KI_L = Tensor_Multiply::x1x2y1_ax1x2_ay1(C_K_KI, V_KL); + for (TAC J : list_J) + { + if (this->filter_atom->filter_for32(Label::ab_ab::a0b0_a2b2, I, J, J)) continue; // restrict (I, J) in the irreducible sector + // [CVC]_{IJKL}+=[CV]_{IK,L}C^L_{LJ} + const Tensor& C_L_LJ = tools.get_Ds_ab(Label::ab::b, L, J); + if (C_L_LJ.empty()) continue; + // (kia) * (alj) = (kilj) -> (ijkl) + // cvc_thread[{I, J, K, L}][R_KL] += einsum("kia, alj -> ijkl", CV_KI_L, C_L_LJ); + // cvc_thread[I][J][K.first][{L.first, R_KL}] += Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_KI_L, C_L_LJ).permute_from({ 1,3,0,2 }); + LRI_Cal_Aux::add_Ds( + Tensor_Multiply::x0x1y1y2_x0x1a_ay1y2(CV_KI_L, C_L_LJ).permute_from({ 1,3,0,2 }), + cvc_thread[I][J][K.first][{L.first, R_KL}]); + } + } + LRI_Cal_Aux::add_Ds_omp_try_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + } + } + LRI_Cal_Aux::add_Ds_omp_wait_map(cvc_thread, cvc, lock_cvc_result_add_map, 1.0); + } // end #pragma omp parallel + + LRI_Cal_Aux::destroy_lock_result(lock_cvc_result_add_map, cvc); + + #ifdef __MKL_RI + mkl_set_num_threads(mkl_threads); + #endif + + malloc_trim(0); + return cvc; +} // end LRI::cal_cvc + +template +std::map>, Tensor>> +LRI::constract_cvc_ds( + const std::map>>>>& cvc) +{ + const Data_Pack_Wrapper data_wrapper(this->data_pool, this->data_ab_name); + const LRI_Cal_Tools tools(this->period, this->data_pool, this->data_ab_name); + + const std::vector list_I = LRI_Cal_Aux::filter_list_map(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).a01, data_wrapper(Label::ab::a).Ds_ab); + const std::vector list_J = LRI_Cal_Aux::filter_list_map(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).b01, data_wrapper(Label::ab::b).Ds_ab); + const std::vector list_K0 = LRI_Cal_Aux::filter_list_set(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).a2, data_wrapper(Label::ab::a).index_Ds_ab[0]); + const std::vector list_L0 = LRI_Cal_Aux::filter_list_set(this->parallel->list_A.at(Label::Aab_Aab::a01b01_a2b2).b2, data_wrapper(Label::ab::b).index_Ds_ab[0]); + + const std::vector list_K = LRI_Cal_Aux::filter_list_map(list_K0, data_wrapper(Label::ab::a2b2).Ds_ab); + const std::vector list_L = LRI_Cal_Aux::filter_list_set(list_L0, data_wrapper(Label::ab::a2b2).index_Ds_ab[0]); + + std::map>> Hs; + for (auto& map_i : cvc) + { + const TA& I = map_i.first; + if (this->filter_atom->filter_for1(Label::ab_ab::a0b0_a2b2, I)) continue; // restrict I in the irreducible sector + for (auto& map_j : map_i.second) + { + const TAC& J = map_j.first; + if (this->filter_atom->filter_for32(Label::ab_ab::a0b0_a2b2, I, J, J)) continue; // restrict (I, J) in the irreducible sector + // init H_IJ + for (auto& map_k : map_j.second) + { + const TA& K = map_k.first; + for (auto& map_l : map_k.second) + { + const TAC& L = map_l.first; + const Tensor& CVC_IJKL = map_l.second; + const Tensor& D_KL = tools.get_Ds_ab(Label::ab::a2b2, K, L); + if (D_KL.empty()) continue; + // Hs[I][J] += Tensor_Multiply::gemv(CVC_IJKL, D_KL) + LRI_Cal_Aux::add_Ds(Tensor_Multiply::gemv(CVC_IJKL, D_KL), Hs[I][J]); + } + } + } + } + return Hs; +} + +} // end namespace RI + diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h index 5cf0c1b..00847bd 100644 --- a/include/RI/ri/LRI.h +++ b/include/RI/ri/LRI.h @@ -64,6 +64,10 @@ class LRI std::map>> &Ds_result, const double fac_add_Ds = 1.0); + std::map>>>> cal_cvc(); + std::map>> constract_cvc_ds( + const std::map>>>>& cvc); + public: std::shared_ptr> parallel = std::make_shared>(); @@ -94,3 +98,4 @@ class LRI #include "LRI.hpp" #include "LRI-set.hpp" #include "LRI-cal_loop3.hpp" +#include "LRI-cal_cvc.hpp" diff --git a/include/RI/ri/LRI_Cal_Aux.h b/include/RI/ri/LRI_Cal_Aux.h index 1708ed1..44dcdad 100644 --- a/include/RI/ri/LRI_Cal_Aux.h +++ b/include/RI/ri/LRI_Cal_Aux.h @@ -160,10 +160,10 @@ namespace LRI_Cal_Aux } } - template + template void add_Ds_omp_try_map( - std::map>> &Ds_result_thread, - std::map>> &Ds_result, + std::map>& Ds_result_thread, + std::map>& Ds_result, std::map &lock_Ds_result_add_map, const double &fac) { @@ -183,10 +183,10 @@ namespace LRI_Cal_Aux } } - template + template void add_Ds_omp_wait_map( - std::map>> &Ds_result_thread, - std::map>> &Ds_result, + std::map>& Ds_result_thread, + std::map>& Ds_result, std::map &lock_Ds_result_add_map, const double &fac) { @@ -338,11 +338,11 @@ namespace LRI_Cal_Aux return list_filter; } - template + template std::map init_lock_result( const std::vector &labels, const std::unordered_map> &list_A, - std::map>> &Ds_result) + std::map>& Ds_result) { std::map lock_Ds_result_add_map; for(const Label::ab_ab &label : labels) diff --git a/include/RI/ri/LRI_Cal_Tools.h b/include/RI/ri/LRI_Cal_Tools.h index e3df5e5..983e34b 100644 --- a/include/RI/ri/LRI_Cal_Tools.h +++ b/include/RI/ri/LRI_Cal_Tools.h @@ -58,7 +58,15 @@ class LRI_Cal_Tools *this->Ds_ab_ptr.at(label), Aa.first, TAC{Ab.first, (Ab.second-Aa.second)%this->period}); } - + inline const Tensor& get_Ds_ab( + const Label::ab& label, + const TAC& Aa, const TA& Ab) const + { + using namespace Array_Operator; + return Global_Func::find( + *this->Ds_ab_ptr.at(label), + Aa.first, TAC{ Ab, (-Aa.second) % this->period }); + } inline std::vector split_b01(const Label::ab_ab &label) const { switch(label)