From 7d8fce17ac955a815eab2eb1761f2bd30e0b883f Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 11 Dec 2024 16:19:30 +0800 Subject: [PATCH 1/4] Refactor: dd TITLE and timer::tick for esolver --- source/module_esolver/esolver_fp.cpp | 2 ++ source/module_esolver/esolver_ks.cpp | 6 ++---- source/module_esolver/esolver_ks_lcao.cpp | 4 ++++ source/module_esolver/esolver_ks_lcao_tddft.cpp | 5 +++++ source/module_esolver/esolver_ks_pw.cpp | 8 ++++++++ source/module_esolver/esolver_of.cpp | 10 ++++++++++ source/module_esolver/esolver_sdft_pw.cpp | 10 ++++++++++ source/module_esolver/lcao_before_scf.cpp | 3 +++ 8 files changed, 44 insertions(+), 4 deletions(-) diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index ba5770414f..94daacd37d 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -129,6 +129,8 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp) //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. void ESolver_FP::after_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_FP", "after_scf"); + // 0) output convergence information ModuleIO::output_convergence_after_scf(this->conv_esolver, this->pelec->f_en.etot); diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index a79ab77493..e644d19e0a 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -427,9 +427,7 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) ModuleBase::timer::tick(this->classname, "runner"); // 2) before_scf (electronic iteration loops) - ModuleBase::timer::tick(this->classname, "before_scf"); this->before_scf(ucell, istep); - ModuleBase::timer::tick(this->classname, "before_scf"); // 3) write charge density if (PARAM.inp.dm_to_rho) @@ -468,9 +466,7 @@ void ESolver_KS::runner(UnitCell& ucell, const int istep) } // end scf iterations // 9) after scf - ModuleBase::timer::tick(this->classname, "after_scf"); this->after_scf(ucell, istep); - ModuleBase::timer::tick(this->classname, "after_scf"); ModuleBase::timer::tick(this->classname, "runner"); return; @@ -740,6 +736,8 @@ void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& i template void ESolver_KS::after_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_KS", "after_scf"); + // 1) call after_scf() of ESolver_FP ESolver_FP::after_scf(ucell, istep); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 63a1201bef..fce56a9d01 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -903,6 +903,8 @@ template void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); + // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 if (PARAM.inp.out_elf[0] > 0) { @@ -1266,6 +1268,8 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); } + + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); } template class ESolver_KS_LCAO; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 880daa8512..525ff5f73f 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -281,6 +281,9 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const i void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_KS_LCAO_TDDFT", "after_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO_TDDFT", "after_scf"); + for (int is = 0; is < PARAM.inp.nspin; is++) { if (module_tddft::Evolve_elec::out_dipole == 1) @@ -306,6 +309,8 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) this->RA); } ESolver_KS_LCAO, double>::after_scf(ucell, istep); + + ModuleBase::timer::tick("ESolver_KS_LCAO_TDDFT", "after_scf"); } void ESolver_KS_LCAO_TDDFT::weight_dm_rho() diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index d57c4f5ffc..385cf27621 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -247,6 +247,7 @@ template void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_PW", "before_scf"); + ModuleBase::timer::tick("ESolver_KS_PW", "before_scf"); //! 1) call before_scf() of ESolver_KS ESolver_KS::before_scf(ucell, istep); @@ -427,6 +428,8 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) this->already_initpsi = true; } } + + ModuleBase::timer::tick("ESolver_KS_PW", "before_scf"); } template @@ -647,6 +650,9 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int template void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_KS_PW", "after_scf"); + ModuleBase::timer::tick("ESolver_KS_PW", "after_scf"); + // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 if (PARAM.inp.out_elf[0] > 0) { @@ -741,6 +747,8 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) auto* onsite_p = projectors::OnsiteProjector::get_instance(); onsite_p->cal_occupations(reinterpret_cast, Device>*>(this->kspw_psi), this->pelec->wg); } + + ModuleBase::timer::tick("ESolver_KS_PW", "after_scf"); } template diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 32d4a90b2e..bdb24829e4 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/source/module_esolver/esolver_of.cpp @@ -200,6 +200,9 @@ void ESolver_OF::runner(UnitCell& ucell, const int istep) */ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) { + ModuleBase::TITLE("ESolver_OF", "before_opt"); + ModuleBase::timer::tick("ESolver_OF", "before_opt"); + //! 1) call before_scf() of ESolver_FP ESolver_FP::before_scf(ucell, istep); @@ -299,6 +302,8 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell) { this->theta_[0] = 0.2; } + + ModuleBase::timer::tick("ESolver_OF", "before_opt"); } /** @@ -483,6 +488,9 @@ bool ESolver_OF::check_exit() */ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) { + ModuleBase::TITLE("ESolver_OF", "after_opt"); + ModuleBase::timer::tick("ESolver_OF", "after_opt"); + // 1) calculate the kinetic energy density if (PARAM.inp.out_elf[0] > 0) { @@ -491,6 +499,8 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell) // 2) call after_scf() of ESolver_FP ESolver_FP::after_scf(ucell, istep); + + ModuleBase::timer::tick("ESolver_OF", "after_opt"); } /** diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index fc18e8e302..abeda884b1 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -94,6 +94,9 @@ void ESolver_SDFT_PW::before_all_runners(UnitCell& ucell, const Input template void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_SDFT_PW", "before_scf"); + ModuleBase::timer::tick("ESolver_SDFT_PW", "before_scf"); + ESolver_KS_PW::before_scf(ucell, istep); delete reinterpret_cast*>(this->p_hamilt); this->p_hamilt = new hamilt::HamiltSdftPW(this->pelec->pot, @@ -110,6 +113,8 @@ void ESolver_SDFT_PW::before_scf(UnitCell& ucell, const int istep) { this->stowf.update_sto_orbitals(PARAM.inp.seed_sto); } + + ModuleBase::timer::tick("ESolver_SDFT_PW", "before_scf"); } template @@ -122,8 +127,13 @@ void ESolver_SDFT_PW::iter_finish(UnitCell& ucell, const int istep, i template void ESolver_SDFT_PW::after_scf(UnitCell& ucell, const int istep) { + ModuleBase::TITLE("ESolver_SDFT_PW", "after_scf"); + ModuleBase::timer::tick("ESolver_SDFT_PW", "after_scf"); + // 1) call after_scf() of ESolver_KS_PW ESolver_KS_PW::after_scf(ucell, istep); + + ModuleBase::timer::tick("ESolver_SDFT_PW", "after_scf"); } template diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 2066b5069b..1cecc6e38c 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -43,6 +43,7 @@ template void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_KS_LCAO", "before_scf"); + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); //! 1) call before_scf() of ESolver_KS ESolver_KS::before_scf(ucell, istep); @@ -347,6 +348,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) 1); } + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); return; } @@ -377,6 +379,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 } + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); return; } From 5eea6d8e6fc83f1b529186f7ef216485589459f2 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 11 Dec 2024 18:47:47 +0800 Subject: [PATCH 2/4] Refactor: remove GlobalC::GridD --- source/module_base/tool_title.cpp | 16 +- .../module_neighbor/sltk_grid_driver.cpp | 4 - .../module_neighbor/sltk_grid_driver.h | 5 - .../module_dm/test/test_dm_R_init.cpp | 2 +- source/module_esolver/esolver_gets.cpp | 14 +- source/module_esolver/esolver_ks_lcao.cpp | 61 +++-- source/module_esolver/esolver_ks_lcao.h | 2 + .../module_esolver/esolver_ks_lcao_tddft.cpp | 1 + source/module_esolver/lcao_before_scf.cpp | 16 +- source/module_esolver/lcao_others.cpp | 22 +- .../module_hamilt_lcao/hamilt_lcaodft/FORCE.h | 2 + .../hamilt_lcaodft/FORCE_STRESS.cpp | 102 ++++---- .../hamilt_lcaodft/FORCE_STRESS.h | 2 + .../hamilt_lcaodft/FORCE_gamma.cpp | 35 +-- .../hamilt_lcaodft/FORCE_k.cpp | 40 ++-- .../hamilt_lcaodft/hamilt_lcao.cpp | 2 +- .../operator_lcao/deepks_lcao.cpp | 7 +- .../operator_lcao/deepks_lcao.h | 2 + .../hamilt_lcaodft/spar_dh.cpp | 25 +- .../hamilt_lcaodft/spar_st.cpp | 2 +- source/module_hamilt_lcao/module_dftu/dftu.h | 86 ++++--- .../module_dftu/dftu_folding.cpp | 56 +++-- .../module_dftu/dftu_force.cpp | 20 +- source/module_io/berryphase.cpp | 8 +- source/module_io/berryphase.h | 6 +- source/module_io/cal_r_overlap_R.cpp | 14 +- source/module_io/cal_r_overlap_R.h | 6 +- source/module_io/output_mat_sparse.cpp | 2 +- source/module_io/output_mulliken.h | 33 ++- source/module_io/td_current_io.cpp | 3 +- source/module_io/td_current_io.h | 1 + source/module_io/to_wannier90_lcao.cpp | 19 +- source/module_io/to_wannier90_lcao.h | 4 +- source/module_io/unk_overlap_lcao.cpp | 18 +- source/module_io/unk_overlap_lcao.h | 2 +- source/module_lr/esolver_lrtd_lcao.cpp | 124 +++++++--- source/module_lr/esolver_lrtd_lcao.h | 1 + .../operator_casida/operator_lr_hxc.cpp | 2 +- source/module_rdmft/rdmft.cpp | 18 +- source/module_rdmft/rdmft.h | 16 +- source/module_rdmft/rdmft_pot.cpp | 226 ++++++++---------- source/module_rdmft/update_state_rdmft.cpp | 4 +- 42 files changed, 544 insertions(+), 487 deletions(-) diff --git a/source/module_base/tool_title.cpp b/source/module_base/tool_title.cpp index e55f966f27..922d144a42 100644 --- a/source/module_base/tool_title.cpp +++ b/source/module_base/tool_title.cpp @@ -17,10 +17,10 @@ namespace ModuleBase void TITLE(const std::string &class_name,const std::string &function_name,const bool disable) { - if (disable) - { - return;//no output - } + // if (disable) + // { + // return;//no output + // } #ifdef __NORMAL std::cout<<" ==> "< "< #endif -namespace GlobalC -{ -Grid_Driver GridD; -} Grid_Driver::Grid_Driver( const int &test_d_in, const int &test_grid_in) diff --git a/source/module_cell/module_neighbor/sltk_grid_driver.h b/source/module_cell/module_neighbor/sltk_grid_driver.h index fbc683cade..32ccf482af 100644 --- a/source/module_cell/module_neighbor/sltk_grid_driver.h +++ b/source/module_cell/module_neighbor/sltk_grid_driver.h @@ -125,9 +125,4 @@ class Grid_Driver : public Grid const short box_y, const short box_z) const; }; - -namespace GlobalC -{ -extern Grid_Driver GridD; -} #endif diff --git a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp index ac2d7161b1..ef4ccef8ce 100644 --- a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp +++ b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp @@ -100,7 +100,7 @@ class DMTest : public testing::Test #endif }; -// test for construct DMR from GlobalC::GridD and UnitCell +// test for construct DMR from GridD and UnitCell TEST_F(DMTest, DMInit1) { // initalize a kvectors diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 3cdfcde99a..cf51f6bbbe 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -94,15 +94,17 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) ucell.infoNL.get_rcutmax_Beta(), PARAM.globalv.gamma_only_local); + Grid_Driver gd; + atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, - GlobalC::GridD, + gd, ucell, search_radius, PARAM.inp.test_atom_input); Record_adj RA; - RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + RA.for_2d(ucell, gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); if (this->p_hamilt == nullptr) { @@ -110,7 +112,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) { this->p_hamilt = new hamilt::HamiltLCAO, std::complex>(ucell, - GlobalC::GridD, + gd, &this->pv, this->kv, *(two_center_bundle_.overlap_orb), @@ -121,7 +123,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) else { this->p_hamilt = new hamilt::HamiltLCAO, double>(ucell, - GlobalC::GridD, + gd, &this->pv, this->kv, *(two_center_bundle_.overlap_orb), @@ -132,13 +134,13 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) const std::string fn = PARAM.globalv.global_out_dir + "SR.csr"; std::cout << " The file is saved in " << fn << std::endl; - ModuleIO::output_SR(pv, GlobalC::GridD, this->p_hamilt, fn); + ModuleIO::output_SR(pv, gd, this->p_hamilt, fn); if (PARAM.inp.out_mat_r) { cal_r_overlap_R r_matrix; r_matrix.init(ucell,pv, orb_); - r_matrix.out_rR(ucell,istep); + r_matrix.out_rR(ucell, gd, istep); } ModuleBase::timer::tick("ESolver_GetS", "runner"); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index fce56a9d01..84a5773285 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -259,8 +259,17 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // 14) initialize rdmft, added by jghan if( PARAM.inp.rdmft == true ) { - rdmft_solver.init( this->GG, this->GK, this->pv, ucell, this->kv, *(this->pelec), - this->orb_, two_center_bundle_, PARAM.inp.dft_functional, PARAM.inp.rdmft_power_alpha); + rdmft_solver.init(this->GG, + this->GK, + this->pv, + ucell, + this->gd, + this->kv, + *(this->pelec), + this->orb_, + two_center_bundle_, + PARAM.inp.dft_functional, + PARAM.inp.rdmft_power_alpha); } ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners"); @@ -296,6 +305,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for PARAM.inp.test_force, PARAM.inp.test_stress, ucell, + this->gd, this->pv, this->pelec, this->psi, @@ -458,7 +468,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) this->kv, orb_.cutoffs(), this->pelec->wg, - GlobalC::GridD + this->gd #ifdef __EXX , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr, @@ -484,7 +494,7 @@ void ESolver_KS_LCAO::after_all_runners(UnitCell& ucell) this->GK, this->kv, this->pelec->wg, - GlobalC::GridD, + this->gd, orb_.cutoffs(), this->two_center_bundle_ #ifdef __EXX @@ -1036,7 +1046,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) this->pelec->klist->kvec_d, ucell, orb_, - GlobalC::GridD, + this->gd, &(this->pv), *(this->psi), dynamic_cast*>(this->pelec)->get_DM(), @@ -1136,14 +1146,23 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) two_center_bundle_, orb_, ucell, - GlobalC::GridD, + this->gd, this->kv, this->p_hamilt); - + //! Perform Mulliken charge analysis if (PARAM.inp.out_mul) { - ModuleIO::cal_mag(&(this->pv), this->p_hamilt, this->kv, this->pelec, this->two_center_bundle_, this->orb_, ucell, istep, true); + ModuleIO::cal_mag(&(this->pv), + this->p_hamilt, + this->kv, + this->pelec, + this->two_center_bundle_, + this->orb_, + ucell, + this->gd, + istep, + true); } } @@ -1189,7 +1208,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) &hR, &ucell, orb_.cutoffs(), - &GlobalC::GridD, + &this->gd, two_center_bundle_.kinetic_orb.get()); const int nspin_k = (PARAM.inp.nspin == 2 ? 2 : 1); @@ -1228,12 +1247,12 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.wannier_spin); myWannier.calculate(ucell, - this->pelec->ekb, - this->pw_wfc, - this->pw_big, - this->sf, - this->kv, - this->psi, + this->pelec->ekb, + this->pw_wfc, + this->pw_big, + this->sf, + this->kv, + this->psi, &(this->pv)); } else if (PARAM.inp.wannier_method == 2) @@ -1247,7 +1266,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.wannier_spin, orb_); - myWannier.calculate(ucell,this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); + myWannier.calculate(ucell, this->gd, this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); } @@ -1259,12 +1278,10 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation"); berryphase bp(&(this->pv)); - bp.lcao_init(ucell, - this->kv, - this->GridT, - orb_); // additional step before calling - // macroscopic_polarization (why capitalize - // the function name?) + bp.lcao_init(ucell, this->gd, this->kv, this->GridT, orb_); + // additional step before calling + // macroscopic_polarization (why capitalize + // the function name?) bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); } diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 9c7e8be612..1730f04a36 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -58,6 +58,8 @@ class ESolver_KS_LCAO : public ESolver_KS { // we will get rid of this class soon, don't use it, mohan 2024-03-28 Record_adj RA; + Grid_Driver gd; + // 2d block-cyclic distribution info Parallel_Orbitals pv; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 525ff5f73f..fc1ba66bc6 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -299,6 +299,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) = dynamic_cast>*>(this->pelec)->get_DM(); ModuleIO::write_current(ucell, + this->gd, istep, this->psi, pelec, diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 1cecc6e38c..9eef0aafce 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -81,7 +81,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, - GlobalC::GridD, + this->gd, ucell, search_radius, PARAM.inp.test_atom_input); @@ -111,7 +111,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->pw_rho->nplane, this->pw_rho->startz_current, ucell, - GlobalC::GridD, + this->gd, dr_uniform, rcuts, psi_u, @@ -128,7 +128,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // (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. - this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -190,7 +190,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), ucell, - GlobalC::GridD, + this->gd, &this->pv, this->pelec->pot, this->kv, @@ -214,15 +214,11 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { const Parallel_Orbitals* pv = &this->pv; // build and save at beginning - GlobalC::ld.build_psialpha(PARAM.inp.cal_force, - ucell, - orb_, - GlobalC::GridD, - *(two_center_bundle_.overlap_orb_alpha)); + GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha)); if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, GlobalC::GridD); + GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd); } } #endif diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index fa82e38ff6..55ee33c593 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -65,7 +65,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) double search_radius = PARAM.inp.search_radius; atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, - GlobalC::GridD, + this->gd, ucell, search_radius, PARAM.inp.test_atom_input, @@ -84,7 +84,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) atom_arrange::search(PARAM.inp.search_pbc, GlobalV::ofs_running, - GlobalC::GridD, + this->gd, ucell, search_radius, PARAM.inp.test_atom_input); @@ -114,7 +114,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->pw_rho->nplane, this->pw_rho->startz_current, ucell, - GlobalC::GridD, + this->gd, dr_uniform, rcuts, psi_u, @@ -131,7 +131,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // (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. - this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell, this->gd, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -193,7 +193,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), ucell, - GlobalC::GridD, + this->gd, &this->pv, this->pelec->pot, this->kv, @@ -217,15 +217,11 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) { const Parallel_Orbitals* pv = &this->pv; // build and save at beginning - GlobalC::ld.build_psialpha(PARAM.inp.cal_force, - ucell, - orb_, - GlobalC::GridD, - *(two_center_bundle_.overlap_orb_alpha)); + GlobalC::ld.build_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, *(two_center_bundle_.overlap_orb_alpha)); if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, GlobalC::GridD); + GlobalC::ld.check_psialpha(PARAM.inp.cal_force, ucell, orb_, this->gd); } } #endif @@ -283,7 +279,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) PARAM.globalv.global_out_dir, GlobalV::ofs_warning, &ucell, - &GlobalC::GridD, + &this->gd, this->kv); } else @@ -312,7 +308,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) PARAM.globalv.global_out_dir, GlobalV::ofs_warning, &ucell, - &GlobalC::GridD, + &this->gd, this->kv, PARAM.inp.if_separate_k, &GlobalC::Pgrid, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h index 171c94a5d8..3f6fc8f01d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h @@ -51,6 +51,7 @@ class Force_LCAO const bool isstress, ForceStressArrays& fsr, // mohan add 2024-06-16 const UnitCell& ucell, + Grid_Driver& gd, const psi::Psi* psi, const elecstate::ElecState* pelec, ModuleBase::matrix& foverlap, @@ -73,6 +74,7 @@ class Force_LCAO // get the ds, dt, dvnl. void allocate(const UnitCell& ucell, + Grid_Driver& gd, const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 43a633daf2..f9dfeed7bb 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -36,6 +36,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const bool istestf, const bool istests, const UnitCell& ucell, + Grid_Driver& gd, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -158,6 +159,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, isforce, isstress, ucell, + gd, fsr, pelec, psi, @@ -181,39 +183,36 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // calculate force and stress for Nonlocal part if(PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { - hamilt::NonlocalNew> tmp_nonlocal( - nullptr, - kv.kvec_d, - nullptr, - &ucell, - orb.cutoffs(), - &GlobalC::GridD, - two_center_bundle.overlap_orb_beta.get() - ); + hamilt::NonlocalNew> tmp_nonlocal(nullptr, + kv.kvec_d, + nullptr, + &ucell, + orb.cutoffs(), + &gd, + two_center_bundle.overlap_orb_beta.get()); const auto* dm_p = dynamic_cast*>(pelec)->get_DM(); - if(PARAM.inp.nspin == 2) + if (PARAM.inp.nspin == 2) { const_cast*>(dm_p)->switch_dmr(1); } const hamilt::HContainer* dmr = dm_p->get_DMR_pointer(1); tmp_nonlocal.cal_force_stress(isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta); - if(PARAM.inp.nspin == 2) + if (PARAM.inp.nspin == 2) { const_cast*>(dm_p)->switch_dmr(0); } } - else if(PARAM.inp.nspin == 4) + else if (PARAM.inp.nspin == 4) { hamilt::NonlocalNew, std::complex>> tmp_nonlocal( - nullptr, - kv.kvec_d, - nullptr, - &ucell, - orb.cutoffs(), - &GlobalC::GridD, - two_center_bundle.overlap_orb_beta.get() - ); + nullptr, + kv.kvec_d, + nullptr, + &ucell, + orb.cutoffs(), + &gd, + two_center_bundle.overlap_orb_beta.get()); // calculate temporary complex DMR for nonlocal force&stress // In fact, only SOC part need the imaginary part of DMR for correct force&stress @@ -225,7 +224,6 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, dm_p->cal_DMR_full(&tmp_dmr); tmp_nonlocal.cal_force_stress(isforce, isstress, &tmp_dmr, fvnl_dbeta, svnl_dbeta); } - //! forces and stress from vdw // Peize Lin add 2014-04-04, update 2021-03-09 @@ -300,13 +298,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } if (PARAM.inp.dft_plus_u == 2) { - GlobalC::dftu.force_stress(ucell, - pelec, - pv, - fsr, // mohan 2024-06-16 - force_dftu, - stress_dftu, - kv); + GlobalC::dftu.force_stress(ucell, gd, pelec, pv, fsr, force_dftu, stress_dftu, kv); } else { @@ -314,7 +306,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, kv.kvec_d, nullptr, // HR are not used for force&stress ucell, - &GlobalC::GridD, + &gd, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), &GlobalC::dftu); @@ -326,7 +318,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // atomic force and stress for DeltaSpin ModuleBase::matrix force_dspin; ModuleBase::matrix stress_dspin; - if(PARAM.inp.sc_mag_switch) + if (PARAM.inp.sc_mag_switch) { if (isforce) { @@ -337,24 +329,22 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, stress_dspin.create(3, 3); } - hamilt::DeltaSpin> tmp_dspin( - nullptr, - kv.kvec_d, - nullptr, - ucell, - &GlobalC::GridD, - two_center_bundle.overlap_orb_onsite.get(), - orb.cutoffs() - ); + hamilt::DeltaSpin> tmp_dspin(nullptr, + kv.kvec_d, + nullptr, + ucell, + &gd, + two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs()); const auto* dm_p = dynamic_cast>*>(pelec)->get_DM(); - if(PARAM.inp.nspin == 2) + if (PARAM.inp.nspin == 2) { const_cast, double>*>(dm_p)->switch_dmr(2); } const hamilt::HContainer* dmr = dm_p->get_DMR_pointer(1); tmp_dspin.cal_force_stress(isforce, isstress, dmr, force_dspin, stress_dspin); - if(PARAM.inp.nspin == 2) + if (PARAM.inp.nspin == 2) { const_cast, double>*>(dm_p)->switch_dmr(0); } @@ -388,12 +378,12 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { if (GlobalC::exx_info.info_ri.real_number) { - exx_lri_double.cal_exx_stress(ucell.omega,ucell.lat0); + 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(ucell.omega,ucell.lat0); + exx_lri_complex.cal_exx_stress(ucell.omega, ucell.lat0); stress_exx = GlobalC::exx_info.info_global.hybrid_alpha * exx_lri_complex.stress_exx; } } @@ -495,7 +485,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // pengfei 2016-12-20 if (ModuleSymmetry::Symmetry::symm_flag == 1) { - this->forceSymmetry(ucell,fcs, symm); + this->forceSymmetry(ucell, fcs, symm); } #ifdef __DEEPKS @@ -503,17 +493,15 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.deepks_out_labels) // not parallelized yet { const std::string file_ftot = PARAM.globalv.global_out_dir + "deepks_ftot.npy"; - LCAO_deepks_io::save_npy_f(fcs, - file_ftot, - ucell.nat, + LCAO_deepks_io::save_npy_f(fcs, file_ftot, ucell.nat, GlobalV::MY_RANK); // Ty/Bohr, F_tot if (PARAM.inp.deepks_scf) { const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - LCAO_deepks_io::save_npy_f(fcs - GlobalC::ld.F_delta, - file_fbase, - ucell.nat, + LCAO_deepks_io::save_npy_f(fcs - GlobalC::ld.F_delta, + file_fbase, + ucell.nat, GlobalV::MY_RANK); // Ry/Bohr, F_base if (!PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now @@ -522,7 +510,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { const std::vector>& dm_gamma = dynamic_cast*>(pelec)->get_DM()->get_DMK_vector(); - GlobalC::ld.cal_gdmx(dm_gamma[0], ucell, orb, GlobalC::GridD, isstress); + GlobalC::ld.cal_gdmx(dm_gamma[0], ucell, orb, gd, isstress); } else { @@ -531,13 +519,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, ->get_DM() ->get_DMK_vector(); - GlobalC::ld.cal_gdmx_k(dm_k, - ucell, - orb, - GlobalC::GridD, - kv.get_nks(), - kv.kvec_d, - isstress); + GlobalC::ld.cal_gdmx_k(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress); } if (PARAM.inp.deepks_out_unittest) { @@ -892,6 +874,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, const bool isforce, const bool isstress, const UnitCell& ucell, + Grid_Driver& gd, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -918,6 +901,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, isstress, fsr, // mohan add 2024-06-15 ucell, + gd, psi, pelec, foverlap, @@ -943,6 +927,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn const bool isforce, const bool isstress, const UnitCell& ucell, + Grid_Driver& gd, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi>* psi, @@ -968,6 +953,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn isstress, fsr, // mohan add 2024-06-16 ucell, + gd, psi, pelec, foverlap, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index 927d08790e..5593fd0afb 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -34,6 +34,7 @@ class Force_Stress_LCAO const bool istestf, const bool istests, const UnitCell& ucell, + Grid_Driver& gd, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -81,6 +82,7 @@ class Force_Stress_LCAO const bool isforce, const bool isstress, const UnitCell& ucell, + Grid_Driver& gd, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi* psi, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index c7aad83123..5d312ad0ce 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -17,6 +17,7 @@ template <> void Force_LCAO::allocate(const UnitCell& ucell, + Grid_Driver& gd, const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, @@ -80,7 +81,7 @@ void Force_LCAO::allocate(const UnitCell& ucell, orb, pv, two_center_bundle, - &GlobalC::GridD, + &gd, nullptr); // calculate dT in LCAP @@ -104,7 +105,7 @@ void Force_LCAO::allocate(const UnitCell& ucell, orb, pv, two_center_bundle, - &GlobalC::GridD, + &gd, nullptr); // calculate asynchronous S matrix to output for Hefei-NAMD @@ -175,6 +176,7 @@ void Force_LCAO::ftable(const bool isforce, const bool isstress, ForceStressArrays& fsr, // mohan add 2024-06-16 const UnitCell& ucell, + Grid_Driver& gd, const psi::Psi* psi, const elecstate::ElecState* pelec, ModuleBase::matrix& foverlap, @@ -206,7 +208,7 @@ void Force_LCAO::ftable(const bool isforce, // allocate DSloc_x, DSloc_y, DSloc_z // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z - this->allocate(ucell,pv, fsr, two_center_bundle, orb); + this->allocate(ucell, gd, pv, fsr, two_center_bundle, orb); const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z }; const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 }; @@ -229,25 +231,24 @@ void Force_LCAO::ftable(const bool isforce, const std::vector>& dm_gamma = dm->get_DMK_vector(); // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm - //GlobalC::ld.cal_projected_DM(dm, ucell, orb, GlobalC::GridD); + // GlobalC::ld.cal_projected_DM(dm, ucell, orb, gd); GlobalC::ld.cal_descriptor(ucell.nat); GlobalC::ld.cal_gedm(ucell.nat); - DeePKS_domain::cal_f_delta_gamma( - dm_gamma, - ucell, - orb, - GlobalC::GridD, - *this->ParaV, - GlobalC::ld.lmaxd, - GlobalC::ld.nlm_save, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, - GlobalC::ld.F_delta, - isstress, - svnl_dalpha); + DeePKS_domain::cal_f_delta_gamma(dm_gamma, + ucell, + orb, + gd, + *this->ParaV, + GlobalC::ld.lmaxd, + GlobalC::ld.nlm_save, + GlobalC::ld.gedm, + GlobalC::ld.inl_index, + GlobalC::ld.F_delta, + isstress, + svnl_dalpha); #ifdef __MPI Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index f00b089a75..6b88139fc4 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -27,6 +27,7 @@ template <> void Force_LCAO>::allocate(const UnitCell& ucell, + Grid_Driver& gd, const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, @@ -98,7 +99,7 @@ void Force_LCAO>::allocate(const UnitCell& ucell, orb, pv, two_center_bundle, - &GlobalC::GridD, + &gd, nullptr); // delete lm.SlocR //----------------------------------------- @@ -129,7 +130,7 @@ void Force_LCAO>::allocate(const UnitCell& ucell, orb, pv, two_center_bundle, - &GlobalC::GridD, + &gd, nullptr); // delete lm.Hloc_fixedR // calculate asynchronous S matrix to output for Hefei-NAMD @@ -148,7 +149,7 @@ void Force_LCAO>::allocate(const UnitCell& ucell, orb, pv, two_center_bundle, - &(GlobalC::GridD), + &(gd), nullptr, // delete lm.SlocR PARAM.inp.cal_syns, PARAM.inp.dmax); @@ -271,6 +272,7 @@ void Force_LCAO>::ftable(const bool isforce, const bool isstress, ForceStressArrays& fsr, // mohan add 2024-06-15 const UnitCell& ucell, + Grid_Driver& gd, const psi::Psi>* psi, const elecstate::ElecState* pelec, ModuleBase::matrix& foverlap, @@ -298,6 +300,7 @@ void Force_LCAO>::ftable(const bool isforce, = dynamic_cast>*>(pelec)->get_DM(); this->allocate(ucell, + gd, pv, fsr, // mohan add 2024-06-16 two_center_bundle, @@ -327,27 +330,26 @@ void Force_LCAO>::ftable(const bool isforce, const std::vector>>& dm_k = dm->get_DMK_vector(); // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm - //GlobalC::ld.cal_projected_DM_k(dm, ucell, orb, GlobalC::GridD); + // GlobalC::ld.cal_projected_DM_k(dm, ucell, orb, gd); GlobalC::ld.cal_descriptor(ucell.nat); GlobalC::ld.cal_gedm(ucell.nat); - DeePKS_domain::cal_f_delta_k( - dm_k, - ucell, - orb, - GlobalC::GridD, - pv, - GlobalC::ld.lmaxd, - kv->get_nks(), - kv->kvec_d, - GlobalC::ld.nlm_save_k, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, - GlobalC::ld.F_delta, - isstress, - svnl_dalpha); + DeePKS_domain::cal_f_delta_k(dm_k, + ucell, + orb, + gd, + pv, + GlobalC::ld.lmaxd, + kv->get_nks(), + kv->kvec_d, + GlobalC::ld.nlm_save_k, + GlobalC::ld.gedm, + GlobalC::ld.inl_index, + GlobalC::ld.F_delta, + isstress, + svnl_dalpha); #ifdef __MPI Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 572c26e4d6..b08da06af8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -192,7 +192,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, // no explicit call yet &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, PARAM.inp.nspin); this->getOperator()->add(veff); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp index 95df523db2..0e8c501f9d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp @@ -30,6 +30,7 @@ DeePKS>::DeePKS(HS_Matrix_K* hsk_in, intor_orb_alpha_(intor_orb_alpha), ptr_orb_(ptr_orb), nks(nks_in) { this->cal_type = calculation_type::lcao_deepks; + this->gd = GridD_in; #ifdef __DEEPKS this->initialize_HR(GridD_in); #endif @@ -158,7 +159,7 @@ void DeePKS>::contributeHR() { ModuleBase::timer::tick("DeePKS", "contributeHR"); const Parallel_Orbitals* pv = this->hsk->get_pv(); - GlobalC::ld.cal_projected_DM(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD); + GlobalC::ld.cal_projected_DM(this->DM, *this->ucell, *ptr_orb_, *(this->gd)); GlobalC::ld.cal_descriptor(this->ucell->nat); GlobalC::ld.cal_gedm(this->ucell->nat); // recalculate the H_V_delta @@ -186,7 +187,7 @@ void DeePKS, double>>::contributeHR() { ModuleBase::timer::tick("DeePKS", "contributeHR"); - GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD); + GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, *(this->gd)); GlobalC::ld.cal_descriptor(this->ucell->nat); // calculate dE/dD GlobalC::ld.cal_gedm(this->ucell->nat); @@ -219,7 +220,7 @@ void DeePKS, std::complex>>::contribut { ModuleBase::timer::tick("DeePKS", "contributeHR"); - GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, GlobalC::GridD); + GlobalC::ld.cal_projected_DM_k(this->DM, *this->ucell, *ptr_orb_, *(this->gd)); GlobalC::ld.cal_descriptor(this->ucell->nat); // calculate dE/dD GlobalC::ld.cal_gedm(this->ucell->nat); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h index f95b18a607..32f7d1cddc 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h @@ -59,6 +59,8 @@ class DeePKS> : public OperatorLCAO const UnitCell* ucell = nullptr; + Grid_Driver* gd = nullptr; + HContainer* H_V_delta = nullptr; // the following variable is introduced temporarily during LCAO refactoring diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp index 9710ec4b26..31be5f4582 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp @@ -34,16 +34,16 @@ void sparse_format::cal_dH(const UnitCell& ucell, const bool cal_deri = true; const bool cal_stress = false; LCAO_domain::build_ST_new(fsr_dh, - 'T', - cal_deri, - cal_stress, - ucell, - orb, - pv, - two_center_bundle, - &GlobalC::GridD, - nullptr, - false); // delete unused parameter lm.Hloc_fixedR + 'T', + cal_deri, + cal_stress, + ucell, + orb, + pv, + two_center_bundle, + &grid, + nullptr, + false); // delete unused parameter lm.Hloc_fixedR LCAO_domain::build_Nonlocal_mu_new(pv, fsr_dh, @@ -52,7 +52,7 @@ void sparse_format::cal_dH(const UnitCell& ucell, ucell, orb, *(two_center_bundle.overlap_orb_beta), - &GlobalC::GridD); + &grid); sparse_format::cal_dSTN_R(ucell,pv, HS_Arrays, fsr_dh, grid, orb.cutoffs(), current_spin, sparse_thr); @@ -60,8 +60,7 @@ void sparse_format::cal_dH(const UnitCell& ucell, delete[] fsr_dh.DHloc_fixedR_y; delete[] fsr_dh.DHloc_fixedR_z; - gint_k - .cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, ucell, GlobalC::GridD); + gint_k.cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, ucell, grid); return; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp index e4fbba3624..1435c91e3e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_st.cpp @@ -79,7 +79,7 @@ void sparse_format::cal_TR(const UnitCell& ucell, orb, pv, two_center_bundle, - &(GlobalC::GridD), + &(grid), HS_Arrays.Hloc_fixedR.data()); sparse_format::set_R_range(HS_Arrays.all_R_coor, grid); diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index 68aae44516..7e8909c96f 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -180,16 +180,15 @@ class DFTU // dim = 1-3 : dS, for force // dim = 4-6 : dS * dR, for stress - void folding_matrix_k( - const UnitCell &ucell, - ForceStressArrays &fsr, - const Parallel_Orbitals &pv, - const int ik, - const int dim1, - const int dim2, - std::complex* mat_k, - const std::vector> &kvec_d); - + void folding_matrix_k(const UnitCell& ucell, + Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const int dim1, + const int dim2, + std::complex* mat_k, + const std::vector>& kvec_d); /** * @brief new function of folding_S_matrix @@ -202,9 +201,9 @@ class DFTU // In dftu_force.cpp // For calculating force and stress fomr DFT+U //============================================================= - public: - + public: void force_stress(const UnitCell& ucell, + Grid_Driver& gd, const elecstate::ElecState* pelec, const Parallel_Orbitals& pv, ForceStressArrays& fsr, @@ -212,43 +211,42 @@ class DFTU ModuleBase::matrix& stress_dftu, const K_Vectors& kv); - private: - - void cal_force_k(const UnitCell &ucell, - ForceStressArrays &fsr, - const Parallel_Orbitals &pv, + private: + void cal_force_k(const UnitCell& ucell, + Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, const int ik, const std::complex* rho_VU, ModuleBase::matrix& force_dftu, const std::vector>& kvec_d); - void cal_stress_k( - const UnitCell &ucell, - ForceStressArrays &fsr, - const Parallel_Orbitals &pv, - const int ik, - const std::complex* rho_VU, - ModuleBase::matrix& stress_dftu, - const std::vector>& kvec_d); - - void cal_force_gamma(const UnitCell &ucell, - const double* rho_VU, - const Parallel_Orbitals &pv, - double* dsloc_x, - double* dsloc_y, - double* dsloc_z, - ModuleBase::matrix& force_dftu); - - void cal_stress_gamma( - const UnitCell &ucell, - const Parallel_Orbitals &pv, - Grid_Driver* gd, - double* dsloc_x, - double* dsloc_y, - double* dsloc_z, - double* dh_r, - const double* rho_VU, - ModuleBase::matrix& stress_dftu); + void cal_stress_k(const UnitCell& ucell, + Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const std::complex* rho_VU, + ModuleBase::matrix& stress_dftu, + const std::vector>& kvec_d); + + void cal_force_gamma(const UnitCell& ucell, + const double* rho_VU, + const Parallel_Orbitals& pv, + double* dsloc_x, + double* dsloc_y, + double* dsloc_z, + ModuleBase::matrix& force_dftu); + + void cal_stress_gamma(const UnitCell& ucell, + const Parallel_Orbitals& pv, + Grid_Driver* gd, + double* dsloc_x, + double* dsloc_y, + double* dsloc_z, + double* dh_r, + const double* rho_VU, + ModuleBase::matrix& stress_dftu); #endif //============================================================= diff --git a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp index aac95478d7..12904bbd72 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp @@ -126,25 +126,33 @@ void DFTU::fold_dSR_gamma( return; } -void DFTU::folding_matrix_k( - const UnitCell &ucell, - ForceStressArrays &fsr, - const Parallel_Orbitals &pv, - const int ik, - const int dim1, - const int dim2, - std::complex* mat_k, - const std::vector> &kvec_d) +void DFTU::folding_matrix_k(const UnitCell& ucell, + Grid_Driver& gd, + ForceStressArrays& fsr, + const Parallel_Orbitals& pv, + const int ik, + const int dim1, + const int dim2, + std::complex* mat_k, + const std::vector>& kvec_d) { ModuleBase::TITLE("DFTU", "folding_matrix_k"); ModuleBase::timer::tick("DFTU", "folding_matrix_k"); ModuleBase::GlobalFunc::ZEROS(mat_k, pv.nloc); double* mat_ptr; - if (dim1 == 1 || dim1 == 4) { mat_ptr = fsr.DSloc_Rx; - } else if (dim1 == 2 || dim1 == 5) { mat_ptr = fsr.DSloc_Ry; - } else if (dim1 == 3 || dim1 == 6) { mat_ptr = fsr.DSloc_Rz; -} + if (dim1 == 1 || dim1 == 4) + { + mat_ptr = fsr.DSloc_Rx; + } + else if (dim1 == 2 || dim1 == 5) + { + mat_ptr = fsr.DSloc_Ry; + } + else if (dim1 == 3 || dim1 == 6) + { + mat_ptr = fsr.DSloc_Rz; + } int nnr = 0; ModuleBase::Vector3 dtau; @@ -161,18 +169,18 @@ void DFTU::folding_matrix_k( for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; - GlobalC::GridD.Find_atom(ucell, tau1, T1, I1); + gd.Find_atom(ucell, tau1, T1, I1); Atom* atom1 = &ucell.atoms[T1]; const int start1 = ucell.itiaiw2iwt(T1, I1, 0); // (2) search among all adjacent atoms. - for (int ad = 0; ad < GlobalC::GridD.getAdjacentNum() + 1; ++ad) + for (int ad = 0; ad < gd.getAdjacentNum() + 1; ++ad) { - const int T2 = GlobalC::GridD.getType(ad); - const int I2 = GlobalC::GridD.getNatom(ad); + const int T2 = gd.getType(ad); + const int I2 = gd.getNatom(ad); Atom* atom2 = &ucell.atoms[T2]; - tau2 = GlobalC::GridD.getAdjacentTau(ad); + tau2 = gd.getAdjacentTau(ad); dtau = tau2 - tau1; double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff_[T1] + orb_cutoff_[T2]; @@ -185,12 +193,12 @@ void DFTU::folding_matrix_k( } else if (distance >= rcut) { - for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum() + 1; ++ad0) + for (int ad0 = 0; ad0 < gd.getAdjacentNum() + 1; ++ad0) { - const int T0 = GlobalC::GridD.getType(ad0); - const int I0 = GlobalC::GridD.getNatom(ad0); + const int T0 = gd.getType(ad0); + const int I0 = gd.getNatom(ad0); - tau0 = GlobalC::GridD.getAdjacentTau(ad0); + tau0 = gd.getAdjacentTau(ad0); dtau1 = tau0 - tau1; dtau2 = tau0 - tau2; @@ -216,9 +224,7 @@ void DFTU::folding_matrix_k( // exp(k dot dR) // dR is the index of box in Crystal coordinates //------------------------------------------------ - ModuleBase::Vector3 dR(GlobalC::GridD.getBox(ad).x, - GlobalC::GridD.getBox(ad).y, - GlobalC::GridD.getBox(ad).z); + ModuleBase::Vector3 dR(gd.getBox(ad).x, gd.getBox(ad).y, gd.getBox(ad).z); const double arg = (kvec_d[ik] * dR) * ModuleBase::TWO_PI; const std::complex kphase = std::complex(cos(arg), sin(arg)); diff --git a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp index 3ab4ef2496..1903a33899 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp @@ -73,6 +73,7 @@ namespace ModuleDFTU { void DFTU::force_stress(const UnitCell& ucell, + Grid_Driver& gd, const elecstate::ElecState* pelec, const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-16 @@ -149,7 +150,7 @@ void DFTU::force_stress(const UnitCell& ucell, { this->cal_stress_gamma(ucell, pv, - &GlobalC::GridD, + &gd, fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z, @@ -208,11 +209,11 @@ void DFTU::force_stress(const UnitCell& ucell, if (PARAM.inp.cal_force) { - cal_force_k(ucell,fsr, pv, ik, &rho_VU[0], force_dftu, kv.kvec_d); + cal_force_k(ucell, gd, fsr, pv, ik, &rho_VU[0], force_dftu, kv.kvec_d); } if (PARAM.inp.cal_stress) { - cal_stress_k(ucell,fsr, pv, ik, &rho_VU[0], stress_dftu, kv.kvec_d); + cal_stress_k(ucell, gd, fsr, pv, ik, &rho_VU[0], stress_dftu, kv.kvec_d); } } // ik } @@ -249,6 +250,7 @@ void DFTU::force_stress(const UnitCell& ucell, } void DFTU::cal_force_k(const UnitCell& ucell, + Grid_Driver& gd, ForceStressArrays& fsr, const Parallel_Orbitals& pv, const int ik, @@ -270,7 +272,7 @@ void DFTU::cal_force_k(const UnitCell& ucell, for (int dim = 0; dim < 3; dim++) { - this->folding_matrix_k(ucell,fsr, pv, ik, dim + 1, 0, &dSm_k[0], kvec_d); + this->folding_matrix_k(ucell, gd, fsr, pv, ik, dim + 1, 0, &dSm_k[0], kvec_d); #ifdef __MPI pzgemm_(&transN, @@ -378,6 +380,7 @@ void DFTU::cal_force_k(const UnitCell& ucell, } void DFTU::cal_stress_k(const UnitCell& ucell, + Grid_Driver& gd, ForceStressArrays& fsr, const Parallel_Orbitals& pv, const int ik, @@ -403,14 +406,7 @@ void DFTU::cal_stress_k(const UnitCell& ucell, { for (int dim2 = dim1; dim2 < 3; dim2++) { - this->folding_matrix_k(ucell, - fsr, // mohan add 2024-06-16 - pv, - ik, - dim1 + 4, - dim2, - &dSR_k[0], - kvec_d); + this->folding_matrix_k(ucell, gd, fsr, pv, ik, dim1 + 4, dim2, &dSR_k[0], kvec_d); #ifdef __MPI pzgemm_(&transN, diff --git a/source/module_io/berryphase.cpp b/source/module_io/berryphase.cpp index a8b5eef5b3..d2ce73230d 100644 --- a/source/module_io/berryphase.cpp +++ b/source/module_io/berryphase.cpp @@ -40,11 +40,15 @@ void berryphase::get_occupation_bands() } #ifdef __LCAO -void berryphase::lcao_init(const UnitCell& ucell, const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb) +void berryphase::lcao_init(const UnitCell& ucell, + Grid_Driver& gd, + const K_Vectors& kv, + const Grid_Technique& grid_tech, + const LCAO_Orbitals& orb) { ModuleBase::TITLE("berryphase", "lcao_init"); lcao_method.init(ucell,grid_tech, kv.get_nkstot(), orb); - lcao_method.cal_R_number(ucell); + lcao_method.cal_R_number(ucell, gd); lcao_method.cal_orb_overlap(ucell); return; } diff --git a/source/module_io/berryphase.h b/source/module_io/berryphase.h index 7a3707c711..c0fbe215be 100644 --- a/source/module_io/berryphase.h +++ b/source/module_io/berryphase.h @@ -37,7 +37,11 @@ class berryphase void get_occupation_bands(); #ifdef __LCAO - void lcao_init(const UnitCell& ucell, const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb); + void lcao_init(const UnitCell& ucell, + Grid_Driver& gd, + const K_Vectors& kv, + const Grid_Technique& grid_tech, + const LCAO_Orbitals& orb); #endif void set_kpoints(const K_Vectors& kv, const int direction); diff --git a/source/module_io/cal_r_overlap_R.cpp b/source/module_io/cal_r_overlap_R.cpp index b43c181576..2758ad131d 100644 --- a/source/module_io/cal_r_overlap_R.cpp +++ b/source/module_io/cal_r_overlap_R.cpp @@ -239,20 +239,20 @@ void cal_r_overlap_R::init(const UnitCell& ucell,const Parallel_Orbitals& pv, co return; } -void cal_r_overlap_R::out_rR(const UnitCell& ucell, const int& istep) +void cal_r_overlap_R::out_rR(const UnitCell& ucell, Grid_Driver& gd, const int& istep) { ModuleBase::TITLE("cal_r_overlap_R", "out_rR"); ModuleBase::timer::tick("cal_r_overlap_R", "out_rR"); int step = istep; // set R coor range - int R_minX = int(-GlobalC::GridD.getTrueCellX()); - int R_minY = int(-GlobalC::GridD.getTrueCellY()); - int R_minZ = int(-GlobalC::GridD.getTrueCellZ()); + int R_minX = int(-gd.getTrueCellX()); + int R_minY = int(-gd.getTrueCellY()); + int R_minZ = int(-gd.getTrueCellZ()); - int R_x = GlobalC::GridD.getCellX(); - int R_y = GlobalC::GridD.getCellY(); - int R_z = GlobalC::GridD.getCellZ(); + int R_x = gd.getCellX(); + int R_y = gd.getCellY(); + int R_z = gd.getCellZ(); std::set> all_R_coor; for (int ix = 0; ix < R_x; ix++) diff --git a/source/module_io/cal_r_overlap_R.h b/source/module_io/cal_r_overlap_R.h index ae23718148..69a4d2c7fb 100644 --- a/source/module_io/cal_r_overlap_R.h +++ b/source/module_io/cal_r_overlap_R.h @@ -9,11 +9,13 @@ #include "module_basis/module_ao/ORB_gaunt_table.h" #include "module_basis/module_ao/ORB_read.h" #include "module_basis/module_ao/parallel_orbitals.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb11.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb21.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb.h" #include "single_R_io.h" -#include "module_cell/unitcell.h" + #include #include #include @@ -31,7 +33,7 @@ class cal_r_overlap_R bool binary = false; void init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); - void out_rR(const UnitCell& ucell, const int& istep); + void out_rR(const UnitCell& ucell, Grid_Driver& gd, const int& istep); void out_rR_other(const UnitCell& ucell, const int& istep, const std::set>& output_R_coor); private: diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index d3a6f6e1b1..f12fd69d3e 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -80,7 +80,7 @@ void output_mat_sparse(const bool& out_mat_hsR, } else { - r_matrix.out_rR(ucell,istep); + r_matrix.out_rR(ucell, grid, istep); } } diff --git a/source/module_io/output_mulliken.h b/source/module_io/output_mulliken.h index 2d78d2fa52..1bcea3a21f 100644 --- a/source/module_io/output_mulliken.h +++ b/source/module_io/output_mulliken.h @@ -93,6 +93,7 @@ void cal_mag(Parallel_Orbitals* pv, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, UnitCell& ucell, + Grid_Driver& gd, const int istep, const bool print) { @@ -134,15 +135,14 @@ void cal_mag(Parallel_Orbitals* pv, auto atomLabels = ucell.get_atomLabels(); if(PARAM.inp.nspin == 2) { - auto sc_lambda = new hamilt::DeltaSpin>( - nullptr, - kv.kvec_d, - nullptr, - ucell, - &GlobalC::GridD, - two_center_bundle.overlap_orb_onsite.get(), - orb.cutoffs() - ); + auto sc_lambda + = new hamilt::DeltaSpin>(nullptr, + kv.kvec_d, + nullptr, + ucell, + &gd, + two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs()); dynamic_cast*>(pelec)->get_DM()->switch_dmr(2); moments = sc_lambda->cal_moment(dmr, constrain); dynamic_cast*>(pelec)->get_DM()->switch_dmr(0); @@ -162,14 +162,13 @@ void cal_mag(Parallel_Orbitals* pv, else if(PARAM.inp.nspin == 4) { auto sc_lambda = new hamilt::DeltaSpin, std::complex>>( - nullptr, - kv.kvec_d, - nullptr, - ucell, - &GlobalC::GridD, - two_center_bundle.overlap_orb_onsite.get(), - orb.cutoffs() - ); + nullptr, + kv.kvec_d, + nullptr, + ucell, + &gd, + two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs()); moments = sc_lambda->cal_moment(dmr, constrain); delete sc_lambda; //const std::vector title = {"Total Magnetism (uB)", "", "", ""}; diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index 96272a70e7..9534107da2 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -118,6 +118,7 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix, double> } void ModuleIO::write_current(const UnitCell& ucell, + Grid_Driver& gd, const int istep, const psi::Psi>* psi, const elecstate::ElecState* pelec, @@ -134,7 +135,7 @@ void ModuleIO::write_current(const UnitCell& ucell, std::vector>*> current_term = {nullptr, nullptr, nullptr}; if (!TD_Velocity::tddft_velocity) { - cal_current = new TD_current(&ucell, &GlobalC::GridD, pv, orb, intor); + cal_current = new TD_current(&ucell, &gd, pv, orb, intor); cal_current->calculate_vcomm_r(); cal_current->calculate_grad_term(); for (int dir = 0; dir < 3; dir++) diff --git a/source/module_io/td_current_io.h b/source/module_io/td_current_io.h index 40aed95214..709a85c4d6 100644 --- a/source/module_io/td_current_io.h +++ b/source/module_io/td_current_io.h @@ -11,6 +11,7 @@ namespace ModuleIO #ifdef __LCAO /// @brief func to output current, only used in tddft void write_current(const UnitCell& ucell, + Grid_Driver& gd, const int istep, const psi::Psi>* psi, const elecstate::ElecState* pelec, diff --git a/source/module_io/to_wannier90_lcao.cpp b/source/module_io/to_wannier90_lcao.cpp index 72651d6991..3ecad5ea36 100644 --- a/source/module_io/to_wannier90_lcao.cpp +++ b/source/module_io/to_wannier90_lcao.cpp @@ -39,6 +39,7 @@ toWannier90_LCAO::~toWannier90_LCAO() } void toWannier90_LCAO::calculate(const UnitCell& ucell, + Grid_Driver& gd, const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi>& psi, @@ -114,7 +115,7 @@ void toWannier90_LCAO::calculate(const UnitCell& ucell, initialize_orb_table(ucell); produce_basis_orb(); - set_R_coor(ucell); + set_R_coor(ucell, gd); count_delta_k(ucell,kv); } @@ -138,7 +139,7 @@ void toWannier90_LCAO::calculate(const UnitCell& ucell, return exp_idkr; }; - FR[i].set_parameters(fr_ptr[i], &ucell, &orb_, &GlobalC::GridD, ParaV, 140, 110); + FR[i].set_parameters(fr_ptr[i], &ucell, &orb_, &gd, ParaV, 140, 110); FR[i].calculate_FR(); } @@ -303,15 +304,15 @@ void toWannier90_LCAO::initialize_orb_table(const UnitCell& ucell) #endif } -void toWannier90_LCAO::set_R_coor(const UnitCell& ucell) +void toWannier90_LCAO::set_R_coor(const UnitCell& ucell, const Grid_Driver& gd) { - int R_minX = int(-GlobalC::GridD.getTrueCellX()); - int R_minY = int(-GlobalC::GridD.getTrueCellY()); - int R_minZ = int(-GlobalC::GridD.getTrueCellZ()); + int R_minX = int(-gd.getTrueCellX()); + int R_minY = int(-gd.getTrueCellY()); + int R_minZ = int(-gd.getTrueCellZ()); - int R_x = GlobalC::GridD.getCellX(); - int R_y = GlobalC::GridD.getCellY(); - int R_z = GlobalC::GridD.getCellZ(); + int R_x = gd.getCellX(); + int R_y = gd.getCellY(); + int R_z = gd.getCellZ(); int R_num = R_x * R_y * R_z; R_coor_car.resize(R_num); diff --git a/source/module_io/to_wannier90_lcao.h b/source/module_io/to_wannier90_lcao.h index 34d0e64dbb..9a28e71d56 100644 --- a/source/module_io/to_wannier90_lcao.h +++ b/source/module_io/to_wannier90_lcao.h @@ -80,12 +80,14 @@ class toWannier90_LCAO : public toWannier90 ~toWannier90_LCAO(); void calculate(const UnitCell& ucell, + Grid_Driver& gd, const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi>& psi, const Parallel_Orbitals* pv); void calculate(const UnitCell& ucell, + Grid_Driver& gd, const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi& psi, @@ -131,7 +133,7 @@ class toWannier90_LCAO : public toWannier90 void initialize_orb_table(const UnitCell& ucell); void produce_basis_orb(); - void set_R_coor(const UnitCell& ucell); + void set_R_coor(const UnitCell& ucell, const Grid_Driver& gd); void count_delta_k(const UnitCell& ucell, const K_Vectors& kv); std::vector delta_k_all; diff --git a/source/module_io/unk_overlap_lcao.cpp b/source/module_io/unk_overlap_lcao.cpp index b784f54a1c..78dbba044e 100644 --- a/source/module_io/unk_overlap_lcao.cpp +++ b/source/module_io/unk_overlap_lcao.cpp @@ -360,7 +360,7 @@ int unkOverlap_lcao::iw2im(const UnitCell& ucell, int iw) } // search for the nearest neighbor atoms -void unkOverlap_lcao::cal_R_number(const UnitCell& ucell) +void unkOverlap_lcao::cal_R_number(const UnitCell& ucell, Grid_Driver& gd) { // The number of overlaps between atomic orbitals 1 and atomic orbitals 2, // or the number of R, is empty when there is no overlap @@ -377,18 +377,18 @@ void unkOverlap_lcao::cal_R_number(const UnitCell& ucell) for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; - GlobalC::GridD.Find_atom(ucell, tau1, T1, I1); + gd.Find_atom(ucell, tau1, T1, I1); - for (int ad = 0; ad < GlobalC::GridD.getAdjacentNum() + 1; ++ad) + for (int ad = 0; ad < gd.getAdjacentNum() + 1; ++ad) { - const int T2 = GlobalC::GridD.getType(ad); - const int I2 = GlobalC::GridD.getNatom(ad); + const int T2 = gd.getType(ad); + const int I2 = gd.getNatom(ad); Atom* atom2 = &ucell.atoms[T2]; - const double R_direct_x = (double)GlobalC::GridD.getBox(ad).x; - const double R_direct_y = (double)GlobalC::GridD.getBox(ad).y; - const double R_direct_z = (double)GlobalC::GridD.getBox(ad).z; + const double R_direct_x = (double)gd.getBox(ad).x; + const double R_direct_y = (double)gd.getBox(ad).y; + const double R_direct_z = (double)gd.getBox(ad).z; - tau2 = GlobalC::GridD.getAdjacentTau(ad); + tau2 = gd.getAdjacentTau(ad); dtau = tau2 - tau1; double distance = dtau.norm() * ucell.lat0; double rcut = rcut_orb_[T1] + rcut_orb_[T2]; diff --git a/source/module_io/unk_overlap_lcao.h b/source/module_io/unk_overlap_lcao.h index 0ba664271e..2554da1142 100644 --- a/source/module_io/unk_overlap_lcao.h +++ b/source/module_io/unk_overlap_lcao.h @@ -54,7 +54,7 @@ class unkOverlap_lcao int iw2iL(const UnitCell& ucell, int iw); int iw2iN(const UnitCell& ucell, int iw); int iw2im(const UnitCell& ucell, int iw); - void cal_R_number(const UnitCell& ucell); + void cal_R_number(const UnitCell& ucell, Grid_Driver& gd); void cal_orb_overlap(const UnitCell& ucell); void prepare_midmatrix_pblas(const UnitCell& ucell, const int ik_L, diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index 502ef01120..d47e962d72 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -326,11 +326,11 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu ucell.infoNL.get_rcutmax_Beta(), PARAM.globalv.gamma_only_local); atom_arrange::search(PARAM.inp.search_pbc, - GlobalV::ofs_running, - GlobalC::GridD, - this->ucell, - search_radius, - PARAM.inp.test_atom_input); + GlobalV::ofs_running, + this->gd, + this->ucell, + search_radius, + PARAM.inp.test_atom_input); this->set_gint(); this->gint_->gridt = &this->gt_; @@ -343,28 +343,28 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu Gint_Tools::init_orb(dr_uniform, rcuts, ucell, orb, psi_u, dpsi_u, d2psi_u); this->gt_.set_pbc_grid(this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pw_big->bx, - this->pw_big->by, - this->pw_big->bz, - this->pw_big->nbx, - this->pw_big->nby, - this->pw_big->nbz, - this->pw_big->nbxx, - this->pw_big->nbzp_start, - this->pw_big->nbzp, - this->pw_rho->ny, - this->pw_rho->nplane, - this->pw_rho->startz_current, - ucell, - GlobalC::GridD, - dr_uniform, - rcuts, - psi_u, - dpsi_u, - d2psi_u, - PARAM.inp.nstream); + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bx, + this->pw_big->by, + this->pw_big->bz, + this->pw_big->nbx, + this->pw_big->nby, + this->pw_big->nbz, + this->pw_big->nbxx, + this->pw_big->nbzp_start, + this->pw_big->nbzp, + this->pw_rho->ny, + this->pw_rho->nplane, + this->pw_rho->startz_current, + ucell, + this->gd, + dr_uniform, + rcuts, + psi_u, + dpsi_u, + d2psi_u, + PARAM.inp.nstream); psi_u.clear(); psi_u.shrink_to_fit(); dpsi_u.clear(); @@ -388,7 +388,7 @@ LR::ESolver_LR::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu this->pw_rho->startz_current, &ucell, &orb); - this->gint_->initialize_pvpR(ucell, &GlobalC::GridD, 1); // always use nspin=1 for transition density + this->gint_->initialize_pvpR(ucell, &this->gd, 1); // always use nspin=1 for transition density // if EXX from scratch, init 2-center integral and calculate Cs, Vs #ifdef __EXX @@ -435,11 +435,26 @@ void LR::ESolver_LR::runner(UnitCell& ucell, const int istep) if (input.lr_solver != "lapack") { pre_op.act(1, offset_is, 1, precondition.data() + offset_is, precondition.data() + offset_is); } } std::cout << "Solving spin-conserving excitation for open-shell system." << std::endl; - HamiltULR hulr(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, *this->psi_ks, this->eig_ks, + HamiltULR hulr(xc_kernel, + nspin, + this->nbasis, + this->nocc, + this->nvirt, + this->ucell, + orb_cutoff_, + this->gd, + *this->psi_ks, + this->eig_ks, #ifdef __EXX - this->exx_lri, this->exx_info.info_global.hybrid_alpha, + this->exx_lri, + this->exx_info.info_global.hybrid_alpha, #endif - this->gint_, this->pot, this->kv, this->paraX_, this->paraC_, this->paraMat_); + this->gint_, + this->pot, + this->kv, + this->paraX_, + this->paraC_, + this->paraMat_); LR::HSolver::solve(hulr, this->X[0].template data(), nloc_per_band, nstates, this->pelec->ekb.c, this->input.lr_solver, this->input.lr_thr, precondition); if (input.out_wfc_lr) { write_states("openshell", this->pelec->ekb.c, this->X[0].template data(), nloc_per_band, nstates); } } @@ -451,12 +466,29 @@ void LR::ESolver_LR::runner(UnitCell& ucell, const int istep) for (int is = 0;is < nspin;++is) { std::cout << "Calculating " << spin_types[is] << " excitations" << std::endl; - HamiltLR hlr(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, *this->psi_ks, this->eig_ks, + HamiltLR hlr(xc_kernel, + nspin, + this->nbasis, + this->nocc, + this->nvirt, + this->ucell, + orb_cutoff_, + this->gd, + *this->psi_ks, + this->eig_ks, #ifdef __EXX - this->exx_lri, this->exx_info.info_global.hybrid_alpha, + this->exx_lri, + this->exx_info.info_global.hybrid_alpha, #endif - this->gint_, this->pot[is], this->kv, this->paraX_, this->paraC_, this->paraMat_, - spin_types[is], input.ri_hartree_benchmark, (input.ri_hartree_benchmark == "aims" ? input.aims_nbasis : std::vector({}))); + this->gint_, + this->pot[is], + this->kv, + this->paraX_, + this->paraC_, + this->paraMat_, + spin_types[is], + input.ri_hartree_benchmark, + (input.ri_hartree_benchmark == "aims" ? input.aims_nbasis : std::vector({}))); // solve the Casida equation LR::HSolver::solve(hlr, this->X[is].template data(), nloc_per_band, nstates, this->pelec->ekb.c + is * nstates, this->input.lr_solver, this->input.lr_thr, precondition/*, @@ -505,10 +537,24 @@ void LR::ESolver_LR::after_all_runners(UnitCell& ucell) auto spin_types = (nspin == 2 && !openshell) ? std::vector({ "singlet", "triplet" }) : std::vector({ "updown" }); for (int is = 0;is < this->X.size();++is) { - LR_Spectrum spectrum(nspin, this->nbasis, this->nocc, this->nvirt, this->gint_, *this->pw_rho, *this->psi_ks, - this->ucell, this->kv, GlobalC::GridD, this->orb_cutoff_, - this->paraX_, this->paraC_, this->paraMat_, - &this->pelec->ekb.c[is * nstates], this->X[is].template data(), nstates, openshell); + LR_Spectrum spectrum(nspin, + this->nbasis, + this->nocc, + this->nvirt, + this->gint_, + *this->pw_rho, + *this->psi_ks, + this->ucell, + this->kv, + this->gd, + this->orb_cutoff_, + this->paraX_, + this->paraC_, + this->paraMat_, + &this->pelec->ekb.c[is * nstates], + this->X[is].template data(), + nstates, + openshell); spectrum.transition_analysis(spin_types[is]); spectrum.optical_absorption(freq, input.abs_broadening, spin_types[is]); } diff --git a/source/module_lr/esolver_lrtd_lcao.h b/source/module_lr/esolver_lrtd_lcao.h index 05e2e7507f..a4ff3f85f5 100644 --- a/source/module_lr/esolver_lrtd_lcao.h +++ b/source/module_lr/esolver_lrtd_lcao.h @@ -49,6 +49,7 @@ namespace LR protected: const Input_para& input; const UnitCell& ucell; + Grid_Driver gd; std::vector orb_cutoff_; // not to use ElecState because 2-particle state is quite different from 1-particle state. diff --git a/source/module_lr/operator_casida/operator_lr_hxc.cpp b/source/module_lr/operator_casida/operator_lr_hxc.cpp index 52ba4fc751..22af349ca1 100644 --- a/source/module_lr/operator_casida/operator_lr_hxc.cpp +++ b/source/module_lr/operator_casida/operator_lr_hxc.cpp @@ -123,7 +123,7 @@ namespace LR // LR_Util::print_HR(*this->gint->get_hRGint(), this->ucell.nat, "VR(grid)"); HR_real_imag.set_zero(); - this->gint->transfer_pvpR(&HR_real_imag, &ucell, &GlobalC::GridD); + this->gint->transfer_pvpR(&HR_real_imag, &ucell, &this->gd); // LR_Util::print_HR(HR_real_imag, this->ucell.nat, "VR(real, 2d)"); LR_Util::set_HR_real_imag_part(HR_real_imag, *this->hR, ucell.nat, type); }; diff --git a/source/module_rdmft/rdmft.cpp b/source/module_rdmft/rdmft.cpp index 3a7f515c9b..7e68635869 100644 --- a/source/module_rdmft/rdmft.cpp +++ b/source/module_rdmft/rdmft.cpp @@ -55,15 +55,16 @@ RDMFT::~RDMFT() } template -void RDMFT::init(Gint_Gamma& GG_in, - Gint_k& GK_in, - Parallel_Orbitals& ParaV_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, + Grid_Driver& gd_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; @@ -74,6 +75,7 @@ void RDMFT::init(Gint_Gamma& GG_in, charge = pelec_in.charge; pelec = &pelec_in; orb = &orb_in; + gd = &gd_in; two_center_bundle = &two_center_bundle_in; XC_func_rdmft = XC_func_rdmft_in; alpha_power = alpha_power_in; diff --git a/source/module_rdmft/rdmft.h b/source/module_rdmft/rdmft.h index 5f4c5848ef..2d9861abf0 100644 --- a/source/module_rdmft/rdmft.h +++ b/source/module_rdmft/rdmft.h @@ -55,7 +55,7 @@ class RDMFT elecstate::ElecState* pelec = nullptr; //! update after ion step - const K_Vectors* kv = nullptr; + const K_Vectors* kv = nullptr; int nk_total = 0; int nbands_total; @@ -81,8 +81,17 @@ class RDMFT // std::vector E_RDMFT(4); //! initialization of rdmft calculation - void 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 init(Gint_Gamma& GG_in, + Gint_k& GK_in, + Parallel_Orbitals& ParaV_in, + UnitCell& ucell_in, + Grid_Driver& gd_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); //! update in ion-step and get V_TV void update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis_in, @@ -189,6 +198,7 @@ class RDMFT // update after ion step const UnitCell* ucell = nullptr; + Grid_Driver* gd = nullptr; const ModulePW::PW_Basis* rho_basis = nullptr; const ModuleBase::matrix* vloc = nullptr; const ModuleBase::ComplexMatrix* sf = nullptr; diff --git a/source/module_rdmft/rdmft_pot.cpp b/source/module_rdmft/rdmft_pot.cpp index f405364abe..06d9c5a3ff 100644 --- a/source/module_rdmft/rdmft_pot.cpp +++ b/source/module_rdmft/rdmft_pot.cpp @@ -50,64 +50,56 @@ template void RDMFT::cal_V_TV() { HR_TV->set_zero(); - - V_ekinetic_potential = new hamilt::EkineticNew>( - hsk_TV, - kv->kvec_d, - HR_TV, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - two_center_bundle->kinetic_orb.get() - ); - - V_nonlocal = new hamilt::NonlocalNew>( - hsk_TV, - kv->kvec_d, - HR_TV, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - two_center_bundle->overlap_orb_beta.get() - ); + + V_ekinetic_potential = new hamilt::EkineticNew>(hsk_TV, + kv->kvec_d, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + two_center_bundle->kinetic_orb.get()); + + V_nonlocal = new hamilt::NonlocalNew>(hsk_TV, + kv->kvec_d, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + two_center_bundle->overlap_orb_beta.get()); if( PARAM.inp.gamma_only ) { - V_local = new rdmft::Veff_rdmft( - GG, - hsk_TV, - kv->kvec_d, - this->pelec->pot, - HR_TV, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "local" - ); + V_local = new rdmft::Veff_rdmft(GG, + hsk_TV, + kv->kvec_d, + this->pelec->pot, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "local"); } else { - V_local = new rdmft::Veff_rdmft( - GK, - hsk_TV, - kv->kvec_d, - this->pelec->pot, - HR_TV, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "local" - ); + V_local = new rdmft::Veff_rdmft(GK, + hsk_TV, + kv->kvec_d, + this->pelec->pot, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "local"); } // update HR_TV in ion-step, now HR_TV has the HR of V_ekinetic + V_nonlcao + V_local @@ -125,42 +117,38 @@ void RDMFT::cal_V_hartree() if( PARAM.inp.gamma_only ) { - V_hartree = new rdmft::Veff_rdmft( - GG, - hsk_hartree, - kv->kvec_d, - this->pelec->pot, - HR_hartree, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "hartree" - ); + V_hartree = new rdmft::Veff_rdmft(GG, + hsk_hartree, + kv->kvec_d, + this->pelec->pot, + HR_hartree, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "hartree"); } else { // this can be optimized, use potHartree.update_from_charge() - V_hartree = new rdmft::Veff_rdmft( - GK, - hsk_hartree, - kv->kvec_d, - this->pelec->pot, - HR_hartree, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "hartree" - ); + V_hartree = new rdmft::Veff_rdmft(GK, + hsk_hartree, + kv->kvec_d, + this->pelec->pot, + HR_hartree, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "hartree"); } // in gamma only, must calculate HR_hartree before HR_local @@ -182,7 +170,7 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) // elecstate::DensityMatrix DM_test(ParaV, nspin, kv->kvec_d, nk_total); // elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test); - // DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + // DM_test.init_DMR(this->gd, &GlobalC::ucell); // DM_test.cal_DMR(); // // compare DM_XC and DM get in update_charge(or ABACUS) @@ -209,46 +197,42 @@ void RDMFT::cal_V_XC(const UnitCell& ucell) if( PARAM.inp.gamma_only ) { // this can be optimized, use potXC.update_from_charge() - V_dft_XC = new rdmft::Veff_rdmft( - GG, - hsk_dft_XC, - kv->kvec_d, - this->pelec->pot, - HR_dft_XC, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "xc", - &etxc, - &vtxc - ); + V_dft_XC = new rdmft::Veff_rdmft(GG, + hsk_dft_XC, + kv->kvec_d, + this->pelec->pot, + HR_dft_XC, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "xc", + &etxc, + &vtxc); } else { // this can be optimized, use potXC.update_from_charge() - V_dft_XC = new rdmft::Veff_rdmft( - GK, - hsk_dft_XC, - kv->kvec_d, - this->pelec->pot, - HR_dft_XC, - &GlobalC::ucell, - orb->cutoffs(), - &GlobalC::GridD, - nspin, - charge, - rho_basis, - vloc, - sf, - "xc", - &etxc, - &vtxc - ); + V_dft_XC = new rdmft::Veff_rdmft(GK, + hsk_dft_XC, + kv->kvec_d, + this->pelec->pot, + HR_dft_XC, + &GlobalC::ucell, + orb->cutoffs(), + this->gd, + nspin, + charge, + rho_basis, + vloc, + sf, + "xc", + &etxc, + &vtxc); } V_dft_XC->contributeHR(); } diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index 4f2d329377..db5f47b6da 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -98,7 +98,7 @@ void RDMFT::update_charge() // calculate DMK and DMR elecstate::DensityMatrix DM_gamma_only(ParaV, nspin); elecstate::cal_dm_psi(ParaV, wg, wfc, DM_gamma_only); - DM_gamma_only.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + DM_gamma_only.init_DMR(this->gd, &GlobalC::ucell); DM_gamma_only.cal_DMR(); for (int is = 0; is < nspin; is++) @@ -128,7 +128,7 @@ void RDMFT::update_charge() // calculate DMK and DMR elecstate::DensityMatrix DM(ParaV, nspin, kv->kvec_d, nk_total); elecstate::cal_dm_psi(ParaV, wg, wfc, DM); - DM.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + DM.init_DMR(this->gd, &GlobalC::ucell); DM.cal_DMR(); for (int is = 0; is < nspin; is++) From 07c39e8986010781f44dd612e33975a4b2337299 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Wed, 11 Dec 2024 18:51:44 +0800 Subject: [PATCH 3/4] convert tool_title.cpp --- source/module_base/tool_title.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/source/module_base/tool_title.cpp b/source/module_base/tool_title.cpp index 922d144a42..6d2716670e 100644 --- a/source/module_base/tool_title.cpp +++ b/source/module_base/tool_title.cpp @@ -17,10 +17,10 @@ namespace ModuleBase void TITLE(const std::string &class_name,const std::string &function_name,const bool disable) { - // if (disable) - // { - // return;//no output - // } + if (disable) + { + return; // no output + } #ifdef __NORMAL std::cout<<" ==> "< "< Date: Thu, 12 Dec 2024 02:02:08 +0800 Subject: [PATCH 4/4] remove useless destructor and flag to enable the default move assignment --- .../module_neighbor/sltk_atom_arrange.cpp | 28 +--------------- .../module_neighbor/sltk_atom_arrange.h | 11 ------- .../module_cell/module_neighbor/sltk_grid.cpp | 29 ++--------------- .../module_cell/module_neighbor/sltk_grid.h | 32 +++---------------- .../module_neighbor/sltk_grid_driver.h | 4 ++- .../test/sltk_atom_arrange_test.cpp | 1 - .../module_neighbor/test/sltk_grid_test.cpp | 4 --- source/module_esolver/esolver_lj.cpp | 9 ------ .../operator_lcao/test/tmp_mocks.cpp | 3 +- source/module_lr/esolver_lrtd_lcao.cpp | 2 ++ 10 files changed, 15 insertions(+), 108 deletions(-) diff --git a/source/module_cell/module_neighbor/sltk_atom_arrange.cpp b/source/module_cell/module_neighbor/sltk_atom_arrange.cpp index d4d020c5ff..8c259794f0 100644 --- a/source/module_cell/module_neighbor/sltk_atom_arrange.cpp +++ b/source/module_cell/module_neighbor/sltk_atom_arrange.cpp @@ -162,30 +162,4 @@ void atom_arrange::search( } return; -} - - -//2015-05-07 -void atom_arrange::delete_vector( - std::ofstream &ofs_in, - const bool pbc_flag, // GlobalV::SEARCH_PBC - Grid_Driver &grid_d, - const UnitCell &ucell, - const double &search_radius_bohr, - const int &test_atom_in) -{ - const double radius_lat0unit2 = search_radius_bohr / ucell.lat0; - - Atom_input at2( - ofs_in, - ucell, - ucell.nat, - ucell.ntype, - pbc_flag, - radius_lat0unit2, - test_atom_in); - - grid_d.delete_vector(at2.getGrid_layerX_minus(),at2.getGrid_layerY_minus(),at2.getGrid_layerZ_minus()); - - grid_d.delete_Cell(); -} +} \ No newline at end of file diff --git a/source/module_cell/module_neighbor/sltk_atom_arrange.h b/source/module_cell/module_neighbor/sltk_atom_arrange.h index 14a2771487..1d263adb3f 100644 --- a/source/module_cell/module_neighbor/sltk_atom_arrange.h +++ b/source/module_cell/module_neighbor/sltk_atom_arrange.h @@ -30,17 +30,6 @@ class atom_arrange const double& rcutmax_Beta, const bool gamma_only_local); - //2015-05-07 - static void delete_vector( - std::ofstream &ofs_in, - const bool pbc_flag, - Grid_Driver &grid_d, - const UnitCell &ucell, - const double &search_radius_bohr, - const int &test_atom_in); - -private: - }; #endif diff --git a/source/module_cell/module_neighbor/sltk_grid.cpp b/source/module_cell/module_neighbor/sltk_grid.cpp index 117a5281e1..52d9cdb760 100644 --- a/source/module_cell/module_neighbor/sltk_grid.cpp +++ b/source/module_cell/module_neighbor/sltk_grid.cpp @@ -16,20 +16,9 @@ CellSet::CellSet() in_grid[2] = 0; } -Grid::Grid(const int& test_grid_in) : test_grid(test_grid_in) -{ - // ModuleBase::TITLE("Grid","Grid"); - //---------------------------------------------------------- - // EXPLAIN : init_cell_flag (use this flag in case memory - // leak) - //---------------------------------------------------------- - init_cell_flag = false; -} +Grid::Grid(const int& test_grid_in) : test_grid(test_grid_in) {} -Grid::~Grid() -{ - this->delete_Cell(); -} +Grid::~Grid() {} void Grid::init(std::ofstream& ofs_in, const UnitCell& ucell, const Atom_input& input) { @@ -49,7 +38,6 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs { ModuleBase::TITLE("SLTK_Grid", "setMemberVariables"); - this->delete_Cell(); // mohan add 2010-09-05 // AdjacentSet::call_times = 0; @@ -96,8 +84,6 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, // output data to ofs Cell[i][j].resize(cell_nz); } } - this->init_cell_flag = true; - this->true_cell_x = input.getGrid_layerX_minus(); this->true_cell_y = input.getGrid_layerY_minus(); this->true_cell_z = input.getGrid_layerZ_minus(); @@ -508,14 +494,3 @@ void Grid::Construct_Adjacent_final(const int i, fatom1.addAdjacent(fatom2); } } -// 2015-05-07 -void Grid::delete_vector(int i, int j, int k) -{ - if (expand_flag) - { - if (this->pbc) - { - this->Cell[i][j][k].atom_map.clear(); - } - } -} diff --git a/source/module_cell/module_neighbor/sltk_grid.h b/source/module_cell/module_neighbor/sltk_grid.h index 56a8b6ee3c..6b6883abe8 100644 --- a/source/module_cell/module_neighbor/sltk_grid.h +++ b/source/module_cell/module_neighbor/sltk_grid.h @@ -41,10 +41,9 @@ class Grid Grid(const int& test_grid_in); virtual ~Grid(); - void init(std::ofstream& ofs, const UnitCell& ucell, const Atom_input& input); + Grid& operator=(Grid&&) = default; - // 2015-05-07 - void delete_vector(int i, int j, int k); + void init(std::ofstream& ofs, const UnitCell& ucell, const Atom_input& input); // Data bool pbc; // periodic boundary condition @@ -64,28 +63,7 @@ class Grid int true_cell_z; std::vector>> Cell; // dx , dy ,dz is cell number in each direction,respectly. - void delete_Cell() // it will replace by container soon! - { - if (this->init_cell_flag) - { - for (int i = 0; i < this->cell_nx; i++) - { - for (int j = 0; j < this->cell_ny; j++) - { - this->Cell[i][j].clear(); - } - } - - for (int i = 0; i < this->cell_nx; i++) - { - this->Cell[i].clear(); - } - - this->Cell.clear(); - this->init_cell_flag = false; - } - } - bool init_cell_flag = false; + // LiuXh add 2019-07-15 double getD_minX() const { @@ -125,8 +103,8 @@ class Grid return true_cell_z; } - private: - const int test_grid; +private: + int test_grid = 0; //========================================================== // MEMBER FUNCTIONS : // Three Main Steps: diff --git a/source/module_cell/module_neighbor/sltk_grid_driver.h b/source/module_cell/module_neighbor/sltk_grid_driver.h index 32ccf482af..4ea18c70de 100644 --- a/source/module_cell/module_neighbor/sltk_grid_driver.h +++ b/source/module_cell/module_neighbor/sltk_grid_driver.h @@ -55,6 +55,8 @@ class Grid_Driver : public Grid ~Grid_Driver(); + Grid_Driver& operator=(Grid_Driver&&) = default; + //========================================================== // EXPLAIN FOR default parameter `adjs = nullptr` // @@ -103,7 +105,7 @@ class Grid_Driver : public Grid private: mutable AdjacentAtomInfo adj_info; - const int test_deconstructor; // caoyu reconst 2021-05-24 + int test_deconstructor = 0; //========================================================== // MEMBER FUNCTIONS : diff --git a/source/module_cell/module_neighbor/test/sltk_atom_arrange_test.cpp b/source/module_cell/module_neighbor/test/sltk_atom_arrange_test.cpp index 369c332bf1..31709e540a 100644 --- a/source/module_cell/module_neighbor/test/sltk_atom_arrange_test.cpp +++ b/source/module_cell/module_neighbor/test/sltk_atom_arrange_test.cpp @@ -52,7 +52,6 @@ Magnetism::~Magnetism() void SetGlobalV() { PARAM.input.test_grid = false; - PARAM.input.test_deconstructor = false; } class SltkAtomArrangeTest : public testing::Test diff --git a/source/module_cell/module_neighbor/test/sltk_grid_test.cpp b/source/module_cell/module_neighbor/test/sltk_grid_test.cpp index af9492e4bc..25b5cf7204 100644 --- a/source/module_cell/module_neighbor/test/sltk_grid_test.cpp +++ b/source/module_cell/module_neighbor/test/sltk_grid_test.cpp @@ -84,7 +84,6 @@ TEST_F(SltkGridTest, Init) Atom_input Atom_inp(ofs, *ucell, ucell->nat, ucell->ntype, pbc, radius, test_atom_in); Grid LatGrid(PARAM.input.test_grid); LatGrid.init(ofs, *ucell, Atom_inp); - EXPECT_TRUE(LatGrid.init_cell_flag); EXPECT_EQ(LatGrid.getCellX(), 11); EXPECT_EQ(LatGrid.getCellY(), 11); EXPECT_EQ(LatGrid.getCellZ(), 11); @@ -126,9 +125,6 @@ TEST_F(SltkGridTest, InitSmall) EXPECT_EQ(LatGrid.cell_nx, 4); EXPECT_EQ(LatGrid.cell_ny, 4); EXPECT_EQ(LatGrid.cell_nz, 4); - // init cell flag - EXPECT_TRUE(LatGrid.init_cell_flag); - ofs.close(); remove("test.out"); } diff --git a/source/module_esolver/esolver_lj.cpp b/source/module_esolver/esolver_lj.cpp index 379aaced17..ccf8d33c30 100644 --- a/source/module_esolver/esolver_lj.cpp +++ b/source/module_esolver/esolver_lj.cpp @@ -88,15 +88,6 @@ void ESolver_LJ::runner(UnitCell& ucell, const int istep) lj_virial(i, j) /= (2.0 * ucell.omega); } } -#ifdef __MPI - atom_arrange::delete_vector( - GlobalV::ofs_running, - PARAM.inp.search_pbc, - grid_neigh, - ucell, - search_radius, - PARAM.inp.test_atom_input); -#endif } double ESolver_LJ::cal_energy() diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp index 38f09cf72c..9f5bb551eb 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/tmp_mocks.cpp @@ -174,7 +174,8 @@ Grid::Grid(const int& test_grid_in) : test_grid(test_grid_in) {} Grid::~Grid() {} Grid_Driver::Grid_Driver(const int& test_d_in, const int& test_grid_in) - : Grid(test_grid_in), test_deconstructor(test_d_in) {} + : Grid(test_grid_in), test_deconstructor(test_d_in) { +} Grid_Driver::~Grid_Driver() {} // filter_adjs delete not adjacent atoms in adjs diff --git a/source/module_lr/esolver_lrtd_lcao.cpp b/source/module_lr/esolver_lrtd_lcao.cpp index d47e962d72..7b362242fe 100644 --- a/source/module_lr/esolver_lrtd_lcao.cpp +++ b/source/module_lr/esolver_lrtd_lcao.cpp @@ -146,6 +146,8 @@ LR::ESolver_LR::ESolver_LR(ModuleESolver::ESolver_KS_LCAO&& ks_sol throw std::invalid_argument("when lr_solver==spectrum, esolver_type must be set to `lr` to skip the KS calculation."); } + this->gd = std::move(ks_sol.gd); + // xc kernel this->xc_kernel = inp.xc_kernel; std::transform(xc_kernel.begin(), xc_kernel.end(), xc_kernel.begin(), tolower);