diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 58a900fe11..60fbc6fae9 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -31,39 +31,38 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c ModuleBase::TITLE("Exx_LRI","init"); ModuleBase::timer::tick("Exx_LRI", "init"); -// if(GlobalC::exx_info.info_global.separate_loop) -// { -// Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::No; -// Hexx_para.mixing_beta = 0; -// } -// else -// { -// if("plain"==GlobalC::CHR.mixing_mode) -// Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Plain; -// else if("pulay"==GlobalC::CHR.mixing_mode) -// Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Pulay; -// else -// throw std::invalid_argument("exx mixing error. exx_separate_loop==false, mixing_mode!=plain or pulay"); -// Hexx_para.mixing_beta = GlobalC::CHR.mixing_beta; -// } - - this->mpi_comm = mpi_comm_in; - this->p_kv = &kv_in; - this->orb_cutoff_ = orb.cutoffs(); + // if(GlobalC::exx_info.info_global.separate_loop) + // { + // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::No; + // Hexx_para.mixing_beta = 0; + // } + // else + // { + // if("plain"==GlobalC::CHR.mixing_mode) + // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Plain; + // else if("pulay"==GlobalC::CHR.mixing_mode) + // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Pulay; + // else + // throw std::invalid_argument("exx mixing error. exx_separate_loop==false, mixing_mode!=plain or pulay"); + // 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->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times ); -// #ifdef __MPI -// Exx_Abfs::Util::bcast( this->info.files_abfs, 0, this->mpi_comm ); -// #endif + // #ifdef __MPI + // Exx_Abfs::Util::bcast( this->info.files_abfs, 0, this->mpi_comm ); + // #endif 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,15 +84,14 @@ 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 ); @@ -106,19 +104,17 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) 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 ); + // 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); + // this->m_abfsabfs.init_radial_table(Rradial); + // 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), @@ -137,7 +133,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) @@ -166,7 +163,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) @@ -185,44 +183,48 @@ void Exx_LRI::cal_exx_ions(const bool write_cv) 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"); const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(pv); + if(p_symrot) + { this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()); } + else + { this->exx_lri.set_symmetry(false, {}); } + 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); + const 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); @@ -245,8 +247,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; } @@ -280,8 +282,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? @@ -304,8 +305,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? @@ -318,16 +318,15 @@ void 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/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp index 8640088130..f7dea18e90 100644 --- a/source/module_ri/LRI_CV.hpp +++ b/source/module_ri/LRI_CV.hpp @@ -37,7 +37,7 @@ LRI_CV::~LRI_CV() template void LRI_CV::set_orbitals( - const LCAO_Orbitals& orb, + const LCAO_Orbitals& orb, const std::vector>> &lcaos_in, const std::vector>> &abfs_in, const std::vector>> &abfs_ccp_in, @@ -47,7 +47,7 @@ void LRI_CV::set_orbitals( ModuleBase::TITLE("LRI_CV", "set_orbitals"); ModuleBase::timer::tick("LRI_CV", "set_orbitals"); - this->orb_cutoff_ = orb.cutoffs(); + this->orb_cutoff_ = orb.cutoffs(); this->lcaos = lcaos_in; this->abfs = abfs_in; this->abfs_ccp = abfs_ccp_in; @@ -109,11 +109,11 @@ auto LRI_CV::cal_datas( if( R_delta.norm()*GlobalC::ucell.lat0 < Rcut ) { const Tresult Data = func_DPcal_data(it0, it1, R_delta, flags); -// if(Data.norm(std::numeric_limits::max()) > threshold) -// { + // if(Data.norm(std::numeric_limits::max()) > threshold) + // { #pragma omp critical(LRI_CV_cal_datas) Datas[list_A0[i0]][list_A1[i1]] = Data; -// } + // } } } }