From f5ea63504617f69bdda79f4e74ae306fcd749336 Mon Sep 17 00:00:00 2001 From: linpz Date: Wed, 5 Aug 2026 11:40:20 +0800 Subject: [PATCH 1/2] 1. add OpenMP in XC_Functional_Libxc::v_xc_libxc() 2. add OpenMP in RI_2D_Comm::split_m2D_ktoR() --- .../module_xc/xc_functional_libxc_vxc.cpp | 44 ++- source/module_ri/RI_2D_Comm.h | 32 ++- source/module_ri/RI_2D_Comm.hpp | 267 +++++++++++++----- 3 files changed, 262 insertions(+), 81 deletions(-) diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp index 2456abe5ba9..f69efc1b94c 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -10,7 +10,7 @@ #include "module_base/tool_title.h" #include - +#include #include std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 @@ -94,16 +94,46 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / switch( func.info->family ) { case XC_FAMILY_LDA: - // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, nrxx, rho.data(), - exc.data(), vrho.data() ); + { + constexpr int nr_batch_size = 1024; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, nr_batch_size) + #endif + for( int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size ) + { + const int ir_end = std::min(ir_start + nr_batch_size, nrxx); + const int nrxx_thread = ir_end - ir_start; + xc_lda_exc_vxc( + &func, + nrxx_thread, + rho.data() + ir_start * nspin, + exc.data() + ir_start, + vrho.data() + ir_start * nspin ); + } break; + } case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: - // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, nrxx, rho.data(), sigma.data(), - exc.data(), vrho.data(), vsigma.data() ); + { + constexpr int nr_batch_size = 1024; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, nr_batch_size) + #endif + for( int ir_start = 0; ir_start < nrxx; ir_start += nr_batch_size ) + { + const int ir_end = std::min(ir_start + nr_batch_size, nrxx); + const int nrxx_thread = ir_end - ir_start; + xc_gga_exc_vxc( + &func, + nrxx_thread, + rho.data() + ir_start * nspin, + sigma.data() + ir_start * ((1==nspin)?1:3), + exc.data() + ir_start, + vrho.data() + ir_start * nspin, + vsigma.data() + ir_start * ((1==nspin)?1:3) ); + } break; + } default: throw std::domain_error("func.info->family ="+std::to_string(func.info->family) +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); diff --git a/source/module_ri/RI_2D_Comm.h b/source/module_ri/RI_2D_Comm.h index c7d023e735d..ca7e2108701 100644 --- a/source/module_ri/RI_2D_Comm.h +++ b/source/module_ri/RI_2D_Comm.h @@ -29,14 +29,30 @@ namespace RI_2D_Comm using TAC = std::pair; //public: - template - extern std::vector>>> - split_m2D_ktoR(const UnitCell& ucell, - const K_Vectors& kv, - const std::vector& mks_2D, - const Parallel_2D& pv, - const int nspin, - const bool spgsym = false); + template + extern std::vector>>> split_m2D_ktoR( + const UnitCell& ucell, + const K_Vectors& kv, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin, + const bool spgsym = false); + + template + extern std::vector>>> split_m2D_ktoR_gamma( + const UnitCell& ucell, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin); + + template + extern std::vector>>> split_m2D_ktoR_k( + const UnitCell& ucell, + const K_Vectors& kv, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin, + const bool spgsym = false); // judge[is] = {s0, s1} extern std::vector, std::set>> diff --git a/source/module_ri/RI_2D_Comm.hpp b/source/module_ri/RI_2D_Comm.hpp index f40e6f59d4a..056d6aff6ef 100644 --- a/source/module_ri/RI_2D_Comm.hpp +++ b/source/module_ri/RI_2D_Comm.hpp @@ -30,93 +30,228 @@ inline RI::Tensor> tensor_conj(const RI::Tensor auto RI_2D_Comm::split_m2D_ktoR(const UnitCell& ucell, - const K_Vectors & kv, - const std::vector&mks_2D, - const Parallel_2D & pv, - const int nspin, + const K_Vectors & kv, + const std::vector&mks_2D, + const Parallel_2D & pv, + const int nspin, const bool spgsym) -> std::vector>>> { ModuleBase::TITLE("RI_2D_Comm","split_m2D_ktoR"); ModuleBase::timer::tick("RI_2D_Comm", "split_m2D_ktoR"); - const TC period = RI_Util::get_Born_vonKarmen_period(kv); + std::vector>>> mRs_a2D + = (period == TC{1, 1, 1}) + ? RI_2D_Comm::split_m2D_ktoR_gamma(ucell, mks_2D, pv, nspin) + : RI_2D_Comm::split_m2D_ktoR_k(ucell, kv, mks_2D, pv, nspin, spgsym); + ModuleBase::timer::tick("RI_2D_Comm", "split_m2D_ktoR"); + return mRs_a2D; +} + +template +auto RI_2D_Comm::split_m2D_ktoR_gamma(const UnitCell& ucell, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin) +-> std::vector>>> +{ + ModuleBase::TITLE("RI_2D_Comm","split_m2D_ktoR_gamma"); + ModuleBase::timer::tick("RI_2D_Comm", "split_m2D_ktoR_gamma"); + const std::map nspin_k = {{1,1}, {2,2}, {4,1}}; const double SPIN_multiple = std::map{ {1,0.5}, {2,1}, {4,1} }.at(nspin); // why? + const TC cell = {0, 0, 0}; std::vector>>> mRs_a2D(nspin); + + #ifdef _OPENMP + // pre-init all outer maps mRs_a2D[is_b][iat] to avoid concurrent std::map rebalancing + for (int is_b = 0; is_b < nspin; ++is_b) + for (int iat0 = 0; iat0 < ucell.nat; ++iat0) + mRs_a2D[is_b][iat0]; + + std::vector locks(ucell.nat); + for (auto& l : locks) + omp_init_lock(&l); + #endif + for (int is_k = 0; is_k < nspin_k.at(nspin); ++is_k) - { - const std::vector ik_list = RI_2D_Comm::get_ik_list(kv, is_k); - for(const TC &cell : RI_Util::get_Born_von_Karmen_cells(period)) - { - RI::Tensor mR_2D; - int ik_full = 0; - for (const int ik : ik_list) + { + using Tdata_m = typename Tmatrix::value_type; + RI::Tensor mk_2D + = RI_Util::Vector_to_Tensor(*mks_2D[is_k], pv.get_col_size(), pv.get_row_size()); + const Tdata_m frac = RI::Global_Func::convert(SPIN_multiple); + RI::Tensor mR_2D = RI::Global_Func::convert(mk_2D * frac); + + #ifdef _OPENMP + #pragma omp parallel for schedule(dynamic) + #endif + for (int iwt0_2D = 0; iwt0_2D != mR_2D.shape[0]; ++iwt0_2D) + { + const int iwt0 = ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) + ? pv.local2global_col(iwt0_2D) + : pv.local2global_row(iwt0_2D); + int iat0, iw0_b, is0_b; + std::tie(iat0, iw0_b, is0_b) = RI_2D_Comm::get_iat_iw_is_block(ucell, iwt0); + const int it0 = ucell.iat2it[iat0]; + for (int iwt1_2D = 0; iwt1_2D != mR_2D.shape[1]; ++iwt1_2D) { - auto set_mR_2D = [&mR_2D](auto&& mk_frac) { - if (mR_2D.empty()) { - mR_2D = RI::Global_Func::convert(mk_frac); - } else { - mR_2D - = mR_2D + RI::Global_Func::convert(mk_frac); - } - }; - using Tdata_m = typename Tmatrix::value_type; - if (!spgsym) + const int iwt1 = ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) + ? pv.local2global_row(iwt1_2D) + : pv.local2global_col(iwt1_2D); + int iat1, iw1_b, is1_b; + std::tie(iat1, iw1_b, is1_b) = RI_2D_Comm::get_iat_iw_is_block(ucell, iwt1); + const int it1 = ucell.iat2it[iat1]; + + const int is_b = RI_2D_Comm::get_is_block(is_k, is0_b, is1_b); + #ifdef _OPENMP + omp_set_lock(&locks[iat0]); + #endif + RI::Tensor& mR_a2D = mRs_a2D[is_b][iat0][{iat1, cell}]; + if (mR_a2D.empty()) { - RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik], pv.get_col_size(), pv.get_row_size()); - const Tdata_m frac = SPIN_multiple - * RI::Global_Func::convert(std::exp( - -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); - if (static_cast(std::round(SPIN_multiple * kv.wk[ik] * kv.get_nkstot_full())) == 2) - { set_mR_2D(mk_2D * (frac * 0.5) + tensor_conj(mk_2D * (frac * 0.5))); } - else { set_mR_2D(mk_2D * frac); } + mR_a2D = RI::Tensor( + {static_cast(ucell.atoms[it0].nw), + static_cast(ucell.atoms[it1].nw)}); } + mR_a2D(iw0_b, iw1_b) = mR_2D(iwt0_2D, iwt1_2D); + #ifdef _OPENMP + omp_unset_lock(&locks[iat0]); + #endif + } + } + } + + #ifdef _OPENMP + for (auto& l : locks) + omp_destroy_lock(&l); + + // prune empty inner maps created by pre-init + for (int is_b = 0; is_b < nspin; ++is_b) + for (auto it = mRs_a2D[is_b].begin(); it != mRs_a2D[is_b].end();) + { + if (it->second.empty()) + it = mRs_a2D[is_b].erase(it); else - { // traverse kstar, ik means ik_ibz - for (auto& isym_kvd : kv.kstars[ik % ik_list.size()]) + ++it; + } + #endif + + ModuleBase::timer::tick("RI_2D_Comm", "split_m2D_ktoR_gamma"); + return mRs_a2D; +} + +template +auto RI_2D_Comm::split_m2D_ktoR_k(const UnitCell& ucell, + const K_Vectors& kv, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin, + const bool spgsym) +-> std::vector>>> +{ + ModuleBase::TITLE("RI_2D_Comm","split_m2D_ktoR_k"); + ModuleBase::timer::tick("RI_2D_Comm", "split_m2D_ktoR_k"); + + const TC period = RI_Util::get_Born_vonKarmen_period(kv); + const std::map nspin_k = {{1,1}, {2,2}, {4,1}}; + const double SPIN_multiple = std::map{ {1,0.5}, {2,1}, {4,1} }.at(nspin); // why? + + std::vector>>> mRs_a2D(nspin); + #ifdef _OPENMP + #pragma omp parallel + #endif + { + std::vector>>> mRs_a2D_thread(nspin); + for (int is_k = 0; is_k < nspin_k.at(nspin); ++is_k) + { + const std::vector ik_list = RI_2D_Comm::get_ik_list(kv, is_k); + const auto cells = RI_Util::get_Born_von_Karmen_cells(period); + #pragma omp for schedule(dynamic) + for (size_t icell = 0; icell < cells.size(); ++icell) + { + const TC& cell = cells[icell]; + RI::Tensor mR_2D; + int ik_full = 0; + for (const int ik : ik_list) + { + using Tdata_m = typename Tmatrix::value_type; + auto set_mR_2D = [&mR_2D](RI::Tensor&& mk_frac) + { + if (mR_2D.empty()) + { mR_2D = RI::Global_Func::convert(mk_frac); } + else + { mR_2D = mR_2D + RI::Global_Func::convert(mk_frac); } + }; + if (!spgsym) { - RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik_full + is_k * kv.get_nkstot_full()], pv.get_col_size(), pv.get_row_size()); + RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik], pv.get_col_size(), pv.get_row_size()); const Tdata_m frac = SPIN_multiple * RI::Global_Func::convert(std::exp( - -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * ((isym_kvd.second * ucell.G) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); - set_mR_2D(mk_2D * frac); - ++ik_full; + -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); + if (static_cast(std::round(SPIN_multiple * kv.wk[ik] * kv.get_nkstot_full())) == 2) + { set_mR_2D(mk_2D * (frac * 0.5) + tensor_conj(mk_2D * (frac * 0.5))); } + else + { set_mR_2D(mk_2D * frac); } } - } - } - for(int iwt0_2D=0; iwt0_2D!=mR_2D.shape[0]; ++iwt0_2D) - { - const int iwt0 =ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) - ? pv.local2global_col(iwt0_2D) - : pv.local2global_row(iwt0_2D); - int iat0, iw0_b, is0_b; - std::tie(iat0,iw0_b,is0_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt0); - const int it0 = ucell.iat2it[iat0]; - for(int iwt1_2D=0; iwt1_2D!=mR_2D.shape[1]; ++iwt1_2D) - { - const int iwt1 =ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) - ? pv.local2global_row(iwt1_2D) - : pv.local2global_col(iwt1_2D); - int iat1, iw1_b, is1_b; - std::tie(iat1,iw1_b,is1_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt1); - const int it1 = ucell.iat2it[iat1]; - - const int is_b = RI_2D_Comm::get_is_block(is_k, is0_b, is1_b); - RI::Tensor &mR_a2D = mRs_a2D[is_b][iat0][{iat1,cell}]; - if (mR_a2D.empty()) { - mR_a2D = RI::Tensor( - {static_cast(ucell.atoms[it0].nw), - static_cast( - ucell.atoms[it1].nw)}); + else + { // traverse kstar, ik means ik_ibz + for (auto& isym_kvd : kv.kstars[ik % ik_list.size()]) + { + RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik_full + is_k * kv.get_nkstot_full()], pv.get_col_size(), pv.get_row_size()); + const Tdata_m frac = SPIN_multiple + * RI::Global_Func::convert(std::exp( + -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * ((isym_kvd.second * ucell.G) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); + set_mR_2D(mk_2D * frac); + ++ik_full; + } + } + } // end for ik + for(int iwt0_2D=0; iwt0_2D!=mR_2D.shape[0]; ++iwt0_2D) + { + const int iwt0 =ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) + ? pv.local2global_col(iwt0_2D) + : pv.local2global_row(iwt0_2D); + int iat0, iw0_b, is0_b; + std::tie(iat0,iw0_b,is0_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt0); + const int it0 = ucell.iat2it[iat0]; + for(int iwt1_2D=0; iwt1_2D!=mR_2D.shape[1]; ++iwt1_2D) + { + const int iwt1 =ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) + ? pv.local2global_row(iwt1_2D) + : pv.local2global_col(iwt1_2D); + int iat1, iw1_b, is1_b; + std::tie(iat1,iw1_b,is1_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt1); + const int it1 = ucell.iat2it[iat1]; + + const int is_b = RI_2D_Comm::get_is_block(is_k, is0_b, is1_b); + RI::Tensor& mR_a2D = mRs_a2D_thread[is_b][iat0][{iat1, cell}]; + if (mR_a2D.empty()) + { + mR_a2D = RI::Tensor( + {static_cast(ucell.atoms[it0].nw), + static_cast(ucell.atoms[it1].nw)}); + } + mR_a2D(iw0_b, iw1_b) = mR_2D(iwt0_2D, iwt1_2D); + } // for iwt1_2D + } // end for iwt0_2D + } // end for icell + } // end for is_k + + #ifdef _OPENMP + #pragma omp critical + #endif + { + for(int is=0; is Date: Wed, 17 Jun 2026 12:42:16 +0800 Subject: [PATCH 2/2] Refactor: simplify Mixing DM --- source/Makefile.Objects | 1 - source/module_ri/CMakeLists.txt | 1 - source/module_ri/Exx_LRI_interface.h | 3 +- source/module_ri/Exx_LRI_interface.hpp | 84 ++++++------ source/module_ri/Mix_DMk_2D.cpp | 110 ++++++++------- source/module_ri/Mix_DMk_2D.h | 65 +++++---- source/module_ri/Mix_Matrix.cpp | 164 ----------------------- source/module_ri/Mix_Matrix.h | 53 -------- source/module_ri/RPA_LRI.hpp | 18 ++- source/module_ri/test/CMakeLists.txt | 4 +- source/module_ri/test/dm_mixing_test.cpp | 163 ++++++---------------- 11 files changed, 189 insertions(+), 477 deletions(-) delete mode 100644 source/module_ri/Mix_Matrix.cpp delete mode 100644 source/module_ri/Mix_Matrix.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index cc023f54d81..0d409aef154 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -601,7 +601,6 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\ Matrix_Orbs22.o\ RI_2D_Comm.o\ Mix_DMk_2D.o\ - Mix_Matrix.o\ symmetry_rotation.o\ symmetry_irreducible_sector.o\ diff --git a/source/module_ri/CMakeLists.txt b/source/module_ri/CMakeLists.txt index a201cab1a31..1f66fe1e6ee 100644 --- a/source/module_ri/CMakeLists.txt +++ b/source/module_ri/CMakeLists.txt @@ -8,7 +8,6 @@ if (ENABLE_LIBRI) Matrix_Orbs22.cpp RI_2D_Comm.cpp Mix_DMk_2D.cpp - Mix_Matrix.cpp ) if(ENABLE_LCAO) diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index 77aaa6c47b4..cc04213f965 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -120,7 +120,8 @@ class Exx_LRI_Interface std::shared_ptr> exx_ptr; private: - Mix_DMk_2D mix_DMk_2D; + + Mix_DMk_2D mix_DMk_2D; bool exx_spacegroup_symmetry = false; ModuleSymmetry::Symmetry_rotation symrot_; diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 0605c7e687d..1858622c0e0 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -62,9 +62,7 @@ void Exx_LRI_Interface::cal_exx_ions(const UnitCell& ucell, const bool ModuleBase::TITLE("Exx_LRI_Interface","cal_exx_ions"); if(!this->flag_finish.init) { throw std::runtime_error("Exx init unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } - this->exx_ptr->cal_exx_ions(ucell, write_cv); - this->flag_finish.ions = true; } @@ -79,7 +77,6 @@ void Exx_LRI_Interface::cal_exx_elec(const std::vectorexx_ptr->cal_exx_elec(Ds, ucell, pv, p_symrot); - this->flag_finish.elec = true; } @@ -93,7 +90,6 @@ void Exx_LRI_Interface::cal_exx_force(const int& nat) { throw std::runtime_error("Exx Hamiltonian unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } this->exx_ptr->cal_exx_force(nat); - this->flag_finish.force = true; } @@ -107,12 +103,14 @@ void Exx_LRI_Interface::cal_exx_stress(const double& omega, const doub { throw std::runtime_error("Exx Hamiltonian unfinished when "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } this->exx_ptr->cal_exx_stress(omega, lat0); - this->flag_finish.stress = true; } template -void Exx_LRI_Interface::exx_before_all_runners(const K_Vectors& kv, const UnitCell& ucell, const Parallel_2D& pv) +void Exx_LRI_Interface::exx_before_all_runners( + const K_Vectors& kv, + const UnitCell& ucell, + const Parallel_2D& pv) { ModuleBase::TITLE("Exx_LRI_Interface","exx_before_all_runners"); // initialize the rotation matrix in AO representation @@ -166,14 +164,15 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, if(GlobalC::exx_info.info_global.cal_exx) { if (this->exx_spacegroup_symmetry) - {this->mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1), PARAM.globalv.gamma_only_local);} + { this->mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1)); } else - {this->mix_DMk_2D.set_nks(kv.get_nks(), PARAM.globalv.gamma_only_local);} + { this->mix_DMk_2D.set_nks(kv.get_nks()); } - if(GlobalC::exx_info.info_global.separate_loop) - { this->mix_DMk_2D.set_mixing(nullptr); } + if (GlobalC::exx_info.info_global.separate_loop) + { this->mix_DMk_2D.set_mixing_plain(GlobalC::exx_info.info_global.mixing_beta_for_loop1); } else { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); } + // for exx two_level scf this->two_level_step = 0; } @@ -190,40 +189,39 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, ModuleBase::TITLE("Exx_LRI_Interface","exx_eachiterinit"); if (GlobalC::exx_info.info_global.cal_exx) { - if (!GlobalC::exx_info.info_global.separate_loop && (this->two_level_step || istep > 0 || PARAM.inp.init_wfc == "file") // non separate loop case - || (GlobalC::exx_info.info_global.separate_loop && PARAM.inp.init_wfc == "file" && this->two_level_step == 0 && iter == 1)) // the first iter in separate loop case + if (!GlobalC::exx_info.info_global.separate_loop + && (this->two_level_step + || istep > 0 + || PARAM.inp.init_wfc == "file") // non separate loop case + || (GlobalC::exx_info.info_global.separate_loop + && PARAM.inp.init_wfc == "file" + && this->two_level_step == 0 + && iter == 1) + ) // the first iter in separate loop case { const bool flag_restart = (iter == 1) ? true : false; auto cal = [this, &ucell,&kv, &flag_restart](const elecstate::DensityMatrix& dm_in) { if (this->exx_spacegroup_symmetry) - { this->mix_DMk_2D.mix(symrot_.restore_dm(kv,dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); } + { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); } else { this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); } - const std::vector>>> - Ds = PARAM.globalv.gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR( - ucell, - *this->exx_ptr->p_kv, - this->mix_DMk_2D.get_DMk_gamma_out(), - *dm_in.get_paraV_pointer(), - PARAM.inp.nspin) - : RI_2D_Comm::split_m2D_ktoR( - ucell, - *this->exx_ptr->p_kv, - this->mix_DMk_2D.get_DMk_k_out(), - *dm_in.get_paraV_pointer(), - PARAM.inp.nspin, - this->exx_spacegroup_symmetry); - - if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + const std::vector>>> Ds = + RI_2D_Comm::split_m2D_ktoR( + ucell, + *this->exx_ptr->p_kv, + this->mix_DMk_2D.get_DMk_out(), + *dm_in.get_paraV_pointer(), + PARAM.inp.nspin, + this->exx_spacegroup_symmetry); + if(this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); } else { this->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer()); } }; if(istep > 0 && flag_restart) - { cal(*dm_last_step); } + { cal(*this->dm_last_step); } else { cal(dm); } } @@ -396,21 +394,23 @@ bool Exx_LRI_Interface::exx_after_converge( // if init_wfc == "file", DM is calculated in the 1st iter of the 1st two-level step, so we mix it here const bool flag_restart = (this->two_level_step == 0 && PARAM.inp.init_wfc != "file") ? true : false; - if (this->exx_spacegroup_symmetry) - {this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), flag_restart);} + if(this->exx_spacegroup_symmetry) + { this->mix_DMk_2D.mix(symrot_.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), flag_restart); } else - {this->mix_DMk_2D.mix(dm.get_DMK_vector(), flag_restart);} - - // GlobalC::exx_lcao.cal_exx_elec(p_esolver->LOC, p_esolver->LOWF.wfc_k_grid); - const std::vector>, RI::Tensor>>> - Ds = std::is_same::value //gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), nspin) - : RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), nspin, this->exx_spacegroup_symmetry); - - if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + { this->mix_DMk_2D.mix(dm.get_DMK_vector(), flag_restart); } + const std::vector>>> Ds = + RI_2D_Comm::split_m2D_ktoR( + ucell, + *this->exx_ptr->p_kv, + this->mix_DMk_2D.get_DMk_out(), + *dm.get_paraV_pointer(), + nspin, + this->exx_spacegroup_symmetry); + if(this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer(), &this->symrot_); } else { this->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer()); } // restore DM but not Hexx + iter = 0; this->two_level_step++; diff --git a/source/module_ri/Mix_DMk_2D.cpp b/source/module_ri/Mix_DMk_2D.cpp index 8abd540a146..bfc610c855b 100644 --- a/source/module_ri/Mix_DMk_2D.cpp +++ b/source/module_ri/Mix_DMk_2D.cpp @@ -4,69 +4,85 @@ //======================= #include "Mix_DMk_2D.h" +#include "module_base/module_mixing/plain_mixing.h" #include "module_base/tool_title.h" -Mix_DMk_2D &Mix_DMk_2D::set_nks(const int nks, const bool gamma_only_in) +#include + +template +Mix_DMk_2D::~Mix_DMk_2D() +{ + if(this->flag_del_mixing) + delete this->mixing; +} + +template +void Mix_DMk_2D::set_nks(const int nks) { - ModuleBase::TITLE("Mix_DMk_2D", "set_nks"); - this->gamma_only = gamma_only_in; - if (this->gamma_only) - this->mix_DMk_gamma.resize(nks); - else - this->mix_DMk_k.resize(nks); - return *this; + this->mix_DMk.clear(); + this->mix_DMk.resize(nks); } -Mix_DMk_2D &Mix_DMk_2D::set_mixing(Base_Mixing::Mixing* mixing_in) +template +void Mix_DMk_2D::set_mixing(Base_Mixing::Mixing* mixing_in) { - ModuleBase::TITLE("Mix_DMk_2D","set_mixing"); - if(this->gamma_only) - for (Mix_Matrix>& mix_one : this->mix_DMk_gamma) - mix_one.init(mixing_in); - else - for (Mix_Matrix>>& mix_one : this->mix_DMk_k) - mix_one.init(mixing_in); - return *this; + if(this->flag_del_mixing) + delete this->mixing; + this->mixing = mixing_in; + this->flag_del_mixing = false; } -Mix_DMk_2D &Mix_DMk_2D::set_mixing_beta(const double mixing_beta) +template +void Mix_DMk_2D::set_mixing_plain(const double& mixing_beta) { - ModuleBase::TITLE("Mix_DMk_2D","set_mixing_beta"); - if(this->gamma_only) - for (Mix_Matrix>& mix_one : this->mix_DMk_gamma) - mix_one.mixing_beta = mixing_beta; - else - for (Mix_Matrix>>& mix_one : this->mix_DMk_k) - mix_one.mixing_beta = mixing_beta; - return *this; + if(this->flag_del_mixing) + delete this->mixing; + this->mixing = new Base_Mixing::Plain_Mixing(mixing_beta); + this->flag_del_mixing = true; } -void Mix_DMk_2D::mix(const std::vector>& dm, const bool flag_restart) +template +void Mix_DMk_2D::mix(const std::vector>& dm, const bool flag_restart) { - ModuleBase::TITLE("Mix_DMk_2D","mix"); - assert(this->mix_DMk_gamma.size() == dm.size()); - for(int ik=0; ikmix_DMk_gamma[ik].mix(dm[ik], flag_restart); + ModuleBase::TITLE("Mix_DMk_2D", "mix"); + if (flag_restart) + { this->restart_all(dm); } + else + { this->mix_all(dm); } } -void Mix_DMk_2D::mix(const std::vector>>& dm, const bool flag_restart) + +template +std::vector*> Mix_DMk_2D::get_DMk_out() const { - ModuleBase::TITLE("Mix_DMk_2D","mix"); - assert(this->mix_DMk_k.size() == dm.size()); - for(int ik=0; ikmix_DMk_k[ik].mix(dm[ik], flag_restart); + std::vector*> DMk_out(this->mix_DMk.size()); + for (int ik = 0; ik < this->mix_DMk.size(); ++ik) + { DMk_out[ik] = &this->mix_DMk[ik].data_out; } + return DMk_out; } -std::vector*> Mix_DMk_2D::get_DMk_gamma_out() const +template +void Mix_DMk_2D::restart_all(const std::vector>& data_in) { - std::vector*> DMk_out(this->mix_DMk_gamma.size()); - for(int ik=0; ikmix_DMk_gamma.size(); ++ik) - DMk_out[ik] = &this->mix_DMk_gamma[ik].get_data_out(); - return DMk_out; + assert(this->mix_DMk.size() == data_in.size()); + assert(this->mixing != nullptr); + for (int ik = 0; ik < data_in.size(); ++ik) + { + this->mix_DMk[ik].data_out = data_in[ik]; + this->mixing->init_mixing_data(this->mix_DMk[ik].mixing_data, data_in[ik].size(), sizeof(Tdata)); + } } -std::vector>*> Mix_DMk_2D::get_DMk_k_out() const + +template +void Mix_DMk_2D::mix_all(const std::vector>& data_in) { - std::vector>*> DMk_out(this->mix_DMk_k.size()); - for(int ik=0; ikmix_DMk_k.size(); ++ik) - DMk_out[ik] = &this->mix_DMk_k[ik].get_data_out(); - return DMk_out; -} \ No newline at end of file + assert(this->mix_DMk.size() == data_in.size()); + assert(this->mixing != nullptr); + for (int ik = 0; ik < data_in.size(); ++ik) + { + this->mixing->push_data(this->mix_DMk[ik].mixing_data, this->mix_DMk[ik].data_out.data(), data_in[ik].data(), nullptr, false); + this->mixing->mix_data(this->mix_DMk[ik].mixing_data, this->mix_DMk[ik].data_out.data()); + } +} + +template class Mix_DMk_2D; +template class Mix_DMk_2D>; diff --git a/source/module_ri/Mix_DMk_2D.h b/source/module_ri/Mix_DMk_2D.h index 9d640fa7f12..f0ff9bcfa91 100644 --- a/source/module_ri/Mix_DMk_2D.h +++ b/source/module_ri/Mix_DMk_2D.h @@ -6,67 +6,62 @@ #ifndef MIX_DMK_2D_H #define MIX_DMK_2D_H -#include "Mix_Matrix.h" #include "module_base/module_mixing/mixing.h" -#include "module_base/matrix.h" -#include "module_base/complexmatrix.h" +#include #include +template class Mix_DMk_2D { public: + ~Mix_DMk_2D(); + /** - * @brief Sets the number of k-points and gamma_only flag. + * @brief Sets the number of k-points. * @param nks Number of k-points. - * @param gamma_only_in Flag indicating if only gamma point is mixed. - * @return Reference to the current object. */ - Mix_DMk_2D &set_nks(const int nks, const bool gamma_only_in); + void set_nks(const int nks); /** * @brief Sets the mixing mode. * @param Mixing Mixing pointer. - * @return Reference to the current object. */ - Mix_DMk_2D &set_mixing(Base_Mixing::Mixing* mixing_in); + void set_mixing(Base_Mixing::Mixing* mixing_in); /** - * @brief Sets the mixing beta. - * @param mixing_beta Mixing beta. - * @return Reference to the current object. - */ - Mix_DMk_2D &set_mixing_beta(const double mixing_beta); - - /** - * @brief Mixes the double density matrix. - * @param dm Double Density matrix. - * @param flag_restart Flag indicating whether restart mixing. + * @brief Sets Base_Mixing::Plain_Mixing. + * @param mixing_beta mixing beta for plain mixing. */ - void mix(const std::vector>& dm, const bool flag_restart); + void set_mixing_plain(const double& mixing_beta); /** - * @brief Mixes the complex density matrix. - * @param dm Complex density matrix. + * @brief Mixes the density matrix. + * @param dm Density matrix. * @param flag_restart Flag indicating whether restart mixing. */ - void mix(const std::vector>>& dm, const bool flag_restart); + void mix(const std::vector>& dm, const bool flag_restart); /** - * @brief Returns the double density matrix. - * @return Double density matrices for each k-points. + * @brief Returns the density matrix. + * @return Density matrices for each k-points. */ - std::vector*> get_DMk_gamma_out() const; - /** - * @brief Returns the complex density matrix. - * @return Complex density matrices for each k-points. - */ - std::vector>*> get_DMk_k_out() const; + std::vector*> get_DMk_out() const; private: - std::vector>> mix_DMk_gamma; - std::vector>>> mix_DMk_k; - bool gamma_only; + struct DMk_Mix_Data + { + std::vector data_out; + Base_Mixing::Mixing_Data mixing_data; + }; + + void restart_all(const std::vector>& data_in); + + void mix_all(const std::vector>& data_in); + + std::vector mix_DMk; + Base_Mixing::Mixing* mixing = nullptr; + bool flag_del_mixing = false; }; -#endif \ No newline at end of file +#endif diff --git a/source/module_ri/Mix_Matrix.cpp b/source/module_ri/Mix_Matrix.cpp deleted file mode 100644 index 21d5b5e3c20..00000000000 --- a/source/module_ri/Mix_Matrix.cpp +++ /dev/null @@ -1,164 +0,0 @@ -//======================= -// AUTHOR : Peize Lin -// DATE : 2023-05-09 -//======================= - -#ifndef MIX_DATA_HPP -#define MIX_DATA_HPP - -#include "Mix_Matrix.h" -#include "module_base/matrix.h" -#include "module_base/complexmatrix.h" -#include "module_base/tool_title.h" - -template<> -void Mix_Matrix::mix(const ModuleBase::matrix& data_in, const bool flag_restart) -{ - if(separate_loop) - { - this->mixing = new Base_Mixing::Plain_Mixing(this->mixing_beta); - } - - if(flag_restart) - { - this->data_out = data_in; - this->mixing->init_mixing_data(this->matrix_data, data_in.nc*data_in.nr, sizeof(*data_in.c)); - } - else - { - this->mixing->push_data(this->matrix_data, data_out.c, data_in.c, nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.c); - } - - if(separate_loop) - { - delete this->mixing; - this->mixing = nullptr; - } -} -template<> -void Mix_Matrix::mix(const ModuleBase::ComplexMatrix& data_in, const bool flag_restart) -{ - if (separate_loop) - { - this->mixing = new Base_Mixing::Plain_Mixing(this->mixing_beta); - } - - if (flag_restart) - { - this->data_out = data_in; - this->mixing->init_mixing_data(this->matrix_data, data_in.nc * data_in.nr, sizeof(*data_in.c)); - } - else - { - this->mixing->push_data(this->matrix_data, data_out.c, data_in.c, nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.c); - } - - if (separate_loop) - { - delete this->mixing; - this->mixing = nullptr; - } -} - -template<> -void Mix_Matrix>::mix(const std::vector& data_in, const bool flag_restart) -{ - if (separate_loop) - { - this->mixing = new Base_Mixing::Plain_Mixing(this->mixing_beta); - } - - if (flag_restart) - { - this->data_out = data_in; - this->mixing->init_mixing_data(this->matrix_data, data_in.size(), sizeof(*data_in.data())); - } - else - { - this->mixing->push_data(this->matrix_data, data_out.data(), data_in.data(), nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.data()); - } - - if (separate_loop) - { - delete this->mixing; - this->mixing = nullptr; - } -} -template<> -void Mix_Matrix>>::mix(const std::vector>& data_in, const bool flag_restart) -{ - if (separate_loop) - { - this->mixing = new Base_Mixing::Plain_Mixing(this->mixing_beta); - } - - if (flag_restart) - { - this->data_out = data_in; - this->mixing->init_mixing_data(this->matrix_data, data_in.size(), sizeof(*data_in.data())); - } - else - { - this->mixing->push_data(this->matrix_data, data_out.data(), data_in.data(), nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.data()); - } - - if (separate_loop) - { - delete this->mixing; - this->mixing = nullptr; - } -} - -/* -// mix for ct::Tensor if useful in the future. comment out for UT-coverage -template<> -void Mix_Matrix::mix(const ct::Tensor& data_in, const bool flag_restart) -{ - ModuleBase::TITLE("Mix_Matrix", "mix"); - if (separate_loop) - { - this->mixing = new Base_Mixing::Plain_Mixing(this->mixing_beta); - } - - if (flag_restart) - { - this->data_out = data_in; - this->mixing->init_mixing_data(this->matrix_data, data_in.NumElements(), sizeof(data_in.SizeOfType(data_in.data_type()))); - } - else - { - switch (data_in.data_type()) - { - case ct::DataType::DT_DOUBLE: - this->mixing->push_data(this->matrix_data, data_out.data(), data_in.data(), nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.data()); - break; - case ct::DataType::DT_COMPLEX_DOUBLE: - this->mixing->push_data(this->matrix_data, data_out.data>(), data_in.data>(), nullptr, false); - this->mixing->mix_data(this->matrix_data, data_out.data>()); - break; - default: - throw std::invalid_argument("Mix_Matrix: data type not supported"); - break; - } - } - - if (separate_loop) - { - delete this->mixing; - this->mixing = nullptr; - } -} -*/ - -template class Mix_Matrix; -template class Mix_Matrix; -template class Mix_Matrix>; -template class Mix_Matrix>>; -// template class Mix_Matrix; - -#endif \ No newline at end of file diff --git a/source/module_ri/Mix_Matrix.h b/source/module_ri/Mix_Matrix.h deleted file mode 100644 index 820a6c47f2f..00000000000 --- a/source/module_ri/Mix_Matrix.h +++ /dev/null @@ -1,53 +0,0 @@ -//======================= -// AUTHOR : Peize Lin -// DATE : 2023-05-09 -//======================= - -#ifndef MIX_MATRIX_H -#define MIX_MATRIX_H - -#include - -#include "module_base/module_mixing/mixing.h" -#include "module_base/module_mixing/plain_mixing.h" -template -class Mix_Matrix -{ - public: - /** - * @brief init Mix_Matrix with Mixing pointer - * - * @param mixing_in Mixing pointer - */ - void init(Base_Mixing::Mixing* mixing_in) - { - this->mixing = mixing_in; - if (this->mixing == nullptr) - this->separate_loop = true; - } - /** - * @brief Get the data out object - * - * @return const Tdata& - */ - const Tdata& get_data_out() const - { - return this->data_out; - } - - /** - * @brief Mixes the input data according to the set mixing mode. - * @param data_in Input data to be mixed. - * @param flag_restart Flag indicating whether restart mixing. - */ - void mix(const Tdata& data_in, const bool flag_restart); - - double mixing_beta = 1.0; - - private: - Tdata data_out; - Base_Mixing::Mixing* mixing = nullptr; - Base_Mixing::Mixing_Data matrix_data; - bool separate_loop = false; -}; -#endif \ No newline at end of file diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 79c851fa3a7..7967b67baa4 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -79,14 +79,14 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix const K_Vectors& kv, const LCAO_Orbitals& orb) { - Mix_DMk_2D mix_DMk_2D; + Mix_DMk_2D mix_DMk_2D; bool exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1); if (exx_spacegroup_symmetry) - {mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1), PARAM.globalv.gamma_only_local);} + {mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1));} else - {mix_DMk_2D.set_nks(kv.get_nks(), PARAM.globalv.gamma_only_local);} + {mix_DMk_2D.set_nks(kv.get_nks());} - mix_DMk_2D.set_mixing(nullptr); + mix_DMk_2D.set_mixing_plain(1.0); ModuleSymmetry::Symmetry_rotation symrot; if (exx_spacegroup_symmetry) { @@ -99,9 +99,13 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix else { mix_DMk_2D.mix(dm.get_DMK_vector(), true); } const std::vector>>> - Ds = PARAM.globalv.gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR(ucell,kv, mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin) - : RI_2D_Comm::split_m2D_ktoR(ucell,kv, mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin, exx_spacegroup_symmetry); + Ds = RI_2D_Comm::split_m2D_ktoR( + ucell, + kv, + mix_DMk_2D.get_DMk_out(), + *dm.get_paraV_pointer(), + PARAM.inp.nspin, + exx_spacegroup_symmetry); // set parameters for bare Coulomb potential GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; diff --git a/source/module_ri/test/CMakeLists.txt b/source/module_ri/test/CMakeLists.txt index 9573a3827eb..701e585b439 100644 --- a/source/module_ri/test/CMakeLists.txt +++ b/source/module_ri/test/CMakeLists.txt @@ -4,11 +4,11 @@ remove_definitions(-D__ROCM) AddTest( TARGET dm_mixing_test LIBS parameter base ${math_libs} device - SOURCES dm_mixing_test.cpp ../Mix_DMk_2D.cpp ../Mix_Matrix.cpp + SOURCES dm_mixing_test.cpp ../Mix_DMk_2D.cpp ) AddTest( TARGET ri_cv_io_test LIBS base ${math_libs} device parameter SOURCES ri_cv_io_test.cpp ) -install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) \ No newline at end of file +install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/source/module_ri/test/dm_mixing_test.cpp b/source/module_ri/test/dm_mixing_test.cpp index 9831f5101dd..dfe539b9d62 100644 --- a/source/module_ri/test/dm_mixing_test.cpp +++ b/source/module_ri/test/dm_mixing_test.cpp @@ -2,7 +2,6 @@ #include "gtest/gtest.h" #include "module_base/module_mixing/broyden_mixing.h" #include "module_ri/Mix_DMk_2D.h" -#include "module_ri/Mix_Matrix.h" /************************************************ * unit test of charge_mixing.cpp & Mix_DMk_2D.cpp @@ -10,8 +9,6 @@ /** * - Tested Functions: - * - Mix_Matrix::mix: - * mix the matrix data according to the set mixing mode. * - Mix_DMk_2D::mix: * mix the density matrix data according to the set mixing mode. * @@ -23,30 +20,6 @@ class DM_Mixing_Test : public ::testing::Test DM_Mixing_Test() { mixing = new Base_Mixing::Broyden_Mixing(ndim, mixing_beta); - mix_data = std::vector(2); - mix_complexdata = std::vector(3); - mix_data[0].create(nr, nc); - mix_data[1].create(nr, nc); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - mix_data[0](i, j) = i * nc + j; - mix_data[1](i, j) = i * nc + j + 0.2; - } - } - mix_complexdata[0].create(nr, nc); - mix_complexdata[1].create(nr, nc); - mix_complexdata[2].create(nr, nc); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - mix_complexdata[0](i, j) = std::complex{double(i), double(j)}; - mix_complexdata[1](i, j) = std::complex{double(i), double(j) + 0.2}; - mix_complexdata[2](i, j) = std::complex{double(i) + 0.8, double(j)}; - } - } mix_data_vector = std::vector>(2); mix_complexdata_vector = std::vector>>(3); mix_data_vector[0].resize(nr * nc); @@ -83,99 +56,16 @@ class DM_Mixing_Test : public ::testing::Test const double mixing_beta = 0.3; protected: - std::vector mix_data; - std::vector mix_complexdata; std::vector> mix_data_vector; std::vector>> mix_complexdata_vector; }; -TEST_F(DM_Mixing_Test, Mix_Matrix) -{ - // Separete loop - // matrix - Mix_Matrix mix_matrix; - mix_matrix.init(nullptr); - mix_matrix.mixing_beta = mixing_beta; - for (int istep = 0; istep < 2; ++istep) - { - mix_matrix.mix(mix_data[istep], (istep == 0)); - } - - ModuleBase::matrix data_out = mix_matrix.get_data_out(); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - EXPECT_DOUBLE_EQ(data_out(i, j), (1 - mixing_beta) * mix_data[0](i, j) + mixing_beta * mix_data[1](i, j)); - } - } - // vector - Mix_Matrix> mix_vector; - mix_vector.init(nullptr); - mix_vector.mixing_beta = mixing_beta; - for (int istep = 0; istep < 2; ++istep) - { - mix_vector.mix(mix_data_vector[istep], (istep == 0)); - } - - std::vector data_out_vector = mix_vector.get_data_out(); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - EXPECT_DOUBLE_EQ(data_out_vector[i * nc + j], (1 - mixing_beta) * mix_data_vector[0][i * nc + j] + mixing_beta * mix_data_vector[1][i * nc + j]); - } - } - //vector> - Mix_Matrix>> mix_vector_complex; - mix_vector_complex.init(nullptr); - mix_vector_complex.mixing_beta = mixing_beta; - for (int istep = 0; istep < 2; ++istep) - { - mix_vector_complex.mix(mix_complexdata_vector[istep], (istep == 0)); - } - - std::vector> data_out_vector_complex = mix_vector_complex.get_data_out(); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - EXPECT_DOUBLE_EQ(data_out_vector_complex[i * nc + j].real(), ((1 - mixing_beta) * mix_complexdata_vector[0][i * nc + j] + mixing_beta * mix_complexdata_vector[1][i * nc + j]).real()); - EXPECT_DOUBLE_EQ(data_out_vector_complex[i * nc + j].imag(), ((1 - mixing_beta) * mix_complexdata_vector[0][i * nc + j] + mixing_beta * mix_complexdata_vector[1][i * nc + j]).imag()); - } - } - - // Broyden mix - Mix_Matrix mix_complexmatrix; - mix_complexmatrix.init(mixing); - mixing->coef = { 1.1, -0.1 }; - for (int istep = 0; istep < 3; ++istep) - { - mix_complexmatrix.mix(mix_complexdata[istep], (istep == 0)); - } - ModuleBase::ComplexMatrix com_data_out = mix_complexmatrix.get_data_out(); - for (int i = 0; i < nr; ++i) - { - for (int j = 0; j < nc; ++j) - { - std::complex first_step_result - = (1 - mixing_beta) * mix_complexdata[0](i, j) + mixing_beta * mix_complexdata[1](i, j); - std::complex second_step_result - = (1 - mixing_beta) * first_step_result + mixing_beta * mix_complexdata[2](i, j); - std::complex ref = second_step_result * mixing->coef[1] + first_step_result * mixing->coef[0]; - EXPECT_DOUBLE_EQ(com_data_out(i, j).real(), ref.real()); - EXPECT_DOUBLE_EQ(com_data_out(i, j).imag(), ref.imag()); - } - } -} - TEST_F(DM_Mixing_Test, Mix_DMk_2D) { //Gamma only - Mix_DMk_2D mix_dmk_gamma; - mix_dmk_gamma.set_nks(1, true); - mix_dmk_gamma.set_mixing(nullptr); - mix_dmk_gamma.set_mixing_beta(mixing_beta); + Mix_DMk_2D mix_dmk_gamma; + mix_dmk_gamma.set_nks(1); + mix_dmk_gamma.set_mixing_plain(1.0); std::vector>> dm_gamma(2); dm_gamma[0] = std::vector>(1); dm_gamma[0][0] = mix_data_vector[0]; @@ -185,20 +75,19 @@ TEST_F(DM_Mixing_Test, Mix_DMk_2D) { mix_dmk_gamma.mix(dm_gamma[istep], (istep == 0)); } - std::vector*> dm_gamma_out = mix_dmk_gamma.get_DMk_gamma_out(); + std::vector*> dm_gamma_out = mix_dmk_gamma.get_DMk_out(); for (int i = 0; i < nr; ++i) { for (int j = 0; j < nc; ++j) { - EXPECT_DOUBLE_EQ(dm_gamma_out[0][0][i * nc + j], (1 - mixing_beta) * mix_data_vector[0][i * nc + j] + mixing_beta * mix_data_vector[1][i * nc + j]); + EXPECT_DOUBLE_EQ(dm_gamma_out[0][0][i * nc + j], mix_data_vector[1][i * nc + j]); } } // not Gamma only - Mix_DMk_2D mix_dmk; - mix_dmk.set_nks(1, false); - mix_dmk.set_mixing(nullptr); - mix_dmk.set_mixing_beta(mixing_beta); + Mix_DMk_2D> mix_dmk; + mix_dmk.set_nks(1); + mix_dmk.set_mixing_plain(1.0); std::vector>>> dm(2); dm[0] = std::vector>>(1); dm[0][0] = mix_complexdata_vector[0]; @@ -208,15 +97,41 @@ TEST_F(DM_Mixing_Test, Mix_DMk_2D) { mix_dmk.mix(dm[istep], (istep == 0)); } - std::vector>*> dm_out = mix_dmk.get_DMk_k_out(); + std::vector>*> dm_out = mix_dmk.get_DMk_out(); + for (int i = 0; i < nr; ++i) + { + for (int j = 0; j < nc; ++j) + { + EXPECT_DOUBLE_EQ(dm_out[0][0][i * nc + j].real(), mix_complexdata_vector[1][i * nc + j].real()); + EXPECT_DOUBLE_EQ(dm_out[0][0][i * nc + j].imag(), mix_complexdata_vector[1][i * nc + j].imag()); + } + } + + // Shared Broyden mix + Mix_DMk_2D> mix_dmk_broyden; + mix_dmk_broyden.set_nks(1); + mix_dmk_broyden.set_mixing(mixing); + mixing->coef = { 1.1, -0.1 }; + std::vector>>> dm_broyden(3); + for (int istep = 0; istep < 3; ++istep) + { + dm_broyden[istep] = std::vector>>(1); + dm_broyden[istep][0] = mix_complexdata_vector[istep]; + mix_dmk_broyden.mix(dm_broyden[istep], (istep == 0)); + } + std::vector>*> dm_broyden_out = mix_dmk_broyden.get_DMk_out(); for (int i = 0; i < nr; ++i) { for (int j = 0; j < nc; ++j) { std::complex first_step_result - = (1 - mixing_beta) * mix_complexdata_vector[0][i * nc + j] + mixing_beta * mix_complexdata_vector[1][i * nc + j]; - EXPECT_DOUBLE_EQ(dm_out[0][0][i * nc + j].real(), first_step_result.real()); - EXPECT_DOUBLE_EQ(dm_out[0][0][i * nc + j].imag(), first_step_result.imag()); + = (1 - mixing_beta) * mix_complexdata_vector[0][i * nc + j] + + mixing_beta * mix_complexdata_vector[1][i * nc + j]; + std::complex second_step_result + = (1 - mixing_beta) * first_step_result + mixing_beta * mix_complexdata_vector[2][i * nc + j]; + std::complex ref = second_step_result * mixing->coef[1] + first_step_result * mixing->coef[0]; + EXPECT_DOUBLE_EQ(dm_broyden_out[0][0][i * nc + j].real(), ref.real()); + EXPECT_DOUBLE_EQ(dm_broyden_out[0][0][i * nc + j].imag(), ref.imag()); } } -} \ No newline at end of file +}