From 1f9779ded0bc91414219fcef587fad4880af172f Mon Sep 17 00:00:00 2001 From: Yu Liu <77716030+YuLiu98@users.noreply.github.com> Date: Thu, 12 Dec 2024 16:41:22 +0800 Subject: [PATCH] Refactor: remove GlobalC::GridD (#5720) * Refactor: dd TITLE and timer::tick for esolver * Refactor: remove GlobalC::GridD * convert tool_title.cpp * remove useless destructor and flag to enable the default move assignment --------- Co-authored-by: maki49 <1579492865@qq.com> --- source/module_base/tool_title.cpp | 4 +- .../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.cpp | 4 - .../module_neighbor/sltk_grid_driver.h | 9 +- .../test/sltk_atom_arrange_test.cpp | 1 - .../module_neighbor/test/sltk_grid_test.cpp | 4 - .../module_dm/test/test_dm_R_init.cpp | 2 +- source/module_esolver/esolver_fp.cpp | 2 + source/module_esolver/esolver_gets.cpp | 14 +- source/module_esolver/esolver_ks.cpp | 6 +- source/module_esolver/esolver_ks_lcao.cpp | 196 ++++++++------- source/module_esolver/esolver_ks_lcao.h | 2 + .../module_esolver/esolver_ks_lcao_tddft.cpp | 6 + source/module_esolver/esolver_ks_pw.cpp | 8 + source/module_esolver/esolver_lj.cpp | 9 - source/module_esolver/esolver_of.cpp | 10 + source/module_esolver/esolver_sdft_pw.cpp | 10 + source/module_esolver/lcao_before_scf.cpp | 19 +- source/module_esolver/lcao_others.cpp | 22 +- .../module_hamilt_lcao/hamilt_lcaodft/FORCE.h | 2 + .../hamilt_lcaodft/FORCE_STRESS.cpp | 172 ++++++------- .../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 | 31 +-- .../operator_lcao/deepks_lcao.h | 2 + .../operator_lcao/test/tmp_mocks.cpp | 3 +- .../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 | 126 +++++++--- 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 +- 55 files changed, 708 insertions(+), 707 deletions(-) diff --git a/source/module_base/tool_title.cpp b/source/module_base/tool_title.cpp index e55f966f27..6d2716670e 100644 --- a/source/module_base/tool_title.cpp +++ b/source/module_base/tool_title.cpp @@ -19,7 +19,7 @@ void TITLE(const std::string &class_name,const std::string &function_name,const { if (disable) { - return;//no output + return; // no output } #ifdef __NORMAL std::cout<<" ==> "< "<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.cpp b/source/module_cell/module_neighbor/sltk_grid_driver.cpp index 8d770ff1f9..b4f3616fed 100644 --- a/source/module_cell/module_neighbor/sltk_grid_driver.cpp +++ b/source/module_cell/module_neighbor/sltk_grid_driver.cpp @@ -8,10 +8,6 @@ #include #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..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 : @@ -125,9 +127,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_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_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_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_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.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 3d089e568e..44e8105dc9 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -33,6 +33,7 @@ //--------------temporary---------------------------- #include "module_base/global_function.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_elecstate/cal_ux.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag @@ -40,7 +41,6 @@ #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" -#include "module_elecstate/cal_ux.h" #include #ifdef __EXX @@ -65,6 +65,7 @@ // test RDMFT #include "module_rdmft/rdmft.h" + #include namespace ModuleESolver @@ -164,7 +165,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // 5) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell , pv, this->kv.get_nks(), orb_); + LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell, pv, this->kv.get_nks(), orb_); #ifdef __EXX // 6) initialize exx @@ -178,12 +179,12 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // initialize 2-center radial tables for EXX-LRI if (GlobalC::exx_info.info_ri.real_number) { - this->exx_lri_double->init(MPI_COMM_WORLD, ucell,this->kv, orb_); + this->exx_lri_double->init(MPI_COMM_WORLD, ucell, this->kv, orb_); this->exd->exx_before_all_runners(this->kv, ucell, this->pv); } else { - this->exx_lri_complex->init(MPI_COMM_WORLD, ucell,this->kv, orb_); + this->exx_lri_complex->init(MPI_COMM_WORLD, ucell, this->kv, orb_); this->exc->exx_before_all_runners(this->kv, ucell, this->pv); } } @@ -198,7 +199,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } // 8) initialize ppcell - this->ppcell.init_vloc(ucell,this->pw_rho); + this->ppcell.init_vloc(ucell, this->pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); // 9) inititlize the charge density @@ -257,10 +258,19 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } // 14) initialize rdmft, added by jghan - if( PARAM.inp.rdmft == true ) + 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 +306,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 +469,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 +495,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 @@ -647,7 +658,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const GlobalC::dftu.set_dmr(dynamic_cast*>(this->pelec)->get_DM()); } // Calculate U and J if Yukawa potential is used - GlobalC::dftu.cal_slater_UJ(ucell,this->pelec->charge->rho, this->pw_rho->nrxx); + GlobalC::dftu.cal_slater_UJ(ucell, this->pelec->charge->rho, this->pw_rho->nrxx); } #ifdef __DEEPKS @@ -706,21 +717,21 @@ void ESolver_KS_LCAO::hamilt2density_single(UnitCell& ucell, int istep, if (PARAM.inp.sc_mag_switch) { spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); - if(!sc.mag_converged() && this->drho>0 && this->drho < PARAM.inp.sc_scf_thr) + if (!sc.mag_converged() && this->drho > 0 && this->drho < PARAM.inp.sc_scf_thr) { // optimize lambda to get target magnetic moments, but the lambda is not near target - sc.run_lambda_loop(iter-1); + sc.run_lambda_loop(iter - 1); sc.set_mag_converged(true); skip_solve = true; } - else if(sc.mag_converged()) + else if (sc.mag_converged()) { // optimize lambda to get target magnetic moments, but the lambda is not near target - sc.run_lambda_loop(iter-1); + sc.run_lambda_loop(iter - 1); skip_solve = true; } } - if(!skip_solve) + if (!skip_solve) { hsolver::HSolverLCAO hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver); hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge); @@ -772,7 +783,6 @@ void ESolver_KS_LCAO::update_pot(UnitCell& ucell, const int istep, const { this->pelec->cal_converged(); } - } //------------------------------------------------------------------------------ @@ -800,9 +810,14 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& { const std::vector>& tmp_dm = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - ModuleDFTU::dftu_cal_occup_m(iter, ucell,tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt); + ModuleDFTU::dftu_cal_occup_m(iter, + ucell, + tmp_dm, + this->kv, + this->p_chgmix->get_mixing_beta(), + this->p_hamilt); } - GlobalC::dftu.cal_energy_correction(ucell,istep); + GlobalC::dftu.cal_energy_correction(ucell, istep); } GlobalC::dftu.output(ucell); } @@ -903,12 +918,14 @@ 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) { this->pelec->cal_tau(*(this->psi)); } - + //! 2) call after_scf() of ESolver_KS ESolver_KS::after_scf(ucell, istep); @@ -978,27 +995,27 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) if (PARAM.inp.out_mat_hs[0]) { ModuleIO::save_mat(istep, - h_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "H", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); + h_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "H", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); ModuleIO::save_mat(istep, - s_mat.p, - PARAM.globalv.nlocal, - bit, - PARAM.inp.out_mat_hs[1], - 1, - PARAM.inp.out_app_flag, - "S", - "data-" + std::to_string(ik), - this->pv, - GlobalV::DRANK); + s_mat.p, + PARAM.globalv.nlocal, + bit, + PARAM.inp.out_mat_hs[1], + 1, + PARAM.inp.out_app_flag, + "S", + "data-" + std::to_string(ik), + this->pv, + GlobalV::DRANK); } } } @@ -1019,25 +1036,24 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) #ifdef __DEEPKS if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) { - hamilt::HamiltLCAO* p_ham_deepks - = dynamic_cast*>(this->p_hamilt); + hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); std::shared_ptr ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {}); LCAO_Deepks_Interface LDI(ld_shared_ptr); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); LDI.out_deepks_labels(this->pelec->f_en.etot, - this->pelec->klist->get_nks(), - ucell.nat, - PARAM.globalv.nlocal, - this->pelec->ekb, - this->pelec->klist->kvec_d, - ucell, - orb_, - GlobalC::GridD, - &(this->pv), - *(this->psi), - dynamic_cast*>(this->pelec)->get_DM(), - p_ham_deepks); + this->pelec->klist->get_nks(), + ucell.nat, + PARAM.globalv.nlocal, + this->pelec->ekb, + this->pelec->klist->kvec_d, + ucell, + orb_, + this->gd, + &(this->pv), + *(this->psi), + dynamic_cast*>(this->pelec)->get_DM(), + p_ham_deepks); ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); } @@ -1045,29 +1061,28 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) //! 9) Perform RDMFT calculations /******** test RDMFT *********/ - if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17 + if (PARAM.inp.rdmft == true) // rdmft, added by jghan, 2024-10-17 { ModuleBase::matrix occ_number_ks(this->pelec->wg); - for(int ik=0; ik < occ_number_ks.nr; ++ik) - { - for(int inb=0; inb < occ_number_ks.nc; ++inb) + for (int ik = 0; ik < occ_number_ks.nr; ++ik) + { + for (int inb = 0; inb < occ_number_ks.nc; ++inb) { occ_number_ks(ik, inb) /= this->kv.wk[ik]; } } - this->rdmft_solver.update_elec(ucell,occ_number_ks, *(this->psi)); + this->rdmft_solver.update_elec(ucell, occ_number_ks, *(this->psi)); - //! initialize the gradients of Etotal with respect to occupation numbers and wfc, - //! and set all elements to 0. + //! initialize the gradients of Etotal with respect to occupation numbers and wfc, + //! and set all elements to 0. ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true); - psi::Psi dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); + psi::Psi dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); dE_dWfc.zero_out(); double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc); } /******** test RDMFT *********/ - #ifdef __EXX // 10) Write RPA information. if (PARAM.inp.rpa) @@ -1081,7 +1096,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) this->kv, orb_); rpa_lri_double.init(MPI_COMM_WORLD, this->kv, orb_.cutoffs()); - rpa_lri_double.out_for_RPA(ucell,this->pv, *(this->psi), this->pelec); + rpa_lri_double.out_for_RPA(ucell, this->pv, *(this->psi), this->pelec); } #endif @@ -1134,19 +1149,28 @@ 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); } } //! 14) Print out atomic magnetization only when 'spin_constraint' is on. - if (PARAM.inp.sc_mag_switch) + if (PARAM.inp.sc_mag_switch) { spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.cal_mi_lcao(istep); @@ -1154,7 +1178,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) sc.print_Mag_Force(GlobalV::ofs_running); } - //! 15) Clean up RA. + //! 15) Clean up RA. //! this should be last function and put it in the end, mohan request 2024-11-28 if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress) { @@ -1187,7 +1211,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); @@ -1226,12 +1250,12 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.wannier_spin); myWannier.set_tpiba_omega(ucell.tpiba, ucell.omega); 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) @@ -1245,27 +1269,25 @@ 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"); } //! 19) berry phase calculations, added by jingan - if (PARAM.inp.calculation == "nscf" && - berryphase::berry_phase_flag && - ModuleSymmetry::Symmetry::symm_flag != 1) + if (PARAM.inp.calculation == "nscf" && berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != 1) { 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.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + 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"); } + + ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf"); } template class ESolver_KS_LCAO; 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 880daa8512..fc1ba66bc6 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) @@ -296,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, @@ -306,6 +310,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_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_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..9eef0aafce 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); @@ -80,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); @@ -110,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, @@ -127,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 @@ -189,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, @@ -213,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 @@ -347,6 +344,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) 1); } + ModuleBase::timer::tick("ESolver_KS_LCAO", "before_scf"); return; } @@ -377,6 +375,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; } 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 4a262da0e4..c740a96550 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -1,26 +1,26 @@ #include "FORCE_STRESS.h" -#include "module_parameter/parameter.h" #include "module_hamilt_lcao/module_dftu/dftu.h" //Quxin add for DFT+U on 20201029 #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/output_log.h" +#include "module_parameter/parameter.h" // new #include "module_base/timer.h" -#include "module_parameter/parameter.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_elecstate/potentials/efield.h" // liuyu add 2022-05-18 #include "module_elecstate/potentials/gatefield.h" // liuyu add 2022-09-13 #include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10 #include "module_hamilt_general/module_vdw/vdw.h" +#include "module_parameter/parameter.h" #ifdef __DEEPKS #include "module_elecstate/elecstate_lcao.h" -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03 -#include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h" // mohan add 2024-07-22 +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03 +#include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h" // mohan add 2024-07-22 #endif +#include "module_elecstate/elecstate_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dspin_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h" -#include "module_elecstate/elecstate_lcao.h" template Force_Stress_LCAO::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : RA(&ra), f_pw(nat_in), nat(nat_in) @@ -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, @@ -179,41 +181,38 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, pv, kv); // calculate force and stress for Nonlocal part - if(PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + 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, ucell, orb, GlobalC::GridD, kv.get_nks(), kv.kvec_d, isstress); + GlobalC::ld.cal_gdmx(dm_gamma, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress); } else { @@ -531,13 +519,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, ->get_DM() ->get_DMK_vector(); - GlobalC::ld.cal_gdmx(dm_k, - ucell, - orb, - GlobalC::GridD, - kv.get_nks(), - kv.kvec_d, - isstress); + GlobalC::ld.cal_gdmx(dm_k, ucell, orb, gd, kv.get_nks(), kv.kvec_d, isstress); } if (PARAM.inp.deepks_out_unittest) { @@ -551,19 +533,17 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } LCAO_deepks_io::save_npy_gvx(ucell.nat, - GlobalC::ld.des_per_atom, - GlobalC::ld.gvx_tensor, - PARAM.globalv.global_out_dir, - GlobalV::MY_RANK); + GlobalC::ld.des_per_atom, + GlobalC::ld.gvx_tensor, + PARAM.globalv.global_out_dir, + GlobalV::MY_RANK); } } else { const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - LCAO_deepks_io::save_npy_f(fcs, - file_fbase, - ucell.nat, - GlobalV::MY_RANK); // no scf, F_base=F_tot + LCAO_deepks_io::save_npy_f(fcs, file_fbase, ucell.nat, + GlobalV::MY_RANK); // no scf, F_base=F_tot } } #endif @@ -732,9 +712,9 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { const std::string file_s = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; LCAO_deepks_io::save_npy_s(scs, - file_s, - ucell.omega, - GlobalV::MY_RANK); // change to energy unit Ry when printing, S_base; + file_s, + ucell.omega, + GlobalV::MY_RANK); // change to energy unit Ry when printing, S_base; } if (PARAM.inp.deepks_scf) { @@ -752,12 +732,11 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } if (PARAM.inp.deepks_out_labels) // not parallelized yet { - const std::string file_s = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_npy_s( - scs, - file_s, - ucell.omega, - GlobalV::MY_RANK); // change to energy unit Ry when printing, S_tot, w/ model + const std::string file_s = PARAM.globalv.global_out_dir + "deepks_stot.npy"; + LCAO_deepks_io::save_npy_s(scs, + file_s, + ucell.omega, + GlobalV::MY_RANK); // change to energy unit Ry when printing, S_tot, w/ model // wenfei add 2021/11/2 if (PARAM.inp.deepks_scf) @@ -767,12 +746,11 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { GlobalC::ld.cal_gvepsl(ucell.nat); - LCAO_deepks_io::save_npy_gvepsl( - ucell.nat, - GlobalC::ld.des_per_atom, - GlobalC::ld.gvepsl_tensor, - PARAM.globalv.global_out_dir, - GlobalV::MY_RANK); // unitless, grad_vepsl + LCAO_deepks_io::save_npy_gvepsl(ucell.nat, + GlobalC::ld.des_per_atom, + GlobalC::ld.gvepsl_tensor, + PARAM.globalv.global_out_dir, + GlobalV::MY_RANK); // unitless, grad_vepsl } } } @@ -869,11 +847,11 @@ void Force_Stress_LCAO::calForcePwPart(const UnitCell& ucell, // local pseudopotential force: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - f_pw.cal_force_loc(ucell,fvl_dvl, rhopw, nlpp.vloc, chr); + f_pw.cal_force_loc(ucell, fvl_dvl, rhopw, nlpp.vloc, chr); //-------------------------------------------------------- // ewald force: use plane wave only. //-------------------------------------------------------- - f_pw.cal_force_ew(ucell,fewalds, rhopw, &sf); // remain problem + f_pw.cal_force_ew(ucell, fewalds, rhopw, &sf); // remain problem //-------------------------------------------------------- // force due to core correlation. @@ -892,6 +870,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 +897,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, isstress, fsr, // mohan add 2024-06-15 ucell, + gd, psi, pelec, foverlap, @@ -943,6 +923,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 +949,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn isstress, fsr, // mohan add 2024-06-16 ucell, + gd, psi, pelec, foverlap, @@ -1009,17 +991,17 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, // local pseudopotential stress: // use charge density; plane wave; local pseudopotential; //-------------------------------------------------------- - sc_pw.stress_loc(ucell,sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr); + sc_pw.stress_loc(ucell, sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr); //-------------------------------------------------------- // hartree term //-------------------------------------------------------- - sc_pw.stress_har(ucell,sigmahar, rhopw, 0, chr); + sc_pw.stress_har(ucell, sigmahar, rhopw, 0, chr); //-------------------------------------------------------- // ewald stress: use plane wave only. //-------------------------------------------------------- - sc_pw.stress_ewa(ucell,sigmaewa, rhopw, 0); // remain problem + sc_pw.stress_ewa(ucell, sigmaewa, rhopw, 0); // remain problem //-------------------------------------------------------- // stress due to core correlation. @@ -1034,7 +1016,7 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, sigmaxc(i, i) = -etxc / ucell.omega; } // Exchange-correlation for PBE - sc_pw.stress_gga(ucell,sigmaxc, rhopw, chr); + sc_pw.stress_gga(ucell, sigmaxc, rhopw, chr); return; } @@ -1042,9 +1024,7 @@ void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, #include "module_base/mathzone.h" // do symmetry for total force template -void Force_Stress_LCAO::forceSymmetry(const UnitCell& ucell, - ModuleBase::matrix& fcs, - ModuleSymmetry::Symmetry* symm) +void Force_Stress_LCAO::forceSymmetry(const UnitCell& ucell, ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm) { double d1, d2, d3; for (int iat = 0; iat < ucell.nat; iat++) 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 38e5be5051..6542f2a0f6 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 47d2f57c9f..476fa43ee7 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 @@ -1,8 +1,8 @@ #include "deepks_lcao.h" -#include "module_parameter/parameter.h" #include "module_base/timer.h" #include "module_base/tool_title.h" +#include "module_parameter/parameter.h" #ifdef __DEEPKS #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #endif @@ -25,11 +25,11 @@ DeePKS>::DeePKS(HS_Matrix_K* hsk_in, const LCAO_Orbitals* ptr_orb, const int& nks_in, elecstate::DensityMatrix* DM_in) - : OperatorLCAO(hsk_in, kvec_d_in, hR_in), - DM(DM_in), ucell(ucell_in), - intor_orb_alpha_(intor_orb_alpha), ptr_orb_(ptr_orb), nks(nks_in) + : OperatorLCAO(hsk_in, kvec_d_in, hR_in), DM(DM_in), ucell(ucell_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 @@ -52,7 +52,7 @@ void hamilt::DeePKS>::initialize_HR(Grid_Driver* Gr ModuleBase::TITLE("DeePKS", "initialize_HR"); ModuleBase::timer::tick("DeePKS", "initialize_HR"); - auto* paraV = this->hR->get_paraV();// get parallel orbitals from HR + auto* paraV = this->hR->get_paraV(); // get parallel orbitals from HR // TODO: if paraV is nullptr, AtomPair can not use paraV for constructor, I will repair it in the future. // this->H_V_delta = new HContainer(paraV); @@ -158,7 +158,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 +186,7 @@ void DeePKS, double>>::contributeHR() { ModuleBase::timer::tick("DeePKS", "contributeHR"); - 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); // calculate dE/dD GlobalC::ld.cal_gedm(this->ucell->nat); @@ -219,7 +219,7 @@ void DeePKS, std::complex>>::contribut { ModuleBase::timer::tick("DeePKS", "contributeHR"); - 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); // calculate dE/dD GlobalC::ld.cal_gedm(this->ucell->nat); @@ -290,9 +290,10 @@ void hamilt::DeePKS>::pre_calculate_nlm( ModuleBase::Vector3 dtau = tau0 - tau1; intor_orb_alpha_->snap(T1, L1, N1, M1, 0, dtau * ucell->lat0, false /*calc_deri*/, nlm); nlm_in[ad].insert({all_indexes[iw1l], nlm[0]}); - if (npol == 2) { + if (npol == 2) + { nlm_in[ad].insert({all_indexes[iw1l + 1], nlm[0]}); -} + } } } } @@ -386,9 +387,10 @@ void hamilt::DeePKS>::calculate_HR() ModuleBase::Vector3& R_index1 = adjs.box[ad1]; auto row_indexes = paraV->get_indexes_row(iat1); const int row_size = row_indexes.size(); - if (row_size == 0) { + if (row_size == 0) + { continue; -} + } std::vector s_1t(trace_alpha_size * row_size); for (int irow = 0; irow < row_size; irow++) @@ -412,9 +414,10 @@ void hamilt::DeePKS>::calculate_HR() hamilt::BaseMatrix* tmp = this->H_V_delta->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); // if not found , skip this pair of atoms - if (tmp == nullptr) { + if (tmp == nullptr) + { continue; -} + } auto col_indexes = paraV->get_indexes_col(iat2); const int col_size = col_indexes.size(); std::vector hr_current(row_size * col_size, 0); 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/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_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..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); @@ -326,11 +328,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 +345,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 +390,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 +437,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 +468,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 +539,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++)