diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 217f14e89d..16f35eea7f 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -32,6 +32,7 @@ #include #ifdef __EXX #include "module_io/restart_exx_csr.h" +#include "module_ri/exx_opt_orb.h" #include "module_ri/RPA_LRI.h" #endif @@ -170,9 +171,19 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell inp.lcao_rmax, ucell, two_center_bundle_, - orb_); + this->orb_); //------------------init Basis_lcao---------------------- + if (PARAM.inp.calculation == "gen_opt_abfs") + { +#ifdef __EXX + Exx_Opt_Orb::generate_matrix(GlobalC::exx_info.info_opt_abfs, this->kv, this->orb_); +#else + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_all_runners", "calculation=gen_opt_abfs must compile __EXX"); +#endif + return; + } + // 5) initialize density matrix // DensityMatrix is allocated here, DMK is also initialized here // DMR is not initialized here, it will be constructed in each before_scf @@ -188,7 +199,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 6) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, pv, this->kv.get_nks(), orb_); + LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, pv, this->kv.get_nks(), this->orb_); #ifdef __EXX // 7) initialize exx @@ -202,11 +213,11 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // initialize 2-center radial tables for EXX-LRI if (GlobalC::exx_info.info_ri.real_number) { - this->exx_lri_double->init(MPI_COMM_WORLD, this->kv, orb_); + this->exx_lri_double->init(MPI_COMM_WORLD, this->kv, this->orb_); } else { - this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, orb_); + this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, this->orb_); } } } @@ -215,7 +226,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // 8) initialize DFT+U if (PARAM.inp.dft_plus_u) { - GlobalC::dftu.init(ucell, &this->pv, this->kv.get_nks(), orb_); + GlobalC::dftu.init(ucell, &this->pv, this->kv.get_nks(), this->orb_); } // 9) initialize ppcell @@ -244,7 +255,7 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell // load the DeePKS model from deep neural network GlobalC::ld.load_model(PARAM.inp.deepks_model); // read pdm from file for NSCF or SCF-restart, do it only once in whole calculation - GlobalC::ld.read_projected_DM((PARAM.inp.init_chg == "file"), PARAM.inp.deepks_equiv, *orb_.Alpha); + GlobalC::ld.read_projected_DM((PARAM.inp.init_chg == "file"), PARAM.inp.deepks_equiv, *this->orb_.Alpha); } #endif @@ -313,7 +324,7 @@ void ESolver_KS_LCAO::cal_force(ModuleBase::matrix& force) this->GG, // mohan add 2024-04-01 this->GK, // mohan add 2024-04-01 two_center_bundle_, - orb_, + this->orb_, force, this->scs, this->sf, @@ -466,7 +477,7 @@ void ESolver_KS_LCAO::after_all_runners() this->GG, this->GK, this->kv, - orb_.cutoffs(), + this->orb_.cutoffs(), this->pelec->wg, GlobalC::GridD #ifdef __EXX @@ -495,7 +506,7 @@ void ESolver_KS_LCAO::after_all_runners() this->kv, this->pelec->wg, GlobalC::GridD, - orb_.cutoffs(), + this->orb_.cutoffs(), this->two_center_bundle_ #ifdef __EXX , @@ -1168,7 +1179,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->pelec->ekb, this->pelec->klist->kvec_d, GlobalC::ucell, - orb_, + this->orb_, GlobalC::GridD, &(this->pv), *(this->psi), @@ -1188,8 +1199,8 @@ void ESolver_KS_LCAO::after_scf(const int istep) rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), MPI_COMM_WORLD, this->kv, - orb_); - rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); + this->orb_); + rpa_lri_double.init(MPI_COMM_WORLD, this->kv, this->orb_.cutoffs()); rpa_lri_double.out_for_RPA(this->pv, *(this->psi), this->pelec); } #endif @@ -1278,7 +1289,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) this->kv.kvec_d, &hR, &GlobalC::ucell, - orb_.cutoffs(), + this->orb_.cutoffs(), &GlobalC::GridD, two_center_bundle_.kinetic_orb.get()); diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 667eac6c20..ea2a533646 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -56,6 +56,10 @@ void ESolver_KS_LCAO::others(const int istep) // return; // use 'return' will cause segmentation fault. by mohan // 2024-06-09 } + else if (cal_type == "gen_opt_abfs") + { + return; + } else if (cal_type == "test_memory") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing memory"); diff --git a/source/module_hamilt_general/module_xc/exx_info.h b/source/module_hamilt_general/module_xc/exx_info.h index d94f5e6375..1c4b75a313 100644 --- a/source/module_hamilt_general/module_xc/exx_info.h +++ b/source/module_hamilt_general/module_xc/exx_info.h @@ -23,20 +23,20 @@ struct Exx_Info struct Exx_Info_Lip { - const Conv_Coulomb_Pot_K::Ccp_Type& ccp_type; - const double& hse_omega; + const Conv_Coulomb_Pot_K::Ccp_Type &ccp_type; + const double &hse_omega; double lambda = 0.3; Exx_Info_Lip(const Exx_Info::Exx_Info_Global& info_global) :ccp_type(info_global.ccp_type), - hse_omega(info_global.hse_omega) {} + hse_omega(info_global.hse_omega) {} }; Exx_Info_Lip info_lip; struct Exx_Info_RI { - const Conv_Coulomb_Pot_K::Ccp_Type& ccp_type; - const double& hse_omega; + const Conv_Coulomb_Pot_K::Ccp_Type &ccp_type; + const double &hse_omega; bool real_number = false; @@ -58,15 +58,28 @@ struct Exx_Info int abfs_Lmax = 0; // tmp Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global) - : ccp_type(info_global.ccp_type), hse_omega(info_global.hse_omega) - { - } + :ccp_type(info_global.ccp_type), + hse_omega(info_global.hse_omega) {} }; Exx_Info_RI info_ri; - Exx_Info() : info_lip(this->info_global), info_ri(this->info_global) + struct Exx_Info_Opt_ABFs { - } + //const Conv_Coulomb_Pot_K::Ccp_Type &ccp_type; + //const double &hse_omega; + + int abfs_Lmax = 0; // tmp + double ecut_exx = 60; + double tolerence = 1E-2; + double kmesh_times = 4; + + //Exx_Info_Opt_ABFs(const Exx_Info::Exx_Info_Global& info_global) + // :ccp_type(info_global.ccp_type), + // hse_omega(info_global.hse_omega) {} + }; + Exx_Info_Opt_ABFs info_opt_abfs; + + Exx_Info() : info_lip(this->info_global), info_ri(this->info_global) {} }; #endif diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 713ce512fc..5f0b60bc72 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -192,7 +192,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 4; use_libxc = false; } - else if( xc_func == "OPT_ORB" || xc_func == "NONE" || xc_func == "NOX+NOC") + else if( xc_func == "NONE" || xc_func == "NOX+NOC") { // not doing anything } diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 732b1e27ec..db5063a71b 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -326,7 +326,6 @@ void Input_Conv::Convert() ModuleBase::GlobalFunc::MAKE_DIR(GlobalC::restart.folder); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" - || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0") { GlobalC::restart.info_save.save_charge = true; GlobalC::restart.info_save.save_H = true; @@ -344,7 +343,6 @@ void Input_Conv::Convert() GlobalC::restart.folder = PARAM.globalv.global_readin_dir + "restart/"; if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" - || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0") { GlobalC::restart.info_load.load_charge = true; GlobalC::restart.info_load.load_H = true; @@ -373,14 +371,11 @@ void Input_Conv::Convert() GlobalC::exx_info.info_global.cal_exx = true; GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; - } else if (dft_functional_lower == "opt_orb") { - GlobalC::exx_info.info_global.cal_exx = false; - Exx_Abfs::Jle::generate_matrix = true; } else { GlobalC::exx_info.info_global.cal_exx = false; } - if (GlobalC::exx_info.info_global.cal_exx || Exx_Abfs::Jle::generate_matrix || PARAM.inp.rpa) + if (GlobalC::exx_info.info_global.cal_exx || PARAM.inp.rpa) { // EXX case, convert all EXX related variables // GlobalC::exx_info.info_global.cal_exx = true; @@ -407,9 +402,9 @@ void Input_Conv::Convert() GlobalC::exx_info.info_ri.cauchy_stress_threshold = PARAM.inp.exx_cauchy_stress_threshold; GlobalC::exx_info.info_ri.ccp_rmesh_times = std::stod(PARAM.inp.exx_ccp_rmesh_times); - Exx_Abfs::Jle::Lmax = PARAM.inp.exx_opt_orb_lmax; - Exx_Abfs::Jle::Ecut_exx = PARAM.inp.exx_opt_orb_ecut; - Exx_Abfs::Jle::tolerence = PARAM.inp.exx_opt_orb_tolerence; + GlobalC::exx_info.info_opt_abfs.abfs_Lmax = PARAM.inp.exx_opt_orb_lmax; + GlobalC::exx_info.info_opt_abfs.ecut_exx = PARAM.inp.exx_opt_orb_ecut; + GlobalC::exx_info.info_opt_abfs.tolerence = PARAM.inp.exx_opt_orb_tolerence; // EXX does not support symmetry for nspin==4 if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1" && PARAM.inp.nspin == 4) diff --git a/source/module_io/read_input_item_system.cpp b/source/module_io/read_input_item_system.cpp index 32ff4bccef..8eccc300da 100644 --- a/source/module_io/read_input_item_system.cpp +++ b/source/module_io/read_input_item_system.cpp @@ -61,7 +61,7 @@ void ReadInput::item_system() } { Input_Item item("calculation"); - item.annotation = "test; scf; relax; nscf; get_wf; get_pchg"; + item.annotation = "scf; relax; md; cell-relax; nscf; get_S; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; item.read_value = [](const Input_Item& item, Parameter& para) { para.input.calculation = strvalue; std::string& calculation = para.input.calculation; @@ -78,7 +78,8 @@ void ReadInput::item_system() "get_S", "get_wf", "get_pchg", - "gen_bessel"}; + "gen_bessel", + "gen_opt_abfs"}; if (!find_str(callist, calculation)) { const std::string warningstr = nofound_str(callist, "calculation"); diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index bb04a3f358..6b8d3a0668 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -21,7 +21,7 @@ #include "module_exx_symmetry/symmetry_rotation.h" class Parallel_Orbitals; - + template class RPA_LRI; @@ -59,19 +59,19 @@ class Exx_LRI void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb); void cal_exx_force(); void cal_exx_stress(); - std::vector> get_abfs_nchis() const; + // std::vector> get_abfs_nchis() const; std::vector< std::map>>> Hexxs; double Eexx; ModuleBase::matrix force_exx; ModuleBase::matrix stress_exx; - + private: const Exx_Info::Exx_Info_RI &info; MPI_Comm mpi_comm; const K_Vectors *p_kv = nullptr; - std::vector orb_cutoff_; + std::vector orb_cutoff_; std::vector>> lcaos; std::vector>> abfs; @@ -99,4 +99,4 @@ class Exx_LRI #include "Exx_LRI.hpp" -#endif +#endif diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 58a900fe11..808627801b 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -47,9 +47,9 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c // Hexx_para.mixing_beta = GlobalC::CHR.mixing_beta; // } - this->mpi_comm = mpi_comm_in; - this->p_kv = &kv_in; - this->orb_cutoff_ = orb.cutoffs(); + this->mpi_comm = mpi_comm_in; + this->p_kv = &kv_in; + this->orb_cutoff_ = orb.cutoffs(); this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times ); @@ -59,11 +59,10 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c const std::vector>> abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); - if(this->info.files_abfs.empty()) { - 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 ); -} + if(this->info.files_abfs.empty()) + { 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::print_orbs_size(this->abfs, GlobalV::ofs_running); auto get_ccp_parameter = [this]() -> std::map @@ -85,21 +84,22 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c throw std::domain_error(std::string(__FILE__)+" line "+std::to_string(__LINE__)); break; } }; - this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times); + this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times); - 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 ); -} + 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->cv.set_orbitals( - orb, + orb, this->lcaos, this->abfs, this->abfs_ccp, this->info.kmesh_times, this->info.ccp_rmesh_times ); ModuleBase::timer::tick("Exx_LRI", "init"); } + + template void Exx_LRI::cal_exx_ions(const bool write_cv) { @@ -112,13 +112,11 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) // this->m_abfslcaos_lcaos.init_radial_table(Rradial); std::vector atoms(GlobalC::ucell.nat); - for(int iat=0; iat atoms_pos; - for(int iat=0; iat latvec = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), RI_Util::Vector3_to_array3(GlobalC::ucell.a2), @@ -128,7 +126,7 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period); // std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour. - const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1+this->info.ccp_rmesh_times, orb_cutoff_); + const std::array period_Vs = LRI_CV_Tools::cal_latvec_range(1+this->info.ccp_rmesh_times, this->orb_cutoff_); const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false); @@ -137,7 +135,8 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) list_As_Vs.first, list_As_Vs.second[0], {{"writable_Vws",true}}); this->cv.Vws = LRI_CV_Tools::get_CVws(Vs); - if (write_cv && GlobalV::MY_RANK == 0) { LRI_CV_Tools::write_Vs_abf(Vs, PARAM.globalv.global_out_dir + "Vs"); } + if (write_cv && GlobalV::MY_RANK == 0) + { LRI_CV_Tools::write_Vs_abf(Vs, PARAM.globalv.global_out_dir + "Vs"); } this->exx_lri.set_Vs(std::move(Vs), this->info.V_threshold); if(PARAM.inp.cal_force || PARAM.inp.cal_stress) @@ -155,7 +154,7 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) } } - const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2, orb_cutoff_); + const std::array period_Cs = LRI_CV_Tools::cal_latvec_range(2, this->orb_cutoff_); const std::pair, std::vector>>>> list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false); @@ -166,7 +165,8 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) {"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}}); std::map>> &Cs = std::get<0>(Cs_dCs); this->cv.Cws = LRI_CV_Tools::get_CVws(Cs); - if (write_cv && GlobalV::MY_RANK == 0) { LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs"); } + if (write_cv && GlobalV::MY_RANK == 0) + { LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs"); } this->exx_lri.set_Cs(std::move(Cs), this->info.C_threshold); if(PARAM.inp.cal_force || PARAM.inp.cal_stress) @@ -183,51 +183,57 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions"); } + + template void Exx_LRI::cal_exx_elec(const std::vector>>>& Ds, - const Parallel_Orbitals& pv, - const ModuleSymmetry::Symmetry_rotation* p_symrot) + const Parallel_Orbitals& pv, + const ModuleSymmetry::Symmetry_rotation* p_symrot) { ModuleBase::TITLE("Exx_LRI","cal_exx_elec"); ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec"); + if (p_symrot) + { this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()); } + else + { this->exx_lri.set_symmetry(false, {}); } + const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(pv); this->Hexxs.resize(PARAM.inp.nspin); this->Eexx = 0; - (p_symrot) ? this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()) : this->exx_lri.set_symmetry(false, {}); for(int is=0; isexx_lri.set_Ds(Ds[is], this->info.dm_threshold, suffix); - this->exx_lri.cal_Hs({ "","",suffix }); - - if (!p_symrot) - { - this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first( - this->mpi_comm, std::move(this->exx_lri.Hs), std::get<0>(judge[is]), std::get<1>(judge[is])); - } - else - { - // reduce but not repeat - auto Hs_a2D = this->exx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs); - // rotate locally without repeat - Hs_a2D = p_symrot->restore_HR(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', Hs_a2D); - // cal energy using full Hs without repeat - this->exx_lri.energy = this->exx_lri.post_2D.cal_energy( - this->exx_lri.post_2D.saves["Ds_" + suffix], - this->exx_lri.post_2D.set_tensors_map2(Hs_a2D)); - // get repeated full Hs for abacus - this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first( - this->mpi_comm, std::move(Hs_a2D), std::get<0>(judge[is]), std::get<1>(judge[is])); - } - this->Eexx += std::real(this->exx_lri.energy); + std::string suffix = ((PARAM.inp.cal_force || PARAM.inp.cal_stress) ? std::to_string(is) : ""); + + this->exx_lri.set_Ds(Ds[is], this->info.dm_threshold, suffix); + this->exx_lri.cal_Hs({ "","",suffix }); + + if (!p_symrot) + { + this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first( + this->mpi_comm, std::move(this->exx_lri.Hs), std::get<0>(judge[is]), std::get<1>(judge[is])); + } + else + { + // reduce but not repeat + auto Hs_a2D = this->exx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs); + // rotate locally without repeat + Hs_a2D = p_symrot->restore_HR(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', Hs_a2D); + // cal energy using full Hs without repeat + this->exx_lri.energy = this->exx_lri.post_2D.cal_energy( + this->exx_lri.post_2D.saves["Ds_" + suffix], + this->exx_lri.post_2D.set_tensors_map2(Hs_a2D)); + // get repeated full Hs for abacus + this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first( + this->mpi_comm, std::move(Hs_a2D), std::get<0>(judge[is]), std::get<1>(judge[is])); + } + this->Eexx += std::real(this->exx_lri.energy); post_process_Hexx(this->Hexxs[is]); } this->Eexx = post_process_Eexx(this->Eexx); this->exx_lri.set_symmetry(false, {}); - ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec"); + ModuleBase::timer::tick("Exx_LRI", "cal_exx_elec"); } template @@ -245,8 +251,8 @@ template double Exx_LRI::post_process_Eexx(const double& Eexx_in) const { ModuleBase::TITLE("Exx_LRI","post_process_Eexx"); - const double SPIN_multiple = std::map{ {1,2}, {2,1}, {4,1} }.at(PARAM.inp.nspin); // why? - const double frac = -SPIN_multiple; + const double SPIN_multiple = std::map{ {1,2}, {2,1}, {4,1} }.at(PARAM.inp.nspin); // why? + const double frac = -SPIN_multiple; return frac * Eexx_in; } @@ -267,6 +273,8 @@ post_process_old } */ + + template void Exx_LRI::cal_exx_force() { @@ -280,8 +288,7 @@ void Exx_LRI::cal_exx_force() for(std::size_t idim=0; idimexx_lri.force[idim]) { this->force_exx(force_item.first, idim) += std::real(force_item.second); -} -} + } } } const double SPIN_multiple = std::map{{1,2}, {2,1}, {4,1}}.at(PARAM.inp.nspin); // why? @@ -291,6 +298,7 @@ void Exx_LRI::cal_exx_force() } + template void Exx_LRI::cal_exx_stress() { @@ -304,8 +312,7 @@ void Exx_LRI::cal_exx_stress() for(std::size_t idim0=0; idim0stress_exx(idim0,idim1) += std::real(this->exx_lri.stress(idim0,idim1)); -} -} + } } } const double SPIN_multiple = std::map{{1,2}, {2,1}, {4,1}}.at(PARAM.inp.nspin); // why? @@ -315,19 +322,22 @@ void Exx_LRI::cal_exx_stress() ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress"); } + +/* template std::vector> Exx_LRI::get_abfs_nchis() const { - std::vector> abfs_nchis; - for (const auto& abfs_T : this->abfs) - { - std::vector abfs_nchi_T; - for (const auto& abfs_L : abfs_T) { - abfs_nchi_T.push_back(abfs_L.size()); -} - abfs_nchis.push_back(abfs_nchi_T); - } - return abfs_nchis; + std::vector> abfs_nchis; + for (const auto& abfs_T : this->abfs) + { + std::vector abfs_nchi_T; + for (const auto& abfs_L : abfs_T) { + abfs_nchi_T.push_back(abfs_L.size()); + } + abfs_nchis.push_back(abfs_nchi_T); + } + return abfs_nchis; } +*/ #endif diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index d1b64c56df..15b1420a36 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -48,7 +48,7 @@ class Exx_LRI_Interface double& get_Eexx() const { return this->exx_ptr->Eexx; } // Processes in ESolver_KS_LCAO - /// @brief in beforescf: set xc type, opt_orb, do DM mixing + /// @brief in beforescf: set xc type, do DM mixing void exx_beforescf(const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const Parallel_2D& pv, const LCAO_Orbitals& orb); /// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 8fa6d8b84a..8d3b3c8428 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -4,7 +4,6 @@ #include "Exx_LRI_interface.h" #include "module_ri/exx_abfs-jle.h" -#include "module_ri/exx_opt_orb.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_base/parallel_common.h" @@ -42,17 +41,16 @@ void Exx_LRI_Interface::exx_beforescf(const K_Vectors& kv, const Charg #ifdef __MPI if (GlobalC::exx_info.info_global.cal_exx) { - if (GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx) { XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); - } else + if (GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx) + { + XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + } + else { if (ucell.atoms[0].ncpp.xc_func == "HF" || ucell.atoms[0].ncpp.xc_func == "PBE0" || ucell.atoms[0].ncpp.xc_func == "HSE") - { - XC_Functional::set_xc_type("pbe"); - } + { XC_Functional::set_xc_type("pbe"); } else if (ucell.atoms[0].ncpp.xc_func == "SCAN0") - { - XC_Functional::set_xc_type("scan"); - } + { XC_Functional::set_xc_type("scan"); } } // initialize the rotation matrix in AO representation this->exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1); @@ -68,30 +66,20 @@ void Exx_LRI_Interface::exx_beforescf(const K_Vectors& kv, const Charg this->exx_ptr->cal_exx_ions(PARAM.inp.out_ri_cv); } - if (Exx_Abfs::Jle::generate_matrix) - { - //program should be stopped after this judgement - Exx_Opt_Orb exx_opt_orb; - exx_opt_orb.generate_matrix(kv, orb); - ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); - return; - } - - // set initial parameter for mix_DMk_2D - 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);} - else - {this->mix_DMk_2D.set_nks(kv.get_nks(), PARAM.globalv.gamma_only_local);} - if(GlobalC::exx_info.info_global.separate_loop) { - this->mix_DMk_2D.set_mixing(nullptr); - } else { - this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); - } - // for exx two_level scf - this->two_level_step = 0; - } + // set initial parameter for mix_DMk_2D + 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); } + else + { this->mix_DMk_2D.set_nks(kv.get_nks(), PARAM.globalv.gamma_only_local); } + if(GlobalC::exx_info.info_global.separate_loop) + { this->mix_DMk_2D.set_mixing(nullptr); } + else + { this->mix_DMk_2D.set_mixing(chgmix.get_mixing()); } + // for exx two_level scf + this->two_level_step = 0; + } #endif // __MPI } @@ -103,14 +91,18 @@ void Exx_LRI_Interface::exx_eachiterinit(const elecstate::DensityMatri if (!GlobalC::exx_info.info_global.separate_loop && this->two_level_step) { const bool flag_restart = (iter == 1) ? 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); } - else { this->mix_DMk_2D.mix(dm.get_DMK_vector(), 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); } const std::vector>,RI::Tensor>>> Ds = PARAM.globalv.gamma_only_local ? RI_2D_Comm::split_m2D_ktoR(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin) : RI_2D_Comm::split_m2D_ktoR(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry); - if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer(), &this->symrot_); } - else { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer()); } + if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer(), &this->symrot_); } + else + { this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer()); } } } } @@ -126,8 +118,8 @@ void Exx_LRI_Interface::exx_hamilt2density(elecstate::ElecState& elec, if (GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx && this->two_level_step == 0 && iter == 1) { - if (GlobalV::MY_RANK == 0) {GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); -} + if (GlobalV::MY_RANK == 0) + { GlobalC::restart.load_disk("Eexx", 0, 1, &this->exx_ptr->Eexx); } Parallel_Common::bcast_double(this->exx_ptr->Eexx); this->exx_ptr->Eexx /= GlobalC::exx_info.info_global.hybrid_alpha; } @@ -150,79 +142,77 @@ bool Exx_LRI_Interface::exx_after_converge( const double& scf_ene_thr) { // only called if (GlobalC::exx_info.info_global.cal_exx) auto restart_reset = [this]() - { // avoid calling restart related procedure in the subsequent ion steps - GlobalC::restart.info_load.restart_exx = true; - this->exx_ptr->Eexx = 0; - }; + { // avoid calling restart related procedure in the subsequent ion steps + GlobalC::restart.info_load.restart_exx = true; + this->exx_ptr->Eexx = 0; + }; - // no separate_loop case - if (!GlobalC::exx_info.info_global.separate_loop) - { - GlobalC::exx_info.info_global.hybrid_step = 1; + if (!GlobalC::exx_info.info_global.separate_loop) // no separate_loop case + { + GlobalC::exx_info.info_global.hybrid_step = 1; - // in no_separate_loop case, scf loop only did twice - // in first scf loop, exx updated once in beginning, - // in second scf loop, exx updated every iter + // in no_separate_loop case, scf loop only did twice + // in first scf loop, exx updated once in beginning, + // in second scf loop, exx updated every iter - if (this->two_level_step) - { - restart_reset(); - return true; - } - else - { - // update exx and redo scf - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); - iter = 0; - std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; - this->two_level_step++; - return false; - } + if (this->two_level_step) + { + restart_reset(); + return true; } else - { // has separate_loop case - const double ediff = std::abs(etot - etot_last_outer_loop) * ModuleBase::Ry_to_eV; - if (two_level_step) { std::cout << FmtCore::format("EDIFF/eV (outer loop): %.8e \n", ediff); } - // exx converged or get max exx steps - if (this->two_level_step == GlobalC::exx_info.info_global.hybrid_step - || (iter == 1 && this->two_level_step != 0) // density convergence of outer loop - || (ediff < scf_ene_thr && this->two_level_step != 0)) //energy convergence of outer loop - { - restart_reset(); - return true; - } - else - { - this->etot_last_outer_loop = etot; - // update exx and redo scf - if (this->two_level_step == 0) - { - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); - } + { + // update exx and redo scf + XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + iter = 0; + std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; + this->two_level_step++; + return false; + } + } + else // has separate_loop case + { + const double ediff = std::abs(etot - etot_last_outer_loop) * ModuleBase::Ry_to_eV; + if (two_level_step) + { std::cout << FmtCore::format("EDIFF/eV (outer loop): %.8e \n", ediff); } + // exx converged or get max exx steps + if (this->two_level_step == GlobalC::exx_info.info_global.hybrid_step + || (iter == 1 && this->two_level_step != 0) // density convergence of outer loop + || (ediff < scf_ene_thr && this->two_level_step != 0)) //energy convergence of outer loop + { + restart_reset(); + return true; + } + else + { + this->etot_last_outer_loop = etot; + // update exx and redo scf + if (this->two_level_step == 0) + { XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); } - std::cout << " Updating EXX " << std::flush; - timeval t_start; gettimeofday(&t_start, nullptr); + std::cout << " Updating EXX " << std::flush; + timeval t_start; gettimeofday(&t_start, nullptr); - const bool flag_restart = (this->two_level_step == 0) ? true : false; + const bool flag_restart = (this->two_level_step == 0) ? 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);} - else - {this->mix_DMk_2D.mix(dm.get_DMK_vector(), 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(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), nspin) - : RI_2D_Comm::split_m2D_ktoR(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), nspin, this->exx_spacegroup_symmetry); + // 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(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), nspin) + : RI_2D_Comm::split_m2D_ktoR(*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), nspin, this->exx_spacegroup_symmetry); - // check the rotation of Ds - // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'D', Ds[0]); + // check the rotation of Ds + // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'D', Ds[0]); - // check the rotation of H(R) before adding exx - // this->symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, this->symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *lm.ParaV)); - // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', *(dynamic_cast*>(&hamilt)->getHR())); - // exit(0); + // check the rotation of H(R) before adding exx + // this->symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, this->symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *lm.ParaV)); + // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', *(dynamic_cast*>(&hamilt)->getHR())); + // exit(0); if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { @@ -236,31 +226,31 @@ bool Exx_LRI_Interface::exx_after_converge( // this->symrot_.print_HR(this->exx_ptr->Hexxs[0], "Hexxs_restore-DM-only"); // test // this->symrot_.print_HR(this->exx_ptr->Hexxs[0], "Hexxs_ref"); // test } - // ======================== test ======================== - // if (this->two_level_step)exit(0); - // check the rotation of S(R) - // this->symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, this->symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *lm.ParaV)); - // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', *(dynamic_cast*>(&hamilt)->getSR())); + // ======================== test ======================== + // if (this->two_level_step)exit(0); + // check the rotation of S(R) + // this->symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, this->symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *lm.ParaV)); + // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', *(dynamic_cast*>(&hamilt)->getSR())); - // check the rotation of D(R): no atom pair? - // symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *this->DM->get_paraV_pointer())); - // symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'D', *(this->DM->get_DMR_pointer(0))); + // check the rotation of D(R): no atom pair? + // symrot_.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, symrot_.get_Rs_from_adjacent_list(GlobalC::ucell, GlobalC::GridD, *this->DM->get_paraV_pointer())); + // symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'D', *(this->DM->get_DMR_pointer(0))); - // check the rotation of Hexx - // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', this->exx_ptr->Hexxs[0]); - // exit(0);// break after test - // ======================== test ======================== - iter = 0; - this->two_level_step++; + // check the rotation of Hexx + // this->symrot_.test_HR_rotation(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, 'H', this->exx_ptr->Hexxs[0]); + // exit(0);// break after test + // ======================== test ======================== + iter = 0; + this->two_level_step++; - timeval t_end; gettimeofday(&t_end, nullptr); - std::cout << "and rerun SCF\t" - << std::setprecision(3) << std::setiosflags(std::ios::scientific) - << (double)(t_end.tv_sec-t_start.tv_sec) + (double)(t_end.tv_usec-t_start.tv_usec)/1000000.0 - << std::defaultfloat << " (s)" << std::endl; + timeval t_end; gettimeofday(&t_end, nullptr); + std::cout << "and rerun SCF\t" + << std::setprecision(3) << std::setiosflags(std::ios::scientific) + << (double)(t_end.tv_sec-t_start.tv_sec) + (double)(t_end.tv_usec-t_start.tv_usec)/1000000.0 + << std::defaultfloat << " (s)" << std::endl; return false; - } } + } restart_reset(); return true; } diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index 76aab8a6e4..b51e9adc19 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -64,13 +64,8 @@ void Matrix_Orbs11::init_radial(const std::vectorMGT))); -} -} -} -} -} -} + Center2_Orb::Orb11(orb_A[TA][LA][NA], orb_B[TB][LB][NB], psb_, this->MGT))); + } } } } } } ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); } @@ -89,13 +84,8 @@ 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_, - this->MGT))); -} -} -} -} -} -} + this->MGT))); + } } } } } } ModuleBase::timer::tick("Matrix_Orbs11", "init_radial"); } @@ -109,13 +99,8 @@ void Matrix_Orbs11::init_radial_table() for (auto& coD: coC.second) { for (auto& coE: coD.second) { for (auto& coF: coE.second) { - coF.second.init_radial_table(); -} -} -} -} -} -} + coF.second.init_radial_table(); + } } } } }} ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); } @@ -123,7 +108,8 @@ void Matrix_Orbs11::init_radial_table(const std::map(position); - for (size_t i = 0; i != 4; ++i) { - radials.insert(iq + i); -} + for (size_t i = 0; i != 4; ++i) + { radials.insert(iq + i); } } for (auto& coC: *center2_orb11_sAB) { for (auto& coD: coC.second) { for (auto& coE: coD.second) { for (auto& coF: coE.second) { - coF.second.init_radial_table(radials); -} -} -} -} + coF.second.init_radial_table(radials); + } } } } } - } -} + } + } ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); } diff --git a/source/module_ri/Matrix_Orbs11.hpp b/source/module_ri/Matrix_Orbs11.hpp index d20aa330df..78b7cdd7f3 100644 --- a/source/module_ri/Matrix_Orbs11.hpp +++ b/source/module_ri/Matrix_Orbs11.hpp @@ -125,14 +125,14 @@ std::array,3> Matrix_Orbs11::cal_grad_overlap_matrix( } template -std::map>>>> Matrix_Orbs11::cal_overlap_matrix_all( - const ModuleBase::Element_Basis_Index::IndexLNM &index_r, +std::map>>>> Matrix_Orbs11::cal_overlap_matrix_all( + const ModuleBase::Element_Basis_Index::IndexLNM &index_r, const ModuleBase::Element_Basis_Index::IndexLNM &index_c ) const { ModuleBase::TITLE("Matrix_Orbs11","cal_overlap_matrix"); - + std::map>>>> matrixes; - + for( const auto &co1 : center2_orb11_s ) { const size_t TA = co1.first; diff --git a/source/module_ri/Matrix_Orbs22.cpp b/source/module_ri/Matrix_Orbs22.cpp index 027f252614..dd772f0348 100644 --- a/source/module_ri/Matrix_Orbs22.cpp +++ b/source/module_ri/Matrix_Orbs22.cpp @@ -48,7 +48,6 @@ void Matrix_Orbs22::init(const int mode, const LCAO_Orbitals& orb, const double this->MGT.init_Gaunt(2 * Lmax + 1); ModuleBase::timer::tick("Matrix_Orbs22", "init"); - std::cout << "Matrix_Orbs22::init()::done" << std::endl; } void Matrix_Orbs22::init_radial(const std::vector>>& orb_A1, @@ -113,7 +112,7 @@ void Matrix_Orbs22::init_radial(const LCAO_Orbitals& orb_A1, void Matrix_Orbs22::init_radial_table() { - ModuleBase::TITLE("Matrix_Orbs22", "init_radial"); + ModuleBase::TITLE("Matrix_Orbs22", "init_radial_table"); ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table"); for (auto& coA: center2_orb22_s) for (auto& coB: coA.second) @@ -132,7 +131,7 @@ void Matrix_Orbs22::init_radial_table() void Matrix_Orbs22::init_radial_table(const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs22", "init_radial_table_Rs"); - ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table"); + ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table_Rs"); for (const auto& RsA: Rs) for (const auto& RsB: RsA.second) { @@ -168,5 +167,5 @@ void Matrix_Orbs22::init_radial_table(const std::map Matrix_Orbs22::cal_overlap_matrix( case Matrix_Order::B2A2A1B1: m(iB2,iA2,iA1,iB1) = overlap; break; case Matrix_Order::B2A2B1A1: m(iB2,iA2,iB1,iA1) = overlap; break; case Matrix_Order::B2B1A1A2: m(iB2,iB1,iA1,iA2) = overlap; break; - case Matrix_Order::B2B1A2A1: m(iB2,iB1,iA2,iA1) = overlap; break; + case Matrix_Order::B2B1A2A1: m(iB2,iB1,iA2,iA1) = overlap; break; default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__)); } } @@ -219,7 +219,7 @@ std::array,3> Matrix_Orbs22::cal_grad_overlap_matrix( { const Tdata overlap = co10.second.cal_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB1, MB2 ); switch(matrix_order) - { + { const std::array grad_overlap = RI_Util::Vector3_to_array3(co10.second.cal_grad_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB1, MB2 )); const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; @@ -280,6 +280,7 @@ std::map < size_t, std::map>>>> matrixes; for( const auto &co1 : center2_orb22_s ) diff --git a/source/module_ri/exx_abfs-abfs_index.h b/source/module_ri/exx_abfs-abfs_index.h index 09cc493f2e..fbfb77af78 100644 --- a/source/module_ri/exx_abfs-abfs_index.h +++ b/source/module_ri/exx_abfs-abfs_index.h @@ -7,13 +7,15 @@ #include "../module_base/element_basis_index.h" #include "../module_basis/module_ao/ORB_atomic_lm.h" -class LCAO_Orbitals; + class LCAO_Orbitals; -class Exx_Abfs::Abfs_Index +namespace Exx_Abfs { -public: - static ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb ); - static ModuleBase::Element_Basis_Index::Range construct_range( const std::vector>> &orb ); -}; +namespace Abfs_Index +{ + extern ModuleBase::Element_Basis_Index::Range construct_range( const LCAO_Orbitals &orb ); + extern ModuleBase::Element_Basis_Index::Range construct_range( const std::vector>> &orb ); +} +} #endif // EXX_ABFS_ABFS_INDEX_H \ No newline at end of file diff --git a/source/module_ri/exx_abfs-construct_orbs.cpp b/source/module_ri/exx_abfs-construct_orbs.cpp index 2452f49f77..5407e69732 100644 --- a/source/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/module_ri/exx_abfs-construct_orbs.cpp @@ -45,16 +45,7 @@ std::vector>> Exx_Abfs::Construct_ } } - for (int T = 0; T < orbs.size() ; T++) - { - for (int L=orbs[T].size()-1; L >= 0 ; L--) - { - if (orbs[T][L].size()>0) - break; - else - orbs[T].resize(L); - } - } + Exx_Abfs::Construct_Orbs::filter_empty_orbs(orbs); return orbs; } @@ -478,6 +469,20 @@ inline const Numerical_Orbital_Lm &Exx_Abfs::Construct_Orbs::get_orbital( } */ +void Exx_Abfs::Construct_Orbs::filter_empty_orbs(std::vector>> &orbs) +{ + for (int T = 0; T < orbs.size() ; T++) + { + for (int L=orbs[T].size()-1; L >= 0 ; L--) + { + if (orbs[T][L].size()>0) + break; + else + orbs[T].resize(L); + } + } +} + void Exx_Abfs::Construct_Orbs::print_orbs_size( const std::vector>> &orbs, std::ostream &os) diff --git a/source/module_ri/exx_abfs-construct_orbs.h b/source/module_ri/exx_abfs-construct_orbs.h index be61fc7ff1..e33136696f 100644 --- a/source/module_ri/exx_abfs-construct_orbs.h +++ b/source/module_ri/exx_abfs-construct_orbs.h @@ -8,56 +8,60 @@ class LCAO_Orbitals; -class Exx_Abfs::Construct_Orbs +namespace Exx_Abfs { -public: - static std::vector>> change_orbs( +namespace Construct_Orbs +{ + extern std::vector>> change_orbs( const LCAO_Orbitals &orb_in, const double kmesh_times ); - static std::vector>> change_orbs( + extern std::vector>> change_orbs( const std::vector>> &orb_in, const double kmesh_times ); - static std::vector>> abfs_same_atom( + extern std::vector>> abfs_same_atom( const LCAO_Orbitals& orb, const std::vector>> &lcaos, const double kmesh_times_mot, const double times_threshold=0); - - static void print_orbs_size( + + extern void print_orbs_size( const std::vector>> &orbs, - std::ostream &os); - -private: - static std::vector>>> psi_mult_psi( + std::ostream &os); + + extern void filter_empty_orbs(std::vector>> &orbs); + + // private: + extern std::vector>>> psi_mult_psi( const std::vector>> &lcaos ); - - static std::vector>>> psir_mult_psir( + + extern std::vector>>> psir_mult_psir( const std::vector>> &lcaos ); - - static std::vector>>> orth( + + extern std::vector>>> orth( const std::vector>>> &psis, const std::vector>> &lcaos, const double norm_threshold = std::numeric_limits::min() ); - static std::vector>>> pca( + extern std::vector>>> pca( const LCAO_Orbitals& orb, const std::vector>> &abfs, const std::vector>> &orbs, const double kmesh_times_mot, const double times_threshold ); - - static std::vector>>> div_r( + + extern std::vector>>> div_r( const std::vector>>> &psirs, - const std::vector &r_radial ); - - static std::vector>> orbital( + const std::vector &r_radial ); + + extern std::vector>> orbital( const std::vector>>> &psis, const std::vector>> &orbs_info, - const double kmesh_times); - - static std::vector>>> get_psi( + const double kmesh_times); + + extern std::vector>>> get_psi( const std::vector>> &orbs ); -}; +} +} #endif // EXX_ABFS_IO_ASA_H diff --git a/source/module_ri/exx_abfs-io.h b/source/module_ri/exx_abfs-io.h index 15dd9303a3..9aed800e4b 100644 --- a/source/module_ri/exx_abfs-io.h +++ b/source/module_ri/exx_abfs-io.h @@ -13,30 +13,31 @@ #include "mpi.h" #endif -class LCAO_Orbitals; + class LCAO_Orbitals; -class Exx_Abfs::IO +namespace Exx_Abfs { -public: - - static std::vector>> construct_abfs( +namespace IO +{ + extern std::vector>> construct_abfs( const LCAO_Orbitals &orbs, const std::vector &files_abfs, - const double kmesh_times=1 ); // close dK, keep Kcut - - static std::vector>> construct_abfs( - const std::vector>> &abfs_pre, + const double kmesh_times=1 ); // close dK, keep Kcut + + extern std::vector>> construct_abfs( + const std::vector>> &abfs_pre, const LCAO_Orbitals &orbs, const std::vector &files_abfs, const double kmesh_times=1 ); // close dK, keep Kcut -private: - static std::vector> construct_abfs_T( + // private: + extern std::vector> construct_abfs_T( const std::string & file_name, const int &T, const int &nk, const double &dk, const double &dr_uniform); -}; +} +} #endif // EXX_ABFS_IO_H diff --git a/source/module_ri/exx_abfs-jle.cpp b/source/module_ri/exx_abfs-jle.cpp index 9a43c6632d..f62fd26006 100644 --- a/source/module_ri/exx_abfs-jle.cpp +++ b/source/module_ri/exx_abfs-jle.cpp @@ -7,27 +7,25 @@ #include "../module_base/mathzone.h" #include "../module_base/math_sphbes.h" // mohan add 2021-05-06 -bool Exx_Abfs::Jle::generate_matrix = false; -int Exx_Abfs::Jle::Lmax = 2; -double Exx_Abfs::Jle::Ecut_exx = 60; -double Exx_Abfs::Jle::tolerence = 1.0e-12; - -void Exx_Abfs::Jle::init_jle( const double kmesh_times, const LCAO_Orbitals& orb ) +std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> +Exx_Abfs::Jle::init_jle( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const double kmesh_times, + const LCAO_Orbitals& orb ) { - jle.resize( GlobalC::ucell.ntype ); - + std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> jle( GlobalC::ucell.ntype ); for (int T = 0; T < GlobalC::ucell.ntype ; T++) { - jle[T].resize( Lmax+1 ); - for (int L=0; L <= Lmax ; ++L) + jle[T].resize( info.abfs_Lmax+1 ); + for (int L=0; L<=info.abfs_Lmax; ++L) { const size_t ecut_number - = static_cast( sqrt( Ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. + = static_cast( sqrt( info.ecut_exx ) * orb.Phi[T].getRcut() / ModuleBase::PI ); // Rydberg Unit. jle[T][L].resize( ecut_number ); std::vector en(ecut_number, 0.0); - ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, tolerence, ModuleBase::GlobalFunc::VECTOR_TO_PTR(en), orb.Phi[T].getRcut()); + ModuleBase::Sphbes::Spherical_Bessel_Roots(ecut_number, L, info.tolerence, ModuleBase::GlobalFunc::VECTOR_TO_PTR(en), orb.Phi[T].getRcut()); for(size_t E=0; E!=ecut_number; ++E) { @@ -56,4 +54,5 @@ void Exx_Abfs::Jle::init_jle( const double kmesh_times, const LCAO_Orbitals& orb } } } + return jle; } diff --git a/source/module_ri/exx_abfs-jle.h b/source/module_ri/exx_abfs-jle.h index 967df31de6..e1b58c8123 100644 --- a/source/module_ri/exx_abfs-jle.h +++ b/source/module_ri/exx_abfs-jle.h @@ -2,25 +2,23 @@ #define EXX_ABFS_JLE_H #include "exx_abfs.h" - +#include "../module_hamilt_general/module_xc/exx_info.h" +#include "../module_basis/module_ao/ORB_atomic_lm.h" #include -#include "../module_basis/module_ao/ORB_read.h" - -class Exx_Abfs::Jle -{ -public: - - std::vector< - std::vector< - std::vector< - Numerical_Orbital_Lm>>> jle; - void init_jle( const double kmesh_times, const LCAO_Orbitals& orb ); + class LCAO_Orbitals; - static bool generate_matrix; - static int Lmax; - static double Ecut_exx; - static double tolerence; -}; +namespace Exx_Abfs +{ +namespace Jle +{ + // jle[T][L][E] + extern std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> + init_jle( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const double kmesh_times, + const LCAO_Orbitals &orb ); +} +} #endif // EXX_ABFS_JLE_H diff --git a/source/module_ri/exx_abfs.h b/source/module_ri/exx_abfs.h index 845eb3d81b..203a412b5e 100644 --- a/source/module_ri/exx_abfs.h +++ b/source/module_ri/exx_abfs.h @@ -1,29 +1,18 @@ #ifndef EXX_ABFS_H #define EXX_ABFS_H -#include -using std::vector; -#include -using std::map; -#include - #include "../module_basis/module_ao/ORB_atomic_lm.h" #include "../module_base/element_basis_index.h" #include "../module_base/matrix.h" #include "../module_base/vector3.h" -class Exx_Abfs +namespace Exx_Abfs { -public: - class Abfs_Index; - class Jle; - class IO; - class Construct_Orbs; - class PCA; - - int rmesh_times = 5; // Peize Lin test - int kmesh_times = 1; // Peize Lin test - -}; + namespace Abfs_Index{} + namespace Jle{} + namespace IO{} + namespace Construct_Orbs{} + namespace PCA{} +} #endif diff --git a/source/module_ri/exx_opt_orb-print.cpp b/source/module_ri/exx_opt_orb-print.cpp index 52505505d9..af38191c83 100644 --- a/source/module_ri/exx_opt_orb-print.cpp +++ b/source/module_ri/exx_opt_orb-print.cpp @@ -3,15 +3,16 @@ #include "exx_abfs-jle.h" void Exx_Opt_Orb::print_matrix( - const K_Vectors &kv, - const std::string& file_name, - const std::vector> &matrix_Q, + const Exx_Info::Exx_Info_Opt_ABFs &info, + const K_Vectors &kv, + const std::string& file_name, + const std::vector> &matrix_Q, const std::vector>> &matrix_S, const RI::Tensor &matrix_V, const size_t TA, const size_t IA, const size_t TB, const size_t IB, - const std::vector& orb_cutoff, - const ModuleBase::Element_Basis_Index::Range &range_jles, - const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const + const std::vector& orb_cutoff, + const ModuleBase::Element_Basis_Index::Range &range_jles, + const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) { auto print_header = [&]( std::ofstream &ofs ) { @@ -20,7 +21,7 @@ void Exx_Opt_Orb::print_matrix( ofs << GlobalC::ucell.latvec.e11 << " " << GlobalC::ucell.latvec.e12 << " " << GlobalC::ucell.latvec.e13 << std::endl; ofs << GlobalC::ucell.latvec.e21 << " " << GlobalC::ucell.latvec.e22 << " " << GlobalC::ucell.latvec.e23 << std::endl; ofs << GlobalC::ucell.latvec.e31 << " " << GlobalC::ucell.latvec.e32 << " " << GlobalC::ucell.latvec.e33 << std::endl; - + if( TA==TB ) { ofs << 1 << " ntype" << std::endl; @@ -28,18 +29,18 @@ void Exx_Opt_Orb::print_matrix( if( IA==IB ) { ofs << 1 << " na" << std::endl; - ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " - << GlobalC::ucell.atoms[TA].tau[IA].y << " " + ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " + << GlobalC::ucell.atoms[TA].tau[IA].y << " " << GlobalC::ucell.atoms[TA].tau[IA].z << std::endl; } else { ofs << 2 << " na" << std::endl; - ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " + ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " << GlobalC::ucell.atoms[TA].tau[IA].y << " " << GlobalC::ucell.atoms[TA].tau[IA].z << std::endl; - ofs << GlobalC::ucell.atoms[TB].tau[IB].x << " " - << GlobalC::ucell.atoms[TB].tau[IB].y << " " + ofs << GlobalC::ucell.atoms[TB].tau[IB].x << " " + << GlobalC::ucell.atoms[TB].tau[IB].y << " " << GlobalC::ucell.atoms[TB].tau[IB].z << std::endl; } } @@ -48,58 +49,58 @@ void Exx_Opt_Orb::print_matrix( ofs << 2 << " ntype" << std::endl; ofs << GlobalC::ucell.atoms[TA].label << " label" << std::endl; ofs << 1 << " na" << std::endl; - ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " - << GlobalC::ucell.atoms[TA].tau[IA].y << " " + ofs << GlobalC::ucell.atoms[TA].tau[IA].x << " " + << GlobalC::ucell.atoms[TA].tau[IA].y << " " << GlobalC::ucell.atoms[TA].tau[IA].z << std::endl; ofs << GlobalC::ucell.atoms[TB].label << " label" << std::endl; ofs << 1 << " na" << std::endl; - ofs << GlobalC::ucell.atoms[TB].tau[IB].x << " " - << GlobalC::ucell.atoms[TB].tau[IB].y << " " + ofs << GlobalC::ucell.atoms[TB].tau[IB].x << " " + << GlobalC::ucell.atoms[TB].tau[IB].y << " " << GlobalC::ucell.atoms[TB].tau[IB].z << std::endl; } - + // ecutwfc_jlq determine the jlq corresponding to plane wave calculation. - ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc" << std::endl; // mohan add 2009-09-08 + ofs << info.ecut_exx << " ecutwfc" << std::endl; // mohan add 2009-09-08 // this parameter determine the total number of jlq. - ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 + ofs << info.ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 if(TA==TB) - ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; + { ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; } else - ofs << orb_cutoff[TA] << " " << orb_cutoff[TB] << " rcut_Jlq" << std::endl; + { ofs << orb_cutoff[TA] << " " << orb_cutoff[TB] << " rcut_Jlq" << std::endl; } // mohan add 'smooth' and 'smearing_sigma' 2009-08-28 ofs << 0 << " smooth" << std::endl; ofs << 0 << " smearing_sigma" << std::endl; - ofs << Exx_Abfs::Jle::tolerence << " tolerence" << std::endl; + ofs << info.tolerence << " tolerence" << std::endl; - ofs << Exx_Abfs::Jle::Lmax << " lmax" << std::endl; + ofs << info.abfs_Lmax << " lmax" << std::endl; ofs << kv.get_nkstot() << " nks" << std::endl; - assert( matrix_V.shape[0] == matrix_V.shape[1] ); - ofs << matrix_V.shape[0] << " nbands" << std::endl; - + assert( matrix_V.shape[0]*matrix_V.shape[1] == matrix_V.shape[2]*matrix_V.shape[3] ); + ofs << matrix_V.shape[0]*matrix_V.shape[1] << " nbands" << std::endl; + auto cal_sum_M = [&range_jles](size_t T) -> size_t { size_t sum_M = 0; for( size_t L = 0; L!=range_jles[T].size(); ++L ) - sum_M += range_jles[T][L].M; + { sum_M += range_jles[T][L].M; } return sum_M; }; const size_t nwfc = (TA==TB && IA==IB) ? cal_sum_M(TA) : cal_sum_M(TA)+cal_sum_M(TB); ofs << nwfc << " nwfc" << std::endl; - - const size_t ecut_numberA = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * orb_cutoff[TA] / ModuleBase::PI ); // Rydberg Unit - const size_t ecut_numberB = static_cast( sqrt( Exx_Abfs::Jle::Ecut_exx ) * orb_cutoff[TB] / ModuleBase::PI ); // Rydberg Unit + + const size_t ecut_numberA = static_cast( sqrt( info.ecut_exx ) * orb_cutoff[TA] / ModuleBase::PI ); // Rydberg Unit + const size_t ecut_numberB = static_cast( sqrt( info.ecut_exx ) * orb_cutoff[TB] / ModuleBase::PI ); // Rydberg Unit if(TA==TB) - ofs << ecut_numberA << " ne" << std::endl; + { ofs << ecut_numberA << " ne" << std::endl; } else - ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; - + { ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; } + ofs << "" << std::endl; - for( int ik=0; ik!=kv.get_nkstot(); ++ik ) + for( int ik=0; ik!=kv.get_nkstot(); ++ik ) { ofs << kv.kvec_c[ik].x << " " << kv.kvec_c[ik].y << " " << kv.kvec_c[ik].z; ofs << " " << kv.wk[ik] * 0.5 << std::endl; @@ -108,44 +109,47 @@ void Exx_Opt_Orb::print_matrix( ofs << std::endl; }; - - + + auto print_Q = [&]( std::ofstream &ofs ) { //--------------------- // < Psi | jY > //--------------------- ofs<< "" << std::endl; - - for( size_t ib=0; ib!=matrix_V.shape[0]; ++ib ) + + for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 ) { - for( size_t iat=0; iat!=matrix_Q.size(); ++iat ) + for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 ) { - const size_t it = (iat==0) ? TA : TB; - for( size_t il=0; il!=range_jles[it].size(); ++il ) + for( size_t iat=0; iat!=matrix_Q.size(); ++iat ) { - for( size_t im=0; im!=range_jles[it][il].M; ++im ) + const size_t it = (iat==0) ? TA : TB; + for( size_t il=0; il!=range_jles[it].size(); ++il ) { - for( size_t iq=0; iq!=range_jles[it][il].N; ++iq ) + for( size_t im=0; im!=range_jles[it][il].M; ++im ) { - ofs<" << std::endl << std::endl; }; - - + + auto print_S = [&]( std::ofstream &ofs, const double scale=1 ) { //--------------------- // < jY | jY > //--------------------- ofs<< "" <" << std::endl << std::endl; }; - - + + auto print_V = [&]( std::ofstream &ofs, const double scale=1 ) { //--------------------- // < Psi | Psi > - //--------------------- + //--------------------- ofs << "" << std::endl; - - for( size_t ib1=0; ib1!=matrix_V.shape[0]; ++ib1 ) + + for( size_t iw0=0; iw0!=matrix_V.shape[0]; ++iw0 ) { - for( size_t ib2=0; ib2!=matrix_V.shape[1]; ++ib2 ) + for( size_t iw1=0; iw1!=matrix_V.shape[1]; ++iw1 ) { - ofs<" << std::endl << std::endl; }; - - std::ofstream ofs(file_name+"_"+ModuleBase::GlobalFunc::TO_STRING(TA)+"_"+ModuleBase::GlobalFunc::TO_STRING(IA)+"_"+ModuleBase::GlobalFunc::TO_STRING(TB)+"_"+ModuleBase::GlobalFunc::TO_STRING(IB)); + + std::ofstream ofs(file_name+"_"+std::to_string(TA)+"_"+std::to_string(IA)+"_"+std::to_string(TB)+"_"+std::to_string(IB)); print_header(ofs); print_Q(ofs); print_S(ofs); diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index d714de1a9d..bb3b60321d 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -10,35 +10,30 @@ #include "module_ri/Matrix_Orbs21.h" #include "module_ri/Matrix_Orbs22.h" #include "module_ri/LRI_CV_Tools.h" +#include -#include "../module_ri/test_code/element_basis_index-test.h" -#include "../module_ri/test_code/test_function.h" -#include - -void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const +void Exx_Opt_Orb::generate_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const K_Vectors &kv, + const LCAO_Orbitals &orb) { -// std::ofstream ofs_mpi(GlobalC::exx_lcao.test_dir.process+"time_"+ModuleBase::GlobalFunc::TO_STRING(GlobalV::MY_RANK),std::ofstream::app); - ModuleBase::TITLE("Exx_Opt_Orb::generate_matrix"); -// ofs_mpi<<"memory:\t"<>> - lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->kmesh_times ); + std::vector>> + lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, info.kmesh_times ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(lcaos); - const std::vector>> - abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + std::vector>> + abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, lcaos, info.kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(abfs); -// ofs_mpi<<"memory:\t"<kmesh_times, orb ); + std::vector< std::vector< std::vector< Numerical_Orbital_Lm>>> + jle = Exx_Abfs::Jle::init_jle( info, info.kmesh_times, orb ); + Exx_Abfs::Construct_Orbs::filter_empty_orbs(jle); -// ofs_mpi<<"memory:\t"<(abfs[T].size())-1 ); -} + GlobalC::exx_info.info_ri.abfs_Lmax = info.abfs_Lmax; + for( 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 ); } const ModuleBase::Element_Basis_Index::Range range_lcaos = Exx_Abfs::Abfs_Index::construct_range( lcaos ); const ModuleBase::Element_Basis_Index::IndexLNM index_lcaos = ModuleBase::Element_Basis_Index::construct_index( range_lcaos ); @@ -46,34 +41,18 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const ModuleBase::Element_Basis_Index::Range range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs ); const ModuleBase::Element_Basis_Index::IndexLNM index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs ); - const ModuleBase::Element_Basis_Index::Range range_jys = Exx_Abfs::Abfs_Index::construct_range( jle.jle ); + const ModuleBase::Element_Basis_Index::Range range_jys = Exx_Abfs::Abfs_Index::construct_range( jle ); const ModuleBase::Element_Basis_Index::IndexLNM index_jys = ModuleBase::Element_Basis_Index::construct_index( range_jys ); -// ofs_mpi<>> radial_R = get_radial_R(); -#if TEST_EXX_RADIAL==2 - { - for(const auto & rA : radial_R) - for(const auto & rB : rA.second) - { - ofs_mpi< - const auto ms_lcaoslcaos_lcaoslcaos = [&]() -> std::map>>>> + const auto ms_lcaoslcaos_lcaoslcaos = [&]() -> std::map>>>> { Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos; - m_lcaoslcaos_lcaoslcaos.init( 1, orb, this->kmesh_times, 1 ); + m_lcaoslcaos_lcaoslcaos.init( 1, orb, info.kmesh_times, 1 ); m_lcaoslcaos_lcaoslcaos.init_radial( lcaos, lcaos, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_lcaoslcaos_lcaoslcaos.init_radial_table(radial_R); @@ -82,15 +61,14 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) #endif return m_lcaoslcaos_lcaoslcaos.cal_overlap_matrix_all( index_lcaos, index_lcaos, index_lcaos, index_lcaos); }(); - -// ofs_mpi<<"memory:\t"< const auto ms_lcaoslcaos_jys = [&]() -> std::map>>>>> { Matrix_Orbs21 m_jyslcaos_lcaos; - m_jyslcaos_lcaos.init( 1, orb, this->kmesh_times, 1 ); - m_jyslcaos_lcaos.init_radial( jle.jle, lcaos, lcaos ); + m_jyslcaos_lcaos.init( 1, orb, info.kmesh_times, 1 ); + m_jyslcaos_lcaos.init_radial( jle, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_jyslcaos_lcaos.init_radial_table(radial_R); #else @@ -99,14 +77,12 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) return m_jyslcaos_lcaos.cal_overlap_matrix_all( index_jys, index_lcaos, index_lcaos ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_jys_jys = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_jys; - m_jys_jys.init( 2, orb, this->kmesh_times, 1 ); - m_jys_jys.init_radial( jle.jle, jle.jle ); + m_jys_jys.init( 2, orb, info.kmesh_times, 1 ); + m_jys_jys.init_radial( jle, jle ); #if TEST_EXX_RADIAL>=1 m_jys_jys.init_radial_table(radial_R); #else @@ -115,13 +91,11 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) return m_jys_jys.cal_overlap_matrix_all( index_jys, index_jys ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_abfs_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_abfs_abfs; - m_abfs_abfs.init( 2, orb, this->kmesh_times, 1 ); + m_abfs_abfs.init( 2, orb, info.kmesh_times, 1 ); m_abfs_abfs.init_radial( abfs, abfs ); #if TEST_EXX_RADIAL>=1 m_abfs_abfs.init_radial_table(radial_R); @@ -131,13 +105,11 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) return m_abfs_abfs.cal_overlap_matrix_all( index_abfs, index_abfs ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_lcaoslcaos_abfs = [&]() -> std::map>>>>> { Matrix_Orbs21 m_abfslcaos_lcaos; - m_abfslcaos_lcaos.init( 1, orb, this->kmesh_times, 1 ); + m_abfslcaos_lcaos.init( 1, orb, info.kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 m_abfslcaos_lcaos.init_radial_table(radial_R); @@ -147,14 +119,12 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) return m_abfslcaos_lcaos.cal_overlap_matrix_all( index_abfs, index_lcaos, index_lcaos ); }(); -// ofs_mpi<<"memory:\t"< const auto ms_jys_abfs = [&]() -> std::map>>>> { Matrix_Orbs11 m_jys_abfs; - m_jys_abfs.init( 2, orb, this->kmesh_times, 1 ); - m_jys_abfs.init_radial( jle.jle, abfs ); + m_jys_abfs.init( 2, orb, info.kmesh_times, 1 ); + m_jys_abfs.init_radial( jle, abfs ); #if TEST_EXX_RADIAL>=1 m_jys_abfs.init_radial_table(radial_R); #else @@ -163,15 +133,6 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) return m_jys_abfs.cal_overlap_matrix_all( index_jys, index_abfs ); }(); -// ofs_mpi<<"memory:\t"<>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, T,I,T,I ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj( + cal_proj_22( ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj( + {cal_proj_21( ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I)[0], ms_lcaoslcaos_abfs.at(T).at(I).at(T).at(I), ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj( + {{cal_proj_11( ms_jys_jys.at(T).at(I).at(T).at(I), {ms_jys_abfs.at(T).at(I).at(T).at(I)}, ms_abfs_abfs_I, {ms_jys_abfs.at(T).at(I).at(T).at(I)})}}; - print_matrix(kv, - "matrix", + print_matrix( + GlobalC::exx_info.info_opt_abfs, + kv, + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", m_lcaoslcaos_jys_proj, m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, TA, IA, TB, IB, - orb.cutoffs(), + orb.cutoffs(), range_jys, index_jys ); } else { - print_matrix(kv, - "matrix", + print_matrix( + GlobalC::exx_info.info_opt_abfs, + kv, + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I), {{ms_jys_jys.at(T).at(I).at(T).at(I)}}, ms_lcaoslcaos_lcaoslcaos.at(T).at(I).at(T).at(I), TA, IA, TB, IB, - orb.cutoffs(), + orb.cutoffs(), range_jys, index_jys ); } } @@ -237,64 +202,68 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const std::vector>> ms_abfs_abfs_I = cal_I( ms_abfs_abfs, TA,IA,TB,IB ); // < lcaos lcaos | lcaos lcaos > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | lcaos lcaos > const RI::Tensor m_lcaoslcaos_lcaoslcaos_proj = - cal_proj( + cal_proj_22( ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB), ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB)); // < lcaos lcaos | jys > - < lcaos lcaos | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector> m_lcaoslcaos_jys_proj = - {cal_proj( + {cal_proj_21( ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[0], ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_21( ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB)[1], ms_lcaoslcaos_abfs.at(TA).at(IA).at(TB).at(IB), ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) })}; // < jys | jys > - < jys | abfs > * < abfs | abfs >.I * < abfs | jys > const std::vector>> m_jys_jys_proj = - {{cal_proj( + {{cal_proj_11( ms_jys_jys.at(TA).at(IA).at(TA).at(IA), { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_11( ms_jys_jys.at(TA).at(IA).at(TB).at(IB), { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }, - {cal_proj( + {cal_proj_11( ms_jys_jys.at(TB).at(IB).at(TA).at(IA), { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TA).at(IA).at(TA).at(IA), ms_jys_abfs.at(TA).at(IA).at(TB).at(IB) }), - cal_proj( + cal_proj_11( ms_jys_jys.at(TB).at(IB).at(TB).at(IB), { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }, ms_abfs_abfs_I, { ms_jys_abfs.at(TB).at(IB).at(TA).at(IA), ms_jys_abfs.at(TB).at(IB).at(TB).at(IB) }) }}; - print_matrix(kv, - "matrix", + print_matrix( + GlobalC::exx_info.info_opt_abfs, + kv, + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", m_lcaoslcaos_jys_proj, m_jys_jys_proj, m_lcaoslcaos_lcaoslcaos_proj, TA, IA, TB, IB, - orb.cutoffs(), + orb.cutoffs(), range_jys, index_jys ); } else { - print_matrix(kv, - "matrix", + print_matrix( + GlobalC::exx_info.info_opt_abfs, + kv, + PARAM.globalv.global_out_dir+"/matrix-opt-abfs", ms_lcaoslcaos_jys.at(TA).at(IA).at(TB).at(IB), {{ms_jys_jys.at(TA).at(IA).at(TA).at(IA), ms_jys_jys.at(TA).at(IA).at(TB).at(IB)}, {ms_jys_jys.at(TB).at(IB).at(TA).at(IA), ms_jys_jys.at(TB).at(IB).at(TB).at(IB)}}, ms_lcaoslcaos_lcaoslcaos.at(TA).at(IA).at(TB).at(IB), TA, IA, TB, IB, - orb.cutoffs(), + orb.cutoffs(), range_jys, index_jys ); } } @@ -305,27 +274,82 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) } // m_big - m_left * m_middle * m_right.T -RI::Tensor Exx_Opt_Orb::cal_proj( - const RI::Tensor & m_big, - const std::vector> & m_left, - const std::vector>> & m_middle, - const std::vector> & m_right ) const +RI::Tensor Exx_Opt_Orb::cal_proj_22( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) { - ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); - -//auto print_nrc = [](const matrix & m){ std::cout<<"\t"< m_proj = m_big; -//print_nrc(m_proj); + ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0y1_x0x1a_y0y1a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +} +RI::Tensor Exx_Opt_Orb::cal_proj_21( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) +{ + ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); + RI::Tensor m_proj = m_big.copy(); for( size_t il=0; il!=m_left.size(); ++il ) { for( size_t ir=0; ir!=m_right.size(); ++ir ) { -//std::cout< m_lm = RI::Tensor_Multiply::x0x1y1_x0x1a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0x1y0_x0x1a_y0a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +} +/*RI::Tensor Exx_Opt_Orb::cal_proj_12( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) +{ + ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0y1_x0a_y0y1a(m_lm, m_right[ir]); + m_proj -= m_lmr; + } + } + return m_proj; +}*/ +RI::Tensor Exx_Opt_Orb::cal_proj_11( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ) +{ + ModuleBase::TITLE("Exx_Opt_Orb::cal_proj"); + RI::Tensor m_proj = m_big.copy(); + for( size_t il=0; il!=m_left.size(); ++il ) + { + for( size_t ir=0; ir!=m_right.size(); ++ir ) + { + // m_proj = m_proj - m_left[il] * m_middle[il][ir] * m_right[ir].T; + const RI::Tensor m_lm = RI::Tensor_Multiply::x0y1_x0a_ay1(m_left[il], m_middle[il][ir]); + const RI::Tensor m_lmr = RI::Tensor_Multiply::x0y0_x0a_y0a(m_lm, m_right[ir]); + m_proj -= m_lmr; } } return m_proj; @@ -333,34 +357,35 @@ RI::Tensor Exx_Opt_Orb::cal_proj( std::vector>> Exx_Opt_Orb::cal_I( const std::map>>>> &ms, - const size_t TA, const size_t IA, const size_t TB, const size_t IB ) const + const size_t TA, const size_t IA, const size_t TB, const size_t IB ) { ModuleBase::TITLE("Exx_Opt_Orb::cal_I"); if( TA==TB && IA==IB ) { - std::vector>> m_I - {{ RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA).shape) }}; - return LRI_CV_Tools::cal_I(m_I); + return {{LRI_CV_Tools::cal_I(RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA)))}}; } else { - std::vector>> m_I - {{ RI::Tensor(ms.at(TA).at(IA).at(TA).at(IA).shape), - RI::Tensor(ms.at(TA).at(IA).at(TB).at(IB).shape) }, - { RI::Tensor(ms.at(TB).at(IB).at(TA).at(IA).shape), - RI::Tensor(ms.at(TB).at(IB).at(TB).at(IB).shape) }}; - return LRI_CV_Tools::cal_I(m_I); + std::vector>> m_in + {{ ms.at(TA).at(IA).at(TA).at(IA), + ms.at(TA).at(IA).at(TB).at(IB) }, + { ms.at(TB).at(IB).at(TA).at(IA), + ms.at(TB).at(IB).at(TB).at(IB) }}; + return LRI_CV_Tools::cal_I(m_in); } } -std::map>> Exx_Opt_Orb::get_radial_R() const +std::map>> Exx_Opt_Orb::get_radial_R() { ModuleBase::TITLE("Exx_Opt_Orb::get_radial_R"); std::map>> radial_R; - for( size_t TA=0; TA!=GlobalC::ucell.ntype; ++TA ) { - for( size_t IA=0; IA!=GlobalC::ucell.atoms[TA].na; ++IA ) { - for( size_t TB=0; TB!=GlobalC::ucell.ntype; ++TB ) { + for( size_t TA=0; TA!=GlobalC::ucell.ntype; ++TA ) + { + for( size_t IA=0; IA!=GlobalC::ucell.atoms[TA].na; ++IA ) + { + for( size_t TB=0; TB!=GlobalC::ucell.ntype; ++TB ) + { for( size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) { const ModuleBase::Vector3 &tauA = GlobalC::ucell.atoms[TA].tau[IA]; @@ -369,8 +394,8 @@ std::map>> Exx_Opt_Orb::get_radial_R() c radial_R[TA][TB].insert( delta_R ); radial_R[TB][TA].insert( delta_R ); } -} -} -} + } + } + } return radial_R; } diff --git a/source/module_ri/exx_opt_orb.h b/source/module_ri/exx_opt_orb.h index ba6910f2ba..157bc016ee 100644 --- a/source/module_ri/exx_opt_orb.h +++ b/source/module_ri/exx_opt_orb.h @@ -1,40 +1,60 @@ #ifndef EXX_OPT_ORB_H #define EXX_OPT_ORB_H +#include "../module_hamilt_general/module_xc/exx_info.h" #include "../module_base/matrix.h" #include "../module_base/element_basis_index.h" #include "module_cell/klist.h" #include "module_basis/module_ao/ORB_read.h" +#include #include #include #include -#include -class Exx_Opt_Orb +namespace Exx_Opt_Orb { -public: - void generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) const; -private: - std::vector>> cal_I( - const std::map>>>> &ms, - const size_t TA, const size_t IA, const size_t TB, const size_t IB ) const; - RI::Tensor cal_proj( - const RI::Tensor & m_big, - const std::vector> & m_left, - const std::vector>> & m_middle, - const std::vector> & m_right ) const; - void print_matrix( + extern void generate_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, + const K_Vectors &kv, + const LCAO_Orbitals& orb); + + // private: + extern std::vector>> cal_I( + const std::map>>>> &ms, + const size_t TA, const size_t IA, const size_t TB, const size_t IB ); + extern RI::Tensor cal_proj_22( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ); + extern RI::Tensor cal_proj_21( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ); + /* + extern RI::Tensor cal_proj_12( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ); + */ + extern RI::Tensor cal_proj_11( + const RI::Tensor & m_big, + const std::vector> & m_left, + const std::vector>> & m_middle, + const std::vector> & m_right ); + extern void print_matrix( + const Exx_Info::Exx_Info_Opt_ABFs &info, const K_Vectors &kv, const std::string& file_name, - const std::vector> &matrix_Q, + const std::vector> &matrix_Q, const std::vector>> &matrix_S, const RI::Tensor &matrix_V, const size_t TA, const size_t IA, const size_t TB, const size_t IB, const std::vector& orb_cutoff, - const ModuleBase::Element_Basis_Index::Range &range_jles, - const ModuleBase::Element_Basis_Index::IndexLNM &index_jles) const; - std::map>> get_radial_R() const; - - int kmesh_times = 4; + const ModuleBase::Element_Basis_Index::Range &range_jles, + const ModuleBase::Element_Basis_Index::IndexLNM &index_jles); + extern std::map>> get_radial_R(); }; #endif