From 66f7bab643468e922ea0fc97b443694724ad5b72 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Sat, 13 Jul 2024 02:14:24 +0800 Subject: [PATCH] remove LM --- source/module_esolver/esolver_ks_lcao.cpp | 100 ++++++++---------- source/module_esolver/esolver_ks_lcao.h | 3 - .../module_esolver/esolver_ks_lcao_elec.cpp | 17 +-- .../module_esolver/esolver_ks_lcao_tddft.cpp | 52 +++++---- source/module_esolver/io_npz.cpp | 13 ++- .../hamilt_lcaodft/LCAO_matrix.h | 35 ------ .../hamilt_lcaodft/hamilt_lcao.cpp | 33 +++--- .../hamilt_lcaodft/hamilt_lcao.h | 27 +++-- .../operator_lcao/op_exx_lcao.cpp | 7 +- .../operator_lcao/op_exx_lcao.h | 10 +- .../operator_lcao/op_exx_lcao.hpp | 45 ++++---- .../hamilt_lcaodft/spar_hsr.cpp | 38 ++++--- .../hamilt_lcaodft/spar_hsr.h | 23 ++-- source/module_io/output_mat_sparse.cpp | 33 +++--- source/module_io/output_mat_sparse.h | 26 ++--- source/module_io/write_HS_R.cpp | 70 ++++++------ source/module_io/write_HS_R.h | 29 ++--- source/module_io/write_Vxc.hpp | 40 +++---- source/module_lr/esolver_lrtd_lcao.cpp | 4 +- 19 files changed, 283 insertions(+), 322 deletions(-) delete mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 06c92693a2..d1bcf6ecde 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -63,22 +63,6 @@ ESolver_KS_LCAO::ESolver_KS_LCAO() { this->classname = "ESolver_KS_LCAO"; this->basisname = "LCAO"; - -// the following EXX code should be removed to other places, mohan 2024/05/11 -#ifdef __EXX - if (GlobalC::exx_info.info_ri.real_number) - { - this->exx_lri_double = std::make_shared>(GlobalC::exx_info.info_ri); - this->exd = std::make_shared>(this->exx_lri_double); - this->LM.Hexxd = &this->exd->get_Hexxs(); - } - else - { - this->exx_lri_complex = std::make_shared>>(GlobalC::exx_info.info_ri); - this->exc = std::make_shared>>(this->exx_lri_complex); - this->LM.Hexxc = &this->exc->get_Hexxs(); - } -#endif } //------------------------------------------------------------------------------ @@ -162,9 +146,6 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) this->init_basis_lcao(inp, ucell); //------------------init Basis_lcao---------------------- - //! pass basis-pointer to EState and Psi - this->LM.ParaV = &(this->ParaV); - // 5) initialize density matrix dynamic_cast*>(this->pelec) ->init_DM(&this->kv, &(this->ParaV), GlobalV::NSPIN); @@ -186,14 +167,20 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) // PLEASE simplify the Exx_Global interface if (GlobalV::CALCULATION == "scf" || GlobalV::CALCULATION == "relax" || GlobalV::CALCULATION == "cell-relax" - || GlobalV::CALCULATION == "md") { - if (GlobalC::exx_info.info_global.cal_exx) { + || GlobalV::CALCULATION == "md") + { + if (GlobalC::exx_info.info_global.cal_exx) + { XC_Functional::set_xc_first_loop(ucell); - if (GlobalC::exx_info.info_ri.real_number) { + // initialize 2-center radial tables for EXX-LRI + if (GlobalC::exx_info.info_ri.real_number) + { + this->exx_lri_double = std::make_shared>(GlobalC::exx_info.info_ri); this->exx_lri_double->init(MPI_COMM_WORLD, this->kv); } else { + this->exx_lri_complex = std::make_shared>>(GlobalC::exx_info.info_ri); this->exx_lri_complex->init(MPI_COMM_WORLD, this->kv); } } @@ -202,7 +189,7 @@ void ESolver_KS_LCAO::before_all_runners(Input& inp, UnitCell& ucell) // 8) initialize DFT+U if (GlobalV::dft_plus_u) { - GlobalC::dftu.init(ucell, this->LM.ParaV, this->kv.get_nks()); + GlobalC::dftu.init(ucell, &this->ParaV, this->kv.get_nks()); } // 9) initialize ppcell @@ -276,7 +263,7 @@ void ESolver_KS_LCAO::init_after_vc(Input& inp, UnitCell& ucell) this->pw_big); dynamic_cast*>(this->pelec) - ->init_DM(&this->kv, this->LM.ParaV, GlobalV::NSPIN); + ->init_DM(&this->kv, &this->ParaV, GlobalV::NSPIN); GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho); @@ -473,21 +460,26 @@ void ESolver_KS_LCAO::after_all_runners() if (PARAM.inp.out_mat_xc) { ModuleIO::write_Vxc(GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::DRANK, - *this->psi, - GlobalC::ucell, - this->sf, - *this->pw_rho, - *this->pw_rhod, - GlobalC::ppcell.vloc, - *this->pelec->charge, - this->GG, - this->GK, - this->LM, - this->kv, - this->pelec->wg, - GlobalC::GridD); + GlobalV::NLOCAL, + GlobalV::DRANK, + &this->ParaV, + *this->psi, + GlobalC::ucell, + this->sf, + *this->pw_rho, + *this->pw_rhod, + GlobalC::ppcell.vloc, + *this->pelec->charge, + this->GG, + this->GK, + this->kv, + this->pelec->wg, + GlobalC::GridD +#ifdef __EXX + , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr + , this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr +#endif + ); } ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners"); @@ -848,7 +840,7 @@ void ESolver_KS_LCAO::hamilt2density(int istep, int iter, double ethr) if (GlobalV::sc_mag_switch) { SpinConstrain& sc = SpinConstrain::getScInstance(); - sc.cal_MW(iter, &(this->LM)); + sc.cal_MW(iter, this->p_hamilt); } // 9) use new charge density to calculate energy @@ -1015,8 +1007,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) for (int ik = 0; ik < this->kv.get_nks(); ++ik) { Hexxk_save.set_zero_hk(); - hamilt::OperatorEXX> opexx_save(&Hexxk_save, - &this->LM, + hamilt::OperatorEXX> opexx_save(&Hexxk_save, nullptr, this->kv); @@ -1069,7 +1060,7 @@ void ESolver_KS_LCAO::iter_finish(int iter) if (GlobalV::sc_mag_switch) { SpinConstrain& sc = SpinConstrain::getScInstance(); - sc.cal_MW(iter, &(this->LM)); + sc.cal_MW(iter, this->p_hamilt); } // 6) calculate the total energy. @@ -1353,18 +1344,17 @@ template ModuleIO::Output_Mat_Sparse ESolver_KS_LCAO::create_Output_Mat_Sparse(int istep) { return ModuleIO::Output_Mat_Sparse(hsolver::HSolverLCAO::out_mat_hsR, - hsolver::HSolverLCAO::out_mat_dh, - hsolver::HSolverLCAO::out_mat_t, - INPUT.out_mat_r, - istep, - this->pelec->pot->get_effective_v(), - this->ParaV, - this->GK, // mohan add 2024-04-01 - two_center_bundle_, - this->LM, - GlobalC::GridD, // mohan add 2024-04-06 - this->kv, - this->p_hamilt); + hsolver::HSolverLCAO::out_mat_dh, + hsolver::HSolverLCAO::out_mat_t, + INPUT.out_mat_r, + istep, + this->pelec->pot->get_effective_v(), + this->ParaV, + this->GK, // mohan add 2024-04-01 + two_center_bundle_, + GlobalC::GridD, // mohan add 2024-04-06 + this->kv, + this->p_hamilt); } //------------------------------------------------------------------------------ diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 6a1860f6e6..b82af80a9c 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -76,9 +76,6 @@ class ESolver_KS_LCAO : public ESolver_KS { // used for gamma only algorithms. Gint_Gamma GG; - // we will get rid of this class soon, don't use it, mohan 2024-03-28 - LCAO_Matrix LM; - Grid_Technique GridT; TwoCenterBundle two_center_bundle_; diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index aa5baf3a02..0bed54ff71 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -100,7 +100,7 @@ void ESolver_KS_LCAO::set_matrix_grid(Record_adj& ra) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - Parallel_Orbitals* pv = this->LM.ParaV; + Parallel_Orbitals* pv = &this->ParaV; ra.for_2d(*pv, GlobalV::GAMMA_ONLY_LOCAL); if (!GlobalV::GAMMA_ONLY_LOCAL) { @@ -183,17 +183,17 @@ void ESolver_KS_LCAO::beforesolver(const int istep) this->p_hamilt = new hamilt::HamiltLCAO( GlobalV::GAMMA_ONLY_LOCAL ? &(this->GG) : nullptr, GlobalV::GAMMA_ONLY_LOCAL ? nullptr : &(this->GK), - &(this->LM), &this->ParaV, this->pelec->pot, this->kv, two_center_bundle_, + DM #ifdef __EXX - DM, - GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step); -#else - DM); + , GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step + , GlobalC::exx_info.info_ri.real_number ? &exx_lri_double->Hexxs : nullptr + , GlobalC::exx_info.info_ri.real_number ? nullptr : &exx_lri_complex->Hexxs #endif + ); } #ifdef __DEEPKS @@ -201,7 +201,7 @@ void ESolver_KS_LCAO::beforesolver(const int istep) // since it depends on ionic positions if (GlobalV::deepks_setorb) { - const Parallel_Orbitals* pv = this->LM.ParaV; + const Parallel_Orbitals* pv = &this->ParaV; // build and save at beginning GlobalC::ld.build_psialpha(GlobalV::CAL_FORCE, GlobalC::ucell, @@ -284,10 +284,12 @@ void ESolver_KS_LCAO::before_scf(const int istep) #ifdef __EXX // set xc type before the first cal of xc in pelec->init_scf if (GlobalC::exx_info.info_ri.real_number) { + this->exd = std::make_shared>(exx_lri_double); this->exd->exx_beforescf(this->kv, *this->p_chgmix); } else { + this->exc = std::make_shared>>(exx_lri_complex); this->exc->exx_beforescf(this->kv, *this->p_chgmix); } #endif // __EXX @@ -583,7 +585,6 @@ void ESolver_KS_LCAO, std::complex>::get_S(void) GlobalV::test_atom_input); this->RA.for_2d(this->ParaV, GlobalV::GAMMA_ONLY_LOCAL); - this->LM.ParaV = &this->ParaV; if (this->p_hamilt == nullptr) { this->p_hamilt = new hamilt::HamiltLCAO, std::complex>( diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 3263c1dd33..45d176ee4a 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -89,18 +89,14 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(Input& inp, UnitCell& ucell) // 5) allocate H and S matrices according to computational resources LCAO_domain::divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, ParaV, kv.get_nks()); - // this part will be updated soon - // pass Hamilt-pointer to Operator - this->LM.ParaV = &(this->ParaV); - // 6) initialize Density Matrix dynamic_cast>*>(this->pelec) - ->init_DM(&kv, this->LM.ParaV, GlobalV::NSPIN); + ->init_DM(&kv, &this->ParaV, GlobalV::NSPIN); // 7) initialize Hsolver if (this->phsol == nullptr) { - this->phsol = new hsolver::HSolverLCAO>(this->LM.ParaV); + this->phsol = new hsolver::HSolverLCAO>(&this->ParaV); this->phsol->method = GlobalV::KS_SOLVER; } @@ -256,28 +252,28 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) if (hsolver::HSolverLCAO>::out_mat_hs[0]) { ModuleIO::save_mat(istep, - h_mat.p, - GlobalV::NLOCAL, - bit, - hsolver::HSolverLCAO>::out_mat_hs[1], - 1, - GlobalV::out_app_flag, - "H", - "data-" + std::to_string(ik), - *this->LM.ParaV, - GlobalV::DRANK); + h_mat.p, + GlobalV::NLOCAL, + bit, + hsolver::HSolverLCAO>::out_mat_hs[1], + 1, + GlobalV::out_app_flag, + "H", + "data-" + std::to_string(ik), + this->ParaV, + GlobalV::DRANK); ModuleIO::save_mat(istep, - s_mat.p, - GlobalV::NLOCAL, - bit, - hsolver::HSolverLCAO>::out_mat_hs[1], - 1, - GlobalV::out_app_flag, - "S", - "data-" + std::to_string(ik), - *this->LM.ParaV, - GlobalV::DRANK); + s_mat.p, + GlobalV::NLOCAL, + bit, + hsolver::HSolverLCAO>::out_mat_hs[1], + 1, + GlobalV::out_app_flag, + "S", + "data-" + std::to_string(ik), + this->ParaV, + GlobalV::DRANK); } } } @@ -311,8 +307,8 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter) } const int nloc = this->ParaV.nloc; - const int ncol_nbands = this->LM.ParaV->ncol_bands; - const int nrow = this->LM.ParaV->nrow; + const int ncol_nbands = this->ParaV.ncol_bands; + const int nrow = this->ParaV.nrow; const int nbands = GlobalV::NBANDS; const int nlocal = GlobalV::NLOCAL; diff --git a/source/module_esolver/io_npz.cpp b/source/module_esolver/io_npz.cpp index 3418bd2527..6edb2bf4a5 100644 --- a/source/module_esolver/io_npz.cpp +++ b/source/module_esolver/io_npz.cpp @@ -458,17 +458,16 @@ void ESolver_KS_LCAO::output_mat_npz(std::string& zipname, const hamilt: } #else - const Parallel_Orbitals* paraV = this->LM->ParaV; - auto row_indexes = paraV->get_indexes_row(); - auto col_indexes = paraV->get_indexes_col(); + auto row_indexes = paraV.get_indexes_row(); + auto col_indexes = paraV.get_indexes_col(); for(int iap=0;iapatom_begin_row[atom_i]; - int start_j = paraV->atom_begin_col[atom_j]; - int row_size = paraV->get_row_size(atom_i); - int col_size = paraV->get_col_size(atom_j); + int start_i = paraV.atom_begin_row[atom_i]; + int start_j = paraV.atom_begin_col[atom_j]; + int row_size = paraV.get_row_size(atom_i); + int col_size = paraV.get_col_size(atom_j); for(int iR=0;iR in 2021-12-2, will be deleted in the future -#include "module_base/abfs-vector3_order.h" -#ifdef __EXX -#include -#endif - -class LCAO_Matrix { - public: - LCAO_Matrix(){}; - ~LCAO_Matrix(){}; - - Parallel_Orbitals* ParaV; - -#ifdef __EXX - using TAC = std::pair>; - std::vector>>>* Hexxd; - std::vector>>>>* - Hexxc; - /// @brief Hexxk for all k-points, only for the 1st scf loop ofrestart load - std::vector> Hexxd_k_load; - std::vector>> Hexxc_k_load; -#endif -}; - -#endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 9d50cf1a85..a59c383e2d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -58,15 +58,18 @@ HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& template HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, - Gint_k* GK_in, - LCAO_Matrix* LM_in, - const Parallel_Orbitals* paraV, - elecstate::Potential* pot_in, - const K_Vectors& kv_in, - const TwoCenterBundle& two_center_bundle, - elecstate::DensityMatrix* DM_in, - int* exx_two_level_step) -{ + Gint_k* GK_in, + const Parallel_Orbitals* paraV, + elecstate::Potential* pot_in, + const K_Vectors& kv_in, + const TwoCenterBundle& two_center_bundle, + elecstate::DensityMatrix* DM_in +#ifdef __EXX + , int* exx_two_level_step + , std::vector>>>* Hexxd + , std::vector>>>>* Hexxc +#endif +) { this->kv = &kv_in; this->classname = "HamiltLCAO"; @@ -360,17 +363,19 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, #ifdef __EXX if (GlobalC::exx_info.info_global.cal_exx) { - Operator*exx = new OperatorEXX>(this->hsk, - LM_in, + // Peize Lin add 2016-12-03 + // set xc type before the first cal of xc in pelec->init_scf + // and calculate Cs, Vs + Operator* exx = new OperatorEXX>(this->hsk, this->hR, *this->kv, - LM_in->Hexxd, - LM_in->Hexxc, + Hexxd, + Hexxc, Add_Hexx_Type::R, exx_two_level_step, !GlobalC::restart.info_load.restart_exx && GlobalC::restart.info_load.load_H); - this->getOperator()->add(exx); + this->getOperator()->add(exx); } #endif // if NSPIN==2, HR should be separated into two parts, save HR into this->hRS2 diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 3e6d387df7..c568b2d447 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -5,11 +5,13 @@ #include "module_elecstate/module_dm/density_matrix.h" #include "module_elecstate/potentials/potential_new.h" #include "module_hamilt_general/hamilt.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" +#ifdef __EXX +#include "module_ri/Exx_LRI.h" +#endif namespace hamilt { @@ -23,15 +25,20 @@ class HamiltLCAO : public Hamilt * @brief Constructor of Hamiltonian for LCAO base * HR and SR will be allocated with Operators */ - HamiltLCAO(Gint_Gamma* GG_in, - Gint_k* GK_in, - LCAO_Matrix* LM_in, - const Parallel_Orbitals* paraV, - elecstate::Potential* pot_in, - const K_Vectors& kv_in, - const TwoCenterBundle& two_center_bundle, - elecstate::DensityMatrix* DM_in, - int* exx_two_level_step = nullptr); + using TAC = std::pair>; + HamiltLCAO(Gint_Gamma* GG_in, + Gint_k* GK_in, + const Parallel_Orbitals* paraV, + elecstate::Potential* pot_in, + const K_Vectors& kv_in, + const TwoCenterBundle& two_center_bundle, + elecstate::DensityMatrix* DM_in +#ifdef __EXX + , int* exx_two_level_step = nullptr + , std::vector>>>* Hexxd = nullptr + , std::vector>>>>* Hexxc = nullptr +#endif + ); /** * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer */ diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp index f4418840b5..56f9520b70 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp @@ -8,18 +8,17 @@ namespace hamilt template<> void OperatorEXX>::add_loaded_Hexx(const int ik) { - BlasConnector::axpy(this->LM->ParaV->get_local_size(), 1.0, this->LM->Hexxd_k_load[ik].data(), 1, this->hsk->get_hk(), 1); + BlasConnector::axpy(this->hR->get_paraV()->get_local_size(), 1.0, this->Hexxd_k_load[ik].data(), 1, this->hsk->get_hk(), 1); } template<> void OperatorEXX, double>>::add_loaded_Hexx(const int ik) { - - BlasConnector::axpy(this->LM->ParaV->get_local_size(), 1.0, this->LM->Hexxc_k_load[ik].data(), 1, this->hsk->get_hk(), 1); + BlasConnector::axpy(this->hR->get_paraV()->get_local_size(), 1.0, this->Hexxc_k_load[ik].data(), 1, this->hsk->get_hk(), 1); } template<> void OperatorEXX, std::complex>>::add_loaded_Hexx(const int ik) { - BlasConnector::axpy(this->LM->ParaV->get_local_size(), 1.0, this->LM->Hexxc_k_load[ik].data(), 1, this->hsk->get_hk(), 1); + BlasConnector::axpy(this->hR->get_paraV()->get_local_size(), 1.0, this->Hexxc_k_load[ik].data(), 1, this->hsk->get_hk(), 1); } } // namespace hamilt 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 3d9b1f9c90..afe27a7733 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 @@ -7,7 +7,6 @@ #include #include "operator_lcao.h" #include "module_cell/klist.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" namespace hamilt { @@ -28,7 +27,6 @@ class OperatorEXX> : public OperatorLCAO using TAC = std::pair>; public: OperatorEXX>(HS_Matrix_K* hsk_in, - LCAO_Matrix* LM_in, hamilt::HContainer* hR_in, const K_Vectors& kv_in, std::vector>>>* Hexxd_in = nullptr, @@ -57,14 +55,16 @@ class OperatorEXX> : public OperatorLCAO bool restart = false; void add_loaded_Hexx(const int ik); - void add_loaded_HexxR() {}; - void clear_loaded_HexxR() {}; + const K_Vectors& kv; - LCAO_Matrix* LM = nullptr; // if k points has no shift, use cell_nearest to reduce the memory cost RI::Cell_Nearest cell_nearest; bool use_cell_nearest = true; + + /// @brief Hexxk for all k-points, only for the 1st scf loop ofrestart load + std::vector> Hexxd_k_load; + std::vector>> Hexxc_k_load; }; } // namespace hamilt 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 90e0d0d968..1d3c029a6a 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 @@ -69,7 +69,6 @@ namespace hamilt template OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, - LCAO_Matrix* LM_in, HContainer*hR_in, const K_Vectors& kv_in, std::vector>>>* Hexxd_in, @@ -77,8 +76,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, Add_Hexx_Type add_hexx_type_in, int* two_level_step_in, const bool restart_in) - : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), - LM(LM_in), + : OperatorLCAO(hsk_in, kv_in.kvec_d, hR_in), kv(kv_in), Hexxd(Hexxd_in), Hexxc(Hexxc_in), @@ -88,6 +86,7 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, { ModuleBase::TITLE("OperatorEXX", "OperatorEXX"); this->cal_type = calculation_type::lcao_exx; + const Parallel_Orbitals* const pv = hR_in->get_paraV(); if (GlobalV::CALCULATION == "nscf") { // if nscf, read HexxR first and reallocate hR according to the read-in HexxR @@ -142,25 +141,25 @@ OperatorEXX>::OperatorEXX(HS_Matrix_K* hsk_in, /// read in Hexx(k) if (std::is_same::value) { - this->LM->Hexxd_k_load.resize(this->kv.get_nks()); + this->Hexxd_k_load.resize(this->kv.get_nks()); for (int ik = 0; ik < this->kv.get_nks(); ik++) { - this->LM->Hexxd_k_load[ik].resize(this->LM->ParaV->get_local_size(), 0.0); + this->Hexxd_k_load[ik].resize(pv->get_local_size(), 0.0); this->restart = GlobalC::restart.load_disk( "Hexx", ik, - this->LM->ParaV->get_local_size(), this->LM->Hexxd_k_load[ik].data(), false); + pv->get_local_size(), this->Hexxd_k_load[ik].data(), false); if (!this->restart) break; } } else { - this->LM->Hexxc_k_load.resize(this->kv.get_nks()); + this->Hexxc_k_load.resize(this->kv.get_nks()); for (int ik = 0; ik < this->kv.get_nks(); ik++) { - this->LM->Hexxc_k_load[ik].resize(this->LM->ParaV->get_local_size(), 0.0); + this->Hexxc_k_load[ik].resize(pv->get_local_size(), 0.0); this->restart = GlobalC::restart.load_disk( "Hexx", ik, - this->LM->ParaV->get_local_size(), this->LM->Hexxc_k_load[ik].data(), false); + pv->get_local_size(), this->Hexxc_k_load[ik].data(), false); if (!this->restart) break; } } @@ -198,8 +197,8 @@ void OperatorEXX>::contributeHR() RI_2D_Comm::add_HexxR( this->current_spin, GlobalC::exx_info.info_global.hybrid_alpha, - this->Hexxd == nullptr ? *this->LM->Hexxd : *this->Hexxd, - *this->LM->ParaV, + *this->Hexxd, + *this->hR->get_paraV(), GlobalV::NPOL, *this->hR, this->use_cell_nearest ? &this->cell_nearest : nullptr); @@ -207,8 +206,8 @@ void OperatorEXX>::contributeHR() RI_2D_Comm::add_HexxR( this->current_spin, GlobalC::exx_info.info_global.hybrid_alpha, - this->Hexxc == nullptr ? *this->LM->Hexxc : *this->Hexxc, - *this->LM->ParaV, + *this->Hexxc, + *this->hR->get_paraV(), GlobalV::NPOL, *this->hR, this->use_cell_nearest ? &this->cell_nearest : nullptr); @@ -233,15 +232,15 @@ void OperatorEXX>::contributeHk(int ik) } else // clear loaded Hexx and release memory { - if (this->LM->Hexxd_k_load.size() > 0) + if (this->Hexxd_k_load.size() > 0) { - this->LM->Hexxd_k_load.clear(); - this->LM->Hexxd_k_load.shrink_to_fit(); + this->Hexxd_k_load.clear(); + this->Hexxd_k_load.shrink_to_fit(); } - else if (this->LM->Hexxc_k_load.size() > 0) + else if (this->Hexxc_k_load.size() > 0) { - this->LM->Hexxc_k_load.clear(); - this->LM->Hexxc_k_load.shrink_to_fit(); + this->Hexxc_k_load.clear(); + this->Hexxc_k_load.shrink_to_fit(); } } } @@ -252,16 +251,16 @@ void OperatorEXX>::contributeHk(int ik) this->kv, ik, GlobalC::exx_info.info_global.hybrid_alpha, - this->Hexxd == nullptr ? *this->LM->Hexxd : *this->Hexxd, - *this->LM->ParaV, + *this->Hexxd, + *this->hR->get_paraV(), this->hsk->get_hk()); else RI_2D_Comm::add_Hexx( this->kv, ik, GlobalC::exx_info.info_global.hybrid_alpha, - this->Hexxc == nullptr ? *this->LM->Hexxc : *this->Hexxc, - *this->LM->ParaV, + *this->Hexxc, + *this->hR->get_paraV(), this->hsk->get_hk()); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index 5e089ef8a9..6420b36d5e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -7,13 +7,17 @@ #include "spar_u.h" void sparse_format::cal_HSR(const Parallel_Orbitals& pv, - LCAO_Matrix& lm, - LCAO_HS_Arrays& HS_Arrays, - Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt>* p_ham) { + LCAO_HS_Arrays& HS_Arrays, + Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int(&nmp)[3], + hamilt::Hamilt>* p_ham +#ifdef __EXX + , const std::vector>>>* Hexxd + , const std::vector>>>>* Hexxc +#endif +) { ModuleBase::TITLE("sparse_format", "cal_HSR"); sparse_format::set_R_range(HS_Arrays.all_R_coor, grid); @@ -94,18 +98,18 @@ void sparse_format::cal_HSR(const Parallel_Orbitals& pv, if (GlobalC::exx_info.info_global.cal_exx) { if (GlobalC::exx_info.info_ri.real_number) { sparse_format::cal_HR_exx(pv, - HS_Arrays, - current_spin, - sparse_thr, - nmp, - *lm.Hexxd); + HS_Arrays, + current_spin, + sparse_thr, + nmp, + *Hexxd); } else { sparse_format::cal_HR_exx(pv, - HS_Arrays, - current_spin, - sparse_thr, - nmp, - *lm.Hexxc); + HS_Arrays, + current_spin, + sparse_thr, + nmp, + *Hexxc); } } #endif // __MPI diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index e4c438763b..497586574e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -1,19 +1,22 @@ #ifndef SPARSE_FORMAT_HSR_H #define SPARSE_FORMAT_HSR_H -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" namespace sparse_format { - -void cal_HSR(const Parallel_Orbitals& pv, - LCAO_Matrix& lm, - LCAO_HS_Arrays& HS_Arrays, - Grid_Driver& grid, - const int& current_spin, - const double& sparse_thr, - const int (&nmp)[3], - hamilt::Hamilt>* p_ham); + using TAC = std::pair>; + void cal_HSR(const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + Grid_Driver& grid, + const int& current_spin, + const double& sparse_thr, + const int(&nmp)[3], + hamilt::Hamilt>* p_ham +#ifdef __EXX + , const std::vector>>>* Hexxd = nullptr + , const std::vector>>>>* Hexxc = nullptr +#endif + ); void cal_HContainer_d( const Parallel_Orbitals& pv, diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index 9d6bf53b0b..0b6be39abc 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -9,24 +9,21 @@ namespace ModuleIO template Output_Mat_Sparse::Output_Mat_Sparse(int out_mat_hsR, - int out_mat_dh, - int out_mat_t, - int out_mat_r, - int istep, - const ModuleBase::matrix& v_eff, - const Parallel_Orbitals& pv, - Gint_k& gint_k, // mohan add 2024-04-01 - const TwoCenterBundle& two_center_bundle, - LCAO_Matrix& lm, - Grid_Driver& grid, // mohan add 2024-04-06 - const K_Vectors& kv, - hamilt::Hamilt* p_ham) + int out_mat_dh, + int out_mat_t, + int out_mat_r, + int istep, + const ModuleBase::matrix& v_eff, + const Parallel_Orbitals& pv, + Gint_k& gint_k, // mohan add 2024-04-01 + const TwoCenterBundle& two_center_bundle, + Grid_Driver& grid, // mohan add 2024-04-06 + const K_Vectors& kv, + hamilt::Hamilt* p_ham) : _out_mat_hsR(out_mat_hsR), _out_mat_dh(out_mat_dh), _out_mat_t(out_mat_t), _out_mat_r(out_mat_r), _istep(istep), - _v_eff(v_eff), _pv(pv), _gint_k(gint_k), // mohan add 2024-04-01 - two_center_bundle_(two_center_bundle), _lm(lm), _grid(grid), // mohan add 2024-04-06 - _kv(kv), _p_ham(p_ham) -{ -} + _v_eff(v_eff), _pv(pv), _gint_k(gint_k), // mohan add 2024-04-01 + two_center_bundle_(two_center_bundle), _grid(grid), // mohan add 2024-04-06 + _kv(kv), _p_ham(p_ham) {} template <> void Output_Mat_Sparse::write() @@ -41,7 +38,7 @@ void Output_Mat_Sparse>::write() //! generate a file containing the Hamiltonian and S(overlap) matrices if (_out_mat_hsR) { - output_HSR(_istep, this->_v_eff, this->_pv, this->_lm, HS_Arrays, this->_grid, _kv, _p_ham); + output_HSR(_istep, this->_v_eff, this->_pv, HS_Arrays, this->_grid, _kv, _p_ham); } //! generate a file containing the kinetic energy matrix diff --git a/source/module_io/output_mat_sparse.h b/source/module_io/output_mat_sparse.h index 9419b35176..4b5bd4498b 100644 --- a/source/module_io/output_mat_sparse.h +++ b/source/module_io/output_mat_sparse.h @@ -3,7 +3,6 @@ #include "module_basis/module_ao/parallel_orbitals.h" #include "module_basis/module_nao/two_center_bundle.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hsolver/hsolver_lcao.h" @@ -15,18 +14,17 @@ class Output_Mat_Sparse { public: Output_Mat_Sparse(int out_mat_hsR, - int out_mat_dh, - int out_mat_t, - int out_mat_r, - int istep, - const ModuleBase::matrix& v_eff, - const Parallel_Orbitals& pv, - Gint_k& gint_k, // mohan add 2024-04-01 - const TwoCenterBundle& two_center_bundle, - LCAO_Matrix& lm, - Grid_Driver& grid, // mohan add 2024-04-06 - const K_Vectors& kv, - hamilt::Hamilt* p_ham); + int out_mat_dh, + int out_mat_t, + int out_mat_r, + int istep, + const ModuleBase::matrix& v_eff, + const Parallel_Orbitals& pv, + Gint_k& gint_k, // mohan add 2024-04-01 + const TwoCenterBundle& two_center_bundle, + Grid_Driver& grid, // mohan add 2024-04-06 + const K_Vectors& kv, + hamilt::Hamilt* p_ham); void write(); @@ -54,8 +52,6 @@ class Output_Mat_Sparse // introduced temporarily for LCAO refactoring const TwoCenterBundle& two_center_bundle_; - LCAO_Matrix& _lm; - // mohan fix bug 2024-04-07, a typical bug!!! Grid_Driver& _grid; // mohan add 2024-04-06 diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 12dc6728bf..380f770421 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -12,18 +12,22 @@ // 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, - const ModuleBase::matrix& v_eff, - const Parallel_Orbitals& pv, - LCAO_Matrix& lm, - LCAO_HS_Arrays& HS_Arrays, - Grid_Driver& grid, // mohan add 2024-04-06 - const K_Vectors& kv, - hamilt::Hamilt>* p_ham, - const std::string& SR_filename, - const std::string& HR_filename_up, - const std::string HR_filename_down, - const bool& binary, - const double& sparse_thr) { + const ModuleBase::matrix& v_eff, + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + Grid_Driver& grid, // mohan add 2024-04-06 + const K_Vectors& kv, + hamilt::Hamilt>* p_ham, +#ifdef __EXX + const std::vector>>>* Hexxd, + const std::vector>>>>* Hexxc, +#endif + const std::string& SR_filename, + const std::string& HR_filename_up, + const std::string HR_filename_down, + const bool& binary, + const double& sparse_thr +) { ModuleBase::TITLE("ModuleIO", "output_HSR"); ModuleBase::timer::tick("ModuleIO", "output_HSR"); @@ -32,26 +36,21 @@ 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, - lm, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham); - } else if (nspin == 2) { + sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham +#ifdef __EXX + , Hexxd, Hexxc +#endif + ); + } + else if (nspin == 2) { int spin_now = 1; // save HR of spin down first (the current spin always be down) - sparse_format::cal_HSR(pv, - lm, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham); + sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham +#ifdef __EXX + , Hexxd, Hexxc +#endif + ); // cal HR of the spin up if (GlobalV::VL_IN_H) { @@ -61,14 +60,11 @@ void ModuleIO::output_HSR(const int& istep, spin_now = 0; } - sparse_format::cal_HSR(pv, - lm, - HS_Arrays, - grid, - spin_now, - sparse_thr, - kv.nmp, - p_ham); + sparse_format::cal_HSR(pv, HS_Arrays, grid, spin_now, sparse_thr, kv.nmp, p_ham +#ifdef __EXX + , Hexxd, Hexxc +#endif + ); } ModuleIO::save_HSR_sparse(istep, diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h index d9b6c21b60..f694b282d0 100644 --- a/source/module_io/write_HS_R.h +++ b/source/module_io/write_HS_R.h @@ -5,25 +5,28 @@ #include "module_basis/module_nao/two_center_bundle.h" #include "module_cell/klist.h" #include "module_hamilt_general/hamilt.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" namespace ModuleIO { + using TAC = std::pair>; void output_HSR(const int& istep, - const ModuleBase::matrix& v_eff, - const Parallel_Orbitals& pv, - LCAO_Matrix& lm, - LCAO_HS_Arrays& HS_Arrays, - Grid_Driver& grid, // mohan add 2024-04-06 - const K_Vectors& kv, - hamilt::Hamilt>* p_ham, - const std::string& SR_filename = "data-SR-sparse_SPIN0.csr", - const std::string& HR_filename_up = "data-HR-sparse_SPIN0.csr", - const std::string HR_filename_down = "data-HR-sparse_SPIN1.csr", - const bool& binary = false, - const double& sparse_threshold = 1e-10); // LiuXh add 2019-07-15, modify in 2021-12-3 + const ModuleBase::matrix& v_eff, + const Parallel_Orbitals& pv, + LCAO_HS_Arrays& HS_Arrays, + Grid_Driver& grid, // mohan add 2024-04-06 + const K_Vectors& kv, + hamilt::Hamilt>* p_ham, +#ifdef __EXX + const std::vector>>>* Hexxd = nullptr, + const std::vector>>>>* Hexxc = nullptr, +#endif + const std::string& SR_filename = "data-SR-sparse_SPIN0.csr", + const std::string& HR_filename_up = "data-HR-sparse_SPIN0.csr", + const std::string HR_filename_down = "data-HR-sparse_SPIN1.csr", + const bool& binary = false, + const double& sparse_threshold = 1e-10); // LiuXh add 2019-07-15, modify in 2021-12-3 void output_dHR(const int& istep, const ModuleBase::matrix& v_eff, diff --git a/source/module_io/write_Vxc.hpp b/source/module_io/write_Vxc.hpp index a87ef6f276..222aeadcd2 100644 --- a/source/module_io/write_Vxc.hpp +++ b/source/module_io/write_Vxc.hpp @@ -270,24 +270,28 @@ void set_gint_pointer>(Gint_Gamma& gint_gamma, /// including terms: local/semi-local XC, EXX, DFTU template void write_Vxc(int nspin, - int nbasis, - int drank, - const psi::Psi& psi, - const UnitCell& ucell, - Structure_Factor& sf, - const ModulePW::PW_Basis& rho_basis, - const ModulePW::PW_Basis& rhod_basis, - const ModuleBase::matrix& vloc, - const Charge& chg, - Gint_Gamma& gint_gamma, // mohan add 2024-04-01 - Gint_k& gint_k, // mohan add 2024-04-01 - LCAO_Matrix& lm, - const K_Vectors& kv, - const ModuleBase::matrix& wg, - Grid_Driver& gd) + int nbasis, + int drank, + const Parallel_Orbitals* pv, + const psi::Psi& psi, + const UnitCell& ucell, + Structure_Factor& sf, + const ModulePW::PW_Basis& rho_basis, + const ModulePW::PW_Basis& rhod_basis, + const ModuleBase::matrix& vloc, + const Charge& chg, + Gint_Gamma& gint_gamma, // mohan add 2024-04-01 + Gint_k& gint_k, // mohan add 2024-04-01 + const K_Vectors& kv, + const ModuleBase::matrix& wg, + Grid_Driver& gd +#ifdef __EXX + , std::vector>>>* Hexxd = nullptr + , std::vector>>>>* Hexxc = nullptr +#endif +) { ModuleBase::TITLE("ModuleIO", "write_Vxc"); - const Parallel_Orbitals* pv = lm.ParaV; int nbands = wg.nc; // 1. real-space xc potential // ModuleBase::matrix vr_xc(nspin, chg.nrxx); @@ -332,9 +336,9 @@ void write_Vxc(int nspin, std::vector> e_orb_locxc; // orbital energy (local XC) std::vector> e_orb_tot; // orbital energy (total) #ifdef __EXX - hamilt::OperatorEXX> vexx_op_ao(&vxc_k_ao, &lm, nullptr, kv, lm.Hexxd, lm.Hexxc, hamilt::Add_Hexx_Type::k); + hamilt::OperatorEXX> vexx_op_ao(&vxc_k_ao, nullptr, 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, &lm, nullptr, kv, lm.Hexxd, lm.Hexxc, hamilt::Add_Hexx_Type::k); + hamilt::OperatorEXX> vexxonly_op_ao(&vexxonly_k_ao, nullptr, 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_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 7f43274fef..1178a1edc0 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -124,8 +124,8 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol this->nbasis = GlobalV::NLOCAL; this->nstates = inp.lr_nstates; LR_Util::setup_2d_division(this->paraMat_, 1, this->nbasis, this->nbasis); - this->paraMat_.atom_begin_row = std::move(ks_sol.LM.ParaV->atom_begin_row); - this->paraMat_.atom_begin_col = std::move(ks_sol.LM.ParaV->atom_begin_col); + this->paraMat_.atom_begin_row = std::move(ks_sol.ParaV.atom_begin_row); + this->paraMat_.atom_begin_col = std::move(ks_sol.ParaV.atom_begin_col); this->paraMat_.iat2iwt_ = ucell.get_iat2iwt(); // move the ground state info