diff --git a/source/source_basis/module_ao/ORB_gaunt_table.cpp b/source/source_basis/module_ao/ORB_gaunt_table.cpp index b930c6bcba..d05f594f7e 100644 --- a/source/source_basis/module_ao/ORB_gaunt_table.cpp +++ b/source/source_basis/module_ao/ORB_gaunt_table.cpp @@ -20,6 +20,7 @@ void ORB_gaunt_table::init_Gaunt(const int &lmax) ModuleBase::TITLE("ORB_gaunt_table", "init_Gaunt"); ModuleBase::timer::tick("ORB_gaunt_table", "init_Gaunt"); + this->Lmax_Gaunt_Coefficients = lmax; const int nlm = (lmax * 2 + 1) * (lmax * 2 + 1); this->Gaunt_Coefficients.create(nlm, nlm, nlm); @@ -184,6 +185,7 @@ void ORB_gaunt_table::init_Gaunt_CH(const int& Lmax) ModuleBase::TITLE("ORB_gaunt_table","init_Gaunt_CH"); ModuleBase::timer::tick("ORB_gaunt_table","init_Gaunt_CH"); + this->Lmax_Gaunt_CH = Lmax; int L = 2*Lmax + 1; int Eff_Np = this->EP_EL(L); diff --git a/source/source_basis/module_ao/ORB_gaunt_table.h b/source/source_basis/module_ao/ORB_gaunt_table.h index 1e295b9125..0e9aefdb70 100644 --- a/source/source_basis/module_ao/ORB_gaunt_table.h +++ b/source/source_basis/module_ao/ORB_gaunt_table.h @@ -93,6 +93,9 @@ class ORB_gaunt_table static int Index_M(const int& m); + int get_Lmax_Gaunt_Coefficients() const { return Lmax_Gaunt_Coefficients; } + int get_Lmax_Gaunt_CH() const { return Lmax_Gaunt_CH; } + private: // Index Function @@ -122,5 +125,8 @@ class ORB_gaunt_table //direct integral ModuleBase::matrix Ylm_Gaunt; + + int Lmax_Gaunt_Coefficients = -1; + int Lmax_Gaunt_CH = -1; }; #endif diff --git a/source/source_hamilt/module_xc/exx_info.h b/source/source_hamilt/module_xc/exx_info.h index 9cbd51c18f..5cca433b94 100644 --- a/source/source_hamilt/module_xc/exx_info.h +++ b/source/source_hamilt/module_xc/exx_info.h @@ -69,8 +69,6 @@ struct Exx_Info double kmesh_times = 4; double Cs_inv_thr = -1; - int abfs_Lmax = 0; // tmp - Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global) : coulomb_param(info_global.coulomb_param) { @@ -80,7 +78,7 @@ struct Exx_Info struct Exx_Info_Opt_ABFs { - int abfs_Lmax = 0; // tmp + int abfs_Lmax = 0; double ecut_exx = 60; double tolerence = 1E-12; std::vector files_jles; diff --git a/source/source_lcao/LCAO_init_basis.cpp b/source/source_lcao/LCAO_init_basis.cpp index f8b60b6298..6d198a3a32 100644 --- a/source/source_lcao/LCAO_init_basis.cpp +++ b/source/source_lcao/LCAO_init_basis.cpp @@ -62,11 +62,6 @@ void init_basis_lcao(Parallel_Orbitals& pv, two_center_bundle.build_beta(ucell.ntype, ucell.infoNL.Beta); } - int Lmax = 0; -#ifdef __EXX - Lmax = GlobalC::exx_info.info_ri.abfs_Lmax; -#endif - #ifdef USE_NEW_TWO_CENTER two_center_bundle.tabulate(); #else diff --git a/source/source_lcao/center2_orb.cpp b/source/source_lcao/center2_orb.cpp index 8df630777b..b5ece1d3e3 100644 --- a/source/source_lcao/center2_orb.cpp +++ b/source/source_lcao/center2_orb.cpp @@ -49,8 +49,7 @@ std::pair Center2_Orb::init_Lmax_2_2(const int& lmax_exx) // used in berryphase by jingan std::pair Center2_Orb::init_Lmax_2_3(const int lmax_orb) { - int Lmax = std::max(-1, lmax_orb); - Lmax++; + const int Lmax = std::max(-1, lmax_orb) + 1; const int Lmax_used = 2 * Lmax + 1; assert(Lmax_used >= 1); return {Lmax_used, Lmax}; @@ -172,15 +171,14 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, // double* integrated_func = new double[kmesh]; - int ll = 0; - if (l != 0) - { - ll = l - 1; - } - - const std::vector>& jlm1 = psb->get_jlx()[ll]; - const std::vector>& jl = psb->get_jlx()[l]; - const std::vector>& jlp1 = psb->get_jlx()[l + 1]; + assert(psb->get_jlx().size()>=l+2); + const int lml = (l>0) ? (l-1) : 0; + const std::vector>& jlm1 = psb->get_jlx().at(lml); + const std::vector>& jl = psb->get_jlx().at(l); + const std::vector>& jlp1 = psb->get_jlx().at(l+1); + assert(jlm1.size()>=rmesh); + assert(jl.size()>=rmesh); + assert(jlp1.size()>=rmesh); #ifdef _OPENMP #pragma omp parallel for schedule(static) @@ -189,6 +187,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, { std::vector integrated_func(kmesh); const std::vector& jl_r = jl[ir]; + assert(jl_r.size()>=kmesh); for (int ik = 0; ik < kmesh; ++ik) { integrated_func[ik] = jl_r[ik] * k1_dot_k2[ik]; @@ -202,6 +201,8 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, // Peize Lin accelerate 2017-10-02 const std::vector& jlm1_r = jlm1[ir]; const std::vector& jlp1_r = jlp1[ir]; + assert(jlm1_r.size()>=kmesh); + assert(jlp1_r.size()>=kmesh); const double fac = l / (l + 1.0); if (l == 0) { @@ -305,10 +306,11 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, std::vector integrated_func(kmesh); - const int lm1 = (l > 0 ? l - 1 : 0); - const std::vector>& jlm1 = psb->get_jlx()[lm1]; - const std::vector>& jl = psb->get_jlx()[l]; - const std::vector>& jlp1 = psb->get_jlx()[l + 1]; + assert(psb->get_jlx().size()>=l+2); + const int lm1 = (l>0) ? (l-1) : 0; + const std::vector>& jlm1 = psb->get_jlx().at(lm1); + const std::vector>& jl = psb->get_jlx().at(l); + const std::vector>& jlp1 = psb->get_jlx().at(l+1); for (const size_t& ir: radials) { @@ -320,6 +322,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, } const std::vector& jl_r = jl[ir]; + assert(jl_r.size()>=kmesh); for (int ik = 0; ik < kmesh; ++ik) { integrated_func[ik] = jl_r[ik] * k1_dot_k2[ik]; @@ -329,8 +332,10 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job, ModuleBase::Integral::Simpson_Integral(kmesh, ModuleBase::GlobalFunc::VECTOR_TO_PTR(integrated_func), dk, temp); rs[ir] = temp * ModuleBase::FOUR_PI; - const std::vector& jlm1_r = jlm1[ir]; - const std::vector& jlp1_r = jlp1[ir]; + const std::vector& jlm1_r = jlm1.at(ir); + const std::vector& jlp1_r = jlp1.at(ir); + assert(jlm1_r.size()>=kmesh); + assert(jlp1_r.size()>=kmesh); const double fac = l / (l + 1.0); if (l == 0) { diff --git a/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp b/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp index aaa1ef2f41..f3efdfbd7e 100644 --- a/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp +++ b/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp @@ -172,29 +172,13 @@ RI::Tensor get_column_mean0_matrix(const RI::Tensor& m) const ModuleBase::Element_Basis_Index::IndexLNM index_abfs = ModuleBase::Element_Basis_Index::construct_index(range_abfs); - const int Lmax_bak = GlobalC::exx_info.info_ri.abfs_Lmax; - GlobalC::exx_info.info_ri.abfs_Lmax = std::numeric_limits::min(); - for (std::size_t T = 0; T != abfs.size(); ++T) - { - GlobalC::exx_info.info_ri.abfs_Lmax - = std::max(GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(abfs[T].size()) - 1); -} - - Matrix_Orbs21 m_abfslcaos_lcaos; - ORB_gaunt_table MGT; - const int Lmax = m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, orb.get_Rmax()); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_abfslcaos_lcaos.init_radial(abfs, lcaos, lcaos, MGT); + Matrix_Orbs21 m_abfslcaos_lcaos; + m_abfslcaos_lcaos.init(abfs, lcaos, lcaos, ucell, orb, kmesh_times, orb.get_Rmax()); std::map>> delta_R; for (std::size_t it = 0; it != abfs.size(); ++it) - { - delta_R[it][it] = {0.0}; -} + { delta_R[it][it] = {0.0}; } m_abfslcaos_lcaos.init_radial_table(delta_R); - - GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak; std::vector, RI::Tensor>>> eig(abfs.size()); for (std::size_t T = 0; T != abfs.size(); ++T) diff --git a/source/source_lcao/module_ri/Exx_LRI.h b/source/source_lcao/module_ri/Exx_LRI.h index 0ca4aff9b6..6f94f1a32b 100644 --- a/source/source_lcao/module_ri/Exx_LRI.h +++ b/source/source_lcao/module_ri/Exx_LRI.h @@ -90,7 +90,6 @@ class Exx_LRI const Exx_Info::Exx_Info_RI &info; MPI_Comm mpi_comm; const K_Vectors *p_kv = nullptr; - ORB_gaunt_table MGT; std::vector orb_cutoff_; std::vector>> lcaos; diff --git a/source/source_lcao/module_ri/Exx_LRI.hpp b/source/source_lcao/module_ri/Exx_LRI.hpp index 5cf3a0072a..ee8c62a718 100644 --- a/source/source_lcao/module_ri/Exx_LRI.hpp +++ b/source/source_lcao/module_ri/Exx_LRI.hpp @@ -39,6 +39,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, this->orb_cutoff_ = orb.cutoffs(); this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos); const std::vector>> abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell, orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); @@ -46,21 +47,18 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, { this->abfs = abfs_same_atom;} else { this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); } + Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs); Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running); - for( size_t T=0; T!=this->abfs.size(); ++T ) - { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(this->abfs[T].size())-1 ); } - this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv); - bool init_MGT = true; + std::shared_ptr MGT = std::make_shared(); for(const auto &settings_list : this->coulomb_settings) { this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, settings_list.second.second, this->info.ccp_rmesh_times); this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb, this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp, - this->info.kmesh_times, this->MGT, init_MGT, settings_list.second.first ); - init_MGT = false; // only init once + this->info.kmesh_times, MGT, settings_list.second.first ); } ModuleBase::timer::tick("Exx_LRI", "init"); @@ -73,11 +71,6 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, ModuleBase::TITLE("Exx_LRI","cal_exx_ions"); ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions"); - // init_radial_table_ions( cal_atom_centres_core(atom_pairs_core_origin), atom_pairs_core_origin ); - - // this->m_abfsabfs.init_radial_table(Rradial); - // this->m_abfslcaos_lcaos.init_radial_table(Rradial); - std::vector atoms(ucell.nat); for(int iat=0; iat::exx_before_all_runners( this->symrot_.find_irreducible_sector( ucell.symm, ucell.atoms, ucell.st, RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat); - this->symrot_.set_abfs_Lmax(GlobalC::exx_info.info_ri.abfs_Lmax); + this->symrot_.set_abfs_Lmax(Exx_Abfs::get_Lmax(this->exx_ptr->abfs)); this->symrot_.cal_Ms(kv, ucell, pv); } } diff --git a/source/source_lcao/module_ri/LRI_CV.h b/source/source_lcao/module_ri/LRI_CV.h index 210fea1afa..aaf00e8656 100644 --- a/source/source_lcao/module_ri/LRI_CV.h +++ b/source/source_lcao/module_ri/LRI_CV.h @@ -40,8 +40,7 @@ class LRI_CV const std::vector>> &abfs_in, const std::vector>> &abfs_ccp_in, const double &kmesh_times, - ORB_gaunt_table& MGT, - const bool& init_MGT, + std::shared_ptr MGT, const bool& init_C); inline std::map>> cal_Vs( diff --git a/source/source_lcao/module_ri/LRI_CV.hpp b/source/source_lcao/module_ri/LRI_CV.hpp index 32e01302e7..60670d5088 100644 --- a/source/source_lcao/module_ri/LRI_CV.hpp +++ b/source/source_lcao/module_ri/LRI_CV.hpp @@ -44,8 +44,7 @@ void LRI_CV::set_orbitals( const std::vector>> &abfs_in, const std::vector>> &abfs_ccp_in, const double &kmesh_times, - ORB_gaunt_table& MGT, - const bool& init_MGT, + std::shared_ptr MGT, const bool& init_C) { ModuleBase::TITLE("LRI_CV", "set_orbitals"); @@ -69,24 +68,17 @@ void LRI_CV::set_orbitals( range_abfs = ModuleBase::Element_Basis_Index::construct_range( abfs ); this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs ); - const int Lmax_v = this->m_abfs_abfs.init(2, ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax); - int Lmax_c = std::numeric_limits::min(); + this->m_abfs_abfs.MGT = this->m_abfslcaos_lcaos.MGT = MGT; + this->m_abfs_abfs.init( + this->abfs_ccp, this->abfs, + ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax); if (init_C) - Lmax_c = this->m_abfslcaos_lcaos.init(1, ucell, orb, kmesh_times, lcaos_rmax); - const int Lmax = std::max(Lmax_v, Lmax_c); + this->m_abfslcaos_lcaos.init( + this->abfs_ccp, this->lcaos, this->lcaos, + ucell, orb, kmesh_times, lcaos_rmax); - if (init_MGT) { - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - } - - this->m_abfs_abfs.init_radial(this->abfs_ccp, this->abfs, MGT); this->m_abfs_abfs.init_radial_table(); if (init_C) { - this->m_abfslcaos_lcaos.init_radial(this->abfs_ccp, - this->lcaos, - this->lcaos, - MGT); this->m_abfslcaos_lcaos.init_radial_table(); } diff --git a/source/source_lcao/module_ri/Matrix_Orbs11.cpp b/source/source_lcao/module_ri/Matrix_Orbs11.cpp index 906a99f008..942f5a4cc6 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs11.cpp +++ b/source/source_lcao/module_ri/Matrix_Orbs11.cpp @@ -9,39 +9,31 @@ #include "source_base/tool_title.h" #include "source_pw/module_pwdft/global.h" -int Matrix_Orbs11::init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, - const double rmax) +void Matrix_Orbs11::init( + const std::vector>>& orb_A, + const std::vector>>& orb_B, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmax) { ModuleBase::TITLE("Matrix_Orbs11", "init"); ModuleBase::timer::tick("Matrix_Orbs11", "init"); this->lat0 = &ucell.lat0; - //const int ntype = orb.get_ntype(); - //int lmax_orb = -1, lmax_beta = -1; - //for (int it = 0; it < ntype; it++) - //{ - // lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - // lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); - //} - int Lmax, Lmax_used; - //if(mode==1) - // { std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_1(lmax_orb, lmax_beta); } - if(mode==2) - { std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_2(GlobalC::exx_info.info_ri.abfs_Lmax); } - //else if(mode==3) - // { std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_3(lmax_orb); } - else - { throw std::invalid_argument("mode = "+std::to_string(mode)+"in file "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } + const int Lmax = std::max({ Exx_Abfs::get_Lmax(orb_A), Exx_Abfs::get_Lmax(orb_B) }); + const int Lmax_used = Exx_Abfs::get_Lmax(orb_A) + Exx_Abfs::get_Lmax(orb_B); //========================================= // (3) make Gaunt coefficients table //========================================= - // this->MGT.init_Gaunt_CH(Lmax); - // this->MGT.init_Gaunt(Lmax); + if(!this->MGT) + { this->MGT = std::make_shared(); } + if(this->MGT->get_Lmax_Gaunt_CH() < Lmax) + { this->MGT->init_Gaunt_CH(Lmax); } + if(this->MGT->get_Lmax_Gaunt_Coefficients() < Lmax) + { this->MGT->init_Gaunt(Lmax); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -55,16 +47,6 @@ int Matrix_Orbs11::init(const int mode, Rmesh, psb_); - ModuleBase::timer::tick("Matrix_Orbs11", "init"); - return Lmax; -} - -void Matrix_Orbs11::init_radial(const std::vector>>& orb_A, - const std::vector>>& orb_B, - const ORB_gaunt_table& MGT) -{ - ModuleBase::TITLE("Matrix_Orbs11", "init_radial"); - ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); for (size_t TA = 0; TA != orb_A.size(); ++TA) { for (size_t TB = 0; TB != orb_B.size(); ++TB) { for (int LA = 0; LA != orb_A[TA].size(); ++LA) { @@ -73,17 +55,14 @@ void Matrix_Orbs11::init_radial(const std::vectorMGT))); + }}}}}} + + ModuleBase::timer::tick("Matrix_Orbs11", "init"); } -void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT) +/* +void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B) { ModuleBase::TITLE("Matrix_Orbs11", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); @@ -98,7 +77,7 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& Center2_Orb::Orb11(orb_A.Phi[TA].PhiLN(LA, NA), orb_B.Phi[TB].PhiLN(LB, NB), psb_, - MGT))); + *this->MGT))); } } } @@ -107,6 +86,7 @@ void Matrix_Orbs11::init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& } ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); } +*/ void Matrix_Orbs11::init_radial_table() { diff --git a/source/source_lcao/module_ri/Matrix_Orbs11.h b/source/source_lcao/module_ri/Matrix_Orbs11.h index a90f0d3fed..069127c355 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs11.h +++ b/source/source_lcao/module_ri/Matrix_Orbs11.h @@ -21,19 +21,13 @@ class Matrix_Orbs11 { public: - // mode: - // 1: - // 2: - int init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, // extend Kcut, keep dK - const double rmax); - - void init_radial(const std::vector>>& orb_A, - const std::vector>>& orb_B, - const ORB_gaunt_table& MGT); - void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B, const ORB_gaunt_table& MGT); + void init( + const std::vector>>& orb_A, + const std::vector>>& orb_B, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, // extend Kcut, keep dK + const double rmax); // extend Rcut, keep dR void init_radial_table(); void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 @@ -67,6 +61,8 @@ class Matrix_Orbs11 const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM& index_r, const ModuleBase::Element_Basis_Index::IndexLNM& index_c) const; + + std::shared_ptr MGT; private: ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; diff --git a/source/source_lcao/module_ri/Matrix_Orbs21.cpp b/source/source_lcao/module_ri/Matrix_Orbs21.cpp index 2814067f72..a11a1333b5 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs21.cpp +++ b/source/source_lcao/module_ri/Matrix_Orbs21.cpp @@ -9,32 +9,33 @@ #include "source_base/tool_title.h" #include "source_pw/module_pwdft/global.h" -int Matrix_Orbs21::init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, - const double rmax) +void Matrix_Orbs21::init( + const std::vector>>& orb_A1, + const std::vector>>& orb_A2, + const std::vector>>& orb_B, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmax) { ModuleBase::TITLE("Matrix_Orbs21", "init"); ModuleBase::timer::tick("Matrix_Orbs21", "init"); this->lat0 = &ucell.lat0; - const int ntype = orb.get_ntype(); - int lmax_orb = -1; - for (int it = 0; it < ntype; it++) - { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); } - int Lmax, Lmax_used; - if(mode==1) - { std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_3_1(GlobalC::exx_info.info_ri.abfs_Lmax, lmax_orb); } - else - { throw std::invalid_argument("mode = "+std::to_string(mode)+"in file "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } + const int Lmax = std::max({ + Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2), + Exx_Abfs::get_Lmax(orb_B) }); + const int Lmax_used = Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2) + Exx_Abfs::get_Lmax(orb_B) + 1; //========================================= // (3) make Gaunt coefficients table //========================================= - // this->MGT.init_Gaunt_CH(2 * Lmax + 1); // why +1 - // this->MGT.init_Gaunt(2 * Lmax + 1); - Lmax = 2 * Lmax + 1; + if(!this->MGT) + { this->MGT = std::make_shared(); } + if(this->MGT->get_Lmax_Gaunt_CH() < Lmax) + { this->MGT->init_Gaunt_CH(Lmax); } + if(this->MGT->get_Lmax_Gaunt_Coefficients() < Lmax) + { this->MGT->init_Gaunt(Lmax); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -48,56 +49,32 @@ int Matrix_Orbs21::init(const int mode, Rmesh, psb_); - ModuleBase::timer::tick("Matrix_Orbs21", "init"); - return Lmax; -} - -void Matrix_Orbs21::init_radial(const std::vector>>& orb_A1, - const std::vector>>& orb_A2, - const std::vector>>& orb_B, - const ORB_gaunt_table& MGT) -{ - ModuleBase::TITLE("Matrix_Orbs21", "init_radial"); - ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); assert(orb_A1.size() == orb_A2.size()); - for (size_t TA = 0; TA != orb_A1.size(); ++TA) - { - for (size_t TB = 0; TB != orb_B.size(); ++TB) - { - for (int LA1 = 0; LA1 != orb_A1[TA].size(); ++LA1) - { - for (size_t NA1 = 0; NA1 != orb_A1[TA][LA1].size(); ++NA1) - { - for (int LA2 = 0; LA2 != orb_A2[TA].size(); ++LA2) - { - for (size_t NA2 = 0; NA2 != orb_A2[TA][LA2].size(); ++NA2) - { - for (int LB = 0; LB != orb_B[TB].size(); ++LB) - { - for (size_t NB = 0; NB != orb_B[TB][LB].size(); ++NB) - { + for (size_t TA = 0; TA != orb_A1.size(); ++TA) { + for (size_t TB = 0; TB != orb_B.size(); ++TB) { + for (int LA1 = 0; LA1 != orb_A1[TA].size(); ++LA1) { + for (size_t NA1 = 0; NA1 != orb_A1[TA][LA1].size(); ++NA1) { + for (int LA2 = 0; LA2 != orb_A2[TA].size(); ++LA2) { + for (size_t NA2 = 0; NA2 != orb_A2[TA][LA2].size(); ++NA2) { + for (int LB = 0; LB != orb_B[TB].size(); ++LB) { + for (size_t NB = 0; NB != orb_B[TB][LB].size(); ++NB) { center2_orb21_s[TA][TB][LA1][NA1][LA2][NA2][LB].insert( - std::make_pair(NB, - Center2_Orb::Orb21(orb_A1[TA][LA1][NA1], - orb_A2[TA][LA2][NA2], - orb_B[TB][LB][NB], - psb_, - MGT))); - } - } - } - } - } - } - } - } - ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); + std::make_pair( + NB, + Center2_Orb::Orb21( + orb_A1[TA][LA1][NA1], + orb_A2[TA][LA2][NA2], + orb_B[TB][LB][NB], + psb_, + *this->MGT))); + }}}}}}}} + ModuleBase::timer::tick("Matrix_Orbs21", "init"); } +/* void Matrix_Orbs21::init_radial(const std::vector>>& orb_A1, const LCAO_Orbitals& orb_A2, - const LCAO_Orbitals& orb_B, - const ORB_gaunt_table& MGT) + const LCAO_Orbitals& orb_B) { ModuleBase::TITLE("Matrix_Orbs21", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); @@ -124,7 +101,7 @@ void Matrix_Orbs21::init_radial(const std::vectorMGT))); } } } @@ -135,6 +112,7 @@ void Matrix_Orbs21::init_radial(const std::vector - int init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, // extend Kcut, keep dK - const double rmax); // extend Rcut, keep dR - - void init_radial(const std::vector>>& orb_A1, - const std::vector>>& orb_A2, - const std::vector>>& orb_B, - const ORB_gaunt_table& MGT); - void init_radial(const std::vector>>& orb_A1, - const LCAO_Orbitals& orb_A2, - const LCAO_Orbitals& orb_B, - const ORB_gaunt_table& MGT); + void init( + const std::vector>>& orb_A1, + const std::vector>>& orb_A2, + const std::vector>>& orb_B, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, // extend Kcut, keep dK + const double rmax); // extend Rcut, keep dR void init_radial_table(); void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 @@ -76,6 +68,8 @@ class Matrix_Orbs21 const ModuleBase::Element_Basis_Index::IndexLNM& index_A1, const ModuleBase::Element_Basis_Index::IndexLNM& index_A2, const ModuleBase::Element_Basis_Index::IndexLNM& index_B) const; + + std::shared_ptr MGT; private: ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; diff --git a/source/source_lcao/module_ri/Matrix_Orbs22.cpp b/source/source_lcao/module_ri/Matrix_Orbs22.cpp index 3c1d0a1e64..69ea735fbc 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs22.cpp +++ b/source/source_lcao/module_ri/Matrix_Orbs22.cpp @@ -9,33 +9,35 @@ #include "source_base/tool_title.h" #include "source_pw/module_pwdft/global.h" -int Matrix_Orbs22::init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, - const double rmax) +void Matrix_Orbs22::init( + const std::vector>>& orb_A1, + const std::vector>>& orb_A2, + const std::vector>>& orb_B1, + const std::vector>>& orb_B2, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmax) { ModuleBase::TITLE("Matrix_Orbs22", "init"); ModuleBase::timer::tick("Matrix_Orbs22", "init"); - this->lat0 = &ucell.lat0; + this->lat0 = &ucell.lat0; - const int ntype = orb.get_ntype(); - int lmax_orb = -1; - for (int it = 0; it < ntype; it++) - { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); } - int Lmax, Lmax_used; - if(mode==1) - { std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_4_1(lmax_orb); } - else - { throw std::invalid_argument("mode = "+std::to_string(mode)+"in file "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); } + const int Lmax = std::max({ + Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2), + Exx_Abfs::get_Lmax(orb_B1) + Exx_Abfs::get_Lmax(orb_B2) }); + const int Lmax_used = Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2) + Exx_Abfs::get_Lmax(orb_B1) + Exx_Abfs::get_Lmax(orb_B2) + 2; //========================================= // (3) make Gaunt coefficients table //========================================= - // this->MGT.init_Gaunt_CH(2 * Lmax + 1); // why +1 - // this->MGT.init_Gaunt(2 * Lmax + 1); - Lmax = 2 * Lmax + 1; + if(!this->MGT) + { this->MGT = std::make_shared(); } + if(this->MGT->get_Lmax_Gaunt_CH() < Lmax) + { this->MGT->init_Gaunt_CH(Lmax); } + if(this->MGT->get_Lmax_Gaunt_Coefficients() < Lmax) + { this->MGT->init_Gaunt(Lmax); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -49,46 +51,37 @@ int Matrix_Orbs22::init(const int mode, Rmesh, psb_); - ModuleBase::timer::tick("Matrix_Orbs22", "init"); - return Lmax; -} - -void Matrix_Orbs22::init_radial(const std::vector>>& orb_A1, - const std::vector>>& orb_A2, - const std::vector>>& orb_B1, - const std::vector>>& orb_B2, - const ORB_gaunt_table& MGT) -{ - ModuleBase::TITLE("Matrix_Orbs22", "init_radial"); - ModuleBase::timer::tick("Matrix_Orbs22", "init_radial"); assert(orb_A1.size() == orb_A2.size()); assert(orb_B1.size() == orb_B2.size()); - for (size_t TA = 0; TA != orb_A1.size(); ++TA) - for (size_t TB = 0; TB != orb_B1.size(); ++TB) - for (int LA1 = 0; LA1 != orb_A1[TA].size(); ++LA1) - for (size_t NA1 = 0; NA1 != orb_A1[TA][LA1].size(); ++NA1) - for (int LA2 = 0; LA2 != orb_A2[TA].size(); ++LA2) - for (size_t NA2 = 0; NA2 != orb_A2[TA][LA2].size(); ++NA2) - for (int LB1 = 0; LB1 != orb_B1[TB].size(); ++LB1) - for (size_t NB1 = 0; NB1 != orb_B1[TB][LB1].size(); ++NB1) - for (int LB2 = 0; LB2 != orb_B2[TB].size(); ++LB2) - for (size_t NB2 = 0; NB2 != orb_B2[TB][LB2].size(); ++NB2) + for (size_t TA = 0; TA != orb_A1.size(); ++TA) { + for (size_t TB = 0; TB != orb_B1.size(); ++TB) { + for (int LA1 = 0; LA1 != orb_A1[TA].size(); ++LA1) { + for (size_t NA1 = 0; NA1 != orb_A1[TA][LA1].size(); ++NA1) { + for (int LA2 = 0; LA2 != orb_A2[TA].size(); ++LA2) { + for (size_t NA2 = 0; NA2 != orb_A2[TA][LA2].size(); ++NA2) { + for (int LB1 = 0; LB1 != orb_B1[TB].size(); ++LB1) { + for (size_t NB1 = 0; NB1 != orb_B1[TB][LB1].size(); ++NB1) { + for (int LB2 = 0; LB2 != orb_B2[TB].size(); ++LB2) { + for (size_t NB2 = 0; NB2 != orb_B2[TB][LB2].size(); ++NB2) { center2_orb22_s[TA][TB][LA1][NA1][LA2][NA2][LB1][NB1][LB2].insert( - std::make_pair(NB2, - Center2_Orb::Orb22(orb_A1[TA][LA1][NA1], - orb_A2[TA][LA2][NA2], - orb_B1[TB][LB1][NB1], - orb_B2[TB][LB2][NB2], - psb_, - MGT))); - ModuleBase::timer::tick("Matrix_Orbs22", "init_radial"); + std::make_pair( + NB2, + Center2_Orb::Orb22( + orb_A1[TA][LA1][NA1], + orb_A2[TA][LA2][NA2], + orb_B1[TB][LB1][NB1], + orb_B2[TB][LB2][NB2], + psb_, + *this->MGT))); + }}}}}}}}}} + ModuleBase::timer::tick("Matrix_Orbs22", "init"); } +/* void Matrix_Orbs22::init_radial(const LCAO_Orbitals& orb_A1, const LCAO_Orbitals& orb_A2, const LCAO_Orbitals& orb_B1, - const LCAO_Orbitals& orb_B2, - const ORB_gaunt_table& MGT) + const LCAO_Orbitals& orb_B2) { ModuleBase::TITLE("Matrix_Orbs22", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs22", "init_radial"); @@ -111,9 +104,10 @@ void Matrix_Orbs22::init_radial(const LCAO_Orbitals& orb_A1, orb_B1.Phi[TB].PhiLN(LB1, NB1), orb_B2.Phi[TB].PhiLN(LB2, NB2), psb_, - MGT))); + *this->MGT))); ModuleBase::timer::tick("Matrix_Orbs22", "init_radial"); } +*/ void Matrix_Orbs22::init_radial_table() { diff --git a/source/source_lcao/module_ri/Matrix_Orbs22.h b/source/source_lcao/module_ri/Matrix_Orbs22.h index 4e9f49f628..e0c3790289 100644 --- a/source/source_lcao/module_ri/Matrix_Orbs22.h +++ b/source/source_lcao/module_ri/Matrix_Orbs22.h @@ -21,24 +21,15 @@ class Matrix_Orbs22 { public: - // mode: - // 1: - int init(const int mode, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const double kmesh_times, // extend Kcut, keep dK - const double rmax); // extend Rcut, keep dR - - void init_radial(const std::vector>>& orb_A1, - const std::vector>>& orb_A2, - const std::vector>>& orb_B1, - const std::vector>>& , - const ORB_gaunt_table& MGT); - void init_radial(const LCAO_Orbitals& orb_A1, - const LCAO_Orbitals& orb_A2, - const LCAO_Orbitals& orb_B1, - const LCAO_Orbitals& orb_B2, - const ORB_gaunt_table& MGT); + void init( + const std::vector>>& orb_A1, + const std::vector>>& orb_A2, + const std::vector>>& orb_B1, + const std::vector>>& orb_B2, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, // extend Kcut, keep dK + const double rmax); // extend Rcut, keep dR void init_radial_table(); void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 @@ -100,6 +91,8 @@ class Matrix_Orbs22 const ModuleBase::Element_Basis_Index::IndexLNM& index_A2, const ModuleBase::Element_Basis_Index::IndexLNM& index_B1, const ModuleBase::Element_Basis_Index::IndexLNM& index_B2) const; + + std::shared_ptr MGT; private: ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; diff --git a/source/source_lcao/module_ri/RPA_LRI.hpp b/source/source_lcao/module_ri/RPA_LRI.hpp index d17479b239..641e07a7f4 100644 --- a/source/source_lcao/module_ri/RPA_LRI.hpp +++ b/source/source_lcao/module_ri/RPA_LRI.hpp @@ -95,7 +95,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, Rs, period, ucell.lat); // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF - symrot.set_abfs_Lmax(GlobalC::exx_info.info_ri.abfs_Lmax); + symrot.set_abfs_Lmax(Exx_Abfs::get_Lmax(exx_lri_rpa.abfs)); symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer()); // output Ts (symrot_R.txt) and Ms (symrot_k.txt) ModuleSymmetry::print_symrot_info_R(symrot, ucell.symm, ucell.lmax, Rs); diff --git a/source/source_lcao/module_ri/exx_abfs.h b/source/source_lcao/module_ri/exx_abfs.h index f08462c94e..42cc41a01d 100644 --- a/source/source_lcao/module_ri/exx_abfs.h +++ b/source/source_lcao/module_ri/exx_abfs.h @@ -23,6 +23,13 @@ class Exx_Abfs int rmesh_times = 5; // Peize Lin test int kmesh_times = 1; // Peize Lin test + static int get_Lmax(const std::vector>>& orb) + { + int Lmax = -1; + for( const auto &orb_T : orb ) + { Lmax = std::max( Lmax, static_cast(orb_T.size())-1 ); } + return Lmax; + } }; #endif diff --git a/source/source_lcao/module_ri/exx_opt_orb.cpp b/source/source_lcao/module_ri/exx_opt_orb.cpp index 4ba6569107..83ea26c18d 100644 --- a/source/source_lcao/module_ri/exx_opt_orb.cpp +++ b/source/source_lcao/module_ri/exx_opt_orb.cpp @@ -47,12 +47,6 @@ void Exx_Opt_Orb::generate_matrix( { jle = Exx_Abfs::IO::construct_abfs( jle, orb, info.files_jles, info.kmesh_times ); } Exx_Abfs::Construct_Orbs::filter_empty_orbs(jle); - // GlobalC::exx_info.info_ri.abfs_Lmax is used for Center2 temporarily - for(const auto &orb_T : abfs) - { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(orb_T.size())-1 ); } - for(const auto &orb_T : jle) - { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast(orb_T.size())-1 ); } - const ModuleBase::Element_Basis_Index::Range range_lcaos = ModuleBase::Element_Basis_Index::construct_range( lcaos ); const ModuleBase::Element_Basis_Index::IndexLNM index_lcaos = ModuleBase::Element_Basis_Index::construct_index( range_lcaos ); @@ -72,11 +66,7 @@ void Exx_Opt_Orb::generate_matrix( { if(judge_orbs_empty(lcaos)) { return {}; } Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos; - ORB_gaunt_table MGT; - const int Lmax = m_lcaoslcaos_lcaoslcaos.init( 1, ucell,orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_lcaoslcaos_lcaoslcaos.init_radial( lcaos, lcaos, lcaos, lcaos, MGT ); + m_lcaoslcaos_lcaoslcaos.init( lcaos, lcaos, lcaos, lcaos, ucell,orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_lcaoslcaos_lcaoslcaos.init_radial_table(radial_R); #else @@ -91,11 +81,7 @@ void Exx_Opt_Orb::generate_matrix( if(judge_orbs_empty(lcaos)) { return {}; } if(judge_orbs_empty(jle)) { return {}; } Matrix_Orbs21 m_jyslcaos_lcaos; - ORB_gaunt_table MGT; - const int Lmax = m_jyslcaos_lcaos.init( 1, ucell , orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_jyslcaos_lcaos.init_radial( jle, lcaos, lcaos, MGT); + m_jyslcaos_lcaos.init( jle, lcaos, lcaos, ucell , orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_jyslcaos_lcaos.init_radial_table( radial_R); #else @@ -109,11 +95,7 @@ void Exx_Opt_Orb::generate_matrix( { if(judge_orbs_empty(jle)) { return {}; } Matrix_Orbs11 m_jys_jys; - ORB_gaunt_table MGT; - const int Lmax = m_jys_jys.init( 2,ucell,orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_jys_jys.init_radial( jle, jle, MGT ); + m_jys_jys.init( jle, jle, ucell,orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_jys_jys.init_radial_table(radial_R); #else @@ -127,11 +109,7 @@ void Exx_Opt_Orb::generate_matrix( { if(judge_orbs_empty(abfs)) { return {}; } Matrix_Orbs11 m_abfs_abfs; - ORB_gaunt_table MGT; - const int Lmax = m_abfs_abfs.init( 2, ucell, orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_abfs_abfs.init_radial( abfs, abfs, MGT ); + m_abfs_abfs.init( abfs, abfs, ucell, orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_abfs_abfs.init_radial_table(radial_R); #else @@ -146,11 +124,7 @@ void Exx_Opt_Orb::generate_matrix( if(judge_orbs_empty(lcaos)) { return {}; } if(judge_orbs_empty(abfs)) { return {}; } Matrix_Orbs21 m_abfslcaos_lcaos; - ORB_gaunt_table MGT; - const int Lmax = m_abfslcaos_lcaos.init( 1, ucell , orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos, MGT ); + m_abfslcaos_lcaos.init( abfs, lcaos, lcaos, ucell , orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_abfslcaos_lcaos.init_radial_table(radial_R); #else @@ -165,11 +139,7 @@ void Exx_Opt_Orb::generate_matrix( if(judge_orbs_empty(jle)) { return {}; } if(judge_orbs_empty(abfs)) { return {}; } Matrix_Orbs11 m_jys_abfs; - ORB_gaunt_table MGT; - const int Lmax = m_jys_abfs.init( 2, ucell,orb, info.kmesh_times, orb.get_Rmax() ); - MGT.init_Gaunt_CH(Lmax); - MGT.init_Gaunt(Lmax); - m_jys_abfs.init_radial( jle, abfs, MGT ); + m_jys_abfs.init( jle, abfs, ucell,orb, info.kmesh_times, orb.get_Rmax() ); #if TEST_EXX_RADIAL>=1 m_jys_abfs.init_radial_table(radial_R); #else