From de2a5563b348293474ea5abe40b8e068742af3e9 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 12:12:09 +0800 Subject: [PATCH 01/24] change ucell in module_exx/exx_lri_interface.hpp --- source/module_ri/Exx_LRI_interface.h | 30 ++++++++++------- source/module_ri/Exx_LRI_interface.hpp | 46 ++++++++++++++++---------- 2 files changed, 47 insertions(+), 29 deletions(-) diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index b3d358a925..42172aada1 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -64,19 +64,25 @@ class Exx_LRI_Interface void exx_hamilt2density(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter); /// @brief in iter_finish: write Hexx, do something according to whether SCF is converged - void exx_iter_finish(const K_Vectors& kv, const UnitCell& ucell, - hamilt::Hamilt& hamilt, elecstate::ElecState& elec, Charge_Mixing& chgmix, - const double& scf_ene_thr, int& iter, const int istep, bool& conv_esolver); + void exx_iter_finish(const K_Vectors& kv, + const UnitCell& ucell, + hamilt::Hamilt& hamilt, + elecstate::ElecState& elec, + Charge_Mixing& chgmix, + const double& scf_ene_thr, + int& iter, + const int istep, + bool& conv_esolver); /// @brief: in do_after_converge: add exx operators; do DM mixing if seperate loop - bool exx_after_converge( - hamilt::Hamilt& hamilt, - const elecstate::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, - const K_Vectors& kv, - const int& nspin, - int& iter, - const int& istep, - const double& etot, - const double& scf_ene_thr); + bool exx_after_converge(const UnitCell& ucell, + hamilt::Hamilt& hamilt, + const elecstate::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, + const K_Vectors& kv, + const int& nspin, + int& iter, + const int& istep, + const double& etot, + const double& scf_ene_thr); int two_level_step = 0; double etot_last_outer_loop = 0.0; elecstate::DensityMatrix* dm_last_step; diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 5eee38d61d..3800ca0419 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -53,14 +53,18 @@ void Exx_LRI_Interface::exx_before_all_runners(const K_Vectors& kv, co } template -void Exx_LRI_Interface::exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb) +void Exx_LRI_Interface::exx_beforescf(const int istep, + const K_Vectors& kv, + const Charge_Mixing& chgmix, + const UnitCell& ucell, + const LCAO_Orbitals& orb) { #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); + if (GlobalC::restart.info_load.load_H_finish && !GlobalC::restart.info_load.restart_exx) { XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } - else if (istep > 0) { XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + else if (istep > 0) { XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } else { @@ -168,9 +172,15 @@ void Exx_LRI_Interface::exx_hamilt2density(elecstate::ElecState& elec, } template -void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, const UnitCell& ucell, - hamilt::Hamilt& hamilt, elecstate::ElecState& elec, Charge_Mixing& chgmix, - const double& scf_ene_thr, int& iter, const int istep, bool& conv_esolver) +void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, + const UnitCell& ucell, + hamilt::Hamilt& hamilt, + elecstate::ElecState& elec, + Charge_Mixing& chgmix, + const double& scf_ene_thr, + int& iter, + const int istep, + bool& conv_esolver) { if (GlobalC::restart.info_save.save_H && (this->two_level_step > 0 || istep > 0) && (!GlobalC::exx_info.info_global.separate_loop || iter == 1)) // to avoid saving the same value repeatedly @@ -194,7 +204,7 @@ void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, const Uni }*/ ////////// for Add_Hexx_Type:R const std::string& restart_HR_path = GlobalC::restart.folder + "HexxR" + std::to_string(GlobalV::MY_RANK); - ModuleIO::write_Hexxs_csr(restart_HR_path, GlobalC::ucell, this->get_Hexxs()); + ModuleIO::write_Hexxs_csr(restart_HR_path, ucell, this->get_Hexxs()); if (GlobalV::MY_RANK == 0) { @@ -214,6 +224,7 @@ void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, const Uni } this->dm_last_step = dynamic_cast*>(&elec)->get_DM(); conv_esolver = this->exx_after_converge( + ucell, hamilt, *dynamic_cast*>(&elec)->get_DM(), kv, @@ -228,6 +239,7 @@ void Exx_LRI_Interface::exx_iter_finish(const K_Vectors& kv, const Uni template bool Exx_LRI_Interface::exx_after_converge( + const UnitCell& ucell, hamilt::Hamilt& hamilt, const elecstate::DensityMatrix& dm, const K_Vectors& kv, @@ -260,7 +272,7 @@ bool Exx_LRI_Interface::exx_after_converge( else { // update exx and redo scf - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); iter = 0; std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; this->two_level_step++; @@ -285,7 +297,7 @@ bool Exx_LRI_Interface::exx_after_converge( // update exx and redo scf if (this->two_level_step == 0) { - XC_Functional::set_xc_type(GlobalC::ucell.atoms[0].ncpp.xc_func); + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); } std::cout << " Updating EXX " << std::flush; @@ -305,11 +317,11 @@ bool Exx_LRI_Interface::exx_after_converge( : 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]); + // this->symrot_.test_HR_rotation(ucell.symm, ucell.atoms, 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())); + // this->symrot_.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, this->symrot_.get_Rs_from_adjacent_list(ucell, GlobalC::GridD, *lm.ParaV)); + // this->symrot_.test_HR_rotation(ucell.symm, ucell.atoms, ucell.st, 'H', *(dynamic_cast*>(&hamilt)->getHR())); // exit(0); if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) @@ -327,15 +339,15 @@ bool Exx_LRI_Interface::exx_after_converge( // ======================== 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())); + // this->symrot_.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, this->symrot_.get_Rs_from_adjacent_list(ucell, GlobalC::GridD, *lm.ParaV)); + // this->symrot_.test_HR_rotation(ucell.symm, ucell.atoms, 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))); + // symrot_.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, symrot_.get_Rs_from_adjacent_list(ucell, GlobalC::GridD, *this->DM->get_paraV_pointer())); + // symrot_.test_HR_rotation(ucell.symm, ucell.atoms, 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]); + // this->symrot_.test_HR_rotation(ucell.symm, ucell.atoms, ucell.st, 'H', this->exx_ptr->Hexxs[0]); // exit(0);// break after test // ======================== test ======================== iter = 0; From 39f62a371b2f3014d7e58289b30ab1c48597b4dc Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 16:49:43 +0800 Subject: [PATCH 02/24] change ucell in exx_lri/init --- source/module_esolver/esolver_ks_lcao.cpp | 5 +++-- source/module_lr/esolver_lrtd_lcao.cpp | 4 ++-- source/module_rdmft/rdmft.cpp | 16 ++++++++++++---- source/module_ri/Exx_LRI.h | 5 ++++- source/module_ri/Exx_LRI.hpp | 5 ++++- source/module_ri/RPA_LRI.h | 1 + source/module_ri/RPA_LRI.hpp | 3 ++- 7 files changed, 28 insertions(+), 11 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index e6a4ab6632..f4cea03d15 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -178,12 +178,12 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // 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, ucell,this->kv, orb_); this->exd->exx_before_all_runners(this->kv, ucell, this->pv); } else { - this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv, orb_); + this->exx_lri_complex->init(MPI_COMM_WORLD, ucell,this->kv, orb_); this->exc->exx_before_all_runners(this->kv, ucell, this->pv); } } @@ -1074,6 +1074,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) RPA_LRI rpa_lri_double(GlobalC::exx_info.info_ri); rpa_lri_double.cal_postSCF_exx(*dynamic_cast*>(this->pelec)->get_DM(), MPI_COMM_WORLD, + ucell, this->kv, orb_); rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 4b3b53a548..25a1fa7ac2 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -228,7 +228,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; } else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; } this->exx_lri = std::make_shared>(exx_info.info_ri); - this->exx_lri->init(MPI_COMM_WORLD, this->kv, ks_sol.orb_); + this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, ks_sol.orb_); this->exx_lri->cal_exx_ions(input.out_ri_cv); } } @@ -398,7 +398,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; } else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; } this->exx_lri = std::make_shared>(exx_info.info_ri); - this->exx_lri->init(MPI_COMM_WORLD, this->kv, orb); + this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, orb); this->exx_lri->cal_exx_ions(input.out_ri_cv); } // else diff --git a/source/module_rdmft/rdmft.cpp b/source/module_rdmft/rdmft.cpp index e50d63fd45..3a7f515c9b 100644 --- a/source/module_rdmft/rdmft.cpp +++ b/source/module_rdmft/rdmft.cpp @@ -55,8 +55,16 @@ RDMFT::~RDMFT() } template -void RDMFT::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& ParaV_in, UnitCell& ucell_in, - K_Vectors& kv_in, elecstate::ElecState& pelec_in, LCAO_Orbitals& orb_in, TwoCenterBundle& two_center_bundle_in, std::string XC_func_rdmft_in, double alpha_power_in) +void RDMFT::init(Gint_Gamma& GG_in, + Gint_k& GK_in, + Parallel_Orbitals& ParaV_in, + UnitCell& ucell_in, + K_Vectors& kv_in, + elecstate::ElecState& pelec_in, + LCAO_Orbitals& orb_in, + TwoCenterBundle& two_center_bundle_in, + std::string XC_func_rdmft_in, + double alpha_power_in) { GG = &GG_in; GK = &GK_in; @@ -155,12 +163,12 @@ void RDMFT::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& Pa if (GlobalC::exx_info.info_ri.real_number) { Vxc_fromRI_d = new Exx_LRI(GlobalC::exx_info.info_ri); - Vxc_fromRI_d->init(MPI_COMM_WORLD, *kv, *orb); + Vxc_fromRI_d->init(MPI_COMM_WORLD, ucell_in,*kv, *orb); } else { Vxc_fromRI_c = new Exx_LRI>(GlobalC::exx_info.info_ri); - Vxc_fromRI_c->init(MPI_COMM_WORLD, *kv, *orb); + Vxc_fromRI_c->init(MPI_COMM_WORLD, ucell_in,*kv, *orb); } } #endif diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index 6c57a1816c..6b990770c4 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -56,7 +56,10 @@ class Exx_LRI void reset_Cs(const std::map>>& Cs_in) { this->exx_lri.set_Cs(Cs_in, this->info.C_threshold); } void reset_Vs(const std::map>>& Vs_in) { this->exx_lri.set_Vs(Vs_in, this->info.V_threshold); } - void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb); + void init(const MPI_Comm &mpi_comm_in, + const UnitCell &ucell, + const K_Vectors &kv_in, + const LCAO_Orbitals& orb); void cal_exx_force(); void cal_exx_stress(); void cal_exx_ions(const bool write_cv = false); diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 58a900fe11..7d7d294965 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -26,7 +26,10 @@ #include template -void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb) +void Exx_LRI::init(const MPI_Comm &mpi_comm_in, + const UnitCell &ucell, + const K_Vectors &kv_in, + const LCAO_Orbitals& orb) { ModuleBase::TITLE("Exx_LRI","init"); ModuleBase::timer::tick("Exx_LRI", "init"); diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h index b169b0884b..554b2d1425 100644 --- a/source/module_ri/RPA_LRI.h +++ b/source/module_ri/RPA_LRI.h @@ -41,6 +41,7 @@ template class RPA_LRI void cal_rpa_cv(); void cal_postSCF_exx(const elecstate::DensityMatrix& dm, const MPI_Comm& mpi_comm_in, + const UnitCell& ucell, const K_Vectors& kv, const LCAO_Orbitals& orb); void out_for_RPA(const Parallel_Orbitals& parav, diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 3096b988f0..8db001f74a 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -73,6 +73,7 @@ void RPA_LRI::cal_rpa_cv() template void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix& dm, const MPI_Comm& mpi_comm_in, + const UnitCell& ucell, const K_Vectors& kv, const LCAO_Orbitals& orb) { @@ -105,7 +106,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix GlobalC::exx_info.info_global.hybrid_alpha = 1; GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times; - exx_lri_rpa.init(mpi_comm_in, kv, orb); + exx_lri_rpa.init(mpi_comm_in, ucell, kv, orb); exx_lri_rpa.cal_exx_ions(PARAM.inp.out_ri_cv); if (exx_spacegroup_symmetry && PARAM.inp.exx_symmetry_realspace) { From 4e7fe0792a191db4ed5f4729a926117d4cd98315 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:00:37 +0800 Subject: [PATCH 03/24] change ucell in exx_abfs-construct.cpp --- source/module_ri/Exx_LRI.hpp | 9 +++++---- source/module_ri/exx_abfs-construct_orbs.cpp | 3 ++- source/module_ri/exx_abfs-construct_orbs.h | 1 + 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 7d7d294965..b7b07b4575 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -53,6 +53,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, this->mpi_comm = mpi_comm_in; this->p_kv = &kv_in; this->orb_cutoff_ = orb.cutoffs(); + const double omega = ucell.omega; this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times ); @@ -66,10 +67,10 @@ 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::print_orbs_size(this->abfs, GlobalV::ofs_running); + } + Exx_Abfs::Construct_Orbs::print_orbs_size(ucell,this->abfs, GlobalV::ofs_running); - auto get_ccp_parameter = [this]() -> std::map + auto get_ccp_parameter = [this,&omega]() -> std::map { switch(this->info.ccp_type) { @@ -79,7 +80,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, { // 4/3 * pi * Rcut^3 = V_{supercell} = V_{unitcell} * Nk const int nspin0 = (PARAM.inp.nspin==2) ? 2 : 1; - const double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full()/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0); + const double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full()/nspin0 * omega / (ModuleBase::PI), 1.0/3.0); return {{"hf_Rcut", hf_Rcut}}; } case Conv_Coulomb_Pot_K::Ccp_Type::Hse: diff --git a/source/module_ri/exx_abfs-construct_orbs.cpp b/source/module_ri/exx_abfs-construct_orbs.cpp index 2452f49f77..dec5e12ef0 100644 --- a/source/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/module_ri/exx_abfs-construct_orbs.cpp @@ -479,6 +479,7 @@ inline const Numerical_Orbital_Lm &Exx_Abfs::Construct_Orbs::get_orbital( */ void Exx_Abfs::Construct_Orbs::print_orbs_size( + const UnitCell &ucell, const std::vector>> &orbs, std::ostream &os) { @@ -486,7 +487,7 @@ void Exx_Abfs::Construct_Orbs::print_orbs_size( const std::vector L_labels = {'s', 'p', 'd'}; for(std::size_t T=0; T>> &orbs, std::ostream &os); From 8192aae254d61d476c29498a9633dd3cb0bf25de Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:14:00 +0800 Subject: [PATCH 04/24] change ucell in module_ri/exx_abfs-jle.cpp --- source/module_ri/Exx_LRI_interface.hpp | 2 +- source/module_ri/exx_abfs-construct_orbs.cpp | 2 +- source/module_ri/exx_abfs-construct_orbs.h | 3 ++- source/module_ri/exx_abfs-jle.cpp | 8 +++++--- source/module_ri/exx_abfs-jle.h | 5 ++++- source/module_ri/exx_opt_orb.cpp | 4 ++-- source/module_ri/exx_opt_orb.h | 2 +- 7 files changed, 16 insertions(+), 10 deletions(-) diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 3800ca0419..72f456649f 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -90,7 +90,7 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, { //program should be stopped after this judgement Exx_Opt_Orb exx_opt_orb; - exx_opt_orb.generate_matrix(kv, orb); + exx_opt_orb.generate_matrix(kv, ucell,orb); ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf"); return; } diff --git a/source/module_ri/exx_abfs-construct_orbs.cpp b/source/module_ri/exx_abfs-construct_orbs.cpp index dec5e12ef0..ddf410a581 100644 --- a/source/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/module_ri/exx_abfs-construct_orbs.cpp @@ -479,7 +479,7 @@ inline const Numerical_Orbital_Lm &Exx_Abfs::Construct_Orbs::get_orbital( */ void Exx_Abfs::Construct_Orbs::print_orbs_size( - const UnitCell &ucell, + const UnitCell& ucell, 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 b0fbda9c34..1d805a362f 100644 --- a/source/module_ri/exx_abfs-construct_orbs.h +++ b/source/module_ri/exx_abfs-construct_orbs.h @@ -4,6 +4,7 @@ #include "exx_abfs.h" #include +#include "module_cell/unitcell.h" #include "../module_basis/module_ao/ORB_atomic_lm.h" class LCAO_Orbitals; @@ -25,7 +26,7 @@ class Exx_Abfs::Construct_Orbs const double times_threshold=0); static void print_orbs_size( - const UnitCell &ucell, + const UnitCell& ucell, const std::vector>> &orbs, std::ostream &os); diff --git a/source/module_ri/exx_abfs-jle.cpp b/source/module_ri/exx_abfs-jle.cpp index 9a43c6632d..f9086fe777 100644 --- a/source/module_ri/exx_abfs-jle.cpp +++ b/source/module_ri/exx_abfs-jle.cpp @@ -12,11 +12,13 @@ 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 ) +void Exx_Abfs::Jle::init_jle(const double kmesh_times, + const UnitCell& ucell, + const LCAO_Orbitals& orb) { - jle.resize( GlobalC::ucell.ntype ); + jle.resize( ucell.ntype ); - for (int T = 0; T < GlobalC::ucell.ntype ; T++) + for (int T = 0; T < ucell.ntype ; T++) { jle[T].resize( Lmax+1 ); for (int L=0; L <= Lmax ; ++L) diff --git a/source/module_ri/exx_abfs-jle.h b/source/module_ri/exx_abfs-jle.h index 967df31de6..d41078248e 100644 --- a/source/module_ri/exx_abfs-jle.h +++ b/source/module_ri/exx_abfs-jle.h @@ -4,6 +4,7 @@ #include "exx_abfs.h" #include +#include "module_cell/unitcell.h" #include "../module_basis/module_ao/ORB_read.h" class Exx_Abfs::Jle @@ -15,7 +16,9 @@ class Exx_Abfs::Jle std::vector< Numerical_Orbital_Lm>>> jle; - void init_jle( const double kmesh_times, const LCAO_Orbitals& orb ); + void init_jle(const double kmesh_times, + const UnitCell& ucell, + const LCAO_Orbitals& orb); static bool generate_matrix; static int Lmax; diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index d714de1a9d..ce5aecb175 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -15,7 +15,7 @@ #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 K_Vectors &kv, const UnitCell& ucell, const LCAO_Orbitals& orb) const { // std::ofstream ofs_mpi(GlobalC::exx_lcao.test_dir.process+"time_"+ModuleBase::GlobalFunc::TO_STRING(GlobalV::MY_RANK),std::ofstream::app); @@ -31,7 +31,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const LCAO_Orbitals& orb) // ofs_mpi<<"memory:\t"<kmesh_times, orb ); + jle.init_jle(this->kmesh_times, ucell , orb, ); // ofs_mpi<<"memory:\t"<>> cal_I( const std::map>>>> &ms, From 00eb36ab0f98892933562d61147e5b866d477ac3 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:14:34 +0800 Subject: [PATCH 05/24] change ucell in module_ri/exx_lip.hpp --- source/module_ri/exx_lip.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_ri/exx_lip.hpp b/source/module_ri/exx_lip.hpp index 94689f5167..1ce25150b1 100644 --- a/source/module_ri/exx_lip.hpp +++ b/source/module_ri/exx_lip.hpp @@ -540,7 +540,7 @@ void Exx_Lip::read_q_pack(const ModuleSymmetry::Symmetry& symm, this->q_pack->kv_ptr->ngk.data(), this->wfc_basis->npwk_max); // mohan update 2021-02-25 // this->q_pack->wf_ptr->init(this->q_pack->kv_ptr->get_nks(),this->q_pack->kv_ptr,this->ucell_ptr,old_pwptr,&ppcell,&GlobalC::ORB,&hm,&Pkpoints); - this->q_pack->wf_ptr->table_local.create(GlobalC::ucell.ntype, GlobalC::ucell.nmax_total, PARAM.globalv.nqx); + this->q_pack->wf_ptr->table_local.create(ucell.ntype, ucell.nmax_total, PARAM.globalv.nqx); // this->q_pack->wf_ptr->table_local.create(this->q_pack->wf_ptr->this->ucell_ptr->ntype, this->q_pack->wf_ptr->this->ucell_ptr->nmax_total, PARAM.globalv.nqx); #ifdef __LCAO Wavefunc_in_pw::make_table_q(GlobalC::ORB.orbital_file, this->q_pack->wf_ptr->table_local); From 0d95907f633657b97442a8479dc078e5bd0133f3 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:23:27 +0800 Subject: [PATCH 06/24] change ucell in module_ri/exx_lri.hpp --- source/module_lr/esolver_lrtd_lcao.cpp | 4 ++-- source/module_rdmft/update_state_rdmft.cpp | 10 ++++++---- source/module_ri/Exx_LRI.h | 2 +- source/module_ri/Exx_LRI.hpp | 3 ++- source/module_ri/Exx_LRI_interface.hpp | 2 +- source/module_ri/RPA_LRI.hpp | 8 ++++---- source/module_ri/exx_opt_orb.cpp | 2 +- 7 files changed, 17 insertions(+), 14 deletions(-) diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 25a1fa7ac2..aa66831cfd 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -229,7 +229,7 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; } this->exx_lri = std::make_shared>(exx_info.info_ri); this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, ks_sol.orb_); - this->exx_lri->cal_exx_ions(input.out_ri_cv); + this->exx_lri->cal_exx_ions(ucell,input.out_ri_cv); } } #endif @@ -399,7 +399,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; } this->exx_lri = std::make_shared>(exx_info.info_ri); this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, orb); - this->exx_lri->cal_exx_ions(input.out_ri_cv); + this->exx_lri->cal_exx_ions(ucell,input.out_ri_cv); } // else #endif diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index 3e6bbad183..358d8f8b05 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -16,8 +16,10 @@ namespace rdmft template -void RDMFT::update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis_in, - ModuleBase::matrix& vloc_in, ModuleBase::ComplexMatrix& sf_in) +void RDMFT::update_ion(UnitCell& ucell_in, + ModulePW::PW_Basis& rho_basis_in, + ModuleBase::matrix& vloc_in, + ModuleBase::ComplexMatrix& sf_in) { ucell = &ucell_in; rho_basis = &rho_basis_in; @@ -31,11 +33,11 @@ void RDMFT::update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis { if (GlobalC::exx_info.info_ri.real_number) { - Vxc_fromRI_d->cal_exx_ions(); + Vxc_fromRI_d->cal_exx_ions(ucell_in); } else { - Vxc_fromRI_c->cal_exx_ions(); + Vxc_fromRI_c->cal_exx_ions(ucell_in); } } #endif diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index 6b990770c4..e23197f1f7 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -62,7 +62,7 @@ class Exx_LRI const LCAO_Orbitals& orb); void cal_exx_force(); void cal_exx_stress(); - void cal_exx_ions(const bool write_cv = false); + void cal_exx_ions(const UnitCell& ucell, const bool write_cv = false); void cal_exx_elec(const std::vector>>>& Ds, const Parallel_Orbitals& pv, const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr); diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index b7b07b4575..7b6058de9d 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -105,7 +105,8 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, } template -void Exx_LRI::cal_exx_ions(const bool write_cv) +void Exx_LRI::cal_exx_ions(const UnitCell& ucell, + const bool write_cv) { ModuleBase::TITLE("Exx_LRI","cal_exx_ions"); ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions"); diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 72f456649f..fcbcadcc3f 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -83,7 +83,7 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, XC_Functional::set_xc_type("pbe"); } } - this->exx_ptr->cal_exx_ions(PARAM.inp.out_ri_cv); + this->exx_ptr->cal_exx_ions(ucell,PARAM.inp.out_ri_cv); } if (Exx_Abfs::Jle::generate_matrix) diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 8db001f74a..84350e60d6 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -89,9 +89,9 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix if (exx_spacegroup_symmetry) { const std::array period = RI_Util::get_Born_vonKarmen_period(kv); - symrot.find_irreducible_sector(GlobalC::ucell.symm, GlobalC::ucell.atoms, GlobalC::ucell.st, - RI_Util::get_Born_von_Karmen_cells(period), period, GlobalC::ucell.lat); - symrot.cal_Ms(kv, GlobalC::ucell, *dm.get_paraV_pointer()); + symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st, + RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat); + symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer()); mix_DMk_2D.mix(symrot.restore_dm(kv, dm.get_DMK_vector(), *dm.get_paraV_pointer()), true); } else { mix_DMk_2D.mix(dm.get_DMK_vector(), true); } @@ -107,7 +107,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix GlobalC::exx_info.info_ri.ccp_rmesh_times = PARAM.inp.rpa_ccp_rmesh_times; exx_lri_rpa.init(mpi_comm_in, ucell, kv, orb); - exx_lri_rpa.cal_exx_ions(PARAM.inp.out_ri_cv); + exx_lri_rpa.cal_exx_ions(ucell,PARAM.inp.out_ri_cv); if (exx_spacegroup_symmetry && PARAM.inp.exx_symmetry_realspace) { exx_lri_rpa.cal_exx_elec(Ds, *dm.get_paraV_pointer(), &symrot); diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index ce5aecb175..28cb5b112e 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -31,7 +31,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co // ofs_mpi<<"memory:\t"<kmesh_times, ucell , orb, ); + jle.init_jle(this->kmesh_times, ucell , orb); // ofs_mpi<<"memory:\t"< Date: Sun, 8 Dec 2024 17:37:01 +0800 Subject: [PATCH 07/24] change ucell in module_ri/exx_lir.hpp/cal_exx_elec --- source/module_esolver/esolver_ks_lcao.cpp | 4 +++- source/module_rdmft/rdmft.h | 4 ++-- source/module_rdmft/rdmft_pot.cpp | 10 +++++----- source/module_rdmft/update_state_rdmft.cpp | 8 +++++--- source/module_ri/Exx_LRI.h | 1 + source/module_ri/Exx_LRI.hpp | 15 ++++++++------- source/module_ri/Exx_LRI_interface.h | 7 +++++-- source/module_ri/Exx_LRI_interface.hpp | 22 ++++++++++++++++------ source/module_ri/RPA_LRI.hpp | 4 ++-- 9 files changed, 47 insertions(+), 28 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index f4cea03d15..dce1ed85f6 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -623,6 +623,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const if (GlobalC::exx_info.info_ri.real_number) { this->exd->exx_eachiterinit(istep, + ucell, *dynamic_cast*>(this->pelec)->get_DM(), this->kv, iter); @@ -630,6 +631,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const else { this->exc->exx_eachiterinit(istep, + ucell, *dynamic_cast*>(this->pelec)->get_DM(), this->kv, iter); @@ -1052,7 +1054,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) occ_number_ks(ik, inb) /= this->kv.wk[ik]; } } - this->rdmft_solver.update_elec(occ_number_ks, *(this->psi)); + this->rdmft_solver.update_elec(ucell,occ_number_ks, *(this->psi)); //! initialize the gradients of Etotal with respect to occupation numbers and wfc, //! and set all elements to 0. diff --git a/source/module_rdmft/rdmft.h b/source/module_rdmft/rdmft.h index 85e07c0613..5f4c5848ef 100644 --- a/source/module_rdmft/rdmft.h +++ b/source/module_rdmft/rdmft.h @@ -90,7 +90,7 @@ class RDMFT //! update in elec-step // Or we can use rdmft_solver.wfc/occ_number directly when optimizing, so that the update_elec() function does not require parameters. - void update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in = nullptr); + void update_elec(const UnitCell& ucell, const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in = nullptr); //! obtain the gradient of total energy with respect to occupation number and wfc double cal_E_grad_wfc_occ_num(); @@ -117,7 +117,7 @@ class RDMFT void cal_V_hartree(); //! construct V_XC based on different XC_functional( i.e. RDMFT class member XC_func_rdmft) - void cal_V_XC(); + void cal_V_XC(const UnitCell& ucell); //! get the total Hamilton in k-space void cal_Hk_Hpsi(); diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index a581910d6d..3f7f20868f 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -175,7 +175,7 @@ void RDMFT::cal_V_hartree() template -void RDMFT::cal_V_XC() +void RDMFT::cal_V_XC(const UnitCell& ucell) { // // //test // DM_XC_pass = DM_XC; @@ -280,11 +280,11 @@ void RDMFT::cal_V_XC() // provide the Ds_XC to Vxc_fromRI(V_exx_XC) if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { - Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV, &this->symrot_exx); + Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, ucell,*ParaV, &this->symrot_exx); } else { - Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV); + Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, ucell,*ParaV); } // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC @@ -308,11 +308,11 @@ void RDMFT::cal_V_XC() // // provide the Ds_XC to Vxc_fromRI(V_exx_XC) if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { - Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV, &this->symrot_exx); + Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, ucell,*ParaV, &this->symrot_exx); } else { - Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV); + Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, ucell,*ParaV); } // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index 358d8f8b05..4f2d329377 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -47,7 +47,9 @@ void RDMFT::update_ion(UnitCell& ucell_in, template -void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in) +void RDMFT::update_elec(const UnitCell& ucell, + const ModuleBase::matrix& occ_number_in, + const psi::Psi& wfc_in, const Charge* charge_in) { // update occ_number, wg, wk_fun_occNum occ_number = (occ_number_in); @@ -76,11 +78,11 @@ void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const p if( this->cal_E_type != 1 ) { // the second cal_E_type need the complete pot to get effctive_V to calEband and so on. - this->pelec->pot->update_from_charge(charge, ucell); + this->pelec->pot->update_from_charge(charge, &ucell); } this->cal_V_hartree(); - this->cal_V_XC(); + this->cal_V_XC(ucell); // this->cal_Hk_Hpsi(); std::cout << "\n******\n" << "update elec in rdmft successfully" << "\n******\n" << std::endl; diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index e23197f1f7..ce40c2de46 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -64,6 +64,7 @@ class Exx_LRI void cal_exx_stress(); void cal_exx_ions(const UnitCell& ucell, const bool write_cv = false); void cal_exx_elec(const std::vector>>>& Ds, + const UnitCell& ucell, const Parallel_Orbitals& pv, const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr); std::vector> get_abfs_nchis() const; diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 7b6058de9d..b0993e9974 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -116,18 +116,18 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, // 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(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), - RI_Util::Vector3_to_array3(GlobalC::ucell.a3)}; + = {RI_Util::Vector3_to_array3(ucell.a1), + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3)}; const std::array period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]}; this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period); @@ -190,6 +190,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, template void Exx_LRI::cal_exx_elec(const std::vector>>>& Ds, + const UnitCell& ucell, const Parallel_Orbitals& pv, const ModuleSymmetry::Symmetry_rotation* p_symrot) { diff --git a/source/module_ri/Exx_LRI_interface.h b/source/module_ri/Exx_LRI_interface.h index 42172aada1..9c18ebad84 100644 --- a/source/module_ri/Exx_LRI_interface.h +++ b/source/module_ri/Exx_LRI_interface.h @@ -57,8 +57,11 @@ class Exx_LRI_Interface void exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb); /// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF - void exx_eachiterinit(const int istep, const elecstate::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, - const K_Vectors& kv, const int& iter); + void exx_eachiterinit(const int istep, + const UnitCell& ucell, + const elecstate::DensityMatrix& dm/**< double should be Tdata if complex-PBE-DM is supported*/, + const K_Vectors& kv, + const int& iter); /// @brief in hamilt2density: calculate Hexx and Eexx void exx_hamilt2density(elecstate::ElecState& elec, const Parallel_Orbitals& pv, const int iter); diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index fcbcadcc3f..8f159be632 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -114,14 +114,18 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, } template -void Exx_LRI_Interface::exx_eachiterinit(const int istep, const elecstate::DensityMatrix& dm, const K_Vectors& kv, const int& iter) +void Exx_LRI_Interface::exx_eachiterinit(const int istep, + const UnitCell& ucell, + const elecstate::DensityMatrix& dm, + const K_Vectors& kv, + const int& iter) { if (GlobalC::exx_info.info_global.cal_exx) { if (!GlobalC::exx_info.info_global.separate_loop && (this->two_level_step || istep > 0)) { const bool flag_restart = (iter == 1) ? true : false; - auto cal = [this, &kv, &flag_restart](const elecstate::DensityMatrix& dm_in) + auto cal = [this, &ucell,&kv, &flag_restart](const elecstate::DensityMatrix& dm_in) { if (this->exx_spacegroup_symmetry) { this->mix_DMk_2D.mix(symrot_.restore_dm(kv,dm_in.get_DMK_vector(), *dm_in.get_paraV_pointer()), flag_restart); } else { this->mix_DMk_2D.mix(dm_in.get_DMK_vector(), flag_restart); } @@ -129,8 +133,14 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, const elecst 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_in.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_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry); - if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer(), &this->symrot_); } - else { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer()); } + if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + { + this->exx_ptr->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); + } + else + { + this->exx_ptr->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer()); + } }; if(istep > 0 && flag_restart) { cal(*dm_last_step); @@ -326,13 +336,13 @@ bool Exx_LRI_Interface::exx_after_converge( 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_); + this->exx_ptr->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer(), &this->symrot_); // this->symrot_.print_HR(this->exx_ptr->Hexxs[0], "Hexxs_irreducible"); // test // this->symrot_.print_HR(this->exx_ptr->Hexxs[0], "Hexxs_restored", 1e-10); // test } else { - this->exx_ptr->cal_exx_elec(Ds, *dm.get_paraV_pointer()); // restore DM but not Hexx + this->exx_ptr->cal_exx_elec(Ds, ucell, *dm.get_paraV_pointer()); // restore DM but not Hexx // 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 } diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 84350e60d6..122879640b 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -110,9 +110,9 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix exx_lri_rpa.cal_exx_ions(ucell,PARAM.inp.out_ri_cv); if (exx_spacegroup_symmetry && PARAM.inp.exx_symmetry_realspace) { - exx_lri_rpa.cal_exx_elec(Ds, *dm.get_paraV_pointer(), &symrot); + exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer(), &symrot); } else { - exx_lri_rpa.cal_exx_elec(Ds, *dm.get_paraV_pointer()); + exx_lri_rpa.cal_exx_elec(Ds, ucell,*dm.get_paraV_pointer()); } // cout<<"postSCF_Eexx: "< Date: Sun, 8 Dec 2024 17:42:11 +0800 Subject: [PATCH 08/24] change ucell in module_ri/exx_lri.hpp --- .../hamilt_lcaodft/FORCE_STRESS.cpp | 8 ++++---- source/module_ri/Exx_LRI.h | 4 ++-- source/module_ri/Exx_LRI.hpp | 18 +++++++++--------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 389119a17a..43a633daf2 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -375,12 +375,12 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { if (GlobalC::exx_info.info_ri.real_number) { - exx_lri_double.cal_exx_force(); + exx_lri_double.cal_exx_force(ucell.nat); force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_double.force_exx; } else { - exx_lri_complex.cal_exx_force(); + exx_lri_complex.cal_exx_force(ucell.nat); force_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_complex.force_exx; } } @@ -388,12 +388,12 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { if (GlobalC::exx_info.info_ri.real_number) { - exx_lri_double.cal_exx_stress(); + exx_lri_double.cal_exx_stress(ucell.omega,ucell.lat0); stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_double.stress_exx; } else { - exx_lri_complex.cal_exx_stress(); + exx_lri_complex.cal_exx_stress(ucell.omega,ucell.lat0); stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_complex.stress_exx; } } diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index ce40c2de46..cbcef837a9 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -60,8 +60,8 @@ class Exx_LRI const UnitCell &ucell, const K_Vectors &kv_in, const LCAO_Orbitals& orb); - void cal_exx_force(); - void cal_exx_stress(); + void cal_exx_force(const int& nat); + void cal_exx_stress(const double& omega, const double& lat0); void cal_exx_ions(const UnitCell& ucell, const bool write_cv = false); void cal_exx_elec(const std::vector>>>& Ds, const UnitCell& ucell, diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index b0993e9974..64da2bcc44 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -219,7 +219,7 @@ void Exx_LRI::cal_exx_elec(const std::vectorexx_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); + Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, 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], @@ -274,20 +274,20 @@ post_process_old */ template -void Exx_LRI::cal_exx_force() +void Exx_LRI::cal_exx_force(const int& nat) { ModuleBase::TITLE("Exx_LRI","cal_exx_force"); ModuleBase::timer::tick("Exx_LRI", "cal_exx_force"); - this->force_exx.create(GlobalC::ucell.nat, Ndim); + this->force_exx.create(nat, Ndim); for(int is=0; isexx_lri.cal_force({"","",std::to_string(is),"",""}); 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? @@ -298,7 +298,7 @@ void Exx_LRI::cal_exx_force() template -void Exx_LRI::cal_exx_stress() +void Exx_LRI::cal_exx_stress(const double& omega, const double& lat0) { ModuleBase::TITLE("Exx_LRI","cal_exx_stress"); ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress"); @@ -310,12 +310,12 @@ 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? - const double frac = 2 * SPIN_multiple / GlobalC::ucell.omega * GlobalC::ucell.lat0; // why? + const double frac = 2 * SPIN_multiple / omega * lat0; // why? this->stress_exx *= frac; ModuleBase::timer::tick("Exx_LRI", "cal_exx_stress"); From b158b78b1889b183f754071ee2a3cd4a53b1a1ee Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:45:44 +0800 Subject: [PATCH 09/24] change ucell in module_ri/exx_opt-print.hpp --- source/module_ri/exx_opt_orb-print.cpp | 45 +++++++++++++------------- source/module_ri/exx_opt_orb.cpp | 12 ++++--- source/module_ri/exx_opt_orb.h | 1 + 3 files changed, 32 insertions(+), 26 deletions(-) diff --git a/source/module_ri/exx_opt_orb-print.cpp b/source/module_ri/exx_opt_orb-print.cpp index 52505505d9..a768d8399b 100644 --- a/source/module_ri/exx_opt_orb-print.cpp +++ b/source/module_ri/exx_opt_orb-print.cpp @@ -3,6 +3,7 @@ #include "exx_abfs-jle.h" void Exx_Opt_Orb::print_matrix( + const UnitCell& ucell, const K_Vectors &kv, const std::string& file_name, const std::vector> &matrix_Q, @@ -15,47 +16,47 @@ void Exx_Opt_Orb::print_matrix( { auto print_header = [&]( std::ofstream &ofs ) { - ofs << GlobalC::ucell.lat0 << std::endl; + ofs << ucell.lat0 << std::endl; - 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; + ofs << ucell.latvec.e11 << " " << ucell.latvec.e12 << " " << ucell.latvec.e13 << std::endl; + ofs << ucell.latvec.e21 << " " << ucell.latvec.e22 << " " << ucell.latvec.e23 << std::endl; + ofs << ucell.latvec.e31 << " " << ucell.latvec.e32 << " " << ucell.latvec.e33 << std::endl; if( TA==TB ) { ofs << 1 << " ntype" << std::endl; - ofs << GlobalC::ucell.atoms[TA].label << " label" << std::endl; + ofs << ucell.atoms[TA].label << " label" << std::endl; if( IA==IB ) { ofs << 1 << " na" << std::endl; - 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 << ucell.atoms[TA].tau[IA].x << " " + << ucell.atoms[TA].tau[IA].y << " " + << ucell.atoms[TA].tau[IA].z << std::endl; } else { ofs << 2 << " na" << std::endl; - 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 << " " - << GlobalC::ucell.atoms[TB].tau[IB].z << std::endl; + ofs << ucell.atoms[TA].tau[IA].x << " " + << ucell.atoms[TA].tau[IA].y << " " + << ucell.atoms[TA].tau[IA].z << std::endl; + ofs << ucell.atoms[TB].tau[IB].x << " " + << ucell.atoms[TB].tau[IB].y << " " + << ucell.atoms[TB].tau[IB].z << std::endl; } } else { ofs << 2 << " ntype" << std::endl; - ofs << GlobalC::ucell.atoms[TA].label << " label" << std::endl; + ofs << 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 << " " - << GlobalC::ucell.atoms[TA].tau[IA].z << std::endl; - ofs << GlobalC::ucell.atoms[TB].label << " label" << std::endl; + ofs << ucell.atoms[TA].tau[IA].x << " " + << ucell.atoms[TA].tau[IA].y << " " + << ucell.atoms[TA].tau[IA].z << std::endl; + ofs << 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 << " " - << GlobalC::ucell.atoms[TB].tau[IB].z << std::endl; + ofs << ucell.atoms[TB].tau[IB].x << " " + << ucell.atoms[TB].tau[IB].y << " " + << ucell.atoms[TB].tau[IB].z << std::endl; } // ecutwfc_jlq determine the jlq corresponding to plane wave calculation. diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 28cb5b112e..f2a96b5a47 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -208,7 +208,8 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co {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, + print_matrix(ucell, + kv, "matrix", m_lcaoslcaos_jys_proj, m_jys_jys_proj, @@ -219,7 +220,8 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co } else { - print_matrix(kv, + print_matrix(ucell, + kv, "matrix", ms_lcaoslcaos_jys.at(T).at(I).at(T).at(I), {{ms_jys_jys.at(T).at(I).at(T).at(I)}}, @@ -276,7 +278,8 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co { 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, + print_matrix(ucell, + kv, "matrix", m_lcaoslcaos_jys_proj, m_jys_jys_proj, @@ -287,7 +290,8 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co } else { - print_matrix(kv, + print_matrix(ucell, + kv, "matrix", 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)}, diff --git a/source/module_ri/exx_opt_orb.h b/source/module_ri/exx_opt_orb.h index fee6d70bb7..1da8e63d1e 100644 --- a/source/module_ri/exx_opt_orb.h +++ b/source/module_ri/exx_opt_orb.h @@ -24,6 +24,7 @@ class Exx_Opt_Orb const std::vector>> & m_middle, const std::vector> & m_right ) const; void print_matrix( + const UnitCell& ucell, const K_Vectors &kv, const std::string& file_name, const std::vector> &matrix_Q, From 9bb0813a6b62aae0b51270dd2aeaaddaa1612420 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 17:48:44 +0800 Subject: [PATCH 10/24] change ucell in module_ri/exx_opt.cpp --- source/module_ri/exx_opt_orb.cpp | 30 +++++++++++++++--------------- source/module_ri/exx_opt_orb.h | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index f2a96b5a47..0bc9dfd231 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -53,7 +53,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co // ofs_mpi<>> radial_R = get_radial_R(); + std::map>> radial_R = get_radial_R(ucell); #if TEST_EXX_RADIAL==2 { for(const auto & rA : radial_R) @@ -172,13 +172,13 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co // ofs_matrixes(GlobalC::exx_lcao.test_dir.matrix+"ms_lcaoslcaos_abfs",ms_lcaoslcaos_abfs); // ofs_matrixes(GlobalC::exx_lcao.test_dir.matrix+"ms_jys_abfs",ms_jys_abfs); - for( size_t TA=0; TA!=GlobalC::ucell.ntype; ++TA ) + for( size_t TA=0; TA!=ucell.ntype; ++TA ) { - for( size_t IA=0; IA!=GlobalC::ucell.atoms[TA].na; ++IA ) + for( size_t IA=0; IA!=ucell.atoms[TA].na; ++IA ) { - for( size_t TB=0; TB!=GlobalC::ucell.ntype; ++TB ) + for( size_t TB=0; TB!=ucell.ntype; ++TB ) { - for( size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) + for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB ) { if( TA==TB && IA==IB ) { @@ -358,23 +358,23 @@ std::vector>> Exx_Opt_Orb::cal_I( } } -std::map>> Exx_Opt_Orb::get_radial_R() const +std::map>> Exx_Opt_Orb::get_radial_R(const UnitCell& ucell) const { 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 IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) + for( size_t TA=0; TA!=ucell.ntype; ++TA ) { + for( size_t IA=0; IA!=ucell.atoms[TA].na; ++IA ) { + for( size_t TB=0; TB!=ucell.ntype; ++TB ) { + for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB ) { - const ModuleBase::Vector3 &tauA = GlobalC::ucell.atoms[TA].tau[IA]; - const ModuleBase::Vector3 &tauB = GlobalC::ucell.atoms[TB].tau[IB]; + const ModuleBase::Vector3 &tauA = ucell.atoms[TA].tau[IA]; + const ModuleBase::Vector3 &tauB = ucell.atoms[TB].tau[IB]; const double delta_R = (-tauA+tauB).norm(); 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 1da8e63d1e..2b68dfa40a 100644 --- a/source/module_ri/exx_opt_orb.h +++ b/source/module_ri/exx_opt_orb.h @@ -34,7 +34,7 @@ class Exx_Opt_Orb 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; + std::map>> get_radial_R(const UnitCell& ucell) const; int kmesh_times = 4; }; From 579e58a9b5fc1afc7bb5403838d8f294ffd342f4 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 19:57:04 +0800 Subject: [PATCH 11/24] change ucell in module_ri/rpa_lri.cpp --- source/module_esolver/esolver_ks_lcao.cpp | 2 +- source/module_ri/Exx_LRI.hpp | 5 +++-- source/module_ri/LRI_CV.h | 4 ++++ source/module_ri/LRI_CV.hpp | 26 +++++++++++++---------- source/module_ri/RPA_LRI.h | 5 +++-- source/module_ri/RPA_LRI.hpp | 17 +++++++++------ 6 files changed, 36 insertions(+), 23 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index dce1ed85f6..9df3070827 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -1080,7 +1080,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) this->kv, orb_); rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); - rpa_lri_double.out_for_RPA(this->pv, *(this->psi), this->pelec); + rpa_lri_double.out_for_RPA(ucell,this->pv, *(this->psi), this->pelec); } #endif diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 64da2bcc44..45d8f17acb 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -138,7 +138,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false); std::map>> - Vs = this->cv.cal_Vs( + Vs = this->cv.cal_Vs(ucell, list_As_Vs.first, list_As_Vs.second[0], {{"writable_Vws",true}}); this->cv.Vws = LRI_CV_Tools::get_CVws(Vs); @@ -148,7 +148,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, if(PARAM.inp.cal_force || PARAM.inp.cal_stress) { std::array>>,3> - dVs = this->cv.cal_dVs( + dVs = this->cv.cal_dVs(ucell, list_As_Vs.first, list_As_Vs.second[0], {{"writable_dVws",true}}); this->cv.dVws = LRI_CV_Tools::get_dCVws(dVs); @@ -166,6 +166,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, std::pair>>, std::array>>,3>> Cs_dCs = this->cv.cal_Cs_dCs( + ucell, list_As_Cs.first, list_As_Cs.second[0], {{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress}, {"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}}); diff --git a/source/module_ri/LRI_CV.h b/source/module_ri/LRI_CV.h index 252259c119..17df63ff0e 100644 --- a/source/module_ri/LRI_CV.h +++ b/source/module_ri/LRI_CV.h @@ -42,17 +42,20 @@ class LRI_CV const double &ccp_rmesh_times_in); inline std::map>> cal_Vs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags); // "writable_Vws" inline std::array>>,3> cal_dVs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags); // "writable_dVws" std::pair>>, std::array>>,3>> cal_Cs_dCs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags); // "cal_dC", "writable_Cws", "writable_dCws", "writable_Vws", "writable_dVws" @@ -91,6 +94,7 @@ class LRI_CV template std::map> cal_datas( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags, diff --git a/source/module_ri/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp index 8640088130..a6c764122a 100644 --- a/source/module_ri/LRI_CV.hpp +++ b/source/module_ri/LRI_CV.hpp @@ -76,6 +76,7 @@ void LRI_CV::set_orbitals( template template auto LRI_CV::cal_datas( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags, @@ -96,17 +97,17 @@ auto LRI_CV::cal_datas( const TA iat0 = list_A0[i0]; const TA iat1 = list_A1[i1].first; const TC &cell1 = list_A1[i1].second; - const int it0 = GlobalC::ucell.iat2it[iat0]; - const int ia0 = GlobalC::ucell.iat2ia[iat0]; - const int it1 = GlobalC::ucell.iat2it[iat1]; - const int ia1 = GlobalC::ucell.iat2ia[iat1]; - const ModuleBase::Vector3 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; - const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; + const int it0 = ucell.iat2it[iat0]; + const int ia0 = ucell.iat2ia[iat0]; + const int it1 = ucell.iat2it[iat1]; + const int ia1 = ucell.iat2ia[iat1]; + const ModuleBase::Vector3 tau0 = ucell.atoms[it0].tau[ia0]; + const ModuleBase::Vector3 tau1 = ucell.atoms[it1].tau[ia1]; const double Rcut = std::min( orb_cutoff_[it0] * rmesh_times + orb_cutoff_[it1], orb_cutoff_[it1] * rmesh_times + orb_cutoff_[it0]); - const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); - if( R_delta.norm()*GlobalC::ucell.lat0 < Rcut ) + const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec); + if( R_delta.norm()*ucell.lat0 < Rcut ) { const Tresult Data = func_DPcal_data(it0, it1, R_delta, flags); // if(Data.norm(std::numeric_limits::max()) > threshold) @@ -124,6 +125,7 @@ auto LRI_CV::cal_datas( template auto LRI_CV::cal_Vs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags) // + "writable_Vws" @@ -134,11 +136,12 @@ auto LRI_CV::cal_Vs( func_DPcal_V = std::bind( &LRI_CV::DPcal_V, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4); - return this->cal_datas(list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_V); + return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_V); } template auto LRI_CV::cal_dVs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags) // + "writable_dVws" @@ -150,11 +153,12 @@ auto LRI_CV::cal_dVs( &LRI_CV::DPcal_dV, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4); return LRI_CV_Tools::change_order( - this->cal_datas(list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV)); + this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV)); } template auto LRI_CV::cal_Cs_dCs( + const UnitCell &ucell, const std::vector &list_A0, const std::vector &list_A1, const std::map &flags) // "cal_dC" + "writable_Cws", "writable_dCws", "writable_Vws", "writable_dVws" @@ -166,7 +170,7 @@ auto LRI_CV::cal_Cs_dCs( &LRI_CV::DPcal_C_dC, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4); std::map, std::array,3>>>> - Cs_dCs_tmp = this->cal_datas(list_A0, list_A1, flags, std::min(1.0,this->ccp_rmesh_times), func_DPcal_C_dC); + Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, std::min(1.0,this->ccp_rmesh_times), func_DPcal_C_dC); std::map>> Cs; std::array>>,3> dCs; diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h index 554b2d1425..8b27c76bfe 100644 --- a/source/module_ri/RPA_LRI.h +++ b/source/module_ri/RPA_LRI.h @@ -38,13 +38,14 @@ template class RPA_LRI } ~RPA_LRI(){}; void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const std::vector& orb_cutoff); - void cal_rpa_cv(); + void cal_rpa_cv(const UnitCell &ucell); void cal_postSCF_exx(const elecstate::DensityMatrix& dm, const MPI_Comm& mpi_comm_in, const UnitCell& ucell, const K_Vectors& kv, const LCAO_Orbitals& orb); - void out_for_RPA(const Parallel_Orbitals& parav, + void out_for_RPA(const UnitCell& ucell, + const Parallel_Orbitals& parav, const psi::Psi& psi, const elecstate::ElecState* pelec); void out_eigen_vector(const Parallel_Orbitals& parav, const psi::Psi& psi); diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 122879640b..35dcb2184f 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -31,10 +31,10 @@ void RPA_LRI::init(const MPI_Comm& mpi_comm_in, const K_Vectors& kv_in } template -void RPA_LRI::cal_rpa_cv() +void RPA_LRI::cal_rpa_cv(const UnitCell& ucell) { - std::vector atoms(GlobalC::ucell.nat); - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) + std::vector atoms(ucell.nat); + for (int iat = 0; iat < ucell.nat; ++iat) { atoms[iat] = iat; } @@ -44,7 +44,8 @@ void RPA_LRI::cal_rpa_cv() const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms(this->mpi_comm, atoms, period_Vs, 2, false); - std::map>> Vs = exx_lri_rpa.cv.cal_Vs(list_As_Vs.first, + std::map>> Vs = exx_lri_rpa.cv.cal_Vs(ucell, + list_As_Vs.first, list_As_Vs.second[0], { {"writable_Vws", true} @@ -57,7 +58,8 @@ void RPA_LRI::cal_rpa_cv() std::pair>>, std::array>>, 3>> - Cs_dCs = exx_lri_rpa.cv.cal_Cs_dCs(list_As_Cs.first, + Cs_dCs = exx_lri_rpa.cv.cal_Cs_dCs(ucell, + list_As_Cs.first, list_As_Cs.second[0], { {"cal_dC", false}, @@ -118,7 +120,8 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix } template -void RPA_LRI::out_for_RPA(const Parallel_Orbitals& parav, +void RPA_LRI::out_for_RPA(const UnitCell& ucell, + const Parallel_Orbitals& parav, const psi::Psi& psi, const elecstate::ElecState* pelec) { @@ -127,7 +130,7 @@ void RPA_LRI::out_for_RPA(const Parallel_Orbitals& parav, this->out_eigen_vector(parav, psi); this->out_struc(); - this->cal_rpa_cv(); + this->cal_rpa_cv(ucell); std::cout << "rpa_pca_threshold: " << this->info.pca_threshold << std::endl; std::cout << "rpa_ccp_rmesh_times: " << this->info.ccp_rmesh_times << std::endl; std::cout << "rpa_lcao_exx(Ha): " << std::fixed << std::setprecision(15) << exx_lri_rpa.Eexx / 2.0 << std::endl; From 14822f837db9824afd4c86ee662f69c64f57922f Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 20:04:25 +0800 Subject: [PATCH 12/24] change ucell in module_ri/rpa_lri_tool.cpp --- source/module_ri/Exx_LRI.hpp | 16 ++++++------ source/module_ri/LRI_CV_Tools.h | 7 +++++- source/module_ri/LRI_CV_Tools.hpp | 41 +++++++++++++++++-------------- source/module_ri/RPA_LRI.hpp | 4 +-- 4 files changed, 39 insertions(+), 29 deletions(-) diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 45d8f17acb..18cff2b2ee 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -133,7 +133,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, 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, ucell, orb_cutoff_); const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false); @@ -141,7 +141,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, Vs = this->cv.cal_Vs(ucell, list_As_Vs.first, list_As_Vs.second[0], {{"writable_Vws",true}}); - this->cv.Vws = LRI_CV_Tools::get_CVws(Vs); + this->cv.Vws = LRI_CV_Tools::get_CVws(ucell,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); @@ -151,16 +151,16 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, dVs = this->cv.cal_dVs(ucell, list_As_Vs.first, list_As_Vs.second[0], {{"writable_dVws",true}}); - this->cv.dVws = LRI_CV_Tools::get_dCVws(dVs); + this->cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs); this->exx_lri.set_dVs(std::move(dVs), this->info.V_grad_threshold); if(PARAM.inp.cal_stress) { - std::array>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(dVs); + std::array>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs); this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold); } } - 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, ucell,orb_cutoff_); const std::pair, std::vector>>>> list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false); @@ -171,18 +171,18 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, {{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress}, {"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); + this->cv.Cws = LRI_CV_Tools::get_CVws(ucell,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) { std::array>>,3> &dCs = std::get<1>(Cs_dCs); - this->cv.dCws = LRI_CV_Tools::get_dCVws(dCs); + this->cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs); this->exx_lri.set_dCs(std::move(dCs), this->info.C_grad_threshold); if(PARAM.inp.cal_stress) { - std::array>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(dCs); + std::array>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs); this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold); } } diff --git a/source/module_ri/LRI_CV_Tools.h b/source/module_ri/LRI_CV_Tools.h index 52085f6bc3..c97f502516 100644 --- a/source/module_ri/LRI_CV_Tools.h +++ b/source/module_ri/LRI_CV_Tools.h @@ -80,20 +80,25 @@ namespace LRI_CV_Tools std::map>> && ds_in); template - extern std::array cal_latvec_range(const double &rcut_times, const std::vector& orb_cutoff); + extern std::array cal_latvec_range(const double &rcut_times, + const UnitCell &ucell, + const std::vector& orb_cutoff); template extern std::map,RI::Tensor>>> get_CVws( + const UnitCell &ucell, const std::map>,RI::Tensor>> &CVs); template extern std::map,std::array,3>>>> get_dCVws( + const UnitCell &ucell, const std::array>,RI::Tensor>>,3> &dCVs); template extern std::array,RI::Tensor>>,3>,3> cal_dMRs( + const UnitCell &ucell, const std::array,RI::Tensor>>,3> &dMs); using TC = std::array; diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp index 3e720ba9f4..31f5a98aa3 100644 --- a/source/module_ri/LRI_CV_Tools.hpp +++ b/source/module_ri/LRI_CV_Tools.hpp @@ -245,14 +245,16 @@ LRI_CV_Tools::change_order(std::map>> template std::array -LRI_CV_Tools::cal_latvec_range(const double &rcut_times, const std::vector& orb_cutoff) +LRI_CV_Tools::cal_latvec_range(const double &rcut_times, + const UnitCell &ucell, + const std::vector& orb_cutoff) { double Rcut_max = 0; - for(int T=0; T proj = ModuleBase::Mathzone::latvec_projection( - std::array,3>{GlobalC::ucell.a1, GlobalC::ucell.a2, GlobalC::ucell.a3}); - const ModuleBase::Vector3 latvec_times = Rcut_max * rcut_times / (proj * GlobalC::ucell.lat0); + std::array,3>{ucell.a1, ucell.a2, ucell.a3}); + const ModuleBase::Vector3 latvec_times = Rcut_max * rcut_times / (proj * ucell.lat0); const ModuleBase::Vector3 latvec_times_ceil = {static_cast(std::ceil(latvec_times.x)), static_cast(std::ceil(latvec_times.y)), static_cast(std::ceil(latvec_times.z))}; @@ -263,23 +265,24 @@ LRI_CV_Tools::cal_latvec_range(const double &rcut_times, const std::vector std::map,RI::Tensor>>> LRI_CV_Tools::get_CVws( + const UnitCell &ucell, const std::map>,RI::Tensor>> &CVs) { std::map,RI::Tensor>>> CVws; for(const auto &CVs_A : CVs) { const TA iat0 = CVs_A.first; - const int it0 = GlobalC::ucell.iat2it[iat0]; - const int ia0 = GlobalC::ucell.iat2ia[iat0]; - const ModuleBase::Vector3 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; + const int it0 = ucell.iat2it[iat0]; + const int ia0 = ucell.iat2ia[iat0]; + const ModuleBase::Vector3 tau0 = ucell.atoms[it0].tau[ia0]; for(const auto &CVs_B : CVs_A.second) { const TA iat1 = CVs_B.first.first; - const int it1 = GlobalC::ucell.iat2it[iat1]; - const int ia1 = GlobalC::ucell.iat2ia[iat1]; + const int it1 = ucell.iat2it[iat1]; + const int ia1 = ucell.iat2ia[iat1]; const std::array &cell1 = CVs_B.first.second; - const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; - const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); + const ModuleBase::Vector3 tau1 = ucell.atoms[it1].tau[ia1]; + const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec); CVws[it0][it1][R_delta] = CVs_B.second; } } @@ -289,6 +292,7 @@ LRI_CV_Tools::get_CVws( template std::map,std::array,3>>>> LRI_CV_Tools::get_dCVws( + const UnitCell &ucell, const std::array>,RI::Tensor>>,3> &dCVs) { std::map,std::array,3>>>> dCVws; @@ -297,17 +301,17 @@ LRI_CV_Tools::get_dCVws( for(const auto &dCVs_A : dCVs[ix]) { const TA iat0 = dCVs_A.first; - const int it0 = GlobalC::ucell.iat2it[iat0]; - const int ia0 = GlobalC::ucell.iat2ia[iat0]; - const ModuleBase::Vector3 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; + const int it0 = ucell.iat2it[iat0]; + const int ia0 = ucell.iat2ia[iat0]; + const ModuleBase::Vector3 tau0 = ucell.atoms[it0].tau[ia0]; for(const auto &dCVs_B : dCVs_A.second) { const TA iat1 = dCVs_B.first.first; - const int it1 = GlobalC::ucell.iat2it[iat1]; - const int ia1 = GlobalC::ucell.iat2ia[iat1]; + const int it1 = ucell.iat2it[iat1]; + const int ia1 = ucell.iat2ia[iat1]; const std::array &cell1 = dCVs_B.first.second; - const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; - const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); + const ModuleBase::Vector3 tau1 = ucell.atoms[it1].tau[ia1]; + const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec); dCVws[it0][it1][R_delta][ix] = dCVs_B.second; } } @@ -320,6 +324,7 @@ LRI_CV_Tools::get_dCVws( template std::array,RI::Tensor>>,3>,3> LRI_CV_Tools::cal_dMRs( + const UnitCell &ucell, const std::array,RI::Tensor>>,3> &dMs) { auto get_R_delta = [&](const TA &iat0, const std::pair &A1) -> std::array diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 35dcb2184f..8169400ae6 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -40,7 +40,7 @@ void RPA_LRI::cal_rpa_cv(const UnitCell& ucell) } const std::array period = {p_kv->nmp[0], p_kv->nmp[1], p_kv->nmp[2]}; - 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, ucell,orb_cutoff_); const std::pair, std::vector>>>> list_As_Vs = RI::Distribute_Equally::distribute_atoms(this->mpi_comm, atoms, period_Vs, 2, false); @@ -52,7 +52,7 @@ void RPA_LRI::cal_rpa_cv(const UnitCell& ucell) }); this->Vs_period = RI::RI_Tools::cal_period(Vs, period); - 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, ucell,orb_cutoff_); const std::pair, std::vector>>>> list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false); From ea871f0e978ff5697296e4ad30a9b7fda9e268ba Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 20:28:40 +0800 Subject: [PATCH 13/24] change ucell in module_ri/ri_2d_comm.cpp --- .../operator_casida/operator_lr_exx.cpp | 12 ++++++++--- source/module_rdmft/rdmft_pot.cpp | 8 +++---- source/module_ri/Exx_LRI.hpp | 2 +- source/module_ri/Exx_LRI_interface.hpp | 8 +++---- source/module_ri/RI_2D_Comm.cpp | 6 +++--- source/module_ri/RI_2D_Comm.h | 11 +++++++--- source/module_ri/RI_2D_Comm.hpp | 21 ++++++++++++------- source/module_ri/RPA_LRI.hpp | 4 ++-- 8 files changed, 44 insertions(+), 28 deletions(-) diff --git a/source/module_lr/operator_casida/operator_lr_exx.cpp b/source/module_lr/operator_casida/operator_lr_exx.cpp index 68ccd403b2..a701e319c9 100644 --- a/source/module_lr/operator_casida/operator_lr_exx.cpp +++ b/source/module_lr/operator_casida/operator_lr_exx.cpp @@ -89,11 +89,17 @@ namespace LR } template - void OperatorLREXX::act(const int nbands, const int nbasis, const int npol, const T* psi_in, T* hpsi, const int ngk_ik, const bool is_first_node)const + void OperatorLREXX::act(const int nbands, + const int nbasis, + const int npol, + const T* psi_in, + T* hpsi, + const int ngk_ik, + const bool is_first_node)const { ModuleBase::TITLE("OperatorLREXX", "act"); // convert parallel info to LibRI interfaces - std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(this->pmat); + std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(ucell,this->pmat); // suppose Cs，Vs, have already been calculated in the ion-step of ground state // and DM_trans has been calculated in hPsi() outside. @@ -107,7 +113,7 @@ namespace LR // if multi-k, DM_trans(TR=double) -> Ds_trans(TR=T=complex) std::vector>>> Ds_trans = aims_nbasis.empty() ? - RI_2D_Comm::split_m2D_ktoR(this->kv, DMk_trans_pointer, this->pmat, 1) + RI_2D_Comm::split_m2D_ktoR(ucell,this->kv, DMk_trans_pointer, this->pmat, 1) : RI_Benchmark::split_Ds(DMk_trans_vector, aims_nbasis, ucell); //0.5 will be multiplied // LR_Util::print_CV(Ds_trans[0], "Ds_trans in OperatorLREXX", 1e-10); // 2. cal_Hs diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index 3f7f20868f..cd600cb933 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -274,8 +274,8 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) // transfer the DM_XC to appropriate format std::vector>,RI::Tensor>>> Ds_XC_d = std::is_same::value //gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR(*kv, DM_XC_pointer, *ParaV, nspin) - : RI_2D_Comm::split_m2D_ktoR(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); + ? RI_2D_Comm::split_m2D_ktoR(ucell,*kv, DM_XC_pointer, *ParaV, nspin) + : RI_2D_Comm::split_m2D_ktoR(ucell,*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); // provide the Ds_XC to Vxc_fromRI(V_exx_XC) if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) @@ -302,8 +302,8 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) // transfer the DM_XC to appropriate format std::vector>,RI::Tensor>>>> Ds_XC_c = std::is_same::value //gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR>(*kv, DM_XC_pointer, *ParaV, nspin) - : RI_2D_Comm::split_m2D_ktoR>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); + ? RI_2D_Comm::split_m2D_ktoR>(ucell,*kv, DM_XC_pointer, *ParaV, nspin) + : RI_2D_Comm::split_m2D_ktoR>(ucell,*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); // // provide the Ds_XC to Vxc_fromRI(V_exx_XC) if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 18cff2b2ee..fa0d41d3c7 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -198,7 +198,7 @@ void Exx_LRI::cal_exx_elec(const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(pv); + const std::vector, std::set>> judge = RI_2D_Comm::get_2D_judge(ucell,pv); this->Hexxs.resize(PARAM.inp.nspin); this->Eexx = 0; diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 8f159be632..cdf3b219c9 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -131,8 +131,8 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, else { this->mix_DMk_2D.mix(dm_in.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_in.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_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry); + ? RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin) + : RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm_in.get_paraV_pointer(), PARAM.inp.nspin, this->exx_spacegroup_symmetry); if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, ucell,*dm_in.get_paraV_pointer(), &this->symrot_); @@ -323,8 +323,8 @@ bool Exx_LRI_Interface::exx_after_converge( // 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); + ? RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), nspin) + : RI_2D_Comm::split_m2D_ktoR(ucell,*this->exx_ptr->p_kv, this->mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), nspin, this->exx_spacegroup_symmetry); // check the rotation of Ds // this->symrot_.test_HR_rotation(ucell.symm, ucell.atoms, ucell.st, 'D', Ds[0]); diff --git a/source/module_ri/RI_2D_Comm.cpp b/source/module_ri/RI_2D_Comm.cpp index db6e6d0244..9a156a7171 100644 --- a/source/module_ri/RI_2D_Comm.cpp +++ b/source/module_ri/RI_2D_Comm.cpp @@ -12,7 +12,7 @@ #include // judge[is] = {s0, s1} -auto RI_2D_Comm::get_2D_judge(const Parallel_2D& pv) +auto RI_2D_Comm::get_2D_judge(const UnitCell& ucell, const Parallel_2D& pv) -> std::vector, std::set>> { ModuleBase::TITLE("RI_2D_Comm","get_2D_judge"); @@ -24,7 +24,7 @@ auto RI_2D_Comm::get_2D_judge(const Parallel_2D& pv) { const int iwt0 = pv.local2global_row(iwt0_2D); int iat0=0;int iw0_b=0;int is0_b=0; - std::tie(iat0,iw0_b,is0_b) = RI_2D_Comm::get_iat_iw_is_block(iwt0); + std::tie(iat0,iw0_b,is0_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt0); iat0_list[is0_b].insert(iat0); } @@ -33,7 +33,7 @@ auto RI_2D_Comm::get_2D_judge(const Parallel_2D& pv) { const int iwt1 = pv.local2global_col(iwt1_2D); int iat1=0;int iw1_b=0;int is1_b=0; - std::tie(iat1,iw1_b,is1_b) = RI_2D_Comm::get_iat_iw_is_block(iwt1); + std::tie(iat1,iw1_b,is1_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt1); iat1_list[is1_b].insert(iat1); } diff --git a/source/module_ri/RI_2D_Comm.h b/source/module_ri/RI_2D_Comm.h index 56e4991472..159a0bf7dc 100644 --- a/source/module_ri/RI_2D_Comm.h +++ b/source/module_ri/RI_2D_Comm.h @@ -31,11 +31,16 @@ namespace RI_2D_Comm //public: template extern std::vector>>> - split_m2D_ktoR(const K_Vectors& kv, const std::vector& mks_2D, const Parallel_2D& pv, const int nspin, const bool spgsym = false); + split_m2D_ktoR(const UnitCell& ucell, + const K_Vectors& kv, + const std::vector& mks_2D, + const Parallel_2D& pv, + const int nspin, + const bool spgsym = false); // judge[is] = {s0, s1} extern std::vector, std::set>> - get_2D_judge(const Parallel_2D& pv); + get_2D_judge(const UnitCell& ucell, const Parallel_2D& pv); template extern void add_Hexx( @@ -72,7 +77,7 @@ namespace RI_2D_Comm //private: extern std::vector get_ik_list(const K_Vectors &kv, const int is_k); - extern inline std::tuple get_iat_iw_is_block(const int iwt); + extern inline std::tuple get_iat_iw_is_block(const UnitCell& ucell,const int& iwt); extern inline int get_is_block(const int is_k, const int is_row_b, const int is_col_b); extern inline std::tuple split_is_block(const int is_b); extern inline int get_iwt(const int iat, const int iw_b, const int is_b); diff --git a/source/module_ri/RI_2D_Comm.hpp b/source/module_ri/RI_2D_Comm.hpp index a168764340..ef8e6e8097 100644 --- a/source/module_ri/RI_2D_Comm.hpp +++ b/source/module_ri/RI_2D_Comm.hpp @@ -29,7 +29,12 @@ inline RI::Tensor> tensor_conj(const RI::Tensor -auto RI_2D_Comm::split_m2D_ktoR(const K_Vectors & kv, const std::vector&mks_2D, const Parallel_2D & pv, const int nspin, const bool spgsym) +auto RI_2D_Comm::split_m2D_ktoR(const UnitCell& ucell, + const K_Vectors & kv, + const std::vector&mks_2D, + const Parallel_2D & pv, + const int nspin, + const bool spgsym) -> std::vector>>> { ModuleBase::TITLE("RI_2D_Comm","split_m2D_ktoR"); @@ -63,7 +68,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const K_Vectors & kv, const std::vector mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik], pv.get_col_size(), pv.get_row_size()); const Tdata_m frac = SPIN_multiple * RI::Global_Func::convert(std::exp( - -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell) * GlobalC::ucell.latvec)))); + -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); if (static_cast(std::round(SPIN_multiple * kv.wk[ik] * kv.get_nkstot_full())) == 2) { set_mR_2D(mk_2D * (frac * 0.5) + tensor_conj(mk_2D * (frac * 0.5))); } else { set_mR_2D(mk_2D * frac); } @@ -75,7 +80,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const K_Vectors & kv, const std::vector mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik_full + is_k * kv.get_nkstot_full()], pv.get_col_size(), pv.get_row_size()); const Tdata_m frac = SPIN_multiple * RI::Global_Func::convert(std::exp( - -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * ((isym_kvd.second * GlobalC::ucell.G) * (RI_Util::array3_to_Vector3(cell) * GlobalC::ucell.latvec)))); + -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * ((isym_kvd.second * ucell.G) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); set_mR_2D(mk_2D * frac); ++ik_full; } @@ -87,7 +92,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const K_Vectors & kv, const std::vector &mR_a2D = mRs_a2D[is_b][iat0][{iat1,cell}]; if (mR_a2D.empty()) { mR_a2D = RI::Tensor( - {static_cast(GlobalC::ucell.atoms[it0].nw), + {static_cast(ucell.atoms[it0].nw), static_cast( - GlobalC::ucell.atoms[it1].nw)}); + ucell.atoms[it1].nw)}); } mR_a2D(iw0_b,iw1_b) = mR_2D(iwt0_2D, iwt1_2D); } @@ -165,7 +170,7 @@ void RI_2D_Comm::add_Hexx( } std::tuple -RI_2D_Comm::get_iat_iw_is_block(const int iwt) +RI_2D_Comm::get_iat_iw_is_block(const UnitCell& ucell,const int& iwt) { const int iat = GlobalC::ucell.iwt2iat[iwt]; const int iw = GlobalC::ucell.iwt2iw[iwt]; diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 8169400ae6..8c576725a9 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -100,8 +100,8 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix const std::vector>>> Ds = PARAM.globalv.gamma_only_local - ? RI_2D_Comm::split_m2D_ktoR(kv, mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin) - : RI_2D_Comm::split_m2D_ktoR(kv, mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin, exx_spacegroup_symmetry); + ? RI_2D_Comm::split_m2D_ktoR(ucell,kv, mix_DMk_2D.get_DMk_gamma_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin) + : RI_2D_Comm::split_m2D_ktoR(ucell,kv, mix_DMk_2D.get_DMk_k_out(), *dm.get_paraV_pointer(), PARAM.inp.nspin, exx_spacegroup_symmetry); // set parameters for bare Coulomb potential GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; From 0107eb5d5eb7d2df5e9909a3a7252eedcf8a5e6b Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 20:53:51 +0800 Subject: [PATCH 14/24] change ucell in module_ri/ri_2d_comm.cpp --- .../hamilt_lcaodft/LCAO_hamilt.hpp | 18 +++++------ .../hamilt_lcaodft/hamilt_lcao.cpp | 1 + .../operator_lcao/op_exx_lcao.h | 5 ++- .../operator_lcao/op_exx_lcao.hpp | 32 +++++++++++-------- .../hamilt_lcaodft/spar_exx.h | 4 +-- .../hamilt_lcaodft/spar_hsr.cpp | 9 ++++-- .../hamilt_lcaodft/spar_hsr.h | 3 +- source/module_io/output_mat_sparse.cpp | 2 +- source/module_io/write_HS_R.cpp | 9 +++--- source/module_io/write_HS_R.h | 3 +- source/module_io/write_vxc.hpp | 4 +-- source/module_rdmft/rdmft_pot.cpp | 2 ++ source/module_ri/RI_2D_Comm.h | 3 +- source/module_ri/RI_2D_Comm.hpp | 26 ++++++++------- 14 files changed, 71 insertions(+), 50 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp index 42317c2df7..10f508c578 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.hpp @@ -21,7 +21,7 @@ // Peize Lin add 2022.09.13 template -void sparse_format::cal_HR_exx( +void sparse_format::cal_HR_exx(const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, @@ -36,15 +36,15 @@ void sparse_format::cal_HR_exx( const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha; std::map> atoms_pos; - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) { + for (int iat = 0; iat < ucell.nat; ++iat) { atoms_pos[iat] = RI_Util::Vector3_to_array3( - GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]] - .tau[GlobalC::ucell.iat2ia[iat]]); + ucell.atoms[ucell.iat2it[iat]] + .tau[ucell.iat2ia[iat]]); } const std::array, 3> latvec - = {RI_Util::Vector3_to_array3(GlobalC::ucell.a1), // too bad to use GlobalC here, - RI_Util::Vector3_to_array3(GlobalC::ucell.a2), - RI_Util::Vector3_to_array3(GlobalC::ucell.a3)}; + = {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here, + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3)}; const std::array Rs_period = {nmp[0], nmp[1], nmp[2]}; @@ -84,7 +84,7 @@ void sparse_format::cal_HR_exx( for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) { - const int iwt0 = RI_2D_Comm::get_iwt(iat0, iw0, is0_b); + const int iwt0 = RI_2D_Comm::get_iwt(ucell,iat0, iw0, is0_b); const int iwt0_local = pv.global2local_row(iwt0); if (iwt0_local < 0) @@ -94,7 +94,7 @@ void sparse_format::cal_HR_exx( for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) { - const int iwt1 = RI_2D_Comm::get_iwt(iat1, iw1, is1_b); + const int iwt1 = RI_2D_Comm::get_iwt(ucell,iat1, iw1, is1_b); const int iwt1_local = pv.global2local_col(iwt1); if (iwt1_local < 0) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index f708ad1310..572c26e4d6 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -403,6 +403,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, // and calculate Cs, Vs Operator* exx = new OperatorEXX>(this->hsk, this->hR, + ucell, *this->kv, Hexxd, Hexxc, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h index 9e1653be1a..cc68a11b06 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h @@ -28,6 +28,7 @@ class OperatorEXX> : public OperatorLCAO public: OperatorEXX>(HS_Matrix_K* hsk_in, hamilt::HContainer* hR_in, + const UnitCell& ucell, const K_Vectors& kv_in, std::vector>>>* Hexxd_in = nullptr, std::vector>>>>* Hexxc_in = nullptr, @@ -58,7 +59,9 @@ class OperatorEXX> : public OperatorLCAO const int istep = 0; // the ion step void add_loaded_Hexx(const int ik); - + + const UnitCell& ucell; + const K_Vectors& kv; // if k points has no shift, use cell_nearest to reduce the memory cost diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp index 56dc6f1f4d..e20ca09b6e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp @@ -49,9 +49,9 @@ namespace hamilt auto* pv = hR->get_paraV(); auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period); bool need_allocate = false; - for (int iat0 = 0;iat0 < GlobalC::ucell.nat;++iat0) + for (int iat0 = 0;iat0 < nat;++iat0) { - for (int iat1 = 0;iat1 < GlobalC::ucell.nat;++iat1) + for (int iat1 = 0;iat1 < nat;++iat1) { // complete the atom pairs that has orbitals in this processor but not in hR due to the adj_list // but adj_list is not enought for EXX, which is more nonlocal than Nonlocal @@ -81,6 +81,7 @@ namespace hamilt template OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, HContainer*hR_in, + const UnitCell& ucell_in, const K_Vectors& kv_in, std::vector>>>* Hexxd_in, std::vector>>>>* Hexxc_in, @@ -89,6 +90,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, int* two_level_step_in, const bool restart_in) : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), + ucell(ucell_in), kv(kv_in), Hexxd(Hexxd_in), Hexxc(Hexxc_in), @@ -115,12 +117,12 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, // Read HexxR in CSR format if (GlobalC::exx_info.info_ri.real_number) { - ModuleIO::read_Hexxs_csr(file_name_exx, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); + ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxd, this->hR); } } else { - ModuleIO::read_Hexxs_csr(file_name_exx, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); + ModuleIO::read_Hexxs_csr(file_name_exx, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); if (this->add_hexx_type == Add_Hexx_Type::R) { reallocate_hcontainer(*Hexxc, this->hR); } } } @@ -154,19 +156,19 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, { // set cell_nearest std::map> atoms_pos; - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) { + for (int iat = 0; iat < ucell.nat; ++iat) { atoms_pos[iat] = RI_Util::Vector3_to_array3( - GlobalC::ucell.atoms[GlobalC::ucell.iat2it[iat]] - .tau[GlobalC::ucell.iat2ia[iat]]); + ucell.atoms[ucell.iat2it[iat]] + .tau[ucell.iat2ia[iat]]); } const std::array, 3> latvec - = { RI_Util::Vector3_to_array3(GlobalC::ucell.a1), - RI_Util::Vector3_to_array3(GlobalC::ucell.a2), - RI_Util::Vector3_to_array3(GlobalC::ucell.a3) }; + = { RI_Util::Vector3_to_array3(ucell.a1), + RI_Util::Vector3_to_array3(ucell.a2), + RI_Util::Vector3_to_array3(ucell.a3) }; this->cell_nearest.init(atoms_pos, latvec, Rs_period); - reallocate_hcontainer(GlobalC::ucell.nat, this->hR, Rs_period, &this->cell_nearest); + reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest); } - else { reallocate_hcontainer(GlobalC::ucell.nat, this->hR, Rs_period); } + else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); } } if (this->restart) @@ -216,10 +218,10 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, { // Read HexxR in CSR format if (GlobalC::exx_info.info_ri.real_number) { - ModuleIO::read_Hexxs_csr(restart_HR_path, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); + ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd); } else { - ModuleIO::read_Hexxs_csr(restart_HR_path, GlobalC::ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); + ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc); } } else @@ -317,6 +319,7 @@ void OperatorEXX>::contributeHk(int ik) if (GlobalC::exx_info.info_ri.real_number) { RI_2D_Comm::add_Hexx( + ucell, this->kv, ik, GlobalC::exx_info.info_global.hybrid_alpha, @@ -325,6 +328,7 @@ void OperatorEXX>::contributeHk(int ik) this->hsk->get_hk()); } else { RI_2D_Comm::add_Hexx( + ucell, this->kv, ik, GlobalC::exx_info.info_global.hybrid_alpha, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h index a9b7b7156f..5c48b4cd63 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_exx.h @@ -11,12 +11,12 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp" #include "module_basis/module_ao/parallel_orbitals.h" - +#include "module_cell/unitcell.h" namespace sparse_format { template -void cal_HR_exx( +void cal_HR_exx(const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index d845c3c8f2..8fa45e7a6e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -7,7 +7,8 @@ #include "spar_exx.h" #include "spar_u.h" -void sparse_format::cal_HSR(const Parallel_Orbitals& pv, +void sparse_format::cal_HSR(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const int& current_spin, @@ -99,7 +100,8 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, if (GlobalC::exx_info.info_global.cal_exx) { if (Hexxd && GlobalC::exx_info.info_ri.real_number) { - sparse_format::cal_HR_exx(pv, + sparse_format::cal_HR_exx(ucell, + pv, HS_Arrays, current_spin, sparse_thr, @@ -108,7 +110,8 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, } else if (Hexxc && !GlobalC::exx_info.info_ri.real_number) { - sparse_format::cal_HR_exx(pv, + sparse_format::cal_HR_exx(ucell, + pv, HS_Arrays, current_spin, sparse_thr, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index 497586574e..a1118f994e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -5,7 +5,8 @@ namespace sparse_format { using TAC = std::pair>; - void cal_HSR(const Parallel_Orbitals& pv, + void cal_HSR(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const int& current_spin, diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index 11cd2caf81..d3a6f6e1b1 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -45,7 +45,7 @@ void output_mat_sparse(const bool& out_mat_hsR, //! generate a file containing the Hamiltonian and S(overlap) matrices if (out_mat_hsR) { - output_HSR(istep, v_eff, pv, HS_Arrays, grid, kv, p_ham); + output_HSR(ucell,istep, v_eff, pv, HS_Arrays, grid, kv, p_ham); } //! generate a file containing the kinetic energy matrix diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 5b59f2ee24..686787af00 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -12,7 +12,8 @@ // The 'sparse_thr' is the accuracy of the sparse matrix. // If the absolute value of the matrix element is less than or equal to the // 'sparse_thr', it will be ignored. -void ModuleIO::output_HSR(const int& istep, +void ModuleIO::output_HSR(const UnitCell& ucell, + const int& istep, const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, @@ -37,7 +38,7 @@ void ModuleIO::output_HSR(const int& istep, if (nspin == 1 || nspin == 4) { const int spin_now = 0; // jingan add 2021-6-4, modify 2021-12-2 - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif @@ -47,7 +48,7 @@ void ModuleIO::output_HSR(const int& istep, int spin_now = 1; // save HR of spin down first (the current spin always be down) - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif @@ -61,7 +62,7 @@ void ModuleIO::output_HSR(const int& istep, spin_now = 0; } - sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham + sparse_format::cal_HSR(ucell,pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham #ifdef __EXX , Hexxd, Hexxc #endif diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h index 6c9bd46768..df4c0251a5 100644 --- a/source/module_io/write_HS_R.h +++ b/source/module_io/write_HS_R.h @@ -11,7 +11,8 @@ namespace ModuleIO { using TAC = std::pair>; -void output_HSR(const int& istep, +void output_HSR(const UnitCell& ucell, + const int& istep, const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, diff --git a/source/module_io/write_vxc.hpp b/source/module_io/write_vxc.hpp index d6fad87dad..fcd03e07a3 100644 --- a/source/module_io/write_vxc.hpp +++ b/source/module_io/write_vxc.hpp @@ -246,10 +246,10 @@ void write_Vxc(const int nspin, std::vector> e_orb_tot; // orbital energy (total) #ifdef __EXX hamilt::OperatorEXX> vexx_op_ao(&vxc_k_ao, - &vxcs_R_ao[0] /*for paraV*/, kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); + &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); hamilt::HS_Matrix_K vexxonly_k_ao(pv, 1); // only hk is needed, sk is skipped hamilt::OperatorEXX> vexxonly_op_ao(&vexxonly_k_ao, - &vxcs_R_ao[0]/*for paraV*/, kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); + &vxcs_R_ao[0],ucell,/*for paraV*/ kv, Hexxd, Hexxc, hamilt::Add_Hexx_Type::k); std::vector> e_orb_exx; // orbital energy (EXX) #endif hamilt::OperatorDFTU> vdftu_op_ao(&vxc_k_ao, kv.kvec_d, nullptr, kv.isk); diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index cd600cb933..f405364abe 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -291,6 +291,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) V_exx_XC = new hamilt::OperatorEXX>( hsk_exx_XC, HR_exx_XC, + ucell, *kv, &Vxc_fromRI_d->Hexxs, nullptr, @@ -319,6 +320,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) V_exx_XC = new hamilt::OperatorEXX>( hsk_exx_XC, HR_exx_XC, + ucell, *kv, nullptr, &Vxc_fromRI_c->Hexxs, diff --git a/source/module_ri/RI_2D_Comm.h b/source/module_ri/RI_2D_Comm.h index 159a0bf7dc..c7d023e735 100644 --- a/source/module_ri/RI_2D_Comm.h +++ b/source/module_ri/RI_2D_Comm.h @@ -44,6 +44,7 @@ namespace RI_2D_Comm template extern void add_Hexx( + const UnitCell& ucell, const K_Vectors& kv, const int ik, const double alpha, @@ -80,7 +81,7 @@ namespace RI_2D_Comm extern inline std::tuple get_iat_iw_is_block(const UnitCell& ucell,const int& iwt); extern inline int get_is_block(const int is_k, const int is_row_b, const int is_col_b); extern inline std::tuple split_is_block(const int is_b); - extern inline int get_iwt(const int iat, const int iw_b, const int is_b); + extern inline int get_iwt(const UnitCell& ucell, const int iat, const int iw_b, const int is_b); } #include "RI_2D_Comm.hpp" diff --git a/source/module_ri/RI_2D_Comm.hpp b/source/module_ri/RI_2D_Comm.hpp index ef8e6e8097..055ea914ef 100644 --- a/source/module_ri/RI_2D_Comm.hpp +++ b/source/module_ri/RI_2D_Comm.hpp @@ -93,7 +93,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const UnitCell& ucell, : pv.local2global_row(iwt0_2D); int iat0, iw0_b, is0_b; std::tie(iat0,iw0_b,is0_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt0); - const int it0 = GlobalC::ucell.iat2it[iat0]; + const int it0 = ucell.iat2it[iat0]; for(int iwt1_2D=0; iwt1_2D!=mR_2D.shape[1]; ++iwt1_2D) { const int iwt1 =ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver) @@ -101,7 +101,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const UnitCell& ucell, : pv.local2global_col(iwt1_2D); int iat1, iw1_b, is1_b; std::tie(iat1,iw1_b,is1_b) = RI_2D_Comm::get_iat_iw_is_block(ucell,iwt1); - const int it1 = GlobalC::ucell.iat2it[iat1]; + const int it1 = ucell.iat2it[iat1]; const int is_b = RI_2D_Comm::get_is_block(is_k, is0_b, is1_b); RI::Tensor &mR_a2D = mRs_a2D[is_b][iat0][{iat1,cell}]; @@ -123,6 +123,7 @@ auto RI_2D_Comm::split_m2D_ktoR(const UnitCell& ucell, template void RI_2D_Comm::add_Hexx( + const UnitCell &ucell, const K_Vectors &kv, const int ik, const double alpha, @@ -146,17 +147,17 @@ void RI_2D_Comm::add_Hexx( const TA &iat1 = Hs_tmpB.first.first; const TC &cell1 = Hs_tmpB.first.second; const std::complex frac = alpha - * std::exp( ModuleBase::TWO_PI*ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec)) ); + * std::exp( ModuleBase::TWO_PI*ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell1)*ucell.latvec)) ); const RI::Tensor &H = Hs_tmpB.second; for(size_t iw0_b=0; iw0_b RI_2D_Comm::get_iat_iw_is_block(const UnitCell& ucell,const int& iwt) { - const int iat = GlobalC::ucell.iwt2iat[iwt]; - const int iw = GlobalC::ucell.iwt2iw[iwt]; + const int iat = ucell.iwt2iat[iwt]; + const int iw = ucell.iwt2iw[iwt]; switch(PARAM.inp.nspin) { case 1: case 2: @@ -212,10 +213,13 @@ RI_2D_Comm::split_is_block(const int is_b) -int RI_2D_Comm::get_iwt(const int iat, const int iw_b, const int is_b) +int RI_2D_Comm::get_iwt(const UnitCell& ucell, + const int iat, + const int iw_b, + const int is_b) { - const int it = GlobalC::ucell.iat2it[iat]; - const int ia = GlobalC::ucell.iat2ia[iat]; + const int it = ucell.iat2it[iat]; + const int ia = ucell.iat2ia[iat]; int iw=-1; switch(PARAM.inp.nspin) { @@ -226,7 +230,7 @@ int RI_2D_Comm::get_iwt(const int iat, const int iw_b, const int is_b) default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__)); } - const int iwt = GlobalC::ucell.itiaiw2iwt(it,ia,iw); + const int iwt = ucell.itiaiw2iwt(it,ia,iw); return iwt; } From d091cad6353942e8f6f34fce7cf02a5d7f3f7054 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 21:17:26 +0800 Subject: [PATCH 15/24] change ucell in module_ri/rpa_lri.cpp --- source/module_ri/RPA_LRI.h | 6 +++--- source/module_ri/RPA_LRI.hpp | 40 ++++++++++++++++++------------------ 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h index 8b27c76bfe..5319fcf76f 100644 --- a/source/module_ri/RPA_LRI.h +++ b/source/module_ri/RPA_LRI.h @@ -49,11 +49,11 @@ template class RPA_LRI const psi::Psi& psi, const elecstate::ElecState* pelec); void out_eigen_vector(const Parallel_Orbitals& parav, const psi::Psi& psi); - void out_struc(); + void out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G); void out_bands(const elecstate::ElecState *pelec); - void out_Cs(); - void out_coulomb_k(); + void out_Cs(const UnitCell &ucell); + void out_coulomb_k(const UnitCell &ucell); // void print_matrix(char *desc, const ModuleBase::matrix &mat); // void print_complex_matrix(char *desc, const ModuleBase::ComplexMatrix &mat); // void init(const MPI_Comm &mpi_comm_in); diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp index 8c576725a9..79c851fa3a 100644 --- a/source/module_ri/RPA_LRI.hpp +++ b/source/module_ri/RPA_LRI.hpp @@ -128,14 +128,14 @@ void RPA_LRI::out_for_RPA(const UnitCell& ucell, ModuleBase::TITLE("DFT_RPA_interface", "out_for_RPA"); this->out_bands(pelec); this->out_eigen_vector(parav, psi); - this->out_struc(); + this->out_struc(ucell.latvec, ucell.G); this->cal_rpa_cv(ucell); std::cout << "rpa_pca_threshold: " << this->info.pca_threshold << std::endl; std::cout << "rpa_ccp_rmesh_times: " << this->info.ccp_rmesh_times << std::endl; std::cout << "rpa_lcao_exx(Ha): " << std::fixed << std::setprecision(15) << exx_lri_rpa.Eexx / 2.0 << std::endl; - this->out_Cs(); - this->out_coulomb_k(); + this->out_Cs(ucell); + this->out_coulomb_k(ucell); std::cout << "etxc(Ha): " << std::fixed << std::setprecision(15) << pelec->f_en.etxc / 2.0 << std::endl; std::cout << "etot(Ha): " << std::fixed << std::setprecision(15) << pelec->f_en.etot / 2.0 << std::endl; @@ -212,7 +212,7 @@ void RPA_LRI::out_eigen_vector(const Parallel_Orbitals& parav, const p } template -void RPA_LRI::out_struc() +void RPA_LRI::out_struc(const ModuleBase::Matrix3& latvec, const ModuleBase::Matrix3& G) { if (GlobalV::MY_RANK != 0) { @@ -221,8 +221,8 @@ void RPA_LRI::out_struc() ModuleBase::TITLE("DFT_RPA_interface", "out_struc"); double TWOPI_Bohr2A = ModuleBase::TWO_PI * ModuleBase::BOHR_TO_A; const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks(); - ModuleBase::Matrix3 lat = GlobalC::ucell.latvec / ModuleBase::BOHR_TO_A; - ModuleBase::Matrix3 G = GlobalC::ucell.G * TWOPI_Bohr2A; + ModuleBase::Matrix3 lat = latvec / ModuleBase::BOHR_TO_A; + ModuleBase::Matrix3 G_RPA = G * TWOPI_Bohr2A; std::stringstream ss; ss << "stru_out"; std::ofstream ofs; @@ -231,9 +231,9 @@ void RPA_LRI::out_struc() ofs << lat.e21 << std::setw(15) << lat.e22 << std::setw(15) << lat.e23 << std::endl; ofs << lat.e31 << std::setw(15) << lat.e32 << std::setw(15) << lat.e33 << std::endl; - ofs << G.e11 << std::setw(15) << G.e12 << std::setw(15) << G.e13 << std::endl; - ofs << G.e21 << std::setw(15) << G.e22 << std::setw(15) << G.e23 << std::endl; - ofs << G.e31 << std::setw(15) << G.e32 << std::setw(15) << G.e33 << std::endl; + ofs << G_RPA.e11 << std::setw(15) << G_RPA.e12 << std::setw(15) << G_RPA.e13 << std::endl; + ofs << G_RPA.e21 << std::setw(15) << G_RPA.e22 << std::setw(15) << G_RPA.e23 << std::endl; + ofs << G_RPA.e31 << std::setw(15) << G_RPA.e32 << std::setw(15) << G_RPA.e33 << std::endl; ofs << p_kv->nmp[0] << std::setw(6) << p_kv->nmp[1] << std::setw(6) << p_kv->nmp[2] << std::setw(6) << std::endl; @@ -286,23 +286,23 @@ void RPA_LRI::out_bands(const elecstate::ElecState* pelec) } template -void RPA_LRI::out_Cs() +void RPA_LRI::out_Cs(const UnitCell& ucell) { std::stringstream ss; ss << "Cs_data_" << GlobalV::MY_RANK << ".txt"; std::ofstream ofs; ofs.open(ss.str().c_str(), std::ios::out); - ofs << GlobalC::ucell.nat << " " << 0 << std::endl; + ofs << ucell.nat << " " << 0 << std::endl; for (auto& Ip: this->Cs_period) { size_t I = Ip.first; - size_t i_num = GlobalC::ucell.atoms[GlobalC::ucell.iat2it[I]].nw; + size_t i_num = ucell.atoms[ucell.iat2it[I]].nw; for (auto& JPp: Ip.second) { size_t J = JPp.first.first; auto R = JPp.first.second; auto& tmp_Cs = JPp.second; - size_t j_num = GlobalC::ucell.atoms[GlobalC::ucell.iat2it[J]].nw; + size_t j_num = ucell.atoms[ucell.iat2it[J]].nw; ofs << I + 1 << " " << J + 1 << " " << R[0] << " " << R[1] << " " << R[2] << " " << i_num << std::endl; @@ -324,14 +324,14 @@ void RPA_LRI::out_Cs() } template -void RPA_LRI::out_coulomb_k() +void RPA_LRI::out_coulomb_k(const UnitCell &ucell) { int all_mu = 0; - vector mu_shift(GlobalC::ucell.nat); - for (int I = 0; I != GlobalC::ucell.nat; I++) + vector mu_shift(ucell.nat); + for (int I = 0; I != ucell.nat; I++) { mu_shift[I] = all_mu; - all_mu += exx_lri_rpa.cv.get_index_abfs_size(GlobalC::ucell.iat2it[I]); + all_mu += exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[I]); } const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks(); std::stringstream ss; @@ -344,7 +344,7 @@ void RPA_LRI::out_coulomb_k() for (auto& Ip: this->Vs_period) { auto I = Ip.first; - size_t mu_num = exx_lri_rpa.cv.get_index_abfs_size(GlobalC::ucell.iat2it[I]); + size_t mu_num = exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[I]); for (int ik = 0; ik != nks_tot; ik++) { @@ -360,7 +360,7 @@ void RPA_LRI::out_coulomb_k() } RI::Tensor> tmp_VR = RI::Global_Func::convert>(JPp.second); - const double arg = 1 * (p_kv->kvec_c[ik] * (RI_Util::array3_to_Vector3(R) * GlobalC::ucell.latvec)) + const double arg = 1 * (p_kv->kvec_c[ik] * (RI_Util::array3_to_Vector3(R) * ucell.latvec)) * ModuleBase::TWO_PI; // latvec const std::complex kphase = std::complex(cos(arg), sin(arg)); if (Vq_k_IJ[J].empty()) @@ -373,7 +373,7 @@ void RPA_LRI::out_coulomb_k() { auto iJ = vq_Jp.first; auto& vq_J = vq_Jp.second; - size_t nu_num = exx_lri_rpa.cv.get_index_abfs_size(GlobalC::ucell.iat2it[iJ]); + size_t nu_num = exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[iJ]); ofs << all_mu << " " << mu_shift[I] + 1 << " " << mu_shift[I] + mu_num << " " << mu_shift[iJ] + 1 << " " << mu_shift[iJ] + nu_num << std::endl; ofs << ik + 1 << " " << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl; From d7a453cfd1082aba72cfdd5da4dcbe5a41345425 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 21:23:16 +0800 Subject: [PATCH 16/24] change ucell in module_ri/lri_cv_tools.cpp --- source/module_ri/LRI_CV_Tools.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp index 31f5a98aa3..39087e6641 100644 --- a/source/module_ri/LRI_CV_Tools.hpp +++ b/source/module_ri/LRI_CV_Tools.hpp @@ -331,13 +331,13 @@ LRI_CV_Tools::cal_dMRs( { const TA iat1 = A1.first; const TC &cell1 = A1.second; - const int it0 = GlobalC::ucell.iat2it[iat0]; - const int ia0 = GlobalC::ucell.iat2ia[iat0]; - const int it1 = GlobalC::ucell.iat2it[iat1]; - const int ia1 = GlobalC::ucell.iat2ia[iat1]; - const ModuleBase::Vector3 tau0 = GlobalC::ucell.atoms[it0].tau[ia0]; - const ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[it1].tau[ia1]; - const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*GlobalC::ucell.latvec); + const int it0 = ucell.iat2it[iat0]; + const int ia0 = ucell.iat2ia[iat0]; + const int it1 = ucell.iat2it[iat1]; + const int ia1 = ucell.iat2ia[iat1]; + const ModuleBase::Vector3 tau0 = ucell.atoms[it0].tau[ia0]; + const ModuleBase::Vector3 tau1 = ucell.atoms[it1].tau[ia1]; + const Abfs::Vector3_Order R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec); return std::array{R_delta.x, R_delta.y, R_delta.z}; }; constexpr int Npos = 3; From 5564184e709b24068a957dd7aad8acfc47388ceb Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 21:36:31 +0800 Subject: [PATCH 17/24] change ucell in module_ri/matirx_orb11.cpp --- source/module_ri/Exx_LRI.hpp | 1 + source/module_ri/LRI_CV.h | 1 + source/module_ri/LRI_CV.hpp | 3 +- source/module_ri/Matrix_Orbs11.cpp | 71 ++++++++++++++++-------------- source/module_ri/Matrix_Orbs11.h | 6 ++- source/module_ri/exx_opt_orb.cpp | 12 ++--- 6 files changed, 52 insertions(+), 42 deletions(-) diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index fa0d41d3c7..0ccaf1fe14 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -97,6 +97,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, } this->cv.set_orbitals( + ucell, orb, this->lcaos, this->abfs, this->abfs_ccp, this->info.kmesh_times, this->info.ccp_rmesh_times ); diff --git a/source/module_ri/LRI_CV.h b/source/module_ri/LRI_CV.h index 17df63ff0e..99b52476ea 100644 --- a/source/module_ri/LRI_CV.h +++ b/source/module_ri/LRI_CV.h @@ -34,6 +34,7 @@ class LRI_CV ~LRI_CV(); void set_orbitals( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &lcaos_in, const std::vector>> &abfs_in, diff --git a/source/module_ri/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp index a6c764122a..0d52d16853 100644 --- a/source/module_ri/LRI_CV.hpp +++ b/source/module_ri/LRI_CV.hpp @@ -37,6 +37,7 @@ LRI_CV::~LRI_CV() template void LRI_CV::set_orbitals( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &lcaos_in, const std::vector>> &abfs_in, @@ -61,7 +62,7 @@ void LRI_CV::set_orbitals( range_abfs = Exx_Abfs::Abfs_Index::construct_range( abfs ); this->index_abfs = ModuleBase::Element_Basis_Index::construct_index( range_abfs ); - this->m_abfs_abfs.init( 2, orb, kmesh_times, (1+this->ccp_rmesh_times)/2.0 ); + this->m_abfs_abfs.init( 2, ucell,orb, kmesh_times, (1+this->ccp_rmesh_times)/2.0 ); this->m_abfs_abfs.init_radial( this->abfs_ccp, this->abfs ); this->m_abfs_abfs.init_radial_table(); diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index 76aab8a6e4..ad01a78a6c 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -9,7 +9,11 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs11::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs11::init(const int mode, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs11", "init"); ModuleBase::timer::tick("Matrix_Orbs11", "init"); @@ -21,7 +25,7 @@ void Matrix_Orbs11::init(const int mode, const LCAO_Orbitals& orb, const double for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -64,13 +68,13 @@ 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 +93,13 @@ 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,17 +113,18 @@ 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"); } -void Matrix_Orbs11::init_radial_table(const std::map>>& Rs) +void Matrix_Orbs11::init_radial_table(const double lat0, + const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs11", "init_radial_table_Rs"); ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); @@ -133,23 +138,23 @@ void Matrix_Orbs11::init_radial_table(const std::map radials; for (const double& R: RsB.second) { - const double position = R * GlobalC::ucell.lat0 / lcao_dr_; + const double position = R * lat0 / lcao_dr_; const size_t iq = static_cast(position); for (size_t i = 0; i != 4; ++i) { - radials.insert(iq + 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.h b/source/module_ri/Matrix_Orbs11.h index 8861198029..a2afe85352 100644 --- a/source/module_ri/Matrix_Orbs11.h +++ b/source/module_ri/Matrix_Orbs11.h @@ -12,7 +12,7 @@ #include "module_basis/module_ao/ORB_gaunt_table.h" #include "module_basis/module_ao/ORB_read.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb11.h" - +#include "module_cell/unitcell.h" #include #include #include @@ -25,6 +25,7 @@ class Matrix_Orbs11 // 1: // 2: void init(const int mode, + const UnitCell& ucell, const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR @@ -34,7 +35,8 @@ class Matrix_Orbs11 void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B); void init_radial_table(); - void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const double lat0, + const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 0bc9dfd231..9b8382f871 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -105,10 +105,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co 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( 2,ucell,orb, this->kmesh_times, 1 ); m_jys_jys.init_radial( jle.jle, jle.jle ); #if TEST_EXX_RADIAL>=1 - m_jys_jys.init_radial_table(radial_R); + m_jys_jys.init_radial_table(ucell.lat0,radial_R); #else m_jys_jys.init_radial_table(); #endif @@ -121,10 +121,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co 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, ucell, orb, this->kmesh_times, 1 ); m_abfs_abfs.init_radial( abfs, abfs ); #if TEST_EXX_RADIAL>=1 - m_abfs_abfs.init_radial_table(radial_R); + m_abfs_abfs.init_radial_table(ucell.lat0,radial_R); #else m_abfs_abfs.init_radial_table(); #endif @@ -153,10 +153,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co 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( 2, ucell,orb, this->kmesh_times, 1 ); m_jys_abfs.init_radial( jle.jle, abfs ); #if TEST_EXX_RADIAL>=1 - m_jys_abfs.init_radial_table(radial_R); + m_jys_abfs.init_radial_table(ucell.lat0,radial_R); #else m_jys_abfs.init_radial_table(); #endif From 00b3046d7ddbb6f38ac03ab043cf7af127f0f626 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 21:49:15 +0800 Subject: [PATCH 18/24] change ucell in module_ri/matirx_orb21.cpp --- source/module_ri/ABFs_Construct-PCA.cpp | 5 +++-- source/module_ri/ABFs_Construct-PCA.h | 2 ++ source/module_ri/Exx_LRI.hpp | 2 +- source/module_ri/LRI_CV.hpp | 2 +- source/module_ri/Matrix_Orbs21.cpp | 13 +++++++++---- source/module_ri/Matrix_Orbs21.h | 6 ++++-- source/module_ri/exx_abfs-construct_orbs.cpp | 6 ++++-- source/module_ri/exx_abfs-construct_orbs.h | 2 ++ source/module_ri/exx_opt_orb.cpp | 10 +++++----- 9 files changed, 31 insertions(+), 17 deletions(-) diff --git a/source/module_ri/ABFs_Construct-PCA.cpp b/source/module_ri/ABFs_Construct-PCA.cpp index 31b1731bce..d5fbfaf477 100644 --- a/source/module_ri/ABFs_Construct-PCA.cpp +++ b/source/module_ri/ABFs_Construct-PCA.cpp @@ -72,6 +72,7 @@ namespace PCA } std::vector, RI::Tensor>>> cal_PCA( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &lcaos, const std::vector>> &abfs, @@ -96,14 +97,14 @@ namespace PCA } Matrix_Orbs21 m_abfslcaos_lcaos; - m_abfslcaos_lcaos.init( 1, orb, kmesh_times, 1 ); + m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); std::map>> delta_R; for( std::size_t it=0; it!=abfs.size(); ++it ) { delta_R[it][it] = {0.0}; } - m_abfslcaos_lcaos.init_radial_table(delta_R); + m_abfslcaos_lcaos.init_radial_table(ucell.lat0,delta_R); GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak; diff --git a/source/module_ri/ABFs_Construct-PCA.h b/source/module_ri/ABFs_Construct-PCA.h index bbd76cb61e..af5390ab9f 100644 --- a/source/module_ri/ABFs_Construct-PCA.h +++ b/source/module_ri/ABFs_Construct-PCA.h @@ -2,6 +2,7 @@ #define ABFS_CONSTRUCT_PCA_H #include "../module_basis/module_ao/ORB_read.h" +#include "module_cell/unitcell.h" #include #include @@ -15,6 +16,7 @@ namespace ABFs_Construct namespace PCA { extern std::vector,RI::Tensor>>> cal_PCA( + const UnitCell& ucell, const LCAO_Orbitals &orb, const std::vector>> &lcaos, const std::vector>> &abfs, // abfs must be orthonormal diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index 0ccaf1fe14..b2d85a4a09 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -62,7 +62,7 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, // #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 ); + abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell, orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold ); if(this->info.files_abfs.empty()) { this->abfs = abfs_same_atom; } else { diff --git a/source/module_ri/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp index 0d52d16853..79175d0d26 100644 --- a/source/module_ri/LRI_CV.hpp +++ b/source/module_ri/LRI_CV.hpp @@ -66,7 +66,7 @@ void LRI_CV::set_orbitals( this->m_abfs_abfs.init_radial( this->abfs_ccp, this->abfs ); this->m_abfs_abfs.init_radial_table(); - this->m_abfslcaos_lcaos.init( 1, orb, kmesh_times, 1 ); + this->m_abfslcaos_lcaos.init( 1, ucell , orb, kmesh_times, 1 ); this->m_abfslcaos_lcaos.init_radial( this->abfs_ccp, this->lcaos, this->lcaos ); this->m_abfslcaos_lcaos.init_radial_table(); diff --git a/source/module_ri/Matrix_Orbs21.cpp b/source/module_ri/Matrix_Orbs21.cpp index a772598ce8..a57c38410d 100644 --- a/source/module_ri/Matrix_Orbs21.cpp +++ b/source/module_ri/Matrix_Orbs21.cpp @@ -9,7 +9,11 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs21::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs21::init(const int mode, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs21", "init"); ModuleBase::timer::tick("Matrix_Orbs21", "init"); @@ -20,7 +24,7 @@ void Matrix_Orbs21::init(const int mode, const LCAO_Orbitals& orb, const double for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -116,7 +120,8 @@ void Matrix_Orbs21::init_radial_table() ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } -void Matrix_Orbs21::init_radial_table(const std::map>>& Rs) +void Matrix_Orbs21::init_radial_table(const double lat0, + const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs21", "init_radial_table_Rs"); ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); @@ -133,7 +138,7 @@ void Matrix_Orbs21::init_radial_table(const std::map radials; for (const double& R: RsB.second) { - const double position = R * GlobalC::ucell.lat0 / lcao_dr_; + const double position = R * lat0 / lcao_dr_; const size_t iq = static_cast(position); for (size_t i = 0; i != 4; ++i) radials.insert(iq + i); diff --git a/source/module_ri/Matrix_Orbs21.h b/source/module_ri/Matrix_Orbs21.h index 680b11b0da..57d85bf1c1 100644 --- a/source/module_ri/Matrix_Orbs21.h +++ b/source/module_ri/Matrix_Orbs21.h @@ -11,7 +11,7 @@ #include "module_basis/module_ao/ORB_gaunt_table.h" #include "module_basis/module_ao/ORB_read.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb21.h" - +#include "module_cell/unitcell.h" #include #include #include @@ -23,6 +23,7 @@ class Matrix_Orbs21 // mode: // 1: void init(const int mode, + const UnitCell& ucell, const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR @@ -35,7 +36,8 @@ class Matrix_Orbs21 const LCAO_Orbitals& orb_B); void init_radial_table(); - void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const double lat0, + const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { diff --git a/source/module_ri/exx_abfs-construct_orbs.cpp b/source/module_ri/exx_abfs-construct_orbs.cpp index ddf410a581..31d9d3a7a7 100644 --- a/source/module_ri/exx_abfs-construct_orbs.cpp +++ b/source/module_ri/exx_abfs-construct_orbs.cpp @@ -88,6 +88,7 @@ std::vector>> Exx_Abfs::Construct_ // P = f * Y std::vector>> Exx_Abfs::Construct_Orbs::abfs_same_atom( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &orbs, const double kmesh_times_mot, @@ -111,7 +112,7 @@ std::vector>> Exx_Abfs::Construct_ #endif const std::vector>>> - abfs_same_atom_pca_psi = pca( orb, abfs_same_atom, orbs, kmesh_times_mot, times_threshold ); + abfs_same_atom_pca_psi = pca(ucell,orb, abfs_same_atom, orbs, kmesh_times_mot, times_threshold ); #if TEST_EXX_LCAO==1 print_orbs(abfs_same_atom_pca_psi,"abfs_same_atom_pca_psi.dat"); @@ -267,6 +268,7 @@ std::vector>>> Exx_Abfs::Construct_O } std::vector>>> Exx_Abfs::Construct_Orbs::pca( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &abfs, const std::vector>> &orbs, @@ -277,7 +279,7 @@ std::vector>>> Exx_Abfs::Construct_O return std::vector>>>(abfs.size()); const std::vector,RI::Tensor>>> - eig = ABFs_Construct::PCA::cal_PCA( orb, orbs, abfs, kmesh_times_mot ); + eig = ABFs_Construct::PCA::cal_PCA(ucell, orb, orbs, abfs, kmesh_times_mot ); const std::vector>>> psis = get_psi( abfs ); std::vector>>> psis_new( psis.size() ); diff --git a/source/module_ri/exx_abfs-construct_orbs.h b/source/module_ri/exx_abfs-construct_orbs.h index 1d805a362f..cd33b19d3d 100644 --- a/source/module_ri/exx_abfs-construct_orbs.h +++ b/source/module_ri/exx_abfs-construct_orbs.h @@ -20,6 +20,7 @@ class Exx_Abfs::Construct_Orbs const double kmesh_times ); static std::vector>> abfs_same_atom( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &lcaos, const double kmesh_times_mot, @@ -43,6 +44,7 @@ class Exx_Abfs::Construct_Orbs const double norm_threshold = std::numeric_limits::min() ); static std::vector>>> pca( + const UnitCell &ucell, const LCAO_Orbitals& orb, const std::vector>> &abfs, const std::vector>> &orbs, diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 9b8382f871..6981b41915 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -26,7 +26,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->kmesh_times ); const std::vector>> - abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom( orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); + abfs = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell,orb, lcaos, this->kmesh_times, GlobalC::exx_info.info_ri.pca_threshold ); // ofs_mpi<<"memory:\t"< std::map>>>>> { Matrix_Orbs21 m_jyslcaos_lcaos; - m_jyslcaos_lcaos.init( 1, orb, this->kmesh_times, 1 ); + m_jyslcaos_lcaos.init( 1, ucell , orb, this->kmesh_times, 1 ); m_jyslcaos_lcaos.init_radial( jle.jle, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 - m_jyslcaos_lcaos.init_radial_table(radial_R); + m_jyslcaos_lcaos.init_radial_table(ucell.lat0, radial_R); #else m_jyslcaos_lcaos.init_radial_table(); #endif @@ -137,10 +137,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co 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, ucell , orb, this->kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 - m_abfslcaos_lcaos.init_radial_table(radial_R); + m_abfslcaos_lcaos.init_radial_table(ucell.lat0,radial_R); #else m_abfslcaos_lcaos.init_radial_table(); #endif From cb4cdd5d5516d4603e33843325c39bac64d68979 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 21:54:25 +0800 Subject: [PATCH 19/24] change ucell in module_ri/matirx_orb2.cpp --- source/module_ri/Matrix_Orbs22.cpp | 13 +++++++++---- source/module_ri/Matrix_Orbs22.h | 5 ++++- source/module_ri/exx_opt_orb.cpp | 4 ++-- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/source/module_ri/Matrix_Orbs22.cpp b/source/module_ri/Matrix_Orbs22.cpp index 027f252614..86bb5a865a 100644 --- a/source/module_ri/Matrix_Orbs22.cpp +++ b/source/module_ri/Matrix_Orbs22.cpp @@ -9,7 +9,11 @@ #include "module_base/tool_title.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Matrix_Orbs22::init(const int mode, const LCAO_Orbitals& orb, const double kmesh_times, const double rmesh_times) +void Matrix_Orbs22::init(const int mode, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const double kmesh_times, + const double rmesh_times) { ModuleBase::TITLE("Matrix_Orbs22", "init"); ModuleBase::timer::tick("Matrix_Orbs22", "init"); @@ -20,7 +24,7 @@ void Matrix_Orbs22::init(const int mode, const LCAO_Orbitals& orb, const double for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -129,7 +133,8 @@ void Matrix_Orbs22::init_radial_table() ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table"); } -void Matrix_Orbs22::init_radial_table(const std::map>>& Rs) +void Matrix_Orbs22::init_radial_table(const double lat0, + const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs22", "init_radial_table_Rs"); ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table"); @@ -152,7 +157,7 @@ void Matrix_Orbs22::init_radial_table(const std::map radials; for (const double& R: RsB.second) { - const double position = R * GlobalC::ucell.lat0 / lcao_dr_; + const double position = R * lat0 / lcao_dr_; const size_t iq = static_cast(position); for (size_t i = 0; i != 4; ++i) radials.insert(iq + i); diff --git a/source/module_ri/Matrix_Orbs22.h b/source/module_ri/Matrix_Orbs22.h index 0c7f08be3f..c30918a741 100644 --- a/source/module_ri/Matrix_Orbs22.h +++ b/source/module_ri/Matrix_Orbs22.h @@ -11,6 +11,7 @@ #include "module_basis/module_ao/ORB_gaunt_table.h" #include "module_basis/module_ao/ORB_read.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb22.h" +#include "module_cell/unitcell.h" #include #include @@ -23,6 +24,7 @@ class Matrix_Orbs22 // mode: // 1: void init(const int mode, + const UnitCell& ucell, const LCAO_Orbitals& orb, const double kmesh_times, // extend Kcut, keep dK const double rmesh_times); // extend Rcut, keep dR @@ -37,7 +39,8 @@ class Matrix_Orbs22 const LCAO_Orbitals& orb_B2); void init_radial_table(); - void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const double lat0, + const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 6981b41915..40bfe6c7f3 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -73,10 +73,10 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co 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, ucell,orb, this->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); + m_lcaoslcaos_lcaoslcaos.init_radial_table(ucell.lat0,radial_R); #else m_lcaoslcaos_lcaoslcaos.init_radial_table(); #endif From 372de4c069b4121c6cdecf7c0a0f0d8afcc6340e Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 22:16:51 +0800 Subject: [PATCH 20/24] change lat0 in Matrix.cpp --- source/module_ri/ABFs_Construct-PCA.cpp | 2 +- source/module_ri/Matrix_Orbs11.cpp | 5 ++--- source/module_ri/Matrix_Orbs11.h | 5 ++--- source/module_ri/Matrix_Orbs11.hpp | 4 ++-- source/module_ri/Matrix_Orbs21.cpp | 4 ++-- source/module_ri/Matrix_Orbs21.h | 5 ++--- source/module_ri/Matrix_Orbs21.hpp | 4 ++-- source/module_ri/Matrix_Orbs22.cpp | 3 +-- source/module_ri/Matrix_Orbs22.h | 5 ++--- source/module_ri/Matrix_Orbs22.hpp | 6 +++--- source/module_ri/exx_opt_orb.cpp | 12 ++++++------ 11 files changed, 25 insertions(+), 30 deletions(-) diff --git a/source/module_ri/ABFs_Construct-PCA.cpp b/source/module_ri/ABFs_Construct-PCA.cpp index d5fbfaf477..933b8c9804 100644 --- a/source/module_ri/ABFs_Construct-PCA.cpp +++ b/source/module_ri/ABFs_Construct-PCA.cpp @@ -104,7 +104,7 @@ namespace PCA for( std::size_t it=0; it!=abfs.size(); ++it ) { delta_R[it][it] = {0.0}; } - m_abfslcaos_lcaos.init_radial_table(ucell.lat0,delta_R); + m_abfslcaos_lcaos.init_radial_table(delta_R); GlobalC::exx_info.info_ri.abfs_Lmax = Lmax_bak; diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index ad01a78a6c..598be1f636 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -19,7 +19,7 @@ void Matrix_Orbs11::init(const int mode, ModuleBase::timer::tick("Matrix_Orbs11", "init"); int Lmax_used, Lmax; - + this->lat0 = ucell.lat0; const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) @@ -123,8 +123,7 @@ void Matrix_Orbs11::init_radial_table() ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); } -void Matrix_Orbs11::init_radial_table(const double lat0, - const std::map>>& Rs) +void Matrix_Orbs11::init_radial_table(const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs11", "init_radial_table_Rs"); ModuleBase::timer::tick("Matrix_Orbs11", "init_radial_table"); diff --git a/source/module_ri/Matrix_Orbs11.h b/source/module_ri/Matrix_Orbs11.h index a2afe85352..02dd663ee8 100644 --- a/source/module_ri/Matrix_Orbs11.h +++ b/source/module_ri/Matrix_Orbs11.h @@ -35,8 +35,7 @@ class Matrix_Orbs11 void init_radial(const LCAO_Orbitals& orb_A, const LCAO_Orbitals& orb_B); void init_radial_table(); - void init_radial_table(const double lat0, - const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { @@ -71,7 +70,7 @@ class Matrix_Orbs11 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - + double lat0 = 0.0; // restore ucell.lat0 std::map Matrix_Orbs11::cal_overlap_matrix( const size_t NB = co6.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const Tdata overlap = co6.second.cal_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA, MB ); + const Tdata overlap = co6.second.cal_overlap( tauA*lat0, tauB*lat0, MA, MB ); const size_t iA = index_A[TA][LA][NA][MA]; const size_t iB = index_B[TB][LB][NB][MB]; switch(matrix_order) @@ -103,7 +103,7 @@ std::array,3> Matrix_Orbs11::cal_grad_overlap_matrix( const size_t NB = co6.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const std::array grad_overlap = RI_Util::Vector3_to_array3(co6.second.cal_grad_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA, MB )); + const std::array grad_overlap = RI_Util::Vector3_to_array3(co6.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA, MB )); const size_t iA = index_A[TA][LA][NA][MA]; const size_t iB = index_B[TB][LB][NB][MB]; for(size_t i=0; ilat0 = ucell.lat0; for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); @@ -120,8 +121,7 @@ void Matrix_Orbs21::init_radial_table() ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } -void Matrix_Orbs21::init_radial_table(const double lat0, - const std::map>>& Rs) +void Matrix_Orbs21::init_radial_table(const std::map>>& Rs) { ModuleBase::TITLE("Matrix_Orbs21", "init_radial_table_Rs"); ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); diff --git a/source/module_ri/Matrix_Orbs21.h b/source/module_ri/Matrix_Orbs21.h index 57d85bf1c1..0e39be9f69 100644 --- a/source/module_ri/Matrix_Orbs21.h +++ b/source/module_ri/Matrix_Orbs21.h @@ -36,8 +36,7 @@ class Matrix_Orbs21 const LCAO_Orbitals& orb_B); void init_radial_table(); - void init_radial_table(const double lat0, - const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { @@ -79,7 +78,7 @@ class Matrix_Orbs21 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - + double lat0 = 0; // restore ucell.lat0 std::map Matrix_Orbs21::cal_overlap_matrix( const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const Tdata overlap = co8.second.cal_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB ); + const Tdata overlap = co8.second.cal_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB ); const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; @@ -139,7 +139,7 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const std::array grad_overlap = RI_Util::Vector3_to_array3(co8.second.cal_grad_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB )); + const std::array grad_overlap = RI_Util::Vector3_to_array3(co8.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB )); const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; diff --git a/source/module_ri/Matrix_Orbs22.cpp b/source/module_ri/Matrix_Orbs22.cpp index 86bb5a865a..2324f16943 100644 --- a/source/module_ri/Matrix_Orbs22.cpp +++ b/source/module_ri/Matrix_Orbs22.cpp @@ -133,8 +133,7 @@ void Matrix_Orbs22::init_radial_table() ModuleBase::timer::tick("Matrix_Orbs22", "init_radial_table"); } -void Matrix_Orbs22::init_radial_table(const double lat0, - const std::map>>& Rs) +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"); diff --git a/source/module_ri/Matrix_Orbs22.h b/source/module_ri/Matrix_Orbs22.h index c30918a741..fa6898a2e3 100644 --- a/source/module_ri/Matrix_Orbs22.h +++ b/source/module_ri/Matrix_Orbs22.h @@ -39,8 +39,7 @@ class Matrix_Orbs22 const LCAO_Orbitals& orb_B2); void init_radial_table(); - void init_radial_table(const double lat0, - const std::map>>& Rs); // unit: ucell.lat0 + void init_radial_table(const std::map>>& Rs); // unit: ucell.lat0 enum class Matrix_Order { @@ -103,7 +102,7 @@ class Matrix_Orbs22 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - + double lat0 = 0; // restore ucell.lat0 std::map< size_t, // TA std::map Matrix_Orbs22::cal_overlap_matrix( const size_t NB2 = co10.first; for( size_t MB2=0; MB2!=2*LB2+1; ++MB2 ) { - const Tdata overlap = co10.second.cal_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB1, MB2 ); + const Tdata overlap = co10.second.cal_overlap( tauA*lat0, tauB*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]; const size_t iB1 = index_B1[TB][LB1][NB1][MB1]; @@ -217,10 +217,10 @@ std::array,3> Matrix_Orbs22::cal_grad_overlap_matrix( const size_t NB2 = co10.first; for( size_t MB2=0; MB2!=2*LB2+1; ++MB2 ) { - const Tdata overlap = co10.second.cal_overlap( tauA*GlobalC::ucell.lat0, tauB*GlobalC::ucell.lat0, MA1, MA2, MB1, MB2 ); + const Tdata overlap = co10.second.cal_overlap( tauA*lat0, tauB*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 std::array grad_overlap = RI_Util::Vector3_to_array3(co10.second.cal_grad_overlap( tauA*lat0, tauB*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]; const size_t iB1 = index_B1[TB][LB1][NB1][MB1]; diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index 40bfe6c7f3..cb03080a52 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -76,7 +76,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_lcaoslcaos_lcaoslcaos.init( 1, ucell,orb, this->kmesh_times, 1 ); m_lcaoslcaos_lcaoslcaos.init_radial( lcaos, lcaos, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 - m_lcaoslcaos_lcaoslcaos.init_radial_table(ucell.lat0,radial_R); + m_lcaoslcaos_lcaoslcaos.init_radial_table(radial_R); #else m_lcaoslcaos_lcaoslcaos.init_radial_table(); #endif @@ -92,7 +92,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_jyslcaos_lcaos.init( 1, ucell , orb, this->kmesh_times, 1 ); m_jyslcaos_lcaos.init_radial( jle.jle, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 - m_jyslcaos_lcaos.init_radial_table(ucell.lat0, radial_R); + m_jyslcaos_lcaos.init_radial_table( radial_R); #else m_jyslcaos_lcaos.init_radial_table(); #endif @@ -108,7 +108,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_jys_jys.init( 2,ucell,orb, this->kmesh_times, 1 ); m_jys_jys.init_radial( jle.jle, jle.jle ); #if TEST_EXX_RADIAL>=1 - m_jys_jys.init_radial_table(ucell.lat0,radial_R); + m_jys_jys.init_radial_table(radial_R); #else m_jys_jys.init_radial_table(); #endif @@ -124,7 +124,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_abfs_abfs.init( 2, ucell, orb, this->kmesh_times, 1 ); m_abfs_abfs.init_radial( abfs, abfs ); #if TEST_EXX_RADIAL>=1 - m_abfs_abfs.init_radial_table(ucell.lat0,radial_R); + m_abfs_abfs.init_radial_table(radial_R); #else m_abfs_abfs.init_radial_table(); #endif @@ -140,7 +140,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_abfslcaos_lcaos.init( 1, ucell , orb, this->kmesh_times, 1 ); m_abfslcaos_lcaos.init_radial( abfs, lcaos, lcaos ); #if TEST_EXX_RADIAL>=1 - m_abfslcaos_lcaos.init_radial_table(ucell.lat0,radial_R); + m_abfslcaos_lcaos.init_radial_table(radial_R); #else m_abfslcaos_lcaos.init_radial_table(); #endif @@ -156,7 +156,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co m_jys_abfs.init( 2, ucell,orb, this->kmesh_times, 1 ); m_jys_abfs.init_radial( jle.jle, abfs ); #if TEST_EXX_RADIAL>=1 - m_jys_abfs.init_radial_table(ucell.lat0,radial_R); + m_jys_abfs.init_radial_table(radial_R); #else m_jys_abfs.init_radial_table(); #endif From ea19823edcacd00e1e0d94c3975542ee6e66c8df Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Sun, 8 Dec 2024 22:21:24 +0800 Subject: [PATCH 21/24] change ucell in Matrix.cpp --- source/module_ri/Matrix_Orbs11.h | 1 + source/module_ri/Matrix_Orbs11.hpp | 9 +++++---- source/module_ri/Matrix_Orbs21.h | 3 ++- source/module_ri/Matrix_Orbs21.hpp | 1 + source/module_ri/Matrix_Orbs22.h | 1 + source/module_ri/Matrix_Orbs22.hpp | 1 + source/module_ri/exx_opt_orb.cpp | 12 ++++++------ 7 files changed, 17 insertions(+), 11 deletions(-) diff --git a/source/module_ri/Matrix_Orbs11.h b/source/module_ri/Matrix_Orbs11.h index 02dd663ee8..f358060d62 100644 --- a/source/module_ri/Matrix_Orbs11.h +++ b/source/module_ri/Matrix_Orbs11.h @@ -63,6 +63,7 @@ class Matrix_Orbs11 template std::map>>>> cal_overlap_matrix_all( + const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM& index_r, const ModuleBase::Element_Basis_Index::IndexLNM& index_c) const; diff --git a/source/module_ri/Matrix_Orbs11.hpp b/source/module_ri/Matrix_Orbs11.hpp index 07d0e5e9cd..49f2089122 100644 --- a/source/module_ri/Matrix_Orbs11.hpp +++ b/source/module_ri/Matrix_Orbs11.hpp @@ -126,6 +126,7 @@ std::array,3> Matrix_Orbs11::cal_grad_overlap_matrix( template std::map>>>> Matrix_Orbs11::cal_overlap_matrix_all( + const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM &index_r, const ModuleBase::Element_Basis_Index::IndexLNM &index_c ) const { @@ -136,16 +137,16 @@ std::map &tauA( GlobalC::ucell.atoms[TA].tau[IA] ); + const ModuleBase::Vector3 &tauA( ucell.atoms[TA].tau[IA] ); for( const auto &co2 : co1.second ) { const size_t TB = co2.first; - for (size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB) + for (size_t IB=0; IB!=ucell.atoms[TB].na; ++IB) { - const ModuleBase::Vector3 &tauB( GlobalC::ucell.atoms[TB].tau[IB] ); + const ModuleBase::Vector3 &tauB( ucell.atoms[TB].tau[IB] ); matrixes[TA][IA][TB][IB] = cal_overlap_matrix( TA, TB, tauA, tauB, index_r, index_c, Matrix_Order::AB ); } diff --git a/source/module_ri/Matrix_Orbs21.h b/source/module_ri/Matrix_Orbs21.h index 0e39be9f69..5d7d858fa8 100644 --- a/source/module_ri/Matrix_Orbs21.h +++ b/source/module_ri/Matrix_Orbs21.h @@ -70,7 +70,8 @@ class Matrix_Orbs21 template std::map>>>>> - cal_overlap_matrix_all(const ModuleBase::Element_Basis_Index::IndexLNM& index_A1, + cal_overlap_matrix_all(const UnitCell& ucell, + 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; diff --git a/source/module_ri/Matrix_Orbs21.hpp b/source/module_ri/Matrix_Orbs21.hpp index 7779ce3f0b..529736575d 100644 --- a/source/module_ri/Matrix_Orbs21.hpp +++ b/source/module_ri/Matrix_Orbs21.hpp @@ -170,6 +170,7 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( template std::map>>>>> Matrix_Orbs21::cal_overlap_matrix_all( + const UnitCell &ucell, 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 diff --git a/source/module_ri/Matrix_Orbs22.h b/source/module_ri/Matrix_Orbs22.h index fa6898a2e3..016a3d303a 100644 --- a/source/module_ri/Matrix_Orbs22.h +++ b/source/module_ri/Matrix_Orbs22.h @@ -93,6 +93,7 @@ class Matrix_Orbs22 template std::map>>>> cal_overlap_matrix_all( + const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM& index_A1, const ModuleBase::Element_Basis_Index::IndexLNM& index_A2, const ModuleBase::Element_Basis_Index::IndexLNM& index_B1, diff --git a/source/module_ri/Matrix_Orbs22.hpp b/source/module_ri/Matrix_Orbs22.hpp index efe3facad5..e4038a0264 100644 --- a/source/module_ri/Matrix_Orbs22.hpp +++ b/source/module_ri/Matrix_Orbs22.hpp @@ -275,6 +275,7 @@ std::array,3> Matrix_Orbs22::cal_grad_overlap_matrix( // 4-parameter interface, for opt_orb template std::map < size_t, std::map>>>> Matrix_Orbs22::cal_overlap_matrix_all( + const UnitCell &ucell, const ModuleBase::Element_Basis_Index::IndexLNM &index_A1, const ModuleBase::Element_Basis_Index::IndexLNM &index_A2, const ModuleBase::Element_Basis_Index::IndexLNM &index_B1, diff --git a/source/module_ri/exx_opt_orb.cpp b/source/module_ri/exx_opt_orb.cpp index cb03080a52..37cd742d7a 100644 --- a/source/module_ri/exx_opt_orb.cpp +++ b/source/module_ri/exx_opt_orb.cpp @@ -80,7 +80,7 @@ void Exx_Opt_Orb::generate_matrix(const K_Vectors &kv, const UnitCell& ucell, co #else m_lcaoslcaos_lcaoslcaos.init_radial_table(); #endif - return m_lcaoslcaos_lcaoslcaos.cal_overlap_matrix_all( index_lcaos, index_lcaos, index_lcaos, index_lcaos); + return m_lcaoslcaos_lcaoslcaos.cal_overlap_matrix_all(ucell,index_lcaos, index_lcaos, index_lcaos, index_lcaos); }(); // ofs_mpi<<"memory:\t"<( index_jys, index_lcaos, index_lcaos ); + return m_jyslcaos_lcaos.cal_overlap_matrix_all(ucell,index_jys, index_lcaos, index_lcaos ); }(); // ofs_mpi<<"memory:\t"<( index_jys, index_jys ); + return m_jys_jys.cal_overlap_matrix_all(ucell,index_jys, index_jys ); }(); // ofs_mpi<<"memory:\t"<( index_abfs, index_abfs ); + return m_abfs_abfs.cal_overlap_matrix_all(ucell,index_abfs, index_abfs ); }(); // ofs_mpi<<"memory:\t"<( index_abfs, index_lcaos, index_lcaos ); + return m_abfslcaos_lcaos.cal_overlap_matrix_all(ucell,index_abfs, index_lcaos, index_lcaos ); }(); // ofs_mpi<<"memory:\t"<( index_jys, index_abfs ); + return m_jys_abfs.cal_overlap_matrix_all(ucell,index_jys, index_abfs ); }(); // ofs_mpi<<"memory:\t"< Date: Sun, 8 Dec 2024 22:25:22 +0800 Subject: [PATCH 22/24] change ucell in module_ri/matirx_orb22,matirx_orb11.cpp --- source/module_ri/Matrix_Orbs21.hpp | 8 ++++---- source/module_ri/Matrix_Orbs22.hpp | 12 ++++++------ .../module_exx_symmetry/symmetry_rotation.cpp | 12 ++++++------ 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/source/module_ri/Matrix_Orbs21.hpp b/source/module_ri/Matrix_Orbs21.hpp index 529736575d..0d8a99049c 100644 --- a/source/module_ri/Matrix_Orbs21.hpp +++ b/source/module_ri/Matrix_Orbs21.hpp @@ -182,16 +182,16 @@ std::map &tauA( GlobalC::ucell.atoms[TA].tau[IA] ); + const ModuleBase::Vector3 &tauA( ucell.atoms[TA].tau[IA] ); for( const auto &co2 : co1.second ) { const size_t TB = co2.first; - for( size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) + for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB ) { - const ModuleBase::Vector3 &tauB( GlobalC::ucell.atoms[TB].tau[IB] ); + const ModuleBase::Vector3 &tauB( ucell.atoms[TB].tau[IB] ); const RI::Tensor &&m = cal_overlap_matrix( TA, TB, tauA, tauB, index_A1, index_A2, index_B, Matrix_Order::A2BA1 ); matrixes[TA][IA][TB][IB].resize(2); diff --git a/source/module_ri/Matrix_Orbs22.hpp b/source/module_ri/Matrix_Orbs22.hpp index e4038a0264..436a952c91 100644 --- a/source/module_ri/Matrix_Orbs22.hpp +++ b/source/module_ri/Matrix_Orbs22.hpp @@ -286,22 +286,22 @@ std::map < size_t, std::map &tauA( GlobalC::ucell.atoms[TA].tau[IA] ); + const ModuleBase::Vector3 &tauA( ucell.atoms[TA].tau[IA] ); for( const auto &co2 : co1.second ) { const size_t TB = co2.first; - for( size_t IB=0; IB!=GlobalC::ucell.atoms[TB].na; ++IB ) + for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB ) { - const ModuleBase::Vector3 &tauB( GlobalC::ucell.atoms[TB].tau[IB] ); + const ModuleBase::Vector3 &tauB( ucell.atoms[TB].tau[IB] ); matrixes[TA][IA][TB][IB] = cal_overlap_matrix( TA, TB, - GlobalC::ucell.atoms[TA].tau[IA], - GlobalC::ucell.atoms[TB].tau[IB], + ucell.atoms[TA].tau[IA], + ucell.atoms[TB].tau[IB], index_A1, index_A2, index_B1, diff --git a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp index bfdef1fab9..9270ef6bf9 100644 --- a/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp +++ b/source/module_ri/module_exx_symmetry/symmetry_rotation.cpp @@ -321,13 +321,13 @@ namespace ModuleSymmetry ofs << gmatc[isym].e21 << " " << gmatc[isym].e22 << " " << gmatc[isym].e23 << std::endl; ofs << gmatc[isym].e31 << " " << gmatc[isym].e32 << " " << gmatc[isym].e33 << std::endl; ofs << "gmatrix_direct=" << std::endl; - ofs << GlobalC::ucell.symm.gmatrix[isym].e11 << " " << GlobalC::ucell.symm.gmatrix[isym].e12 << " " << GlobalC::ucell.symm.gmatrix[isym].e13 << std::endl; - ofs << GlobalC::ucell.symm.gmatrix[isym].e21 << " " << GlobalC::ucell.symm.gmatrix[isym].e22 << " " << GlobalC::ucell.symm.gmatrix[isym].e23 << std::endl; - ofs << GlobalC::ucell.symm.gmatrix[isym].e31 << " " << GlobalC::ucell.symm.gmatrix[isym].e32 << " " << GlobalC::ucell.symm.gmatrix[isym].e33 << std::endl; + ofs << ucell.symm.gmatrix[isym].e11 << " " << ucell.symm.gmatrix[isym].e12 << " " << ucell.symm.gmatrix[isym].e13 << std::endl; + ofs << ucell.symm.gmatrix[isym].e21 << " " << ucell.symm.gmatrix[isym].e22 << " " << ucell.symm.gmatrix[isym].e23 << std::endl; + ofs << ucell.symm.gmatrix[isym].e31 << " " << ucell.symm.gmatrix[isym].e32 << " " << ucell.symm.gmatrix[isym].e33 << std::endl; ofs << "kgmatrix_direct=" << std::endl; - ofs << GlobalC::ucell.symm.kgmatrix[isym].e11 << " " << GlobalC::ucell.symm.kgmatrix[isym].e12 << " " << GlobalC::ucell.symm.kgmatrix[isym].e13 << std::endl; - ofs << GlobalC::ucell.symm.kgmatrix[isym].e21 << " " << GlobalC::ucell.symm.kgmatrix[isym].e22 << " " << GlobalC::ucell.symm.kgmatrix[isym].e23 << std::endl; - ofs << GlobalC::ucell.symm.kgmatrix[isym].e31 << " " << GlobalC::ucell.symm.kgmatrix[isym].e32 << " " << GlobalC::ucell.symm.kgmatrix[isym].e33 << std::endl; + ofs << ucell.symm.kgmatrix[isym].e11 << " " << ucell.symm.kgmatrix[isym].e12 << " " << ucell.symm.kgmatrix[isym].e13 << std::endl; + ofs << ucell.symm.kgmatrix[isym].e21 << " " << ucell.symm.kgmatrix[isym].e22 << " " << ucell.symm.kgmatrix[isym].e23 << std::endl; + ofs << ucell.symm.kgmatrix[isym].e31 << " " << ucell.symm.kgmatrix[isym].e32 << " " << ucell.symm.kgmatrix[isym].e33 << std::endl; ofs << "euler_angle/pi: " << euler_angles_test[isym].x / ModuleBase::PI << " " << euler_angles_test[isym].y / ModuleBase::PI << " " << euler_angles_test[isym].z / ModuleBase::PI << std::endl; for (int l = 0;l <= lmax;++l) From 704f0b892cb3bf047f45ef9a7fa14972002acb64 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Mon, 9 Dec 2024 01:56:31 +0000 Subject: [PATCH 23/24] [pre-commit.ci lite] apply automatic fixes --- source/module_ri/Matrix_Orbs21.cpp | 108 ++++++++++++++++--------- source/module_ri/exx_opt_orb-print.cpp | 13 +-- 2 files changed, 78 insertions(+), 43 deletions(-) diff --git a/source/module_ri/Matrix_Orbs21.cpp b/source/module_ri/Matrix_Orbs21.cpp index 89267d7b04..da363fdc89 100644 --- a/source/module_ri/Matrix_Orbs21.cpp +++ b/source/module_ri/Matrix_Orbs21.cpp @@ -62,21 +62,29 @@ void Matrix_Orbs21::init_radial(const std::vectorMGT))); + this->MGT))); +} +} +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); } @@ -87,21 +95,29 @@ void Matrix_Orbs21::init_radial(const std::vectorMGT))); + this->MGT))); +} +} +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); } @@ -109,15 +125,23 @@ void Matrix_Orbs21::init_radial_table() { ModuleBase::TITLE("Matrix_Orbs21", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); - for (auto& coA: center2_orb21_s) - for (auto& coB: coA.second) - for (auto& coC: coB.second) - for (auto& coD: coC.second) - for (auto& coE: coD.second) - for (auto& coF: coE.second) - for (auto& coG: coF.second) - for (auto& coH: coG.second) - coH.second.init_radial_table(); + for (auto& coA: center2_orb21_s) { + for (auto& coB: coA.second) { + for (auto& coC: coB.second) { + for (auto& coD: coC.second) { + for (auto& coE: coD.second) { + for (auto& coF: coE.second) { + for (auto& coG: coF.second) { + for (auto& coH: coG.second) { + coH.second.init_radial_table(); +} +} +} +} +} +} +} +} ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } @@ -125,7 +149,7 @@ void Matrix_Orbs21::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_orb21_sAB) - for (auto& coD: coC.second) - for (auto& coE: coD.second) - for (auto& coF: coE.second) - for (auto& coG: coF.second) - for (auto& coH: coG.second) - coH.second.init_radial_table(radials); + for (auto& coC: *center2_orb21_sAB) { + for (auto& coD: coC.second) { + for (auto& coE: coD.second) { + for (auto& coF: coE.second) { + for (auto& coG: coF.second) { + for (auto& coH: coG.second) { + coH.second.init_radial_table(radials); +} +} +} +} +} +} } - } + } +} ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } diff --git a/source/module_ri/exx_opt_orb-print.cpp b/source/module_ri/exx_opt_orb-print.cpp index a768d8399b..9c78206b6f 100644 --- a/source/module_ri/exx_opt_orb-print.cpp +++ b/source/module_ri/exx_opt_orb-print.cpp @@ -65,10 +65,11 @@ void Exx_Opt_Orb::print_matrix( // this parameter determine the total number of jlq. ofs << Exx_Abfs::Jle::Ecut_exx << " ecutwfc_jlq" << std::endl;//mohan modify 2009-09-08 - if(TA==TB) + if(TA==TB) { ofs << orb_cutoff[TA] << " rcut_Jlq" << std::endl; - else + } else { 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; @@ -85,8 +86,9 @@ void Exx_Opt_Orb::print_matrix( 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 ) + for( size_t L = 0; L!=range_jles[T].size(); ++L ) { 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); @@ -94,10 +96,11 @@ void Exx_Opt_Orb::print_matrix( 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 - if(TA==TB) + if(TA==TB) { ofs << ecut_numberA << " ne" << std::endl; - else + } else { ofs << ecut_numberA << " " << ecut_numberB << " ne" << std::endl; +} ofs << "" << std::endl; for( int ik=0; ik!=kv.get_nkstot(); ++ik ) From 80eb6a16934d1a70a9d9644434ee801f63de618c Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Tue, 10 Dec 2024 10:50:28 +0800 Subject: [PATCH 24/24] update reference of lat0 in ri --- source/module_ri/Matrix_Orbs11.cpp | 3 +- source/module_ri/Matrix_Orbs11.h | 2 +- source/module_ri/Matrix_Orbs11.hpp | 2 + source/module_ri/Matrix_Orbs21.cpp | 172 +++++++++++++++++------------ source/module_ri/Matrix_Orbs21.h | 2 +- source/module_ri/Matrix_Orbs21.hpp | 2 + source/module_ri/Matrix_Orbs22.cpp | 2 + source/module_ri/Matrix_Orbs22.h | 2 +- source/module_ri/Matrix_Orbs22.hpp | 3 +- 9 files changed, 115 insertions(+), 75 deletions(-) diff --git a/source/module_ri/Matrix_Orbs11.cpp b/source/module_ri/Matrix_Orbs11.cpp index 598be1f636..01572cbc8d 100644 --- a/source/module_ri/Matrix_Orbs11.cpp +++ b/source/module_ri/Matrix_Orbs11.cpp @@ -19,7 +19,7 @@ void Matrix_Orbs11::init(const int mode, ModuleBase::timer::tick("Matrix_Orbs11", "init"); int Lmax_used, Lmax; - this->lat0 = ucell.lat0; + this->lat0 = &ucell.lat0; const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) @@ -127,6 +127,7 @@ void Matrix_Orbs11::init_radial_table(const std::maplat0; for (const auto& RsA: Rs) { for (const auto& RsB: RsA.second) { diff --git a/source/module_ri/Matrix_Orbs11.h b/source/module_ri/Matrix_Orbs11.h index f358060d62..b994f93265 100644 --- a/source/module_ri/Matrix_Orbs11.h +++ b/source/module_ri/Matrix_Orbs11.h @@ -71,7 +71,7 @@ class Matrix_Orbs11 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - double lat0 = 0.0; // restore ucell.lat0 + double* lat0=nullptr; // restore ucell.lat0 std::map Matrix_Orbs11::cal_overlap_matrix( const Matrix_Order &matrix_order) const { RI::Tensor m; + const double lat0 = *this->lat0; const size_t sizeA = index_A[TA].count_size; const size_t sizeB = index_B[TB].count_size; switch(matrix_order) @@ -75,6 +76,7 @@ std::array,3> Matrix_Orbs11::cal_grad_overlap_matrix( const Matrix_Order &matrix_order) const { std::array,3> m; + const double lat0 = *this->lat0; const size_t sizeA = index_A[TA].count_size; const size_t sizeB = index_B[TB].count_size; for(int i=0; ilat0 = ucell.lat0; + this->lat0 = &ucell.lat0; for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); @@ -62,29 +62,37 @@ void Matrix_Orbs21::init_radial(const std::vectorMGT))); -} -} -} -} -} -} -} -} + this->MGT))); + } + } + } + } + } + } + } + } ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); } @@ -95,29 +103,37 @@ void Matrix_Orbs21::init_radial(const std::vectorMGT))); -} -} -} -} -} -} -} -} + this->MGT))); + } + } + } + } + } + } + } + } ModuleBase::timer::tick("Matrix_Orbs21", "init_radial"); } @@ -125,23 +141,31 @@ void Matrix_Orbs21::init_radial_table() { ModuleBase::TITLE("Matrix_Orbs21", "init_radial"); ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); - for (auto& coA: center2_orb21_s) { - for (auto& coB: coA.second) { - for (auto& coC: coB.second) { - for (auto& coD: coC.second) { - for (auto& coE: coD.second) { - for (auto& coF: coE.second) { - for (auto& coG: coF.second) { - for (auto& coH: coG.second) { - coH.second.init_radial_table(); -} -} -} -} -} -} -} -} + for (auto& coA: center2_orb21_s) + { + for (auto& coB: coA.second) + { + for (auto& coC: coB.second) + { + for (auto& coD: coC.second) + { + for (auto& coE: coD.second) + { + for (auto& coF: coE.second) + { + for (auto& coG: coF.second) + { + for (auto& coH: coG.second) + { + coH.second.init_radial_table(); + } + } + } + } + } + } + } + } ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } @@ -149,6 +173,7 @@ void Matrix_Orbs21::init_radial_table(const std::maplat0; for (const auto& RsA: Rs) { for (const auto& RsB: RsA.second) { @@ -164,25 +189,32 @@ void Matrix_Orbs21::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_orb21_sAB) + { + for (auto& coD: coC.second) + { + for (auto& coE: coD.second) + { + for (auto& coF: coE.second) + { + for (auto& coG: coF.second) + { + for (auto& coH: coG.second) + { + coH.second.init_radial_table(radials); + } + } + } + } + } } - for (auto& coC: *center2_orb21_sAB) { - for (auto& coD: coC.second) { - for (auto& coE: coD.second) { - for (auto& coF: coE.second) { - for (auto& coG: coF.second) { - for (auto& coH: coG.second) { - coH.second.init_radial_table(radials); -} -} -} -} -} -} } - } -} + } + } ModuleBase::timer::tick("Matrix_Orbs21", "init_radial_table"); } diff --git a/source/module_ri/Matrix_Orbs21.h b/source/module_ri/Matrix_Orbs21.h index 5d7d858fa8..1d2895b4b2 100644 --- a/source/module_ri/Matrix_Orbs21.h +++ b/source/module_ri/Matrix_Orbs21.h @@ -79,7 +79,7 @@ class Matrix_Orbs21 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - double lat0 = 0; // restore ucell.lat0 + double* lat0 = nullptr; // restore ucell.lat0 std::map Matrix_Orbs21::cal_overlap_matrix( const Matrix_Order &matrix_order) const { RI::Tensor m; + const double lat0 = *this->lat0; const size_t sizeA1 = index_A1[TA].count_size; const size_t sizeA2 = index_A2[TA].count_size; const size_t sizeB = index_B[TB].count_size; @@ -98,6 +99,7 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( const Matrix_Order &matrix_order) const { std::array,3> m; + const double lat0 = *this->lat0; const size_t sizeA1 = index_A1[TA].count_size; const size_t sizeA2 = index_A2[TA].count_size; const size_t sizeB = index_B[TB].count_size; diff --git a/source/module_ri/Matrix_Orbs22.cpp b/source/module_ri/Matrix_Orbs22.cpp index 2324f16943..d2592d5ce5 100644 --- a/source/module_ri/Matrix_Orbs22.cpp +++ b/source/module_ri/Matrix_Orbs22.cpp @@ -19,6 +19,7 @@ void Matrix_Orbs22::init(const int mode, ModuleBase::timer::tick("Matrix_Orbs22", "init"); int Lmax_used, Lmax; + this->lat0 = &ucell.lat0; const int ntype = orb.get_ntype(); int lmax_orb = -1, lmax_beta = -1; for (int it = 0; it < ntype; it++) @@ -137,6 +138,7 @@ void Matrix_Orbs22::init_radial_table(const std::maplat0; for (const auto& RsA: Rs) for (const auto& RsB: RsA.second) { diff --git a/source/module_ri/Matrix_Orbs22.h b/source/module_ri/Matrix_Orbs22.h index 016a3d303a..c3fa5331f8 100644 --- a/source/module_ri/Matrix_Orbs22.h +++ b/source/module_ri/Matrix_Orbs22.h @@ -103,7 +103,7 @@ class Matrix_Orbs22 ModuleBase::Sph_Bessel_Recursive::D2* psb_ = nullptr; ORB_gaunt_table MGT; const double lcao_dr_ = 0.01; - double lat0 = 0; // restore ucell.lat0 + double* lat0 = nullptr; // restore ucell.lat0 std::map< size_t, // TA std::map Matrix_Orbs22::cal_overlap_matrix( const ModuleBase::Element_Basis_Index::IndexLNM &index_B2, const Matrix_Order &matrix_order) const { + const double lat0 = *this->lat0; RI::Tensor m; const size_t sizeA1 = index_A1[TA].count_size; const size_t sizeA2 = index_A2[TA].count_size; @@ -184,7 +185,7 @@ std::array,3> Matrix_Orbs22::cal_grad_overlap_matrix( default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__)); } } - + const double lat0 = *this->lat0; for( const auto &co3 : center2_orb22_s.at(TA).at(TB) ) { const int LA1 = co3.first;